Astronomy & Astrophysics manuscript no. paper \begin{matrix}c \\ ◯ \end{matrix} ESO 2020 November 28, 2020

The effect of a dynamo-generated field on the Parker wind

P. Jakab^{1,2} and A. Brandenburg^{1,3,4,5}

1 Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden 2 Department of Nuclear and Subnuclear Physics, Pavol Jozef Safarik University in Kosice, Srobarova 2, 041 54 Kosice, Slovakia 3 Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden 4 JILA and Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80303, USA 5 McWilliams Center for Cosmology & Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

November 28, 2020, Revision: 1.93

ABSTRACT

Context. Stellar winds are an integral part of the underlying dynamo, the motor of stellar activity. The wind controls the star’s angular momentum loss, which depends on the magnetic field geometry which, in turn, varies significantly in time and latitude. Aims. Here we study basic properties of a self-consistent model that includes simple representations of both the global stellar dynamo in a spherical shell and the exterior in which the wind accelerates and becomes supersonic. Methods. We numerically solve an axisymmetric mean-field model for the induction, momentum, and continuity equations using an isothermal equation of state. The model allows for the simultaneous generation of a mean magnetic field and the development of a Parker wind. The resulting flow is transonic at the critical point, which we arrange to be between the inner and outer radii of the model. The boundary conditions are assumed to be such that the magnetic field is antisymmetric about the equator, i.e., dipolar. Results. At the solar rotation rate, the dynamo is oscillatory and of α^{2} type. In most of the domain, the magnetic field corresponds to that of a split monopole. The magnetic energy flux is largest between the stellar surface and the critical point. The angular momentum flux is highly variable in time and can reach negative values, especially at midlatitudes. At rapid rotation of up to 50 times the solar value, most of the magnetic field is lost along the axis within the inner tangential cylinder of the model. Conclusions. The model reveals unexpected features that are not generally anticipated from models that are designed to reproduce the solar wind: highly variable angular momentum fluxes even from just an α^{2} dynamo in the star. A major caveat of our isothermal models with a magnetic field produced by a dynamo is the difficulty to reach small enough plasma betas without the dynamo itself becoming unrealistically strong inside the star.

Key words. Sun: sunspots – Sun: dynamo – turbulence – magnetohydrodynamics (MHD) – hydrodynamics

1. Introduction

The emergence of a wind around stars is a remarkable and somewhat counter-intuitive phenomenon. The existence of the solar wind was already suggested by the fact that the tails of comets always point away from the Sun (Biermann, 1951). Nevertheless, the wind was thought to be a relatively slow phenomenon associated with an evaporation of the corona (Chamberlain, 1960). The physical nature and mathematical theory of the solar wind was first understood by Parker (1958). His theory showed that the wind starts off as a subsonic flow some distance above the corona. It gradually gains in speed as the gravitational force diminishes and the effective outward pull resulting from the quadratic increase of the cross-sectional area in Bernoulli’s law. This is a purely hydrodynamic phenomenon, unlike what was suggested by the popular notion of the solar corpuscular radiation at the time.

Stellar winds play a crucial role in a star’s life. Without the wind, the Sun would still be spinning rapidly and magnetically superactive. A proper understanding of the rotational evolution of a star through magnetic braking via a wind is important not only for stellar evolution, but it also plays a role in understanding the diversity of magnetic activity as a function of rotation rate and age (van Saders et al., 2016). As the star reaches the age of the Sun, the magnetic field either changes its geometry such that stellar braking is reduced (Metcalfe & van Saders, 2017; See et al., 2019) or it can continue to brake and the star’s differential rotation becomes antisolar-like (Gastine et al., 2014; Käpylä et al., 2014), i.e., the equator spins slower than the poles. Stellar winds can also be important for the dynamo itself in that they can transport magnetic helicity away from the dynamo region, and thereby alleviate what is known as catastrophic quenching; see Mitra et al. (2011) for mean-field models and Del Sordo et al. (2013) for computations of the magnetic helicity flux in simulations in a turbulent wind. Magnetic winds also affect the density and dynamics of cosmic rays in the heliosphere. Computing selfconsistently the dynamo-generated magnetic field evolution in the heliosphere is therefore crucial and for modeling the magnetic shielding of Galactic cosmic rays on the Earth.

The theory of a magnetized stellar wind by Weber & Davis (1967) employes a prescribed and time-independent stellar magnetic field, so any feedback on the underlying dynamics was ignored. This is also true of the recent numerical models of Réville et al. (2015), who compared different magnetic multipoles as initial conditions of their models. This has changed only in recent years. Given that the wind normally dominates over the magnetic field, one can separate the dynamics of the wind from that of the solar dynamo. Pinto et al. (2011), and more recently Perri et al. (2018), modeled this by using two separate codes that are magnetically coupled through a matching condition at the solar surface. In more recent work, Perri et al. (2020) have extended their model to include also a mean-field dynamo solution into the Pluto code, rather than matching the solutions of two separate codes. This allows for a feedback of the wind onto the dynamo.

P. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind

This is therefore similar to the work presented here, except that they still invoke what they call a multilayered boundary condition. This means that different equations are being solved inside and outside the star. The model is therefore still not fully self consistent, but in some ways more realistic than ours.

The purpose of the present paper is to explore some basic properties of stellar winds in the presence of dynamo-generated magnetic fields. It is appropriate to adopt a mean-field model, where we solve the equations for the azimuthally averaged magnetic and velocity fields. In this paper, those mean fields are denoted by an overbar. The effects of turbulence are then parameterized through a turbulent viscosity and a turbulent magnetic diffusivity. In the star’s convection zone, there are also cyclonic convective motions giving rise to kinetic helicity of opposite signs in the two hemispheres. This is modeled through an α effect (Krause & Rädler, 1980). The turbulent magnetic diffusivity is here assumed constant.

The presence of the magnetic field causes the kinetic and magnetic stresses to be different from zero. The turbulent viscosity is itself a result of kinetic and magnetic stresses caused by the fluctuating components of the magnetic and velocity fields. In the theory of turbulent accretion discs (Frank et al., 1992), those stresses are parameterized by the Shakura & Sunyaev (1973) parameter, α_{SS}. It quantifies the stress in terms of the background differential rotation, the sound speed, and the scale height. In accretion discs, where the differential rotation is Keplerian, this amounts to a scaling of the stress by the sound speed squared. In our case, the differential rotation is not related to the sound speed, but the basic mechanism of angular momentum transfer is the same, and we can still express the total stress in a similar fashion.

Unlike the work of Perri et al. (2018), we consider the evolution of the dynamo and the wind within a single code. At this point, our aim is not to produce a realistic model of the Sun, but rather a physically consistent model under conditions where the dynamics of the wind can no longer be separated from that of the dynamo. Our models can also be applied to conditions of rapid rotation, which strongly affects the wind. This can be particularly relevant to young stars in their T Tauri phase. We begin by presenting the basic equations of our model and turn then to the discussion of our results.

The simplest wind solution is the isothermal one that was already found by Parker (1958). Heating is not explicitly invoked. Its physics resembles that of a siphon flow. Once a fluid parcel has moved over the top of the effective gravitational potential, it simply continues to fall and pulls the remaining fluid behind it. The top of the effective potential corresponds to the critical point where the flow speed crosses the sonic point. We arrange this point to be in the middle of the computational domain such that the flow speed becomes supersonic well before the outer point r_{out}. We fit the dynamo-active zone (or stellar envelope) with an α effect different from zero into a spherical shell between the inner point of the computational domain, r_{in}, and a radius R, which models the surface of the star.

The usefulness of an isothermal solution can be justified by considering the fact that the sound speed both at the bottom of the convection zone and in the solar wind is about 100 km s^{−1}, corresponding to a temperature of a million degrees. The lower temperature near the photosphere is obviously ignored. For an isothermal gas, the mean pressure \overline{p} is then simply proportional to the gas density \overline{ρ} with \overline{p} = \overline{ρ}c_{s}^{2}, where c_{s} is the isothermal sound speed. The pressure gradient is then given by (∇\overline{p})/\overline{ρ} = c_{s}^{2}∇ ln \overline{ρ}. The implications of a cool photosphere will be discussed at the end of the paper.

We begin by discussing first the basic equations, boundary conditions, and parameters in Sect. 2. We then present our results in Sect. 3, and draw our conclusions in Sect. 4.

2. The model

We adopt spherical polar coordinates, (r, θ, ϕ), with the origin at the center of the star. The vector r points away from the center, the colatitude θ increases away from the north pole, and ϕ increases in the eastward direction. We assume axisymmetry, i.e., ∂/∂ϕ = 0.

2.1. Basic equations

We write the mean magnetic field as \overline{B} = ∇ × \overline{A}, where \overline{A} is the mean vector potential. This ensures that ∇ · \overline{B} = 0 at all times. The evolution equations for \overline{A}, the mean velocity \overline{U}, and the logarithmic mean density \overline{ρ}, are

\frac{∂\overline{A}}{∂t} = \overline{U} × \overline{B} + α\overline{B} − η_{T}μ_{0}\overline{J},  (1)
\frac{D\overline{U}}{Dt} = −c_{s}^{2}∇ ln \overline{ρ} − \frac{GM}{r^{2}} r̂ + \frac{1}{\overline{ρ}}\overline{J} × \overline{B} − ν_{T}\overline{Q},  (2)
\frac{D ln \overline{ρ}}{Dt} = −∇ · \overline{U},  (3)

where D/Dt = ∂/∂t + \overline{U} · ∇ is the advective derivative, G is Newton’s constant, M is the stellar mass, r̂ = r/r is the radial unit vector, η_{T} and ν_{T} are the sums of turbulent and microphysical values of magnetic diffusivity and kinematic viscosity, respectively, α is the aforementioned coefficient in the α effect, \overline{J} = ∇ × \overline{B}/μ_{0} is the mean current density, μ_{0} is the vacuum permeability,

−\overline{Q} = ∇^{2}\overline{U} + \frac{1}{3}∇∇ · \overline{U} + 2\overline{S} · ∇ ln(ν_{T}\overline{ρ}) (4)

is a term appearing in the viscous force, where \overline{S} is the traceless rate of strain tensor of the mean flow with components \overline{S}_{ij} = \frac{1}{2}(\overline{U}_{i,j} + \overline{U}_{j,i}) − \frac{1}{3}δ_{ij}∇ · \overline{U}. The dot in Eq. (4) denotes the contraction over the free index of ∇ ln(ν_{T}\overline{ρ}).

The mean magnetic field is generated by the α effect. This leads to exponential growth, provided the value of α is above a certain critical value. Eventually, the dynamo must saturate because the Lorentz force from the mean field, \overline{J} × \overline{B}, drives fluid motions that feed back onto the dynamo to limit its growth. This way of achieving saturation is sometimes referred to as Malkus & Proctor (1975) mechanism. In addition, there can be feedback from the small-scale magnetic field that leads to a nonlinear suppression of α, which is referred to as α quenching. We assume here a simple quenching function for α, which is then written in the form

α(r, θ, \overline{B}) = \frac{α_{0}f_{α}(r) cos θ sin^{n} θ}{1 + Q_{α}\overline{B}^{2}/B_{eq}^{2}} , (5)

where n = 6 is chosen to concentrate the α effect to low latitudes (Jabbari et al., 2015; Cole et al., 2016), Q_{α} is a quenching parameter that determines the typical field strength, which is expected to be of the order of Q_{α}^{−1/2}B_{eq}, and

f_{α}(r) = Θ((r − R)/w_{α}) (6)

is a radial profile function with Θ(x) being a smoothened step function from 0 to 1 as x crosses zero. Here, R and w_{α} determine the location and width of the transition. The value of Q_{α} determines the nonlinear equilibration of the dynamo, in addition to the macroscopic feedback from the Lorentz force mentioned above. Our model thus comprises three distinct layers with

r_{in} < R < r_{c} < r_{out}, (7)

where r_{in} < r < R is the dynamo region (modeling the stellar envelope), R < r < r_{c} is the wind acceleration region (modeling the locations of the solar corona and the Alfvén point), and r_{c} < r < r_{out} is the supersonic wind region with r_{c} = GM/2c_{s}^{2} being the critical point.

2.2. Boundary conditions

In most of the cases, we apply a uniform angular velocity Ω_{0} on the inner boundary r = r_{in} by setting u_{ϕ} = r_{in} sin θ Ω_{0}. For the other two velocity components, we adopt “open” boundary conditions by setting the second radial derivative to zero. This condition turns out to be stable in all cases considered in this paper. It allows for a weak inflow to replenish the mass loss on the outer boundary r = r_{out}, where we apply open boundary conditions for all three velocity components. No precautions are taken to ensure that the mass in the computational domain stays constant. It turns out, however, that the total mass remains nearly unchanged. This is, to some extent, also explained by the fact that the total mass loss rate is small compared with other inverse time scales in the problem.

For the magnetic field, we adopt a perfect conductor boundary conditions on the inner radius, i.e.,

\frac{∂\overline{A}_{r}}{∂r} = \overline{A}_{θ} = \overline{A}_{ϕ} = 0 on r = r_{in}, (8)

and a radial field condition on the outer radius, i.e.,

\overline{A}_{r} = \frac{∂\overline{A}_{θ}}{∂r} + \frac{\overline{A}_{θ}}{r} = \frac{∂\overline{A}_{ϕ}}{∂r} + \frac{\overline{A}_{ϕ}}{r} = 0 on r = r_{out}. (9)

On the pole, we assume

\frac{∂\overline{A}_{r}}{∂θ} = \overline{A}_{θ} = \overline{A}_{ϕ} = 0 on θ = 0^{∘} , (10)

while on the equator, we assume

\frac{∂\overline{A}_{r}}{∂θ} = \overline{A}_{θ} = \frac{∂\overline{A}_{ϕ}}{∂θ} = 0 on θ = 90^{∘}. (11)

Since our simulations are axisymmetric, the magnetic field is conveniently represented via \overline{B}_{ϕ} and \overline{A}_{ϕ}. In particular, contours of r sin θ \overline{A}_{ϕ} give the magnetic field lines of the poloidal field, \overline{B}_{pol} = ∇ × (\overline{A}_{ϕ}ϕ̂).

2.3. Wind solution as initial condition

As initial condition for \overline{U} ≡ (u, 0, 0) and \overline{ρ}, we adopt the Parker wind solution. In some cases we also add a finite angular velocity with constant angular momentum, although its effect on the dynamics is ignored in the initial condition. We begin by discussing the Parker wind solution, which can be obtained by solving the Bernoulli equation,

\frac{1}{2}u^{2} + c_{s}^{2} ln \overline{ρ} − GM/r = const, (12)

along with the equation of mass conservation, which states that the mass loss rate is given by Ṁ = 4πr^{2}\overline{ρ}u. We then obtain

\frac{1}{2}u^{2} − c_{s}^{2} ln u − c_{s}^{2} ln r^{2} − GM/r = Φ_{0}, (13) where Φ_{0} = −3/2 is obtained by inserting the values u = r_{c} = 1 for the critical point. We solve the Bernoulli equation iteratively. For r ≤ r_{c}, using u = c_{s}r/r_{c} initially, we iterate

c_{s}^{2} ln u_{i+1}(r) = \frac{1}{2}u_{i}^{2} − c_{s}^{2} ln r^{2} − GM/r − Φ_{0}, (14)

while for r > r_{c}, using u_{0} = 2c_{s} initially, we iterate

\frac{1}{2}u_{i+1}^{2}(r) = c_{s}^{2} ln u_{i} + c_{s}^{2} ln r^{2} + GM/r + Φ_{0}. (15)

This iteration procedure was implemented by Jörn Warnecke and Dhrubaditya Mitra into the Pencil Code^{1} in 2012. We choose the initial value of Ṁ to be Ṁ_{0}.

2.4. Parameters and estimates for the Sun

It is convenient to work with nondimensional units by measuring speeds in units of the isothermal sound speed and lengths in units of the critical radius, r_{c} = GM/2c_{s}^{2}. In the following, we use tildae to denote nondimensional quantities. Using typical numbers for the Sun, we have

c_{s} = 10^{7} cm s^{−1} = 100 km s^{−1}, (16)

GM = GM_{⊙} ≈ 1.3 × 10^{26} cm^{3} s^{−2}, and therefore

r_{c} = GM_{⊙}/2c_{s}^{2} ≈ 7 × 10^{11} cm ≈ 10R_{⊙} ≈ 0.05 AU. (17)

In the Sun, the turbulent viscosity is ν_{T} ≈ u_{rms}ℓ/3 ≈ 10^{13} cm^{2} s^{−1}. The nondimensional viscosity is then

ν̃_{T} ≡ \frac{2ν_{T}c_{s}}{GM_{⊙}} ≈ 2 × 10^{−6}, (18)

which is rather small.

For numerical stability, as already alluded to, we cannot choose the value of ν_{T} to be too small. In practice, for a numerical resolution of 128 × 32 mesh points in the r and θ directions, we can choose ν̃_{T} ≈ 0.01. For 4096 × 1024 mesh points, on the other hand, we can reduce it by a factor of 128 to ν̃_{T} ≈ 8×10^{−5}. This then also means that in the stellar convection zone, we cannot adopt significantly smaller values, as is expected theoretically based on our earlier estimates of u_{rms} and ℓ.

The nondimensional value of the angular velocity is given by

Ω̃ = r_{c}Ω_{0}/c_{s} = GMΩ_{0}/2c_{s}^{3} ≈ 0.2, (19)

where we have used Ω_{0} = 3 × 10^{−6} s^{−1}. The strength of the dynamo is determined by the two dynamo numbers,

C_{α} = α_{0}R/η_{T} and C_{Ω} = ΔΩR^{2}/η_{T}. (20)

The excitation conditions for dipolar and quadrupolar parities are generally fairly close together (Roberts, 1972). This is because the magnetic field is strongest at high latitudes, so the hemispheric coupling is weak. In the following we restrict ourselves to solutions with dipolar parity. We vary the value of C_{α} and focus on values that are about twice supercritical.

In our simulations, we adopt nondimensional units by setting

r_{c} = c_{s} = Ṁ_{0} = μ_{0} = 1, (21)

which implies that GM = 2. Our unit of mass is then [M] = Ṁ_{0}r_{c}/c_{s}. For the Sun, we have Ṁ_{0} ≈ 6 × 10^{12} g s^{−1}, so that our unit of density is [ρ] = Ṁ_{0}/c_{s}r_{c}^{2}, which is about 1.2 × 10^{−18} g cm^{−3} for the Sun. Therefore, our unit of B is [B] = (μ_{0}[ρ])^{1/2}c_{s}, which is about 0.04 G for the Sun. The value of Newton’s constant G never enters on its own. It could be determined a posteriori, if we knew the total stellar mass. In our model, we can compute the mass M_{∗} of the stellar envelope in r_{in} ≤ r ≤ R, but this still leaves the mass of the stellar core undetermined. In the following, it is often convenient to retain the symbols r_{c}, c_{s}, Ṁ_{0}, and μ_{0} to remind ourselves of the normalization.

1 http://github.com/pencil-code  (Pencil  Code Collaboration., 2020), DOI:10.5281/zenodo.2315093

There are a few other parameters of the model that we keep fixed. In all cases we use w_{α} = 0.02 for the transition thickness of α near the surface; see Eq. (6). We always take r_{in} = 0.1 and R = 0.2. This corresponds to a fractional shell thickness of 50% instead of the 30% in the case of the Sun, but we should keep in mind that there are other properties that agree with the Sun only qualitatively. Another example is our smaller choice of R/r_{in} = 5 instead of the solar value of about 10. In all our simulations with 4096 × 1024 meshpoints, we use ν̃_{T} = 8 × 10^{−5}.

2.5. Comparison of characteristic time scales

In our simulations, sound speed and the critical radius are set to unity, so the characteristic sound travel time,

τ_{s} = r_{c}/c_{s} (22)

is therefore also unity. When we adopt the stellar rotation rate, Ω̃ = 0.2, the corresponding rotational time scale

τ_{Ω} = Ω_{0}^{−1} (23)

is then five, and the rotation period is 2π/Ω̃ ≈ 30. The characteristic time scale for the dynamo is the turbulent diffusive time (e.g., Stix, 1974),

τ_{TD} = R^{2}/η_{T}, (24)

which is around 500 in our models. Another interesting time scale for our models is the mass loss time,

τ_{massloss} = M/⟨Ṁ⟩ ≈ M_{0}/Ṁ_{0}. (25)

In our models, M_{0} ≈ 7000 and Ṁ_{0} = 1, so τ_{massloss} ≈ 7000. It turns out that the spindown time is of a similar order of magnitude. It is given by

τ_{spindown} = J_{∗}/⟨J̇⟩, (26)

where J_{∗} = ∫_{∗} \overline{ρ}ϖ^{2}Ω dV is the angular momentum of the stellar envelope, with ϖ = r sin θ being the cylindrical radius, Ω = \overline{U}_{ϕ}/ϖ is the local angular velocity, and J̇ is the angular momentum loss, which we calculate in Sect. 3.5. The asterisk on the integral denotes the volume of the envelope. The mass loss and spindown times are the longest among the time scales considered here, so the mass in the envelope cannot change significantly during the time scales of interest for the wind and the dynamo.

3. Results

After some preliminary studies at low resolution of 128 × 32 meshpoints with ν_{T} = η_{T} = 10^{−2}r_{c}c_{s}, we performed highresolution simulations with 4096 × 1024 meshpoints, where we were able to decrease ν_{T} and η_{T} to 8 × 10^{−5}r_{c}c_{s}. These values are still above the physically motivated value, but for numerical stability reasons, they cannot be decreased further without invoking artificial viscosity and magnetic diffusivity.

Fig. 1. Radial dependence of (a) Ṁ for different latitude ranges, and (b) M_{r} (solid line) for Model A. The dotted line in (b) refers to the mass within the computational domain only, so it vanishes on r = r_{out}.

Fig. 2. Time series of the three magnetic field components at one point for Model A. Here, \overline{B}_{θ} and \overline{B}_{ϕ} are multiplied by 25 to make those components better visible. Note that all three components of \overline{B} are asymmetric about zero. The 12 long tick marks on the lower abscissa show the times for which snapshots will be discussed later on.

Our main model is called Model A, which has the solar value of Ω and a minimal amount of viscosity and magnetic diffusivity that can still be tolerated. Later, we also consider more rapidly rotating models cases (Models B and C).

3.1. Mass loss

In Figure 1a, we show the local mass loss density,

Ṁ(r, θ, t) = 4πr^{2}\overline{ρ}(r, θ, t) \overline{U}_{r}(r, θ, t), (27)

whose average over θ and t, ⟨Ṁ⟩ = ∫_{0}^{π} ∫_{t_{0}}^{t_{0}+T} Ṁ dt sin θ dθ, is close to the initial value Ṁ_{0}. This is not too surprising, but it should be emphasized that this is not enforced as a condition.

The good agreement suggests that the open boundary condition at the bottom draws in a similar amount of mass at the inner boundary as what is lost at the outer boundary.

To get a sense of the radial mass distribution in our model, we plot in Figure 1b the cumulative mass,

M_{r}(r, θ, t) = ∫_{r}^{∞} 4πr^{′2}\overline{ρ}(r′, θ, t) dr′, (28)

for different values of θ at t = 858. We see that the total mass at r = r_{in} is about 7000 mass units; one mass unit here is Ṁ_{0}r_{c}/c_{s}. The mass above the surface is about 10, so 99.9% of the total mass in the computational domain is contained in the stellar envelope in r_{in} ≤ r ≤ R. Thus, if no mass was replenished on the inner boundary, the time it would take to lose all mass at the initial rate would be τ_{massloss} = M/Ṁ = 7000.

We emphasize at this point that the full stellar mass is undetermined, because the value of Newton’s constant G never enters on its own. We could, in principle, constrain it by assuming, for example, that the density in the core is constant and equal to that at r = r_{in}. This would give for the minimal core mass M_{core} ≫ 36000, which is five times the mass in the envelope. Using GM_{core} = 2, we find G ≪ 6 × 10^{−5}c_{s}^{3}/Ṁ_{0}, which is satisfied by a large margin for the values quoted above. We stress, however, that this estimate was done only for illustrative purposes.

3.2. Oscillatory model at solar rotation rate

We focus on a simulation with the solar value of the angular velocity, i.e., Ω̃ = 0.2 (Model A). In these cases, the magnetic field is oscillatory, but in a rather nonlinear fashion; see Figure 2, where we plot the time dependence of the three magnetic field components at one point in the wind. The \overline{B}_{r} component is positive most of the time and much smoother than the \overline{B}_{θ} and \overline{B}_{ϕ} components. The period T is about 41 time units. This corresponds to about 0.1 yr, which is short compared with the actual solar 22 year cycle, but still about five times longer than the cycle period in the model of Perri et al. (2020). Their parameters are otherwise comparable to ours: η_{T} ≈ 4 × 10^{14} cm^{2} s^{−1} in both models, Ṁ = 3 × 10^{−14}M_{⊙} yr^{−1} (a third of our value), an Alfvén radius of about two stellar radii, and a domain size of 20 solar radii (twice our value).

3.3. Magnetic field geometry

In Figure 3 we show a sequence of magnetic field visualizations at different times. To make the magnetic field in the outer parts better visible, we multiply B_{ϕ} by r^{2}. Here, we show the time span from t̃ = 814 to 858, covering just a little over a period. We overplot the surfaces where \overline{U}_{r} is transalfvénic (solid white lines), i.e., where \overline{U}_{r} exceeds the Alfvén speed v_{A} = (\overline{B}^{2}/μ_{0}\overline{ρ})^{1/2}. The surface is corrugated, but its mean radius is around 0.4 r_{c}. We also shows the surfaces where \overline{U}_{r} is transmagnetosonic, i.e., where \overline{U}_{r} exceeds the fast magnetosonic speed c_{ms} (dashed white line), which obeys c_{ms}^{2} = c_{s}^{2} + v_{A}^{2}. The mean radius of the magnetosonic surface is close to r_{c}.

Butterfly diagrams of \overline{B}_{r}(θ, t) and \overline{B}_{ϕ}(θ, t) are shown in Figure 4. The field in the wind does not show any migration in latitude, as is expected from models of the solar dynamo. Figure 5 shows only the inner part of the domain. We see regions with open and closed field lines at different times. However, there is no clear magnetic field migration that manifests itself in the Sun in a Maunder’s butterfly diagram of sunspot locations versus time and latitude.

It is interesting to note the appearance of V-shaped field lines in the panels for t = 822–834 and perhaps also for t = 846. This means that there are magnetic field lines in the wind that are not anchored in the star. This may be a bit surprising, but we have to remember that the magnetic field is time-dependent and the medium electrically conducting. The time-varying magnetic field can therefore induce toroidal currents in the stellar wind, which then produce poloidal field lines that are closed outside the star. This phenomenon may be similar to what is known as “switchbacks” in the solar wind (Bale et al., 2019; Squire et al., 2020).

3.4. Poynting flux

The wind carries with it not only mass, but also kinetic and magnetic energies. The latter is quantified by the mean Poynting flux,

F_{Poy}(r, t) = ∮ (\overline{E} × \overline{B}/μ_{0}) · dS, (29)

where \overline{E} = η_{T}μ_{0}\overline{J} − α\overline{B} − \overline{U} × \overline{B} is the mean electric field. The magnetic energy loss is then Ė_{M} = 4πr^{2}F_{Poy}. In the steady state, ⟨Ė_{M}⟩ would be independent of r if there was no Ohmic dissipation and no conversion between kinetic and magnetic energies in the wind.

As a good estimate for the magnetic energy loss of the solar wind, Brandenburg et al. (2011) computed Ė_{M}(r) ≈ 4πr^{2}⟨(\overline{B}^{2}/2μ_{0}) \overline{U}_{r}⟩, which they found to be of the order of 10^{18} W and slowly decreasing with radius. Estimating the total magnetic energy content within the convection zone based on a mean field of 300 G over the convection zone of volume 4π(R^{3} − r_{in}^{3})/3, we find a time scale of about 10 years, which is comparable with the solar cycle period.

Figure 6 shows the latitudinal dependence of Ė_{M} at different times for Model A. It depends not only on latitude and time, but also somewhat on radius. There is a window at high latitudes where it is almost constant in θ, but the width of this window changes with time. It can have a width of over 45^{∘} (e.g., at t = 818 and 858), but it can also be almost nonexistent (e.g., at t = 842). Comparing with Figure 3, we see that this window of nearly constant Ė_{M} corresponds to regions where the radial field in the wind ist mostly negative. The dips in Ė_{M} correspond to regions where the radial field is weak and changes sign. Near the equator, Ė_{M} shows a sharp drop for most times, except for t = 826. Again, comparing with Figure 3, we see that nothing special happens near those dips, except that for t = 826 the field is a bit weaker. These dips are probably a consequence of the radial field reversal in the equatorial plane and the existence of a field component that is purely vertical to the equatorial plane, thus inhibiting the wind.

Next, we look at the radial dependence of the kinetic and magnetic energy losses for different times and latitudes. The result is shown in Figure 7, where we define compute them as

Ė_{K} = 4πr^{2}(\overline{ρ}\overline{U}^{2}/2) u_{r}, (30)

Ė_{M} = 4πr^{2}(\overline{B}^{2}/2μ_{0}) u_{r}, (31)

respectively. It turns out that Ė_{M} is much smaller than Ė_{K}. To accommodate both quantities in the same plot, we have multiplied Ė_{M} by a factor of 20.

Fig. 3. Color representation of r^{2}B_{ϕ}(r, θ) for different times for Model A. The nearly concentric red solid lines show the surfaces where \overline{U}_{r} is transalfvénic and the red dashed ones show the surfaces where it is transmagnetosonic. The times correspond to the long tick marks of Figure 2.

Fig. 4. Butterfly diagrams of B_{r}(r, θ) and B_{ϕ}(r, θ) for Model A at r/r_{c} = 1.9. Note again the asymmetry of those components with respect to zero, which is different from the properties of the solar magnetic field.

We see that Ė_{K} increases with radius. This is a peculiar feature of isothermal models which is absent both in isentropic models with constant specific entropy and in nonisentropic models with variable specific entropy; see Figs. 9.18 and 9.20 of Brandenburg (2003), respectively. This is mainly because in those models the sound speed decreases with radius in such a way that the Mach number still increases, just as in the isothermal models. Thus, the basic dynamics is similar in that the flow becomes supersonic. In isothermal models, where the sound speed is constant, this transition must always be accompanied by a radial increase of the wind speed. In this sense, a polytropic model would seem more realistic, but it would still ignore the internal energy or entropy equation, which would be even more important for making our models more realistic, as is discussed below; see Sect. 3.8.

Fig. 5. Similar to Figure 3, but this time with a color representation of B_{ϕ}(r, θ) showing only the region close to the center. Note the occurrence of V-shaped field lines during certain times at 822 ≤ t ≤ 834, and 846. The field shows radial outward migration during certain times: negative B_{ϕ} at low latitudes for 814 ≤ t ≤ 826, and positive B_{ϕ} at midlatitudes for 834 ≤ t ≤ 854.

Fig. 6. Latitudinal dependence of the magnetic energy loss at different times for Model A. Note the occurrence of a plateau for small values of θ for 814 ≤ t ≤ 830 and after t = 854.

Fig. 7. Radial dependence of Ė_{K} and Ė_{M} for different latitude ranges at different times for Model A. Note that Ė_{M} has been multiplied by a factor of 20. Ė_{K} shows only little variability and always increases radially outward, while Ė_{M} has a maximum near the Alfvén surface at r/r_{c} ≈ 0.4. The maxima are particularly high for 822 ≤ t ≤ 826.

Fig. 8. Similar to Figure 7, for a semilogarithmic representation, without having rescaled Ė_{M}. Blue (red) lines indicate kinetic (magnetic) energy losses. Note that Ė_{M} ≈ Ė_{K} near the Alfvén surface at r/r_{c} ≈ 0.4, which is marked by a vertical line.

We also see that for t = 826, when the field was a bit weaker further out in the wind (Figure 3), the kinetic energy loss is particularly strong around the Alfvén surface; see Figure 7. At other times, especially for t ≤ 842 ≤ 850, the kinetic energy loss is generally much weaker. Comparing again with Figure 3, this corresponds to times when the radial field near the equator is strong.

In Figure 8, we show Ė_{M} and Ė_{K} for r ≤ r_{c} as a semilogarithmic representation. We see that Ė_{M} ≈ Ė_{K} at r/r_{c} ≈ 0.4. The radial profiles of Ė_{K} are fairly independent of θ and t. This is because the wind is rather powerful and not much affected by rotation or magnetic fields, which are the main factors that provide non-spherically symmetric contributions to the system.

It is interesting to note that Ė_{M}(r) has a maximum at r ≈ r_{c}/2. This radius is a certain distance above the stellar surface

Fig. 9. Latitudinal dependence of the angular momentum loss J̇(r, θ, t) for 1.5 ≤ r/r_{c} ≤ 2 at different times for Model A. The blue (red) lines refer to kinetic (magnetic) contributions, and the black lines denote the turbulent viscous contribution. Positive (negative) values are shown as solid (dotted) lines. Note the strong latitudinal variability of J̇.

and still below the critical point. This radius coincides with the Alfvén radius; see Figure 3. This is the point where most of the star’s magnetic energy has been deposited into the wind. In the Sun, we expect that this energy deposition occurs in the corona. One may tentatively associate the location of the maximum of Ė_{M}(r) with some representation of the star’s corona, although it is unclear whether there is any relation to the real corona of the Sun.

At large radii, r ≫ r_{c}, the magnetic energy loss declines slowly with radius. Such a decline has also been seen for the solar wind (Brandenburg et al., 2011). In the Sun, it may be connected with the conversion of magnetic energy into heat. In Figure 8, we see that Ė_{M} ≈ Ė_{K} at r/r_{c} ≈ 0.4.

3.5. Angular momentum flux

There are no sinks or sources to the angular momentum, to the angular momentum density, \overline{ρ}ϖ^{2}Ω, satisfies a conservation equation of the form (Mestel, 1968, 1999)

\frac{∂}{∂t} (\overline{ρ}ϖ^{2}Ω) = −∇ · F^{AM}, (32)

where

F^{AM} = \overline{ρ}ϖ^{2}Ω\overline{U} − ϖ\overline{B}_{ϕ}\overline{B}/μ_{0} − \overline{ρ}ν_{T}ϖ^{2}∇Ω (33)

is the angular momentum flux. Analogously to the energy loss, the expression for the angular momentum loss is J̇ = 4πr^{2}F_{r}^{AM}, which is shown in Figure 9 for 1.5 ≤ r/r_{c} ≤ 2 for the kinetic, magnetic, and viscous contributions, J̇_{K}, J̇_{M}, and J̇_{ν}, respectively. We see that the angular momentum flux is highly structured, with positive and negative contributions at different latitudes and times. At these radii, the kinetic term proportional to \overline{U}_{ϕ}\overline{U}_{r} dominates over the magnetic term proportional to \overline{B}_{ϕ}\overline{B}_{r}, and the turbulent viscous term is negligible.

The strongly negative contributions to the angular momentum flux are unexpected and may be connected with the time

Fig. 10. Time averaged radial profiles of latitudinally averaged J̇ for θ < 70^{∘} (upper panel) and θ > 70^{∘} (lower panel). The blue (red) lines refer to kinetic (magnetic) contributions, and the black lines denote the turbulent viscous contribution. Positive (negative) values are shown as solid (dotted) lines. The kinetic contributions from different times are shown in gray.

dependence of the solution. It may be of interest to study angular momentum fluxes along magnetic field lines; see the work of Pantolmos & Matt (2017), who compare flow speeds along different field lines. For our unsteady wind solutions, this procedure may no longer be particularly advantageous. However,

Fig. 11. Angular velocity contours superimposed on a color representation of \overline{U}_{r}(r, θ) for Model A. Positive (negative) values of Ω are shown as solid (dotted) black contours. Note the extended region at midlatitudes were Ω < 0.

to get some idea about the latitudes contributing to the negative angular momentum flux, we show in Figure 10 the radial dependence of the time- and latitude-averaged profiles of J̇ separately for the cones θ < 70^{∘} (away from the equator) and θ ≥ 70^{∘} (around the equator). We see that negative angular momentum fluxes dominate and originate mainly from regions away from the equator. Nevertheless, in the range 0.2 ≤ r/r_{c} ≤ 0.6, J̇_{K} and J̇_{M} can be of comparable magnitude, as is expected from the theory of Weber & Davis (1967). This range agrees well with the Alfvén radius; see Figure 3.

To understand the variability of Ω and the occurrence of negative values at certain times, we show in Figure 11 angular velocity contours superimposed on a color representation of \overline{U}_{r}(r, θ). Note that Ω is often negative over an extended range of midlatitudes. As we have seen above, this is chiefly responsible for the inward angular momentum transport discussed above. This could be related to our rather primitive modeling of the hydrodynamics inside the star, which lacks realistic differential rotation, for example. We return to this question briefly in the conclusions. Note also that \overline{U}_{r}(r, θ) shows clear latitudinal variations. The occurrence of regions with negative angular momentum transport is interesting in view of the recent discovery of fast wind episodes observed with Parker Solar Probe at certain longitudes (Finley et al., 2020a,b). Our model is of course axisymmetric and cannot address longitudinal variations, but it reminds us that negative angular momentum transport is not impossible.

We should point out that J̇ is given here in standard units where Ṁ = c_{s} = 1. Therefore, Figure 9 can be directly interpreted as a plot of the mean-field (MF) analogoue of the Shakura & Sunyaev (1973) parameter,

α_{SS}^{MF} = (\overline{ρ}\overline{U}_{r}\overline{U}_{ϕ} − \overline{B}_{r}\overline{B}_{ϕ}/μ_{0})/c_{s}^{2}. (34)

Here, the superscript MF indicates that this expression is applied to the two-dimensional mean fields rather than to the fluctuations, as in the usual turbulent case. This parameter is also frequently used in solar wind studies (see Eq. (2) of Finley et al., 2019); see also Keppens & Goedbloed (1999), Réville et al. (2015), and Pantolmos & Matt (2017) for earlier twodimensional stellar wind models.

The angular momentum in the dynamo zone is J_{∗} ≈ 68 in our units. Owing to cancelation, it is difficult to determine reliable values of J̇ and α_{SS}^{MF}, but for the purpose of a preliminary assessment, it suffices to estimate J̇ ≈ 0.01. As we discuss below in more detail, there can be certain periods where J̇ can even be negative. This then implies spindown or spinup at a rate τ_{spin} ≈ 7000, which is indeed similar to the value of τ_{massloss} quoted in Sect. 2.5. It may well be that α_{SS}^{MF} is much less than 0.01. This would then imply an even larger value of τ_{spindown}.

Fig. 12. Angular velocity contours superimposed on a color representation of \overline{U}_{r}(r, θ) for Model B (a) with Ω̃ = 1 and Model C (b) with Ω̃ = 10. The nearly concentric red solid lines show the surfaces where \overline{U}_{r} is transalfvénic and the red dashed ones show the surfaces where it is transmagnetosonic.

3.6. Resulting dynamo parameters

In our model, differential rotation is automatically established as a result of magnetic braking. Since our turbulent viscosity is assumed to be purely isotropic, differential rotation can only result from the torque on the star established by the magnetized wind (Mestel, 1968). This leads to a nearly constant angular momentum per unit mass, i.e., ϖ^{2}Ω ≈ const. The contours of constant angular velocity tends to approach a pattern that is close to cylindrical, as will be discussed below in the context of rapid rotation. Given that Ω ∝ ϖ^{−2}, the angular velocity difference between the r_{in} and R is ΔΩ = (1 − r_{in}^{2}/R^{2})Ω_{0} = 0.75 Ω_{0}. Therefore, we have for the second dynamo parameter in Eq. (20) the values C_{Ω} = 75, 375, and 3750 for Ω̃ = 0.2, 1, and 10, respectively. The first dynamo parameter in Eq. (20) is C_{α} = 125, where we have used α̃_{0} ≡ α_{0}/c_{s} = 0.05 for Model A, and η_{T} = 8 × 10^{−5}r_{c}c_{s}.

3.7. Rapid rotation

The study of models at rapid rotation is motivated by the interest in understanding the evolution of magnetic activity of young stars, i.e., before they have slowed down to the solar rotation rate. For us, there is also another motivation in that all our models were of α^{2} type, i.e., the Ω effect was weak and C_{Ω} was not much larger than C_{α}, as required for an αΩ dynamo (Brandenburg and Subramanian, 2005). To increase C_{Ω}, the rotation rate could be increased. Another possibility is to lower C_{α}. However, to prevent the dynamo from decaying, one would need to decrease η_{T} even further, but this is computationally difficult.

For rapid rotation, the magnetic field lines and contours of the toroidal magnetic field are much more concentrated to the bottom of the dynamo region, r ≈ r_{in}. At faster rotation, the contours become more cylindrical. This is an effect of the Taylor–Proudman theorem and results generally in small variations along the rotation axis.

Fig. 13. Magnetic field lines superimposed on a color representation of \overline{B}_{ϕ}(r, θ) for Model B with Ω̃ = 1. Strong fields only occur near r = r_{in}; the weak-field regions elsewhere cannot be seen.

The Taylor–Proudman theorem applies primarily the angular velocity contours. This can be seen by writing the relevant part of the \overline{U} · ∇\overline{U} nonlinearity of Eq. (2) in the form

ϕ̂ · ∇ × (−\overline{U} · ∇\overline{U})_{p} = ϖ \frac{∂}{∂z} Ω^{2} + ..., (35)

where Ω = \overline{U}_{ϕ}/ϖ is the local angular velocity, and the dots indicate the presence of other terms not relevant here. In Figure 12a we show contours of Ω together with a color-coded representation of \overline{U}_{r}. We see that the Ω contours are already strongly cylindrical for Ω̃ = 1. As we increase the value of Ω̃ to 10, the cylindrical contours begin to extent much further out along the rotation axis; see Figure 12b.

For Ω̃ = 10, the radial velocity develops a marked indentation inside of what is known as the inner tangent cylinder where

ϖ ≥ r_{in} (inner tangent cylinder); (36)

see Figure 12b. Here the outflow is suppressed and supersonic flows occur only for z ≥ 2r_{c} ≈ r_{out}, i.e., near the outer boundary of the computational domain. For Ω̃ = 1, by comparison, the contours of \overline{U}_{r}(r, θ) are almost perfect spherically symmetric – much more so than even for the case with Ω̃ = 0.2; cf. Figure 11. Similar results have also been found by Washimi & Shibata (1993) in their rotating models where a central dipole magnetic field was assumed.

It turns out that our models are now no longer oscillatory and are thus still not of αΩ type, contrary to what was originally hoped for. Visualizations of the toroidal and poloidal fields for Models B and C are shown in Figures 13 and 14. The fields are strong only inside the star, where the dynamo is active. Outside the star, the field is much weaker and not visible in our graphical representation, but it is never vanishing.

To discuss the nonoscillatory nature of these two models, it is useful to consider the dynamo parameters C_{α} and C_{Ω}. We find

Fig. 14. Similar to Figure 13, but for Model C with Ω̃ = 10. Stronger fields now extend along the axis outside the star.

C_{α} = 125 and C_{Ω} = 3750 for η_{T} = ν_{T} = 8×10^{−5}; see Table 1. To get an idea about the latitudinal variation of the magnetic field in the wind, the plot Ė_{M} as a function of θ for different radii. The result is shown in Figure 15. It turns out that the magnetic activity is confined to a narrow cone with an opening angle of about 15^{∘}.

Noticeable magnetic energy losses are found only near the rotation axis. As a function of radius, similarly to the case of slow rotation, Ė_{M}(r) has a maximum somewhere in R < r < r_{c}, which is where the Alfvén point lies. Furthermore, Model B has a much smaller magnetic energy loss at large radii than Model A.

Table 1. Summary of the simulations discussed in this paper.

Model  α̃ Q_{α}  Ω̃ C_{α} C_{Ω}  \overline{B}_{max} P_{cyc}
A  0.05 10^{−2} 0.2 125 75  6–13 41.0
B  0.1 10^{−2} 1 250 375 16.0  —
C  0.1 10^{−1} 10 250 3750 8.8  —

Fig. 15. Latitudinal dependence of Ė_{M} for different radius ranges for Model B. Note that Ė_{M} is large only near the axis.

Fig. 16. Radial profiles of the latitudinally averaged J̇ for Models A–C in panels (a)–(c). The blue (red) lines refer to kinetic (magnetic) contributions, and the black lines denote the turbulent viscous contribution. Positive (negative) values are shown as solid (dotted) lines. The total (kinetic, magnetic, and viscous) angular momentum transport is dominated by the kinetic contribution, except for Model A, where the magnetic contribution is rather strong, but negative in the outer parts.

The model shows similarities with earlier simulations of outflows emanating from stellar accretion disc dynamos (von Rekowski et al., 2003, 2004), but there the opening angle was closer to 30^{∘}. In the present simulations, the opening angle is essentially zero. It corresponds to a cylinder in which most of the magnetic fields are ejected, although the flow speed is here strongly reduced.

For these rapidly rotating models, we expect significant outward angular momentum transport. To demonstrate this in more detail, we show in Figure 16 the radial profiles of the latitu-

Fig. 17. Radial dependence of the plasma β for Models A (dotted black lines), B (dashed red line), and C (solid blue line).

dinally averaged J̇ for Models A–C for the kinetic, magnetic, and viscous contributions, just as we did in Figure 10. Since Models B and C are steady, time averaging is only needed for Model A.

Figure 16 shows that in Models B and C, the angular momentum transport is outward and J̇ is independent of r throughout most of the wind. For Model A, however, the time-averaged angular momentum transport becomes negative some distance away from the Alfvén point. Note also that J̇ is more than ten times larger in Model C than in the ten times more slowly rotating Model B. For Model B, the viscous contribution exceeds the magnetic one at all radii, while in Model C, the magnetic contribution exceeds the viscous one for r/r_{c} > 1.5. Inside the star, the angular momentum transport is negative and caused by a strong poleward circulation. The viscous contribution is also rather strong, but positive.

3.8. Comparison of the plasma betas for our models

We have seen that in Model A with the slowest rotation, the angular momentum flux was occasionally inward, especially at midlatitudes. We then considered Models B and C with faster rotation in the hope that not only the outward angular momentum flux would be outward, but also that the dynamo in the star would be in the αΩ regime. We found that the angular momentum flux was then indeed outward, but the dynamo was still in the α^{2} regime. In the introduction, we did already emphasize that the lack of a cool photosphere just beneath the corona was ignored. This makes it generally very difficult to reach low plasma betas, which we define as

β = 2ρc_{s}^{2}/B^{2}. (37)

In Figure 17, we plot the radial dependencies of the minimum value of β, β_{min}, for Models A–C. We see that the largest values of β_{min} occur for Model B with an intermediate angular velocity. Increasing the angular velocity further (Model C) increases the field strength and does therefore also lead to a smaller value of β_{min}. The smallest values occur for Model A. This is mainly because Model A is the only model where the magnetic field in the wind is of comparable strength at all latitudes. For faster rotation, the field in the wind is strongly concentrated around the axis.

Let us now return to the potential role of the photosphere. The photosphere of a star is the region where it cools and loses specific entropy. Everywhere else in the wind, the specific entropy does not change much, and therefore the potential enthalpy

Fig. 18. \overline{U}_{r} (red solid lines) and \overline{U}_{ϕ} (blue dashed lines) at the stellar surface at r = R as a function of colatitude for the same times as in Figure 10. \overline{U}_{r} is usually positive, but |\overline{U}_{ϕ}| is much larger and most of the time in the prograde direction, but sometimes it is retrograde, which is a consequence of the low low moment of inertia of the stellar envelope in our model.

must be approximately constant (von Rekowski et al., 2003). The potential enthalpy is defined as H = h + Φ, where h = c_{p}T is the specific enthalpy with c_{p} being the specific heat at constant pressure and T the temperature, and Φ = −GM/r is the potential energy. Hydrostatic equilibrium requires that

0 = −∇H + h∇s/c_{p}, (38)

where s is the specific entropy. For the corona, this implies T = GM/rc_{p} ≈ 2 × 10^{6} K, which is realistic and agrees also with our model. Toward the photosphere, T decreases abruptly because of surface cooling, and therefore the density increases abruptly. Thus, the density would then be much larger than what was possible in our models. This, in turn, would allow us to reach much larger field strengths and therefore smaller plasma betas. Another important consequence of having larger densities in the stellar envelope would be that the angular velocity at the stellar surface would always be in the prograde direction. In our present models, this is not always the case, as can be seen from Figure 18, where we show the radial and azimuthal velocities at the stellar surface. We see that the local rotational velocity is there occasionally in the retrograde direction, especially near the equator.

4. Conclusions

Our work has shown that a simplified realization of a dynamo with a stellar wind can easily be treated self-consistently in one and the same model, provided certain compromises are being made. The assumption of an isothermal equation of state has simplified matters conceptionally. Relaxing this restriction would allow us to include the energy deposition in the corona and to model the effects of a sharp density drop at the stellar surface. This might require a significant increase in resolution near the surface, which in turn requires the use of a nonuniform mesh. Another restriction has been the use of a relatively large turbulent magnetic diffusivity and viscosity. This was mainly needed to resolve shocks that develop within the wind. Those typically emerged in response to rapid changes in the magnetic field. This could probably be avoided by allowing for an additional shock viscosity, but this has been avoided in the present work. On the other hand, the angular momentum flux associated with turbulent viscosity was already negligible, so its presence may not have caused any artifacts.

Future work might involve the inclusion of a Λ effect (Rüdiger, 1980, 1989), which would allow for the development of differential rotation in the stellar envelope. Without including the effects of stellar winds, such models with combined α and Λ effects were studied by Brandenburg et al. (1990, 1991), who found significant alignment of the Ω contours with the rotation axis unless the baroclinic term was also included (Brandenburg et al., 1992). But this may change when their boundary condition on r = R is replaced by a continuous transition to the solar exterior; see Warnecke et al. (2013) for spherical convection simulations with a simplified representation of a stellar corona.

The inclusion of the Λ effect might allow us to model the stellar dynamo more realistically. It would be interesting to see how this affects the angular momentum transport and whether it could help in producing predominantly outward angular momentum transport in cases of slow rotation. It might then allow us to study dynamos in the αΩ regime. This has not been possible in the present model for reasons that we are not entirely clear, because the value of C_{Ω} was thought to be already large enough. There could have been other side effects arising from the coupling to the outflow that are not yet fully understood. Nevertheless, it is interesting to note that the inward angular momentum transport occurs even in the Sun within fast-wind regions at certain longitudes; see Finley et al. (2020a,b).

Another important aspect requiring further attention is the study of angular momentum losses from mean-field stresses. Our work has shown that the angular momentum loss can be quantified in terms of a nondimensional Shakura–Sunyaev parameter. This is a somewhat unusual concept in the context of stellar winds, but it may help putting the theories of turbulent stellar winds and accretion disk on a common footing.

Acknowledgements. We thank the referee for many useful remarks and suggestions that have significantly improved the manuscript. This work was supported in part through the Erasmus+ Programme of the European Union (P.J.) and the Swedish Research Council, grant 2019-04234 (A.B.). We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm.

Software and Data Availability. The source code used for the simulations of this study, the Pencil Code (Pencil Code Collaboration., 2020), is freely avail-

able on https://github.com/pencil-code/. The DOI of the code is

https://doi.org/10.5281/zenodo.2315093 (Brandenburg, 2018). The simulation setups and corresponding data are freely available on https://doi.org/ 10.5281/zenodo.4284439 (Jakab & Brandenburg, 2020).

References

Bale, S. D., Badman, S. T., Bonnell, J. W., Bowen, T. A., Burgess, D., Case, A. W., Cattell, C. A., Chandran, B. D. G., Chaston, C. C., Chen, C. H. K., Drake, J. F., de Wit, T. Dudok, Eastwood, J. P., Ergun, R. E., Farrell, W. M., Fong, C., Goetz, K., Goldstein, M., Goodrich, K. A., Harvey, P. R., et al. 2019, Nature, 576, 237

Biermann, L. 1951, Z. f. Ap., 29, 274 Brandenburg, A. 2003, in Advances in nonlinear dynamos (The Fluid Mechanics

of Astrophysics and Geophysics, Vol. 9), ed. A. Ferriz-Mas & M. Núñez (Taylor & Francis, London and New York), 269

Brandenburg, A., on behalf of the Pencil Code Collaboration, 2018, Pencil

Code, v2018.12.16, Zenodo, DOI:10.5281/zenodo.2315093 Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 Brandenburg, A., Moss, D., Rüdiger, G., & Tuominen, I. 1990, Solar Phys., 128,

243251 Brandenburg, A., Moss, D., Rüdiger, G., Tuominen, I. 1991, Geophys.

Astrophys. Fluid Dyn., 61, 179 Brandenburg, A., Moss, D., & Tuominen, I. 1992, A&A, 265, 328 Brandenburg, A., Subramanian, K., Balogh, A., & Goldstein, M. L. 2011, ApJ,

734, 9 Chamberlain, J. W. 1960, ApJ, 131, 47 Cole, E., Brandenburg, A., Käpylä, P. J., & Käpylä, M. J. 2016, A&A, 593, A134 Del Sordo, F., Guerrero, G., & Brandenburg, A. 2013, MNRAS, 429, 1686 Finley, A. J., Hewitt, A. L., Matt, S. P., Owens, M., Pinto, R. F., & Réville, V.

2019, ApJ, 885, L30 Finley, A. J., Matt, S. P., Réville, V., Pinto, R. F., Owens, M., Kasper, J. C.,

Korreck, K. E., Case, A. W., Stevens, M. L., Whittlesey, P., Larson, D., & Livi, R. 2020a, ApJ, 902, L4

Finley, A. J., McManus, M. D., Matt, S. P., Kasper, J. C., Korreck, K. E., Case, A. W., Stevens, M. L., Whittlesey, P., Larson, D., Livi, R., Bale, S. D., Dudok de Wit, T., Goetz, K., Harvey, P. R., MacDowall, R. J., Malaspina, D. M., & Pulupa, M. 2020b, A&A, in press, arXiv:2011.00016, https://doi.org/10.1051/0004-6361/202039288

Frank, J., King, A. R., & Raine, D. J. 1992, Accretion power in astrophysics

(Cambridge: Cambridge Univ. Press) Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS,

438, L76 Jabbari, S., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2015,

ApJ, 805, 166 Jakab, J. & Brandenburg, A., 2020, Datasets for “The effect of a

dynamo-generated field on the Parker wind,” v2020.11.22, Zenodo, DOI:10.5281/zenodo.4284439

Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 Keppens, R., & Goedbloed, J. P. 1999, A&A, 343, 251 Krause, F., & Rädler, K.-H. 1980, Mean-field Magnetohydrodynamics and

Dynamo Theory (Oxford: Pergamon Press) Malkus, W. V. R., & Proctor, M. R. E. 1975, J. Fluid Mech., 67, 417 Mestel, L. 1968, MNRAS, 138, 359 Mestel, L. 1999, Stellar Magnetism (Clarendon Press, Oxford) Metcalfe, T. S., & van Saders, J. 2017, Solar Phys., 292, 126 Mitra, D., Moss, D., Tavakol, R., & Brandenburg, A. 2011, A&A, 526, A138 Parker, E. N. 1958, ApJ, 128, 664 Pantolmos, G., & Matt, S. P. 2017, ApJ, 849, 83 Pencil Code Collaboration: Brandenburg, A., Johansen, A., Bourdin, P. A.,

Dobler, W., Lyra, W., Rheinhardt, M., Bingert, S., Haugen, N. E. L., Mee, A., Gent, F., Babkovskaia, N., Yang, C.-C., Heinemann, T., Dintrans, B., Mitra, D., Candelaresi, S., Warnecke, J., Käpylä, P. J., Schreiber, A., Chatterjee, P., Käpylä, M. J., Li, X.-Y., Krüger, J., Aarnes, J. R., Sarson, G. R., Oishi, J. S., Schober, J., Plasson, R., Sandin, C., Karchniwy, E., Rodrigues, L. F. S., Hubbard, A., Guerrero, G., Snodin, A., Losada, I. R., Pekkilä, J., & Qian, C. 2020, Journal of Open Source Software, 5, 2684, arXiv:2009.08231

Perri, B., Brun, A. S., Réville, V., & Strugarek, A. 2018, J. Plasma Phys., 84,

765840501 Perri, B., Brun, A. S., Réville, V., & Strugarek, A. 2020, in 223, ed. Magnetic

field evolution in solar-type stars (IAUS 354: Solar and Stellar Magnetic Fields: Origins and Manifestations), 215 A. Kosovichev, K. Strassmeier & M. Jardine Proc. IAU Symp., Vol. 354

Pinto, R. F., Brun, A. S., Jouve, L., & Grappin, R. 2011, ApJ, 737, 72 Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798,

116 Roberts, P. H. 1972, Phil. Trans. R. Soc., A272, 663 Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dyn., 16, 239 Rüdiger, G. 1989, Differential rotation and stellar convection: Sun and solar-type

stars (Gordon & Breach, New York) See, V., Matt, S. P., Finley, A. J., Folsom, C. P., Boro Saikia, S., Donati, J.-

F., Fares, R., Hébrard, É. M., Jardine, M. M., Jeffers, S. V., Marsden, S. C., Mengel, M. W., Morin, J., Petit, P., Vidotto, A. A., Waite, I. A., the BCool Collaboration 2019, ApJ, 886, 120

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 Shore, S. N. 1992, An introduction to astrophysical hydrodynamics (San Diego:

Academic Press) Stix, M. 1974, A&A, 37, 121 Squire, J., Chandran, B. D. G., & Meyrand, R. 2020, ApJ, 891, 2 van Saders, J. L., Ceillier, T., Metcalfe, T. S., Silva Aguirre, V., Pinsonneault, M.

H., García, R. A., Mathur, S., & Davies, G. R. 2016, Nature, 529, 181 von Rekowski, B., Brandenburg, A., Dobler, W., & Shukurov, A. 2003, A&A,

398, 825 von Rekowski, B., & Brandenburg, A. 2004, A&A, 420, 17 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013, ApJ, 778,

141 Washimi, H., & Shibata, S. 1993, MNRAS, 262, 936 Weber, E. J., & Davis, L., Jr. 1967, ApJ, 148, 217