# Stellar winds, dead zones, and coronal mass ejections

## 1 Introduction

The solar wind outflow presents a major challenge to numerical modeling since it is a fully three-dimensional (3D), time-dependent physical environment, where regions of supersonic and subsonic speeds coexist in a tenuous, magnetized plasma. Ulysses observations (McComas et al. swoops (1998)) highlighted again that the solar wind about the ecliptic plane is fundamentally dynamic in nature, while the fast speed wind across both solar poles is on the whole stationary and uniform. Recent SOHO measurements (Hassler et al. hassler (1999)) demonstrated how the fast wind emanating from coronal holes is rooted to the ‘honeycomb’ structure of the chromospheric magnetic network, making the outflow truly 3D, while the daily coronal mass ejections are in essence highly time-varying. Moreover, one really needs to study these time-dependent, multi-dimensional aspects in conjunction with the coronal heating puzzle (Holzer & Leer holzer (1997)).

Working towards that goal, Wang et al. (wang (1998)) recently modeled the solar wind using a two-dimensional, time-dependent, magnetohydrodynamic (MHD) description with heat and momentum addition as well as thermal conduction. Their magnetic topology shows both open (polar) and closed (equatorial) field line regions. When heating the closed field region, a sharp streamer-like cusp forms at its tip as the region continuously expands and evaporates. A quasi-stationary wind model results where the emphasis is on reaching a qualitative and quantitative agreement with the observed latitudinal variation (reproducing, in particular, the sharp transition at roughly $`\pm 20^{}`$ latitude between fast and slow solar wind) by tuning the spatial dependence of the artificial volumetric heating and momentum sources.

We follow another route towards global solar wind modeling, working our way up stepwise from stationary 1D to 3D MHD configurations. In a pure ideal, stationary, axisymmetric MHD approach, numerical simulations can benefit greatly from analytical theory. This is demonstrated in Ustyugova et al. (love (1998)), where stationary magneto-centrifugally driven winds from rotating accretion disks were calculated numerically and critically verified by MHD theory.

In this paper, we extend the Wang et al. (wang (1998)) modeling efforts to 2.5D, by including toroidal vector components while remaining axisymmetric. This allows us to explore stellar wind regimes where rotation is also important. The magnetic field still has open and closed field line regions, but in ideal MHD, the closed field region is a ‘dead’ zone from where no plasma can escape. The unknown coronal heating is avoided by assuming a polytropic equation of state and dropping the energy equation all together. The stationary, axisymmetric, polytropic MHD models are analysed as in Ustyugova et al. (love (1998)).

In particular, we investigate the effects of (i) having both open and closed field line regions in axisymmetric stellar winds; and of (ii) time-dependent perturbations within these transonic outflows. While we still ignore the basic question of why there should be a hot corona in the first place, we make significant progress towards fully 3D, dynamic models. The advantages of a gradual approach towards such ‘final’ model were pointed out in Keppens & Goedbloed (1999a ). There, we initiated our effort to numerically model stellar outflows by gradually relaxing the assumptions inherent in the most well-known solar wind model: the isothermal Parker wind (Parker parker (1958)). In a sequence of stationary, 1D, 1.5D, and 2.5D, hydrodynamic and magnetohydrodynamic stellar wind models, all obtained with the Versatile Advection Code (VAC, see Tóth vacapjl (1996, 1997); Tóth & Keppens vac3 (1998); Keppens & Tóth 1999a , and http://www.phys.uu.nl/t̃oth), we demonstrated that we can now routinely calculate axisymmetric magnetized wind solutions for (differentially) rotating stars. An important generalization of previous modeling efforts (Sakurai sakuraiAA (1985, 1990)) is that the field topology can have both open and closed field line regions, so we model ‘wind’ and ‘dead’ zones self-consistently. In essence, our work extends the early model efforts by Pneuman & Kopp (pneu (1971)) in (i) going from an isothermal to a polytropic equation of state; (ii) allowing for stellar rotation; and (iii) including time-dependent phenomena. While we get qualitatively similar solutions for solar-like conditions, we differ entirely in the numerical procedure employed and in the way boundary conditions are specified. Keppens & Goedbloed (1999a ) contained one such MHD model for fairly solar-like parameter values. In this paper, we start with a critical examination of this ‘reference’ model. The obtained transonic outflow, accelerating from subslow to superfast speeds, must obey the conservation laws predicted by theory, by conserving various physical quantities along streamlines. This will be checked in Sect. 2. Section 3 continues with a physical analysis of the model and investigates the influence of the magnetic field strength and of the latitudinal extent of the ‘dead’ zone. These parameters have a clear influence on the global wind structure, especially evident in the appearance and location of its critical surfaces where the wind speed equals the slow, Alfvén, and fast magnetosonic speeds. We also present one such wind solution for a star which rotates twenty times faster than our sun. Finally, Sect. 4 relaxes the stationarity of the wind pattern, by forcing coronal mass ejections on top of the wind pattern. Conclusions are given in Sect. 5.

## 2 Reference model and conservation laws

### 2.1 Solution procedure

We recall from Keppens & Goedbloed (1999a ) that we solve the following conservation laws for the density $`\rho `$, the momentum vector $`\rho \text{v}`$, and the magnetic field B:

$$\frac{\rho }{t}+(\rho \text{v})=0,$$
(1)   
$$\frac{(\rho \text{v})}{t}+[\rho \text{v}\text{v}+p_{tot}\text{I}\text{B}\text{B}]=\rho \text{g},$$
(2)   
$$\frac{\text{B}}{t}+(\text{v}\text{B}\text{B}\text{v})=0.$$
(3)   
Here, $`p_{tot}=p+\frac{1}{2}B^2`$ is the total pressure, I is the identity tensor, $`\text{g}=(GM_{}/r^2)\widehat{𝐞}_r`$ is the external stellar (mass $`M_{}`$) gravitational field with $`r`$ indicating radial distance. We assume $`p\rho ^\gamma `$ (dimensionless, we take $`p=\rho ^\gamma /\gamma `$), where in this paper we only construct models for specified polytropic index $`\gamma =1.13`$. This compares to the value 1.05 used in recent work by Wu, Guo, & Dryer (wusp97 (1997)) and an empirically determined value of 1.46 derived from Helios 1 data by Totten, Freeman, & Arya (totten (1995)).

The discretized Eqs. (1)–(3) are solved on a radially stretched polar grid in the poloidal plane using a Total Variation Diminishing Lax-Friedrich discretization (see e.g. Tóth & Odstrčil vac1 (1996)) with Woodward limiting (Collela & Woodward wood (1984)). Stationary ($`/t=0`$) solutions are identified when the relative change in the conservative variables from subsequent time steps drops below a chosen tolerance (sometimes down to $`10^7`$). We explained in Keppens & Goedbloed (1999a ) how we benefitted greatly from implicit time integration (see Tóth, Keppens, & Botchev implvac2 (1998); Keppens et al. implvac1 (1999); van der Ploeg, Keppens, Tóth auke (1997)) for obtaining axisymmetric ($`/\phi =0`$) hydrodynamic ($`\text{B}=0`$) stellar outflows characterized by $`\rho (R,Z)`$ and $`\text{v}(R,Z)`$, where $`(R,Z)`$ are Cartesian coordinates in the poloidal plane. Denoting the base radius by $`r_{}`$, these hydrodynamic models cover $`r[1,50]r_{}`$ and have as escape speed $`v_{\mathrm{esc}}=\sqrt{2GM_{}/r_{}}=3.3015c_s`$, with $`c_s`$ the base sound speed. They are also characterized by a rotational parameter $`\zeta =\mathrm{\Omega }_{}r_{}/c_s=0.0156`$ (if not specified otherwise), and impose boundary conditions at the base such that (i) $`v_\phi =\zeta R`$; and (ii) the poloidal base speed $`\text{v}_p`$ is in accordance with a prescribed radial mass flux $`\rho \text{v}_p=f_{\mathrm{mass}}\widehat{𝐞}_r/r^2`$. The value of the mass loss rate parameter $`f_{\mathrm{mass}}`$ is taken from a 1D polytropic, rotating Parker wind valid for the equatorial regions under identical parameter values. For $`\zeta =0.0156`$, we get $`f_{\mathrm{mass}}=0.01377`$. We clarify below in which way the values for the dimensionless quantities $`v_{\mathrm{esc}}/c_s`$, $`\zeta `$, and $`f_{\mathrm{mass}}`$ relate to the prevailing solar conditions.

To arrive at a ‘reference’ MHD wind solution, two more parameters enter the description which quantify the initial field strength and the desired extent of the ‘dead’ zone. A stationary, axisymmetric, magnetized stellar wind is the end result of a time stepping process which has the initial density $`\rho (R,Z)`$ and toroidal velocity component $`v_\phi (R,Z)`$ from the HD solution with identical $`\gamma `$, $`v_{\mathrm{esc}}`$, and $`\zeta `$ parameters. The poloidal velocity is also copied from the HD solution in a polar ‘wind’ zone where $`\theta <\theta _{\mathrm{wind}}`$ (upper quadrant with $`\theta =0`$ at pole), quantified by its polar angle $`\theta _{\mathrm{wind}}`$. The ‘dead’ zone is appropriately initialized by a zero poloidal velocity. The field is initially set to a monopole field in the ‘wind’ zone, where

$$B_R(R,Z;t=0)=B_0R/r^3,B_Z(R,Z;t=0)=B_0Z/r^3,$$
(4)   
where $`r^2=R^2+Z^2`$, coupled to a dipole field in the ‘dead’ zone with

$$B_R(R,Z;t=0)=3a_d\frac{ZR}{r^5},B_Z(R,Z;t=0)=a_d\frac{(2Z^2R^2)}{r^5}.$$
(5)   
The strength of the dipole is taken $`a_d=B_0/(2\mathrm{cos}\theta _{\mathrm{wind}})`$ to keep the radial field component $`B_r`$ continuous at $`\theta =\theta _{\mathrm{wind}}`$. The initial $`B_\phi `$ component is zero throughout. Keppens & Goedbloed (1999a ) took $`B_0=3.69`$ and $`\theta _{\mathrm{wind}}=60^{}`$, so that the corresponding dead zone covered only a $`\pm 30^{}`$ latitudinal band about the stellar equator. For the sun at minimum activity, the extent of the coronal hole is typically such that $`\theta _{\mathrm{wind}}=30^{}`$, so it will be useful to vary this parameter in what follows (Sect. 3). We use a resolution of $`300\times 40`$ in the full poloidal halfplane, and impose symmetry conditions at both poles, and free outflow at the outer radius $`50r_{}`$ (where all quantities are extrapolated linearly in ghost cells). At the stellar base, we similarly extrapolate density and all magnetic field components from their initial values, but let these quantities adjust in value while keeping this initial gradient in the ghost cells. This implies that the density and the magnetic field at the base is determined during the time stepping process to arrive at steady-state. We enforce the $`\text{B}=0`$ condition using a projection scheme (Brackbill & Barnes barnes (1980)), to end up with a physically realistic magnetic configuration (despite the ‘monopolar’ field in the wind zone). The stellar boundary condition for the momentum equation allows us to specify a differential rotation rate $`\zeta (\theta )`$ and latitudinally varying mass flux through $`f_{\mathrm{mass}}(\theta )`$. We set

$$\rho \text{v}_p=f_{\mathrm{mass}}(\theta )\widehat{𝐞}_r/r^2,v_\phi =\zeta (\theta )R+B_\phi v_p/B_p.$$
(6)   
The reference model has a rigid rotation rate according to $`\zeta =0.0156`$, while $`f_{\mathrm{mass}}=0.01377`$ in the wind region and zero in the equatorial dead zone.

As emphasized in Keppens & Goedbloed (1999a ), our choice of boundary conditions is motivated by the variational principle governing all axisymmetric, stationary, ideal MHD equilibria (see Sect. 2.2). The analytic treatment shows that the algebraic Bernoulli equation, together with the cross-field momentum balance, really determine the density profile and the magnetic flux function concurrently. In keeping with this formalism, we impose a base mass flux and a stellar rotation, and let the density and the magnetic field configuration adjust freely at the base. In prescribing the stellar rotation, we exploit the freedom available in the variational principle by setting a flux function at the base. Noteworthy, the Pneuman & Kopp (pneu (1971)) model, as well as many more recent modeling efforts for stellar MHD winds, fix the base normal component of the magnetic field together with the density. Below, we demonstrate that our calculated meridional density structure compares well with recent observations by Gallagher et al. (gall (1999)).

The values for the dimensionless parameters $`v_{\mathrm{esc}}/c_s`$, $`\zeta `$ and $`B_0`$ (actually the ratio of the coronal Alfvén speed to $`c_s`$) are solar-like in the following sense. At a reference radius $`r_{}=1.25R_{}`$, we take values for the number density $`N_o10^8\mathrm{cm}^3`$, temperature $`T_o=1.5\times 10^6\mathrm{K}`$, coronal field strength $`B_o2\mathrm{G}`$, and rotation rate $`\mathrm{\Omega }_{}=2.998\times 10^6\mathrm{s}^1`$. For $`\gamma =1.13`$ and assuming a mean molecular weight $`\stackrel{~}{\mu }=0.5`$, the base sound speed then turns out to be $`c_s=167.241\mathrm{km}/\mathrm{s}`$, with all dimensionless ratios as used in the reference model. Further, the value for the mass loss rate parameter $`f_{\mathrm{mass}}=0.01377`$ is then in units of $`1.06\times 10^{13}\mathrm{g}/\mathrm{s}`$, so that a split-monopole magnetic configuration leads to a realistic mass loss rate $`\dot{M}4\pi f_{\mathrm{mass}}2.9\times 10^{14}M_{}\mathrm{yr}^1`$. Since the reference model has a constant mass flux in its wind zone, the presence of the dead zone reduces this value by exactly $`(1\mathrm{cos}\theta _{\mathrm{wind}})=1/2`$. Units enter through the reference radius $`r_{}`$, the base sound speed $`c_s`$, and the base density $`\rho _{}=N_om_p\stackrel{~}{\mu }`$ (with proton mass $`m_p`$).

### 2.2 Streamfunctions

The final stationary wind pattern is shown below in Fig. 3 (see also Fig. 5 in Keppens & Goedbloed 1999a ). The physical correctness of this numerically obtained ideal MHD solution can be checked as follows. All axisymmetric stationary ideal MHD equilibria are derivable from a single variational principle $`\delta L=\delta 𝑑V=0`$ with Lagrangian density (Goedbloed, Keppens, Lifschitz eps98 (1998); Keppens & Goedbloed 1999b ):

$$(M^2,\psi ,\psi ;R,Z)=\frac{1}{2R^2}(1M^2)\psi ^2\frac{\mathrm{\Pi }_1}{M^2}+\frac{\mathrm{\Pi }_2}{\gamma M^{2\gamma }}\frac{\mathrm{\Pi }_3}{1M^2}.$$
(7)   
To obtain an analytic ideal MHD solution, the minimizing Euler-Lagrange equations need to be solved simultaneously for the poloidal flux function $`\psi (R,Z)`$ and the squared poloidal Alfvén Mach number $`M^2(R,Z)\rho v_p^2/B_p^2`$. Here, $`\text{B}_p=(1/R)\widehat{𝐞}_\phi \times \psi `$. In contrast with the translationally symmetric case (Goedbloed & Lifschitz hans1 (1997); Lifschitz & Goedbloed hans2 (1997)), the governing variational principle contains factors $`R^2`$, while the profiles $`\mathrm{\Pi }_1`$ and $`\mathrm{\Pi }_3`$ are no longer flux functions. In particular,

$$\mathrm{\Pi }_1\chi ^2\left(H+\frac{R^2\mathrm{\Omega }^2}{2}+\frac{GM_{}}{r}\right),$$
(8)   
$$\mathrm{\Pi }_2\frac{\gamma }{\gamma 1}\chi ^{2\gamma }S,$$
(9)   
$$\mathrm{\Pi }_3\frac{\chi ^2}{2}\left(R\mathrm{\Omega }\frac{\mathrm{\Lambda }}{R}\right)^2.$$
(10)   
where five flux functions $`H,\mathrm{\Omega },S,\mathrm{\Lambda },\chi ^{}`$ enter. These direct integrals of the axisymmetric, stationary ideal MHD equations are:

* the Bernoulli function ($``$ energy)

$$H(\psi )\frac{1}{2}v^2+\frac{\rho ^{\gamma 1}\gamma S}{\gamma 1}\frac{GM_{}}{r}v_\phi ^2+v_\phi B_\phi \frac{v_p}{B_p},$$
(11)
* the derivative of the stream function $`\chi ^{}\chi /\psi `$. Indeed, the poloidal stream function $`\chi (R,Z)`$ necessarily obeys $`\chi (\psi )`$, provided that the toroidal component of the electric field vanishes. These are immediate checks on the numerical solution, namely $`v_RB_Z=v_ZB_R`$, or the fact that streamlines and field lines in the poloidal plane must be parallel (easily seen in Fig. 3).
* the entropy $`S`$, which for our polytropic numerical solutions is constant by construction: $`S1/\gamma `$,
* a quantity related to the angular momentum flux $`\text{F}_{\mathrm{AM}}=\rho \text{v}_pRv_\phi \text{B}_pRB_\phi \rho \text{v}_p\mathrm{\Lambda }`$, defined as

$$\mathrm{\Lambda }(\psi )Rv_\phi RB_\phi \frac{B_p}{\rho v_p},$$
(12)
* and the derivative of the electric field potential

$$\mathrm{\Omega }(\psi )\frac{1}{R}\left(v_\phi \frac{v_p}{B_p}B_\phi \right).$$
(13)

Various combinations of these flux functions can be made, for instance Goedbloed & Lifschitz (hans1 (1997)) used the following flux function (instead of $`\mathrm{\Lambda }`$)

$$I(\psi )RB_\phi Rv_\phi \frac{\rho v_p}{B_p}=\chi ^{}\mathrm{\Lambda }.$$
(14)   
Figure 1 presents gray-scale contour plots of three streamfunctions $`I`$ (actually $`\mathrm{log}I`$), $`\mathrm{\Lambda }`$, and $`H`$, calculated from the reference solution, with its poloidal field lines overlaid. Ideally, these contours must match the field line structure exactly: all deviations are due to numerical errors. Inspection shows that the agreement is quite satisfactory. A quantitative measure of the errors is given in the lower two frames. At bottom left, we plotted the relative deviation of $`\mathrm{\Omega }`$ from the value enforced at the base $`\zeta (\theta )`$ (constant to 0.0156 in this model). The solid line marks the 10% deviation, actual values range from \[0.0092, 0.0187\]. It should be noted that this is a very stringent test of the solution, since for the chosen parameters, the wind is purely thermally driven, and the stellar rotation is dynamically unimportant. The largest deviations are apparent at the rotation axis (symmetry axis) in $`\mathrm{\Omega }`$, which is not unexpected due to its $`1/R`$ dependence. Other inconsistencies are in the region which has drastically changed from its initial zero velocity, purely dipole magnetic field structure: open field lines coming from the polar regions are now draped around a dipolar ‘dead’ zone of limited radial extent. This dead zone simply corotates with the base angular velocity, and has a vanishing poloidal velocity $`v_p`$ and toroidal field $`B_\phi `$. Around that zone, the stellar wind traces the open field lines. The final bottom right frame shows $`E_\phi `$, virtually vanishing everywhere, and only at the very base are values of order $`𝒪(10^2)`$. For completeness, we show the plasma beta $`\beta =2p/B^2=1`$ contour which exceeds unity in an hourglass pattern that stretches out from the dead zone to large radial distances.

Overall, the obtained stationary numerical solution passes all criteria for being physically acceptable. We expect that most errors disappear when using a higher resolution. We already exploited a radial grid accumulated near the stellar surface, necessary for resolving the near-surface acceleration. However, we could benefit also from a higher resolution in polar angle, now only 40 points for the full half circle, by, for instance, using the up-down symmetry.

## 3 Extensions of the reference model

With the accuracy of the numerical solutions confirmed by inspection of the streamfunctions, we can start the discussion of the influence of the physical parameters $`B_0`$, $`\theta _{\mathrm{wind}}`$, and $`\zeta `$ on the global wind structure. First, we present a more detailed analysis of the reference solution itself.

Fig. 2 shows the density structure at left, where we plot number density as a function of polar angle for three fixed radial distances, namely at the base $`1.27R_{}`$, at $`11.9R_{}`$ and at $`12.7R_{}`$. Keppens & Goedbloed (1999a ) already demonstrated the basic effect visible here: the equatorial density is higher than the polar density. A recent determination of the $`(r,\theta )`$ dependence of the coronal electron densities within $`1R_{}r1.2R_{}`$ by Gallagher et al. (gall (1999)) concluded that the density falloff is faster in the equatorial region than at the poles, and that the equatorial densities within the observed region are a factor of three larger than in the polar coronal hole. Their study lists number densities of order $`10^8\mathrm{cm}^3`$, as in our reference model. Our base density at $`1.27R_{}`$ (dotted in Fig. 2, left panel) has a distinct latitude variation reflecting the combined open-closed field line structure. Interestingly, the observations in Gallagher et al. (gall (1999)) show a similar structure, with quoted values of $`8.3\times 10^7\mathrm{cm}^3`$ at $`1.2R_{}`$ above the pole, increasing to $`1.6\times 10^8\mathrm{cm}^3`$ at the same distance along the equator. In fact, a dip in the density variation was present due to an active region situated above the equator. Qualitatively, we recover this variation at the boundary of the dead zone. Again, we stress that the base density is calculated self-consistently, hence not imposed as a base boundary condition. The other two radial cuts situated beyond the dead zone agree quite well with the conclusions drawn by the observational study.

The reference wind solution also conforms with some well-known studies in MHD wind modeling. Suess & Nerney (suess (1973)) and Nerney & Suess (nern (1975)) pointed out how consistent axisymmetric stellar wind modeling which include magnetic fields and rotation automatically lead to a meridional flow away from the equator. At large radial distances, the flow profile should be of the form $`v_\theta \mathrm{sin}(2\theta )`$, with a poleward collimation of the magnetic field. This variation in polar angle is general and independent of the precise base field structure. In the middle panel of Fig. 2, we show the latitude dependence of $`v_\theta `$ at $`50R_{}`$ for the reference model. Note the perfect agreement with the predicted variation.

Since the solar rotation rate, quantified by the parameter $`\zeta `$, is low, the calculated wind solution in the poloidal plane should be similar to the one presented by Pneuman & Kopp (pneu (1971)). They constructed purely poloidal, isothermal and axisymmetric models of the solar wind including a helmet streamer (or ‘dead’ zone). An iterative technique was used to solve for the steady coronal expansion, while the density and the radial magnetic field were fixed at the base. They enforced a dipolar $`B_r`$ with a strength of $`1\mathrm{G}`$ at the poles, half the value we use at the initialization. Their uniform coronal temperature was taken to be $`1.56\times 10^6\mathrm{K}`$, almost identical to our base temperature $`T_o`$. Their base number density was imposed to be $`1.847\times 10^8\mathrm{cm}^3`$, independent of latitude, and they assumed a slightly higher value for the mean molecular weigth, namely $`\stackrel{~}{\mu }=0.608`$. This leads to a base density which is a factor of 2.246 higher than the one used in our model. With these differences in mind (together with our polytropic equation of state and the rotational effects), we show in the right panel of Fig. 2 the magnetic structure and the location of the sonic (where $`v_p=c_s`$) and the Alfvénic surface (where $`v_p=B_p/\sqrt{\rho }`$) in a manner used in the original publication of Pneuman & Kopp (pneu (1971)), their Fig. 4. In the $`(1/r,\mathrm{cos}\theta )`$ projection, the Alfvénic transition on the equator is at the cusp of the helmet structure before the sonic point, while the sonic surface is closer to the solar surface at the poles. The qualitative agreement is immediately apparent, although our solution method is completely different, most notably in the prescription of the boundary conditions. By calculating the base density and magnetic field configuration self-consistently, we generalize the solution procedure employed by Pneuman & Kopp (pneu (1971)) as we gain control of the size of the dead zone through our parameter $`\theta _{\mathrm{wind}}`$. This allows us to study the influence of the base topology of the magnetic field on the global wind acceleration pattern in what follows.

Fig. 3 confronts three steady-state wind solutions with our reference model, which differ in the latitudinal extent of the dead zone and/or in the magnetic field strength. With $`\theta _{\mathrm{wind}}=60^{}`$ and $`B_0=3.69`$ (corresponding to a $`2\mathrm{G}`$ base coronal field strength) for the reference case A, we increased the latitudinal extent of the dead zone by taking $`\theta _{\mathrm{wind}}=30^{}`$ in case B, doubled the field strength parameter $`B_0`$ in model C, and took both $`B_0=7.4`$ and $`\theta _{\mathrm{wind}}=30^{}`$ to arrive at model D. We recall that $`B_0`$ specifies only the initial field strength used in the time-stepping process towards a stationary solution. The final base field strength turns out to be of roughly the same magnitude, but differs in its detailed latitudinal variation. The changes in the gobal wind pattern are qualified by the resulting deformations of the critical surfaces (hourglass curves) where the wind speed equals the slow, Alfvén, and fast speeds. The plotted region stretches out to $`18r_{}22.5R_{}`$.

By enlarging the dead zone under otherwise identical conditions (from A to B), the polar, open field lines are forced to fan out more rapidly with radial distance. As a result, the acceleration of the plasma occurs closer to the stellar surface, and the critical curves become somewhat more isotropic in polar angle. The Alfvén surface moves inwards at the poles, and shifts outwards above the now larger dead zone at the equator, approaching a circle with an equatorial imprint of the dead zone. If we keep the dead zone small, but double the initial field strength $`B_0`$ (from A to C), the opposite behaviour occurs: the critical curves, hence the entire acceleration behaviour of the wind, become much more anisotropic. The most pronounced change is an inward shift of the polar slow transition, and an outward shift of the Alfvén and fast polar transition. This behaviour is in agreement with what a Weber-Davis model (Weber & Davis wd (1967)) predicts to happen when the field strength is increased (note that the Weber-Davis model only applies to the equatorial region). When both the field strength and the dead zone are doubled (from A to D), the resulting Alfvén and fast critical curves are rather isotropic due to the influence of the dead zone. The polar slow transition is displaced inward while the polar Alfvén and fast curve are shifted outward, as expected for the higher field. The detailed equatorial behaviour is clearly modulated by the existing dead zone. We note that all wind solutions presented are still thermally driven, since the solar-like rotation rate is rather low and the field strengths are very modest. The changes are entirely due to reasonable variations in magnetic field topology and only a factor of two in field strength. One could tentatively argue that such variations occur in the solar wind pattern within its 11-year magnetic cycle. In gray-scale, Fig. 3 shows the absolute value of the toroidal field component $`B_\phi `$ (this field changes sign across the equator). Note that the stellar rotation has wound up the field lines in a zone midway between the poles and the dead zone. For higher rotation rates (see below), the associated magnetic pressure build up due to rotation can influence the wind pattern and cause collimation (Trussoni, Tsinganos, & Sauty trusso (1997)). Due to the four-lobe structure, one can expect parameter regimes which lead to both poleward collimation (as in the monopole-field models of Sakurai sakuraiAA (1985)) and equatorward streamline bending.

Figure 4 compares the radial dependence of the poloidal velocity for the four models (A, B, C, D) at the pole (left panel) and the equator (middle panel). For comparison, we overplotted the same quantities in each panel for a solution with a split-monopole base field at the same parameter values. This monopolar field solution is identical in nature to the Sakurai (sakuraiAA (1985, 1990)) models and was shown in Keppens & Goedbloed (1999a , Fig. 4). At the pole, a faster acceleration to higher speeds as compared to the reference model results from either increasing the field strength or enlarging the dead zone. Moreover, all four models show a faster initial acceleration than the corresponding monopolar field model. Radio-scattering measurements of the polar solar wind speed (Grall et al. grall (1996)) indicated that the polar wind acceleration is almost complete by $`10R_{}`$, much closer than expected. Our model calculations show that a fast acceleration can result from modest increases in the coronal field strength and dead zone extent (model D has a solar-like dead zone of $`\pm 60^{}`$).

The middle panel of Fig. 4 shows the distinct decrease in equatorial wind speed due to the dead zone, when compared to a split-monopole solution. The equatorial velocities are reduced by $`10`$ to $`40\mathrm{km}/\mathrm{s}`$, depending on the size of the dead zone and the base field strength. Enlarging the dead zone reduces the wind speed significantly (compare A to B, and C to D). To a lesser degree, the same effect is true for an increase in coronal field strength (compare A to C, and B to D). Keppens & Goedbloed (1999a , Figure 6) contained a polar plot of the velocity and the density at fixed radial distance for the reference model, where at least qualitatively, a transition from high density, low speed equatorial wind to lower density, high speed polar wind is noticable. As evidenced by Fig. 4, this difference in equatorial and polar wind is even more pronounced for larger dead zones. The velocities reached are too low for explaining the solar wind speeds – the Weber-Davis wind solution of identical parameters reaches $`263\mathrm{km}/\mathrm{s}`$ at 1 AU. However, this is a well-known shortcoming of a polytropic MHD description for modeling the solar wind. Wu et al. (wujgr99 (1999)) therefore resort to an ad hoc procedure where the polytropic index is an increasing function of radial distance, $`\gamma (r)`$, to attain a more realistic $`420\mathrm{km}/\mathrm{s}`$ wind speed at 1 AU, corresponding to the ‘slow’ solar wind. A more quantitative agreement with the observations at these distances must await models where we take the energy equation into account and/or model extra momentum addition as in Wang et al. (wang (1998)). The equatorial toroidal velocity profile is shown at right in Fig. 4. Note that increasing $`B_0`$ or enlarging the dead zone both negatively affect the degree to which the corona corotates with the star.

Keppens & Goedbloed (1999a ) also contained a hydrodynamic solution for a much faster rotation rate quantified by $`\zeta =0.3`$, or twenty times the solar rotation rate. The additional centrifugal acceleration moves the sonic transition closer to the star along the equator, and induces an equatorward streamline bending at the base at higher latitudes (see also Tsinganos & Sauty tsinhd (1992)). One could meaningfully ask what remains of this effect when a two-component field structure is present as well. Therefore, we calculate an MHD wind for this rotation rate, with $`B_0=3.69`$ and $`\theta _{\mathrm{wind}}=60^{}`$ as in the reference case. The corresponding mass loss rate parameter (only used in the wind zone) is $`f_{\mathrm{mass}}=0.01553`$. Fig. 5 displays the wind structure as in Fig. 3, with the gray-scale indicating the logarithm of the density pattern. For this solution, we used the up-down symmetry to double the resolution in polar angle at the same computational cost. The shape of the critical curves has changed dramatically, with a significant outward poleward shift of the Alfvén curve, and a clear separation between Alfvén and fast critical curves – in agreement with a 1.5D Weber-Davis prediction. The actual position of the critical surfaces may be influenced by an interaction with the outer boundary at $`50r_{}`$ in the time-stepping towards a stationary solution. In fact, the combined Alfvén-fast polar transition has shifted outside the computational domain, and the residual could not be decreased to arbitrary small values but stagnated at $`𝒪(10^7)`$ following this interaction. Within the plotted region of $`37R_{}`$, the solution is acceptable as explained in Sect. 2. Note how the density structure shows an increase towards the equator, causing a very effective thermo-centrifugal acceleration of the equatorial wind above the dead zone. The equatorward streamline bending occuring in the purely hydrodynamic wind is still important, but now clearly affected by the presence of the dead zone. The toroidal magnetic pressure built up by the stellar rotation along the mid-latitude open field lines is shaping the wind structure as a whole. In those regions, we have $`\rho v_p^2/B_\phi ^2<1`$ together with $`2p/B_\phi ^2<1`$. Thereby, it also leads to streamline bending, both poleward as clearly seen in the high latitude field lines, and equatorwards in the vicinity of the stellar surface. In this way, the magnetic topology consisting of a dead and a wind zone, combined with fast rotation, leads to magnetically dominated collimation along the stellar poles, together with magneto-rotational deflections along the equator. The latter leads to enhanced densities in the equatorial plane.

Figure 6 shows the radial dependence at a polar angle $`\theta =41.6^{}`$ of the radial velocity $`v_r`$, sound speed $`c_s`$, radial Alfvén speed $`A_r=B_r/\sqrt{\rho }`$, and azimuthal speed $`v_\phi `$. Note that the radial velocity reaches up to $`400\mathrm{km}/\mathrm{s}`$ (compare with the $`200\mathrm{km}/\mathrm{s}`$ velocities reached under solar conditions as shown in Fig. 4), as a result of the additional centrifugal acceleration. The rotation plays a significant dynamical role here, in contrast to the ‘solar-like’ models discussed earlier and displayed in Fig. 3. The toroidal speed $`v_\phi `$ reaches above $`100\mathrm{km}/\mathrm{s}`$, a factor of $`10`$ higher than along the ecliptic as shown in Fig. 4 (left panel). The corotation obtained by the numerical procedure to find the stationary state can again be quantified by the relative error $`\mathrm{\Omega }\zeta /\zeta `$: the bottom panel of Fig. 6 proves that it is less than 3 % along that same radial cut.

## 4 Triggering coronal mass ejections

In the process of generating a stationary wind solution, various dynamic phenomena take place which may have physical relevance. For instance, the equatorial conic section delineated by $`\theta [\theta _{\mathrm{wind}},\pi \theta _{\mathrm{wind}}]`$ which was initially static ($`v_p=0`$) and dipolar throughout, first gets ‘invaded’ by plasma emanating from the wind zone. The open field structure is dragged in towards the equator, and most of the dipolar field is moved out of the domain, except for the remaining ‘dead’ zone. One could qualitatively relate some of these changes in the global magnetic topology with observed coronal phenomena.

In reality however, coronal mass ejections represent major disturbances which happen on top of the stationary transonic solar wind. They are associated with sudden, significant mass loss and cause violent disruptions of the global field pattern. Most notably, one frequently observes the global coronal wind structure to return to its previous stationary state, after the passage of the CME. Within the realms of our stellar wind models, we can trigger CMEs on top of the outflow pattern, study their motion, and at the same time demonstrate that the numerical solutions indeed are stable to such violent perturbations by returning to a largely unchanged stationary state. We still restrict ourselves to axisymmetric calculations, so the geometry of our ‘CME’ events is rather artificial. In future work, we intend to model these CMEs in their true 3D setting.

As background stellar wind, we use a slightly modified model B from the previous section. Model B had a large dead zone, $`\theta _{\mathrm{wind}}=30^{}`$, $`B_0=3.69`$, and a rigid rotation with $`\zeta =0.0156`$ corresponding to $`\mathrm{\Omega }_{}3\times 10^6\mathrm{s}^1`$. We changed the boundary condition on $`v_\phi `$ to mimick a ‘solar-like’ differential rotation, by taking

$$\zeta (\theta )=\zeta _0+\zeta _2\mathrm{cos}^2\theta +\zeta _4\mathrm{cos}^4\theta ,$$
(15)   
with $`\zeta _0=0.0156`$, $`\zeta _2=0.00197`$, and $`\zeta _4=0.00248`$. This enforces the equator to rotate faster than the poles in accord with the observations. As expected for this low rotation rate, this has no significant influence on the wind acceleration pattern. The coronal mass ejection is an equally straigthforward modification of the boundary condition imposed on the poloidal momentum equation, namely $`\rho \text{v}_p=f_{\mathrm{mass}}(\theta ,t)\widehat{𝐞}_r/r^2`$ with

$$f_{\mathrm{mass}}(\theta ,t)=f_{\mathrm{wind}}(\theta )+g_{\mathrm{CME}}\mathrm{sin}\left(\frac{\pi t}{\tau _{\mathrm{CME}}}\right)\mathrm{cos}^2\left(\frac{\pi }{2}\frac{\theta \theta _{\mathrm{CME}}}{a_{\mathrm{CME}}}\right),$$
(16)   
for $`0t\tau _{\mathrm{CME}}`$ and $`\theta _{\mathrm{CME}}a_{\mathrm{CME}}\theta \theta _{\mathrm{CME}}+a_{\mathrm{CME}}`$, and otherwise

$$f_{\mathrm{mass}}(\theta ,t)=f_{\mathrm{wind}}(\theta ).$$
(17)   
The wind related mass loss rate $`f_{\mathrm{wind}}(\theta )`$ contains the polar angle dependence due to the dead zone, as before. The extra four parameters control the magnitude of the CME mass loss rate $`g_{\mathrm{CME}}`$, the duration $`\tau _{\mathrm{CME}}`$, and the location $`\theta _{\mathrm{CME}}`$ and extent $`0a_{\mathrm{CME}}\pi /2`$ in polar angle for the mass ejection. We only present one CME scenario for parameter values $`g_{\mathrm{CME}}=2`$, $`\tau _{\mathrm{CME}}=0.5`$, $`\theta _{\mathrm{CME}}=60^{}`$, and $`a_{\mathrm{CME}}=30^{}`$. Note that the up-down symmetry is hereby deliberately broken. This scenario mimicks a mass ejection which detaches from the coronal base within 45 minutes and which has an associated mass flux of about $`2\times 10^{13}\mathrm{g}/\mathrm{s}`$. In fact, the total amount of mass lost due to the CME can be evaluated from

$$M_{\mathrm{lost}}^{\mathrm{CME}}=2g_{\mathrm{CME}}\tau _{\mathrm{CME}}\frac{\pi ^2}{\pi ^2a_{\mathrm{CME}}^2}\left\{\mathrm{cos}(\theta _{\mathrm{CME}}a_{\mathrm{CME}})\mathrm{cos}(\theta _{\mathrm{CME}}+a_{\mathrm{CME}})\right\}.$$
(18)   
For the chosen parameter values, this works out to be $`M_{\mathrm{lost}}^{\mathrm{CME}}=\frac{36}{35}\sqrt{3}`$, corresponding to $`0.98\times 10^{17}\mathrm{g}`$, a typical value for a violent event.

Figure 7 shows the density difference between the evolving mass ejection and the background stellar wind at left, the magnetic field structure (middle), and the toroidal velocity component $`v_\phi `$ (right) at times $`t=1`$ (1hr 27’ after onset) and $`t=3`$ (4hr 20’ after onset). The region is plotted up to $`15r_{}18.75R_{}`$. Although the event is triggered in the upper quadrant dead zone only, its violent character also disturbs the overlying open field (or wind) zone. The added plasma, trapped in the dead zone, even perturbs the lower quadrant wind zone at later times. Note that the CME induces global, abrupt changes in the toroidal velocity component. The outermost closed field lines get stretched out radially, pulling the dead zone along (see Fig. 7-8). In ideal MHD calculations, they can never detach through reconnection, although numerical diffusion can cause it to happen. We observed the outermost field lines of the dead zone to travel outwards without noticable reconnection. The overall wind pattern in the first 15 base radii thereby approximately returns to its original stationary state, as shown in Figure 8 which gives the solution at $`t=5`$ (7hr 13’ after onset) and $`t=30`$ (more than 43 hrs after onset). Figure 9 shows how a hypothetical spacecraft at $`21r_{}26.23R_{}`$, close to the ecliptic, would record the CME-passage as a sudden increase in density and poloidal velocity which eventually relax to their pre-event levels. The event is followed by an increased azimuthal flow regime and shows large amplitude oscillations in magnetic field strength and orientation.

As we assumed axisymmetry, this simulation serves as a crude model for CME-type phenomena. Interestingly, axisymmetric numerical simulations of toroidal flux ‘belts’ launched from within the dead zone of a purely meridional, polytropic MHD wind can relate favourably to satellite magnetic cloud measurements at 1AU (Wu et al. wujgr99 (1999)). Note that we triggered a ‘CME’ by prescribing a time and space dependent mass flux at the stellar base, where the density and the magnetic field components could adjust freely. Alternatively, as used in studies by Mikić & Linker (mikic (1994)), global coronal restructuring can be triggered by shearing a coronal arcade. Parameter studies of axisymmetric, but ultimately 3D solutions, could investigate the formation and appearance of various MHD shock fronts depending on plasma beta, Mach numbers, etc.

## 5 Conclusions and outlook

Continuing our gradual approach towards dynamic stellar wind simulations in three dimensions, we studied the influence of the magnetic field strength and topology (allowing for wind and dead zones), of the stellar rotation, and of sudden mass ejecta on axisymmetric MHD winds.

We demonstrated how reasonable changes in the coronal magnetic field (factor of two in field strength and in dead zone extent) influence the detailed acceleration behaviour of the wind. Larger dead zones cause effective, fairly isotropic acceleration to super-Alfvénic velocities since the polar, open field lines are forced to fan out rapidly with radial distance. The Alfvén transition moves outwards when the coronal field strength increases. The equatorial wind outflow is in these models sensitive to the presence and extent of the dead zone, but has, by construction, a vanishing $`B_\phi `$ and a $`\beta >1`$ zone from the tip of the dead zone to large radial distances. The parameter values for these models are solar-like, hence the winds are mostly thermally driven and, in particular, emanated from slowly rotating stars.

For a twenty times faster than solar rotation rate, the wind structure changes dramatically, with a clear separation of the Alfvén and fast magnetosonic critical curves. At these rotation rates, a pure hydrodynamic model predicts equatorward streamline bending from higher latitudes. Our MHD models show how this is now mediated by the dipolar dead zone. An equatorial belt of enhanced density stretches from above the dead zone outwards, where effective thermo-centrifugally driven outflow occurs. The magnetic field structure shows signs of a strong poleward collimation, due to the significant toroidal field pressure build up at these spin rates. As pointed out by Tsinganos & Bogovalov (tsing99 (1999)), this situation could apply to our own sun in an earlier evolutionary phase.

It could be of interest to make more quantitative parameter studies of the interplay between field topology, rotation rate, etc. in order to apprehend transonic stellar outflows driven by combinations of thermal, magnetic, and centrifugal forces. A systematic study of the angular momentum loss rates as a function of dead zone extent, magnetic field strength, and rotation rate can aid in stellar rotational evolution modeling (Keppens, Charbonneau, & McGregor rotpaper (1995); Keppens binary (1997)). Specifically, Li (jianke (1999)) pointed out that the present solar magnetic breaking rate is consistent with either one of two magnetic topologies: (i) one with the standard coronal field strength of $``$ 1 $`G`$ and a small $`<2R_{}`$ dead zone; or (ii) one with a larger $``$ 5 $`G`$ dipole strength and a sizeable dead zone. When we calculate the torque exerted on the star by the magnetized winds A, B, C and D shown in Fig. 3 as

$$\tau _{\mathrm{wind}}=4\pi _0^{\frac{\pi }{2}}𝑑\theta \mathrm{\Lambda }\rho r^2v_R,$$
(19)   
(noting the axi- and up-down symmetry), we find $`\tau _{\mathrm{wind}}^A0.139\times 10^{31}\mathrm{dyne}\mathrm{cm}`$, $`\tau _{\mathrm{wind}}^B0.062\times 10^{31}\mathrm{dyne}\mathrm{cm}`$, $`\tau _{\mathrm{wind}}^C0.246\times 10^{31}\mathrm{dyne}\mathrm{cm}`$, and $`\tau _{\mathrm{wind}}^D0.123\times 10^{31}\mathrm{dyne}\mathrm{cm}`$. This confirms Li’s result, since a simultaneous doubling of the coronal field strength and the dead zone extent (from model A to D) hardly changes the torque magnitude. As could be expected, only enlarging the dead zone lowers the breaking efficiency (as pointed out in Solanki, Motamen, & Keppens samipap (1997)), while only raising the field strength leads to faster spin-down. Interestingly, a Weber-Davis prediction with governing parameters identical to the reference model A (as presented in Keppens & Goedbloed 1999a ) gives a value $`\tau _{\mathrm{wind}}=4\pi \rho _Ar_A^2v_{rA}\frac{2}{3}\mathrm{\Omega }_{}r_A^22.387\times 10^{31}\mathrm{dyne}\mathrm{cm}`$ (all quantities evaluated at the Alfvén radius $`r_A`$), one order of magnitude larger! The same conclusion was reached by Priest & Pneuman (priest74 (1974)) by estimating the angular momentum loss rate from the purely meridional Pneuman & Kopp (pneu (1971)) model. Although the latter model does not include rotation (so that $`\mathrm{\Lambda }`$ in Eq. (19) is strictly zero for this model), Priest & Pneuman (priest74 (1974)) could estimate the torque for a solar rotation rate from the obtained variation of the poloidal Alfvén radius as a function of latitude (our Fig. 2, right panel). The resulting estimate was only 15 % of that for a monopole base field. Our exactly evaluated spin-down rates are 2.6 % to 10.3 % of a split-monopole case. The large difference arises due to the presence of the dead zone and the fact that $`B_\phi `$ vanishes across the equator for the wind solutions from Fig. 3. Indeed, evaluating the torque from Eq. 19 for the monopolar wind solution from Keppens & Goedbloed (1999a , Figure 4) gives $`\tau _{\mathrm{wind}}2.326\times 10^{31}\mathrm{dyne}\mathrm{cm}`$, in agreement with the Weber-Davis estimate. Hence, it should be clear that full MHD modeling is a useful tool to further evaluate and constrain different magnetic braking mechanisms.

We showed how CME events can be simulated on top of these transonic outflows. The detailed wind structure is stable to violent mass dumps, even when ejected in the dead zone. Note that we restricted ourselves to axisymmetric perturbations, and it will be of interest to show whether the axisymmetric solutions are similarly stable to non-axisymmetric perturbations (as recently investigated for shocked accretion flows on compact objects in Molteni, Tóth, & Kuznetsov molt (1999)). One could then focus on truly 3D mass ejecta and their parametric dependence (possibly allow for direct comparison with LASCO observations of coronal mass ejections), or even experiment with unaligned rotation and magnetic axes. A 3D time-dependent analytic model by Gibson & Low (gibson (1998)) can be used as a further check on the numerics. Alternatively, we may decide to zoom in on (3D) details of the wind structure at the boundaries of open and closed field line regions or about the ecliptic plane, to see whether shear flow driven Kelvin-Helmholtz instabilities (Keppens et al. kh2d (1999); Keppens & Tóth 1999b ) develop in these regions.

The Versatile Advection Code was developed as part of the project on ‘Parallel Computational Magneto-Fluid Dynamics’, funded by the Dutch Science Foundation (NWO) Priority Program on Massively Parallel Computing, and coordinated by JPG. Computer time on the Cray C90 was sponsored by the Dutch ‘Stichting Nationale Computerfaciliteiten’ (NCF). We thank Keith MacGregor for suggesting to compare torque magnitudes for different models and an anonymous referee for making several useful comments.