A GRID OF MHD MODELS FOR STELLAR MASS LOSS AND SPIN-DOWN RATES OF SOLAR ANALOGS

In order to produce a grid of solutions for the coronal and wind structure of solar analogs, we follow the approach by CO09 to calculate the mass loss rate and the angular momentum loss rate. We use the BATS-R-US MHD code (Powell et al. 1999) in its version developed for the solar corona (Cohen et al. 2007), which solves the set of MHD equations for the conservation of mass, momentum, magnetic induction, and energy. The model is driven by boundary conditions for the radial surface magnetic field, where the initial condition for the 3D field in the whole domain is potential (Altschuler & Newkirk 1969). The acceleration of the stellar wind is done via a source term in the energy equation that is scaled semi-empirically with the topology of the initial potential field and the magnetic flux-tube expansion factor according to the relations found for the solar wind by Wang & Sheeley (1990) and Arge & Pizzo (2000) as implemented by Cohen et al. (2007). The code then solves the MHD equations with this source term until the solution is relaxed to a steady state, which resembles the observed solar wind structure (when using solar observations to drive the model). A detailed description of the model and the equations it solves can be found in Cohen et al. (2009) and Cohen et al. (2010a).

The unique feature that makes this model more "realistic" is the semi-empirical relation between the magnetic topology and the wind speed, as well as the base density. The latter is imposed to resemble the observed density difference between the slow and fast wind, and turns out to be crucial in determining the mass loss rate of each solution. In Section 4, we discuss how this density difference in a bi-modal stellar wind is the reason for the discrepancy between a solution for a "dipole solar case" and a solution based on actual solar minimum (nearly a dipole) data. In other models with an assumed spherical wind (and density) distribution imposed on a certain magnetic field topology, the MHD solution relaxes to a steady state only due to pressure balance between the wind pressure and the magnetic pressure. However, the stellar wind is assumed to simply "be there," and no mechanism to drive the wind appears in the equations themselves. In our model, a steady state is obtained due not only to the same pressure balance, but also to a volumetric energy source which is dictated by the magnetic field topology. This way, we obtain a realistic, bi-modal steady-state solution for the wind, with faster wind originating from coronal holes and slower wind coming from the boundaries of the helmet streamers as observed in the solar wind (see, e.g., McComas 2007).

As noted, the advantage of our model is that it successfully reproduces the solar wind in a realistic manner. By assuming that stellar winds of solar analogs are driven by the same process as the solar wind, we use the model to calculate different steady-state solutions for a grid of models, where we employ three free parameters: (1) the stellar dipole field strength, (2) the stellar rotation period, and (3) the value of the coronal base density (i.e., the inner boundary condition for the density in the simulation domain). For the stellar magnetic field strength, B_⋆, we use dipole polar values of 5, 10, 50, 100, and 500 G; for the stellar rotation period, Ω_⋆, we use values of 0.5, 2, 5, 10, and 25 days; and for the coronal base density, n_⋆, we use values of 2e8, 5e8, 1e9, 5e9, 1e10, and 5e10 cm⁻³.

We point out that the values we use for the base density are significantly lower than the value used in Matt & Pudritz (2008), that formed the basis of the calculations in Matt et al. (2012), and in Vidotto et al. (2011), which were in the range 1e9–1e13 cm⁻³. While these high values were obtained from the literature and very active stars do have traces of plasma at such high densities (e.g., Ness et al. 2004; Testa et al. 2004) these measurements are for plasma in closed coronal loops, where the density is higher. Testa et al. (2004) also notes that the filling factor—the surface area covered by such high density plasma—is only of the order of a percent or so. The ambient coronal density, which is the actual source for the wind that leaves the star along open field lines, is much lower, unless the lower boundary of the numerical model is set at the chromosphere, which is not the case in any of the models discussed here. To support this argument, we show that we obtain the solar mass loss rate, $\dot{M}_\odot =2\hbox{--}4\times 10^{-14}\;M_\odot \;{\rm yr}^{-1}$, for typical solar parameters and a base density value at the lower end of our density range, whereas solar active regions that dominate the X-ray emission tend to have plasma densities in the range 1e9–1e10 cm⁻³ (e.g., Doschek 1997; Del Zanna & Mason 2003; Young et al. 2009).

For each BATS-R-US model solution, we calculate the non-spherical and arbitrary Alfvén surface at which the flow speed, u, equals the Alfvén speed,

$$\begin{equation} u_A=B/\sqrt{4\pi \rho }, \end{equation} \tag{ 1 }$$

where B is the magnetic field strength and ρ = nm_p is the mass density (n is the number density and m_p is the proton mass). Beyond this surface, the wind velocity exceeds the Alfvén speed, and the wind is no longer in contact with the star via magnetic fields. The Alfvén surface then represents the surface at which the wind effectively escapes and exerts no more torque on the star. We then integrate the mass flux, ρu, over all the Alfvén surface elements, da_i. This procedure allows us to obtain the total mass loss rate, $\dot{M}$:

$$\begin{equation} \dot{M}=\frac{\partial M}{\partial t}=\int \nabla \cdot (\rho \mathbf {u})dV=\int \rho \mathbf {u}\cdot da=\sum \limits _{i}\rho _i \mathbf {u}_i\cdot \mathbf {da}_i. \end{equation} \tag{ 2 }$$

In a similar manner, we integrate the differential version of the formula for the angular momentum loss rate (Equation (35) in Weber & Davis 1967) to obtain the total angular momentum loss rate, $\dot{J}$:

$$\begin{equation} \dot{J}=\frac{2}{3}\Omega _\star \dot{M}r_A^2=\frac{2}{3}\sum \limits _{i} \Omega _\star \sin ^2{\theta }_i\;r^2_{i} \;\rho _i \mathbf {u}_i\cdot \mathbf {da}_i. \end{equation} \tag{ 3 }$$

Here, r_i is the Alfvén radius for a particular magnetic field line, or the distance as measured from the center of the star to a given point on the Alfvén surface. The sin θ_i term is due to the fact that the torque applied on the star depends on the angle between the lever arm and the rotation axis.

Note that in Equations (3) and (4), r_i is not constant but is the differential distance at each point of integration over the Alfvén surface. The fact that we integrate the differential r_is allows us to remove any assumption of geometrical symmetries and to perform calculations for any complex topology.

Equations (2) and (3) tell us that the mass and angular momentum loss rates depend on the interplay between the actual mass flux through the Alfvén surface (dictated by the density; CO09) and the size of the surface. The size of the surface determines both the total area of integration, as well as the size of the lever arm that applies a torque on a star to spin it down. This is why determining the stellar angular momentum loss rate for a complex field topology is not trivial (CO09).

Figure 1 shows the Alfvén surface and a meridional slice of the simulation domain for selected test cases. The color contours represent the local distribution of the mass loss rate. It can be seen that the size of the Alfvén surface increases with an increase of the magnetic field and decreases with the increase of the base density—this is a direct consequence of the dependence of the Alfvén speed on these quantities (Equation (1)) and the distance from the star at which the Alfvén speed is exceeded by the wind. It can also be seen that the overall mass loss rate increases with increasing base density. This is partly a consequence of mass continuity in that the base density essentially defines the lower wind boundary condition. The effect of the rotation period on the Alfvén surface is mostly geometrical. When moving from slow to fast rotation, the Alfvén surface becomes stretched toward the poles, so that the total surface area increases.

Zoom In Zoom Out Reset image size

We note that some of the Alfvén surfaces in Figure 1 are not fully rotationally symmetric. This is due to artifacts originating from the use of a Cartesian grid of finite resolution for models with rotational symmetry. We have performed a number of test runs for representative cases with much higher resolution that resulted in rotationally symmetric shapes of the Alfvén surfaces. However, comparisons of mass and angular momentum loss rates from these high-resolution solutions with those of the lower resolution calculations have shown differences of not more than 5%. Due to the large number of simulations required for the study and the currently very high computational cost of higher resolution calculations, we have retained the lower resolution approach and consequently some of the results presented here may display residual asymmetric artifacts.

Figure 2 summarizes the results of our grid of models. Each pair of plots shows the mass loss rate and the angular momentum loss rate as a function of magnetic field strength for a given base density value. Each curve on the plots represents a particular stellar rotation period. A small number of test cases for very high density and low magnetic field values are considered unphysical because the equatorial rotation velocity on the stellar surface is already greater than the Alfvén speed (parts of the Alfvén surface are inside the star). These cases are marked with a red spot in the plots of Figure 2.

Zoom In Zoom Out Reset image size

The most notable trend is the rather linear (in log space, i.e., power law) dependence of the mass and angular momentum loss rates on the dipole field strength. The slope of this trend decreases as we increase the base density, and the difference between the mass loss rate for the weakest and strongest field strength changes from about two orders of magnitude for a base density of 2e8 cm⁻³ to about one order of magnitude for base density of 5e10 cm⁻³. A linear dependence of the angular momentum and mass loss rates on the base density explicitly appears in Equations (2) and (3).

To illustrate the effect of each parameter, Figure 3 shows the Alfvén surface, the local mass loss rate distribution, and selected magnetic field lines for the most extreme cases (lowest and highest values for each of the parameters). The parameters of these cases are summarized in Table 1. While we consider case E unphysical, we still find the display of the solution valuable for understanding the general trends of the solutions.

Zoom In Zoom Out Reset image size

First, it can be seen that for slow stellar rotation, the coronal magnetic field lines are stretched and opened up radially by the stellar wind. As we move to very fast rotation, with periods of about 2 days or shorter, the coronal field is wrapped in the azimuthal direction as the azimuthal component of the wind velocity increases with respect to the radial component. For low base density, this tangling does not affect the mass loss rate much. However, for a strong magnetic field and high base density, the tangling leads to an increase in the coronal density due to the capture of dense plasma in large tangled coronal loops (case G is about an order of magnitude higher than case F near the equator). Second, the increase in magnetic field strength modifies the distribution of the stellar wind speed and, as a result, it also modifies the distribution of the coronal density. The combination of these dictates the mass loss rate. Figure 4 shows meridional slices of cases B and D colored with the contours of the stellar wind radial speed and number density, along with selected magnetic field lines. It can be seen that for a low magnetic field strength of 5 G, a fast, less dense stellar wind occupies the majority of the domain except for a narrow region near the equator. The low coronal density is due to the fact that most field lines are open so the stellar wind plasma can escape. In contrast, when the stellar field is strong, a large amount of the magnetic flux remains closed, and only weak and slow wind is developed at lower latitudes. Even at higher latitude, the wind does not exceed a speed of about 500 km s⁻¹, which is still considered "slow" and dense in terms of the fast solar wind. For stronger fields, then, the expansion of the magnetic flux tubes that carry the mass flux from the surface into space is greater than for the cases with weaker fields, and the overall wind speed is dominated by the slow, more dense plasma (Wang & Sheeley 1990; Arge & Pizzo 2000). As a result, the overall mass loss rate is greater as B increases.

Zoom In Zoom Out Reset image size

In order to quantify our results, we obtain a scaling law for the stellar mass loss rate, $\dot{M}_\star$, and the stellar angular momentum loss rate, $\dot{J}_\star$, as a function of the stellar base density, stellar dipole field, and stellar rotation period. These scaling laws are as follows:

$$\begin{equation} \frac{\dot{M}_\star }{\dot{M}_\odot }= K \left(\frac{n_\star }{n_\odot } \right)^\alpha \left(\frac{B_\star }{B_\odot } \right)^{\left(\frac{n_\odot }{n_\star } \right)^\beta } \left(\frac{P_\odot }{P} \right)^{\left(1-\frac{n_\odot }{n} \right)^\gamma } \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{\dot{J}_\star }{\dot{J}_\odot } = \left(\frac{P_\odot }{P_\star }\right) \left(\frac{\dot{M}_\star }{\dot{M}_\odot } \right). \end{equation} \tag{ 5 }$$

We use $\dot{M}_\odot =3\times 10^{-14}\;M_\odot \;{\rm yr}^{-1}$ and $\dot{J}_\odot =2\times 10^{29}\;{\rm g}\;{\rm cm}^2\;{\rm s}^{-2}$, which are the average solar mass loss rate and angular momentum loss rate between solar minimum and solar maximum periods, n_☉ = 2 × 10⁸ cm⁻³, B_☉ = 10 G, and P_☉ = 27 days. For lower values of the base density and for the solar case, the above scaling laws are reduced to

$$\begin{equation} \frac{\dot{M}_\star }{\dot{M}_\odot }= K \left(\frac{B_\star }{B_\odot } \right) \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{\dot{J}_\star }{\dot{J}_\odot } = \left(\frac{P_\odot }{P_\star }\right) K \left(\frac{B_\star }{B_\odot } \right). \end{equation} \tag{ 7 }$$

As stated in Section 2.2, based on our simulations of the Sun, we believe that the coronal base density should be at the lower end of the density range investigated here. Unlike the study of Matt et al. (2012), since we have only investigated solar mass stars, our scaling law does not depend on surface gravity nor the solar mass.

The scaling laws for rotation periods of 0.5 and 25 days are displayed in the plots of Figure 2 as dashed and dotted lines, respectively. The set of parameters that provides the best fit for all curves, simultaneously, is K = 3, α = 0.8, β = 0.2, and γ = 0.1. We also show a secondary best fit with K = 2, α = 0.8, β = 0.1, and γ = 0.1. In general, these equations describe the dependence of the stellar loss rates on the magnetic field strength, with a slope that is reduced with increasing density, and with the modification caused by rotation. The deviations from these lines are not more than a factor of two for most points, except for the cases with very high (and probably unrealistic) base density, or for the cases of very weak magnetic field, for which the Alfvén surface gets very close to the stellar surface.

In Figure 5, we directly compare our scaling law with the scaling law in Equation (9) from Matt et al. (2012). We also compare the two scaling laws for the magnetic field depending on rotation as B∝P⁻¹ or B∝P⁻². For each rotation period, we obtain the corresponding value of the magnetic field, and then obtain the torque from Equation (7) (assuming a base density of n_☉). We also use the corresponding magnetic field to estimate $\dot{M}$ from Equation (6), and use these values of B(P) and $\dot{M(B)}$ to estimate the torque via Equation (9) from Matt et al. (2012). We use the same values described in Matt et al. (2012) for K1, K2, and m, along with the solar rotation, Ω_⋆ = 2.7 × 10⁻⁶, the solar radius, R_⋆ = R_☉, and V_e = 600 km s⁻¹. Overall, the trends of the torques using the two scaling laws are quite similar. However, the overall torque in Matt et al. (2012) is about 12 times higher for a field strength of 10 G and 10 times higher for a field strength of 5 G. Looking at the equation for the torque as a function of the average Alfvén radius:

$$\begin{equation} \tau _w=\dot{M}\Omega _\star r^2_A, \end{equation} \tag{ 8 }$$

and adopting a solar mass loss rate of 2 × 10⁻¹⁴ M_☉ yr⁻¹ (1.8 × 10¹² g s⁻¹), and Ω_☉ = 2.7 × 10⁻⁶ s⁻¹, our predicted torque is obtained for r_A = 2–10 R_☉, while the torque predicted by Matt et al. (2012) is obtained for r_A = 15 R_☉, despite the fact that Table 1 in that paper predicts r_A to be between 8 and 10 R_☉ for the solar parameters of f = 0.004 and ϒ = 10²–10³. It is possible that the different in the torque between our scaling and that in Matt et al. (2012) is simply based on a different value assumed for the solar case.

Zoom In Zoom Out Reset image size

We note that the angular momentum loss rate, $\dot{J_\star }$ is well represented by a simple scaling of the mass loss rate proportionally to the rotation rate. This differs from the simple analytical relation developed for spherically symmetric mass loss that also includes a term involving the Alfvén radius, r_A, to some power depending on the field geometry (e.g., Weber & Davis 1967; Mestel 1984; Kawaler 1988; Wood et al. 2004; Reiners & Mohanty 2012). As we emphasized earlier, for our 3D models r_A is a more complex 3D entity, and more properly the average radius of a local wind density-weighted surface as defined in Equation (3). Nevertheless, what we find is that this Alfvén "radius" term is essentially constant in our angular momentum loss results, at least to the precision to which Equation (5) represents them. This might be somewhat surprising given the large range in magnetic field strength investigated, but it has an interesting implication. The magnetic field at the Alfvén radius, $B_{r_A}$ can be expressed in terms of the surface magnetic field, B₀, as

$$\begin{equation} B_{r_A}=B_0\left(\frac{R}{r_A}\right)^b, \end{equation} \tag{ 9 }$$

where the index b depends on the magnetic field configuration. Combining this with Equation (1), we can write

$$\begin{equation} u_A=\frac{B_0}{\sqrt{4\pi \rho _A}}\left(\frac{R}{r_A}\right)^b, \end{equation} \tag{ 10 }$$

where R is the stellar radius and ρ_A is the wind density at the Alfvén radius. Since in our models r_A is invariant, we find

$$\begin{equation} \rho _A u_A^2\propto B_0^2, \end{equation} \tag{ 11 }$$

that is, the kinetic energy density of the escaping wind responsible for the angular momentum loss is proportional to the magnetic field energy density. We also find that the wind speed does not vary strongly between models, with a maximum of about 800 km s⁻¹, similar to the solar wind speed and the escape velocity. The bulk of the kinetic energy variations between models are instead due to wind density variations, and so for the equatorial regions that dominate the angular momentum loss, we can write, very approximately, $\rho _A \propto B_0^2$. This simple dependence arises because, in order to escape through the magnetic field, the wind pressure must first overcome the magnetic pressure, and in this way the field pressure near the equator, B²/8π, acts as a density regulator.

CO09 have shown how active regions on the stellar surface can modify the mass loss and angular momentum loss rates. We should therefore ask ourselves if a dipole approximation is a good proxy for the stellar field. In order to answer this question, we perform a similar calculation of the mass loss and angular momentum loss rates but this time, we use actual solar magnetograms to drive the model. These Michelson Doppler Imager magnetograms¹ were obtained during a solar minimum period (Carrington Rotation 1922, 1997 May) where the solar magnetic field was nearly dipolar, and solar maximum period (Carrington Rotation 1962, 2000 April). Figure 6 shows a meridional cut colored with contours of local mass loss rate distribution for the solar minimum and solar maximum cases, as well as a case with a dipole field of 10 G. The polar field for the solar minimum case is roughly 8.5 G which is not far from that of the dipole case. The rotation period for all three runs is P = 25 days and the base density is n₀ = 2e8 cm⁻³.

Zoom In Zoom Out Reset image size

It can be clearly seen that the solar minimum and dipole cases are similar in orientation, but the Alfvén surface of the dipole case is much smaller than the one of the real Sun. The reason for this deviation is, again, the dependence of the wind structure on the flux-tube expansion. The expansion in reality in the polar regions seems to be much higher than in the dipole case. As a result, coronal polar regions in the dipole case are occupied by slow, more dense wind than the fast, less dense wind in the real case. It is also possible that for the real case, non-dipole components of the field contribute to the strength of the coronal magnetic field, leading to an increase in the size of the Alfvén surface.

Overall, the dipole case seems to overestimate the solar minimum mass loss rate, as well as the angular momentum loss rate by a factor of about three.

It is well established, on both theoretical and observational grounds, that the magnetic activity level of a star of a given mass is primarily dependent on the rotation period. We can derive the expected mass and angular momentum losses as a function of the rotation period if we can assume, or make reasonable guesses for, the values of the other two free parameters in our models, the magnetic field strength and base density. As noted in Section 2.2, the relevant magnetic field for driving the wind is the large-scale dipolar field, rather than the smaller scale and much stronger magnetic field associated with active regions and sunspots (see also Cohen 2011). From the perspective of understanding stellar winds, this is unfortunate because this dipolar, or poloidal, field is essentially never observed in late-type stars. Instead, stellar measurements of magnetic fields tend to be dominated by the closed field regions more analogous to solar active regions, where the magnetic fields are up to two orders of magnitude stronger (see also the discussion of Cranmer & Saar 2011). While previous studies have utilized these measurements of surface magnetic flux to infer wind properties and relations (e.g., Ivanova & Taam 2003; Schrijver et al. 2003; Holzwarth & Jardine 2007; Cranmer & Saar 2011), our model is somewhat different in that is it explicitly based on the dipolar field component that is expected to accompany, but be distinct from, the smaller-scale active regions. We therefore do not follow the same approach.

Lacking direct measurements, we identify the large-scale stellar dipole fields we have adopted here for our wind model driving with the poloidal field of an αΩ-type dynamo. (e.g., Durney & Robinson 1982) used elementary α − Ω dynamo theory to show that the poloidal magnetic field strength for an unsaturated dynamo should vary approximately linearly with the stellar angular velocity, B_p∝Ω. Although the true situation could be considerably more complicated, this relation offers a starting point. If we assume a linear Ω dependence of the dipole field strength, adopt a value of 10 G for a solar rotation period of 27 days, and fix the coronal base density at a value of 2 × 10⁸ cm⁻³, we can interpolate among our model solutions and use Equation (4) to derive the mass loss rate as a function of rotation period. This curve is illustrated in Figure 7, overplotted with our grid model results for the selected solar parameters. Also illustrated is a curve assuming a rotation dependence of B_p∝Ω² and the corresponding points from our model, together with the mass loss rate predicted for a solar mass star in the model of Cranmer & Saar (2011). The one remaining free parameter is the base density. Ivanova & Taam (2003) and Holzwarth & Jardine (2007) assessed evidence from X-ray luminosity versus rotation to relate the base plasma density to rotation rate adopting $n_\star \propto \Omega _\star ^{0.6}$. While, again, the X-ray emission underlying this assessment is dominated by closed field measurements, some increase in base density with activity level is probably not unreasonable. For the sake of comparison, we also illustrate in Figure 7 a curve corresponding to a case of base density increasing linearly with rotation rate, n_⋆∝Ω_⋆.

Zoom In Zoom Out Reset image size

We obtain mass loss rates quite similar to those predicted by Cranmer & Saar (2011) for a fixed base density and a field strength proportional to the square of the rotation velocity. Perhaps the most interesting region of Figure 7 from this perspective is for longer rotation periods than the Sun. Our model including a rotation-dependent density term, together with a linear dependence of magnetic field strength on rotation velocity, levels off toward a much higher asymptotic mass loss rate than the Cranmer & Saar (2011) model which continues to decline. This behavior arises in the density dependent terms in Equation (4). There are unfortunately no observations of mass loss rates for inactive stars with significantly longer rotation periods than the Sun with which to determine whether mass loss continues to decline strongly from the solar rate toward longer rotation periods, or declines much more slowly or levels off.

While our scaling law for the mass loss rate covers a range which is similar to that in Wood et al. (e.g., 2005), our mass loss rate drops much more quickly with rotation rate. The mass loss rate is still close to 100 times higher than solar for 1 Gyr (10–12 days period) in Wood et al. (2005) while it is less than 10 times the solar value in our results. This is quite significant in the context of the mass loss rate of the young Sun during the early evolution of the Earth and the faint young Sun paradox (see, e.g., Sagan & Mullen 1972; Kasting 1993). This difference between our results and the mass loss rate derived from scaling the stellar activity level strengthens the argument made by Cohen (2011) that only a weak dependence of the stellar mass loss on the observed stellar activity level is expected based on solar observations.

Once a relation between the magnetic field strength and the rotation period is known, the angular momentum loss relation in Equation (5) can also be used to examine the expected spin-down relation of main-sequence solar-like stars. Writing the stellar angular momentum as $J=kMR_\star ^2\Omega _\star$, where k is a constant depending on the density profile of the stellar interior, and assuming constant mass and radius, we can write $\dot{J}\propto \dot{\Omega _\star }$. Since $\dot{J}\propto B_p\Omega _\star$, and assuming $B_p\propto \Omega _\star ^a$, we have

$$\begin{equation} \frac{d\Omega _\star }{dt}=-\Omega ^{(a+1)}, \end{equation} \tag{ 12 }$$

implying Ω_⋆∝t^−1/a. We therefore retrieve the Skumanich (1972) spin-down relation for solar-type stars, Ω_⋆∝t^−1/2, for $B_p\propto \Omega ^2_\star$.