A systematic description of wind-driven protoplanetary discs
AAstronomy & Astrophysics manuscript no. selfsimilar © ESO 2021January 27, 2021
A systematic description of wind-driven protoplanetary discs (cid:63)
Geo ff roy R. J. Lesur Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, Francee-mail: [email protected]
Received ——–; accepted January 24, 2021
ABSTRACT
Aims.
Planet forming discs are believed to be very weakly turbulent in the regions outside of 1 AU. For this reason, it is now believedthat magnetised winds could be the dominant mechanism driving accretion in these systems. However, there is today no self-consistentway to describe discs subject to a magnetised wind, in a way similar to the α disc model. In this article, I explore in a systematic waythe parameter space of wind-driven protoplanetary discs and present scaling laws which can be used in reduced models à la α disc. Methods.
I compute a series of self-similar wind solutions, assuming the disc is dominated by ambipolar and Ohmic di ff usions. Thesesolution are obtained by looking for stationary solutions in the finite-volume code PLUTO using a relaxation method and continuation. Results.
Self similar solutions are obtained for values of plasma β ranging from 10 to 10 , for several Ohmic and ambipolar di ff usionstrengths. Mass accretions rates of the order of 10 − M (cid:12) / yr are obtained for poloidal field strength β = O (10 ) or equivalently 1 mGat 10 AU. In addition, the ejection e ffi ciency is always close to 1, implying that wind mass loss rate can be larger than the inner massaccretion rate if the wind-emitting region is large. The resulting magnetic lever arms are typically lower than 2, possibly reaching1.5 in weakest field cases. Remarkably, the mean transport properties (accretion rate, mass loss rate) depend mostly on the fieldstrength and much less on the disc di ff usivities or surface density. The disc internal structure is nevertheless strongly a ff ected byOhmic resistivity, strongly resistive discs being subject to accretion at the surface while ambipolar only models lead to mid-planeaccretion. Finally, I provide a complete set of scaling laws and semi-analytical wind solutions, which can be used to fit and interpretobservations. Conclusions.
Magnetised winds are unavoidable in protoplanetary discs as soon as they are embedded in an ambient poloidal magneticfield. Very detailed disc microphysics are not always needed to describe them, and simplified models such as self-similar solutionsmanage to capture most of the physics seen in full 3D simulations. The remaining di ffi culty to get a complete theory of wind-drivenaccretion lies in the transport of the large scale field, which remains poorly constrained and not well understood. Key words. accretion discs, outflows
1. Introduction
Protoplanetary discs are relatively cold and dense objects, whichtypically last for a few million years around young stellar ob-jects (YSOs). However, it is well known that these discs are pri-marily accretion discs, in which matter is slowly falling ontothe central star. This accretion rate has now been measured indozens of objects, and typically lies in the range of 10 − to10 − M (cid:12) / yr. Accretion in astrophysical discs is usually thoughtto be due to a magnetohydrodynamic (MHD) instability, themagneto-rotational instability (MRI, Balbus & Hawley 1991),which transports angular momentum outwards and mass in-wards. It was quickly realised that applying the MRI in the con-text of protoplanetary discs is a di ffi cult task because of thedramatically low ionisation fraction expected in these objects.Indeed, the inclusion of Ohmic (Gammie 1996) and ambipolardi ff usions (Perez-Becker & Chiang 2011) led to the conclusionthat MRI cannot account for the observed accretion rates in thesediscs by two to three orders of magnitude.The lack of a proper mechanism to trigger turbulence andangular momentum transport in protoplanetary discs revived anold idea mostly abandoned since the finding of the MRI: mag-netised winds. Indeed, while the disc bulk might be too di ff usiveto sustain MHD turbulence, it is still weakly coupled to the am- (cid:63) The self-similar solutions presented in this paper are available fordownload on github: https: // github.com / glesur / PPDwind bient large scale field originating from the parent cloud of theYSO. This ambient field can in principle be enough to triggera magnetic breaking of the disc, leading to mass accretion andthe formation of a wind. This idea was originally proposed inthe context of YSOs by Wardle & Koenigl (1993). It was laterrevived via numerical simulations, first in local shearing box se-tups (Bai & Stone 2013; Simon et al. 2013; Lesur et al. 2014) andthen in global simulations (Gressel et al. 2015; Béthune et al.2017; Bai 2017; Wang et al. 2019; Gressel et al. 2020). All ofthese simulations include Ohmic and ambipolar di ff usion, whilesome of them also include the Hall e ff ect which is believed tobe important in the densest regions of the disc (around 1 AU).The strength of these di ff usions is in turn computed from ioni-sation and chemical models of various complexity. In addition,almost all of these models include some sort of heating of thedisc atmosphere, designed to mimic non-thermal radiation heat-ing expected in these systems (e.g. Thi et al. 2019). In the endof the day, each published model is of very high complexity,and only a handful of simulations can be performed exploring a(very) limited subspace of parameter space.In this context, the purpose of this work is to take a step back,simplifying the physics in order to explore in a more or less sys-tematic manner the coupling between protoplanetary discs andmagnetised winds. This work follows in essence the approachinitiated by Ferreira (1997) on globally self-similar models andlater extended to viscous (Casse & Ferreira 2000a) and weakly Article number, page 1 of 13 a r X i v : . [ a s t r o - ph . S R ] J a n & A proofs: manuscript no. selfsimilar magnetised (Jacquemin-Ide et al. 2019) discs. I however use arelatively di ff erent technic to get the wind solutions, which Ipresent in the next section, along with the physical model. Ithen focus on a fiducial series of wind solution and explore theirphysical properties. Finally, I vary the disc di ff usivities to get anoverview of the more exotic configurations before concluding.
2. Model
In the following, I consider a magnetised disc orbiting a cen-tral object of mass M . The disc is subject to the gravitationalpull of the central object only (self-gravity being neglected), andthe Lorentz force due to electrical currents. Disc is assumed tobe weakly ionised, hence I consider a generalised Ohm’s lawwhich include Ohmic and ambipolar di ff usivities. These di ff u-sivities are prescribed in section 2.3. In the following, I will use either spherical ( r , θ, ϕ ) coordinatesor cylindrical ( R , ϕ, z ) coordinates, depending on the context.Note that the cylindrical radius is written in upper case, whilethe spherical radius is lower case. I solve the non-ideal MHDequations, here in spherical coordinates: ∂ t ρ + ∇ · ρ v = , (1) ∂ t ρ v + ∇ · ρ vv = − ∇ ρ c s + J × B c − GM ρ r e r , (2) ∂ t B = − ∇ × E (3) E = − v × B + π c (cid:16) η O J + η H J × ˆ B − η A J × ˆ B × ˆ B (cid:17) (4)where I have defined the isothermal sound speed c s ( r ), assumingthe flow was locally isothermal, the plasma current J = c ∇ × B / π , the gravitational constant G the magnetic field direction ˆ B , the electromotive field E and the Ohmic, Hall and ambipolardi ff usivities η O , η H and η A . Note that the di ff usivities are alsofunctions of space and magnetic field strength. I stress that noturbulence is assumed in this model, only molecular magneticdi ff usivities which results from the low ionisation fraction of theplasma. If accretion is occurring, it is only the result of torqueswhich are self-consistently computed in the model. In the context of protoplanetary discs, the strong di ff usivitiesare known to almost suppress all non-axisymmetric structures(Béthune et al. 2017; Bai 2017) in the regions outside 1 AU. Forthat reason, I assume that the flow is 2.5D: I conserve three com-ponents for v and u but neglect ϕ derivatives.In order to simplify the problem even more, I will assumethe flow is globally self-similar. This implies that as one movesaway from the central object, the flow "looks" the same. This ap-proach has several advantages: first it avoids issues with the in-ner boundary conditions which are usually problematic in globalnumerical models. Second, it allows us to explore systematicallythe parameter space at exquisite resolution with a limited nu-merical cost, as the problem becomes essentially a 1D problem.I follow Ferreira & Pelletier (1993) and define the self-similarscaling for any field Q as Q ( r , θ ) = (cid:32) rr (cid:33) γ Q ˜ Q ( θ ) , (5) where γ Q is the self-similar exponent and ˜ Q is a 1D functionwhich completely determines the flow. The scaling of the grav-itational force as 1 / r imposes the self-similar scaling to all ofthe other components, i.e γ v = −
12 ; γ B = −
54 ; γ ρ = −
32 ; γ η = . Similarly to the velocity field, the sound speed c s is proportionalto r − / , and ˜ c s is a prescribed function of θ . Note also that theself-similar scaling is here written in spherical coordinates, butcan be transformed into cylindrical ( R , z , ϕ ) coordinates by not-ing that r = R sin( θ ) and z / R = tan − ( θ ): Q ( R , z ) = (cid:32) RR (cid:33) γ Q ˆ Q ( z / R )where ˆ Q ( z / R ) = sin γ Q ( θ ) ˜ Q ( θ ).In the following, it will be useful to define the Keplerian ve-locity v K ( R ) ≡ √ GM / R and the Keplerian angular frequency Ω K ≡ v K / R . I will also use the disc geometrical thickness h ≡ c s ( R , z = / Ω K ( R ). Note that by construction from the self-similar scaling, h / R ≡ ε is a constant of the problem. The di ff usivities η are usually computed from complex thermo-chemical network (e.g. Thi et al. 2019). It is customary to usedimensionless numbers to quantify the di ff usivities: the OhmicReynolds number Rm and the ambipolar Elsasser number Λ A defined asRm ≡ Ω K h η O (6) Λ A ≡ V Ω K η A (7)these dimensionless number are more useful than the more tra-ditionally used Elsasser numbers since in most cases they donot depend on the magnetic field strength (Wardle & Ng 1999).Hence, they only depend on the gas properties (density, temper-ature, composition).The use of the self-similar ansatz implies that Rm and Λ A must be functions of z / R (or θ ) only. This is not true in full chem-ical models, except for ambipolar di ff usion, which is known tobe of order of unity across a wide range of scales (see e.g. Thiet al. 2019, Fig. 8-middle panel). In the following, I will there-fore mostly focus on models dominated by ambipolar di ff usion,and introduce Ohmic di ff usion subsequently.Following Thi et al. (2019) models, I prescribe the ambipolardi ff usion profile to be Λ A ( z / R ) = Λ A0 exp (cid:34) z ( λ h ) (cid:35) (8)where λ is a free parameter quantifying the thickness of the non-ideal layer in units of h .In models including Ohmic di ff usion, I use the following dif-fusivity profileRm( z / R ) = Rm exp (cid:32) z ( λ h ) (cid:33)(cid:32) ˆ ρ ( z / h )ˆ ρ (0) (cid:33) − . This profile, and in particular the dependency on the density ra-tio, is chosen to be consistent with the ambipolar di ff usion pro-file, assuming the plasma is made of two types of charged species(such as electrons and molecular ions). Article number, page 2 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs In contrast to usual self-similar approaches (Casse & Ferreira2000a; Jacquemin-Ide et al. 2019) where stationary equation aresolved with a shooting methods through the critical points ofthe flow, here I solve the time-dependent equations (1)—(3) us-ing the code PLUTO, a finite-volume, shock-capturing scheme(Mignone et al. 2007). I use a spherical domain, with the shapeof a shell with only one grid point in the radial direction and 2048points distributed homogeneously in the θ direction. I choose thegrid to be centred on r = r =
1, and the shell extends from θ = .
15 to θ = π − .
15. The radial boundary conditions are setenforcing the self-similar relations described in section 2.2. Forthe boundary in the θ direction, I use standard outflow boundaryconditions.The initial condition is a disc in hydrostatic equilibriumthreaded by a large scale vertical ( z ) magnetic field whose ini-tial value is set to have β = in the disc mid-plane. This initialcondition is strongly unstable and the disc very quickly (in lessthan 10 Ω − ) launches an outflow before reaching a quasi steady-state.Given that the code is time-dependent and our choice ofboundary conditions, the disc mass is not necessarily conserved.Indeed, a fraction of the accreted material can be lost in the wind,leading to a slow decrease of the disc mass. In historical models,this is taken into account by slightly adjusting the power expo-nent γ ρ to ensure that mass is constant in the disc. This is notpossible in our case since the solution is dynamically evolving:one would need to change γ ρ as a function of time. Instead, Itherefore choose to renormalise the mass at each time step tokeep the total mass inside the domain constant.For similar reasons, the total magnetic flux threading the disccan evolve on secular timescales in this approach. This is notallowed in usual self-similar approaches, which usually assumethat the total toroidal electromotive field is null. Here, I choosenot to enforce such a constraint. This allows me to measure thetransport of magnetic flux, but it also implies that the magneticflux evolves with time. Hence, as for the density, I multiply theeach component of the field by a fixed factor at each time step,adjusted so that the magnetic flux threading the disc is equal tothe desired value. Several diagnostics can be derived from self-similar solution.The most useful diagnostics are related to accretion theory, i.e.how the disc surface density and mass accretion rate evolve withspace and time. Let me therefore develop the mass, angular mo-mentum and magnetic flux conservation equations as: ∂ Σ ∂ t + R ∂∂ R ˙ M acc π = (cid:104) ρ v z (cid:105) z − z , (9)˙ M acc v K π = ∂∂ R (cid:32) R (cid:90) z − z d z T R ϕ (cid:33) + R (cid:104) T z ϕ (cid:105) z − z , (10) ∂ B z ∂ t = − R ∂∂ R R E ϕ , (11)where the overline denotes an azimuthal and ensemble average(in the case of self-similar stationary solution, this average is notstrictly needed), Σ = (cid:82) z − z d z ρ is the disc surface density, z isthe disc surface, ˙ M acc ≡ − π R (cid:82) z − z d z ρ v R is the mass accretionrate, T x φ = ρ v x ( v ϕ − v K ) − B x B ϕ / π is the stress ϕ componentwhile B z = B z ( z =
0) and E ϕ are the vertical magnetic field strength and azimuthal EMF in the disc plane, respectively. Itshould be noted that in principle, accretion theory only requirethe first two relations. However, it is now well established thatin a wind-driven disc, stresses and mass loss rates are also func-tion of the mean poloidal field strength, hence, in this context,accretion theory has to be supplemented by (11).Although equations (9)—(11) fully describes the secularevolution of any disc at hand, it should be realised that the righthand-side terms are a priori unknown. This closure problem iswell known in the disc community, and is usually solved byusing the α disc paradigm in cases where only the T r ϕ term ispresent. In the case of wind-driven disc, 4 terms are actuallypresent. Following the α disc idea, I therefore define 4 param-eters: ζ ± ≡ ± ρ v z ( ± z ) ΣΩ K , (12) α ≡ (cid:82) z − z d z T R ϕ (cid:82) z − z d zP = (cid:82) z z d z T R ϕ ΣΩ h , (13) υ ± ≡ ± T z ϕ ( ± z ) ΣΩ h , (14) v B ≡ E ϕ Ω K hB z , (15)where I have assumed that the disc was isothermal to define α .With these definitions, a positive v B implies that the field is trans-ported outwards . Note that the mass loss rate ζ is comparable tothe definition of Scepi et al. (2018), up to a factor of order unity.In numerical applications, I use z = h .It is also customary to measure the magnetic field strength asa function of the plasma β parameter. In this manuscript, I willdefine the plasma β on the strength of the mean vertical fieldthreading the disc mid-plane, i.e. β ≡ πρ ( z = c s B z . (16)The transport coe ffi cients ζ , α , υ and v B being dimensionlesscoe ffi cients, they are expected to depend only on dimensionlessnumbers. We therefore expect them to depend on the magneticfield strength β , but also on the di ff usivity of the disc Rm and Λ A .In principle, once all of these dependencies are known, accretiontheory is said to be complete and one can compute the evolutionof a wind-driven protoplanetary disc in a way similar to the his-torical α disc. Therefore, determining these dependencies is themain objective of this article.
3. Fiducial simulation
In the following, I will define the fiducial run as a simulationwith ambipolar di ff usion only, Λ A = λ = h above the disc mid-plane. Thedisc aspect ratio is fixed to ε = . c s ( R , Z ) = c s ( R , z = β = . Note that I have also ran this simulation at halfresolution (i.e. 1024 points in θ ) to check the convergence ofthe solution. The two solutions di ff er by less than 1%, so I amconfident that the full resolution simulations presented here areperfectly resolved. Article number, page 3 of 13 & A proofs: manuscript no. selfsimilar
Fig. 1.
Space-time diagram showing the evolution for the first 1000 Ω − of B θ (top), B ϕ (middle) and v θ (bottom) for the fiducial run. The sys-tem quickly reaches a steady-state. In order to continue the solutionas a function of the field strength, the field is slowly increased every300 Ω − . While most of the solutions presented in this article arequasi-steady state, let me briefly discuss the space-time evolu-tion of the simulations used to compute the solutions. Such aspace-time evolution is shown in Fig. 1, for the fiducial run. Asit can be seen, the system rapidly relaxes into a steady-state. Dur-ing the first 100 Ω − , I apply a linear damping to the equations ofmotion to damp epicyclic oscillations more rapidly than the sys-tem would naturally do. These oscillations are due to the suddenlaunching of the wind, and are a spurious result of my wind-freeinitial conditions (a similar procedure was used in shearing boxby Lesur et al. 2013).After the first 100 Ω − , the system reaches a quasi-steadystate, from which I can compute the information needed: angularmomentum transport, mass flux, etc. In order to continue the so-lution as a function of β , I slowly increase the total field strengthbetween 300 Ω − and 400 Ω − , and let the system reach its newequilibrium for the next 200 Ω − . I then repeat this procedureover and over, until the code stops because of singular points inthe solution. In doing so, I reached β =
35 in the fiducial run,by steps of 10% decrease in β . In the fiducial setup, I have alsocontinued the solutions from β = up to β = by stepsof 10% increase to study the domain of validity of the solutionsobtained by this procedure. Hence, the domain explored for thefiducial setup is wider than for the other solutions (section 4).In the following, I average the quasi-steady flow obtained forthe last 100 Ω − of "relaxation" for each magnetisation and usethe result to compute the flow properties. The 2D topology of the flow can be deduced from the 1D "shell"solutions of the simulation by reconstructing the 2D fields usingthe self-similar scalings (5). Using this procedure, I compute theflow topology for the two extreme cases ( β = and β = β = ) shown in Fig. 2. In this fidu-cial simulation, the flow is symmetric with respect to the discmid-plane, as seen here. The strongly magnetised solution pro-duces a super-sonic accretion flow, while the weaker magnetisedcase tend to produce subsonic accretion ( v p / c s < . ffi cient, thanksto ambipolar di ff usion. As the field gets stronger case however,the field lines are clearly pinched around the mid-plane. This in-dicates a strong field dragging by the accretion flow.The wind is in all cases super fast-magnetosonic. The fastsurface is closer to the disc surface when the magnetisation getslarger, while the Alfvén surface tends to stay at the same alti-tude (around 4 . h ). Generally speaking, stronger magnetisationleads to faster and more massive winds (for an identical disc).Interestingly, while the shape of the stream and field lines di ff ersignificantly in the disc, their aspect is very similar in the windregion starting above the disc surface. Note that in all of thesecases, the flow show signs of collimation towards the axis. The transport properties defined in section 2.5 can be evaluatedas a function of the magnetisation for the fiducial run. This givesthe dependencies shown in Fig. 3. The resulting transport co-e ffi cients are remarkably close to power laws of β (except for v B which shows a shallower dependence on β than a power law).Several important trends can be already deduced from this result.First, it should be pointed out that since the solutions are allsteady-state and laminar, the α measured here is a pure laminarstress . It is not due to any form of turbulence, as it is just the re-sult of the large scale wind which transports angular momentum both radially and vertically. Still, it can be compared to what isfound in turbulent discs. The dependence of α on β − is steeperthan the one usually found in ideal-MHD shearing box simu-lations, which usually give α ∝ β − / (e.g. Salvesen et al. 2016;Scepi et al. 2018). Note however that the ideal and non-ideal val-ues for α match for β ∼
10. Hence, my values of α are generallysmaller than their ideal MHD counterpart by a ratio ∼ β / / ff usion in thedisc.In general, the values for α , υ and ζ are compatible with thevalues found in the shearing box literature for β ∼ − andequivalent di ff usivity regimes dominated by Λ A = β areslightly steeper than Bai (2013).Interestingly, v B > di ff using outwards . This isqualitatively consistent with Bai & Stone (2017), but the di ff u-sion speeds are reduced by an order of magnitude in my case.For β = , I get v B = . × − while Bai & Stone (2017)get v B = × − . A careful examination of the ideal and non-ideal contributions to E ϕ (Fig. 4) shows that advection almostcompletely balance di ff usion, di ff usion winning only by a 1%excess. It is also worth noting that while the ideal and ambipo-lar EMFs vary strongly in the disc, their sum is approximatelyconstant with θ , indicating that poloidal field lines are “moving”radially without much vertical deformation in the disc.Using eq. 10, the safe-similar scalings (eq. 5) and the defini-tions of section 2.5, it is straightforward to show that the massaccretion rate is directly related to α and υ ± ˙ M acc π Σ R Ω K ≡ ˙ m = ε (cid:16) αε (cid:124)(cid:123)(cid:122)(cid:125) ˙ m R + υ + + υ − ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) ˙ m z (cid:17) Note that I measure v B in units of c s while Bai & Stone (2017) quote v B in units of v K , hence a factor of ε should be added when comparingthe two results.Article number, page 4 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs Fig. 2.
Flow topology in the fiducial run (ambipolar di ff usion only Λ A = R = β = ; 10 ; 35. The green lines denotes critical lines ofthe flow: Alfvénic (plain) and fast magnetosonic (dot-dashed). The green dashed line represents the disc "surface" where the flow becomes ideal,arbitrarily located at z = . h for all of the solutions. where I have split ˙ M acc in a contribution from the radial ( ˙ m R )and vertical ( ˙ m z ) stresses. Even though these are usually referredto as "turbulent" and "wind" contributions (even in cases wherethe flow is laminar), I emphasise that there is no turbulent stress.Hence both contributions are intrinsically due to the wind. Usingthe power law scalings in Fig. 3, I find that˙ m z ˙ m R = . β . . Hence whenever β (cid:29)
1, the vertical stress is largely dominantand the radial transport of angular momentum can be neglectedaltogether, which is consistent with previous 3D numerical sim-ulations (Béthune et al. 2017). Neglecting ˙ m R , I can derive the mass accretion rate in a real-istic disc using only the scaling for υ ± , which gives˙ M acc = . × − (cid:32) Σ
10 g . cm − (cid:33) (cid:18) R
10 A . U . (cid:19) . × (cid:32) MM (cid:12) (cid:33) . (cid:18) β (cid:19) − . M (cid:12) / yr , (17)or substituting β with the field strength and the surface density˙ M acc = . × − (cid:32) Σ
10 g . cm − (cid:33) . (cid:18) R
10 A . U . (cid:19) . (cid:32) MM (cid:12) (cid:33) − . × (cid:18) ε . (cid:19) − . (cid:18) B z (cid:19) . M (cid:12) / yr . (18)It is worth noting that the mass accretion rate does not de-pend strongly on the disc surface density, in contrast to the usualviscous disc model. Instead, I find that it is the field strength that Article number, page 5 of 13 & A proofs: manuscript no. selfsimilar β − − − − − ζ + ( w i nd m a ss l o ss ) fiducial0 . β − . β − − − − − − α (r a d i a l a n g . m o m . tr a n s p o rt) fiducial40 β − . β − − − − v B ( fl u x tr a n s p o rt) fiducial7 × − β − . + 4 β − . β − − − − − − υ + ( v e rt i c a l a n g . m o m . tr a n s p o rt) fiducial1 . β − . Fig. 3.
Transport coe ffi cient as a function of the disc magnetisation for the fiducial run (ambipolar di ff usion only Λ A0 = ζ + and υ + ) are shown. The best fits using power laws are shown in plain orange line. is the dominant control parameter in these solutions. More quan-titatively, I find accretion rates compatible with expected onesfrom observations ( ∼ − —10 − M (cid:12) / yr) when β ∈ [10 , ]or B z ∼ a few mG. A word of caution though: accretion rates aretypically inferred from the material which falls onto the centralstar, while here the mass accretion rate is computed in the discat 10 AU. If the wind extract a lot material between the star and10 AU, the accretion rates have no reason to match.A common way to quantify what fraction of the accretedmass is ejected is through the ejection e ffi ciency parameter ξ de- fined as ξ ≡ π R (cid:104) ρ v z (cid:105) z − z ˙ M acc = d log ˙ M acc d log R in steady statewhere the last equality results from the continuity equation as-suming ∂ t Σ =
0. It is easy to show that the ejection e ffi ciency isdirectly connected to the transport coe ffi cients: ξ = ζ + + ζ − ε ( υ + + υ − ) = . β . where the first equality assumes ˙ m R is negligible while the sec-ond uses the scaling found in the fiducial run (Fig. 3). This im- Article number, page 6 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs − . − . − .
25 0 .
00 0 .
25 0 .
50 0 . π/ − θ − − − − − − − − − E ϕ / B z IdealAmbipolarTotal
Fig. 4.
Zoom on the toroidal electromotive field of the disc, splittingthe ideal MHD and ambipolar di ff usion contributions for β = inthe fiducial case. Note that the two electromotive fields almost exactlybalance. The total EMF at θ = π/ E ϕ / B z (cid:39) × − and is approxi-mately constant throughout the disc. This gives v B = × − . plies that ξ = O (1) for β ∼ . Therefore, the accretion rate isapproximately a linearly increasing function of radius in steady-state in these models!Finally, while most of the scalings presented here use β andtherefore the vertical field strength as the main control parame-ter, these solutions are actually dominated by the azimuthal fieldcomponent B ϕ . Using the definition of υ it easy to show that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) B ϕ B z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z = (cid:114) π βυ (cid:39) . β . , (19)where the last equality comes from the scaling of the fiducial run.Hence, for β = , | B ϕ | is more than 14 times larger than B z atthe disc surface, indicating that the field is strongly wrapped atthat location. Note finally that since the υ + = − υ − , the sign of B φ is reversed on the top and bottom side of the disc. Although it is possible to interpret accretion in terms of angularmomentum budget (the disc angular momentum is transportedinto the wind), this gives us little information about how the discis actually accreting. This question can be addressed by lookingat the forces acting on the disc. In particular, because the out-flow is magnetised, the Lorentz force F L = J × B / c is key. If F L ,ϕ <
0, the Lorentz force breaks down the disc rotation by cre-ating a negative torque, driving accretion. One easily sees thatthe magnetic torque is closely related to the poloidal current F L , ϕ = J p × B p c . Because the poloidal field is more homogeneously distributed inthe disc (Fig. 2), the magnetic torque is strongest in the regionsof strong J p . Hence the poloidal current is a direct indicator ofthe torque exerted by the wind on the disc. The distribution ofpoloidal currents in the fiducial β = case is shown in Fig. 5.I find a strongly focussed current sheet in the disc mid-plane,with currents directed outwards. It is this current which inducesa strong negative torque, and which, in turn, is responsible for Fig. 5.
Poloidal current in the fiducial simulation at β = . The cur-rent intensity is represented in colour. Note the presence of a strongoutwards current in the disc mid-plane. The green lines represents thesame characteristics as Fig. 2. − . − . . . . z/h − − π J / c J R z − z − / Fig. 6.
Zoom on the radial current in the fiducial simulation at β = .The current sheet is strongly focused, with a scaling compatible with J r ∝ z − (see text). accretion in the disc mid-plane. This current then closes back inthe disc surface and in the wind.It might be surprising at first sight to find such a narrow cur-rent sheet in the system given that the disc is strongly di ff usive( Λ A = ff usion isa self-focusing di ff usion process. It’s been argued that in "free"reconnection sites, the current scales like z − / , where z is thedistance from the null point (Zweibel & Brandenburg 1997). Inthe particular case presented here, I find that the current is even Article number, page 7 of 13 & A proofs: manuscript no. selfsimilar more focused with J r ∝ z − close to the mid-plane. Note that itis not so surprising since in the present case, we’re dealing with areconnection site "forced" by the outflow launched from the sur-face, hence the scaling need not be identical to that of Zweibel& Brandenburg (1997). Winds in protoplanetary discs are often refereed to as magneto-thermal winds. This is because most (if not all) of the globalsimulations published up to now have been using a "hot atmo-sphere". These hot atmospheres are either achieved by a pre-scribed heating function (Béthune et al. 2017), or with a morecomplete treatment of the thermodynamics including varioussources of heating and cooling (Wang et al. 2019; Gressel et al.2020). In any case, the energetics of these winds is partially dom-inated by thermal heating, hence their name.In the fiducial solution presented here, the equation of stateis locally isothermal. Hence the temperature of the atmosphere isthat of the disc at the same cylindrical radius R . For that reason,my fiducial wind is not magneto-thermal, as we shall see.Since the wind is a steady-state ideal MHD flow, several con-served quantities can be derived from the equations of motion.Let me therefore consider a streamline passing through the mid-plane at R . I define Ω as the Keplerian angular velocity at R : Ω ≡ Ω K ( R ). I then follow Blandford & Payne (1982) defining: – The mass loading parameter κ = πρ u p Ω R B p B , where v p and B p are the poloidal velocity and field strength,quantifies the amount of mass loaded onto the field lines. – The rotation parameter ω ≡ ΩΩ − κ B B ϕ πρ RR Ω which measures the rotation speed of magnetic surfaces – The magnetic lever arm λ ≡ Ω R Ω R − RB ϕ R κ B which measures the amoung of angular momentum extractedby the wind – The Bernoulli invariant B = v Ω O R (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) Kinetic − R √ R + z (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) Gravitational − ω RB ϕ κ R B (cid:124) (cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32) (cid:125) Magnetic + w Ω R (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) Thermal measures the energy content of the flow. In this definition, Ihave included the thermal contribution to the flow energeticsas the work done by thermal pressure forces along the fieldline w ≡ (cid:90) ∞ s − d (cid:96) · ∇ P ρ . If the flow is adiabatic, this is simply the enthalpy. However,the locally isothermal approximation does not lead to an adi-abatic flow, hence I have to use the integral form. − Z/R . . . . . . ω , λ κ ωλκ Fig. 7.
MHD invariants plotted along one given field line for the fiducialrun at β = . As expected, the invariants are approximately constantonce in the ideal-MHD region above the disc. The asymptotic valuesare κ = . , ω = .
99 and λ = . I have computed the mass loading, rotation and lever arm pa-rameters along a magnetic field line for the fiducial simulation at β = (Fig. 7). As expected, these parameters are constant onceabove the non-ideal region of the disc. Surprisingly, the leverarm is very small, reaching only λ = .
56. As is well known,cold MHD winds require λ > /
2, so we are just above thislimit value. A close inspection at the Bernoulli invariant (Fig. 8)shows that thermal driving is indeed negligible in the energybudget of the outflow, confirming that the solution is a cold MHDwind, not a magneto-thermal wind. It also shows that B > − Z/R − . − . − . − . . . . . . B KineticGravitationalMagneticThermalTotal
Fig. 8.
Bernoulli invariant plotted along one given field line for the fidu-cial run at β = . The Bernoulli invariant becomes constant once inthe ideal MHD region, as expected. The thermal contribution to the en-ergetics is negligible compared to the magnetic contribution, hence thewind is classified as a "cold" MHD wind. A systematic exploration of the κ, λ parameters as a functionof β in the fiducial simulation gives Fig. 9. This figure can di- Article number, page 8 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs rectly be compared to figure 2 in Blandford & Payne (1982). Itis interesting to note that these new solutions at low magnetisa-tion have a much larger κ and much smaller λ than the historicalsolutions of Blandford & Payne (1982). Such a trend was alreadyobserved by Jacquemin-Ide et al. (2019) in ideal MHD for β > λ , with an asymp-tote in the high β limit close to λ = .
4. This is smaller thanthe λ = / z =
0, while here the wind is emitted from z = H , so the lever arm is physically allowed to be slightlysmaller than 3 /
2. This is confirmed by analysing the Bernoulliinvariant at β = , which still shows that the magnetic energyis dominating the energetics at the disc surface (not shown). κ . . . . . λ l og ( β ) Fig. 9.
Dependence on κ and λ of the fiducial solutions as a functionof the disc magnetisation β . As the field strength increases, λ increaseswhile the mass loading parameter decreases.
4. Parameter space exploration
In order to explore the impact of the disk di ff usivity on the be-haviour of the outflow, I have varied the strength of ambipolardi ff usion in the disc plane Λ A0 (eq. 8), keeping the same verticalprofile. The resulting transport properties are shown in figure. 10.The general trend is that as di ff usion is increased, the mass lossrate and angular momentum fluxes decrease, while magnetic fluxtransport increases. This overall suggest that wind-driven accre-tion tends to be reduced when ambipolar di ff usion gets stronger,as expected naively. It should however be noted that all transportcoe ffi cients do not vary in the same proportion. Taking β = as the reference case, I find that a decrease of Λ A0 by a factor16 leads to a decrease of ζ , α and υ by respectively a factor 3 . . .
25, while v ψ increases by a factor 4 .
0. This thereforesuggests a relatively weak sensitivity of the transport coe ffi cienton Λ A0 (shallower than Λ . ), especially for the vertical angularmomentum coe ffi cient υ . The weak dependence of υ on the discdissipation properties was already observed in shearing box sim-ulations (see Lesur 2020, figure 39 and related text) and is hereconfirmed in global geometry. Overall this confirms that wind-driven angular momentum extraction is only moderately a ff ectedby the strength of ambipolar di ff usion, and that the mass accre-tion rate (17) is approximately valid (up to a factor of a few) forthe range of ambipolar di ff usion considered here. In addition to the strength of ambipolar di ff usion, I have also ex-plored the impact of Ohmic di ff usion on the flow topology. Inprotoplanetary discs, this leads to wind solutions valid in Gam-mie (1996) historical “dead zone”, found around 1 AU in typi-cal chemical models. To limit the computation time required bysuch models, I have chosen to limit myself to Rm >
1. In theliterature, one often finds Ohmic di ff usion evaluated in terms ofOhmic Elsasser number Λ O , which depends on the field strength.One typically has Λ = Rm β so Rm (cid:38) R (cid:46) ff usion fixed to its fiducial value Λ A0 =
1. Hence, Ohmic dif-fusion is added to the di ff usivity tensor of the fiducial run.As for ambipolar di ff usion, let me start with the transport co-e ffi cients as a function of Rm (fig. 11). Here, the angular mo-mentum transport coe ffi cients α and υ are barely a ff ected byOhmic di ff usion, even when Rm =
1. This suggests that theaccretion rate in Gammie (1996) dead zone should not be sodi ff erent from that observed in regions dominated by ambipolardi ff usion. I note however that field di ff usion rate is increased bypossibly several orders of magnitude. This is probably the moststringent e ff ect of adding Ohmic di ff usion to the system. WhileRm ≥
10 solutions are top-down symmetric (so only ζ + and υ + are plotted), the Rm = β (cid:39) . The direction of this symmetry breaking israndom. Starting from slightly di ff erent initial condition, I couldget ζ + > ζ − or the opposite, so there is no physically preferred direction for this symmetry breaking.The fact that the transport coe ffi cients are approximately un-a ff ected by Ohmic di ff usion should not lead to the wrong impres-sion that the solutions are physically identical. Indeed, a carefulexamination of the solution at Rm = ff erent from the fiducial solu-tion. First, the current sheet which was initially localised in themid-plane (Fig. 5) is now localised on the disc surfaces. Thiscan be easily interpreted as a result of Ohmic di ff usion. Indeed,while ambipolar di ff usion can have a focusing e ff ect, Ohmic dif-fusion is really only a di ff usive process. Hence, the current sheetcannot form in the strongly di ff usive mid-plane, and is thereforepushed away, at the disc surface. This change of localisation ofthe current sheet in turn impacts the accretion flow, since ac-cretion is driven by the magnetic torque created by the poloidalcurrent (see section 3.4). This is more evidently seen in Fig. 13,where the accretion flow clearly moves from the disc mid-planetowards the disc surface as Rm is decreased, thereby follow-ing the localisation of the current layer. In the end, poloidal cur-rent and accretion flows are both localised at the disc surface.However, since the current is set globally by the electrical cir-cuit enforced by the wind, its intensity is not really a ff ected bythis change in localisation, implying that total torques and accre-tion rates are not dramatically modified, as indicated by transportcoe ffi cients.Having understood that the current sheet gets divided intotwo parts and pushed towards the surface, one can easily imag-ine that such a process might not be totally symmetric for somerange of parameters. This is what happens in the Rm = β = . In this case, the current sheet is morepronounced in the southern hemisphere (Fig. 14), leading to anaccretion stream more pronounced on the bottom side of the disc,and inclined field lines in the disc bulk. I should stress that whileaccretion occurs in the southern side of the disc in this particular Article number, page 9 of 13 & A proofs: manuscript no. selfsimilar β − − − ζ + ( w i nd m a ss l o ss ) Λ A0 = 4Λ A0 = 2Λ A0 = 1Λ A0 = 0 . A0 = 0 . β − − − − α (r a d i a l a n g . m o m . tr a n s p o rt) Λ A0 = 4Λ A0 = 2Λ A0 = 1Λ A0 = 0 . A0 = 0 .
25 10 β − − − v B ( fl u x tr a n s p o rt) Λ A0 = 4Λ A0 = 2Λ A0 = 1Λ A0 = 0 . A0 = 0 . β − − − − υ + ( v e rt i c a l a n g . m o m . tr a n s p o rt) Λ A0 = 4Λ A0 = 2Λ A0 = 1Λ A0 = 0 . A0 = 0 . Fig. 10.
Transport coe ffi cients as a function of the disc magnetisation and mid-plane Elsasser number Λ A0 . example, the wind is more pronounced on the northern side ofthe disc since | ζ + | > | ζ − | and υ + > υ − . Hence angular momentumand mass flow away on the side opposite to that of accretion.Note that such a behaviour was also observed by Béthune et al.(2017) in dissymmetric solutions.
5. Discussion and Conclusions
In this work, I have used a finite volume code to derive self-similar wind solutions applicable to the di ff usive regime of pro-toplanetary discs. In contrast to previous self-similar solutionspublished to date, no assumption is made on flux stationarityor top / down symmetry. It also circumvent the di ffi culties of full-blown 3D numerical simulations which are often limited by theirintegration time and inner boundary conditions. Using this tech- nic, it is possible to explore systematically a wide range of pa-rameters at a reduced numerical cost.Using this tool, I have presented a series of wind solutionsranging from β = to β =
35 valid in the non-ideal regions ofprotoplanetary discs R (cid:38) always exist andare unavoidable as soon as a large scale magnetic field threadsthe disc. This statement is true even for β = fields whichcorrespond to a few µ G at 10 AU for typical surface densities.They give accretion rates ˙ M acc (eq. 18) which depends essen- Article number, page 10 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs β − − − − ζ ( w i nd m a ss l o ss ) fiducialRm = 10 Rm = 10Rm = 1 (top)Rm = 1 (bottom) β − − − α (r a d i a l a n g . m o m . tr a n s p o rt) fiducialRm = 10 Rm = 10Rm = 1 β − − − v B ( fl u x tr a n s p o rt) fiducialRm = 10 Rm = 10Rm = 1 β − − − − υ ( v e rt i c a l a n g . m o m . tr a n s p o rt) fiducialRm = 10 Rm = 10Rm = 1 (top)Rm = 1 (bottom) Fig. 11.
Transport coe ffi cients as a function of the disc magnetisation and mid-plane Reynolds number Rm . Note that these solutions also includeambipolar di ff usion with Λ A0 =
1. The fiducial run is physically equivalent to Rm = ∞ . tially on the field strength, and which are broadly speaking com-patible with observed accretion rates onto T-Tauri stars providedthat 10 < β < or equivalently B z ∼ a few mG. Interestingly,˙ M acc depends only weakly on the disc surface density, in contrastto viscous disc models. However, the direct comparison betweenthe accretion rate in the disc bulk (as measured by eq. 18) andthe mass accretion measured onto the star is probably mislead-ing. The latter is probably significantly smaller than the formerbecause the ejection index (see below) is of order one, hence the disc mass accretion rates could be significantly larger (possiblyby an order of magnitude) than the ones quoted here.The ejection index ξ , which essentially quantifies the ratio ofejected to accreted mass, depends only weakly on the magneticfield strength and is of the order of (but slightly less than) unity for the range of magnetisation I have explored. This means that the mass loss rate is comparable to the mass accretion rate ata given radius. Consequently, if one assumes steady state, theaccretion rate is expected to increase significantly with radius.As expected from a high ξ wind (Lesur 2020, eq. 11.22), thewind lever arm is always relatively small with typically λ < β > , without requiring any heating. Such a high ξ and low λ was already found in previous numerical work (Béthune et al.2017; Bai 2017), and it was initially thought that atmosphericheating, included in all these models, was the main reason forthis result, following the argument of Casse & Ferreira (2000b).I show here that heating is actually not key to get these high ξ -low λ solutions, a result also reported recently by Jacquemin-Ide et al. (2019) in the context of fully ionised discs. Indeed, the Article number, page 11 of 13 & A proofs: manuscript no. selfsimilar
Fig. 12.
Flow structure for the solution Rm = β = . Left: streamlines and sonic Mach number. Middle: field lines and density map. Right:Poloidal current lines and current density. Compared to the fiducial case, the current layer is now localised at the disc surface, leading to a surfaceaccretion flow. − . − . . . . z/h − − − − − ρ v R / ρ m i d c s Rm = 1Rm = 10fiducial − . − . . . . z/h − − − − − v R / c s Rm = 1Rm = 10fiducial Fig. 13.
Inward mass flux (top) and radial velocity Mach number (bot-tom) as a function of Ohmic di ff usion strength at β = . Note thatthe accretion flow moves towards the disc surface as Ohmic di ff usion isincreased, following the poloidal current sheet (see text). solutions presented here are locally isothermal and I have shownthat thermal heating is negligible in the wind energy budget evenfor β = . Hence, low λ and high ξ are quite clearly signaturesof weakly magnetised ( β (cid:29)
1) outflows.The vertical field strength predicted to get accretion ratescompatible to typical T-Tauri stars (about a mG at 10 AU, seeeq. 18) are compatible with recent upper limits obtained fromZeeman measurements (e.g. Vlemmings et al. 2019). When itcomes to field detection, it is quite clear that the toroidal compo-nent is a better candidate than the vertical one since the formeris expected to be much stronger (by a factor 10 or so at β ∼ )than the latter (eq. 19). Note however that the toroidal compo-nent is expected to change sign across the mid-plane, so this de-tection is possible only for tracers which do not average out thefield through the disc thickness.When it comes to the disc microphysics, I find that the disc-averaged mass accretion & ejection rates depend only weakly onambipolar and Ohmic di ff usivities. This implies that one doesnot need very elaborated ionisation models in the disc to get theright accretion / ejection bulk properties. However, the disc ver-tical structure actually strongly depends on di ff usion. Solutionswith large Ohmic di ff usion (i.e. valid closer to 1 AU) tend to ex-hibit accretion at the disc surface while weaker Ohmic di ff usion( R >
10 AU) show accretion in the disc mid-plane. Some of thestrong Ohmic di ff usion solutions also exhibit top / down dissym-metry, where accretion mostly occurs mostly on one side of thedisc. Such a dissymmetry was also found in simulations includ-ing Ohmic and ambipolar di ff usion (Béthune et al. 2017; Gresselet al. 2020) in similar range of parameters, so this result is notunheard of. These topology properties will likely have a strongimpact on dust growth and planet migration theory, but again,they barely a ff ect the vertically-averaged transport properties.Finally, In all of these models, the large scale magnetic fieldis found to be transported outwards , in agreement with Bai &Stone (2017). However, the transport velocity is about an orderof magnitude lower than that measured by Bai & Stone (2017) insimilar conditions (ambipolar di ff usion only). Moreover, I findthat the transport rate varies significantly with the strength ofOhmic di ff usion and to a lesser extent with Ambipolar di ff usion.Overall, This dependency might explain why simulations and Article number, page 12 of 13eo ff roy R. J. Lesur : A systematic description of wind-driven protoplanetary discs Fig. 14.
Flow structure for the solution Rm = β = . Left: streamlines and sonic Mach number. Middle: field lines and density map. Right:Poloidal current lines and current density. . This solution exemplify the dissymmetric solutions found around β = at Rm = models tend to disagree on the magnetic transport rate (Bai &Stone 2017; Leung & Ogilvie 2019; Gressel et al. 2020).Clearly, the large scale field transport is the main bottleneckof a complete theory of wind-driven accretion since the strengthof the large scale poloidal field is the main control parameter.This cannot be ignored for the secular evolution of discs sincethe field strength is expected to vary with time. Indeed, whileclass II objects seem to require field strengths of a few mG at 10sof AU to get the right accretion rates (see above), core collapsecalculations tend to suggest a field strength of 100 mG when thedisc forms (e.g. Masson et al. 2016). Hence, the field strengthmust be reduced by at least one order of magnitude (more prob-ably two) between class 0 and class 2, indicating that flux trans-port must be relatively e ffi cient. The problem of field transportshould therefore deserve more attention in the future if winds areto be used in secular evolution models. Acknowledgements.
I am thankful to Jonatan Jacquemin-Ide, Etienne Martel,Jonathan Ferreira and Cornelis Dullemond for fruitful discussions regarding thephysics of self-similar solutions. I acknowledge support from the European Re-search Council (ERC) under the European Union Horizon 2020 research andinnovation programme (Grant agreement No. 815559 (MHDiscs)). All of thecomputations presented in this paper were performed using the GRICAD infras-tructure (https: // gricad.univ-grenoble-alpes.fr), which is supported by Grenobleresearch communities. References