On the maximum grain size entrained by photoevaporative winds
MMNRAS , 1–11 (2016) Preprint 1 October 2018 Compiled using MNRAS L A TEX style file v3.0
On the maximum grain size entrained by photoevaporativewinds
Mark A. Hutchison, (cid:63) Guillaume Laibe, and Sarah T. Maddison Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We model the behaviour of dust grains entrained by photoevaporation-driven windsfrom protoplanetary discs assuming a non-rotating, plane-parallel disc. We obtain ananalytic expression for the maximum entrainable grain size in extreme-UV radiation-driven winds, which we demonstrate to be proportional to the mass loss rate of the disc.When compared with our hydrodynamic simulations, the model reproduces almost allof the wind properties for the gas and dust. In typical turbulent discs, the entrainedgrain sizes in the wind are smaller than the theoretical maximum everywhere but theinner disc due to dust settling.
Key words: planets and satellites: atmospheres – protoplanetary discs – circumstellarmatter – stars: pre-main-sequence
Small dust grains in the upper atmosphere dominate theopacity for the disc at many wavelengths, thereby shieldingthe bulk of the disc from energetic radiation from the starand controlling the disc’s thermal/geometric structure (Cal-vet et al. 1991; Chiang & Goldreich 1997; D’Alessio et al.1998). Imaging scattered starlight and thermal re-emissionof absorbed stellar radiation from dust in these upper lay-ers is still a vital diagnostic tool used to characterise discsand their structure (Watson et al. 2007; Andrews 2015). Ittherefore follows that physical processes that affect the dy-namics of these grains (e.g. settling, grain growth, and discwinds) may have an impact in the way that we interpretobservations of protoplanetary discs (Testi et al. 2014).Aerodynamic drag from disc winds can loft dust into theatmospheres surrounding discs. Using order-of-magnitudeforce balance arguments, Takeuchi et al. (2005) estimate themaximum grain size that can be carried out by photoevapo-rative winds. Better estimates were obtained by Owen et al.(2011) who test a posteriori whether dust can be entrainedalong gas streamlines in single-phase photoevaporation sim-ulations. More recently, we have performed fully coupled, gasand dust hydrodynamic simulations of protoplanetary discsundergoing dust settling and extreme ultraviolet (EUV) in-duced photoevaporation (Hutchison et al. 2016, hereafter,HPLM16). Based on the suite of simulations for that study,we concluded that only micron sized dust grains and smallerare entrained by photoevaporative winds in typical discs (cid:63)
E-mail: [email protected] found around T Tauri stars. The exact cutoff, however, wasfound to depend on stellar mass, stellar irradiation flux, gasdensity at the base of the flow, and distance from the centralstar.Numerical simulations of gas and dust potentially pro-vide one of the best windows on dust dynamics in discs, butowing to the numerical difficulty associated with simulatingsmall dust grains in such steeply stratified atmospheres, theyare still too unwieldy to use in a practical sense for globaldisc studies across multiple systems. A nice alternative tousing numerical simulations is the self-similar solution forthermal disc winds derived by Clarke & Alexander (2016).In this study, we provide another alternative by deriving aneasy to use (semi-)analytic solution that recovers the major-ity of the results from our hydrodynamic simulations.The paper is organised as follows: in Section 2 we derivethe equations for our semi-analytic model and compare themodel with hydrodynamic simulations; in Section 3 we useour model to explore different parameters that affect dustentrainment in disc winds; in Section 4 we discuss the effectsof settling; and in Section 6 we summarise our findings.
Previously, Hutchison & Laibe (2016) derived an analytic so-lution for EUV-driven winds assuming a non-rotating, plane-parallel atmosphere. The simple geometry makes the prob-lem tractable, retains the vertical disc structure, and repro-duces the vertical winds near the ionisation front. Later,HPLM16 showed using hydrodynamic simulations that the c (cid:13) a r X i v : . [ a s t r o - ph . S R ] A ug Hutchison, Laibe, & Maddison back reaction on outflowing gas due to entrained dust grainsis negligible due to the small dust-to-gas ratios in the up-per atmospheres of discs. We exploit this fact to extend ourmodel to two fluids by directly inserting the analytic windsolution for the gas into the fluid equations for dust.
We assume an isothermal thin disc supported by pressure-gravity balance due to the vertical component of gravityfrom a central star g = − G Mz ( R + z ) / ˆz , (1)where G is the gravitational constant, M is the mass ofthe central star, and z is the height above the midplane.In this geometry, the parameter R and the variable z makeup a quasi-2D coordinate system centred on the star. Thewind speed for isothermal photoevaporation can be writtenin closed form using the Lambert W function (Corless et al.1996; Veberiˇc 2012), v g = c s (cid:115) − W (cid:20) − exp (cid:18) − G Mc √ R + z − (cid:19)(cid:21) , (2)where v g is the gas velocity and c s ≈
10 km s − is the isother-mal sound speed of the wind. The gas density is related tothe velocity via the relation˙ m g = ρ g v g , (3)where ˙ m is the constant mass loss rate per unit area of theoutflow in a stationary regime.The value of ˙ m is best determined by the fluid quanti-ties at the base of the flow. The gas velocity in the wind iswell constrained by Equation (2), so this amounts to deter-mining ρ g,i . For simplicity, we will assume the density in thedisc and wind is piecewise continuous, but the reality is thatcollisional heating from the ionised wind will distort the den-sity structure of the neutral disc near the ionisation front,causing the density to decrease faster than if photoevapo-ration was not present. We have not performed an exhaus-tive study of how ρ g,i changes with stellar and disc param-eters, but the suite of simulations performed by HPLM16show that collisional heating from the wind causes the ini-tial outflow to level off at ∼
40% of the assumed ionisationfront density. Neglecting this density offset equates to over-estimating dust entrainment in the wind (see Section 3.3).However, HPLM16 also showed that a non-rotating, plane-parallel atmosphere underestimates dust entrainment by ap-proximately the same amount. As a result, we only worryabout scaling the density when comparing our model di-rectly with numerical simulations.
In the limit of small dust-to-gas ratio, the steady-state fluidequations for the dust are simplified to the following equa-tions: ∇ · ( ρ d v d ) = 0 , (4) v d · ∇ v d = 1 t s ( v g − v d ) + g , (5)where t s is the Epstein drag stopping time (Epstein 1924), t s = (cid:114) πγ ρ grain sc s ρ g = ρ eff sc s ρ g , (6)and s and ρ grain are respectively the intrinsic size anddensity of the individual dust grains. For convenience, wesimplify the expression in the second equality by defining ρ eff ≡ ρ grain (cid:112) πγ/ m d = ρ d v d . (7)Meanwhile, the momentum equation (5) reduces to a singleordinary differential equation v d d v d d z = 1 t s ( v g − v d ) − G Mz ( R + z ) / , (8)where v g is a function of z and is given by Equation (2).Equation (8) can be written in dimensionless form using theparameters ¯ v ≡ v/v K , ¯ z ≡ z/R , and v K = (cid:112) G M/R :¯ v d d¯ v d d¯ z = S t − (cid:18) − ¯ v d ¯ v g (cid:19) − ¯ z (1 + ¯ z ) / , (9)where, by analogy to dust dynamics in a disc, we have de-fined the Stokes number of the wind to be S t ≡ ρ eff sv c s ˙ m g R = v /Rv g,i /t s,i , (10)where v g,i , ρ g,i , and t s,i = ρ eff s/c s ρ g,i denote the gas velocity,the gas density and the stopping time at the base of the flow,respectively. We emphasise that the Stokes number in thewind, denoted here as S t , is a priori different to the Stokesnumber in the disc (St). Equation (10) can be seen as theratio between the gravitational force and the force requiredto keep the grains entrained at the base of the flow. When ¯ z → ∞ , the term − ¯ z/ (1 + ¯ z ) / → v g → v g / d¯ z →
0. Upon ap-plying these limits to Equation (9), the solution ¯ v d (cid:39) ¯ v g with a vanishing derivative can readily by seen by inspec-tion. Note that this holds for all Stokes numbers. On theother hand, when ¯ z → ¯ z i , where ¯ z i is the location of theinitial flow, the Stokes number determines the behaviour ofthe solution. Possible solutions can be categorised into twomain classes based on whether the dust velocity is initiallyincreasing or decreasing. Increasing:
For S t → v d / d¯ z > z . However, to keepthe drag term from becoming unbounded, the gas and dustvelocities must be approximately equal, to zeroth order in S t . Two distinct subclasses of grains are possible:(i) perfectly -entrained grains that adhere to the zerothorder approximation, and(ii) well -entrained grains that do not. MNRAS , 1–11 (2016) n the maximum entrainable grain size Decreasing:
For S t → ∞ (i.e. low drag), the solution sat-isfies, to zeroth order in S t − ,¯ v d ∂ ¯ v d ∂ ¯ z (cid:39) − ¯ z (1 + ¯ z ) / . (11)Fully entrained flows require ¯ v d >
0, but allowing ¯ v d andd¯ v d / d¯ z to be general yields a total of three new subclasses:(iii) weakly -entrained dust grains with positive velocitiesthroughout the flow (i.e. away from the midplane),(iv) partially -entrained dust grains whose velocitieschange sign in the flow, and(v) non -entrained dust grains whose initial velocities arenegative (i.e. toward the midplane).Note that these latter two subclasses are best interpretedwith an Eulerian perspective, since at large ¯ z the velocitieswill always become positive and converge to the gas veloc-ity. Although a steady state flow is impossible to achieveif ¯ v d reverses direction, this would imply a pile-up of somekind for partially-entrained dust grains. Moreover, the dis-tinction between weakly- and partially-entrained dust grainssuggests there exists a critical Stokes number S t c for whichthe minimum of ¯ v d is zero. This ensures that when S t < S t c particles are entrained by the photoevaporative wind. We denote ¯ z c as the height at which ¯ v d = d¯ v d / d¯ z = 0. Sub-stituting these values into Equation (9) yields the followingrelation, S t c ¯ z c (1 + ¯ z ) / = 1 . (12)Taking the derivative of Equation (12) with respect to ¯ z c removes the S t c dependence and allows us to solve for ¯ z c ,¯ z c = ± √ . (13)This is also the location for the peak gravitational force.Substituting ¯ z c back into Equation (12) gives us the criticalStokes number, S t c = 3 √ (cid:39) . . (14)From Equation (10), we can then solve for the maximumentrainable grain size in the wind, s max = 3 √ c s ˙ m g Rρ eff v , (15)where ˙ m g is obtained from Equation (3) using initial condi-tions at the base of the flow. The critical height ¯ z c is closeenough to the disc surface for the plane-parallel approxima-tion to remain valid, but is subject to the validity of theassumptions made about the underlying disc.We can alternatively derive this limit from a Lagrangianperspective by retaining the d v d / d t term in Equation (8) andrewriting the equation using an effective potential, V eff ≡ − (cid:18) ¯ z S t + 1 √ z (cid:19) , (16)such thatD ¯ z D¯ t = − S t ¯ v g D¯ z D¯ t − d V eff d¯ z , (17) where the convective derivative (D/D¯ t ≡ ∂/∂ ¯ t + ¯v · ∇ ) istaken with respect to the dimensionless time variable, ¯ t ≡ Ω K t , with Ω K = v K /R . Except when S t > S t c , the neweffective potential monotonically decreases to infinity (i.e.dust grains are entrained and escape the system). Above S t c a local minimum forms that keeps grains with s > s max bound to the disc. The important point here is that V eff isindependent of the functional form of v g as long as v g ∝ /ρ g , or equivalently, as long as the continuity equation isvalid. We must integrate Equation (9) numerically to find the dustvelocity, but this requires initial conditions for the flow, in-cluding the gas and dust structure in the underlying disc.Equation (2) specifies the gas velocity for all z , but doesnot identify where the flow begins. Physically, the ionisationfront location marking the base of the flow is set by the in-tensity of the impinging radiation field and its correspondingoptical depth. For maximum flexibility, we leave ρ g,i as aninput parameter for our model and parameterise the pene-tration depth of the EUV radiation by defining ξ ≡ ρ g,i /ρ g,0 ,where ρ g,0 is the local midplane density of the disc. The lo-cation of the ionisation front is then obtained by solving forthe height z i at which the disc density is equal to ρ g,i .For the density profile of the disc we use the isothermalthin disc approximation ρ g ( z ) = ρ g , exp (cid:20) − z H (cid:21) , (18)where H is the local scale height of the disc. We specify theentire R - z disc structure using a power-law parameterisation(see, e.g., Laibe et al. 2012)Σ g = Σ g , (cid:18) R (cid:19) − p , (19) H = H (cid:18) R (cid:19) / − q/ , (20) ρ g,0 = Σ g √ πH , (21)where Σ g is the local surface density for the gas while quan-tities with the subscript 1 AU are reference values measuredat R = 1 AU. The parameters p and q are power-law expo-nents controlling the density and temperature (i.e. flaring)of the disc, respectively. Observations and simulations indi-cate that p and q can cover a range of values – q being thebetter constrained out of the two (e.g. Dutrey et al. 1996;Andrews & Williams 2005, 2007; Laibe et al. 2012; Pinte& Laibe 2014). In keeping with these studies, we adopt theranges p ∈ [0 , .
5] and q ∈ [0 . , . ε = 0 .
01. While un-physical, distributing the dust in this manner makes it mucheasier to compare different grain sizes, identify trends, andgain valuable insight about how dust behaves in discs. Theonly caveat is that we partially lose the ability to make pre-dictions about dust in real outflows, particularly for largergrains which we would expect to be concentrated closer tothe midplane and away from the ionisation front. In fact,
MNRAS000
MNRAS000 , 1–11 (2016)
Hutchison, Laibe, & Maddison from HPLM16 we know that settling can keep dust grainsthat are at least a few times smaller than s max from everentering the wind. Correctly accounting for this behaviourrequires hydrodynamic simulations or assuming a steady-state stratified disc structure for the dust (see Section 4).One may intuitively expect the initial dust velocity tostart from rest, but this is physically inconsistent with therelation ρ d = ˙ m d /v d , which would create an unphysical den-sity spike at the surface of the disc. The numerical modelsfrom HPLM16 show that the ionisation front is a very nar-row, dynamically complicated transition region where nei-ther phase is particularly well represented by the solutionsabove. Once the flow has settled, the dust already has a fi-nite velocity and the solution is better represented using azero derivative at ¯ z = ¯ z i . Substituting this into Equation (9)and solving for ¯ v d results in the following initial conditionfor the dust¯ v d,i = ¯ v g,i (cid:34) − S t ¯ z (1 + ¯ z ) / (cid:35) . (22) s max We verify that s max is valid by comparing the semi-analyticmodel here with the hydrodynamic model in HPLM16. Indoing so, we adopt the following values for our disc: R =5 AU, M = 1 M (cid:12) , Σ g , = 100 g cm − , H = 0 .
05 AU, p = 1, and q = 0 .
5. We assume the EUV penetration depthis ξ = 10 − , which is consistent with thermo-chemical mod-els of typical T Tauri discs (Woitke et al. 2016). In theionised wind, gas particles are held isothermally (i.e. γ = 1)at T = 10 K such that c s ≈
10 km s − . Finally, we as-sume that the intrinsic dust density for all our dust grains is ρ grain = 3 g cm − . Plugging these values into Equation (15)give s max (cid:39) . µ m. Except where noted otherwise, we willuse the fiducial values above for the remainder of the paper.Although the outflow in this geometry is 1D in nature,the transition between disc and outflow is better captured inmulti-dimensional simulations. This is because ionisation ofparticles at the disc surface causes compression of adjacentneutral particles, which in 1D produces sporadic bursts ofoutflow. Thus for our numerical simulation, we place 200 028particles on a uniform (staggered) lattice inside a Cartesianbox, ( x, z ) ∈ [ − . , .
9] AU and set the gas/dust masses anddust fraction using the method described in HPLM16 forunequal-mass, one-fluid particles. We use periodic boundaryconditions in x and dynamic boundaries in z , consistent witha steady-state wind flowing away from the disc (Hutchison& Laibe 2016). We create the dynamic boundaries by con-verting the ionised particles at t = 0 into ghost particlesand forcing them to move in the vertical direction at thelocal wind speed prescribed by Equation (2). The numberof ghost particles produced by this setup is 5372 on eachside of the disc. Photoevaporation is created by heating gasparticles to 10 K when ρ g ≤ ρ g,i = ξρ g,0 . The disc beginsthe simulation in isothermal hydrostatic equilibrium, but isevolved adiabatically from t = 0 to capture the collisionalheating at the ionisation front.We measure the reduction in ρ g,i caused by ionisationfront heating in the simulation to be ∼ s max from 0 .
82 to 0 . µ m. Note, we need not worryabout specifying ρ d,i because it has no influence in deter- v [ k m / s ] s = 0.01 µ m s = 0.05 µ m s = 0.1 µ m s = 0.2 µ m s = 0.3 µ m -202-10 -5 0 5 10 z [AU] s = 0.4 µ m Figure 1.
Dust velocities as a function of z for the 6 indicatedgrain sizes. The coloured points are results from hydrodynamicsimulations at t = 80 yr while the black lines are the correspond-ing semi-analytic results. The simulations confirm that s max is be-tween 0 . . µ m, the predicted value being 0 . µ m. Note tran-sient oscillations appear in the numerical solution, but alwaysoscillate closely about the semi-analytic solution. (An animatedversion is playable in the online pdf.) mining the value of s max . Figure 1 shows a snapshot ofthe hydrodynamic velocities at t = 80 yr overlaid with thesemi-analytic curves assuming a 40% reduction in gas den-sity in the flow. The hydrodynamic solutions are somewhatnoisy due to fluctuating motions in the flow produced bythe stochastic ionisation of gas particles at the ionisationfront, but they always oscillate about their respective semi-analytic solution. The trend of v g → s → s max is a clearindicator that Equation (15) is correct. As further proof,when we try s = 4 µ m, no dust is entrained in the outflow.We emphasise that s max is a robust, physical limit setby gas properties in the wind and is not affected by anyproperties of the disc. In fact, the above simulations wererun with dust settling enabled to show that dust density inthe disc has no effect on the entrainment properties in thewind (unless of course the density is zero, see Section 4). Be-cause it does not matter what physical mechanisms supplythe dust to the disc surface (e.g. turbulence, migration, ac-cretion jets), s max is model independent. This is importantbecause any observational constraint on s in the wind canbe converted into a strict lower bound on ˙ m g by invertingEquation (15). In terms of surface density, this translatesinto a photoevaporative mass loss rate of˙Σ g,photo ≥ π √ ρ eff s obs v c s , (23)where s obs is the largest grain size observed in the flow. In this section, we use Equations (9) and (15) to investigatehow s max and dust entrainment depend on the model’s discand stellar parameters by systematically varying the grainsize, disc radius, base flow density, and stellar mass. MNRAS , 1–11 (2016) n the maximum entrainable grain size Using the values above, we solve Equation (9) for six differ-ent grain sizes, s = [0 . , . , . , . , . , . µ m, and plotthe resulting velocity and density profiles for both gas anddust at R = 5 AU as a function of z in Figure 2. There isa steady decline in entrainment with increasing grain size.We have selected these sizes in order to have at least onerepresentative sample from each of our earlier defined en-trainment subclasses: Perfect:
The 0 . µ m grains mirror the gas velocity and den-sity profiles almost exactly. Well:
The 0 . µ m grains always have a positive accelera-tion, but noticeably deviate from the gas solution. Weak:
Both the 0 . . µ m grains exhibit a tell-taledip (peak) in their velocity (density) profiles. Partial:
The velocity for the 0 . µ m grains goes negative,causing them to stall above the disc surface. Theycannot escape or set up a steady-state outflow. Non:
The initial velocity of the 1 . µ m grains is negativeso that not even partial ejection can be achieved.Also evident in Figure 2 is the fact that all dust velocitiesconverge to that of the gas for large z , regardless of grainsize. This is even true – albeit unphysical – for partially-and non-entrained dust grains at sufficiently large z (seethe dotted lines). Finally, at this radius, s max ≈ . µ mand we have verified that v d = d v d / d z = 0 at z c = R/ √ z c where aerodynamic drag takes over as thedominant force, the larger grains experience a greater accel-eration due to their large differential velocities with the gas.As a result, their dust densities drop more rapidly than smallgrains and a seesaw pattern develops with a common pivotpoint. Small, perfectly-entrained grains that trace the den-sity profile for the gas represent the smallest (largest) possi-ble density that can be achieved by the dust at small (large) The common pivot point is actually a coincidence. In generalcrossings are well localised, but occur at different locations. z [AU] v [ k m / s ] s = 0.01 µ m s = 0.1 µ m s = 0.4 µ m s = 0.7 µ m s = 0.9 µ m s = 1.2 µ m Gas0 5 10 15 20 25 z [AU] -20-19-18-17 l og ( ρ ) [ g / c m ] Figure 2.
Velocity (top) and density (bottom) profiles for gas(dashed) and dust (solid) in a photoevaporative wind at R =5 AU. Grains with s (cid:38) . µ m cannot be entrained in the wind,but may lead to a pile-up near the surface of the disc. Largevelocity differences between phases is an indication of the lack ofdust entrainment. The velocity-density relation in Equation (7)leads to a seesaw pattern in the density. z . Thus, assuming (cid:15) is constant in the wind is a good, fair,and poor approximation for perfectly-, well-, and weakly-entrained dust grains, respectively.At the other extreme, the sign reversal in the velocityfor grains s ≈ . µ m suggests a pile-up occurs just abovethe disc’s surface (below 1 µ m the initial velocities are neg-ative). The opacity created by the structure in the densityprofile of the wind will affect the flux of radiation throughthe disc’s atmosphere. The feedback that this will have onphotoevaporation is complicated and requires proper radia-tive transfer calculations. Furthermore, the thermal emissionand/or scattered light from the dust grains in the flow canhave observational signatures unique to photoevaporatingdiscs, as shown by Owen et al. (2011). We leave this forfuture study as this goes beyond the scope of this paper. MNRAS000
Velocity (top) and density (bottom) profiles for gas(dashed) and dust (solid) in a photoevaporative wind at R =5 AU. Grains with s (cid:38) . µ m cannot be entrained in the wind,but may lead to a pile-up near the surface of the disc. Largevelocity differences between phases is an indication of the lack ofdust entrainment. The velocity-density relation in Equation (7)leads to a seesaw pattern in the density. z . Thus, assuming (cid:15) is constant in the wind is a good, fair,and poor approximation for perfectly-, well-, and weakly-entrained dust grains, respectively.At the other extreme, the sign reversal in the velocityfor grains s ≈ . µ m suggests a pile-up occurs just abovethe disc’s surface (below 1 µ m the initial velocities are neg-ative). The opacity created by the structure in the densityprofile of the wind will affect the flux of radiation throughthe disc’s atmosphere. The feedback that this will have onphotoevaporation is complicated and requires proper radia-tive transfer calculations. Furthermore, the thermal emissionand/or scattered light from the dust grains in the flow canhave observational signatures unique to photoevaporatingdiscs, as shown by Owen et al. (2011). We leave this forfuture study as this goes beyond the scope of this paper. MNRAS000 , 1–11 (2016)
Hutchison, Laibe, & Maddison z [ AU ] s = 0.01 µ ms = 0.7 µ m5 10 15 20 25 R [AU] z [ AU ] -20-19.5-19-18.5-18 l og ( ρ du s t ) [ g / c m ] Figure 3.
Density contours for s = 0 .
01 and 0 . µ m dust grainsare created by horizontally stacking 1D calculations at differentradii. Empty regions, like the one in the bottom panel at R (cid:46) s > s max . When s ∼ s max , the maximumdust density in the wind occurs at z c (dashed orange line) ratherthan along the ionisation front (thick grey line). This suggeststhat the opacity in the wind may not monotonically decreasewith z . The solid black line ( z = R ) and hash marks indicatethe region where lack of radial pressure gradients and centrifugalmotion cause our approximation to break down. The assumed underlying 2D disc structure makes our cal-culations radially consistent at the base of the flow so wecan approximate the 2D dust density in the wind by hori-zontally stacking vertical density maps from different radii.Figure 3 compares density contours for two different grainsizes, s = 0 .
01 and 0 . µ m. However, beyond z ∼ R differen-tial pressure and centrifugal effects become significant andour approximation breaks down (Hollenbach et al. 1994).The hash marks on the black line, z = R , indicate the re-gion where our approximation likely breaks down. In con-trast to smaller grains, whose density in the wind monotoni-cally declines with z , large grains with s ∼ s max tend to havestrongly peaked dust densities with maxima near z ∼ z c .Figure 3 shows that the entrainment region is not thesame for all dust grains. This radial size sorting of dust,first pointed out by Owen et al. (2011), is nicely picked upby the analytic expression for the maximum grain size inEquation (15). A better illustration of this radial sortingis obtained by plotting s max as a function of R , as shownin Figure 4 for five disc profiles varying p ∈ [0 , .
5] whileholding q = 0 . q ∈ [0 . , . p = 1. Physically, the exact shape of s max isdetermined by the relative rate of decline between gravity R [AU]
10 100 m a x s [ µ m ] (0.0,0.5)(0.4,0.5)(0.8,0.5)(1.0,0.8)(1.0,0.7)(1.0,0.6)(1.0,0.5)(1.0,0.4)(1.1,0.5)(1.5,0.5) s = 0.5 µ m s = 0.4 µ m s = 0.3 µ m s = 0.2 µ m Figure 4.
The maximum entrainable grain size plotted as a func-tion of radius for the different surface density and temperaturepower law exponents ( p, q ) labeled next to each curve. Arrowsindicate the width of the entrainment region for a few grains inthe ( p, q ) = (1 . , .
5) disc. and density as a function of R and z . Thus it comes as nosurprise that s max has a strong dependence on both p and q ,particularly at large radii where they affect the disc structurethe most.Usually grains with s < max( s max ) have a unique in-ner and outer entrainment radius beyond which they can-not be dragged into the flow, but there are a few exceptions.First, because photoevaporation cannot operate below theso-called ‘critical-radius’ (Liffman 2003; Adams et al. 2004;Alexander et al. 2014), R c,EUV (cid:39) . (cid:18) MM (cid:12) (cid:19) AU , (24)all grains with s < s max will share the same inner entrain-ment radius. A similar situation occurs for the outer entrain-ment radius if the disc is truncated at the outer edge (e.g.by external photoevaporation; see Facchini et al. 2016). Forcompleteness’ sake, we also mention that s max has no ex-tremum when q ≥ p , suggesting there is no outer entrain-ment radius for grains. However, the reality is that discs arefinite and s max will eventually drop to zero regardless. The penetration depth has a strong effect on dust entrain-ment by determining ρ g,i and z i . Rigorously, ξ also influ-ences v g,i because the initial velocity is obtained by inte-grating backwards from the sonic point to z i . However, thedisc density profile is so steep and the velocity profile is so MNRAS , 1–11 (2016) n the maximum entrainable grain size -12.2 -12 -11.8 -11.6 -11.4 log( ; g ; ) [g = cm ] s m a x [ m ] log( ) s m a x [ m ] Figure 5.
The maximum entrainable grain size at R = 5 AU asa function of ρ g,i (holding Σ g , = 100 g/cm ; blue) and ρ g,0 (holding ξ = 10 − ; orange). This confirms the linear density-entrainment relationship found by HPLM16 and shows that it isrobust against changes to the disc profile. shallow that the actual impact of ξ on v g,i is small. Because ξ is directly proportional to the initial base flow density –and since ρ g and s are inversely proportional in the Epsteindrag regime – ξ scales the entrainment properties linearlyin s , as shown by the orange dashed line and axes in Fig-ure 5. Note this scaling has no effect on the actual shapeof the velocity/density profiles, just a vertical offset in thedensities.Varying Σ g,1AU has no other effect than to scale ρ g,0 ,and hence ξ . The exponents p and q have a similar effect.The blue shaded region and axes in Figure 5 illustrate theentire range of midplane densities produced by p ∈ [0 . , . q ∈ [0 . , .
8] at R = 5 AU and their associated valuesof s max assuming ξ = 10 − . The range of ( p, q ) pairs givesthe data a slight spread, but the average slope of ρ g,0 vs s max is almost identical to that obtained for ξ . This is nottrue for H g,1AU and H , however. In fact, H has the oppositeeffect because it is inversely proportional to ρ g,0 . Smearingthe same density over a larger volume results in a diminishedbase flow density and a smaller maximum grain size s max . Massive stars tend to have more massive discs and higherluminosities (Andrews et al. 2013; Andrews 2015), therebyaffecting Σ g and ξ , respectively. However, before addressingthese complexities, it is instructive to look at how gravityalone affects dust entrainment in winds. From Equation (3)we can see that the overall velocity profile for the gas de-creases as stellar mass increases. This is illustrated by thedashed lines in the top panel of Figure 6 using the followingstellar masses, M = [0 . , . , , . , . M (cid:12) . The result-ing dust velocities for s = 0 . µ m (solid lines in top panel)are shown in each case and go from being well-entrained at0 . M (cid:12) to being completely non-entrained at 1 . M (cid:12) . Fur-thermore, in Figure 7 we find that s max globally decreasesas the stellar mass increases. Thus, dust entrainment is in-versely proportional to the stellar mass. z [AU] v [ k m = s ] M = 0 : M - M = 0 : M - M = 1 M - M = 1 : M - M = 1 : M - z [AU]
0 5 10 15 20 25 l og ( ; ) [ g = c m ] -20-19-18-17 Figure 6.
Similar to Figure 2, but varying stellar mass instead ofgrain size. Gas and dust properties are represented using dashedand solid lines, respectively. The grain size for all dust curves is s = 0 . µ m. Increasing the stellar mass reduces the initial outflowvelocities of both gas and dust. As a result, fewer grain sizes canbe entrained in the flow. R [AU]
1 10 100 s m a x [ m ] : M - : M - M - : M - : M - Figure 7.
The maximum entrainable grain size plotted as a func-tion of radius for the stellar masses listed in the plot. The largergravitational well produced by higher mass stars require more en-ergy to escape, thus prohibiting the larger, slower moving grainsfrom becoming entrained in the outflow.MNRAS000
The maximum entrainable grain size plotted as a func-tion of radius for the stellar masses listed in the plot. The largergravitational well produced by higher mass stars require more en-ergy to escape, thus prohibiting the larger, slower moving grainsfrom becoming entrained in the outflow.MNRAS000 , 1–11 (2016)
Hutchison, Laibe, & Maddison
The gas and dust densities in the bottom panel of Fig-ure 6 are phenomenologically similar to those in Figure 2,but illustrate two new points of interest. First, the densitygradient in the wind gets steeper as the mass of the star in-creases. This suggests that the gas/dust envelope surround-ing photoevaporating discs is smaller for host stars thatare more massive. Second, placing identical dust distribu-tions around stars of different masses will produce differentdensity profiles in the wind. The weakened entrainment athigher stellar masses will result in diminished outflow ve-locities and higher wind densities. Therefore, although weexpect a narrower grain size distribution in the wind, theirdensities will be enhanced compared to the same grains atlower stellar mass.We can now see that increasing the stellar mass weak-ens dust entrainment while increasing disc mass and/or stel-lar luminosity strengthens it. It is possible to tune thesethree parameters to minimise their combined effect on en-trainment, but there is always a residule density signature.For example, compensating for the weakened entrainmentaround a high mass star by increasing ρ g,0 and/or ξ would,according to Figure 5, increase the already inflated winddensity. Thus, even in the case where dust entrainment isminimally affected, higher mass stars should exhibit higherdust densities in their outflows. In this section we use a turbulent disc model to derive thedust density for a settled disc. With this density relation,we derive a new constraint on the maximum grain size inwinds that reflects how some entrainable dust grains settlewell below the launch point for the flow. Finally, we discusshow settling affects the results in Section 3.
We saw in Section 2.3, that settling does not affect the accu-racy of the semi-analytic model in determining s max or thevelocity of entrained dust grains – as long as there is a non-zero dust density at the ionisation front. With the velocitiespinned down, the continuity equation ensures that the shapeof the density profile is known too. However, Figure 8 showsthat dust settling within the disc will result in a size depen-dent ρ d,i that requires external calibration. In principle, thisis similar to the 40% shift we applied earlier to ρ g,i , but ismade complicated by the unique interaction each grain sizehas with the disc. Calibrating ρ d,i from simulations wouldrender the model superfluous, so we approximate the densityusing a turbulent disc model.The continuity approximation used throughout Sec-tion 3 can already be thought of as extreme turbulence inthe disc. Here we just switch to a more physical descriptionby assuming a finite turbulent viscosity, ν t = αc / Ω K where α is a dimensionless constant (Shakura & Sunyaev 1973).Following Dubrulle et al. (1995), the vertical distribution ofdust density in a turbulent disc governed by ∂ρ d ∂t + ∇ ( ρ d ∆ v ) = ∇ (cid:20) ρ g κ t ∇ (cid:18) ρ d ρ g (cid:19)(cid:21) , (25) -10 -5 0 5 10-30-20 z [AU] l og ( ρ ) [ g / c m ] s = 0.01 µ m s = 0.05 µ m s = 0.1 µ m s = 0.2 µ m s = 0.3 µ m -10 Figure 8.
Dust densities for the entrained dust grains in thehydrodynamic simulations run in Section 2.3. The semi-analyticmodel can predicted the shape of the density profile because ρ d ∝ /v d , but it cannot account for the initial density offset for eachgrain size in a settled disc without first knowing the density ofeach grain size at the ionisation front. with the turbulent diffusivity, κ t , is given by κ t = Ω K αH √ , (26)where Γ = 5 / K t s (cid:28)
1, which is true forall grain sizes considered in this study. Noting that (i) thevertical timescale is much shorter than the radial timescale;(ii) assuming that the dust is always small enough to set-tle at the local terminal velocity (i.e. ∆ v z = − z Ω t s ); and(iii) restricting ourselves to stationary solutions allows us torewrite Equation (25) as a separable first-order differentialequation, − z Ω t s ρ d = κ t (cid:18) zH ρ d + d ρ d d z (cid:19) , (27)which has the solution ρ d = ε ρ g exp (cid:20) − √ α ρ grain s Ω K ρ g,0 c s (cid:18) ρ g,0 ρ g − (cid:19)(cid:21) . (28) Our hydrodynamic simulations are non-turbulent, but theycan provide an order of magnitude estimate of the α neededin this model to produce a realistic stratified dusty disc. Us-ing α = 0 .
05, Equation (28) reproduces a density structurethat is similar to that shown in Figure 8. Although the dustdensity in the model and simulation can drop to arbitrarilylow values, it is natural to believe that below some thresholdthe density becomes physically irrelevant. Alternatively, thethreshold could be caused by a technological limitation inthe observing power of a particular telescope. Either way,limiting ρ d allows us to invert Equation (28) and solve for MNRAS , 1–11 (2016) n the maximum entrainable grain size Table 1.
Coefficients C mn for the fit of (cid:15) i in Equation (31). m n − .
046 00 0 .
852 90 − .
021 26 − .
001 0921 − .
808 50 0 .
089 92 0 .
004 97 02 − .
068 33 − .
004 08 0 0 the maximum grain size, s tur , allowed/observed to enter thewind due to settling in the disc: s tur = α (cid:112) π (1 + Γ) (cid:18) ξ − ξ (cid:19) (cid:18) Σ g ρ grain (cid:19) log (cid:18) ε ε i (cid:19) , (29)where we have defined ε i ≡ ρ d ( z i ) /ρ g ( z i ).In this study, we are more interested in the fundamentalprocess by which settling limits the grains in the flow thanby observational feasibility; therefore, we define the abovethreshold as the point where 100 × (1 − β )% of the grainsare contained within | z | < z i . For this to be meaningful,we assume 0 < β (cid:28)
1. The role of β in this analysis is todetermine the likelihood that a given grain size will havesettled below the ionisation front and not be entrained inthe wind. Mathematically, we do this in three steps: (i) wenormalise the density to get the probability density function(PDF), (ii) we integrate the PDF over the interval [ − z i , z i ]and set it equal to the desired threshold fraction of grainsin the disc: (cid:82) z i − z i ρ d ( z, s ) d z (cid:82) ∞−∞ ρ d ( z, s ) d z = 1 − β, (30)and (iii) we solve for the critical value s that makes thisrelation true. The family of solutions for Equation (30) isindeed given by Equation (29), but ε i is a non-linear functionof both ξ and β . In the interest of keeping s tur analytic, wefit ε i with a polynomial surface ε i (cid:39) (cid:88) m =0 3 (cid:88) n =0 C mn ξ m β n , (31)with coefficients C mn given in Table 1. The root-mean-square error for this fit is ∼ s tur in this model, the more likely it will applyin other turbulent disc models, albeit as a softer limit. Nu-merically, the smallest value of β we can use before experi-encing roundoff errors is 10 − . This corresponds to a (cid:38) σ confidence level that s tur has settled below the ionisationfront. The analytic form of s tur allows us to decrease β muchlower, revealing that s tur asymptotically reaches a maximumvalue. However, extrapolating far beyond the fitted data of-ten leads to erroneous results, so we feel it is sufficient toassume β = 10 − .Although we believe these predictions are representa-tive of discs undergoing photoevaporation, the exact valuesand likelihoods are biased by the assumptions of our tur-bulent disc model. Turbulence in the upper disc is not wellunderstood and could be different to what we assume here.Furthermore, recent non-ideal magneto-hydrodynamic sim-ulations suggest that the accretion stress may be largelylaminar in the inner (cid:46)
30 AU of the disc (e.g. Bai 2013, 2014; Lesur et al. 2014; Gressel et al. 2015; Simon et al.2015), thereby reducing values of α controlling turbulentmixing in these regions. This would result in even smallervalues of s tur , i.e. a smaller reservoir of grains to be car-ried out by the wind.. Finally, our analysis does not takeinto account the collisional evolution and redistribution ofsmall dust grains, which can potentially trap small grainsin a layer significantly thinner than the gas (Krijt & Ciesla2016). This could similarly reduce s tur below that impliedby Equation (29). Given this new constraint, the true maximum grain size inthe wind is given by the minimum of s max and s tur . Fig-ure 9 shows the effect this new relation has on our earliercalculations. The left panel varies the disc’s surface densityexponent p while the right panel varies stellar mass (note s tur does not depend on q ). The dashed and dotted linesshow s max and s tur respectively, while the solid lines indi-cate the true maximum grain size entrained at each radius.In the right panel we have assumed that both ξ and Σ g,1AU are multiplied by a mass factor, (cid:112) M/M (cid:12) . The mass factorhere is simply illustrative; without a mass dependence, thecurves would be colinear. As mentioned earlier, any posi-tive correlation like this between stellar mass and ξ and/orΣ g,1AU will increase the dust density in the wind. Thus, de-spite both s max and s tur being inversely proportional to stel-lar mass, when taking into account that telescopes have alimited resolution, the higher wind densities at higher stellarmasses may result in a larger range of observable grain sizesin the wind (as opposed to entrainable grain sizes).One of the most striking features of Figure 9 is thatdust settling is the dominant mechanism determining themaximum grain size in the wind at almost all radii inthe disc. Moreover, since photoevaporation only operates at R > R c,EUV , the region dominated by s max can be verysmall. This is an important result because it means thatonly dust in the inner few AU of the disc will experienceweak entrainment in the winds. Consequently, processes re-liant on slow moving dust grains, such as mass transport tothe outer disc and pile-ups at the disc surface, are restrictedto the inner disc. Note, this does not rule out photoevapo-ration as the transport mechanism for the earlier mentionedcrystalline grains since these grains are formed in situ in theinner disc.The values used for ξ and Σ g,1AU are fairly conservativeand increasing either one will lead to larger grains entrainedin the outflow. However, since both s max and s tur are affectedsimilarly, the transition from wind-limited to settling-limitedflow remains largely unaffected. The slight change that doesoccur is because the dust’s density profile is no longer aperfect gaussian like that of the gas. Ultimately, this distortsthe linear density-entrainment relationship discussed earlier,but the deviation is small and decreases as ξ increases. Strictly speaking, the results above are only applicable toEUV photoevaporation, but the flows produced by EUV,
MNRAS000
MNRAS000 , 1–11 (2016) Hutchison, Laibe, & Maddison
1 10 100 R [AU] s m a x [ µ m ]
1 10 100 R [AU] . M (cid:31) . M (cid:31) M (cid:31) . M (cid:31) . M (cid:31) p = 0 p = 0 . p = 0 . p = 1 . p = 1 . s max s tur s max s tur Figure 9.
The maximum grain size in the wind as a function of radius while varying p (left) and M (right) assuming a minimumdust-to-gas ratio, ε i = 10 − . The dotted and dot-dashed lines show the maximum entrainable grain size with and without settling,respectively. The true value of s max is the minimum of these two curves. In the right panel we assume that both ξ and Σ g,1AU have a (cid:112) M/M (cid:12) dependence. The stark cutoff in the inner disc is due to quenching of photoevaporation inside of R c,EUV in Equation (24). Inevery case, settling severely limits the value of s max in the wind. X-ray, and FUV radiation are all hydrodynamically simi-lar. Therefore, we expect dust grains entrained by X-rayor FUV driven winds to behave qualitatively like the dustgrains described in this paper. However, there are some im-portant differences to consider that can influence the waydust grains behave in winds (see Alexander et al. 2014; Gortiet al. 2016). EUV radiation accelerates winds faster thanthe other energy regimes because it almost instantaneouslyheats winds to ∼ K. Heating due to X-rays and FUV ra-diation is more gradual and cooler. Maximum X-ray temper-atures reach ∼ > (cid:46)
100 K in the outerdisc. In contrast, the energy regimes rank in opposite or-der in terms of penetration as evidenced by their associatedpenetration columns of neutral hydrogen: N H,EUV ∼ –10 cm − , N H,X-ray ∼ –10 cm − , and N H,FUV ∼ –10 cm − (Ercolano et al. 2009; Alexander et al. 2014)).While the lower outflow velocities (due to lower tem-peratures) achieved by X-ray and FUV radiation reduce en-trainment efficiency by at least a factor of a few (HPLM16),according to Section 3.3 the significantly larger penetrationdepths likely overcompensate and produce better entrain-ment overall, raising the value of s max . At the same time,by penetrating deeper into the disc, X-ray and FUV windshave access to larger, more settled grains than EUV winds.As a result, s tur also increases and we may expect the samequalitative behaviour as shown in Figure 9, but shifted tolarger grain sizes.The indirect thermal coupling with the dust should beminimal for X-rays since they are almost unaffected by smalldust grains in the upper atmospheres of discs (see Gorti et al.2016). However, X-rays are primarily absorbed by heavy el-ements in the disc – many of which can be found in dustgrains – and the effects of dust settling on X-ray driven pho-toevaporation have not been studied in detail (Alexanderet al. 2014). FUV winds are much more likely to be affectedby changes in the dust phase. Gorti et al. (2015) report no-ticeable reductions in disc lifetimes when dust evolution (i.e. settling, migration, and coagulation/fragmentation) is con-sidered alongside FUV photoevaporation. While increasedmass-loss rates imply better entrainment, the radial varia-tion in FUV wind temperatures and fragmentation efficiency(i.e. replenishment of small grains), plus the tendency fordust to migrate inwards, suggests the increase does not occuruniformly across all radii. Furthermore, the opacity from en-trained dust grains has never been modeled self-consistentlyand could impede FUV photoevaporation rates, especiallyin the outer disc. Further studies need to be conducted inorder to determine the relative importance of each of theseprocesses and how they affect the radial profile of grain sizescarried into the winds. We have developed a simple but powerful model using a non-rotating, plane-parallel, photoevaporating atmosphere to es-timate the conditions by which dust grains can be carriedinto a photoevaporative wind. Equation (9) gives the maxi-mum grain size, s max , that can be entrained by the outflow.The model accurately recovers almost all of the results fromour more rigorous hydrodynamic simulations in HPLM16 fordifferent stellar and disc parameters. This implies that anyobservational constraint on the grain size in the wind can betranslated into a strict lower bound on the mass loss rate ofthe disc, Equation (23).This relation is consistent with, but more versatile thanearlier estimates of s max by Takeuchi et al. (2005) and Owenet al. (2011). We show in particular that s max varies withdisc radius with max( s max ) < µ m for typical T Tauristars. However, the largest grain size entrained in the flowmay be much smaller than s max , since dust settling pre-vents the replenishment of large grains in the wind (ex-cept in the inner few AU of the disc). In addition to de-termining the maximum entrained grain size, we uncoverfive distinct behavioural classes of dust grains in photoe- MNRAS , 1–11 (2016) n the maximum entrainable grain size vaporative winds: perfectly-, well-, weakly-, partially-, andnon-entrained grains. These five classes exhibit different out-flow velocities, which may lead to stratified dust layers in thewind and recapture of weakly-entrained dust grains by thedisc at large radii.The stellar mass has a non-linear relationship with s max that alters the shape of the velocity and density profiles forthe gas and dust. We find that more massive stars will tendto host winds with a more compact envelope, higher winddensity, and higher dust-to-gas ratio. Although the maxi-mum entrainable grain size decreases with increasing stellarmass, any positive correlation between luminosity and/ordisc mass with the mass of the central star may result in alarger range of observable grain sizes at higher stellar mass. ACKNOWLEDGEMENTS
We thank Daniel Price for useful discussions and the anony-mous referee whose comments helped improve this paper.M.H. acknowledges funding from a Swinburne UniversityPostgraduate Research Award (SUPRA) and support fromthe Astronomical Society of Australia. G.L. acknowledgesfunding from the European Research Council for the FP7ERC advanced grant project ECOGAL. S.T.M. acknowl-edges partial support from PALSE (Programme Avenir LyonSaint-Etienne). Simulations were performed on the swin-STAR supercomputer at Swinburne University of Technol-ogy.
REFERENCES
Adams F. C., Hollenbach D., Laughlin G., Gorti U., 2004, ApJ,611, 360Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L.,2014, Protostars and Planets VI, pp 475–496Andrews S. M., 2015, PASP, 127, 961Andrews S. M., Williams J. P., 2005, ApJ, 631, 1134Andrews S. M., Williams J. P., 2007, ApJ, 671, 1800Andrews S. M., Rosenfeld K. A., Kraus A. L., Wilner D. J., 2013,ApJ, 771, 129 Bai X.-N., 2013, ApJ, 772, 96Bai X.-N., 2014, ApJ, 791, 137Calvet N., Patino A., Magris G. C., D’Alessio P., 1991, ApJ, 380,617Chiang E. I., Goldreich P., 1997, ApJ, 490, 368Clarke C. J., Alexander R. D., 2016, MNRAS,Corless R., Gonnet G., Hare D., Jeffrey D., Knuth D., 1996, Ad-vances in Computational Mathematics, 5, 329D’Alessio P., Cant¨o J., Calvet N., Lizano S., 1998, ApJ, 500, 411Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237Dutrey A., Guilloteau S., Duvert G., Prato L., Simon M., SchusterK., Menard F., 1996, A&A, 309, 493Epstein P. S., 1924, Physical Review, 23, 710Ercolano B., Clarke C. J., Drake J. J., 2009, ApJ, 699, 1639Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593Gorti U., Hollenbach D., Dullemond C. P., 2015, ApJ, 804, 29Gorti U., Liseau R., S´andor Z., Clarke C., 2016, Space Sci. Rev.,Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ,801, 84Hansen B. M. S., 2014, MNRAS, 440, 3545Hollenbach D., Johnstone D., Lizano S., Shu F., 1994, ApJ, 428,654Hutchison M. A., Laibe G., 2016, Publ. Astron. Soc. Australia,33, e014Hutchison M. A., Price D. J., Laibe G., Maddison S. T., 2016,Monthly Notices of the Royal Astronomical Society, 461, 742Krijt S., Ciesla F. J., 2016, ApJ, 822, 111Laibe G., Gonzalez J.-F., Maddison S. T., 2012, A&A, 537, A61Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56Liffman K., 2003, Publ. Astron. Soc. Australia, 20, 337Miyake T., Suzuki T. K., Inutsuka S.-i., 2016, ApJ, 821, 3Owen J. E., Ercolano B., Clarke C. J., 2011, MNRAS, 411, 1104Pinte C., Laibe G., 2014, A&A, 565, A129Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Simon J. B., Lesur G., Kunz M. W., Armitage P. J., 2015, MN-RAS, 454, 1117Takeuchi T., Clarke C. J., Lin D. N. C., 2005, ApJ, 627, 286Testi L., et al., 2014, Protostars and Planets VI, pp 339–361Veberiˇc D., 2012, Computer Physics Communications, 183, 2622Watson A. M., Stapelfeldt K. R., Wood K., M´enard F., 2007,Protostars and Planets V, pp 523–538Woitke P., et al., 2016, A&A, 586, A103This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000