The Effect of Combined Magnetic Geometries on Thermally Driven Winds I: Interaction of Dipolar and Quadrupolar Fields
Cool stars with outer convective envelopes are observed to have magnetic fields with a variety of geometries, which on large scales are dominated by a combination of the lowest order fields such as the dipole, quadrupole and octupole modes. Magnetised stellar wind outflows are primarily responsible for the loss of angular momentum from these objects during the main sequence. Previous works have shown the reduced effectiveness of the stellar wind braking mechanism with increasingly complex, but singular, magnetic field geometries. In this paper, we quantify the impact of mixed dipolar and quadrupolar fields on the spin-down torque using 50 MHD simulations with mixed field, along with 10 of each pure geometries. The simulated winds include a wide range of magnetic field strength and reside in the slow-rotator regime. We find that the stellar wind braking torque from our combined geometry cases are well described by a broken power law behaviour, where the torque scaling with field strength can be predicted by the dipole component alone or the quadrupolar scaling utilising the total field strength. The simulation results can be scaled and apply to all main-sequence cool stars. For Solar parameters, the lowest order component of the field (dipole in this paper) is the most significant in determining the angular momentum loss.
Subject headings:magnetohydrodynamics (MHD) - stars: low-mass - stars: stellar winds, outflows - stars: magnetic field- stars: rotation, evolution
The spin down of cool stars () is a complex function of mass and age, as shown by the increasing number of rotation period measurements for large stellar populations (Barnes, 2003; Irwin & Bouvier, 2009; Barnes, 2010; Agüeros et al., 2011; Meibom et al., 2011; McQuillan et al., 2013; Bouvier et al., 2014; Stauffer et al., 2016; Davenport, 2017). Observed properties of these stars show a wide range of mass loss rates, coronal temperatures, field strengths and geometries, which all connect with stellar rotation to control the loss of angular momentum (Reiners & Mohanty, 2012; Gallet & Bouvier, 2013; Van Saders & Pinsonneault, 2013; Brown, 2014; Matt et al., 2015; Gallet & Bouvier, 2015; Amard et al., 2016; Blackman & Owen, 2016; See et al. in prep). Despite the wide range of interlinking stellar properties an overall trend of spin down with an approximately Skumanich law is observed at late ages; (Skumanich, 1972; Soderblom, 1983).
For Sun-like stars on the main sequence, the spin-down process is governed primarily by their magnetised stellar winds which remove angular momentum over the star’s lifetime. Parker (1958) originally posited that stellar winds must exist due to the thermodynamic pressure gradient between the high temperature corona and interplanetary space. Continued solar observations have constrained theoretical models for the solar wind to a high degree of accuracy (van der Holst et al., 2014; Usmanov et al., 2014; Oran et al., 2015). Recent models of the solar wind are beginning to accurately reproduce the energetics within the corona and explain the steady outflow of plasma into the Heliosphere (e.g. Grappin et al., 1983; Van der Holst et al., 2010; Pinto et al., 2016). The wind driving is now known to be much more complex than a thermal pressure gradient, with authors typically heating the wind through the dissipation of Alfvén waves in the corona. Other cool stars are observed with x-ray emissions indicating hot stellar coronae like that of the Sun (Rosner et al., 1985; Hall et al., 2007; Wright et al., 2004; Wolk et al., 2005). Similar stellar winds and wind heating mechanisms are therefore expected to exist across a range of Sun-like stars. Assuming equivalent mass loss mechanisms, results from the Solar wind are incorporated into more general stellar wind modelling efforts (e.g. Cohen & Drake, 2014; Alvarado-Gómez et al., 2016).
Detailed studies of wind driving physics remain computationally expensive to run, so are usually applied on a case-by-case basis. How applicabile the heating physics gained from modelling the Solar wind is to other stars still in question. With the reliability of such results even for the global properties of a given star in question, large parameter studies with simpler physics remain useful. A more general method can allow for parametrisations which are more appropriate to the variety of stellar masses and rotation periods found in observed stellar populations. Parker-type solutions remain useful for this due to their simplicity and versatility (Parker, 1965; Mestel, 1968; Sakurai, 1990; Keppens & Goedbloed, 1999). In these solutions, wind plasma is accelerated from the stellar surface and becomes transonic at the sonic surface. With the addition of magnetic fields the wind also become trans-alfvénic, i.e faster than the Alfvén speed, at the Alfvén surface. Weber & Davis (1967) showed for a one-dimensional magnetised wind that the Alfvén radius represented a lever arm for the spin-down torque. Since the introduction of this result, many researchers have produced scaling laws for the Alfvén radius (Mestel, 1984; Kawaler, 1988; Matt & Pudritz, 2008; Matt et al., 2012; Ud-Doula et al., 2009; Pinto et al., 2011; Réville et al., 2015a; Pantolmos. in prep) all of which highlight the importance of the magnetic field strength and mass loss rate in correctly parametrising a power law dependence. In such formulations, the mass loss rate is incourporated as a free parameter as the physical mechanisms which determines it are not yet completely understood. Measuring the mass loss rate from Sun-like stars is particularly difficult due to the wind’s tenuous nature and poor emission. Wood (2004) used Lyman- absorption from the interaction of stellar winds and their local interstellar medium to measure mass loss rates, but the method is model-dependent and only available for a few stars. Theoretical work from Cranmer & Saar (2011) predicts the mass loss rates from Sun-like stars, but it is uncertain if the physics used within the model scales correctly between stars. Therefore, parameter studies where the mass loss rate is an unknown parameter are needed.
In addition to the mass loss rate, the angular momentum loss rate is strongly linked with the magnetic properties of a given star. Frequently researchers assume the dipole component of the field to be the most significant in governing the global wind dynamics (e.g. Ustyugova et al., 2006; Zanni & Ferreira, 2009; Gallet & Bouvier, 2013; Cohen & Drake, 2014; Gallet & Bouvier, 2015; Matt et al., 2015; Johnstone et al., 2015). Zeeman Doppler Imaging (ZDI) studies (e.g. Morin et al., 2008; Petit et al., 2008; Fares et al., 2009; Vidotto et al., 2014b; Jeffers et al., 2014; See et al., 2015, 2016; Folsom et al., 2016; Hébrard et al., 2016; See et al., 2017), provide information on the large scale surface magnetic fields of active stars. Observations have shown stellar magnetic fields to be much more complex than simple dipoles, containing combinations of many different field modes. ZDI is a topographic technique typically decomposes the field at the stellar surface into individual spherical harmonic modes. The 3D field geometry can then be recovered with field extrapolation techniques using the ZDI map as an inner boundary. Several studies have considered how these observed fields affect the global wind properties. Typically used to determine an initial 3D field solution, then a magnetohydrodynamics code evolves this initial state in time until a steady state solution for the wind and magnetic field geometry is attained (e.g. Vidotto et al., 2011; Cohen et al., 2011; Garraffo et al., 2016b; Réville et al., 2016; Alvarado-Gómez et al., 2016; Nicholson et al., 2016; do Nascimento Jr et al., 2016). These works are less conducive to the production of semi-analytical formulations, as the principle drivers of the spin-down process are hidden within complex field geometries, rotation and wind heating physics.
A few studies show systematically how previous torque formulations depend on magnetic geometry using single modes. Réville et al. (2015a) explored thermally driven stellar winds with dipolar, quadrupolar and octupolar field geometries. They concluded that higher order field modes produce a weaker torque for the same field strength and mass loss, which is supported by results from Garraffo et al. (2016a). Despite these studies and works like them, only one study has systematic scaled the mass loss rate for a mixed field geometry field (Strugarek et al., 2014a). However, the aforementioned studies of the angular momentum loss from Sun-like stars have yet to address the systematic addition of individual spherical harmonic field modes.
Mixed geometry fields are observed within our closest star, the Sun, which undergoes a 11 year cycle oscillating between dipolar and quadrupolar field modes from cycle minimum to maximum respectively (DeRosa et al., 2012). Observed Sun-like stars also exhibit a range of spherical harmonic field combinations. Simple magnetic cycles are observed using ZDI, both HD 201091 (Saikia et al., 2016) and HD 78366 (Morgenthaler et al., 2012) show combinations of the dipole, quadrupole and octupole field modes oscillating similarly to the solar field. Other cool stars exist with seemingly stochastic changing field combinations (Petit et al., 2009; Morgenthaler et al., 2011). Observed magnetic geometries all contain combinations of different spherical harmonic modes with a continuous range of mixtures, it is unclear what impact this will have on the braking torque.
In this study we will investigate the significance of the dipole field when combined with a quadrupolar mode. We focus on these two field geometries, which are thought to contribute in anti-phase to the solar cycle and perhaps more generally to stellar cycles in cool stars. Section 2 covers the numerical setup with a small discussion of the magnetic geometries for which we develop stellar wind solutions. Section 3 presents the main simulation results, including discussion of the qualitative wind properties and field structure, along with quantitative parametrisations for the stellar wind torque. Here we also highlight the dipole’s importance in the braking, and introduce an approximate scaling relation for the torque. Finally in Section 4 we focus on the magnetic field in the stellar wind, first a discussion of the overall evolution of the flux, then a discussion of the open flux and opening radius within our simulations. Conclusions and thoughts for further work can then be found in Section 5. The Appendix contains a short note on the wind acceleration profiles of our wind solutions.
2. Simulation Method
2.1. Numerical Setup
This work uses the magnetohydrodynamics (MHD) code PLUTO (Mignone et al., 2007; Mignone, 2009), a finite-volume code which solves Riemann problems at cell boundaries in order to calculate the flux of conserved quantities through each cell. PLUTO is modular by design, capable of interchanging solvers and physics during setup. The present work uses a diffusive numerical scheme, the solver of Harten, Lax, and van Leer, HLL (Einfeldt, 1988), which allows for greater numerical stability in the higher strength magnetic field cases. The magnetic field solenoidality condition () is maintained using the Constrained Transport method (See Tóth (2000) for discussion).
The MHD equations are solved in a conservative form, with each equation relating to the conservation of mass, momentum and energy, plus the induction equation for magnetic field,
Here is the mass density, is the velocity field, is the gravitational acceleration, is the magnetic field
We assume the wind profiles to be axisymmetric and solve the MHD equations using a spherical geometry in 2.5D, i.e. our domain contains two spatial dimensions () but allows for 3D axisymetric solutions for the fluid flow and magnetic field using three vector components (). The domain extends from one stellar radius () out to with a uniform grid spacing in and a geometrically stretched grid in , which grows from an initial spacing of to at the outer boundary. The computational mesh contains grid cells. These choices allow for the highest resolution near the star, where we set the boundary conditions that govern the wind profile in the rest of the domain.
Initially a polytropic parker wind (Parker, 1965; Keppens & Goedbloed, 1999) with fills the domain, along with a super-imposed background field corresponding to our chosen magnetic geometry and strength. During the time-evolution, the plasma pressure, density, and poloidal components of the magnetic field () are held fixed at the stellar surface, whilst the poloidal components of the velocity () are allowed to evolve in response to the magnetic field (the boundary is held with and ). We then enforce the flow at the surface to be parallel to the magnetic field (). The star rotates as a solid body, with linearly extrapolated into the boundary and set using the stellar rotation rate ,
where the subscript “p” denotes the poloidal components () of a given vector. This condition enforces an effective rotation rate for the field lines which, in steady state ideal MHD, should be equal to the stellar rotation rate and conserved along field lines (Zanni & Ferreira, 2009; Réville et al., 2015a). This ensures the footpoints of the stellar magnetic field are correctly anchored into the surface of the star. The final boundary conditions are applied to the outer edges of the simulation, a simple outflow (zero derivative) is set at allowing for the outward transfer of mass, momenta and magnetic field, along with an axisymmetric condition along the rotation axis ( and ). Due to the supersonic flow properties at the outer boundary and its large radial extent compared with the location of the fast magnetosonic surface, any artefacts from the outer boundary cannot propagate upwind into the domain.
The code is run, following the MHD equations above, until a steady state solution is found. The magnetic fields modify the wind dynamics compared to the spherically symmetric initial state, with regions of high magnetic pressure shutting off the radial outflow. In this way, the applied boundary conditions allow for closed and open regions of flow to form (e.g Washimi & Shibata, 1993; Keppens & Goedbloed, 2000), as observed within the solar wind. In some cases of strong magnetic field small reconnection events are seen, caused by the numerical diffusivity of our chosen numerical scheme. Reconnection events are also seen in Pantolmos & Matt (in prep) and discussed within their Appendix. We adopt a similar method for deriving flow quantities in cases exhibiting periodic reconnection events. In such cases, once a quasi-steady state is established a temporal average of quantities such as the torque and mass loss are used.
Inputs for the simulations are given as ratios of characteristic speeds which control key parameters such as the wind temperature (), field strength () and rotation rate (). Where is the sound speed at the surface, is the Alfvén speed at the north pole, is the rotation speed at the equator, is the surface escape speed and is the keplerian speed at the equator. In this way, all simulations represent a family of solutions for stars with a range of gravities. As this work focuses on the systematic addition of dipolar and quadrupolar geometries, we fix the rotation rate for all our simulations. Matt et al. (2012) showed that the non-linear effects of rotation on their torque scaling can be neglected for slow rotators. They defined velocities as a fraction of the breakup speed,
The Alfvén radius remains independent of the stellar spin rate until , after which the effects of fast rotation start to be important. For this study a solar rotation rate is chosen (), which is well within the slow rotator regime. We set the temperature of the wind with , higher than used previosuly in Réville et al. (2015a). This choice of higher sound speed drives the wind to slightly higher terminal speeds, which are more consistent with observed solar wind speeds. Each geometry is studied with 10 different field strengths controlled by the input parameter , which is defined here with the Alfvén speed on the stellar north pole (see following Section). Table 1 lists all our variations of for each geometry.
Due to the use of characteristic speeds as simulation inputs, our results can be scaled to any stellar parameters. For example, using solar parameters, the wind is driven by a coronal temperature of 1.4MK and our parameter space covers a range of stellar magnetic field strengths from 0.9G to 87G over the pole. Changing these normalisations will modify this range.
2.2. Magnetic Field Configuration
Within this work, we consider magnetic field geometries that encompass a range of dipole and quadrupole combinations with different relative strengths. We represent the mixed fields using the ratio, , of dipolar field to the total combined field strength.
In this study the magnetic fields of the dipole and quadrupole are described in the formalism of Gregory et al. (2010) using polar field strengths,
The total field, comprised of the sum of the two geometries,
where the total polar field , is controlled by the parameter,
This work considers aligned magnetic moments such that ranges from 1 to 0, corresponding to all the field strength in the dipolar or quadrupolar mode respectively. As with , is calculated at the north pole. This sets the relative strengths of the dipole and quadrupole fields,
Alternative parametrisations are commonly used in the analysis of ZDI observations and dynamo modelling. These communities use the surface averaged field strengths, , or the ratio of magnetic energy density () stored within each of the dipole and quadrupole field modes at the stellar surface. During the solar magnetic cycle, values of can range from at solar maximum to at solar minimum (DeRosa et al., 2012). A transformation from our parameter to the ratio of energies is simply given by:
where the numerical pre-factor accounts for the integration of magnetic energy in each mode over the stellar surface.
Initial field configurations are displayed in Figure 1. The pure dipolar and quadrupolar cases are shown in comparison to two mixed cases (). These combined geometry fields add in one hemisphere and subtract in the other. This effect is due to the different symmetry families each geometry belongs to, with the dipole’s polarity reversing over the equator unlike the equatorially symmetric quadrupole. Continuing the use of “primary” and “secondary” families as in McFadden et al. (1991) and DeRosa et al. (2012), we refer to the dipole as primary and quadrupole as secondary. The fields are chosen such that they align in polarity in the northern hemisphere. This choice has no impact on the derived torque or mass loss rate due to the symmetry of the quadrupole about the equator. Either aligned or anti-aligned, these fields will always create one additive hemisphere and one subtracting; swapping their relative orientations simply switches the respective hemispheres. This is in contrast to combining dipole & octupole fields, where the aligned and anti-aligned cases cause subtraction at the equator or poles respectively (Gregory et al., 2016; Finley & Matt. in prep).
Figure 1 indicates that even with equal quadrupole and dipole polar field strengths, , the overall dipole topology will remain. In this case the magnetic energy density in the dipolar mode is 1.5 times greater than the quadrupolar mode and with the more rapid radial decay of the quadrupolar field, this explains the overall dipolar topology. A higher fraction of quadrupole is required to produce a noticeable deviation from this configuration, which is shown at . More than half of the parameter space that we explore lies in the range where the energy density of the quadrupole mode is greater than that of the dipole (). For this study both the pure dipolar and quadrupolar fields are used as controls (both of which were studied in detail within Réville et al. (2015a)), and 5 mixed cases parametrised by values ( = 0.8, 0.5, 0.3, 0.2, 0.1). We include to demonstrate the dominance of the dipole at higher values. Each value is given a unique identifying colour which is maintained in all figures throughout this paper. Table 1 contains a complete list of parameters for all cases, which are numbered by increasing and quadrupole fraction.
3. Simulation results
3.1. Morphology of the Field and Wind Outflow
Figure 1 shows the topological changes in field structure from the addition of dipole and quadrupole fields. It is evident in these initial magnetic field configurations that the global magnetic field becomes asymmetric about the equator for mixed cases, as does the magnetic boundary condition which is maintained fixed at the stellar surface. It is not immediately clear how this will impact the torque scaling from Réville et al. (2015a), who studied only single geometries.
Results for these field configurations using our PLUTO simulations are displayed in Figure 2. The dipole and quadrupole cases are shown in conjunction with the mixed field cases, . The Figure displays for a comparable value of polar magnetic field strength, the different sizes of Alfvén surface that are produced. The mixed magnetic geometries modify the size and morphology of the Alfvén and sonic surfaces. Due to the slow rotation, the fast and slow magnetosonic surfaces are co-located with the sonic and Alfvén surfaces (the fast magnetosonic surface being always the larger of the two surfaces).
The field geometry is found to imprint itself onto the stellar wind velocity with regions of closed magnetic field confining the flow creating areas of co-rotating plasma, referred to as deadzones (Mestel, 1968). Steady state wind solutions typically have regions of open field where a faster wind and most of the torque is contained, along with these deadzone(s) around which a slower wind is produced. Similarly to the solar wind, slower wind can be found on the open field lines near the boundary of closed field (Feldman et al., 2005; Riley et al., 2006; Fisk et al., 1998). Observations of the Sun reveal the fast wind component emerging from deep within coronal holes, typically over the poles, and the slow wind component originating from the boundary between coronal holes and close field regions. Due to the polytropic wind used here, we do not capture the different heating and acceleration mechanisms required to create a true fast and slow solar-like wind (as seen with the Ulysses spacecraft e.g. McComas et al., 2000; Ebert et al., 2009). Our models produce an overall wind speed consistent with slow solar wind component, which we assume to represent the average global flow. More complex wind driving and coronal heating physics are required to recover a multi-speed wind, as observed from the Sun (Cranmer et al., 2007; Pinto et al., 2016).
Figure 3 displays a grid of simulations with a range of magnetic field strengths and values ( ranges from 3.6 to 54; values consistent with the solar cycle maximum), where the mixing of the fields plays a clear role in the changing dynamics of the flow. Regions of closed magnetic field cause significant changes to the morphology of the wind. A single deadzone is established on the equator by the dipole geometry whereas the quadrupole creates two over mid latitudes. Mixed cases have intermediate states between the pure regimes. Within our simulations the deadzones are accompanied by streamers which form above closed field regions and drive slower speed wind than from the open field regions. The dynamics of these streamers, their location and size are an interesting result of the changing topology of the flow.
The dashed coloured lines within Figure 3 show where the field polarity reverses using , which traces the location of the streamers. The motion of the streamers through the grid of simulations is then observed. With increasing quadrupole field, the single dipolar streamer moves into the northern hemisphere and with continued quadrupole addition a second streamer appears from the southern pole and travels towards the northern hemisphere until the quadrupolar streamers are recovered both sitting at mid latitudes. This motion can also be seen for fixed cases as the magnetic field strength is decreased. For a given value the current sheets sweep towards the southern hemisphere with increased polar field strength, in some cases (36 and 38) moving onto the axis of rotation. This is the opposite behaviour to decreasing the value, i.e. the streamer configuration is seen to take a more dipolar morphology as the field strength is increased. Additionally within Figure 3, for low field strengths each produces a comparable Alfvén surface with very similar morphology, all dominated by the quadrupolar mode.
3.2. Global Flow Quantities
Our simulations produce steady state solutions for the density, velocity and magnetic field structure. To compute the wind torque on the star we calculate , a quantity related directly to the angular momentum flux (Keppens & Goedbloed, 2000),
Within axisymmetric steady state ideal MHD, is conserved along any given field line. However we find variations from this along the open-closed field boundary due to numerical diffusion across the sharp transition in quantities found there. The spin-down torque, , due to the transfer of angular momentum in the wind is then given by the area integral,
where is the area of any surface enclosing the star. For illustrative purposes, Figure 3 shows the Alfvén surface coloured by angular momentum flux (thick multi-coloured line), which is seen to be strongly focused around the equatorial region. The angular momentum flux is calculated normal to the Alfvén surface,
where is the normal unit vector to the Alfvén surface. The mass loss rate from our wind solutions is calculated similarly to the torque,
Both expressions for the mass loss and torque are evaluated using spherical shells of area which are outside the closed field regions. This allows for the calculation of an average Alfvén radius (which is cylindrical from the rotation axis) in terms of the torque, mass flux and rotation rate,
Throughout this work, is used as a normalised torque which accounts for the mass loss rates which we do not control. Values of the average Alfvén radius are tabulated within Table 1. is shown in Figure 3 using a grey vertical dashed line. For each case, the cylindrical Alfvén radius is offset inwards of the maximum Alfvén radius from the simulation, a geometrical effect as this corresponds to the average cylindrical and includes variations in flow quantities as well. Exploring Figure 3, the motion of the deadzones/current sheets have little impact on the overall torque. For example, no abrupt increase in the Alfvén radius is seen from case 34 to 36 (where the southern streamer is forced onto the rotation axis) compared to cases 44 and 46. The torque is instead governed by the magnetic field strength in the wind which controls the location of the Alfvén surface.
We parametrise the magnetic and mass loss properties using the “wind magnetisation” defined by,
where is the combined field strength at the pole. Previous studies that used this parameter defined it with the equatorial field strength (e.g. Matt & Pudritz, 2008; Matt et al., 2012; Réville et al., 2015a; Pantolmos & Matt. in prep). We use polar values unlike previous authors due to the additive property of the radial field at the pole, for aligned axisymmetric fields. Note that selecting one value of the field on the surface will not always produce a value which describes the field as a whole. The polar strength works for these aligned fields, but will easily break down for un-aligned fields and anti-aligned axisymmetric odd fields, thus it suits the present study, but a move away from this parameter in future is warranted.
During analysis, the wind magnetisation, , is treated as an independent parameter that determines the Alfvén radius and thus the torque, . We increase by setting a larger , creating a stronger global magnetic field. Table 1 displays all the input values of and as well as the resulting global outflow properties from our steady state solutions, which are used to formulate the torque scaling relations within this study. Figure 4 displays all 70 simulations in space. Cases are colour-coded here by their value, a convention which is continued throughout this work.
3.3. Single Mode Torque Scalings
The efficiency of the magnetic braking mechanism is known to be dependent on the magnetic field geometry. This has been previsously shown for single mode geometries (e.g. Réville et al., 2015a; Garraffo et al., 2016a). We first concider two pure gemetries, dipole and quadrupole, using the formulation from Matt & Pudritz (2008),
where and are fitting parameters for the pure dipole and quadrupole cases, using the surface field strength. Here we empirically fit ; the interpretation of is discussed in Matt & Pudritz (2008), Réville et al. (2015a) and Pantolmos & Matt (in prep), where it is determined to be dependant on magnetic geometry and the wind acceleration profile. The Appendix contains further discussion of the wind acceleration profile and its impact on this power law relationship.
The left panel of Figure 5 shows the Alfvén radii vs the wind magnetisations for all cases (colour-coded with their value). Solid lines show scaling relations for dipolar (red) and quadrupolar (blue) geometries, as first shown in Réville et al. (2015a). We calculate best fit values for and for the dipole and quadrupole, tabulated in Table 2. Values here differ due to our hotter wind ( than their ), using polar , and we do not account for our low rotation rate. As previously shown, the dipole field is far more efficient at transferring angular momentum than the quadrupole. In this study we concider the effect of combined geometries, within Figure 5 these cases lie between the dipole and quadrupole slopes, with no single power law of this form to describe them.
Pantolmos & Matt (in prep) have shown the role of the velocity profile in the power law dependence of the torque. In our simulations, the acceleration of the flow from the base wind velocity to its terminal speed is primarily governed by the thermal pressure gradient, however magnetic topologies can all modify the radial velocity profile (as can changes in wind temperature, , and rapid rotation, not included in our study). Effects on the torque formulations due to these differences in acceleration can be removed via the multiplication of with . In their work, the authors determine the theoretical power law dependence, , from one-dimensional analysis. In this formulation the slope of the power law is controlled only by the order of the magnetic geometry, , which is and for the dipole and quadrupole respectively,
where and are fit parameters to our wind solutions, tabulated in Table 2. The value of is calculated as an average of the velocity at all points on the Alfvén surface in the meridional plane.
Equation (22) is able to predict accurately the power law dependence for the two pure modes using the order of the spherical harmonic field, . We show this in the right panel of Figure 5, where the Alfvén radii are plotted against the new parameter, . A similar qualitative behaviour is shown to the scaling with in the left panel. Using the theoretical power law dependencies, the dipolar (red) and quadrupolar (blue) slopes are plotted with and respectively. Using a single fit constant for both sloes within this figure shows good agreement with the simulation results.
More accurate values of and are fit for each mode independently. These values produce a better fit and are compared with the theoretical values in Table 2. The mixed simulations show a similar qualitative behaviour to the plot against .
Obvious trends are seen within the mixed case scatter. A saturation to quadrupolar Alfvén radii values for lower and values is observed, along with a power law trend with a dipolar gradient for higher and values. This indicates that both geometries play a role in governing the lever arm, with the dipole dominating the braking process at higher wind magnetisations.
3.4. Broken Power Law Scaling For Mixed Field Cases
Observationally the field geometries of cool stars are, at large scales, dominated by the dipole mode with higher order modes playing smaller roles in shaping the global field. It is the global field which controls the spin-down torque in the magnetic braking process. Higher order modes (such as the quadrupole) decay radially much faster than the dipole and as such they have a reduced contribution to setting the Alfvén speed at distances larger than a few stellar radii.
We calculate , which only takes into account the dipole’s field strength,
Taking as a hypothesis that the field controlling the location of the Alfvén radius is the dipole component, a power law scaling using can be constructed in the same form as Matt & Pudritz (2008),
Substitution of the dipole component into equation (22) similarly gives,
where , , , and will be parameters fit to simulations.
A comparison of these approximations can be seen in Figure 5, where equations (24) (left panel) and (25) (right panel) are plotted with dashed lines for all the values used in our simulations. Mixed cases which lie above the quadrupolar slope are shown to agree with the dashed-lines in both forms. Such cases are dominated by the dipole component of the field only, irrespective of the quadrupolar component.
The role of the dipole is even more clear in Figure 6 where only the dipole component of is plotted for each simulation. The solid red line in Figure 6, given by equation (24), shows agreement at a given with deviation from this caused by a regime change onto the quadrupolar slope (shown in dashed colour).
The behaviour of our simulated winds, despite using a combination of field geometries, simply follow existing scaling relations with this modification. In general, the dipole () prediction shows good agreement with the simulated wind models, except in cases where the Alfvén surface is close-in to the star. In these cases, the quadrupole mode still has magnetic field strength able to control the location of the Alfvén surface. Interestingly, and in contrast to the dipole-dominated regime, the quadrupole dominated regime behaves as if all the field strength is within the quadrupolar mode. This is visible within Figure 5 for low values of and .
The mixed field scaling can be described as a broken power law, set by the maximum of either the dipole component or the pure quadrupolar relation. With the break in the power law given by ,
where is the location of the intercept for the dipole component and pure quadrupole scalings,
The solid lines in Figure 4 show the value of , equation (27), diving the two regimes. Specifically, the solutions above the solid black line behave as if only the dipole component () is governing the Alfvén radius.
Transitioning from regimes is not perfectly abrupt. Therefore producing an analytical solution for the mixed cases which includes this behaviour would increase the accuracy for stars near the regime change. E.g. we have formulated a slightly better fit, using a relationship based on the quadrature addition of different regions of field. However it provides no reduction to the error on this simpler form and is not easily generalised to higher topologies. For practical purposes, the scaling of equation (26) and (27) predict accurately the simulation torque with increasing magnetic field strength for a variety of dipole fractions. We therefore present the simplest available solution, leaving the generalised form to be developed within future work.
4. The impact of geometry on the magnetic flux in the wind
4.1. Evolution of the Flux
The magnetic flux in the wind is a useful diagnostic tool. The rate of the stellar flux decay with distance is controlled by the overall magnetic geometry. We calculate the magnetic flux as a function of radial distance by evaluating the integral of the magnetic field threading closed spherical shells, where we take the absolute value of the flux to avoid field polarity cancellations,
Considering the initial potential fields of the two pure modes this is simply a power law in field order ,
where dipole and quadrupole, we denote the flux with “” for the potential field. Figure 7 displays the flux decay of all values of for each value, grey lines. The behaviour is qualitatively identical to that observed within previous works (e.g. Schrijver et al., 2003; Johnstone et al., 2010; Vidotto et al., 2014a; Réville et al., 2015a), where the field decays as the potential field does until the pressure of the wind forces the field into a purely radial configuration with a constant magnetic flux, referred to as the open flux. The power law dependence of equation (29) indicates for higher mode magnetic fields, the decay will be faster. We therefore expect the more quadrupolar dominated fields studied in this work to have less open flux.
In the case of mixed geometries a simple power law is not available for the initial potential configurations, instead we evaluate the flux using equation (28), where is the initial potential field for each mixed geometry. This allows us to calculate the radial evolution of the flux for a given which we compare to the simulated cases. Figure 7 shows the flux normalised by the surface flux versus radial distance from the star. For each value, the magnetic flux decay of the potential field (black solid line) is shown with the different strength simulations (grey solid lines). A comparison of the flux decay for all potential magnetic geometries is available in the bottom right panel showing, as expected, the increasingly quadrupolar fields decaying faster.
In this study we control which, for a given surface density, sets the polar magnetic field strength for our simulations. The stellar flux for different topologies and the same will differ and must be taken into account in order to describe the dipole and quadrupolar components (dashed red and blue) in Figure 7. We plot the magnetic flux of the potential field quadrupole component alone in dotted blue for each value,
and similarly the potential field dipole component of the magnetic flux,
where in both equations the surface flux of a pure dipole/quadrupole (, ) field is required to match our normalised flux representation.
Due to the rapid decay of the quadrupolar mode, the flux at large radial distances for all simulations containing the dipole mode is described by the dipolar component. The quadrupole component decay sits below and parallel to the potential field prediction for small radii, becoming indistinguishable for the lowest values as the flux stored in the dipole is decreased. Importantly for small radii, simulations containing a quadrupolar component are dominated by the quadrupolar decay following a power law decay, which can be seen by shifting the blue dashed line upwards to intercept at the stellar surface.
This result for the flux decay is reminiscent of the broken power law description for the Alfvén radius in Section 3.4. The field acts as a quadrupole using the total field for small radii and the dipole component only for large radii. There is a transition between these two regimes that is not described by either approximation. But is shown by the potential solution in solid black.
4.2. Topology Independent Open Flux Formulation
The magnetic flux within the wind decays following the potential field solution closely until the magnetic field geometry is opened by the pressures of the stellar wind and the field lines are forced into a nearly radial configuration with constant flux, shown in Figure 7 for all simulations. The importance of this open flux is discussed by Réville et al. (2015a). These authors showed a single power law dependence for the Alfvén radius, independent of magnetic geometry, when parametrised in terms of the open flux, ,
which, ignoring the effects of rapid rotation, can be fit with,
where, and are fitting parameters for the open flux formulation.
Using the open flux parameter, Figure 8 shows a collapse towards a single power law dependence as in Réville et al. (2015a). However our wind solutions show a systematic difference in power law dependence from dipole to quadrupole. On careful inspection of the result from Figure 6 of Réville et al. (2015a), the same systematic trend between their topologies and the fit scaling is seen.
Pantolmos & Matt (in prep) find solutions for thermally driven winds with different coronal temperatures, from these they find the wind acceleration profiles of a given wind to very significantly alter the slope in - space. From this work our trend with geometry indicates that each geometry must have a slightly different wind acceleration profile. This is most likely due to difference in the super radial expansion of the flux tubes for each geometry, which is not taken into account with equation (33). The field geometry is imprinted onto the wind as it accelerates out to the Alfvén surface. As such, this scaling relation is not entirely independent of topology. Further details on the wind acceleration profile within our study is available in the Appendix. Pantolmos (in prep) are able to include the effects of acceleration in their scaling through multiplication of with . The expected semi-analytic solution from Pantolmos & Matt (in prep) is given,
where the fit parameters are derived from one-dimensional theory as constants, and .
We are able to reproduce this power law fit of with the wind acceleration effects removed, on the right panel of Figure 8. Including all simulations in the fit, we arrive at values of and for the constants of proportionality and power law dependence. However a systematic difference is still seem from one value to another. More precise fits can be found for each geometry independently, but the systematic difference appearing in the right panel implies a modification to our semi-analytic formulations is required to describe the torque fully in terms of the open flux.
Here we show the scaling law from Réville et al. (2015a) is improved with the modification from Pantlomos (in prep). This formulation is able to describe the Alfvén radius scaling with changing open flux and mass loss. However with the open flux remaining an unknown from observations and difficult to predict, scaling laws that incorporate known parameters (such as those of equations (26) and (27)) are still needed for rotational evolution calculations.
4.3. The Relationship Between the Opening and Alfvén Radii
The location of the field opening is an important distance. It is both critical for determining the torque and for comparison to potential field source surface (PFSS) models (Altschuler & Newkirk, 1969), which set the open flux with a tunable free parameter . The opening radius, , we define is the radial distance at which the potential flux reaches the value of the open flux (). This definition is chosen because it relates to the 1D analysis employed to describe the power law dependences of our torque scaling relations. Specifically, a known value of allows for a precise calculation of the open flux (a priori from the potential field equations), which then gives the torque on the star within our simulations. The physical opening of the simulation field takes place at slightly larger radii than this with the field becoming non-potential due to its interaction with the wind (which explains why the closed field regions seen in Figure 3 typically extend slightly beyond ). A similar smooth transition is produced with PFSS modelling.â
is marked for each simulation in Figure 7 and again for comparative purposes in the bottom right panel. It is clear that smaller opening radii are found for lower cases. Due to their more rapidly decaying flux, they tend to have a smaller fraction of the stellar flux remaining in the open flux. From the radial decay of the magnetic field, the open flux and opening radii are observed to be dependent on the available stellar flux and topology. Pantolmos & Matt (in prep) have recently shown these to also be dependent on the wind acceleration profile. This complex dependence makes it difficult to predict the open flux for a given system.
Our simulations produce values for the average Alfvén radius, , and the opening radius, , for the 7 different geometries studied. It is interesting to consider the relative size of these radii as they both characterise key dynamic properties for each stellar wind solution. For all cases shown in Figure 3, the opening radii are plotted in dashed red, allowing for the relative size to be compared with the cylindrical Alfvén radius, shown in dash grey. With increasing magnetic field strength (), both radii are seen to grow from case to case, however with increasing , the cylindrical Alfvén radius generally grows faster than the opening radius. To quantify this, Figure 9 shows a plot of the Alfvén radii vs the opening radii for all cases. Linear trends of and are indicated with dashed lines. For each value, the relationship between the Alfvén and opening radius () is seen to systematically decrease with increasing higher order field component. In all cases, for small radii a shallower slope is observed which then steepens with increasing radial extent.
The dependence of the Alfvén radius and opening radius on field geometry and magnetisation is a constraint on PFSS models, which are readily used with ZDI observations as a less computationally expensive alternative to MHD modelling (Jardine et al., 1999, 2002; Dunstone et al., 2008; Cohen et al., 2010; Johnstone et al., 2010; Rosén et al., 2015; Réville et al., 2015b). PFSS models are a useful tool, however require the source surface radius, , as an input. Authors often set a source surface and change the geometry and strength of the field freely (Fares et al., 2010; See et al., 2015, 2017). We find however for a given value there exists a differing relation for the opening radius, as we define it here, to the Alfvén radius and magnetisation. These trends are observed to continue for higher mode fields (Finley & Matt. in prep), with decreasing overall with increased field complexity. As such, our results confirm that the opening radius should not remain fixed when changing geometries or increasing the wind magnetisation. We find the relationship of to change in both cases. With fixed magnetisation, the opening radius should move towards the star for higher order fields to maintain a constant thermal driving. Maintaining the opening radius whilst increasing the field complexity infers that the wind has a reduced acceleration. Similarly with increased wind magnetisation the opening radius should move further from the star. The value of as we have defined it, is directly related to the source surface radius, and for a given magnetic geometry, the two should scale approximately together. For example, for a dipole field, comparing our definition of to the PFSS model shows that equals an approximately constant value of 3/2 . Thus conclusions made about the opening radii are constraints on future PFSS modelling.
A method for predicting within our simulations remains unknown, however it is understood that is key to predicting the torque from our simulated winds. We do however find the ratio of to be roughly constant for a given geometry, deviations from which may be numerical or suggest additional physics which we do not explore here.
We undertake a systematic study of the two simplest magnetic geometries, dipolar and quadrupolar, and for the first time their combinations with varying relative strengths. We parametrise the study using the ratio, , of dipolar to total combined field strength, which is shown to be a key variable in our new torque formulation.
We have shown that a large proportion of the magnetic field energy needs to be in the quadrupole for any significant morphology changes to be seen in the wind. All cases above 50% dipole field show a single streamer and are dominated by dipolar behaviour. Even in cases of small we observe the dipole field to be the key parameter controlling the morphology of the flow, with the quadrupolar field rapidly decaying away for most cases leaving the dipole component behind. For smaller field strengths the Alfvén radii appears close to the star, where the quadrupolar field is still dominant, and thus a quadrupolar morphology is established. Increasing the fraction of quadrupolar field strength allows this behaviour to continue for larger Alfvén radii.
The morphology of the wind can be concidered in the context of star-planet or disk interactions. Our findings suggest that the connectivity, polarity and strength of the field within the orbital plane depend in a simple way on the relative combination of dipole and quadrupole fields. Different combinations of these two field modes change the location of the current sheet(s) and the relative orientation of the stellar wind magnetic field with respect to any planetary or disk magnetic field. Asymmetries such as these can modify the poynting flux exchange for close-in planets (Strugarek et al., 2014b) or the strength of magnetospheric driving and geomagnetic storms on Earth-like exoplanets. Cohen et al. (2014) use observed magnetic fields to simulate the stellar wind environment surrounding the planet hosting star EV Lac. They calculate the magnetospheric joule heating on the exoplanets orbiting the M dwarf, finding significant changes to atmospheric properties such as thickness and temperature. Additionally, transient phenomena in the Solar wind such as coronal mass ejections are shown to deflect towards streamer belts (Kay et al., 2013). This has been applied to mass ejections around M dwarfs stars (Kay & Opher, 2014), and could similarly be applied here using knowledge of the streamer locations from our model grid.
If the host star magnetic field can be observed and decomposed into constituent field modes, containing dominant dipole and quadrupole components, a qualitative assessment of the stellar wind environment can be made. We find the addition of these primary and secondary fields to create an asymmetry which may shift potentially habitable exoplanets in and out of volatile wind streams. Observed planet hosting stars such as Bootis have already been shown to have global magnetic fields which are dominated by combinations of these low order field geometries (Donati et al., 2008). With further investigation it is possible to qualitatively approximate the conditions for planets in orbit of such stars. For dipole and quadrupole dominated host stars with a given magnetic field strength our grid of models provide an estimate of the location of the streamers and open field regions.
Within this work we build on the scaling relations from, Matt et al. (2012), Réville et al. (2015a) and Pantolmos & Matt (in prep). We confirm existing scaling laws and explore a new mixed field parameter space with similar methods. From our wind solutions we fit the variables, , , and (see Table 4), which describe the torque scaling for the pure dipole and quadrupole modes. From the 50 mixed case simulations, we produce an approximate scaling relation which takes the form of a broken power law, as a single power law fit is not available for the mixed geometries cases in space.
For low and low dipole fraction, the Alfvén radius behaves like a pure quadrupole,
At higher and dipole fractions, the torque is only dependent on the dipolar component of the field,
The later formulation is used when the Alfvén radius of a given dipole & quadrupole mixed field is greater than the pure quadrupole case for the same , i.e. the maximum of our new formula or the pure quadrupole. We define to separate the two regimes (see Figure 4).
The importance of the relative radial decay of both modes and the location of the opening and Alfvén radii appear to play a key role, and deserve further follow up investigation. This work analytically fits the decay of the magnetic flux, but a parametric relationship for the field opening remains uncertain. The relation of the relative sizes of the Alfvén and opening radii are found to be dependent on geometry, which can be used to inform potential field source surface modelling, where by the source surface must be specified when changing the field geometry.
Paper II includes the addition of octupolar field geometries, another primary symmetry family which introduces an additional complication in the relative orientation of the octupole to the dipole. It is shown however, that the mixing of any two axisymmetric geometries will follow a similar behaviour, especially if each belongs to different symmetry families (Finley & Matt. in prep). The lowest order mode largely dominates the dynamics of the torque until the Alfvén radii and opening radii are sufficiently close to the star for the higher order modes to impact the field strength.
Appendix A Wind Acceleration
The creation of a semi-analytic formulation for the Alfvén radius for a variety of stellar parameters has been the goal of many studies proceeding this (e.g. Matt & Pudritz, 2008; Matt et al., 2012; Réville et al., 2015a; Pantolmos & Matt. in prep). Using a one-dimensional approximation based on work by Kawaler (1988), previous studies have aimed to predict the power law dependence, , of the torque formulations used within this work.
Using the one-dimensional framework, the field strength is assumed to decay as a power law , which in this study is only valid for the pure cases. Pantolmos & Matt (in prep) show the effect of wind acceleration can be removed from the torque scaling relations through the multiplication of and with . The power law dependences then becomes,
The modified dependent parameter, , is used throughout this work (see Figures 5 and 8), and the analytic predictions for the power law slopes are shown to have good agreement with our simulations. This dependent variable however, requires additional information about the wind speed at the Alfvén surface which is often unavailable.
Typically, rotation evolution models use the available stellar surface parameters e.g. . Therfore knowledge of the flow speed at the Alfvén radius, , is required for the semi-analytic formulations. is shown by Pantolmos & Matt. (in prep) and Réville et al. (2015a) to share a similar profile to a one-dimensional thermal wind, . Figure 10 displays the average Alfvén speed vs the Alfvén radius for all 70 simulations (coloured points). The parker wind solution (Parker, 1965) used in the initial condition is displayed for comparison (dashed line). Nearly all simulations follow the hydrodynamic solution, with a behaviour mostly independent of . Towards higher values of the Alfvén radius, a noticeable separation starts to develop between geometries. This range is accessed less by the higher order geometries as the range of Alfvén radii is much smaller than that for the pure dipole mode.
In order to include the effects of wind acceleration in the simplified 1D analysis to explain the simulation scalings between and , Réville et al. (2015a) introduced a parametrisation for the acceleration of the wind to the Alfvén radius with a power law dependence in radial distance using ,
A single power law with is fit to the simulation data, which is chosen for simplicity within the 1D formalism. The use of this parameter is approximate if is a power law in , which we show over the parameter space has a significant deviation. Using the semi-analytic theory, Réville et al. (2015a) then derived the power law dependence for the scaling (Equation (21)),
which includes geometric and wind acceleration parameters in the form of and respectively. Using this result, is computed for both the dipole () and quadrupole () geometries in Table 4, and compared to the simulation results with good agreement.
Pantolmos & Matt. (in prep) explain the power-law dependence, so long as remains constant and the wind acceleration profile is known. Reiners & Mohanty (2012), Réville et al. (2015a) and Pantolmos & Matt. (in prep) all analytically describe the power law dependence of the open flux formulation (Equation (33)) using the power law dependence ,
The result is independent of geometry, . As before the parameter approximates the wind driving as a power law in radius, which is fit with a single power law for both geometries such that should be the same for both the dipole and quadrupole. This prediction is tabulated in Table 4, however the simulation slopes are shown to no longer agree with the result. It is suggested that the open flux slope is much more sensitive to the wind acceleration than the formulation, therefore slight changes in flow acceleration modify the result. Slightly different slopes can be fit for the dipole and quadrupole cases which can recover the different values, however this is seemingly just a symptom of the power law approximation breaking down.
We conclude that the approximate power law of equation (A3) give a reasonable adjustment to the torque prediction for known wind velocity profiles, despite the badness of fit to the simulations points. Even though the power-law approximation to the wind velocity profile (equation A3) is not a precise fit to the data in Figure 10, the value of does provide a way to approximately include the contribution of the wind acceleration to the fit power-law exponents and . A more precise formulation could be derived based on a Parker-like wind profile without the use of a power law, however the torque scaling with is relative insensitive to the chosen approximate velocity profile.
- slugcomment: Draft: March 6, 2018
- The PLUTO code operates with a factor of absorbed into the normalisation of B. Tabulated parameters are given in cgs units with this factor incorporated.
- It could be argued that this should be weighted by the total area of the Alfvén surface, but for simplicity we calculate the un-weighted average.
- A choice in our parameter space may have made this clearer to see in Figure 8, due to the increased heating and therefore larger range of acceleration allowing the topology to impact the velocity profile.
- Agüeros, M. A., Covey, K. R., Lemonias, J. J., et al. 2011, The Astrophysical Journal, 740, 110
- Altschuler, M. D., & Newkirk, G. 1969, Solar Physics, 9, 131
- Alvarado-Gómez, J., Hussain, G., Cohen, O., et al. 2016, Astronomy & Astrophysics, 594, A95
- Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, Astronomy & Astrophysics, 587, A105
- Barnes, S. A. 2003, The Astrophysical Journal, 586, 464
- —. 2010, The Astrophysical Journal, 722, 222
- Blackman, E. G., & Owen, J. E. 2016, Monthly Notices of the Royal Astronomical Society, 458, 1548
- Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, Protostars and Planets VI, 433
- Brown, T. M. 2014, The Astrophysical Journal, 789, 101
- Cohen, O., Drake, J., Glocer, A., et al. 2014, The Astrophysical Journal, 790, 57
- Cohen, O., Drake, J., Kashyap, V., Hussain, G., & Gombosi, T. 2010, The Astrophysical Journal, 721, 80
- Cohen, O., & Drake, J. J. 2014, The Astrophysical Journal, 783, 55
- Cohen, O., Kashyap, V., Drake, J., et al. 2011, The Astrophysical Journal, 733, 67
- Cranmer, S. R., & Saar, S. H. 2011, The Astrophysical Journal, 741, 54
- Cranmer, S. R., Van Ballegooijen, A. A., & Edgar, R. J. 2007, The Astrophysical Journal Supplement Series, 171, 520
- Davenport, J. R. A. 2017, ApJ, 835, 16
- DeRosa, M., Brun, A., & Hoeksema, J. 2012, The Astrophysical Journal, 757, 96
- do Nascimento Jr, J.-D., Vidotto, A., Petit, P., et al. 2016, The Astrophysical Journal Letters, 820, L15
- Donati, J.-F., Moutou, C., Fares, R., et al. 2008, Monthly Notices of the Royal Astronomical Society, 385, 1179
- Dunstone, N., Hussain, G., Collier Cameron, A., et al. 2008, Monthly Notices of the Royal Astronomical Society, 387, 481
- Ebert, R., McComas, D., Elliott, H., Forsyth, R., & Gosling, J. 2009, Journal of Geophysical Research: Space Physics, 114
- Einfeldt, B. 1988, SIAM Journal on Numerical Analysis, 25, 294
- Fares, R., Donati, J.-F., Moutou, C., et al. 2009, Monthly Notices of the Royal Astronomical Society, 398, 1383
- —. 2010, Monthly Notices of the Royal Astronomical Society, 406, 409
- Feldman, U., Landi, E., & Schwadron, N. 2005, Journal of Geophysical Research: Space Physics, 110
- Fisk, L., Schwadron, N., & Zurbuchen, T. 1998, Space Science Reviews, 86, 51
- Folsom, C. P., Petit, P., Bouvier, J., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 580
- Gallet, F., & Bouvier, J. 2013, Astronomy & Astrophysics, 556, A36
- —. 2015, Astronomy & Astrophysics, 577, A98
- Garraffo, C., Drake, J. J., & Cohen, O. 2016a, Astronomy & Astrophysics, 595, A110
- —. 2016b, The Astrophysical Journal Letters, 833, L4
- Grappin, R., Leorat, J., & Pouquet, A. 1983, Astronomy and Astrophysics, 126, 51
- Gregory, S., Jardine, M., Gray, C., & Donati, J. 2010, Reports on Progress in Physics, 73, 126901
- Gregory, S. G., Donati, J.-F., & Hussain, G. A. 2016, arXiv preprint arXiv:1609.00273
- Hall, J. C., Lockwood, G., & Skiff, B. A. 2007, The Astronomical Journal, 133, 862
- Hébrard, É., Donati, J.-F., Delfosse, X., et al. 2016, Monthly Notices of the Royal Astronomical Society, 461, 1465
- Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
- Irwin, J., & Bouvier, J. 2009, in IAU Symp, Vol. 258
- Jardine, M., Barnes, J. R., Donati, J.-F., & Cameron, A. C. 1999, Monthly Notices of the Royal Astronomical Society, 305, L35
- Jardine, M., Collier Cameron, A., & Donati, J.-F. 2002, Monthly Notices of the Royal Astronomical Society, 333, 339
- Jeffers, S., Petit, P., Marsden, S., et al. 2014, Astronomy & Astrophysics, 569, A79
- Johnstone, C., Güdel, M., Lüftinger, T., Toth, G., & Brott, I. 2015, Astronomy & Astrophysics, 577, A27
- Johnstone, C., Jardine, M., & Mackay, D. 2010, Monthly Notices of the Royal Astronomical Society, 404, 101
- Kawaler, S. D. 1988, The Astrophysical Journal, 333, 236
- Kay, C., & Opher, M. 2014, in American Astronomical Society Meeting Abstracts# 224, Vol. 224
- Kay, C., Opher, M., & Evans, R. M. 2013, The Astrophysical Journal, 775, 5
- Keppens, R., & Goedbloed, J. 1999, Astron. Astrophys, 343, 251
- —. 2000, The Astrophysical Journal, 530, 1036
- Matt, S., & Pudritz, R. E. 2008, The Astrophysical Journal, 678, 1109
- Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, The Astrophysical Journal Letters, 799, L23
- Matt, S. P., MacGregor, K. B., Pinsonneault, M. H., & Greene, T. P. 2012, The Astrophysical Journal Letters, 754, L26
- McComas, D., Barraclough, B., Funsten, H., et al. 2000, Journal of Geophysical Research: Space Physics, 105, 10419
- McFadden, P., Merrill, R., McElhinny, M., & Lee, S. 1991, Journal of Geophysical Research: Solid Earth, 96, 3923
- McQuillan, A., Aigrain, S., & Mazeh, T. 2013, Monthly Notices of the Royal Astronomical Society, 432, 1203
- Meibom, S., Mathieu, R. D., Stassun, K. G., Liebesny, P., & Saar, S. H. 2011, The Astrophysical Journal, 733, 115
- Mestel, L. 1968, Monthly Notices of the Royal Astronomical Society, 138, 359
- —. 1984, in Cool Stars, Stellar Systems, and the Sun (Springer), 49
- Mignone, A. 2009, Memorie della Societa Astronomica Italiana Supplementi, 13, 67
- Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228
- Morgenthaler, A., Petit, P., Morin, J., et al. 2011, Astronomische Nachrichten, 332, 866
- Morgenthaler, A., Petit, P., Saar, S., et al. 2012, Astronomy & Astrophysics, 540, A138
- Morin, J., Donati, J.-F., Petit, P., et al. 2008, Monthly Notices of the Royal Astronomical Society, 390, 567
- Nicholson, B., Vidotto, A., Mengel, M., et al. 2016, Monthly Notices of the Royal Astronomical Society, 459, 1907
- Oran, R., Landi, E., van der Holst, B., et al. 2015, The Astrophysical Journal, 806, 55
- Parker, E. 1965, Space Science Reviews, 4, 666
- Parker, E. N. 1958, The Astrophysical Journal, 128, 664
- Petit, P., Dintrans, B., Morgenthaler, A., et al. 2009, Astronomy & Astrophysics, 508, L9
- Petit, P., Dintrans, B., Solanki, S., et al. 2008, Monthly Notices of the Royal Astronomical Society, 388, 80
- Pinto, R., Brun, A., & Rouillard, A. 2016, Astronomy & Astrophysics, 592, A65
- Pinto, R. F., Brun, A. S., Jouve, L., & Grappin, R. 2011, The Astrophysical Journal, 737, 72
- Reiners, A., & Mohanty, S. 2012, The Astrophysical Journal, 746, 43
- Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015a, The Astrophysical Journal, 798, 116
- Réville, V., Brun, A. S., Strugarek, A., et al. 2015b, The Astrophysical Journal, 814, 99
- Réville, V., Folsom, C. P., Strugarek, A., & Brun, A. S. 2016, The Astrophysical Journal, 832, 145
- Riley, P., Linker, J., Mikić, Z., et al. 2006, The Astrophysical Journal, 653, 1510
- Rosén, L., Kochukhov, O., & Wade, G. A. 2015, The Astrophysical Journal, 805, 169
- Rosner, R., Golub, L., & Vaiana, G. 1985, Annual review of astronomy and astrophysics, 23, 413
- Saikia, S. B., Jeffers, S., Morin, J., et al. 2016, Astronomy & Astrophysics, 594, A29
- Sakurai, T. 1990, Computer Physics Reports, 12, 247
- Schrijver, C. J., DeRosa, M. L., et al. 2003, The Astrophysical Journal, 590, 493
- See, V., Jardine, M., Vidotto, A., et al. 2015, Monthly Notices of the Royal Astronomical Society, 453, 4301
- —. 2016, Monthly Notices of the Royal Astronomical Society, 462, 4442
- —. 2017, Monthly Notices of the Royal Astronomical Society, stw3094
- Skumanich, A. 1972, The Astrophysical Journal, 171, 565
- Soderblom, D. 1983, The Astrophysical Journal Supplement Series, 53, 1
- Stauffer, J., Rebull, L., Bouvier, J., et al. 2016, The Astronomical Journal, 152, 115
- Strugarek, A., Brun, A., Matt, S., et al. 2014a, in SF2A-2014: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 279
- Strugarek, A., Brun, A. S., Matt, S. P., & Réville, V. 2014b, The Astrophysical Journal, 795, 86
- Tóth, G. 2000, Journal of Computational Physics, 161, 605
- Ud-Doula, A., Owocki, S. P., & Townsend, R. H. 2009, Monthly Notices of the Royal Astronomical Society, 392, 1022
- Usmanov, A. V., Goldstein, M. L., & Matthaeus, W. H. 2014, The Astrophysical Journal, 788, 43
- Ustyugova, G., Koldoba, A., Romanova, M., & Lovelace, R. 2006, The Astrophysical Journal, 646, 304
- Van der Holst, B., Manchester IV, W., Frazin, R., et al. 2010, The Astrophysical Journal, 725, 1373
- van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, The Astrophysical Journal, 782, 81
- Van Saders, J. L., & Pinsonneault, M. H. 2013, The Astrophysical Journal, 776, 67
- Vidotto, A., Jardine, M., Morin, J., et al. 2014a, Monthly Notices of the Royal Astronomical Society, 438, 1162
- Vidotto, A., Jardine, M., Opher, M., Donati, J., & Gombosi, T. 2011, in 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, Vol. 448, 1293
- Vidotto, A., Gregory, S., Jardine, M., et al. 2014b, Monthly Notices of the Royal Astronomical Society, 441, 2361
- Washimi, H., & Shibata, S. 1993, Monthly Notices of the Royal Astronomical Society, 262, 936
- Weber, E. J., & Davis, L. 1967, The Astrophysical Journal, 148, 217
- Wolk, S., Harnden Jr, F., Flaccomio, E., et al. 2005, The Astrophysical Journal Supplement Series, 160, 423
- Wood, B. E. 2004, Living Reviews in Solar Physics, 1, 1
- Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2004, The Astrophysical Journal Supplement Series, 152, 261
- Zanni, C., & Ferreira, J. 2009, Astronomy & Astrophysics, 508, 1117