The Effect of Combined Magnetic Geometries on Thermally Driven Winds. I. Interaction of Dipolar and Quadrupolar Fields

This work uses the MHD code PLUTO (Mignone et al. 2007; Mignone 2009), a finite-volume code that 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 (${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$) 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 the magnetic field:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot \rho {\boldsymbol{v}}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{m}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{mv}}-{\boldsymbol{BB}}+{\boldsymbol{I}}{p}_{T})=\rho {\boldsymbol{a}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot ((E+{p}_{T}){\boldsymbol{v}}-{\boldsymbol{B}}({\boldsymbol{v}}\cdot {\boldsymbol{B}}))={\boldsymbol{m}}\cdot {\boldsymbol{a}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{vB}}-{\boldsymbol{Bv}})=0.\end{eqnarray} \tag{ 4 }$$

Here ρ is the mass density, v is the velocity field, a is the gravitational acceleration, B is the magnetic field¹ , p_T = p + B²/2 is the combined thermal and magnetic pressure, and m = ρv is the momentum density. The total energy density is written as $E=\rho \epsilon$+ m²/(2ρ) + B²/2, with $\epsilon$ representing the internal energy per unit mass of the fluid. In addition, I is the identity matrix. A polytropic wind is used for this study, such that the closing equation of state takes the form $\rho \epsilon =p/(\gamma -1)$, where γ represents the polytropic index.

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 (r, θ) but allows for 3D axisymmetric solutions for the fluid flow and magnetic field using three vector components (r, θ, ϕ). The domain extends from one stellar radius (R_*) out to 60 R_* with a uniform grid spacing in θ and a geometrically stretched grid in r, which grows from an initial spacing of 0.01 to 1.08 R_* at the outer boundary. The computational mesh contains N_r × N_θ = 256 × 512 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 γ = 1.05 fills the domain, along with a superimposed 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 (B_r, B_θ) are held fixed at the stellar surface, while the poloidal components of the velocity (v_r, v_θ) are allowed to evolve in response to the magnetic field (the boundary is held with dv_r/dr = 0 and dv_θ/dr = 0). We then enforce the flow at the surface to be parallel to the magnetic field (${\boldsymbol{v}}| | {\boldsymbol{B}}$). The star rotates as a solid body, with B_ϕ linearly extrapolated into the boundary and v_ϕ set using the stellar rotation rate Ω_*,

$$\begin{eqnarray}&&{v}_{\phi }={{\rm{\Omega }}}_{* }r\,\sin \theta +\displaystyle \frac{{{\boldsymbol{v}}}_{{\rm{p}}}\cdot {{\boldsymbol{B}}}_{{\rm{p}}}}{| {{{\boldsymbol{B}}}_{{\rm{p}}}| }^{2}}{B}_{\phi },\end{eqnarray} \tag{ 5 }$$

where the subscript "p" denotes the poloidal components (r, θ) of a given vector. This condition enforces an effective rotation rate for the field lines that, 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 that 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 60 R_* allowing for the outward transfer of mass, momenta, and magnetic field, along with an axisymmetric condition along the rotation axis (θ = 0 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 artifacts 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 G. Pantolmos & S. Matt (2017, in preparation) and discussed in 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 torque and mass loss are used.

Inputs for the simulations are given as ratios of characteristic speeds that control key parameters such as the wind temperature (c_s/v_esc), field strength (v_A/v_esc), and rotation rate (v_rot/v_kep). Where ${c}_{s}=\sqrt{\gamma p/\rho }$ is the sound speed at the surface, ${v}_{A}={B}_{* }/\sqrt{4\pi \rho }$ is the Alfvén speed at the north pole, v_rot is the rotation speed at the equator, ${v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{* }/{R}_{* }}$ is the surface escape speed, and ${v}_{\mathrm{kep}}=\sqrt{{{GM}}_{* }/{R}_{* }}$ 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 of our simulations. Matt et al. (2012) showed that the nonlinear effects of rotation on their torque scaling can be neglected for slow rotators. They defined velocities as a fraction of the breakup speed,

$$\begin{eqnarray}&&{\left.f=\displaystyle \frac{{v}_{\mathrm{rot}}}{{v}_{\mathrm{kep}}}\right|}_{r={R}_{* },\theta =\pi /2}=\displaystyle \frac{{{\rm{\Omega }}}_{* }{R}_{* }^{3/2}}{{({{GM}}_{* })}^{1/2}}.\end{eqnarray} \tag{ 6 }$$

The Alfvén radius remains independent of the stellar spin rate until f ≈ 0.03, after which the effects of fast rotation start to be important. For this study, a solar rotation rate is chosen (f = 4.46 × 10⁻³) that is well within the slow-rotator regime. We set the temperature of the wind with c_s/v_esc = 0.25, higher than the c_s/v_esc = 0.222 used 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 v_A/v_esc, which is defined here with the Alfvén speed on the stellar north pole (see next section). Table 1 lists all of our variations of v_A/v_esc for each geometry.

Table 1.  Input Parameters and Results from the 70 Simulations

Case	${{ \mathcal R }}_{\mathrm{dip}}$	vA/vesc	$\langle {R}_{A}\rangle /{R}_{* }$	ϒ	Ro/R*	ϒopen	$\langle v({R}_{A})\rangle /{v}_{\mathrm{esc}}$	Case	${{ \mathcal R }}_{\mathrm{dip}}$	vA/vesc	$\langle {R}_{A}\rangle /{R}_{* }$	ϒ	Ro/R*	ϒopen	$\langle v({R}_{A})\rangle /{v}_{\mathrm{esc}}$
1	1	0.1	3.06	11.1	1.31	294	0.123	36	0.3	2	5.66	2930	2.61	2040	0.242
2	1	0.3	4.19	73.2	1.88	819	0.183	37	0.3	3	6.76	6850	3.01	3460	0.283
3	1	0.5	5.05	192	2.33	1450	0.221	38	0.3	6	9.41	31200	3.8	8840	0.360
4	1	1	6.88	773	2.95	3550	0.287	39	0.3	12	13	137000	5.05	21600	0.432
5	1	1.5	8.56	1880	3.41	6530	0.334	40	0.3	24	15.7	360000	6.18	37300	0.476
6	1	2	10	3660	3.8	9970	0.367	41	0.2	0.1	2.43	10.7	1.2	120	0.078
7	1	3	12.6	9280	4.54	18100	0.414	42	0.2	0.3	2.96	72.4	1.54	245	0.109
8	1	6	18.2	43900	6.07	47000	0.463	43	0.2	0.5	3.33	190	1.76	368	0.129
9	1	12	25.1	178000	8	109000	0.544	44	0.2	1	4.04	729	2.1	701	0.163
10	1	24	29.6	452000	9.75	180000	0.543	45	0.2	1.5	4.61	1630	2.39	1070	0.187
11	0.8	0.1	2.51	11.2	1.2	245	0.114	46	0.2	2	5.09	2930	2.56	1480	0.205
12	0.8	0.3	3.89	73.5	1.76	651	0.168	47	0.2	3	5.92	6840	2.9	2390	0.240
13	0.8	0.5	4.64	192	2.1	1120	0.203	48	0.2	6	7.93	31600	3.58	5890	0.301
14	0.8	1	6.19	751	2.73	2620	0.261	49	0.2	12	10.4	129000	4.54	13500	0.392
15	0.8	1.5	7.6	1780	3.12	4690	0.305	50	0.2	24	12.6	359000	5.56	24500	0.439
16	0.8	2	8.88	3390	3.46	7210	0.339	51	0.1	0.1	2.44	10.5	1.2	121	0.079
17	0.8	3	11.1	8660	4.14	13100	0.386	52	0.1	0.3	2.95	71.3	1.54	243	0.110
18	0.8	6	16.3	41700	5.67	35000	0.463	53	0.1	0.5	3.29	188	1.76	358	0.129
19	0.8	12	22.9	183000	7.77	84500	0.531	54	0.1	1	3.92	722	2.16	652	0.164
20	0.8	24	27.9	475000	9.07	147000	0.560	55	0.1	1.5	4.41	1620	2.44	964	0.190
21	0.5	0.1	2.63	11.4	1.14	168	0.095	56	0.1	2	4.81	2890	2.61	1290	0.208
22	0.5	0.3	3.38	74.1	1.54	407	0.140	57	0.1	3	5.53	6840	2.9	2050	0.244
23	0.5	0.5	3.94	191	1.82	674	0.169	58	0.1	6	7.13	33900	3.52	4850	0.317
24	0.5	1	5.11	736	2.33	1500	0.223	59	0.1	12	8.96	149000	4.31	10300	0.376
25	0.5	1.5	6.11	1660	2.67	2510	0.259	60	0.1	24	10.5	452000	5.16	17700	0.408
26	0.5	2	7.03	3050	2.95	3740	0.289	61	0	0.1	2.47	10.2	1.2	127	0.081
27	0.5	3	8.65	7500	3.46	6670	0.334	62	0	0.3	2.98	70.3	1.59	256	0.113
28	0.5	6	12.6	36600	4.6	17900	0.407	63	0	0.5	3.33	185	1.82	377	0.134
29	0.5	12	18.3	172000	6.3	46000	0.464	64	0	1	3.96	715	2.22	682	0.168
30	0.5	24	23	485000	7.49	83300	0.519	65	0	1.5	4.44	1600	2.5	1010	0.196
31	0.3	0.1	2.46	11	1.14	124	0.077	66	0	2	4.83	2890	2.67	1350	0.214
32	0.3	0.3	3.04	73.4	1.48	268	0.109	67	0	3	5.54	6950	2.95	2150	0.252
33	0.3	0.5	3.46	191	1.71	420	0.130	68	0	6	6.98	34900	3.63	4910	0.326
34	0.3	1	4.3	733	2.1	870	0.171	69	0	12	8.46	158000	4.43	9970	0.390
35	0.3	1.5	5.03	1630	2.39	1420	0.215	70	0	24	9.65	584000	5.16	16400	0.421

Download table as:  ASCIITypeset image

Due to the use of characteristic speeds as simulation inputs, our results can be scaled to any stellar parameter. For example, using solar parameters, the wind is driven by a coronal temperature of ≈1.4 MK, and our parameter space covers a range of stellar magnetic field strengths from 0.9 to 87 G over the pole. Changing these normalizations will modify this range.

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, ${{ \mathcal R }}_{\mathrm{dip}}$, of the 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:

$$\begin{eqnarray}&&{B}_{r,\mathrm{dip}}(r,\theta )={B}_{* }^{l=1}{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{3}\cos \theta ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{B}_{\theta ,\mathrm{dip}}(r,\theta )=\displaystyle \frac{1}{2}{B}_{* }^{l=1}{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{3}\sin \theta ,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{B}_{r,\mathrm{quad}}(r,\theta )=\displaystyle \frac{1}{2}{B}_{* }^{l=2}{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{4}(3{\cos }^{2}\theta -1),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{B}_{\theta ,\mathrm{quad}}(r,\theta )={B}_{* }^{l=2}{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{4}\cos \theta \sin \theta .\end{eqnarray} \tag{ 10 }$$

The total field, comprised of the sum of the two geometries,

$$\begin{eqnarray}&&{\boldsymbol{B}}(r,\theta )={{\boldsymbol{B}}}_{\mathrm{dip}}(r,\theta )+{{\boldsymbol{B}}}_{\mathrm{quad}}(r,\theta ),\end{eqnarray} \tag{ 11 }$$

where the total polar field, ${B}_{* }={B}_{* }^{l=1}+{B}_{* }^{l=2}$, is controlled by the ${{ \mathcal R }}_{\mathrm{dip}}$ parameter,

$$\begin{eqnarray}&&{\left.{{ \mathcal R }}_{\mathrm{dip}}=\displaystyle \frac{{B}_{r,\mathrm{dip}}}{{B}_{r,\mathrm{dip}}+{B}_{r,\mathrm{quad}}}\right|}_{r={R}_{* },\theta =0}=\displaystyle \frac{{B}_{* }^{l=1}}{{B}_{* }}.\end{eqnarray} \tag{ 12 }$$

This work considers aligned magnetic moments such that ${{ \mathcal R }}_{\mathrm{dip}}$ ranges from 1 to 0, corresponding to all the field strengths in the dipolar or quadrupolar mode, respectively. As with v_A/v_esc, ${{ \mathcal R }}_{\mathrm{dip}}$ is calculated at the north pole. This sets the relative strengths of the dipole and quadrupole fields,

$$\begin{eqnarray}&&{B}_{* }^{l=1}={{ \mathcal R }}_{\mathrm{dip}}{B}_{* },\qquad {B}_{* }^{l=2}=(1-{{ \mathcal R }}_{\mathrm{dip}}){B}_{* },\end{eqnarray} \tag{ 13 }$$

Alternative parameterizations are commonly used in the analysis of ZDI observations and dynamo modeling. These communities use the surface-averaged field strengths, $\langle | B| \rangle$, or the ratio of magnetic energy density (E_m ∝ B²) stored within each of the dipole and quadrupole field modes at the stellar surface. During the solar magnetic cycle, values of ${B}_{\mathrm{quad}}^{2}/{B}_{\mathrm{dip}}^{2}$ can range from ≈10–100 at solar maximum to ≈10⁻² at solar minimum (DeRosa et al. 2012). A transformation from our parameter to the ratio of energies is simply given by

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\mathrm{quad}}^{2}}{{B}_{\mathrm{dip}}^{2}}=\displaystyle \frac{2}{3}\displaystyle \frac{{(1-{{ \mathcal R }}_{\mathrm{dip}})}^{2}}{{{ \mathcal R }}_{\mathrm{dip}}^{2}},\end{eqnarray} \tag{ 14 }$$

where the numerical prefactor 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 (${{ \mathcal R }}_{\mathrm{dip}}=0.5,0.1$). 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 antialigned, 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 and octupole fields, where the aligned and antialigned cases cause subtraction at the equator or poles, respectively (Gregory et al. 2016; A. Finley & S. Matt 2017, in preparation).

Zoom In Zoom Out Reset image size

Figure 1 indicates that even with equal quadrupole and dipole polar field strengths, ${{ \mathcal R }}_{\mathrm{dip}}=0.5$, the overall dipole topology will remain. In this case, the magnetic energy density in the dipolar mode is 1.5 times greater than that in the quadrupolar mode coupled 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 ${{ \mathcal R }}_{\mathrm{dip}}=0.1$. 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 (${B}_{\mathrm{quad}}^{2}/{B}_{\mathrm{dip}}^{2}\gt 1.0$). For this study, the pure dipolar and quadrupolar fields are used as controls (both of which were studied in detail within Réville et al. 2015a), and five mixed cases are parameterized by ${{ \mathcal R }}_{\mathrm{dip}}$ values (${{ \mathcal R }}_{\mathrm{dip}}=0.8$, 0.5, 0.3, 0.2, 0.1). We include ${{ \mathcal R }}_{\mathrm{dip}}=0.8$ to demonstrate the dominance of the dipole at higher values. Each ${{ \mathcal R }}_{\mathrm{dip}}$ value is given a unique identifying color that is maintained in all figures throughout this paper. Table 1 contains a complete list of parameters for all cases, which are numbered by increasing v_A/v_esc and quadrupole fraction.

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 that 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 ${{ \mathcal R }}_{\mathrm{dip}}=0.5,0.1$. The figure displays the different sizes of Alfvén surface that are produced for a comparable value of polar magnetic field strength. 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 colocated with the sonic and Alfvén surfaces (the fast magnetosonic surface always being the larger of the two surfaces).

Zoom In Zoom Out Reset image size

The field geometry is found to imprint itself onto the stellar wind velocity with regions of closed magnetic field confining the flow and creating areas of corotating plasma, referred to as dead zones (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 dead zone(s), around which a slower wind is produced. Similar to the solar wind, slower wind can be found on the open field lines near the boundary of the closed field (Fisk et al. 1998; Feldman et al. 2005; Riley et al. 2006). 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 the coronal holes and closed 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 a 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 multispeed 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 ${{ \mathcal R }}_{\mathrm{dip}}=0.3,0.2,0.1$ values (${B}_{\mathrm{quad}}^{2}/{B}_{\mathrm{dip}}^{2}$ 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 dead zone is established on the equator by the dipole geometry, whereas the quadrupole creates two over the midlatitudes. Mixed cases have intermediate states between the pure regimes. Within our simulations, the dead zones are accompanied by streamers that form above the closed field regions and drive a slower-speed wind than that from the open field regions. The dynamics of these streamers and their location and size are an interesting result of the changing topology of the flow.

Zoom In Zoom Out Reset image size

The dashed colored lines in Figure 3 show where the field polarity reverses using B_r = 0, 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 toward the northern hemisphere until the quadrupolar streamers are recovered, both sitting at midlatitudes. This motion can also be seen for fixed ${{ \mathcal R }}_{\mathrm{dip}}$ cases as the magnetic field strength is decreased. For a given ${{ \mathcal R }}_{\mathrm{dip}}$ value, the current sheets sweep toward the southern hemisphere with increased polar field strength, in some cases (36 and 38) moving onto the axis of rotation. This is the opposite behavior to decreasing the ${{ \mathcal R }}_{\mathrm{dip}}$ value; i.e., the streamer configuration is seen to take a more dipolar morphology as the field strength is increased. Additionally in Figure 3, for low field strengths, each ${{ \mathcal R }}_{\mathrm{dip}}$ produces a comparable Alfvén surface with very similar morphology, all dominated by the quadrupolar mode.

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 ${{\boldsymbol{F}}}_{\mathrm{AM}}={\rm{\Lambda }}\rho {\boldsymbol{v}}$ (Keppens & Goedbloed 2000),

$$\begin{eqnarray}&&{\rm{\Lambda }}(r,\theta )=r\,\sin \theta \left({v}_{\phi }-\displaystyle \frac{{B}_{\phi }}{\rho }\displaystyle \frac{| {{\boldsymbol{B}}}_{{\rm{p}}}{| }^{2}}{{{\boldsymbol{v}}}_{{\rm{p}}}\cdot {{\boldsymbol{B}}}_{{\rm{p}}}}\right).\end{eqnarray} \tag{ 15 }$$

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,

$$\begin{eqnarray}&&\tau ={\int }_{A}{\rm{\Lambda }}\rho {\boldsymbol{v}}\cdot d{\boldsymbol{A}},\end{eqnarray} \tag{ 16 }$$

where A is the area of any surface enclosing the star. For illustrative purposes, Figure 3 shows the Alfvén surface colored by angular momentum flux (thick multicolored lines), which is seen to be strongly focused around the equatorial region. The angular momentum flux is calculated normal to the Alfvén surface,

$$\begin{eqnarray}&&\displaystyle \frac{d\tau }{{dA}}={\rm{\Lambda }}\rho {\boldsymbol{v}}\cdot \hat{{\boldsymbol{A}}}={{\boldsymbol{F}}}_{\mathrm{AM}}\cdot \hat{{\boldsymbol{A}}},\end{eqnarray} \tag{ 17 }$$

where $\hat{{\boldsymbol{A}}}$ is the normal unit vector to the Alfvén surface. The mass-loss rate from our wind solutions is calculated similarly to the torque,

$$\begin{eqnarray}&&\dot{M}={\int }_{A}\rho {\boldsymbol{v}}\cdot d{\boldsymbol{A}}.\end{eqnarray} \tag{ 18 }$$

Both expressions for the mass loss and torque are evaluated using spherical shells of area A that 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,

$$\begin{eqnarray}&&\langle {R}_{A}\rangle =\sqrt{\displaystyle \frac{\tau }{\dot{M}{{\rm{\Omega }}}_{* }}}.\end{eqnarray} \tag{ 19 }$$

Throughout this work, $\langle {R}_{A}\rangle$ is used as a normalized torque that accounts for the mass-loss rates that we do not control. Values of the average Alfvén radius are tabulated in Table 1, and $\langle {R}_{A}\rangle$ is shown in Figure 3 using dashed gray lines. For each case, the cylindrical Alfvén radius is offset inward of the maximum Alfvén radius from the simulation, a geometrical effect, as this corresponds to the average cylindrical R_A and includes variations in flow quantities as well. Exploring Figure 3, the motion of the dead zones/current sheets has little impact on the overall torque. For example, no abrupt increase in the Alfvén radius is seen from cases 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 that controls the location of the Alfvén surface.

We parameterize the magnetic and mass-loss properties using the "wind magnetization" defined by

$$\begin{eqnarray}&&{\rm{\Upsilon }}=\displaystyle \frac{{B}_{* }^{2}{R}_{* }^{2}}{\dot{M}{v}_{\mathrm{esc}}},\end{eqnarray} \tag{ 20 }$$

where B_* 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; G. Pantolmos & S. Matt 2017, in preparation). 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 that describes the field as a whole. The polar strength works for these aligned fields but will easily break down for unaligned fields and antialigned axisymmetric odd l fields; thus, it suits the present study, but a move away from this parameter in future is warranted.

During analysis, the wind magnetization, ϒ, is treated as an independent parameter that determines the Alfvén radius $\langle {R}_{A}\rangle$ and thus the torque τ. We increase ϒ by setting a larger v_A/v_esc, creating a stronger global magnetic field. Table 1 displays all the input values of ${{ \mathcal R }}_{\mathrm{dip}}$ and v_A/v_esc, 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 ${\rm{\Upsilon }}\mbox{--}{{ \mathcal R }}_{\mathrm{dip}}$ space. Cases are color-coded by their ${{ \mathcal R }}_{\mathrm{dip}}$ value, a convention that is continued throughout this work.

Zoom In Zoom Out Reset image size

The efficiency of the magnetic braking mechanism is known to be dependent on the magnetic field geometry. This has been previously shown for single-mode geometries (e.g., Réville et al. 2015a; Garraffo et al. 2016b). We first consider two pure geometries, dipole and quadrupole, using the formulation from Matt & Pudritz (2008),

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{s}{{\rm{\Upsilon }}}^{{m}_{s}},\end{eqnarray} \tag{ 21 }$$

where K_s and m_s are fitting parameters for the pure dipole and quadrupole cases using the surface field strength. Here we empirically fit m_s; the interpretation of m_s is discussed in Matt & Pudritz (2008), Réville et al. (2015a), and G. Pantolmos & S. Matt (2017, in preparation), where it is determined to be dependent 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 versus the wind magnetizations for all cases (color-coded with their ${{ \mathcal R }}_{\mathrm{dip}}$ value). The solid lines show the scaling relations for dipolar (red) and quadrupolar (blue) geometries, as first shown in Réville et al. (2015a). We calculate best-fit values for K_s and m_s for the dipole and quadrupole, tabulated in Table 2. Values here differ due to our hotter wind (c_s/v_esc = 0.25 versus their c_s/v_esc = 0.222) using polar B_*, 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 consider 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.

Zoom In Zoom Out Reset image size

G. Pantolmos & S. Matt (2017, in preparation) 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 ${v}_{\mathrm{esc}}/\langle v({R}_{A})\rangle$. In their work, the authors determined the theoretical power-law dependence, m_l,th = 1/(2l + 2), from one-dimensional analysis. In this formulation, the slope of the power law is controlled only by the order of the magnetic geometry, l, which is l = 1 and l = 2 for the dipole and quadrupole, respectively,

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{l}{\left[{\rm{\Upsilon }}\displaystyle \frac{{v}_{\mathrm{esc}}}{\langle v({R}_{A})\rangle }\right]}^{{m}_{l}},\end{eqnarray} \tag{ 22 }$$

where K_l and m_l are fit parameters to our wind solutions, tabulated in Table 2. The value of $\langle v({R}_{A})\rangle$ 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 accurately predict the power-law dependence for the two pure modes using the order of the spherical harmonic field, l. We show this in the right panel of Figure 5, where the Alfvén radii are plotted against the new parameter, ${\rm{\Upsilon }}{v}_{\mathrm{esc}}/\langle v({R}_{A})\rangle$. A similar qualitative behavior is shown in the scaling with ϒ in the left panel. Using the theoretical power-law dependencies, the dipolar (red) and quadrupolar (blue) slopes are plotted with m_l,th = 1/4 and 1/6, respectively. Using a single-fit constant K_l = 1 for both slopes within this figure shows good agreement with the simulation results.

More accurate values of K_l and m_l 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 behavior to the plot against ϒ.

Obvious trends are seen within the mixed-case scatter. A saturation to quadrupolar Alfvén radii values for lower ϒ and ${{ \mathcal R }}_{\mathrm{dip}}$ values is observed, along with a power-law trend with a dipolar gradient for higher ϒ and ${{ \mathcal R }}_{\mathrm{dip}}$ values. This indicates that both geometries play a role in governing the lever arm, with the dipole dominating the braking process at higher wind magnetizations.

Observationally, the field geometries of cool stars are, at large scales, dominated by the dipole mode, with higher-order l modes playing smaller roles in shaping the global field. It is the global field that controls the spin-down torque in the magnetic braking process. Higher-order modes (such as the quadrupole) radially decay 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 ϒ_dip, which only takes into account the dipole's field strength,

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}_{\mathrm{dip}}={\left(\displaystyle \frac{{B}_{* }^{l=1}}{{B}_{* }}\right)}^{2}\displaystyle \frac{{B}_{* }^{2}{R}_{* }^{2}}{\dot{M}{v}_{\mathrm{esc}}}={{ \mathcal R }}_{\mathrm{dip}}^{2}{\rm{\Upsilon }}.\end{eqnarray} \tag{ 23 }$$

Taking as a hypothesis that the field controlling the location of the Alfvén radius is the dipole component, a power-law scaling using ϒ_dip can be constructed in the same form as that of Matt & Pudritz (2008),

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{s,\mathrm{dip}}{[{{\rm{\Upsilon }}}_{\mathrm{dip}}]}^{{m}_{s,\mathrm{dip}}}={K}_{s,\mathrm{dip}}{[{{ \mathcal R }}_{\mathrm{dip}}^{2}{\rm{\Upsilon }}]}^{{m}_{s,\mathrm{dip}}}.\end{eqnarray} \tag{ 24 }$$

Substitution of the dipole component into Equation (22) similarly gives

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{l,\mathrm{dip}}{\left[{{ \mathcal R }}_{\mathrm{dip}}^{2}{{\rm{\Upsilon }}}_{}\displaystyle \frac{{v}_{\mathrm{esc}}}{\langle v({R}_{A})\rangle }\right]}^{{m}_{l,\mathrm{dip}}},\end{eqnarray} \tag{ 25 }$$

where K_s,dip, m_s,dip, K_l,dip, and m_l,dip 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 ${{ \mathcal R }}_{\mathrm{dip}}$ values used in our simulations. Mixed cases that 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 ${{ \mathcal R }}_{\mathrm{dip}}$, with deviation from this caused by a regime change onto the quadrupolar slope (shown by colored dashed lines).

Zoom In Zoom Out Reset image size

The behavior of our simulated winds, despite using a combination of field geometries, simply follows existing scaling relations with this modification. In general, the dipole (ϒ_dip) prediction shows good agreement with the simulated wind models, except in cases where the Alfvén surface is close to the star. In these cases, the quadrupole mode still has a 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 in Figure 5 for low values of ϒ and ${{ \mathcal R }}_{\mathrm{dip}}$.

The mixed field $\langle {R}_{A}\rangle$ 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 ϒ_crit,

$$\begin{eqnarray}\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}=\left\{\begin{array}{ll}{K}_{s,\mathrm{dip}}{[{{ \mathcal R }}_{\mathrm{dip}}^{2}{\rm{\Upsilon }}]}^{{m}_{s,\mathrm{dip}}}, & \mathrm{if}\ {\rm{\Upsilon }}\gt {{\rm{\Upsilon }}}_{\mathrm{crit}}({{ \mathcal R }}_{\mathrm{dip}}),\\ {K}_{s,\mathrm{quad}}{[{\rm{\Upsilon }}]}^{{m}_{s,\mathrm{quad}}}, & \mathrm{if}\ {\rm{\Upsilon }}\leqslant {{\rm{\Upsilon }}}_{\mathrm{crit}}({{ \mathcal R }}_{\mathrm{dip}})\end{array}\right.,\end{eqnarray} \tag{ 26 }$$

where ϒ_crit is the location of the intercept for the dipole component and pure quadrupole scalings,

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}_{\mathrm{crit}}({{ \mathcal R }}_{\mathrm{dip}})={\left[\displaystyle \frac{{K}_{s,\mathrm{dip}}}{{K}_{s,\mathrm{quad}}}{{ \mathcal R }}_{\mathrm{dip}}^{2{m}_{s,\mathrm{dip}}}\right]}^{\tfrac{1}{{m}_{s,\mathrm{quad}}-{m}_{s,\mathrm{dip}}}}.\end{eqnarray} \tag{ 27 }$$

The solid lines in Figure 4 show the value of ϒ_crit (Equation (27)), dividing the two regimes. Specifically, the solutions above the black line behave as if only the dipole component (ϒ_dip) is governing the Alfvén radius.

Transitioning from regimes is not perfectly abrupt. Therefore, producing an analytical solution for the mixed cases that includes this behavior would increase the accuracy for stars near the regime change. For example, 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 in the error on this simpler form and is not easily generalized to higher topologies. For practical purposes, the scaling of Equations (26) and (27) accurately predicts the simulation torque with increasing magnetic field strength for a variety of dipole fractions. We therefore present the simplest available solution, leaving the generalized form to be developed in future work.

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,

$$\begin{eqnarray}&&{\rm{\Phi }}(r)={\oint }_{r}| {\boldsymbol{B}}\cdot d{\boldsymbol{A}}| .\end{eqnarray} \tag{ 28 }$$

Considering the initial potential fields of the two pure modes, this is simply a power law in field order l,

$$\begin{eqnarray}&&{\rm{\Phi }}{(r)}_{P}={{\rm{\Phi }}}_{* }{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{l},\end{eqnarray} \tag{ 29 }$$

where l = 1 dipole and l = 2 quadrupole, and we denote the flux with P for the potential field. Figure 7 displays the flux decay of all values of v_A/v_esc for each ${{ \mathcal R }}_{\mathrm{dip}}$ value (gray lines). The behavior is qualitatively identical to that observed in previous works (e.g., Schrijver et al. 2003; Johnstone et al. 2010; Vidotto et al. 2014b; 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 that, for higher l 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.

Zoom In Zoom Out Reset image size

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 B is the initial potential field for each mixed geometry. This allows us to calculate the radial evolution of the flux for a given ${{ \mathcal R }}_{\mathrm{dip}}$, which we compare to the simulated cases. Figure 7 shows the flux normalized by the surface flux versus radial distance from the star. For each ${{ \mathcal R }}_{\mathrm{dip}}$ value, the magnetic flux decay of the potential field (black line) is shown with the different-strength v_A/v_esc simulations (gray lines). A comparison of the flux decay for all potential magnetic geometries is given in the bottom right panel, showing, as expected, the increasingly quadrupolar fields decaying faster.

In this study, we control v_A/v_esc, which, for a given surface density, sets the polar magnetic field strength for our simulations. The stellar flux for different topologies and the same B_* will differ and must be taken into account in order to describe the dipole and quadrupolar components (dashed red and blue lines) in Figure 7. We plot the magnetic flux of the potential field quadrupole component alone with a dashed blue line for each ${{ \mathcal R }}_{\mathrm{dip}}$ value,

$$\begin{eqnarray}&&{\rm{\Phi }}{(r)}_{P,\mathrm{quad}}=(1-{{ \mathcal R }}_{\mathrm{dip}}){{\rm{\Phi }}}_{* ,\mathrm{quad}}{\left(\displaystyle \frac{{R}_{* }}{r}\right)}^{2},\end{eqnarray} \tag{ 30 }$$

and, similarly, the potential field dipole component of the magnetic flux,

$$\begin{eqnarray}&&{\rm{\Phi }}{(r)}_{P,\mathrm{dip}}={{ \mathcal R }}_{\mathrm{dip}}{{\rm{\Phi }}}_{* ,\mathrm{dip}}\left(\displaystyle \frac{{R}_{* }}{r}\right),\end{eqnarray} \tag{ 31 }$$

where in both equations the surface flux of a pure dipole/quadrupole (Φ_*,dip, Φ_*,quad) field is required to match our normalized 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 ${{ \mathcal R }}_{\mathrm{dip}}$ 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 an l = 2 power-law decay, which can be seen by shifting the blue dashed line upward to intercept Φ/Φ_* = 1 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 it is shown by the potential solution (black lines).

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 parameterized in terms of the open flux, Φ_open,

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}_{\mathrm{open}}=\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{open}}^{2}/{R}_{* }^{2}}{\dot{M}{v}_{\mathrm{esc}}},\end{eqnarray} \tag{ 32 }$$

which, ignoring the effects of rapid rotation, can be fit with

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{o}{[{{\rm{\Upsilon }}}_{\mathrm{open}}]}^{{m}_{o}},\end{eqnarray} \tag{ 33 }$$

where m_o and K_o are fitting parameters for the open flux formulation.

Using the open flux parameter, Figure 8 shows a collapse toward 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.³ We calculate the best fits for each pure mode separately, i.e., the dipole and quadrupole, tabulated in Table 3.

Zoom In Zoom Out Reset image size

G. Pantolmos & S. Matt (2017, in preparation) find solutions for thermally driven winds with different coronal temperatures. From these, they find that the wind acceleration profiles of a given wind very significantly alter the slope in R_A–ϒ_open 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 differences in the superradial 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 in our study are available in the Appendix. G. Pantolmos (2017, in preparation) are able to include the effects of acceleration in their scaling through multiplication of ϒ_open with ${v}_{\mathrm{esc}}/\langle v({R}_{A})\rangle$. The expected semianalytical solution from G. Pantolmos & S. Matt (2017, in preparation) is given as

$$\begin{eqnarray}&&\displaystyle \frac{\langle {R}_{A}\rangle }{{R}_{* }}={K}_{c}{\left[{{\rm{\Upsilon }}}_{\mathrm{open}}\displaystyle \frac{{v}_{\mathrm{esc}}}{\langle v({R}_{A})\rangle }\right]}^{{m}_{c}},\end{eqnarray} \tag{ 34 }$$

where the fit parameters are derived from 1D theory as constants, K_c,th = 1/4π and m_c,th = 1/2.

We are able to reproduce this power-law fit of ϒ_open, with the wind acceleration effects removed, in the right panel of Figure 8. Including all simulations in the fit, we arrive at values of ${K}_{c}=1.01{K}_{c,\mathrm{th}}\pm 0.07$ and ${m}_{c}=0.942{m}_{c,\mathrm{th}}\pm 0.009$ for the constants of proportionality and power-law dependence. However, a systematic difference is still seem from one ${{ \mathcal R }}_{\mathrm{dip}}$ value to another. More precise fits can be found for each geometry independently, but the systematic difference appearing in the right panel implies that a modification to our semianalytical formulations is required to describe the torque fully in terms of the open flux.

Here we show that the scaling law from Réville et al. (2015a) is improved with the modification from G. Pantolmos (2017, in preparation). 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.

The location of the field opening is an important distance. It is critical both 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 R_ss. The opening radius, R_o, we define as the radial distance at which the potential flux reaches the value of the open flux (Φ_P(R_o) = Φ_open). 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 R_o 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 nonpotential due to its interaction with the wind (which explains why the closed field regions seen in Figure 3 typically extend slightly beyond R_o). A similar smooth transition is produced with PFSS modeling.

In Figure 7, R_o is marked for each simulation and for comparative purposes in the bottom right panel. It is clear that smaller opening radii are found for lower ${{ \mathcal R }}_{\mathrm{dip}}$ 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. G. Pantolmos & S. Matt (2017, in preparation) 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, $\langle {R}_{A}\rangle$, and the opening radius, R_o, for the seven different geometries studied. It is interesting to consider the relative size of these radii, as they both characterize key dynamic properties for each stellar wind solution. For all cases shown in Figure 3, the opening radii are plotted with dashed red lines, allowing for the relative size to be compared with the cylindrical Alfvén radius, shown with dashed gray lines. With increasing magnetic field strength (ϒ), both radii are seen to grow from case to case; however, with increasing ${{ \mathcal R }}_{\mathrm{dip}}$, 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 versus the opening radii for all cases. Linear trends of R_A/R_o = 3.2 and 1.7 are indicated with dashed lines. For each ${{ \mathcal R }}_{\mathrm{dip}}$ value, the relationship between the Alfvén and opening radius ($\langle {R}_{A}\rangle /{R}_{o}$) 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.

Zoom In Zoom Out Reset image size

The dependence of the Alfvén radius and opening radius on field geometry and magnetization is a constraint on PFSS models, which are readily used with ZDI observations as a less computationally expensive alternative to MHD modeling (Jardine et al. 1999, 2002; Dunstone et al. 2008; Cohen et al. 2010; Johnstone et al. 2010; Réville et al. 2015b; Rosén et al. 2015). PFSS models are a useful tool; however, they require the source surface radius, R_ss, 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, that for a given ${{ \mathcal R }}_{\mathrm{dip}}$ value there exists a differing relation for the opening radius, as we define it here, to the Alfvén radius and magnetization. These trends are observed to continue for higher l mode fields (A. Finley & S. Matt 2017, in preparation), with $\langle {R}_{A}\rangle /{R}_{o}$ 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 magnetization. We find that the relationship of $\langle {R}_{A}\rangle /{R}_{o}$ changes in both cases. With fixed magnetization, the opening radius should move toward the star for higher-order fields to maintain a constant thermal driving. Maintaining the opening radius while increasing the field complexity infers that the wind has a reduced acceleration. Similarly, with increased wind magnetization, the opening radius should move further from the star. The value of R_o 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 R_o to the PFSS model shows that R_ss equals an approximately constant value of 3/2 R_o. Thus, conclusions made about the opening radii are constraints on future PFSS modeling.

A method for predicting R_o within our simulations remains unknown; however, it is understood that R_o is key to predicting the torque from our simulated winds. We do, however, find the ratio of $\langle {R}_{A}\rangle /{R}_{o}$ to be roughly constant for a given geometry, deviations from which may be numerical or suggest additional physics that we do not explore here.