The effect of a dynamo-generated field on the Parker wind
aa r X i v : . [ a s t r o - ph . S R ] J un Astronomy&Astrophysicsmanuscript no. paper c (cid:13)
ESO 2020June 5, 2020
The effect of a dynamo-generated field on the Parker wind
P. Jakab , and A. Brandenburg , , , Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Nuclear Physics and Subnuclear Physics, Pavol Jozef Safarik University in Kosice, Jesenna 5, 041 54 Kosice,Slovakia Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden JILA and Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80303, USA McWilliams Center for Cosmology & Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USAJune 5, 2020, Revision: 1.69
ABSTRACT
Context.
Stellar winds are an integral part of the underlying dynamo, the motor of stellar activity. The wind controls the star’s angularmomentum 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 dynamoin 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 usingan isothermal equation of state. The model allows for the simultaneous generation of a mean magnetic field and the development ofa Parker wind. The resulting flow is transonic at the critical point, which we arrange to be between the inner and outer radii of themodel. 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 α type. In most of the domain, the magnetic field correspondsto that of a split monopole. The magnetic energy flux is largest between the stellar surface and the critical point. At rapid rotation ofup 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 reproducethe solar wind: highly variable angular momentum flux fluxes even from just an α dynamo in the star. For rapid rotation, magneticfields are ejected mostly along the axis, where the wind speed is reduced. Key words.
Sun: sunspots – Sun: dynamo – turbulence – magnetohydrodynamics (MHD) – hydrodynamics
1. Introduction
The emergence of a wind around stars is a remarkable and some-what counter-intuitive phenomenon. The existence of the solarwind was already suggested by the fact that the tails of comets al-ways point away from the Sun (Biermann, 1951). Nevertheless,the wind was thought to be a relatively slow phenomenon asso-ciated with an evaporation of the corona (Chamberlain, 1960).The physical nature and mathematical theory of the solar windwas first understood by Parker (1958). His theory showed thatthe wind starts off as a subsonic flow some distance above thecorona. It gradually gains in speed as the gravitational forcediminishes and the effective outward pull resulting from thequadratic increase of the cross-sectional area in Bernoulli’s law.This is a purely hydrodynamic phenomenon, unlike what wassuggested by the popular notion of the solar corpuscular radia-tion at the time.Stellar winds play a crucial role in a star’s life. Without thewind, the Sun would still be spinning rapidly and magneticallysuperactive. A proper understanding of the rotational evolutionof a star through magnetic breaking via a wind is importantnot only for stellar evolution, but it also plays a role in under-standing the diversity of magnetic activity as a function of rota-tion rate and age (van Saders et al., 2016). As the star reachesthe age of the Sun, the magnetic field either changes its ge-ometry such that stellar breaking is reduced (See et al., 2019;Metcalfe & van Saders, 2017) or it can continue to break and thestar’s differential rotation becomes antisolar-like (Gastine et al., 2014; K¨apyl¨a et al., 2014), i.e., the equator spins slower than thepoles. Stellar winds can also be important for the dynamo it-self in that they can transport magnetic helicity away from thedynamo region, and thereby alleviate what is known as catas-trophic quenching; see Mitra et al. (2011) for mean-field mod-els and Del Sordo et al. (2013) for computations of the magnetichelicity flux in simulations in a turbulent wind. Magnetic windsalso affect the density and dynamics of cosmic rays in the helio-sphere. Computing selfconsistently the dynamo-generated mag-netic field evolution in the heliosphere is therefore crucial andfor modeling the magnetic shielding of Galactic cosmic rays onthe Earth.The theory of a magnetized stellar wind by Weber & Davis(1967) employes a prescribed and time-independent stellar mag-netic field, so any feedback on the underlying dynamics wasignored. This is also true of the recent numerical models ofR´eville et al. (2015), who compared different magnetic multi-poles as initial conditions of their models. This has changed onlyin recent years. Given that the wind normally dominates over themagnetic field, one can separate the dynamics of the wind fromthat of the solar dynamo. In recent work of Perri et al. (2018),this was modeled using two separate codes that are magneticallycoupled through a matching condition at the solar surface.The purpose of the present paper is to explore some basicproperties of stellar winds in the presence of dynamo-generatedmagnetic fields. It is appropriate to adopt a mean-field model,where we solve the equations for the azimuthally averaged mag-
1. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind netic and velocity fields. In this paper, those mean fields are de-noted by an overbar. The effects of turbulence are then param-eterized through a turbulent viscosity and a turbulent magneticdiffusivity. In the star’s convection zone, there are also cyclonicconvective motions giving rise to kinetic helicity of oppositesigns in the two hemispheres. This is modeled by an α effect(Krause & R¨adler, 1980). The turbulent magnetic diffusivity ishere assumed constant.The presence of the magnetic field causes the kinetic andmagnetic stresses to be different from zero. The turbulent viscos-ity is itself a result of kinetic and magnetic stresses caused by thefluctuating components of the magnetic and velocity fields. Inthe theory of turbulent accretion discs (Frank et al., 1992), thosestresses are parameterized by the Shakura & Sunyaev (1973) pa-rameter, α SS . It quantifies the stress in terms of the backgrounddifferential rotation, the sound speed, and the scale height. Inaccretion discs, where the differential rotation is Keplerian, thisamounts to a scaling of the stress by the sound speed squared.In our case, the differential rotation is not related to the soundspeed, but the basic mechanism of angular momentum transferis the same, and we can still express the total stress in a similarfashion.Unlike the work of Perri et al. (2018), we consider the evo-lution of the dynamo and the wind with a single code. At thispoint, our aim is not to produce a realistic model of the Sun, butrather a physically consistent model under conditions where thedynamics of the wind can no longer be separated from that of thedynamo. Our models can also be applied to conditions of rapidrotation, which strongly affects the wind. This can be particu-larly relevant to young stars in their T Tauri phase. We begin bypresenting the basic equations of our model and turn then to thediscussion of our results.The simplest wind solution is the isothermal one that was al-ready found by Parker (1958). Heating is not explicitly invoked.Its physics resembles that of a siphon flow. Once a fluid parcelhas moved over the top of the effective gravitational potential, itsimply continues to fall and pulls the remaining fluid behind it.The top of the effective potential corresponds to the critical pointwhere the flow speed crosses the sonic point (see Appendix A foran illustration). We arrange this point to be in the middle of thecomputational domain such that the flow speed becomes super-sonic well before the outer point r out . We fit the dynamo-activezone (or stellar envelope) with an α effect different from zerointo a spherical shell between the inner point of the computa-tional domain, r in , and a radius R , which models the surface ofthe star.The usefulness of an isothermal solution can be justified byconsidering the fact that the sound speed both at the bottom ofthe convection zone and in the solar wind is about
100 km s − ,corresponding to a temperature of a million degrees. The lowertemperature near the photosphere is obviously ignored. For anisothermal gas, the mean pressure p is then simply proportionalto the gas density ρ with p = ρc , where c s is the sound speed.The pressure gradient is then given by ( ∇ p ) /ρ = c ∇ ln ρ .We begin by discussing first the basic equations, boundaryconditions, and parameters in Sect. 2. We then present our resultsin Sect. 3, and draw our conclusions in Sect. 4.
2. The model
We adopt spherical polar coordinates, ( r, θ, φ ) , with the origin atthe center of the star. The vector r points away from the center,the colatitude θ increases away from the North pole, and φ in- creases in the eastward direction. We assume axisymmetry, i.e., ∂/∂φ = 0 . We write the mean magnetic field as B = ∇ × A , where A isthe mean vector potential. This ensures that ∇ · B = 0 at alltimes. The evolution equations for A , the mean velocity U , andthe logarithmic mean density ρ , are ∂ A ∂t = U × B + α B − η T J , (1) D U D t = − c ∇ ln ρ − GMr ˆ r + 1 ρ J × B − ν T Q , (2) D ln ρ D t = − ∇ · U , (3)where D / D t = ∂/∂t + U · ∇ is the advective derivative, c s isthe isothermal sound speed, G is Newton’s constant, M is thestellar mass, ˆ r = r /r is the radial unit vector, η T and ν T are thesums of turbulent and microphysical values of magnetic diffusiv-ity and kinematic viscosity, respectively, α is the aforementionedcoefficient in the α effect, J = ∇ × B /µ is the mean currentdensity, µ is the vacuum permeability, − Q = ∇ U + ∇∇ · U + 2 S · ∇ ln( ν T ρ ) (4)is a term appearing in the viscous force, where S is the tracelessrate of strain tensor of the mean flow with components S ij = ( U i,j + U j,i ) − δ ij ∇ · U . The dot in Eq. (4) denotes thecontraction over the free index of ∇ ln( ν T ρ ) .The mean magnetic field is generated by the α effect. Thisleads to exponential growth, provided the value of α is abovea certain critical value. Eventually, the dynamo must saturatebecause the Lorentz force from the mean field, J × B , drivesfluid motions that feed back onto the dynamo to limit its growth.This way of achieving saturation is sometimes referred to asMalkus & Proctor (1975) mechanism. In addition, there can befeedback from the small-scale magnetic field that leads to a non-linear suppression of α , which is referred to as α quenching. Weassume here a simple quenching function for α , which is thenwritten in the form α ( r, θ, B ) = α f α ( r ) cos θ sin n θ Q α B /B , (5)where n = 6 is chosen to concentrate the α effect to low lati-tudes (Jabbari et al., 2015; Cole et al., 2016), Q α is a quenchingparameter that determines the typical field strength, which is ex-pected to be of the order of Q − / α B eq , and f α ( r ) = Θ (cid:16) ( r − R ) /w α (cid:17) (6)is a radial profile function with Θ( x ) being a smoothened stepfunction from 0 to 1 as x crosses zero. Here, R and w α determinethe location and width of the transition. The value of Q α de-termines the nonlinear equilibration of the dynamo, in additionto the macroscopic feedback from the Lorentz force mentionedabove. Our model thus comprises three distinct layers with r in < R < r ∗ < r out , (7)where r in < r < R is the dynamo region (modeling the stellarenvelope), R < r < r ∗ is the wind acceleration region (model-ing the locations of the solar corona and the Alfv´en point), and r ∗ < r < r out is the supersonic wind region with r ∗ = GM/ c being the critical point.
2. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
In most of the cases, we apply a uniform angular velocity Ω onthe inner boundary r = r in by setting u φ = r in sin θ Ω . Forthe other two velocity components, we adopt “open” boundaryconditions by setting the second radial derivative to zero. Thiscondition turns out to be stable in all cases considered in thispaper. It allows for a weak inflow to replenish the mass loss onthe outer boundary r = r out , where we apply open boundaryconditions for all three velocity components. No precautions aretaken to ensure that the mass in the computational domain staysconstant. It turns out, however, that the total mass remains nearlyunchanged. This is, to some extent, also explained by the factthat the total mass loss rate is small compared with other inversetime scales in the problem.For the magnetic field, we adopt a perfect conductor bound-ary conditions on the inner radius, i.e., ∂A r ∂r = A θ = A φ = 0 on r = r in , (8)and a radial field condition on the outer radius, i.e., A r = ∂A θ ∂r + A θ r = ∂A φ ∂r + A φ r = 0 on r = r out . (9)On the pole, we assume ∂A r ∂θ = A θ = A φ = 0 on θ = 0 ◦ , (10)while on the equator, we assume ∂A r ∂θ = A θ = ∂A φ ∂θ = 0 on θ = 90 ◦ . (11)Since our simulations are axisymmetric, the magnetic field isconveniently represented via B φ and A φ . In particular, contoursof r sin θ A φ give the magnetic field lines of the poloidal field, B pol = ∇ × ( A φ ˆ φ ) . As initial condition for U ≡ ( u, , and ρ , we adopt the Parkerwind solution. In some cases we also add a finite angular ve-locity with constant angular momentum, although its effect onthe dynamics is ignored in the initial condition. We begin bydiscussing the Parker wind solution, which can be obtained bysolving the Bernoulli equation, u + c ln ρ − GM/r = const , (12)along with the equation of mass conservation, which states thatthe mass loss rate is given by ˙ M = 4 πr ρu . We then obtain u − c ln u − c ln r − GM/r = Φ , (13)where Φ = − / is obtained by inserting the values u = r ∗ =1 for the critical point. We solve the Bernoulli equation itera-tively. For r ≤ r ∗ , using u = c s r/r ∗ initially, we iterate c ln u i +1 ( r ) = u i − c ln r − GM/r − Φ , (14)while for r > r ∗ , using u = 2 c s initially, we iterate u i +1 ( r ) = c ln u i + c ln r + GM/r + Φ . (15)This iteration procedure was implemented by J¨orn Warnecke andDhrubaditya Mitra into the P ENCIL C ODE in 2012. We choosethe initial value of ˙ M to be ˙ M . http://github.com/pencil-code ,DOI:10.5281/zenodo.2315093 It is convenient to work with nondimensional units by measur-ing speeds in units of the isothermal sound speed and lengthsin units of the critical radius, r ∗ = GM/ c . In the following,we use tildae to denote nondimensional quantities. Using typicalnumbers for the Sun, we have c s = 10 cm s − = 100 km s − , (16) GM = GM ⊙ ≈ . × cm s − , and therefore r ∗ = GM ⊙ / c ≈ × cm ≈ R ⊙ ≈ .
05 AU . (17)In the Sun, the turbulent viscosity is ν T ≈ u rms ℓ/ ≈ cm s − . The nondimensional viscosity is then ˜ ν T ≡ ν T c s GM ⊙ ≈ × − , (18)which is rather small.For numerical stability, as already alluded to, we cannotchoose the value of ν T to be too small. In practice, for a numeri-cal resolution of × mesh points in the r and θ directions,we can choose ˜ ν T ≈ . . For × mesh points, on theother hand, we can reduce it by a factor of 128 to ˜ ν T ≈ × − .This then also means that in the stellar convection zone, we can-not adopt significantly smaller values, as is expected theoreti-cally based on our earlier estimates of u rms and ℓ .The nondimensional value of the angular velocity is given by ˜Ω = r ∗ Ω /c s = GM Ω / c ≈ . , (19)where we have used Ω = 3 × − s − . The strength of thedynamo is determined by the two dynamo numbers, C α = α R/η T and C Ω = ∆Ω R /η T . (20)The excitation conditions for dipolar and quadrupolar paritiesare generally fairly close together (Roberts, 1972). This is be-cause the magnetic field is strongest at high latitudes, so thehemispheric coupling is weak. In the following we restrict our-selves 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 s = ˙ M = µ = 1 , (21)which implies that GM = 2 . Our unit of mass is then [ M ] =˙ M r ∗ /c s . For the Sun, we have ˙ M ≈ × g s − , sothat our unit of density is [ ρ ] = ˙ M /c s r ∗ , which is about . × − g cm − for the Sun. Therefore, our unit of B is [ B ] = ( µ [ ρ ]) / c s , which is about .
04 G for the Sun. The valueof Newton’s constant G never enters on its own. It could be de-termined a posteriori, if we knew the total stellar mass. In ourmodel, we can compute the mass M ∗ of the stellar envelope in r in ≤ r ≤ R , but this still leaves the mass of the stellar coreundetermined. In the following, it is often convenient to retainthe symbols r ∗ , c s , ˙ M , and µ to remind ourselves of the nor-malization.There are a few other parameters of the model that we keepfixed. In all cases we use w α = 0 . for the transition thicknessof α near the surface; see Eq. (6). We always take r in = 0 . and R = 0 . . This corresponds to a fractional shell thickness of 50%instead of the 30% in the case of the Sun, but we should keep inmind that there are other properties that agree with the Sun onlyqualitatively. Another example is our smaller choice of R/r in =5 instead of the solar value of about 10. In all our simulationswith × meshpoints, we use ˜ ν T = 8 × − r ∗ c s .
3. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 1.
Radial dependence of ˙ M (a) and M r (b) for Model A. In our simulations, sound speed and the critical radius are set tounity, so the characteristic sound travel time, τ s = r ∗ /c s (22)is therefore also unity. When we adopt the stellar rotation rate, ˜Ω = 0 . , the corresponding rotational time scale τ Ω = Ω − (23)is then five, and the rotation period is π/ ˜Ω ≈ . The charac-teristic time scale for the dynamo is the turbulent diffusive time(e.g., Stix, 1974), τ TD = R /η T , (24)which is around in our models. Another interesting timescale for our models is the mass loss time, τ massloss = M/ h ˙ M i ≈ M / ˙ M . (25)In our models, M ≈ and ˙ M = 1 , so τ massloss ≈ . Itturns out that the spindown time is of a similar order of magni-tude. It is given by τ spindown = J ∗ / h ˙ J i , (26)where J ∗ = R ∗ ρ̟ Ω dV is the angular momentum of the stel-lar envelope, and ˙ J is the angular momentum loss. The asteriskon the integral denotes the volume of the envelope. The massloss and spindown times are the longest among the time scalesconsidered here, so the mass in the envelope cannot change sig-nificantly during the time scales of interest for the wind and thedynamo. Fig. 2.
Time series of the three magnetic field components atone point for Model A.
3. Results
After some preliminary studies at low resolution of × meshpoints with ν T = η T = 10 − r ∗ c s , we performed high-resolution simulations with × meshpoints, where wewere able to decrease ν T and η T to × − r ∗ c s . These valuesare still above the physically motivated value, but for numeri-cal stability reasons, they cannot be decreased further withoutinvoking artificial viscosity and magnetic diffusivity.Our main model is called Model A, which has the solar valueof Ω and a minimal amount of viscosity and magnetic diffusivitythat can still be tolerated. Later, we also consider more rapidlyrotating models cases. In Figure 1(a), we show the local mass loss density, ˙ M ( r, θ, t ) = 4 πr ρ ( r, θ, t ) U r ( r, θ, t ) , (27)whose average over θ and t , h ˙ M i = R π R t + Tt ˙ M dt sin θ dθ ,is close to the initial value ˙ M . This is not too surprising, butit should be emphasized that this is not enforced as a condition.The good agreement suggests that the open boundary conditionat the bottom draws in a similar amount of mass at the innerboundary as what is lost at the outer boundary.To get a sense of the radial mass distribution in our model,we plot in Figure 1(b) the cumulative mass, M r ( r, θ, t ) = Z ∞ r πr ′ ρ ( 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 ˙ M r ∗ /c s .The mass above the surface is about 10, so 99.9% of the totalmass in the computational domain is contained in the stellar en-velope in r in ≤ r ≤ R . Thus, if no mass was replenished onthe inner boundary, the time it would take to lose all mass at theinitial rate would be τ massloss = M/ ˙ M = 7000 .We emphasize at this point that the full stellar mass is unde-termined, because the value of Newton’s constant G never en-ters on its own. We could, in principle, constrain it by assuming,for example, that the density in the core is constant and equalto that at r = r in . This would give for the minimal core mass M core ≫ , which is five times the mass in the envelope.
4. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 3.
Color representation of r B φ ( r, θ ) for different times for Model A. The nearly concentric white solid lines show the surfaceswhere U r is transonic and the dashed ones show the surfaces where it is transmagnetosonic. Fig. 4.
Butterfly diagrams of B φ ( r, θ ) and B r ( r, θ ) for Model A at r/r ∗ = 1 . .Using GM core = 2 , we find G ≪ × − c / ˙ M , which is sat-isfied by a large margin for the values quoted above. We stress,however, that this estimate was done only for illustrative pur-poses. We focus on a simulation with the solar value of the angularvelocity, i.e., ˜Ω = 0 . (Model A). In these cases, the magneticfield is oscillatory, but in a rather nonlinear fashion; see Figure 2,where we plot the time dependence of the three magnetic fieldcomponents at one point in the wind. The B r component is pos-itive most of the time and much smoother than the B θ and B φ components. The period T is about 41 time units.
5. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 5.
Similar to Figure 3, but this time with a color representation of B φ ( r, θ ) showing only the region close to the center. In Figure 3 we show a sequence of magnetic field visualizationsat different times. To make the magnetic field in the outer partsbetter visible, we multiply B φ by r . Here, we show the timespan from ˜ t = 814 to , covering just a little over a period. Weoverplot the surfaces where U r is transonic (solid green lines),i.e., where U r exceeds the Alfv´en speed v A = ( B /µ ρ ) / .The surface is corrugated, but its mean radius is around . r ∗ .We also shows the surfaces where U r is transmagnetosonic, i.e.,where U r exceeds the fast magnetosonic speed c ms (dash greenline), which obeys c = c + v . The mean radius of the mag-netosonic surface is close to r ∗ .Butterfly diagrams of B r ( θ, t ) and B φ ( θ, t ) are shown inFigure 4. The field in the wind does not show any migrationin 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 differenttimes. However, there is no clear magnetic field migration thatmanifests itself in the Sun in a Maunder’s butterfly diagram ofsunspot locations versus time and latitude. The wind carries with it not only mass, but also kinetic and mag-netic energies. The latter is quantified by the mean Poynting flux, F Poy ( r, t ) = I ( E × B /µ ) · d S , (29)where E = µ η T J − α B − U × B is the mean electric field.The magnetic energy loss is then ˙ E M = 4 πr F Poy . In the steadystate, h ˙ E M i would be independent of r if there was no Ohmicdissipation and no conversion between kinetic and magnetic en-ergies in the wind.As a good estimate for the magnetic energy loss of thesolar wind, Brandenburg et al. (2011) computed ˙ E M ( r ) ≈ πr h ( B / µ ) U r i , which they found to be of the order of W and slowly decreasing with radius. Estimating the to-tal magnetic energy content within the convection zone basedon a mean field of
300 G over the convection zone of volume π ( R − r ) / , we find a time scale of about 10 years, whichis comparable with the solar cycle period.Figure 6 shows the latitudinal dependence of ˙ E M at differenttimes for Model A. It depends not only on latitude and time, butalso somewhat on radius. There is a window at high latitudeswhere it is almost constant in θ , but the width of this windowchanges with time. It can have a width of over ◦ (e.g., at t =
6. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 6.
Latitudinal dependence of the magnetic energy loss at different times for Model A.
Fig. 7.
Radial dependence of the magnetic and kinetic energy losses at different times for Model A. Note that ˙ E M has beenmultiplied by a factor of 20.
7. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 8.
Similar to Figure 7, for a semilogarithmic representation, without having rescaled ˙ E M . Blue (red) lines indicate kinetic(magnetic) energy losses. Note that ˙ E M ≈ ˙ E K at r/r ∗ ≈ . . Fig. 9.
Latitudinal dependence of the angular momentum loss ˙ J ( r, θ, t ) at different times for Model A. The blue (red) lines refer tokinetic (magnetic) contributions, and the black lines denote the turbulent viscous contribution. Positive (negative) values are shownas solid (dotted) lines. and 858), but it can also be almost nonexistent (e.g., at t =842 ). The kinetic energy in the wind increases with radius at alltimes. It must get its energy from the thermal energy, but if thisis not included in our isothermal model.
8. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Next, we look at the radial dependence of the kinetic andmagnetic energy losses for different times and latitudes. The re-sult is shown in Figure 7, where we define compute them as ˙ E K = 4 πr ( ρ U / u r , (30) ˙ E M = 4 πr ( B / µ ) u r , (31)respectively. It turns out that ˙ E M is much smaller than ˙ E K . Toaccommodate both quantities in the same plot, we have multi-plied ˙ E M by a factor of 20.We see that ˙ E K increases with radius. The radial profiles of ˙ E K are fairly independent of θ and t . This is because the windis rather powerful and not much affected by rotation or magneticfields, which are the main factors that provide non-sphericallysymmetric contributions to the system.It is interesting to note that ˙ E M ( r ) has a maximum at r ≈ r ∗ / . This radius is a certain distance above the stellar surfaceand still below the critical point. This radius coincides with theAlfv´en radius; see Figure 3. This is the point where most of thestar’s magnetic energy has been deposited into the wind. In theSun, we expect that this energy deposition occurs in the corona.One may tentatively associate the location of the maximum of ˙ E M ( r ) with some representation of the star’s corona, although itis unclear whether there is any relation to the real corona of theSun.At large radii, r ≫ r ∗ , the magnetic energy loss declinesslowly with radius. Such a decline has also been seen for thesolar wind (Brandenburg et al., 2011). In the Sun, it may beconnected with the conversion of magnetic energy into heat. InFigure 8, we show ˙ E M and ˙ E K for r ≤ r ∗ as a semilogarithmicrepresentation. We see that ˙ E M ≈ ˙ E K at r/r ∗ ≈ . . There are no sinks or sources to the angular momentum, tothe angular momentum density, ρ̟ Ω , satisfies a conservationequation of the form (Mestel, 1968, 1999) ∂∂t (cid:0) ρ̟ Ω (cid:1) = − ∇ · F AM , (32)where F AM = ρ̟ Ω U − ̟B φ B /µ − ρν T ̟ ∇ Ω (33)is the angular momentum flux. Analogously to the energyloss, the expression for the angular momentum loss is ˙ J =4 πr F AM r , which is shown in Figure 9. We see that the angularmomentum flux is highly structured, with positive and negativecontributions at different latitudes and times. The kinetic termproportional to U φ U r dominates over the magnetic term pro-portional to B φ B r , and the turbulent viscous term is negligible.We should point out that ˙ J is given here in standard unitswhere ˙ M = c s = 1 . Therefore, Figure 9 can be directly inter-preted as a plot of the Shakura & Sunyaev (1973) parameter, α MFSS = ( ρU r U φ − B r Bφ/µ ) /c . (34)Here, the superscript MF indicates that this expression is appliedto the two-dimensional mean fields rather than to the fluctua-tions, as in the usual turbulent case. The angular momentum in the dynamo zone is J ∗ ≈ inour units. Owing to cancelation, it is difficult to determine re-liable values of ˙ J and α MFSS , but for the purpose of a prelimi-nary assessment, it suffices to estimate ˙ J ≈ . . This then im-plies τ spindown = 6800 , which is indeed similar to the value of τ massloss quoted in Sect. 2.5. It may well be that α MFSS is muchless than 0.01. This would then imply an even larger value of τ spindown . In our model, differential rotation is automatically establishedas a result of magnetic breaking. Since our turbulent viscosityis assumed to be purely isotropic, differential rotation can onlyresult from the torque on the star established by the magnetizedwind (Mestel, 1968). This leads to a nearly constant angular mo-mentum per unit mass, i.e., ̟ Ω ≈ const . The contours of con-stant angular velocity tends to approach a pattern that is closeto cylindrical, as will be discussed below in the context of rapidrotation. Given that Ω ∝ ̟ − , the angular velocity differencebetween the r in and R is ∆Ω = (1 − r /R )Ω = 0 .
75 Ω .Therefore, we have for the second dynamo parameter in Eq. (20)the values C Ω = 75 , , and for ˜Ω = 0 . , , and , re-spectively. The first dynamo parameter in Eq. (20) is C α = 125 ,where we have used ˜ α ≡ α /c s = 0 . for Model A, and η T = 8 × − r ∗ c s . The study of models at rapid rotation is motivated by the in-terest in understanding the evolution of magnetic activity ofyoung stars, i.e., before they have slowed down to the solar ro-tation rate. For us, there is also another motivation in that allour models were of α 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 ro-tation rate could be increased. Another possibility is to lower C α .However, to prevent the dynamo from decaying, one would needto decrease η T even further, but this is computationally difficult.For rapid rotation, the magnetic field lines and contoursof the toroidal magnetic field are much more concentrated tothe bottom of the dynamo region, r ≈ r in . At faster rotation,the contours become more cylindrical. This is an effect of theTaylor–Proudman theorem and results generally in small varia-tions along the rotation axis.The Taylor–Proudman theorem applies primarily the angularvelocity contours. This can be seen by writing the relevant partof the U · ∇ U nonlinearity of Eq. (2) in the form ˆ φ · ∇ × (cid:0) − U · ∇ U (cid:1) p = ̟ ∂∂z Ω + ..., (35)where Ω = U φ /̟ is the local angular velocity, and the dots indi-cate the presence of other terms not relevant here. In Figure 10(a)we show contours of Ω together with a color-coded represen-tation of U r . We see that the Ω contours are already stronglycylindrical for ˜Ω = 1 . As we increase the value of ˜Ω to 10, thecylindrical contours begin to extent much further out along therotation axis; see Figure 10(b).For ˜Ω = 10 , the radial velocity develops a marked indenta-tion inside of what is known as the inner tangent cylinder where ̟ ≥ r in (inner tangent cylinder) ; (36)
9. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 10.
Angular velocity contours superimposed on a color rep-resentation of U r ( r, θ ) for Model B (a) with ˜Ω = 1 and Model C(b) with ˜Ω = 10 . Fig. 11.
Magnetic field lines superimposed on a color represen-tation of B φ ( r, θ ) for Model B with ˜Ω = 1 .see Figure 10(b). Here the outflow is suppressed and supersonicflows occur only for z ≥ r ∗ ≈ r out , i.e., near the outer bound-ary of the computational domain.It turns out that our models are now no longer oscillatoryand are thus still not of α Ω type, contrary to what was originallyhoped for; see Figures 11 and 12. We find C α = 125 and C Ω =3750 for η T = ν T = 8 × − ; see Table 1. To get an idea aboutthe latitudinal variation of the magnetic field in the wind, the plot ˙ E M as a function of θ for different radii. The result is shown in Fig. 12.
Similar to Figure 11, but for Model C with ˜Ω = 10 .Figure 13. It turns out that the magnetic activity is confined to anarrow cone with an opening angle of about ◦ .Noticeable magnetic energy losses are found only near therotation axis; see the blue dotted lines in Figure 14. As a functionof radius, similarly to the case of slow rotation, ˙ E M ( r ) has amaximum somewhere in R < r < r ∗ , which is where the Alfv´enpoint lies. Furthermore, Model B has a much smaller magneticenergy loss at large radii than Model A. Here, we have multiplied ˙ E M ( r ) by a factor of to make it comparable to the values of ˙ E K ( r ) .The model shows similarities with earlier simulationsof outflows emanating from stellar accretion disc dynamos(von Rekowski et al., 2003, 2004), but there the opening anglewas closer to ◦ . In the present simulations, the opening angle Table 1.
Summary of the simulations discussed in this paper.
Model ˜ α Q α ˜Ω C α C Ω B max P cyc A .
05 10 − . − . −
10 250 3750 8.8 —10. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. 13.
Latitudinal dependence of the magnetic energy loss forModel B.
Fig. 14.
Radial dependence of ˙ M (a) and ˙ E M (b) for ≤ θ ≤ (close to the axis) as blue dotted lines, ≤ θ ≤ , as blacksolid lines, and ≤ θ ≤ (close to the equator) as red dashedlines for Model B. The radial dependence of ˙ E K is shown asgreen solid lines.is essentially zero. It corresponds to a cylinder in which most ofthe magnetic fields are ejected, although the flow speed is herestrongly reduced.
4. Conclusions
Our work has shown that a simplified realization of a dynamowith a stellar wind can easily be treated self-consistently inone and the same model, provided certain compromises are be-ing made. The assumption of an isothermal equation of statehas simplified matters conceptionally. Relaxing this restrictionwould allow us to include the energy deposition in the coronaand to model the effects of a sharp density drop at the stellar sur-face. 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 turbu-lent magnetic diffusivity and viscosity. This was mainly neededto resolve shocks that develop within the wind. Those typicallyemerged in response to rapid changes in the magnetic field. Thiscould probably be avoided by allowing for an additional shockviscosity, but this has been avoided in the present work. On theother hand, the angular momentum flux associated with turbu-lent viscosity was already negligible, so its presence may nothave caused any artifacts.Future work might be the inclusion of a Λ effect (R¨udiger,1980, 1989), which would allow for the development of differen-tial rotation in the stellar envelope. Without including the effectsof stellar winds, such models with combined α and Λ effectswere studied by Brandenburg et al. (1990, 1991), who found sig-nificant alignment of the Ω contours with the rotation axis unlessthe 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; seeWarnecke et al. (2013) for spherical convection simulations witha simplified representation of a stellar corona.The inclusion of the Λ effect might allow us to model thestellar dynamo more realistically. It might then allow us to studydynamos in the α Ω regime. This has not been possible in thepresent model for reasons that we are not entirely clear, becausethe value of C Ω was thought to be already large enough. Therecould have been other side effects arising from the coupling tothe outflow that are not yet fully understood.Another important aspect requiring further attention is thestudy of angular momentum losses. Our work has shown that theangular momentum loss can be quantified in terms of the nondi-mensional Shakura–Sunyaev parameter. This is a somewhat un-usual concept in the context of stellar winds, but it may helpputting the theories of turbulent stellar winds and accretion diskon a common footing. Acknowledgements.
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 re-sources provided by the Swedish National Allocations Committee at the Centerfor Parallel Computers at the Royal Institute of Technology in Stockholm.
References
Biermann, L. 1951, Z. f. Ap., 29, 274Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1Brandenburg, A., Moss, D., R¨udiger, G., & Tuominen, I. 1990, Solar Phys., 128,243251Brandenburg, A., Moss, D., R¨udiger, G., Tuominen, I. 1991, Geophys.Astrophys. Fluid Dyn., 61, 179Brandenburg, A., Moss, D., & Tuominen, I. 1992, A&A, 265, 328Brandenburg, A., Subramanian, K., Balogh, A., & Goldstein, M. L. 2011, ApJ,734, 9Chamberlain, J. W. 1960, ApJ, 131, 47Cole, E., Brandenburg, A., K¨apyl¨a, P. J., & K¨apyl¨a, M. J. 2016, A&A, 593, A134Del Sordo, F., Guerrero, G., & Brandenburg, A. 2013, MNRAS, 429, 1686Frank, 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, L76Jabbari, S., Brandenburg, A., Kleeorin, N., Mitra, D., & Rogachevskii, I. 2015,ApJ, 805, 166K¨apyl¨a, P. J., K¨apyl¨a, M. J., & Brandenburg, A. 2014, A&A, 570, A43Krause, F., & R¨adler, K.-H. 1980, Mean-field Magnetohydrodynamics andDynamo Theory (Oxford: Pergamon Press)Malkus, W. V. R., & Proctor, M. R. E. 1975, J. Fluid Mech., 67, 417Mestel, L. 1968, MNRAS, 138, 359Mestel, L. 1999, Stellar Magnetism (Clarendon Press, Oxford)Metcalfe, T. S., & van Saders, J. 2017, Solar Phys., 292, 126
11. Jakab and A. Brandenburg: The effect of a dynamo-generated field on the Parker wind
Fig. A.1.
Effective potential Φ eff = − c ln r − GM/r (black),together with the Bernoulli term − c ln r (red), and the gravi-tational potential − GM/r (blue).
Mitra, D., Moss, D., Tavakol, R., & Brandenburg, A. 2011, A&A, 526, A138Parker, E. N. 1958, ApJ, 128, 664Perri, B., Brun, A. S., R´eville, V., & Strugarek, A. 2018, J. Plasma Phys., 84,765840501R´eville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798,116Roberts, P. H. 1972, Phil. Trans. R. Soc., A272, 663R¨udiger, G. 1980, Geophys. Astrophys. Fluid Dyn., 16, 239R¨udiger, G. 1989, Differential rotation and stellar convection: Sun and solar-typestars (Gordon & Breach, New York)See, V., Matt, S. P., Finley, A. J., Folsom, C. P., Boro Saikia, S., Donati, J.-F., Fares, R., H´ebrard, ´E. M., Jardine, M. M., Jeffers, S. V., Marsden, S. C.,Mengel, M. W., Morin, J., Petit, P., Vidotto, A. A., Waite, I. A., the BCoolCollaboration 2019, ApJ, 886, 120Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Shore, S. N. 1992, An introduction to astrophysical hydrodynamics (San Diego:Academic Press)Stix, M. 1974, A&A, 37, 121van 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, 181von Rekowski, B., Brandenburg, A., Dobler, W., & Shukurov, A. 2003, A&A,398, 825von Rekowski, B., & Brandenburg, A. 2004, A&A, 420, 17Warnecke, J., K¨apyl¨a, P. J., Mantere, M. J., & Brandenburg, A. 2013, ApJ, 778,141Weber, E. J., & Davis, L., Jr. 1967, ApJ, 148, 217
Appendix A: Effective solar wind potential
The purpose of this appendix is to illustrate the effective po-tential experienced by the siphon flow in the isothermal Parkerwind. By comparison, in the case of Roche-lobe overflow, theone-dimensional Bernoulli equation can be written as (cid:0) u − c (cid:1) dd r ln u = − dΦ eff d r , (A.1)where u = U r and Φ eff ( r ) is just the usual gravitational po-tential in a binary system between the two components. For theParker wind, this term takes the form Φ eff = − c ln r − GM/r, (A.2)which is plotted in Figure A.1 along with its two contributions.Evidently, − GM/r steadily increases. Thus, to overcome thegravitational pull, the first term is needed. It becomes important at large radii and is then negative. It characterizes the rapid ex-pansion of the cross-sectional area, πr , and has then the sameeffect as the increase of the cross-sectional area in a Laval noz-zle; see Shore (1992)., and has then the sameeffect as the increase of the cross-sectional area in a Laval noz-zle; see Shore (1992).