{'fraction_non_alphanumeric': 0.030527511559815042, 'fraction_numerical': 0.019855682309083054, 'mean_word_length': 4.500968139412075, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


INTRODUCTION
The spin down of cool stars (M * 1.3M ) 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; Ω * ∝ τ −0.5 (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; *af472@exeter.ac.uk 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. 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;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. 2015See et al. , 2016Folsom 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.

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 (∇ · 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 magnetic field, Here ρ is the mass density, v is the velocity field, a is the gravitational acceleration, B is the magnetic field 1 , p T = p + B 2 /8π is the combined thermal and magnetic pressure, m is the momentum density given by ρv and E is the total energy density. The energy of the system is written as E = ρ + m 2 /(2ρ) + B 2 /2, with representing the internal energy per unit mass of the fluid. I is the identity matrix. A polytropic wind is used for this study, such that the closing equation of state takes the form ρ = p/(γ − 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 axisymetric solutions for the fluid flow and magnetic field using three vector components (r, θ, φ). The domain extends from one stellar radius (R * ) out to 60R * with a uniform grid spacing in θ and a geometrically stretched grid in r, which grows from an initial spacing of 0.01R * to 1.08R * 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 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 (B r , B θ ) are held fixed at the stellar surface, whilst 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 (v||B). The star rotates as a solid body, with B φ linearly extrapolated into the boundary and v φ set using the stellar rotation rate Ω * , where the subscript "p" denotes the poloidal components (r, θ) 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 60R * 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, 1 The PLUTO code operates with a factor of 1/ √ 4π absorbed into the normalisation of B. Tabulated parameters are given in cgs units with this factor incorporated. 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 (c s /v esc ), field strength (v A /v esc ) and rotation rate (v rot /v kep ). Where c s = γp/ρ is the sound speed at the surface, v A = B * / √ 4πρ is the Alfvén speed at the north pole, v rot is the rotation speed at the equator, v esc = 2GM * /R * is the surface escape speed and v kep = 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 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 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 −3 ), which is well within the slow rotator regime. We set the temperature of the wind with c s /v esc = 0.25, higher than c s /v esc = 0.222 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 v A /v esc , which is defined here with the Alfvén speed on the stellar north pole (see following Section). Table 1 lists all our variations of v A /v esc 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. . Initial magnetic configurations for a dipolar field, quadrupolar field and two mixed cases (red, green, magenta and blue for the dipole fractions of 100%, 50%, 10% and purely quadrupole respectively). Mixed cases have the dominant pure field geometry over-plotted in dashed colour. The combined fields add in the northern hemisphere and subtract in southern hemisphere because they belong to opposite field symmetry families. With as much as half the field strength in the quadrupole, shown in green, the topology of the field is still dominated by the dipole field.

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, R dip , 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 B * = B l=1 * +B l=2 * , is controlled by the R dip parameter, This work considers aligned magnetic moments such that R dip ranges from 1 to 0, corresponding to all the field strength in the dipolar or quadrupolar mode respectively. As with v A /v esc , R dip 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, |B| , or the ratio of magnetic energy density (E m ∝ B 2 ) stored within each of the dipole and quadrupole field modes at the stellar surface. During the solar magnetic cycle, values of B 2 quad /B 2 dip can range from ≈ 10 − 100 at solar maximum to ≈ 10 −2 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 (R 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 Figure 2. Logarithm of density normalised by the surface value for dipolar, quadrupolar and mixed magnetic fields for cases 7, 27, 57, 67 (see Table 1). The winds are initialised using the same initial polytropic parker wind solution with γ = 1.05 and cs/vesc = 0.25. Stellar rotation rate and magnetic field strength are set with f = 4.46 × 10 −3 and v A /vesc = 3.0. The Alfvén and sonic Mach surfaces are shown in blue and black respectively, in addition the fast and slow magnetosonic surfaces are indicated with dot-dash and dashed white lines. A transition from one to two streamers is seen with increasing quadrupolar field (decreasing R dip ), and the two combined field cases exhibit the top bottom asymmetry from the field addition and subtraction. 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, R 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 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 R 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 2 quad /B 2 dip > 1.0). 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 R dip values (R dip = 0.8, 0.5, 0.3, 0.2, 0.1). We include R dip = 0.8 to demonstrate the dominance of the dipole at higher values. Each R dip 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 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 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.

Morphology of the Field and Wind Outflow
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, R dip = 0.5, 0.1. 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 respectively), coloured by poloidal wind speed, with field lines in white. The current sheets are indicated by dashed lines, whose colour corresponds to their R dip value in future figures. The streamer configuration is modified by both changes to the field strength and mixing ratio. Increased field strength or R dip value tends to revolve the southern hemisphere streamer towards the south pole. The Alfvén surfaces have been coloured to show the flux of angular momentum normal to the surface [Units normalised by: 8 × 10 −6 ρ * v kep R * ]. The average Alfvén radius, R A , from equation (19) is shown in dashed grey. The sonic surface and opening radius are shown in solid black and dashed red respectively. The morphology and properties of the lower field cases are nearly indistinguishable, with only slight differences to the streamer locations. The reduction in torque with increasing quadrupolar fraction can be visually seen by moving down the grid. The most dipolar field sits in the top right panel and the most quadrupolar in the bottom left, these models are chosen to emphasis the transition in field dominance.
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). The dashed coloured lines within 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 towards the northern hemisphere until the quadrupolar streamers are recovered both sitting at mid latitudes. This motion can also be seen for fixed R dip cases as the magnetic field strength is decreased. For a given R dip 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 R dip 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 R dip produces a comparable Alfvén surface with very similar morphology, all dominated by the quadrupolar mode.

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 F AM = Λρv (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 A 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 multicoloured 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 A 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, R A 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. R A 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 R A 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 B * is the combined field strength at the pole. Previous studies that used this parameter defined it with the   . 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 l 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 R A 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 R 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 Υ − R dip space. Cases are colour-coded here by their R dip value, a convention which is continued throughout this work.

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 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 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 R dip 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 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 than 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 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 veloc-ity 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 esc / v(R A ) . In their work, the authors determine 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, where K l and m l are fit parameters to our wind solutions, tabulated in Table 2. The value of v(R A ) is calculated as an average of the velocity at all points on the Alfvén surface in the meridional plane. 2 Equation (22) is able to predict accurately 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, Υv esc / v(R A ) . 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 m l,th = 1/4 and m l,th = 1/6 respectively. Using a single fit constant K l = 1 for both sloes 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 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 R dip values is observed, along with a power law trend with a dipolar gradient for higher Υ and R 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 magnetisations.

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 l 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 Υ dip , 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  (21) and (22) Topology(l) Ks ms K l m l m l,th (l) Dipole (1) 1.49 ± 0.03 0.231 ± 0.003 0.92 ± 0.04 0.258 ± 0.005 0.250 Quadrupole (2) 1.72 ± 0.03 0.132 ± 0.003 1.11 ± 0.04 0.156 ± 0.004 0.167 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Υ 10 0 10 1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Υv esc / v(R A ) 10 0 10 1   power law scaling using Υ dip can be constructed in the same form as Matt & Pudritz (2008), Substitution of the dipole component into equation (22) similarly gives, 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  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 R dip 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 (Υ dip ) 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 R dip .
The mixed field R A 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 , (26) where Υ crit is the location of the intercept for the dipole component and pure quadrupole scalings, . (27) The solid lines in Figure 4 show the value of Υ crit , equation (27), diving the two regimes. Specifically, the solutions above the solid 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 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.

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 l, where l = 1 dipole and l = 2 quadrupole, 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 R dip 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 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.
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 R dip 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 R dip value, the magnetic flux decay of the potential field (black solid line) is shown with the different strength v A /v esc 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 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) in Figure 7. We plot the magnetic flux of the potential field quadrupole component alone in dotted blue for each R dip value, and similarly the potential field dipole component of the magnetic flux, where in both equations the surface flux of a pure dipole/quadrupole (Φ * ,dip , Φ * ,quad ) 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 R 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 a l = 2 power law decay, which can be seen by shifting the blue dashed line upwards 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 is shown by the potential solution in solid black.

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 Figure 7. Magnetic flux vs Radial distance for all cases studied within this work, compared with analytical predictions. Solid grey lines show the 10 simulation fields for each field geometry. Solutions of equation (28) for the potential field magnetic flux are shown in black solid lines for each R dip value. In each case, the flux of dipole and quadrupole components using a potential field are plotted with dashed red and blue respectively, equations (30) which, ignoring the effects of rapid rotation, can be fit with, 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. 3 We calculate best fits for each pure mode separately i.e. the dipole and quadrupole, tabulated in Table 3.
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 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 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 Υ open with v esc / v(R A ) . The expected semi-analytic solution from Pantolmos & Matt (in prep) is given, where the fit parameters are derived from onedimensional 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, on the right panel of Figure 8. Including all simulations in the fit, we arrive at values of K c = 1.01K c,th ± 0.07 and m c = 0.942m c,th ± 0.009 for the constants of proportionality and power law dependence. However a systematic difference is still seem from one R dip 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 Table 3 Open Flux Best Fit Parameters to equations (33) and (34) Topology(l) Ko mo Dipole (1) 0.37 ± 0.05 0.360 ± 0.006 Quadrupole (2) 0.62 ± 0.01 0.283 ± 0.002 Kc K c,th mc m c,th Topology Independent 0.08 ± 0.03 0.0796 0.471 ± 0.003 0.500 10 1 10 2 10 3 10 4 10 5 10 6 Υ open v esc / v(R A ) 10 0 10 1 10 1 10 2 10 3 10 4 10 5 10 6 Υ open 10 0 10 1  (34) is shown to reduce the scatter for all simulations. A systematic discrepancy is still seen from the single power law with changing geometry.
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.

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 R ss . The opening radius, R o , we define is 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 non-potential due to its interaction with the wind (which explains why the closed field regions seen in Figure 3 typically extend slightly beyond   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, R A , and the opening radius, R o , 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 R 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 vs 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 R dip value, the relationship between the Alfvén and opening radius ( R A /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.
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(Jardine et al. , 2002Dunstone 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, 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. 2015See et al. , 2017. We find however for a given R dip 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 l mode fields (Finley & Matt. in prep), with R A /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 magnetisation. We find the relationship of R A /R o 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 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 modelling.
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 R A /R o to be roughly constant for a given geometry, deviations from which may be numerical or suggest additional physics which we do not explore here.

CONCLUSION
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, R dip , 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 R dip 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.  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 . 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, K s,dip , m s,dip , K s,quad and m s,quad (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, (36) At higher Υ and dipole fractions, the torque is only dependent on the dipolar component of the field, (38) 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 Υ crit 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.  (Hunter 2007). Table 4 Predicting ms and mo using q = 0.8 ± 0.1 Topology(l) ms m s,th (l, q) mo m o,th (q) Dipole (1) 0.231 ± 0.003 0.21 ± 0.01 0.360 ± 0.006 0.36 ± 0.02 Quadrupole (2) 0.132 ± 0.003 0.15 ± 0.01 0.283 ± 0.002 0.36 ± 0.02

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, m, of the torque formulations used within this work.
Using the one-dimensional framework, the field strength is assumed to decay as a power law B(r) = B * (R * /r) l+2 , 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 Υ open with v esc / v(R A ) . The power law dependences then becomes, m l,th = 1/(2l + 2), and similarly, m c,th = 1/2.
The modified dependent parameter, Υv esc / v(R A ) , 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, v(R A ), is required for the semi-analytic formulations. v(R A ) is shown by Pantolmos   Figure 10. Scatter of the average Alfvén speed at the Alfvén surface as a function of the average Alfvén radius. The dashed line shows a hydrodynamic parker wind with cs/vesc = 0.25, the solid line shows a fit to all our simulation data. Variation is seen between the dipolar and quadrupolar data towards the extreme values of the Alfvén radius. The combined average wind acceleration profile in black gives a q = 0.84. The winds in our simulations are set with a higher coronal temperature than Réville et al. (2015a), thus show a larger acceleration (they produce q ≈ 0.7). & Matt. (in prep) and Réville et al. (2015a) to share a similar profile to a one-dimensional thermal wind, v(r). 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 R dip . Towards higher values of the Alfvén radius, a noticeable separation starts to develop between geometries. This range is accessed less by the higher l 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 R A 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 q, v(R A )/v esc = (R A /R * ) q . (A3) A single power law with q = 0.84 is fit to the simulation data, which is chosen for simplicity within the 1D formalism. The use of this q parameter is approximate if v(R A ) is a power law in R A , 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)), m s,th = 1/(2l + 2 + q), which includes geometric and wind acceleration parameters in the form of l and q respectively. Using this result, m s,th is computed for both the dipole (l = 1) and quadrupole (l = 2) geometries in Table 4, and compared to the simulation results with good agreement.
The result is independent of geometry, l. As before the q parameter approximates the wind driving as a power law in radius, which is fit with a single power law for both geometries such that m o,th 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 m o 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 q does provide a way to approximately include the contribution of the wind acceleration to the fit power-law exponents m o and m s . 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.