TY - JOUR
AU - Clarke, Cathie, J
AB - ABSTRACT We model the gas and dust dynamics in a turbulent protoplanetary disc undergoing extreme-UV photoevaporation in order to better characterize the dust properties in thermal winds (e.g. size distribution, flux rate, trajectories). Our semi-analytic approach allows us to rapidly calculate these dust properties without resorting to expensive hydrodynamic simulations. We find that photoevaporation creates a vertical gas flow within the disc that assists turbulence in supplying dust to the ionization front. We examine both the delivery of dust to the ionization front and its subsequent entrainment in the overlying wind. We derive a simple analytic criterion for the maximum grain size that can be entrained and show that this is in good agreement with the results of previous simulations where photoevaporation is driven by a range of radiation types. We show that, in contrast to the case for magnetically driven winds, we do not expect large-scale dust transport within the disc to be effected by photoevaporation. We also show that the maximum size of grains that can be entrained in the wind (smax) is around an order of magnitude larger than the maximum size of grains that can be delivered to the front by advection alone (⁠|$s_{\mathrm{crit}}\lesssim 1 \,\, \mu {\mathrm{m}}$| for Herbig Ae/Be stars and |$\lesssim 0.01 \,\, \mu {\mathrm{m}}$| for T Tauri stars). We further investigate how larger grains, up to a limiting size slimit, can be delivered to the front by turbulent diffusion alone. In all cases, we find smax ≳ slimit so that we expect that any dust that is delivered to the front can be entrained in the wind and that most entrained dust follows trajectories close to that of the gas. protoplanetary discs, circumstellar matter, planetary systems, stars: pre-main sequence, dust, extinction 1 INTRODUCTION Planets are born in the nurturing cocoon of gas and dust created by single, or potentially binary (Quintana et al. 2007), parent stars early in their own development. While the circumstances surrounding the birth of each planet may be unique – from fragmentation via gravitational instability (e.g. Kuiper 1951; Boss 1997; Inutsuka, Machida & Matsumoto 2010; Nayakshin 2010) to streaming instability of large dust reservoirs (e.g. Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007) followed by core accretion (e.g. Safronov 1969; Pollack et al. 1996; Ikoma, Nakazawa & Emori 2000) – all planets, regardless of how they were born, inevitably experience the death of the protoplanetary disc phase early in their evolution (≲ 5–10 Myr; e.g. Haisch, Lada & Lada 2001; Hernández et al. 2007; Mamajek 2009). These discs are crucial to the early development of planets because they provide the fodder for growth and a driver of radial migration (e.g. Goldreich & Tremaine 1980; Ward 1986, 1997; Tanaka, Takeuchi & Ward 2002). Therefore, putting constraints on the dispersal mechanisms and/or dispersal rates of gas and dust in these discs is important for understanding the architectures of planetary systems we observe today (e.g Alexander & Pascucci 2012; Mordasini et al. 2012). Many dispersal mechanisms are likely to play a role in the evolution of protoplanetary discs: viscous accretion (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), planet–disc interactions (Calvet et al. 2005), grain growth (Dullemond & Dominik 2005), photophoresis (Krauss et al. 2007), MRI-driven winds (Suzuki & Inutsuka 2009), binary star interactions (Marsh & Mahoney 1992), MRI-driven dust depleted flows (Chiang & Murray-Clay 2007), magneto-thermal winds (Bai et al. 2016), and photoevaporation (Clarke, Gendrin & Sotomayor 2001). It is however often assumed that photoevaporation is responsible for the final clearing of the disc, in part because of its success in explaining the rapid transition between disc-bearing and disc-less stars (Owen, Ercolano & Clarke 2011b; Ercolano et al. 2015). Predictions of gas mass-loss rates (e.g. Ercolano et al. 2008; Gorti, Dullemond & Hollenbach 2009; Owen, Clarke & Ercolano 2012; Picogna et al. 2019) and their effect on disc evolution (Ercolano & Rosotti 2015; Ercolano & Pascucci 2017; Carrera et al. 2017; Jennings, Ercolano & Rosotti 2018) have been studied in great detail. While dust was initially assumed to be blown out with the gas (Johnstone, Hollenbach & Bally 1998; Matsuyama, Johnstone & Hartmann 2003; Adams et al. 2004), Throop & Bally (2005) provided the first study of the differential removal of dust and gas in a photoevaporative wind. With the majority of dust mass in discs locked up in large grains near the disc mid-plane and traced by its sub-mm thermal emission, the best prospect of using dust to trace photoevaporative winds is through scattered light imaging above the disc photosphere. Owen, Ercolano & Clarke (2011a) found that dust pulled from the disc can leave a characteristic imprint on the surrounding environment that can be detected in older edge-on discs. They further found that the morphology of their scattered light images depended on the maximum grain size pulled into the wind at different radii. Takeuchi, Clarke & Lin (2005) provided the first order-of-magnitude estimate for the maximum grain size in winds. Later, Hutchison, Laibe & Maddison (2016b, hereafter HLM16b) found an analytic relation for the maximum entrainable grain size in vertical winds and postulated that the equation could be inverted to give a lower limit on the gas mass-loss rate if observations could constrain the grain sizes that were being pulled from the disc. While these examples are illustrative of how dust may be used to constrain disc winds, each of these studies suffers from some limitation in modelling the delivery and/or subsequent entrainment of solids in the wind. To date, there have been five two-phase (gas + dust) photoevaporation wind models presented in the literature – none of which provide the robust understanding that we need to confidently set constraints on dispersal mechanisms using dust measurements: Owen et al. (2011a) used a full grain-size distribution for the dust (Mathis, Rumpl & Nordsieck 1977), a self-consistent ionization front prescription for the gas (Hollenbach et al. 1994), and checked for dust entrainment along gas streamlines in a 2D extreme-ultraviolet (EUV)-driven wind. The disadvantage with their model is that the flow originates from an infinitely-thin, perfectly-mixed disc and the dust in the wind was required to have the same streamlines as the gas. Facchini, Clarke & Bisbas (2016) also used a realistic grain-size distribution and, with the aid of the 3D-PDR code (Bisbas et al. 2012), self-consistently solved the hydrodynamics of gas and dust in a radial flow from the outer disc. In addition to only being 1D, their model is restricted to external photoevaporation caused by far-UV (FUV) sources. Hutchison et al. (2016a, hereafter HPLM16a) used smoothed particle hydrodynamics to self-consistently model the two-phase hydrodynamics with a single grain size in thin vertical atmospheres undergoing EUV photoevaporation. However, their 1D model neglects rotational effects in the wind, the ionization front location is a free parameter, and the disc is laminar (i.e. settling dust never reaches a steady state). HLM16b developed a fast semi-analytic model based on HPLM16a that included turbulent mixing within the disc but treated the gas as being static below the ionization front with no mechanism for advecting the dust upwards. Their model also suffers from having a 1D flow and lacks a self-consistently set ionization front. Franz et al. (2020) used a 2D hydrodynamical disc model irradiated by X-ray and EUV spectra (XEUV) from a central T Tauri star to measure entrainment in winds using tracer particles of various sizes inserted at the disc surface. Since this study only tracks the trajectories of dust particles within the wind, it cannot assess dust mass-loss rates, since the delivery of dust into the wind, and how this depends on the grain size and the properties of the underlying turbulent disc, is not addressed. One of the common themes emerging from the work of HPLM16a and HLM16b is that the local conditions in the disc beneath the flow can limit the dust distribution in the wind, despite the wind being able to entrain larger grains. Therefore, accurate diagnostics of dust in winds (such as size distribution, mass-loss, and density) require models that can self-consistently link the gas and dust properties back to their respective reservoirs in the disc. Two-phase, 3D radiation hydrodynamic simulations can potentially provide the best constraints on dusty winds, but as such models are not readily available, we focus on combining and improving upon the following elements from previous studies: We set the ionization front location as a function of incident EUV flux (Hollenbach et al. 1994). Using an adaptation of the self-similar 2D wind solution of Clarke & Alexander (2016) that accounts for the finite height of the disc surface, we set the gas velocity above the ionization front so as to enable the 2D ionized wind solution to pass smoothly through its critical point. Below the ionization front, we solve for the vertical profile of dust species of various sizes including gravitational settling, turbulent diffusion, and hydrodynamic drag (from the vertical gas flow feeding the base of the wind). We approximate and spatially resolve the physical width of the ionization front to provide a smooth transition between the turbulent disc and the laminar wind. The combination of the above features allows us to track the flow of dust from the disc mid-plane into the wind, thereby coupling the dusty wind to the disc. Our goal is to investigate the size range, flux, and subsequent trajectory of dust passing through the ionization front. In what follows, we will refer to processes by which grains arrive at the ionization front as delivery and the subsequent dynamical evolution of the dust in the ionized wind flow as entrainment. The structure of this paper is as follows. In Section 2, we describe the vertical structure of the gas disc and how we solve for the steady-state dust distribution and corresponding dust flux, as a function of grain size. In Section 3, we describe results relating to both dust delivery and entrainment. Section 4 discusses the implications of our results for dusty photoevaporating discs, and Section 5 summarizes our conclusions. 2 METHODS We proceed by building a gas model for the disc and wind. Then we add dust to the disc mid-plane and solve the fluid equations for the dust assuming the backreaction on the gas is negligible (for justification, see HPLM16a). In what follows, we use the subscripts g and d to distinguish between gas and dust variables with the same name. We consider a disc whose internal structure is approximately described by the following power-law parametrizations (see, e.g. Laibe, Gonzalez & Maddison 2012): $$\begin{eqnarray*} \rho _{\mathrm{g,0}}& = \frac{\Sigma _{\mathrm{g}}}{\sqrt{2\pi } H}, \end{eqnarray*}$$(1) $$\begin{eqnarray*} \Sigma _{\mathrm{g}}& \propto \left(\frac{R}{R_{\mathcal {G}}} \right)^{-p}, \end{eqnarray*}$$(2) $$\begin{eqnarray*} H& \propto \left(\frac{R}{R_{\mathcal {G}}} \right)^{3/2- q/2}, \end{eqnarray*}$$(3) $$\begin{eqnarray*} c_{\mathrm{s}}& = \Omega _{\mathrm{K}}H\propto \left(\frac{R}{R_{\mathcal {G}}} \right)^{-q/2}, \end{eqnarray*}$$(4) where ρg, 0 is the mid-plane density, Σg is the local surface density, |$\Omega _{\mathrm{K}}= (\mathcal {G}M/R^3)^{1/2}$| is the Keplerian orbital frequency, and p and q are power-law exponents respectively controlling the surface density and temperature profiles of the disc. In this paper, we assume the following generic reference values for our disc: p = 1, q = 0.5, M = M⊙, Σg, 1 au = 100 g cm−3, H1 au = 0.05 au, and |$c_{\mathrm{s}}{}_{,1\!\,\, {\rm{\small AU}}}= \Omega _{\mathrm{K}}{}_{,1\!\,\, {\rm{\small AU}}}\, H{}_{1\!\,\, {\rm{\small AU}}}\approx 1.5 \,\, {\mathrm{km}} \,\, {\mathrm{s}}^{-1}$|⁠. Later, however, we will modify the radial power-law profiles for the sound speed and disc scale height in order to be consistent with the height of the ionization front and the vertical flow it sets up within the disc.1 2.1 Gas flow We begin by identifying the physical constraints at the ionization front that are needed to connect the gas flow in the disc to the wind. 2.1.1 Jump conditions at the ionization front Conservation of mass and momentum across the ionization front allows us to relate the density (ρg) and the perpendicular velocity (u⊥) in the neutral disc with corresponding quantities in the ionized wind: $$\begin{eqnarray*} \rho _{\mathrm{g,n}}u_{\mathrm{n}\perp }& = \rho _{\mathrm{g,i}}u_{\mathrm{i}\perp }, \end{eqnarray*}$$(5) $$\begin{eqnarray*} \rho _{\mathrm{g,n}}c_{\mathrm{s,n}}^2 + \rho _{\mathrm{g,n}}u_{\mathrm{n}\perp }^2 & = \rho _{\mathrm{g,i}}c_{\mathrm{s,i}}^2 + \rho _{\mathrm{g,i}}u_{\mathrm{i}\perp }^2, \end{eqnarray*}$$(6) where cs is the sound speed of the gas, and the subscripts n and i refer to the neutral and ionized sides, respectively. Setting the sound speed on the neutral side to be cs, n = cs (which is assumed to be a function of radius only), we can solve equations (5) and (6) for ρg, n and un⊥. After some manipulation, we find $$\begin{eqnarray*} u_{\mathrm{n}\perp }& = \dfrac{c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \pm \sqrt{\left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \right)^2 - 4 c_{\mathrm{s}}^2 u_{\mathrm{i}\perp }^2}}{2 u_{\mathrm{i}\perp }}, \end{eqnarray*}$$(7) $$\begin{eqnarray*} \rho _{\mathrm{g,n}}& = \dfrac{\rho _{\mathrm{g,i}}u_{\mathrm{i}\perp }}{u_{\mathrm{n}\perp }} = \frac{\rho _{\mathrm{g,i}}\left[ c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \mp \sqrt{\left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \right)^2 - 4 c_{\mathrm{s}}^2 u_{\mathrm{i}\perp }^2} \right] }{2 c_{\mathrm{s}}^2}, \end{eqnarray*}$$(8) where we take the negative sign in equation (7) and the positive sign in equation (8) to be the physical solutions because the velocity in the disc must be smaller than the wind speed. If the gas flow is perpendicular to the ionization front (i.e. if ui = ui⊥ and un = un⊥), then equations (7) and (8) are sufficient to constrain the physical properties of our flow. However, when the ionization front is tilted by angle iIF with respect to the flow, as assumed in this study, there is an additional constraint on the gas velocity; namely that the velocity parallel to the front remains constant across the discontinuity: ui∥ = un∥ = u∥. Assuming that we know the wind density and total velocity at the base of the ionized flow as well as the sound speed in the disc and wind, which leaves us with five unknown variables. The neutral quantities un⊥ and ρg, n continue to depend on ui⊥ according to equations (7) and (8), while the remaining unknown velocities (ui⊥, un, and u∥) can all be related using trigonometric relations, thereby closing the system of equations. With some effort, it can be shown that u∥ satisfies the following polynomial equation: $$\begin{eqnarray*} c_6 u_\parallel ^6 + c_4 u_\parallel ^4 + c_2 u_\parallel ^2 + c_0 = 0, \end{eqnarray*}$$(9) with coefficients $$\begin{eqnarray*} c_6 & = \cot ^2{i_{\rm IF}} \csc ^2{i_{\rm IF}}, \end{eqnarray*}$$(10) $$\begin{eqnarray*} c_4 & = \left[ 2 \left(c_{\mathrm{s}}^2 - c_{\mathrm{s,i}}^2 \right) - u_{\mathrm{i}}^2 \left(1 + \csc ^2{i_{\rm IF}} \right) \right] \cot ^2{i_{\rm IF}}, \end{eqnarray*}$$(11) $$\begin{eqnarray*} c_2 & = c_{\mathrm{s}}^4 + \left[ \left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}}^2 \right)^2 - 2 c_{\mathrm{s}}^2 u_{\mathrm{i}}^2 \right] \cot ^2{i_{\rm IF}}, \end{eqnarray*}$$(12) $$\begin{eqnarray*} c_0 & = - c_{\mathrm{s}}^4 u_{\mathrm{i}}^2. \end{eqnarray*}$$(13) Although the roots of equation (9) can be expressed analytically in closed form, they are too extensive to be of practical use. Instead, we solve for u∥ numerically, after which ui⊥ and the remaining quantities are easily calculated. For a full derivation of equations (7)–(13), see Appendix A. 2.1.2 Determination of the velocity at the base of the ionized flow In order to determine ui, we follow the approach of Clarke & Alexander (2016), who computed the structure of thermally driven axisymmetric winds under the assumption of self-similarity. Clarke & Alexander demonstrated that there is a maximum launch velocity for which the flow structure contains a critical point at the sonic point and that this maximum velocity agrees well with the launch velocity found in hydrodynamical solutions of disc winds. While Clarke & Alexander only treated the case that the wind is launched from the disc mid-plane, we here consider the more realistic scenario of a wind launched from an inclined surface. In order to preserve self-similarity in the wind, we assume this surface makes a constant angle iIF = 0.3 rad to the horizontal, as detailed in Section 2.1.4. In this study, we keep iIF fixed for all parameter combinations in order to make comparisons between simulations easier. We find the launch velocity by sampling |$u_{\mathrm{i}}\in [0,\, c_{\mathrm{s}}]$| and calculating the parallel and perpendicular components ui∥ and ui⊥ according to Section 2.1.1. These components determine the initial launch angle of the flow. Using these initial conditions, we then integrate the equations describing the topology, density, and velocity structure of a self-similar flow.2 Finally, we iterate on this procedure until we find the maximum velocity ui that permits a transonic flow that extends to infinity. 2.1.3 Ionized boundary The ionized density at the base of the wind is independently determined by solving the ionization-recombination balance for an isothermal disc irradiated by EUV radiation. Using the weak stellar wind model from Hollenbach et al. (1994), the density of the ionized wind in the inner few au of the disc can be expressed analytically by $$\begin{eqnarray*} \rho _{\mathrm{g,i}}& = m_{\mathrm{H}}C_0 \left(\dfrac{3 \Phi _{\mathrm{i}}}{4 \pi \alpha _2 R^3} \right)^{1/2}, \end{eqnarray*}$$(14) where R is the cylindrical disc radius measured from a central star of mass M, mH is the mass of hydrogen, α2 = 2.6 × 10−13 cm3 s−1 is the recombination coefficient for all states except to the ground state (case B), Φi is the stellar EUV luminosity, and C0 = 0.2 is the numerical factor used by Hollenbach et al. (1994) to bring their analytical and numerical results into agreement. In this paper, we will only consider two extremes for the EUV luminosity, |$\Phi _{\mathrm{i}}= [10^{41},\, 10^{45}] \,\, {\mathrm{photons}} \,\, {\mathrm{s}}^{-1}$| (hereafter referred to as Φ41 and Φ45), and assume that the resulting wind has an isothermal temperature of |$T = 10~000\, {\rm K}$| or, equivalently, a sound speed cs, i ∼ 10 km s−1. Note Φ41 corresponds to a typical T Tauri star whereas Φ45 would represent an early-type Herbig Be star (Gorti et al. 2009) with a larger mass. Because changing the stellar mass clouds the interpretation of the results, we choose to keep M constant as we vary Φi. Although equation (14) formally only applies to disc radii between where the disc scale height equals the radius of the star and the so-called gravitational radius |$R_{\mathcal {G}}= \mathcal {G}M/ c_{\mathrm{s}}^2$|⁠, we have adopted it as a general relation for our disc. We caution that the resulting mass-loss rate diverges to large radii as |$\sqrt{R}$| and that, in real systems, the radial range over which this approximation holds is limited by the finite energy input of the star. Nevertheless, our reasons for this choice are two-fold: First, we have found that similarity solutions for the wind only exist for flows where the base density has a radial power-law slope ≲ −2, thereby preventing us from using the other relations in the Hollenbach et al. (1994) model. Secondly, the R−3/2 scaling in equation (14) is consistent with XEUV disc winds (Picogna et al. 2019), which will later aid in comparing our results to Franz et al. (2020). 2.1.4 Ionization front location The location of the ionization front was indirectly set when we chose a constant surface inclination iIF while solving for the initial wind velocity; now we must ensure that the gas properties within the disc are consistent with the location we have chosen. We assume that the flow in the neutral region is perpendicular to the disc mid-plane and thus the momentum equation in this region is $$\begin{eqnarray*} u\frac{{\rm d} u}{{\rm d} z} = - \frac{c_{\mathrm{s}}^2 }{\rho _{\mathrm{g}}} \frac{{\rm d} \rho _{\mathrm{g}}}{{\rm d} z} - \frac{ \mathcal {G}Mz}{(R^2 + z^2)^{3/2}}, \end{eqnarray*}$$(15) whose solution is given by3 $$\begin{eqnarray*} u& = c_{\mathrm{s}}\sqrt{- \, \mathrm{\it W}_0 \! \! \left[ -\exp { \left(-\frac{2 \mathcal {G}M}{c_{\mathrm{s}}^2 \sqrt{R^2 + z^2}} - C_1 \right) } \right]}, \end{eqnarray*}$$(16) $$\begin{eqnarray*} \rho _{\mathrm{g}}& = \frac{\rho _{\mathrm{g,0}}u_0}{u}, \end{eqnarray*}$$(17) where W0 is the Lambert W function. To find C1, we look for solutions that satisfy the boundary conditions given by equations (1), (7), and (8). For an inclined ionization front, mass conservation requires ρg, 0u0 = ρg, nun/cosiIF, thereby fixing the initial flow velocity at the mid-plane, u0. Then, using the transcendental form of equation (16) (evaluated at z = 0), the expression for C1 is $$\begin{eqnarray*} C_1 = \mathcal {M}_0^2 - \ln \mathcal {M}_0^2 - 2 \mathcal {M}_{\mathrm{K}}^2. \end{eqnarray*}$$(18) For convenience here and in what follows, we have defined the Mach numbers |$\mathcal {M}_{\mathrm{K}}\equiv v_{\mathrm{K}}/ c_{\mathrm{s}}$| (where vK = RΩK is the mid-plane Keplerian velocity) and |$\mathcal {M}\equiv u/ c_{\mathrm{s}}$|⁠. Using the latter definition, |$\mathcal {M}_0$| and |$\mathcal {M}_{\mathrm{n}}$| refer to the Mach number at the mid-plane and neutral side of the ionization front, respectively. We can use a similar approach to obtain an analytic relation for the height of the ionization front, zIF. At z = zIF, the transcendental form of equation (16) gives $$\begin{eqnarray*} \mathcal {M}_{\mathrm{n}}^2 - \ln {\mathcal {M}_{\mathrm{n}}^2} = \frac{2 \mathcal {M}_{\mathrm{K}}^2}{\sqrt{1+ \left(z_{\mathrm{IF}}/R\right)^2}} + C_1. \end{eqnarray*}$$(19) Substituting in the expression for C1 and rearranging to find zIF, we obtain $$\begin{eqnarray*} z_{\mathrm{IF}}= R\sqrt{4 \mathcal {M}_{\mathrm{K}}^4 \left[ \mathcal {M}_{\mathrm{n}}^2 - \mathcal {M}_0^2 + \ln {\left(\dfrac{\mathcal {M}_0^2}{\mathcal {M}_{\mathrm{n}}^2} \right)} + 2 \mathcal {M}_{\mathrm{K}}^2 \right]^{-2} - 1}. \end{eqnarray*}$$(20) Inserting the generic disc parameters from equations (1)–(4) gives an opening angle that is not constant with radius. We therefore invert equation (20) to find the required disc sound speed that produces tan −1(zIF/R) = iIF. While the deviation in cs from equation (4) is ≲ 2 at all radii, the fact that cs is no longer independent of the ionization front means that we must include this inversion calculation in our iterative scheme to find ui. Note that changing the sound speed also requires that we alter the disc scale height, which is later used in determining the dust flow. In essence, we sacrifice the simplicity of the radial power-law profiles in equations (3) and (4) to maintain consistency with our inclined ionization front and gas flow. 2.2 Dust flow With the gas solution fully determined up to and immediately above the ionization front, we can now derive the corresponding solutions for the dust. We omit the dynamical effect of dust back-reaction on the gas for the following reasons. In the disc, the gas pressure from the near hydrostatic equilibrium restores any perturbations caused by the dust motion. In the wind, the total dust-to-gas ratio from a realistic grain-size distribution is orders of magnitude smaller than the canonical 0.01, a case already shown by (HPLM16a) to exhibit negligible back-reaction on the gas for individual grain sizes. By solving for dynamical equilibria in the z-direction with fixed mid-plane dust densities, we are implicitly assuming that the vertical flow time (on which this equilibrium is established) is much shorter than the time-scale on which the mid-plane density changes (either as a result of the wind or of radial drift in the disc mid-plane). We find that this condition is readily fulfilled in practice. 2.2.1 Dust dynamics in the turbulent disc To obtain the dust velocities in the disc, we consider the momentum equation,4 $$\begin{eqnarray*} \frac{\partial v}{\partial t} + v\frac{\partial v}{\partial z} \simeq - \frac{ \left(v- u\right)}{t_{\mathrm{s}}} - \frac{\mathcal {G}Mz}{\left(R^2 + z^2 \right)^{3/2}}, \end{eqnarray*}$$(21) where ts is the stopping time in the Epstein drag regime (Epstein 1924), $$\begin{eqnarray*} t_{\mathrm{s}}= \sqrt{\frac{\pi \gamma }{8}} \frac{\rho _{\mathrm{grain}}s }{c_{\mathrm{s}}\rho _{\mathrm{g}}} = \frac{\rho _{\mathrm{eff}}s }{c_{\mathrm{s}}\rho _{\mathrm{g}}}, \end{eqnarray*}$$(22) with ρgrain the intrinsic density of individual dust grains. For convenience, we simplify the expression in the second equality by defining |$\rho _{\mathrm{eff}}\equiv \rho _{\mathrm{grain}}\sqrt{\pi \gamma /8}$| as an effective grain density. Because we are only interested in the steady-state flow, we can omit the time derivative on the left-hand side of equation (21). Furthermore, for grains relevant to photoevaporation, the advection term is typically much smaller than the drag and gravitational components in the disc interior and v in this region can be approximated by the local terminal velocity of the dust: $$\begin{eqnarray*} v= u- \frac{t_{\mathrm{s}}\mathcal {G}Mz}{\left(R^2 + z^2\right)^{3/2}}. \end{eqnarray*}$$(23) With the grain size fixed, the stopping time increases substantially along the flow on account of the exponential decline in gas density with z. This causes the terminal velocity approximation to break down (particularly near the ionization front) and failure to use the correct velocity leads to an overestimate of the dust flux leaving the disc. We emphasize that this is a general phenomenon that affects even the smallest grains due to the step-like transition at the ionization front being comparable in width to the stopping distance of the dust. We therefore obtain v numerically for all dust grains by solving equation (21) using Matlab’s variable-step, variable-order ordinary differential equation solver ODE15s (Shampine & Reichelt 1997; Shampine, Reichelt & Kierzenka 1999). Because v = 0 is a removable singularity of equation (21), large grains whose velocity changes signs during the flow require special attention. In these exceptional cases, we make the substitution |$\operatorname{sign}{(v)}\sqrt{|w|} = v$| in equation (21), which allows us to (i) remove the singularity by absorbing the multiplicative factor of v into the derivative and (ii) maintain a real dust velocity inside of the drag term when the dust velocity is negative. Once we have the dust velocity we proceed with finding the dust density in the disc. There is a rich body of work in the literature investigating the vertical dust distribution in turbulent discs (e.g. Takeuchi & Lin 2002; Schräpler & Henning 2004; Johansen & Klahr 2005; Fromang & Papaloizou 2006; Fromang & Nelson 2009; Charnoz et al. 2011; Birnstiel, Fang & Johansen 2016; HLM16b), the vast majority of which are based on the seminal work of Dubrulle, Morfill & Sterzik (1995). However, as pointed out by Riols & Lesur (2018), the net vertical flow induced by disc winds can alter the dust distribution in the disc. To capture this effect in the Dubrulle et al. (1995) model, we use the dust velocity derived in equation (23) to compute the dust density in the advection-diffusion equation, $$\begin{eqnarray*} \frac{\partial \rho _{\mathrm{d}}}{\partial t} + \frac{\partial }{\partial z} \left[ \rho _{\mathrm{d}}v- \rho _{\mathrm{g}}\mathcal {D}_{\mathrm{d}}\frac{\partial }{\partial z} \left(\frac{\rho _{\mathrm{d}}}{\rho _{\mathrm{g}}}\right) \right] = 0, \end{eqnarray*}$$(24) where we approximate the diffusion coefficient for the dust, |$\mathcal {D}_{\mathrm{d}}$|⁠, using the expression from Charnoz et al. (2011): $$\begin{eqnarray*} \mathcal {D}_{\mathrm{d}}\approx \frac{\alpha c_{\mathrm{s}}H}{\mathrm{Sc}}. \end{eqnarray*}$$(25) Here α is the Shakura–Sunyaev turbulence parameter (Shakura & Sunyaev 1973) and Sc is the dimensionless Schmidt number. There are a number of definitions for Sc in the literature (see Youdin & Lithwick 2007; Laibe 2014, for a discussion), but we have opted to use Sc = 1 + St since we are only considering motion in the vertical direction (Laibe 2014). A shortcoming of this model is that it assumes uniform turbulence in the disc, which may, in reality, vary depending on the physical mechanism generating the turbulence. For example, Shi & Chiang (2014) show that gravito-turbulence is relatively uniform vertically (only differing by a factor of ∼2 from mid-plane to surface), but turbulence induced by the magnetorotational instability is more variable (fluctuating by a factor of ∼15). Another shortcoming is that it assumes that the Stokes number (St ≡ ΩKts) is much smaller than unity. In practice, we find that all grains that enter the wind satisfy this condition in the turbulent region of the disc, but not necessarily in the wind. In Section 3.1, we discuss some limitations to our model that arise when the stopping time becomes comparable with the time-scale on which dust crosses the ionization front. A significant issue in coupling the wind to a turbulent disc is that we need to specify the value of |$\mathcal {D}_{\mathrm{d}}$| in the wind (i.e. above the ionization front). In default of a model for turbulence in the wind, we will assume that this region of the flow is laminar so that |$\mathcal {D}_{\mathrm{d}}$| tends to zero. While it is reasonable to assume that the physical scale on which the disc makes the transition from a neutral to ionized state is governed by recombination physics, the relevant length scale on which |$\mathcal {D}_{\mathrm{d}}$| declines at the ionization front is not obvious. Clearly, step functions for |$\mathcal {D}_{\mathrm{d}}$| and ρg should be avoided because the derivatives in equation (24) would be dependent on the resolution of our numerical grid. For this reason, we smooth the disc quantities |$\mathcal {D}_{\mathrm{d}}$|⁠, ts, cs, and u (smoothing of ρg is indirectly set by u through mass conservation) on to their respective values in the wind5 using a tanh function parametrized in terms of a length scale W. For example, the smoothed functional form for the Stokes number is given by $$\begin{eqnarray*} \mathrm{St}(z) = \frac{1}{2} \left[ (\tilde{\mathrm{St}}(z) + \mathrm{St}_{\mathrm{i}}) - (\tilde{\mathrm{St}}(z) -\mathrm{St}_{\mathrm{i}}) \tanh {\left(\frac{z - z_{\mathrm{IF}}}{W/3} \right)} \right], \end{eqnarray*}$$(26) where |$\tilde{\mathrm{St}}(z)$| is obtained from equation (22) using non-smoothed disc variables and Sti is the value of the Stokes number in the ionized flow above the front. Throughout the discussion of the results we characterize the dust properties in terms of St(zIF); since the Stokes number increases strongly as the density drops at the front, St(zIF) can be seen to be around half of the Stokes number in the fully ionized wind above the ionization front. By applying these smoothed forms for variables at the ionization front, we can test the sensitivity of our results to the assumptions we have made about the relevant length scales while ensuring that W is sufficiently above the grid scale such that spatial derivatives can be reliably computed. When not otherwise specified, we will estimate the width by $$\begin{eqnarray*} W = \frac{m_{\mathrm{H}}}{\sigma \rho _{\mathrm{g}}(z_{\mathrm{IF}})}, \end{eqnarray*}$$(27) where σ = 6.3 × 10−18 cm2 is the photoionization cross-section for neutral hydrogen and ρg(zIF) is the smoothed density at the ionization front. The conceptually simplest approach is to evolve equation (24) on an Eulerian grid for a specific grain size in a background gas velocity and density field that smoothly connects on to the wind. In this implementation, we need only to fix the dust density at the mid-plane; an upper boundary condition is unnecessary as long as the grid extends into the region where the flow becomes laminar and equation (24) effectively becomes first order. In practice, as the solution evolves towards a steady state, it ‘finds’ a steady-state flux, $$\begin{eqnarray*} \mathcal {F}=-{\mathcal {D}_{\mathrm{d}}}\biggl (\frac{\mathrm{d} \rho _{\mathrm{d}}}{\mathrm{d} z} - \frac{\mathrm{d} \ln {\rho _{\mathrm{g}}}}{\mathrm{d} z} \rho _{\mathrm{d}}\biggr) + \rho _{\mathrm{d}}v, \end{eqnarray*}$$(28) and corresponding gradient of the dust-to-gas ratio at the mid-plane with a smoothly varying steady-state dust density profile. With the lack of a second fully prescribed boundary condition, it is not immediately obvious as to why the Eulerian simulation prefers any one solution over another. Fruitfully, we can investigate the family of analytic solutions to equation (24) in a steady-state characterized by constant |$\mathcal {F}$| (or equivalently by choosing the gradient of the dust-to-gas ratio at z = 0). Assuming |$\mathcal {D}_{\mathrm{d}}\ne 0$|⁠, these solutions take the general form $$\begin{eqnarray*} \rho _{\mathrm{d}}= \varepsilon \rho _{\mathrm{g}}e^{I} \left(1 - \int ^z_0 \dfrac{\mathcal {F}e^{-I}}{\varepsilon \rho _{\mathrm{g}}\mathcal {D}_{\mathrm{d}}} \mathrm{d} z^{\prime } \right), \end{eqnarray*}$$(29) where the mid-plane dust-to-gas ratio ε is a function of grain size, and, for convenience, we have defined $$\begin{eqnarray*} I = \int ^z_0 \frac{v}{\mathcal {D}_{\mathrm{d}}} \mathrm{d} z^{\prime }. \end{eqnarray*}$$(30) After experimenting with a range of values of |$\mathcal {F}$| it is readily apparent that |$\mathcal {F}$| only affects the density profile in the region near the ionization front itself, specifically where |$\mathcal {D}_{\mathrm{d}}$| is close to zero. Small values of |$\mathcal {F}$| cause the density to blow up to ∞ while large values cause the density to plummet to −∞. The explanation for this behaviour is as follows. Since |$\mathcal {F}$| effectively sets the total dust flux through the disc, selecting an arbitrarily large flux implies that a correspondingly large amount of dust is being lost to the wind. With the advective velocity fixed, the only way to impose such a condition is to adjust the dust profile so as deliver this flux; in practice, this involves establishing a large negative or positive gradient in the dust-to-gas ratio in that vicinity. These requirements cause the density to become unbounded at the ionization front for all but an extremely narrow range of fluxes that only differ beyond their six or seventh significant digit – for all intents and purposes, a single value of |$\mathcal {F}$|⁠. With this insight, it can be seen that the value of |$\mathcal {F}$| selected is such that the disc and wind solutions connect smoothly at the ionization front (i.e. no maxima, minima, or kinks). Thus, it comes as no surprise that this ‘unique’ value of |$\mathcal {F}$| reproduces the steady-state solution obtained by evolving equation (24) numerically. We therefore conclude that our selection of |$\mathcal {F}$| in the analytic model is both unique and physically motivated. The obvious benefit to using equation (29) over conventional numerical schemes is that, even after sweeping over possible |$\mathcal {F}$| values, we can obtain the steady-state dust density in seconds – a speed-up of ∼104 times over the numerical method detailed below. Furthermore, the analytic solution scales well to very high resolutions, requiring only a few minutes to compute the vertical dust velocity and density at a given disc radius on a grid with 10 million points. This latter quality is particularly important in this study so that we can ensure that the ionization front is always adequately resolved as we lower W (the smallest of which requires 60 million points). As a final remark, it is important to note that equation (29) breaks down as |$\mathcal {D}_{\mathrm{d}}\rightarrow 0$|⁠, which, due to the tanh smoothing, always occurs within ∼3W of the ionization front, regardless of the underlying disc parameters. Likewise, when α ≲ 10−3 equation (29) can also predict unphysical behaviour, even in the bulk of the disc. In these regions where the turbulent diffusion is small, we can (in most cases) circumvent these issues by transitioning to a purely advective solution where the two solutions smoothly connect. In practice, we do this by calculating the derivative of |$\rho _{\mathrm{d}}= \dot{m_{\rm d}}/v$| for each potential transition point, where |$\dot{m_{\rm d}}$| is a constant representing the mass flux at each of these points. We start at the last defined point in the diffusive solution and work backwards toward the disc mid-plane until we find the point where the slope from the two solutions are equal. We have found that the above procedure agrees well with numerical simulations (as long as sufficient resolution is used), thus extending the parameter space that we can model. 2.3 Comparing the semi-analytic model to time-dependent calculations To assess the validity of our disc model, we solved the full time-dependent hydrodynamic equations (21) and (24) using a modified version of the grid code described in appendix A2 of Hutchison, Price & Laibe (2018). Assuming Φi = Φ45 and R = 10 au, we discretize the region |$z \in [0,\, 1.5z_{\mathrm{IF}}]$| with Nz = 40 002 cell-centred grid points, including ghost points, such that we span the width of the ionization front W ≈ 0.002 au with ∼17 grid points. We specify the inner boundary conditions at z = 0 to be du/dz = 0 and ρg = ρg, 0 while leaving the outer boundary conditions at z = zIF open (i.e. discretized versions of equations (15) and (21) appropriate for describing the mid-point between adjacent grid points Nz − 1 and Nz). We directly import the gas velocity, gas density, and diffusion coefficient calculated from our model (including the smoothing between disc and wind), but start the dust from rest (v = 0) with a hydrostatic density profile that extends smoothly to the outer edge of our numerical grid. We then allow the velocity and density to evolve in time until dynamic equilibrium is reached. The coloured lines in Fig. 1 show the steady-state Eulerian density and fluxes for a model with α = 0.05, dust with size |$s = 0.811\,\, \mu {\mathrm{m}}$|⁠, intrinsic grain density ρgrain = 3 g cm−3, and a mid-plane dust-to-gas ratio per local logarithmic size bin of ε = 2.53 × 10−6 (consistent with a MRN grain-size distribution spanning |$s \in [1 \,\, {\mathrm{nm}},\, 10\,\, {\mathrm{cm}}]$|⁠, discretized into 100 equally spaced logarithmic bins, and an over all mid-plane dust-to-gas ratio of 0.016). Corresponding results from our semi-analytic model calculated on a grid with Nz = 106 points are overlaid in black and show excellent agreement with the hydrodynamic solution. Note that the grain size shown in this example is such that the dust velocity changes sign below the ionization front and the net upward flux of dust at the front is driven by diffusion. Small deviations can be seen in the transition region near the ionization front (see inset panels), but are most likely caused by the different resolutions used for each model (the higher resolution in the analytic case helps to find a better fit for the purely advective solution in the wind). Even with the difference in resolution, the semi-analytic solution takes seconds to compute while the steady-state hydrodynamic solution takes days to reach a true dynamic equilibrium. Figure 1. Open in new tabDownload slide Comparison between our semi-analytic model (black lines) and the steady-state hydrodynamical Eulerian solution (coloured lines) for dust grains with size |$s = 0.811\,\, \mu {\mathrm{m}}$| at R = 10 au. Left-hand panel: the dust density as a function of height above the mid-plane. The inset panel is a zoom-in of the region z ∈ [zIF ± 3W], where zIF = 3.093 36 au, and shows that the ionization front is well resolved. Right-hand panel: the corresponding advective, diffusive, and total dust flux as a function of z. Negative fluxes are drawn with dotted lines in the Eulerian solution and dot–dashed lines in the semi-analytic solution. The inset panel again zoom-ins on the ionization front and shows the complex behaviour of the advective and diffusive fluxes near zIF. Resolving this region is crucial to finding the correct dust flux leaving the disc. Figure 1. Open in new tabDownload slide Comparison between our semi-analytic model (black lines) and the steady-state hydrodynamical Eulerian solution (coloured lines) for dust grains with size |$s = 0.811\,\, \mu {\mathrm{m}}$| at R = 10 au. Left-hand panel: the dust density as a function of height above the mid-plane. The inset panel is a zoom-in of the region z ∈ [zIF ± 3W], where zIF = 3.093 36 au, and shows that the ionization front is well resolved. Right-hand panel: the corresponding advective, diffusive, and total dust flux as a function of z. Negative fluxes are drawn with dotted lines in the Eulerian solution and dot–dashed lines in the semi-analytic solution. The inset panel again zoom-ins on the ionization front and shows the complex behaviour of the advective and diffusive fluxes near zIF. Resolving this region is crucial to finding the correct dust flux leaving the disc. 3 RESULTS We start by defining a few quantities/regimes that will aid in the analysis of our results. We then present a suite of simulations that vary grain size, turbulence, and width of the ionization front at a single radial location (R = 10 au) and mid-plane gas density (Σg = 10 g cm−2). We examine the possible dust trajectories allowed by the wind at this location and compare the range of delivered grains to the ionization front to the maximum entrainable grain size by the wind. Finally, we show how the grain properties change with disc radius and ionizing flux. Unless otherwise specified, we assume the following set of fiducial parameters: a solar-mass star with an EUV luminosity of Φ45 and an intrinsic grain density of 3 g cm−3. 3.1 Preliminaries Much of the analysis and discussion that follows is framed in terms of the dimensionless Stokes number because it both highlights the important physics of the problem and is relatively unaffected by changes to the ionizing luminosity. On the other hand, St depends on the background gas density and varies substantially from the disc mid-plane to the wind. To avoid this latter ambiguity, our discussion of St is confined almost exclusively to the the ionization front, particularly the smoothed value defined in equation (26). These values can be related back to physical grain sizes by using the secondary axes provided in Figs 2, 3, and 6. Figure 2. Open in new tabDownload slide The normalized flux (ratio of the flux of given grain size to its value if advected with the gas) as a function of grain size (upper horizontal axis) and Stokes number at the ionization front (lower horizontal axis). Each line denotes the relation for fixed value of the width, W, of the ionization front (inset legend). The vertical dashed line indicates the critical grain size (scrit) or Stokes number (Stcrit) for which the dust velocity is zero at the ionization front. The left- and right-hand panels represent α = 0.005 and 0.05, respectively. Note that the conversion between Stokes number and grain size is not constant with W, but variations for the cases shown are small (e.g. the left edge of each curve corresponds to |$s = 10^{-3}\,\, \mu {\mathrm{m}}$|⁠). Figure 2. Open in new tabDownload slide The normalized flux (ratio of the flux of given grain size to its value if advected with the gas) as a function of grain size (upper horizontal axis) and Stokes number at the ionization front (lower horizontal axis). Each line denotes the relation for fixed value of the width, W, of the ionization front (inset legend). The vertical dashed line indicates the critical grain size (scrit) or Stokes number (Stcrit) for which the dust velocity is zero at the ionization front. The left- and right-hand panels represent α = 0.005 and 0.05, respectively. Note that the conversion between Stokes number and grain size is not constant with W, but variations for the cases shown are small (e.g. the left edge of each curve corresponds to |$s = 10^{-3}\,\, \mu {\mathrm{m}}$|⁠). Figure 3. Open in new tabDownload slide The flux of dust per logarithmic grain size interval for the case of dust following an MRN size distribution in the disc mid-plane. The curves and vertical dotted line have the same meaning as in Fig. 2, while the oblique dashed line indicates the flux value in the advection limit (f = 1). The total mass flux that escapes to the wind at this radius (R = 10 au) is listed for each curve in the inset box on the right of each figure, assuming a maximum grain size of sMRN = 10 cm and a total mid-plane dust-to-gas ratio of 0.01. Figure 3. Open in new tabDownload slide The flux of dust per logarithmic grain size interval for the case of dust following an MRN size distribution in the disc mid-plane. The curves and vertical dotted line have the same meaning as in Fig. 2, while the oblique dashed line indicates the flux value in the advection limit (f = 1). The total mass flux that escapes to the wind at this radius (R = 10 au) is listed for each curve in the inset box on the right of each figure, assuming a maximum grain size of sMRN = 10 cm and a total mid-plane dust-to-gas ratio of 0.01. To help with our analysis, we define the normalized flux, $$\begin{eqnarray*} f\equiv \frac{\mathcal {F}}{\mathcal {F}_{\rm adv}} = 1 + \frac{\mathcal {F}_{\rm diff}}{\mathcal {F}_{\rm adv}}, \end{eqnarray*}$$(31) where |$\mathcal {F}_{\rm adv}$| and |$\mathcal {F}_{\rm diff}$| are the advective and diffusive components of the total flux. Importantly, f is effectively bound between [0,1] under realistic conditions. The upper bound f = 1 may correspond to the case that the dust and gas kinematics are the same (i.e. small grains with small ts) so that the dust-to-gas ratio for that grain size remains constant along the flow trajectory. However, f = 1 could also correspond to a purely advective solution (i.e. when diffusion is negligible) but with the dust and gas velocity diverging with increasing z and the dust-to-gas ratio varying so as to maintain constant flux (provided that the dust velocity remains positive at all z). Either way, we will refer to the case f = 1 as the advection limit. In the opposite extreme of f → 0, we approach what we call the diffusive limit. Here |$\mathcal {F}_{\rm diff}\sim -\mathcal {F}_{\rm adv}$| and, similar to the advective limit, can correspond to different scenarios. The first scenario is relevant to dust grains that are large enough to decouple from the gas flow in the disc such that the dust velocity becomes negative. The resulting sign reversal of |$\mathcal {F}_{\rm adv}$| first occurs at the ionization front for a critical grain size scrit (or in terms of Stokes numbers Stcrit). In a purely advective flow, grains with s > scrit would experience a runaway pile-up before ever reaching the ionization front. However, diffusion considerably smooths the advective pile-ups in these situations and takes over as the sole delivery mechanism to the ionization front. A second scenario occurs when |$\mathcal {F}_{\rm adv}$| is still positive at the ionization front, but is either weak (e.g. due to large St) or the turbulent diffusivity is strong (e.g. due to a steep gradient in the dust-to-gas ratio or large α). The fact that the dust-to-gas ratio increases with z for the purely advective solution means that the diffusive flux tends to be negative and |$\mathcal {F}\lt \mathcal {F}_{\rm adv}$|⁠. In other words, diffusion actually opposes the delivery of dust to the ionization front. In either scenario, we see the diffusion limit is characterized by a low normalized flux f. More quantitatively, we can obtain a limiting form for the diffusive limit by considering the inability of dust to accelerate as rapidly as the gas over the finite width of the ionization front. As the gas accelerates to ui across width W, the dust is accelerated by Δvi ∼ uiδt/ts, where δt ∼ W/Δvi is the time-scale for dust to cross the front. Thus, in this limit, |$\Delta v_{\mathrm{i}}\sim \sqrt{u_{\mathrm{i}}W/t_{\mathrm{s}}}$|⁠; in the absence of diffusion, this would mean that the dust-to-gas ratio was boosted by a ratio |${u_{\mathrm{i}}}/{\Delta v_{\mathrm{i}}}$| across the front in order to satisfy continuity for the dust. However, if diffusion is strong enough to iron out such a strong increase in dust-to-gas ratio across the front, this limits the value of f to $$\begin{eqnarray*} f_{\rm diff} \sim \frac{\Delta v_{\mathrm{i}}}{u_{\mathrm{i}}} = \sqrt{\frac{W \Omega _{\mathrm{K}}}{\mathrm{St}(z_{\mathrm{IF}})u_{\mathrm{i}}}}. \end{eqnarray*}$$(32) Interestingly, the solutions can approach this limit without imposing strong turbulence in the disc, provided that the front is sufficiently thin. This is because the sluggish acceleration of the dust produces a steep positive gradient in the dust-to-gas ratio at the ionization front, which in turn excites strong diffusive mixing back into the disc. It is however worth raising a caveat about the reality of the solutions in the diffusion limit that we cannot investigate further within the framework of the Dubrulle et al. (1995) formulation. When we solve the advection diffusion equation (equation 24), in conjunction with equation (21), we are assuming that, whereas the Stokes number controls the coupling of the grain motion to the mean flow, the turbulent motion of the grains is equal to that prescribed for the gas. Although we allow the Schmidt number (equation 25) to vary with Stokes number as 1 + St, this means that since St < 1 for grains entering the wind, in practice Sc ∼ 1 and thus it is assumed that the grains experience the same diffusivity as the gas. This is a reasonable assumption for many applications where the relevant time-scale associated with the turbulence is the local dynamical time, Ω−1, but in the present case, we are concerned with the value of the diffusivity over a very small spatial region through which the grains are passing on a much shorter time-scale. We thus caution that we are likely to be overestimating the effective diffusivity of the dust at the front. If the effective diffusivity should indeed be lower, we can expect the ironing out of small-scale variations in the dust-to-gas ratio at the front to be less severe (i.e. the process that is responsible for driving the solution to the diffusion limit). Exploring this issue is beyond the scope of this paper and would require an explicit modeling of dust particle motions that are only partially coupled to the turbulent gaseous background. 3.2 Delivering dust to the ionization front Fig. 2 displays the normalized flux f(s) defined in the previous section for the cases |$\alpha = [0.005,\, 0.05]$| and a range of values of W, the scale length over which the diffusivity and gas density are smoothed at the ionization front. For a fixed gas profile, the grain size (upper horizontal axis) scales linearly with the Stokes number at the ionization front, St(zIF) (lower horizontal axis), the smoothed value at z = zIF (equation 26); the vertical dotted line represents Stcrit (the corresponding grain size is denoted scrit), the maximum value of St(zIF) for which the dust velocity is upwardly directed for all z < zIF in the limit of an infinitely thin ionization front. Deviations from this infinitely thin limit are small as long as W ≪ zIF, a condition we expect to hold for real systems (e.g. equation 27). However, larger widths can experience critical sizes that are a factor of a few larger (e.g. |$s_{\mathrm{crit}}\simeq 2 \,\, \mu {\mathrm{m}}$| for W = 1 au at R = 10 au, a factor of ∼2.5 higher than our fiducial case). Fig. 2 shows that for large W and low α, f(s) is close to the advection limit and declines steeply for St(zIF) > Stcrit. In line with our earlier prediction, as diffusion increases in importance (lower W and higher α), the ironing out of the positive gradient in the dust-to-gas ratio at the ionization front acts to suppress the flow of dust. At the same time, diffusion facilitates dust transport across the front for grains with St(zIF) ≳ Stcrit, with the upper limit on grain size increasing weakly with W. We will use slimit and Stlimit to distinguish this upper size limit set by the existence of solutions with positive flux from the maximum entrainable size limit set by the wind (see Section 3.4). We provide a quantitative demonstration of the way that the fluxes vary between the advection and diffusive limits in the following section but, for now, note qualitatively how this affects the behaviour seen in Fig. 2. For high W, dust with St(zIF) < Stcrit has f ∼ 1 (efficient delivery to the ionization front). In the limit of low W, grains are not efficiently delivered even for St(zIF) ≪ Stcrit because in the diffusive limit ,|$f\propto 1/\sqrt{\mathrm{St}(z_{\mathrm{IF}})}$|⁠. In practice, therefore, the value of W influences whether grains with St(zIF) up to a few times Stcrit (i.e. of size |$s \gtrsim s_{\mathrm{crit}}= 0.8 \,\, \mu {\mathrm{m}}$| for these parameters) can reach the wind. It can be seen from Fig. 2 that the decline in f as W is decreased is stronger for higher α because this means that fluxes reach the diffusive limit at higher W. Thus, the ability of dust to reach the wind is jointly controlled by the width of the acceleration region and also by whether there are microphysical processes at the front that can mix dust across this region. 3.3 Dust delivery for an MRN dust size distribution We now consider the case that the dust in the mid-plane follows a MRN distribution (number of grains in size range s to s + ds scaling as s−3.5ds up to a maximum grain size sMRN). Assuming spherical dust grains with uniform intrinsic density, the mid-plane dust-to-gas ratio per logarithmic size bin, ε, scales as |$\sqrt{s}$|⁠. Fig. 3 depicts the flux of dust per logarithmic size interval (⁠|$\mathcal {F}$|⁠, i.e. equation 28) as a function of grain size (upper scale) and St(zIF) (lower scale) when |$\alpha = [0.005,\, 0.05]$|⁠. The heavy dashed line represents the advection limit (i.e. where all grains with St(zIF) < Stcrit are advected with a flux that is the product of the gas flux and the mid-plane dust-to-gas ratio for each size bin: |$\mathcal {F}= \varepsilon u_0\rho _{\mathrm{g,0}}\propto \sqrt{s}$|⁠). The main utility of this plot is that it shows (i) how far below the maximum (advection) limit the dust flux falls and (ii) the grain size distribution of the dust delivered to the ionization front. In each panel the various model lines correspond to the same values of W in Fig. 2. Fig. 3 confirms the behaviour described in Section 3.2 that, in the limit of large W, the dust loss is close to the advection limit and then rolls over at a grain size corresponding to a few times Stcrit (⁠|${\sim }1\,\, \mu {\mathrm{m}}$| for Φ45 and R = 10 au). As W is reduced, the mass-loss decreases across all grain sizes, but with a greater reduction for larger grains. The limiting behaviour at small W is that the mass-loss per logarithmic size bin is roughly constant for grains up to Stcrit and then cuts off abruptly for larger grain sizes. The behaviour is similar for the two α values but the α = 0.05 profiles start to deviate from the advection limit at higher W values, due to more efficient mixing across the ionization front. In Fig. 4, |$\mathcal {F}$| is plotted as a function of W for five different grain sizes. For each size, a horizontal dotted line marks the corresponding advection limit. Again it is evident that |$\mathcal {F}$| tends to the advection limit at high W and that the W value at which this transition occurs increases for higher α. In the α = 0.05 case, all grain sizes converge to the same relation at low W. Noting that in the diffusion limit (equation 32), |$f(s) \propto 1/\sqrt{\mathrm{St}(z_{\mathrm{IF}})} \propto 1/\sqrt{s}$| and that, for an MRN distribution, |$\mathcal {F}\propto \sqrt{s} f$|⁠, it follows that |$\mathcal {F}$| is independent of s in the diffusion limit and (from equation 32), scales as |$\sqrt{W}$|⁠. The convergence of |$\mathcal {F}$| on to this locus (shown as the dashed line) at low W is thus confirmation of the way that diffusion limits the dust flux for St(zIF) < Stcrit and enhances the dust flux for St(zIF) > Stcrit (e.g. |$s = 0.811 \,\, \mu {\mathrm{m}}$|⁠, distinguished from the others by a dot–dashed line). For α = 0.005, we see qualitatively similar behaviour in that |$\mathcal {F}$| falls below the advection limit as W is decreased, but has not converged to the diffusion limit for the values of W shown. The behaviour reported here, where diffusion can in some cases suppress the dust loss in the wind, is a consequence of solving equation (24) (Dubrulle et al. 1995) in the case where diffusion is effective across the narrow ionization front; we however draw attention to the discussion at the end of Section 3.1 concerning the realism of such solutions in cases where the dust is unable to couple to turbulent motions in the gas on the length scale of the front. Figure 4. Open in new tabDownload slide The flux of dust per logarithmic grain size interval for a variety of grain sizes as a function of the ionization front width W. For each grain size, the flux value corresponding to the advection limit (f = 1) is shown by a horizontal dotted line of the same colour. The oblique dashed line indicates the diffusion limit (equation 32). As expected, the flux converges to the diffusive limit at small W and large α (right-hand panel). For smaller α (left-hand panel), the flux at small W is caught in an intermediate regime that is strongly influenced by both advection and diffusion. A notable exception is |$s = 0.811 \,\, \mu {\mathrm{m}}$| (dot–dashed line), which converges to the diffusive limit regardless of α because it has a Stokes number larger than Stcrit (i.e. no advective flux at the ionization front). Figure 4. Open in new tabDownload slide The flux of dust per logarithmic grain size interval for a variety of grain sizes as a function of the ionization front width W. For each grain size, the flux value corresponding to the advection limit (f = 1) is shown by a horizontal dotted line of the same colour. The oblique dashed line indicates the diffusion limit (equation 32). As expected, the flux converges to the diffusive limit at small W and large α (right-hand panel). For smaller α (left-hand panel), the flux at small W is caught in an intermediate regime that is strongly influenced by both advection and diffusion. A notable exception is |$s = 0.811 \,\, \mu {\mathrm{m}}$| (dot–dashed line), which converges to the diffusive limit regardless of α because it has a Stokes number larger than Stcrit (i.e. no advective flux at the ionization front). 3.4 Connecting to the wind solution We examine the trajectories of dust grains that reach the ionization front by considering their entrainment in an ionized wind whose structure and kinematics is given by a variant of the self-similar wind solution of Clarke & Alexander (2016). Specifically, whereas Clarke & Alexander (2016) considered pressure-driven winds that launch from the disc mid-plane; here we consider winds that launch from the finite height above the disc at which ionization balance is achieved. In order to preserve the assumption of self-similarity it is necessary to describe the launching surface as an axisymmetric inclined plane, specified by the angle iIF between its normal and the normal to the disc mid-plane. Here we adopt iIF = 0.3 and find that the maximum launch velocity of gas above the ionization front, which allows the flow to make a smooth transition between subsonic and supersonic flow is 0.44cs, i = 4.26 km s−1. Above the ionization front, we consider the dynamical evolution of dust grains that are subject to the combination of acceleration due to gravitational, centrifugal, and drag forces. We do not assume that the grains are always traveling at their local terminal velocity (equation 23) as the lower densities in the wind do not necessarily allow application of this ‘short friction time’ assumption. The Lagrangian equation of motion in the (R, z, ϕ) directions can be written as $$\begin{eqnarray*} \frac{D v_{\rm R}}{D t} &= - \frac{v_{R}-u_{R}}{t_{\mathrm{s}}} - \frac{\mathcal {G}MR}{\left(R^2 + z^2 \right)^{3/2}} + \frac{l_\phi ^2}{R^3}, \end{eqnarray*}$$(33) $$\begin{eqnarray*} \frac{D l_\phi }{D t} &= - \frac{l_\phi -Ru_\phi }{t_{\mathrm{s}}} , \end{eqnarray*}$$(34) $$\begin{eqnarray*} \frac{D v_{\rm z}}{D t} &= - \frac{v_z-u_z}{t_{\mathrm{s}}} - \frac{\mathcal {G}Mz}{\left(R^2 + z^2 \right)^{3/2}}, \end{eqnarray*}$$(35) where lϕ is the specific angular momentum of the dust. It is assumed that each gas streamline is characterized by its specific angular momentum at the wind base. As described in Section 2.1.2, the gas emerges nearly perpendicularly to the ionization front and thus at an angle iIF to the vertical. The gas velocity at the base of the flow is a significant fraction of the sound speed in the ionized gas while the dust crosses the ionization front with a speed close to the gas flow below the front (which is subsonic with respect to the cool neutral disc gas). Consequently, the initial drag terms acting on the dust just above the front are close to those experienced by a stationary grain. Assuming vR = uR = vz = 0, the ratio of the vertical and radial components of the momentum equation is given by $$\begin{eqnarray*} \frac{D v_z}{D v_{\rm R}} = \frac{\sec ^3{i_{\rm IF}} - \chi \tan {i_{\rm IF}}}{\chi (\sec ^3{i_{\rm IF}} - 1)}, \end{eqnarray*}$$(36) where $$\begin{eqnarray*} \chi \equiv \frac{\mathcal {G}M}{R^2}\frac{t_{\rm s}}{|u|} = \mathrm{St}\frac{v_{\mathrm{K}}}{|u|}. \end{eqnarray*}$$(37) If Dvz/DvR < tan iIF, or equivalently |$\chi \gt \cot {i_{\rm IF}}$|⁠, grains promptly re-intersect the ionization front and are not entrained. For R = 10 au and iIF = 0.3, this would imply that the maximum Stokes number for which grains can be entrained in the wind occurs a little above unity (Sti = Sti, max ∼ 1.47, corresponding to a physical size of approximately 0.07 and |$7 \,\, \mu {\mathrm{m}}$| for Φ41 and Φ45, respectively). However, in the context of our smoothed ionization front, we must scale this value by ∼1/2 (since the gas velocity is approximately half of the ionized wind velocity) in order to be consistent with the Stokes numbers we measure at zIF in our model. Thus, in everything that follows, we use Stmax to refer to the scaled maximum measured at the mid-point between the disc and wind. Importantly, the flux delivered to the ionization front always cuts off below this limit, only approaching Stmax at very large W (see Figs 2 and 3). We therefore expect all grains that pass through the ionization front at R = 10 au to become entrained by the wind (we explore radial variations using equation 27 to set W in Section 3.5). For grains that are not promptly returned to the disc, we integrate equations (33)–(35) to obtain the 3D trajectories in the wind. At every time-step the value of ϕ = tan−1(z/R) is used to map the dust grain on to the appropriate gas streamline passing through this point. Pre-tabulated values of the self-similar gas streamline solution are used to identify ϕ with the distance along the streamline normalized to the base radius. Once we have the base radius and density normalization we can then calculate the local gas density, streamline inclination, and poloidal gas velocity from the self-similar solution. The local azimuthal velocity of the gas is simply calculated using conservation of specific angular momentum, assuming the gas is Keplerian at the flow base. These gas properties are then used to calculate the drag terms in the equations of motion. The left-hand panel in Fig. 5 depicts the projected R-z gas streamline and two example dust trajectories originating at R = 10 au. In order to isolate the entrainment capabilities of the wind, we assume the gas starts at the ionization front with the ionized velocity and density without smoothing (although, for consistency with the discussion of dust delivery to the front, we will continue to label dust particles in terms of their smoothed Stokes number St(zIF) defined in equation 26). The yellow dot–dashed curve represents the largest dust particle that can be launched without immediately re-intersecting the disc (⁠|$s_{\mathrm{max}}= 6.85 \,\, \mu {\mathrm{m}}$|⁠). The Stokes number at the base of the yellow curve is St(zIF) = Stmax ∼ 0.8, slightly higher than our estimate based on equations (36) and (37). The initial path skims just above the disc surface (see the inset panel) before turning upward on a trajectory that is angled about halfway between the disc surface and the gas streamline. The orange solid line is for a dust particle that is 10 times smaller (i.e. near the critical Stokes number Stcrit ∼ 0.08) and whose trajectory is more similar to the gas streamline. It then follows that the great majority of dust grains leaving the disc (St(zIF) < Stcrit) follow trajectories that are close to that of the gas. The right-hand panel shows the full 3D trajectories. The azimuthal drag provides a minor correction to the dust trajectory at each time-step that accumulates over time, but major differences do not appear until after the dust has already become unbound. This suggests that small deviations from our assumed gas flow may alter the detailed structure of the dust flow without compromising our overall conclusions. Figure 5. Open in new tabDownload slide Left-hand panel: gas streamline (blue dashed line) and trajectories for two dust sizes at R = 10 au. The orange solid line shows the trajectory of a grain near Stcrit ∼ 0.08 while the yellow dot–dashed line is the trajectory for the maximum entrainable grain size Stmax ∼ 0.8. The cross marks the sonic point for the gas and the open circles where the dust grains become unbound. The inset panel shows a square zoomed-in region of width 0.25 au near the base of the wind (black dotted line). Grains larger than Stmax immediately intersect with the disc rather than raining down at larger radii. Right-hand panel: the corresponding 3D streamline/trajectories for the gas and dust in the left-hand panel. The grey dotted lines show how the dust trajectories would change if the azimuthal drag were neglected. Figure 5. Open in new tabDownload slide Left-hand panel: gas streamline (blue dashed line) and trajectories for two dust sizes at R = 10 au. The orange solid line shows the trajectory of a grain near Stcrit ∼ 0.08 while the yellow dot–dashed line is the trajectory for the maximum entrainable grain size Stmax ∼ 0.8. The cross marks the sonic point for the gas and the open circles where the dust grains become unbound. The inset panel shows a square zoomed-in region of width 0.25 au near the base of the wind (black dotted line). Grains larger than Stmax immediately intersect with the disc rather than raining down at larger radii. Right-hand panel: the corresponding 3D streamline/trajectories for the gas and dust in the left-hand panel. The grey dotted lines show how the dust trajectories would change if the azimuthal drag were neglected. We see that the dust grains of all sizes move monotonically away from the wind base. In particular, there are no grain sizes for which grains begin to ascend and then re-descend to join the disc at larger radius. This result is to be contrasted with what is found in the case of magneto-centrifugal winds by Giacalone et al. (2019), who find rain out of grains even in the case of dust species with low Stokes number (where the short friction time assumption is valid). The reason for this difference is likely to be the very different profiles of poloidal velocity (normalized to the local escape velocity) in the two cases. In magnetocentrifugal winds, even tightly coupled grains would not achieve the escape velocity until they attained a substantial fraction of the Alfven radius, having then traversed ∼5–10 times the radius of the flow base. Grains that are somewhat less well coupled (Stokes number of order 1) cannot reach this point either because their trajectories intersect the tilted surface of the disc or because they acquire a downward velocity component. In the present case of an EUV-driven thermal wind, the gas launches more rapidly from the ionization front and a tightly coupled grain would become unbound after traversing a distance of a few 10s of |${{\ \rm per\ cent}}$| of the radius of the flow base. Even for more moderately coupled grains (Sti ∼ 1) for which the initial acceleration does not immediately return them to the disc, the significantly more rapid flow close to the flow base results in them rapidly becoming unbound. For yet larger grains (Sti ∈ [1, Sti, max]), the net acceleration vector returns them promptly below the ionization front; thus, these also do not contribute to radial transport of grains in the disc. The above calculations assume that the dust is immediately exposed to the gas properties in the ionized wind (i.e. an infinitely thin ionization front). If instead we start the dust at rest from zIF using the smoothed gas properties from our disc model, then the method of calculation has to be adjusted because the imposed structure of the front is not self-similar like the wind. To accommodate for this, at each time-step, we smooth the appropriate gas streamline in the wind to its corresponding disc value at the base using the same general formalism as our disc model (e.g. equation 26), only now the smoothing is done as a function of path length, the disc value is constant, and the wind value is variable. Depending on the assumed ionizing luminosity Φi and front width W, we find that Stmax is reduced by ∼20–|$40{{\ \rm per\ cent}}$| but that flow trajectories are otherwise little changed. One type of trajectory that is not readily resolved in our calculations are those grains that are immediately returned to the disc. These grains are the generalization of the hovering grains seen in 1D wind simulations (HPLM16a, HLM16b). It is plausible that once these grains re-enter the disc, the increased gas density will again propel them upwards, leading to an oscillation about the ionization front. However, in contrast to the hovering grains in 1D, the non-zero radial velocity would lead to radial migration along the disc surface until they either become fully entrained in the wind or sink back down into the disc interior. Although surface dust transport within the disc is potentially of interest (as in the case of transport of crystalline grains; e.g. van Boekel et al. 2005), our conclusion that dust delivery shuts off at Stokes numbers somewhat below those for which dust entrainment in the wind is ineffective leads us to expect that this effect is not significant and we do not explore it further here. 3.5 Dependence of dust wind properties on launch radius and ionizing flux We have argued above that all of the dust delivered through the ionization front at R = 10 au can be entrained in the wind. Based on the fact that Stcrit (an indicative value for the Stokes number above which dust is not efficiently transported across the ionization front) is an order of magnitude smaller than Stmax (the maximum Stokes number at the base of the ionized region for which grains are entrained in the flow), we further concluded that the the majority of dust in the wind follows trajectories that are similar to that of the gas. We now consider the radial scaling of Stmax and Stcrit to see whether these conclusions change with radius. From equations (36) and (37), we see that for a self-similar wind (where iIF and ui are independent of radius), the entrainment criterion depends only on χ and hence on the product St(zIF)vK. Since the upper limit to χ is constant, the maximum Stokes number at the ionization front scales as |$\mathrm{St}_{\mathrm{max}}\propto \sqrt{R}$|⁠. Meanwhile, we can estimate Stcrit by first calculating scrit using the short friction time limit within the disc (i.e. setting v(zIF) = 0 in equation 23 and solving for s), $$\begin{eqnarray*} s_{\mathrm{crit}}= u_{\mathrm{n}}\frac{c_{\mathrm{s}}\rho _{\mathrm{g,n}}}{\rho _{\mathrm{eff}}} \frac{\left(R^2 + z_{\mathrm{IF}}^2 \right)^{3/2}}{\mathcal {G}Mz_{\mathrm{IF}}} \propto R^{0.25} M^{-0.5}, \end{eqnarray*}$$(38) and then using scrit to calculate the mid-point between the neutral and ionized Stokes numbers: St(zIF) = 0.5(ts, n + ts, i)ΩK. Conveniently, for the case that the density profile of the ionized wind base is given by equation (14) and where ts, i ≫ ts, n, the radial scaling for Stcrit is equal to that of scrit, allowing us to plot the Stokes numbers and grain sizes together in the same figure. We take advantage of this in Fig. 6, where we compare the radial scaling of Stcrit and scrit for different ionizing luminosities Φi and turbulence levels α. Note that provided ts, i ≫ ts, n, Stcrit is independent of Φi, an approximation that is well borne out in Fig. 6. We also plot Stlimit, which is the maximum Stokes number at zIF for which our solutions in Section 3.2 yield a non-zero flux at the ionization front. In general, we find that Stcrit < Stlimit ≲ Stmax (the latter inequality being approximately equal if smoothing is factored into the maximum limit), supporting our earlier conclusion that the majority of grains leaving the disc follow similar trajectories to that of the gas for all radii relevant to photoevaporation. Figure 6. Open in new tabDownload slide Radial profiles for the three listed Stokes numbers (left-hand axis) and associated grain sizes (right-hand axis) when |$\Phi _{\mathrm{i}}= [10^{41},\, 10^{45}] \,\, {\mathrm{photons}} \,\, {\mathrm{s}}^{-1}$| and |$\alpha = [0.005,\, 0.05]$|⁠. The thick orange lines mark the critical Stokes number/grain size above which advection can no longer deliver grains to the ionization front and is adequately approximated by equation (38) (black dotted line with triangular markers), despite the simplifying assumptions that go into its derivation. The thin blue lines represent the upper limit to the Stokes number/grain size that can be delivered to the wind by diffusion in the disc. Finally, the black dashed line with square markers is the theoretical maximum entrainable Stokes number/grain size based on properties in the wind, as derived from equations (36) and (37). Figure 6. Open in new tabDownload slide Radial profiles for the three listed Stokes numbers (left-hand axis) and associated grain sizes (right-hand axis) when |$\Phi _{\mathrm{i}}= [10^{41},\, 10^{45}] \,\, {\mathrm{photons}} \,\, {\mathrm{s}}^{-1}$| and |$\alpha = [0.005,\, 0.05]$|⁠. The thick orange lines mark the critical Stokes number/grain size above which advection can no longer deliver grains to the ionization front and is adequately approximated by equation (38) (black dotted line with triangular markers), despite the simplifying assumptions that go into its derivation. The thin blue lines represent the upper limit to the Stokes number/grain size that can be delivered to the wind by diffusion in the disc. Finally, the black dashed line with square markers is the theoretical maximum entrainable Stokes number/grain size based on properties in the wind, as derived from equations (36) and (37). We see from Fig. 6 that there is a very weak dependence of Stlimit on Φi and α. The flux dependence derives from its effect on the density of the wind base and hence the value of W, (equation 27). As we saw previously in Figs 2 and 3, the width W affects Stlimit and the Stokes number corresponding to the peak dust flux (varying α produces similar effects). In the latter case, the peak flux can occur at Stokes numbers larger than Stcrit when Φi is small and R is large. As Stlimit approaches Stmax, the deviation between dust and gas trajectories in the wind becomes more significant; however, the fact that this occurs at very low Φi, where dust fluxes are in any case low, means that this is unlikely to be of observational consequence.7 Since we have shown that both Stlimit and Stmax are very insensitive to α and Φi, it follows that the corresponding grain sizes depend on Φi only via its influence on the base density of the ionized flow. This implies that both smax and scrit scale with |$\sqrt{\Phi _{\mathrm{i}}}$|⁠, which explains the |$\sqrt{\Phi _{\mathrm{i}}/\Phi _{41}}$| dependence in the right-hand axis of Fig. 6. Finally, in real systems Φi is correlated with the stellar mass. While we do not vary M in this study, we note there is an approximate |$1/\sqrt{M}$| scaling in equation (38) that results from cs ≈ ΩKH and unρg, n = constant. 4 IMPLICATIONS FOR DUSTY PHOTOEVAPORATIVE WINDS We have shown that any grain that enters the ionized flow with a Stokes number less than a threshold value Stmax will be entrained in the flow. At R = 10 au, this includes grains that are |$\lesssim 7 \,\, \mu {\mathrm{m}}$| for Φ45 (⁠|$\lesssim 0.07 \,\, \mu {\mathrm{m}}$| for Φ41). Above this threshold, dust re-enters the disc locally, although it may be possible for Stokes numbers very close to Stmax to skim along the disc surface until conditions are more favourable for entrainment (a process that we do not model here). Grains that are less than around |$10 {{\ \rm per\ cent}}$| of Stmax follow streamlines that are close to that of the gas, but all grains that are initially entrained in the wind rapidly achieve escape velocity and do not later rain out on to the disc at larger radii. We have found that the diffusive description typically used to model the turbulent transport of dust in the disc interior leads to some non-intuitive results when applied to a near-discontinuous ionization front. Namely, both the flux and size distribution of dust leaving the disc at a particular radius are sensitive to processes occurring on the length scale of the front. We presented a suite of simulations at R = 10 au corresponding to different assumed front widths that show that the efficacy of mixing at the ionization front is bracketed by two limits, the advection limit and the diffusive limit, as set out in Section 3.1. In order to focus more on parameters relevant to real discs, in Section 3.5, we switched to using equation (27) to approximate the physical width of the ionization front. At R = 10 au, we found W ≃ 0.002 au for Φ45 (appropriate for a very luminous Herbig Ae/Be star), which corresponds most closely to the purple lines in Fig. 3. Comparing the two turbulence cases reveals flux variations on the order of a few (except near the size slimit where the flux drops sharply to zero) and that slimit, which shifts from ∼2 to |$3 \,\, \mu {\mathrm{m}}$|⁠, is not strongly dependent on the level of turbulence. In both cases, the dust flux is well below the scaled value of the gas flux even for grains sizes substantially less than |$s \sim s_{\mathrm{crit}}= 0.72 \,\, \mu {\mathrm{m}}$|⁠. Further comparison with Fig. 4 suggests that, at this value of W, most grains are in the diffusive limit when α = 0.05 and an intermediate state between the diffusive and advective limits for α = 0.005. We emphasize that the suppression of the dust flux in these calculations is a consequence of the strong diffusion that results from applying the advection-diffusion equation across a very narrow front (see discussion in Section 3.1). Of course not all stars encircled by protoplanetary discs have such large EUV luminosities. Lowering the ionizing flux to Φ41 (appropriate for a low-mass T Tauri star) increases W by two orders of magnitude (e.g. W ≃ 0.15 au at R = 10 au). Importantly, however, there is a commensurate shift in the width at which the crossover between the advective and diffusive limits takes place. Therefore, the general trends and conclusions for Φ45 and Φ41 are similar, although there is some minor, order unity, variation in the radial profiles for Stlimit (see Fig. 6) and the Stokes number corresponding to the peak flux. In contrast, mapping these Stokes numbers on to grain sizes results in a |$\sqrt{\Phi _{\mathrm{i}}}$| dependence (see the right-hand scale of Fig. 6). It is also interesting to note that W inherits its radial dependence from ρg(zIF) ∝ R−1.5. It then follows that |$W/z_{\mathrm{IF}}\propto \sqrt{R}$| and that the ionization front in the inner disc is more prone to the limiting influences of diffusion than the outer disc. One caveat to the arguments above is that equation (27) is calculated using the smoothed gas density, which is much closer to ρg, i than ρg, n (a natural consequence of mass conservation caused by the rapid acceleration of the gas across the ionization front). However, since ρg(zIF) and ρg, n bracket nearly the entire drop in density across the ionization front, it stands to reason that the correct ‘physical’ interpretation is also bracketed by the effects predicted by each density. If we instead insert ρg, n into equation (27), we find that W shrinks by almost two orders of magnitude (e.g. W ≃ 2 × 10−5 au at R = 10 au) and the diffusive effects at the ionization front are magnified. Moreover, the radial dependence of W/zIF becomes a decreasing function of radius (W/zIF ∝ R−0.25), suggesting that diffusion becomes more important with R. Clearly, ρg, n predicts a more pessimistic limit on the dust flux and ρg(zIF) the more optimistic, but both show evidence of W regulating the dust flow through the ionization front. We caution once again that our results have been obtained using a simple prescription for diffusive mixing that is parametrized in terms of a turbulent mixing parameter α. Furthermore, we have assumed that α transitions between its disc value to zero over a length scale W. Clearly the application of a diffusive description to a situation where the transport properties are changing over a length scale that is much smaller than any length scales associate with disc turbulence is questionable. It is for this reason that we de-emphasize a quantitative comparison. Nevertheless, whatever the appropriate microphysical description for dust mixing in the narrow region over which the gas is accelerated at the ionization front, we argue that the diffusive limit would represent the appropriate limit if mixing prevented steep changes in the dust-to-gas ratio over the width of the front. As mentioned above, however, the attainment of this limit depends on the ability of dust grains to couple to the turbulent motions of the gas over the narrow front width, which is an imposed assumption in the formulation of Dubrulle et al. (1995). Given our ignorance of the details of microphysical mixing at the ionization front, we will focus the rest of our discussion on the results that are more robust against these uncertainties. Fig. 3 illustrates the general point that St(zIF) = Stcrit (i.e. the Stokes number above which advection can no longer supply dust to the ionization front) provides an order of magnitude indication of the Stokes number where the dust delivery turns down sharply, this limit being insensitive to the value of W.8 This allows us to make general statements about the range of dust sizes that are delivered through the ionization front. In contrast, the actual value of the dust flux (relative to the gas) is very sensitive to the details of mixing at the front. We therefore focus more on what we can say about the dust sizes that enter the wind. The lower limit to the size distribution is relatively unimportant since we expect these particles to simply be advected away with the gas. The upper limit, on the other hand, can be defined in one of two ways. In an absolute sense, the upper limit at a particular launch radius is set by the competition between delivery and entrainment (i.e. min[Stlimit, Stmax]). In practice, we find these limits are comparable – particularly when the acceleration through the ionization front is factored into the entrainment calculation. However, from Figs 2 and 3, we see that the dust flux is rapidly extinguished above Stcrit, which leads us instead to define Stcrit as an effective maximum to our grain-size distribution in our wind. When combined with the trajectory calculations performed in Fig. 5, we have reason to believe that the majority of grains delivered through the ionization front are small enough to be well entrained in the wind. Pushing this simplification one step further to assume that the dust is stuck to the gas streamlines, then gradients in the grain-size distribution and associated optical properties of the wind purely derive from the variation in maximum grain size as a function of launch radius (roughly the grain size corresponding to St(zIF) = Stcrit), as depicted in Fig. 6. The common turnover in flux near Stcrit also allows us to place some constraints on the evolution of the dust-to-gas ratio in photoevaporating discs. The lack of dust loss for St(zIF) > Stcrit means that, per unit mass of gas lost in the wind, the minimum fraction of dust that remains in the disc is the fraction of dust in the mid-plane with s > scrit; for an MRN size distribution spreading over many magnitudes in size, this fraction is given by |$1-\sqrt{s_{\mathrm{crit}}/s_{\mathrm{MRN}}}$|⁠. Given that scrit is of (sub-)micron scale and that there is strong evidence from multiwavelength sub-mm studies that grain growth has proceeded to mm sizes and above in protoplanetary discs (e.g. Testi et al. 2014), it follows that a negligible fraction of dust mass is lost to the wind. While this seems to imply that photoevaporation could be a promising mechanism for driving up the dust-to-gas ratio in protoplanetary discs (e.g. Throop & Bally 2005), it should be stressed that this conclusion only holds in regions where sMRN is large (e.g. has not been reduced by strong radial drift) and photoevaporation can remove a significant portion of the local gas mass (see Sellek, Booth & Clarke 2020 for a similar conclusion in the case of dust entrained in winds driven by external FUV radiation in the outermost regions of protoplanetary discs). Enhancements in the dust-to-gas ratio by photoevaporation, if any, may therefore be limited to localized regions of the disc, such as dust traps, that can retain large dust grains over long time-scales. Having discussed the results from our own work, it is interesting to see how these results compare with previous studies of dusty EUV winds. At a radius of |$R_{\mathcal {G}}= 9.46 \,\, {\rm{\small AU}}$|⁠, Owen et al. (2011a) found a maximum grain size in the wind of around |$2 \,\, \mu {\mathrm{m}}$|⁠. This value is within |${\sim }3 {{\ \rm per\ cent}}$| of the maximum size implied by |$\chi = \cot {i_{\rm IF}}$| (the dashed line in Fig. 6), after correcting for the factor of 3 difference in grain density and factor of 100 difference in ionizing luminosity (ρgrain = 1 g cm−3 and Φi = 1043 photons s−1, respectively). This agreement is unsurprising considering both calculations are based on the entrainment properties of the wind alone and we both source Hollenbach et al. (1994) for our ionization front density. On the other hand, the curtailment of the dust flux between scrit and slimit is more analogous to the settling limit proposed by HLM16b, only here we have generalized the effect in two key ways. First, we account for the vertical drag from the flow feeding the wind and, secondly, we use a physical model to approximate the density, location, and finite extent of the ionization front. Repeating the calculations of HLM16b using our disc parameters and ionization front location, their model predicts a maximum grain size comparable to our scrit in the inner disc and about an order of magnitude smaller than our scrit in the outer disc (e.g. |$s_{\mathrm{max}}= [0.75, \,\,0.097] \,\, \mu {\mathrm{m}}$| at R = [1, 100] au assuming Φ45). These differences show the importance of accounting for the vertical gas flow in the disc and accurately modelling the ionization front when obtaining the grain-size distribution in the wind. Ultimately, our calculations predict that the maximum grain size in photoevaporative winds is intermediate to the sizes proposed by Owen et al. (2011a) and HLM16b. 4.1 Application to winds driven by non-ionizing radiation This paper focuses on winds driven by EUV radiation where the gas at the disc-wind interface is ionized and rapidly accelerated to near sonic speeds (∼0.44cs, i = 4.26 km s−1) over a spatially thin ionization front. Using a diffusive model for turbulent diffusion, we have found that the extent to which small dust (St(zIF) ≲ 0.1) can be advected with the wind is sensitive to the treatment of mixing at the ionization front (see Fig. 3). We now consider what these calculations say about dust entrainment in other types of thermally driven photoevaporative winds. It is first important to stress that the radiation type ‘driving’ the wind is that which provides the heating in the region where the gas transitions from subsonic to supersonic flow. Therefore, for example, a disc exposed to a mixture of EUV and FUV radiation can be driven, in this sense, by FUV radiation but also be mainly heated by EUV radiation at a point far out in the flow (Johnstone et al. 1998), a phenomenon that gives rise to the well-known offset ionization fronts observed in the proplyds in the Orion Nebula Cluster and elsewhere (O’dell, Wen & Hu 1993; Kim et al. 2016). In such cases, dust entrainment and acceleration to beyond the escape velocity is achieved at points in the flow that lie far interior to the ionization front. Thus the effects observed in our simulations, where the ionization front can throttle back the entrainment of dust, are least relevant to FUV driven winds. In contrast, winds that are driven by XEUV radiation still exhibit a spatially sharp transition at the disc-wind interface (Picogna et al. 2019), potentially sharp enough to experience some of the effects observed in this study. For winds that lack a sharply defined front, we can expect the gas temperature to vary smoothly along the flow trajectory (Facchini et al. 2016; Owen et al. 2012). The resulting dust entrainment may therefore be similar to our solutions with large W where the gas is accelerated gradually over a substantial fraction of the total disc height. The larger densities and faster velocities in the acceleration region lead to an increase in both Stcrit and the flux, the latter being close to the advective limit nearly up to Stcrit before being rapidly quenched. The combination of having a larger grain-size distribution at higher fluxes may indicate that, given the same gas flux, FUV winds have the potential to be more dust rich than winds that exhibit a sharp disc–wind interface. The X-ray case has been investigated by Franz et al. (2020) for a mass-loss rate at R = 10 au that, in our model, corresponds to an EUV luminosity of ∼1043 photons s−1. The focus of their study is on the entrainment properties of the wind, which they obtain by calculating trajectories of dust particles introduced at the base of the heated flow. Since trajectories alone are insufficient to make statements regarding flux, their results are best compared with smax in our model. Accounting for the different parameters employed by Franz et al. (2020) (stellar mass, gas density and velocity, sound speed, local inclination angle of the disc surface, and intrinsic grain density), setting |$\chi = \cot {i_{\rm IF}}$| gives a maximum size of |$s_{\mathrm{max}}\sim 7.6 \,\, \mu {\mathrm{m}}$| at R = 10 au – only a little smaller than the |$9\,\, \mu {\mathrm{m}}$| they report at the same radius. The approximate agreement of smax to XEUV calculations can be attributed to the fact that (i) XEUV winds have a sharp disc-wind interface with clearly defined wind properties9 at the base of the flow and (ii) smax is a dynamical limit, derived without reference to the heating mechanism generating the wind. As a final point of interest, Franz et al. (2020) observed some of their larger grains reconnecting with the disc at large radii (R > 100 au), following a rapid decline in the base density of the wind. It is likely we do not observe such trajectories in our wind model because of the fixed power-law relation we assume for ρg, i in equation (14). Rainout of grains in photoevaporative winds may therefore only occur in regions where the radial base density profile is very steep. 5 CONCLUSIONS In this study, we have coupled a turbulent gas disc with an inclined self-similar EUV-driven photoevaporative wind and attempted to track the flow of dust as it travels from the disc mid-plane into the wind. By solving the fluid equations for the dust in the disc and the equations of motion for the dust in the wind, we are able to explore aspects of both delivery and entrainment of dust in EUV winds. HLM16b previously argued that delivery of dust to the ionization front (in their case via diffusion in a static disc) controlled the upper size limit of dust grains in photoevaporative winds for a large fraction of the disc. In this study, we find that including the advection of dust within the gas flow feeding the wind helps to improve the diffusive delivery of large grains to the ionization front, bringing the maximum deliverable size from the disc and maximum entrainable size from the wind nearly into agreement. However, while the deliverable size limit is increased, the exiting flux of grains near this limit remains small. Our simulations instead point to an effective maximum set by the steep turnover in dust flux near what we call the critical grain size, which is intermediate to the sizes proposed earlier in the literature for EUV winds (Owen et al. 2011a; HLM16b). Importantly, this turnover at the critical size limit holds over a wide range of ionizing luminosities, turbulence strengths, and ionization front widths. The critical grain size corresponds to the size at which the advective flow of dust is first interrupted by the dust velocity reversing directions (usually near the ionization front). For the fiducial case modeled here (a solar-mass star with an EUV luminosity of 1045 photons s−1, corresponding to a very luminous Herbig Ae/Be star, and an intrinsic grain density of 3 g cm−3), the critical grain size is |$\lesssim 1 \,\, \mu {\mathrm{m}}$|⁠. This is sufficiently small compared to typical range of grain sizes that are present in the disc mid-plane such that only a very small fraction of the dust mass is lost in the wind. Thus, EUV-driven photoevaporation provides a mechanism for driving up the dust-to-gas ratio in discs (for a discussion on how our results are likely to affect estimates of dust transport in photoevaporative winds driven by other types of energetic radiation, see Section 4.1). Alternatively, for EUV luminosities typical of T Tauri stars (∼1041 photons s−1) the critical grain size delivered to the ionization front is |$\lesssim 10^{-2}\,\, \mu {\mathrm{m}}$|⁠, implying that EUV winds generated by T Tauri stars should be essentially dust free. Therefore, any observations of dust in haloes around T Tauri stars would likely be due to alternative mechanisms (e.g. magneto-centrifugal winds or infall). To summarize other key results: We find that essentially all grain sizes delivered to the ionization front are entrained in the wind and eventually escape to infinity. The maximum entrainable grain size in EUV winds is set by the wind properties at the surface of the disc rather than at some later point in the wind. Thus, in contrast to hydro-magnetic winds (Giacalone et al. 2019), we do not find evidence for fall-back of grains on to the disc at larger radii for our assumed wind-base profile. However, it may be possible for steeper base density profiles, particularly in regions where photoevaporation experiences a sudden drop in efficiency (as seen in Franz et al. 2020). Taking advantage of the fact that the maximum size limit is set at the surface of the disc, we have derived a simple entrainment criterion based on whether the ratio of the vertical to radial forces result in the dust immediately intersecting with the inclined surface of the disc (⁠|$\chi \gt \cot {i_{\rm IF}}$|⁠; see Section 3.4). Importantly, this criterion is derived without reference to the heating mechanism in the wind and we find that it is in good agreement with all previous predictions in the literature. Large grains near the maximum limit do not vertically sink back into the disc but re-enter with an outward radial velocity, potentially allowing them to skim radially along the surface of the neutral wind until they are either entrained or become permanently trapped by the disc. The drag from the vertical flow created by photoevaporation aids in lofting grains into the surface layers of the disc. Even if these grains are not ultimately lost to the wind, they could enhance mechanisms such as the one recently proposed by Owen & Kollmeier (2019) regarding the removal of surface grains by radiation pressure in transition discs. For grains whose velocity remains positive throughout the flow (i.e. are always delivered to the ionization front), we find that the flux is sensitive to the efficiency of diffusive mixing in the vicinity of the ionization front. We have identified two analytical limits that explain the observed flux variation in our model: Diffusive limit: characterized by strong diffusive mixing that leads to low dust delivery rates to the ionization front. Although more prominent for large grains, small grains are also susceptible. Advective limit: characterized by weak diffusive mixing, allowing advection to deliver grains to the ionization front in the ratio in which they are present in the disc mid-plane. In deriving these limits, we have uncovered a subtle inconsistency in using standard disc turbulence models (where the relevant time-scale is the local dynamical time of the disc) to model the diffusivity of dust that passes through a spatially thin ionization front on much shorter time-scales. Further modelling of the microphysics at the ionization front is required to determine which limit is more physically plausible. However, since the fraction of the mid-plane dust that is delivered to the ionization front is low even in the advection limit (the maximum attainable flux), this uncertainty does not affect the conclusions drawn above. ACKNOWLEDGEMENTS This work has been carried out within the framework of the National Centre for Competence in Research PlanetS, supported by the Swiss National Science Foundation. We thank Gavin Coleman, Rainer Schräpler, Richard Booth, Sebastien Krijt, and James Sutherland for fruitful discussions. We also thank the anonymous referee whose insightful comments helped improve this paper. MAH acknowledges the financial support of the SNSF. CJC is grateful for funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (DUSTBUSTERS) and for support from the STFC consolidated grant ST/S000623/1. CJC also gratefully acknowledges hospitality from LMU, funded by the CAM-LMU initiative, during completion of this project. DATA AVAILABILITY The data underlying this paper will be shared on reasonable request to the corresponding author. Footnotes 1 Also note that by including the flow in the disc, the mid-plane density, surface density, and the disc scale height all deviate very slightly from their respective hydrostatic relations. However, given that the changes to the disc structure occur almost exclusively in the low-density regions near the disc surface, these deviations are negligible. 2 In connecting the disc to the wind, self-similarity is not strictly maintained. Because the disc is not globally isothermal, the radial dependence in equation (4) weakly influences the launching angle of the wind. Fortunately, the large temperature jump at the ionization front ensures that ui∥ ≪ ui⊥ and that the ionized wind emerges nearly perpendicular to the ionization front. The departure from self-similarity is small (variations of |$\lt 1{{\ \rm per\ cent}}$|⁠) and can safely be ignored. 3 Note that this also describes the plane-parallel wind solution given by Hutchison & Laibe (2016; see also HLM16b) but with different boundary conditions. 4 The approximate equality in equation (21) is due to a small simplification we have made to the drag term, which should contain a factor of 1/(1 + ε), which we set to unity on the basis of the very small dust-to-gas ratios ε in our solutions. 5 Because the solutions in this section are strictly 1D, we map the 2D gas streamline on to the vertical coordinate assuming z is equal to the distance along the streamline. In so doing, we inherently assume there are no radial changes to the gravitational force and that the dust follows the same trajectory as the gas, both of which are approximately true near the base of the wind. 6 Note that the normalization of the dust density does not affect the form of the resulting profiles, given the neglect of back-reaction of dust on the gas. 7 Another parameter modulated by Φi is the height of the ionization front zIF, but the steep, near-Gaussian density gradient at four to five disc scale heights above the mid-plane renders this effect insignificant: since |$\rho _{\mathrm{g,i}}\propto \sqrt{\Phi _{\mathrm{i}}}$|⁠, a factor of 100 reduction in ionizing luminosity only requires an order of magnitude change in ρg, i, which can be achieved by a very modest (order |$10 {{\ \rm per\ cent}}$|⁠) increase in iIF. Note that in order to maintain our assumption of constant iIF, we have absorbed the variation with Φi into the disc sound speed and scale height, as detailed in Section 2.1.4. 8 In our calculations, the peak dust flux sometimes occurs slightly above Stcrit (e.g. at small Φi and large R) but never by more than a factor of ∼3–4. 9 It would be more difficult to apply smax to FUV winds where the acceleration of the wind is more gradual. REFERENCES Adams F. C. , Hollenbach D., Laughlin G., Gorti U., 2004 , ApJ , 611 , 360 10.1086/421989 Crossref Search ADS Alexander R. D. , Pascucci I., 2012 , MNRAS , 422 , L82 10.1111/j.1745-3933.2012.01243.x Crossref Search ADS Bai X.-N. , Ye J., Goodman J., Yuan F., 2016 , ApJ , 818 , 152 10.3847/0004-637X/818/2/152 Crossref Search ADS Birnstiel T. , Fang M., Johansen A., 2016 , Space Sci. Rev. , 205 , 41 Bisbas T. G. , Bell T. A., Viti S., Yates J., Barlow M. J., 2012 , MNRAS , 427 , 2100 10.1111/j.1365-2966.2012.22077.x Crossref Search ADS Boss A. P. , 1997 , Science , 276 , 1836 10.1126/science.276.5320.1836 Crossref Search ADS Calvet N. et al. , 2005 , ApJ , 630 , L185 10.1086/491652 Crossref Search ADS Carrera D. , Gorti U., Johansen A., Davies M. B., 2017 , ApJ , 839 , 16 10.3847/1538-4357/aa6932 Crossref Search ADS Charnoz S. , Fouchet L., Aleon J., Moreira M., 2011 , ApJ , 737 , 33 10.1088/0004-637X/737/1/33 Crossref Search ADS Chiang E. , Murray-Clay R., 2007 , Nat. Phys. , 3 , 604 10.1038/nphys661 Crossref Search ADS Clarke C. J. , Alexander R. D., 2016 , MNRAS , 460 , 3044 10.1093/mnras/stw1178 Crossref Search ADS Clarke C. J. , Gendrin A., Sotomayor M., 2001 , MNRAS , 328 , 485 10.1046/j.1365-8711.2001.04891.x Crossref Search ADS Dubrulle B. , Morfill G., Sterzik M., 1995 , Icarus , 114 , 237 10.1006/icar.1995.1058 Crossref Search ADS Dullemond C. P. , Dominik C., 2005 , A&A , 434 , 971 10.1051/0004-6361:20042080 Crossref Search ADS Epstein P. S. , 1924 , Phys. Rev. , 23 , 710 10.1103/PhysRev.23.710 Crossref Search ADS Ercolano B. , Pascucci I., 2017 , R. Soc. Open Sci. , 4 , 170114 10.1098/rsos.170114 Crossref Search ADS PubMed Ercolano B. , Rosotti G., 2015 , MNRAS , 450 , 3008 10.1093/mnras/stv833 Crossref Search ADS Ercolano B. , Drake J. J., Raymond J. C., Clarke C. C., 2008 , ApJ , 688 , 398 10.1086/590490 Crossref Search ADS Ercolano B. , Koepferl C., Owen J., Robitaille T., 2015 , MNRAS , 452 , 3689 10.1093/mnras/stv1528 Crossref Search ADS Facchini S. , Clarke C. J., Bisbas T. G., 2016 , MNRAS , 457 , 3593 10.1093/mnras/stw240 Crossref Search ADS Franz R. , Picogna G., Ercolano B., Birnstiel T., 2020 , A&A , 635 , A53 10.1051/0004-6361/201936615 Crossref Search ADS Fromang S. , Nelson R. P., 2009 , A&A , 496 , 597 10.1051/0004-6361/200811220 Crossref Search ADS Fromang S. , Papaloizou J., 2006 , A&A , 452 , 751 10.1051/0004-6361:20054612 Crossref Search ADS Giacalone S. , Teitler S., Königl A., Krijt S., Ciesla F. J., 2019 , ApJ , 882 , 33 10.3847/1538-4357/ab311a Crossref Search ADS Goldreich P. , Tremaine S., 1980 , ApJ , 241 , 425 10.1086/158356 Crossref Search ADS Gorti U. , Dullemond C. P., Hollenbach D., 2009 , ApJ , 705 , 1237 10.1088/0004-637X/705/2/1237 Crossref Search ADS Haisch K. E. Jr., Lada E. A., Lada C. J., 2001 , ApJ , 553 , L153 10.1086/320685 Crossref Search ADS Hernández J. et al. , 2007 , ApJ , 671 , 1784 10.1086/522882 Crossref Search ADS Hollenbach D. , Johnstone D., Lizano S., Shu F., 1994 , ApJ , 428 , 654 10.1086/174276 Crossref Search ADS Hutchison M. A. , Laibe G., 2016 , Publ. Astron. Soc. Aust. , 33 , e014 10.1017/pasa.2016.10 Crossref Search ADS Hutchison M. A. , Price D. J., Laibe G., Maddison S. T., 2016a , MNRAS , 461 , 742 10.1093/mnras/stw1126 (HPLM16a) Crossref Search ADS Hutchison M. A. , Laibe G., Maddison S. T., 2016b , MNRAS , 463 , 2725 10.1093/mnras/stw2191 (HPLM16b) Crossref Search ADS Hutchison M. , Price D. J., Laibe G., 2018 , MNRAS , 476 , 2186 10.1093/mnras/sty367 Crossref Search ADS Ikoma M. , Nakazawa K., Emori H., 2000 , ApJ , 537 , 1013 10.1086/309050 Crossref Search ADS Inutsuka S.-i. , Machida M. N., Matsumoto T., 2010 , ApJ , 718 , L58 10.1088/2041-8205/718/2/L58 Crossref Search ADS Jennings J. , Ercolano B., Rosotti G. P., 2018 , MNRAS , 477 , 4131 10.1093/mnras/sty964 Crossref Search ADS Johansen A. , Klahr H., 2005 , ApJ , 634 , 1353 10.1086/497118 Crossref Search ADS Johansen A. , Youdin A., 2007 , ApJ , 662 , 627 10.1086/516730 Crossref Search ADS Johnstone D. , Hollenbach D., Bally J., 1998 , ApJ , 499 , 758 10.1086/305658 Crossref Search ADS Kim J. S. , Clarke C. J., Fang M., Facchini S., 2016 , ApJ , 826 , L15 10.3847/2041-8205/826/1/L15 Crossref Search ADS Krauss O. , Wurm G., Mousis O., Petit J.-M., Horner J., Alibert Y., 2007 , A&A , 462 , 977 10.1051/0004-6361:20066363 Crossref Search ADS Kuiper G. P. , 1951 , 50th Anniversary of the Yerkes Observatory and Half a Century of Progress in Astrophysics . McGraw-Hill , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Laibe G. , 2014 , MNRAS , 437 , 3037 10.1093/mnras/stt1928 Crossref Search ADS Laibe G. , Gonzalez J.-F., Maddison S. T., 2012 , A&A , 537 , A61 10.1051/0004-6361/201015349 Crossref Search ADS Lynden-Bell D. , Pringle J. E., 1974 , MNRAS , 168 , 603 10.1093/mnras/168.3.603 Crossref Search ADS Mamajek E. E. , 2009 , in Usuda T., Tamura M., Ishii M., eds, AIP Conf. Ser. Vol. 1158, Exoplanets and Disks: Their Formation and Diversity . Am. Inst. Phys , New York , p. 3 10.1063/1.3215910 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Marsh K. A. , Mahoney M. J., 1992 , ApJ , 395 , L115 10.1086/186501 Crossref Search ADS Mathis J. S. , Rumpl W., Nordsieck K. H., 1977 , ApJ , 217 , 425 10.1086/155591 Crossref Search ADS Matsuyama I. , Johnstone D., Hartmann L., 2003 , ApJ , 582 , 893 10.1086/344638 Crossref Search ADS Mordasini C. , Alibert Y., Georgy C., Dittkrist K.-M., Klahr H., Henning T., 2012 , A&A , 547 , A112 10.1051/0004-6361/201118464 Crossref Search ADS Nayakshin S. , 2010 , MNRAS , 408 , L36 10.1111/j.1745-3933.2010.00923.x Crossref Search ADS O’dell C. R. , Wen Z., Hu X., 1993 , ApJ , 410 , 696 Crossref Search ADS Owen J. E. , Kollmeier J. A., 2019 , MNRAS , 487 , 3702 10.1093/mnras/stz1591 Crossref Search ADS Owen J. E. , Ercolano B., Clarke C. J., 2011a , MNRAS , 411 , 1104 10.1111/j.1365-2966.2010.17750.x Crossref Search ADS Owen J. E. , Ercolano B., Clarke C. J., 2011b , MNRAS , 412 , 13 10.1111/j.1365-2966.2010.17818.x Crossref Search ADS Owen J. E. , Clarke C. J., Ercolano B., 2012 , MNRAS , 422 , 1880 10.1111/j.1365-2966.2011.20337.x Crossref Search ADS Picogna G. , Ercolano B., Owen J. E., Weber M. L., 2019 , MNRAS , 487 , 691 10.1093/mnras/stz1166 Crossref Search ADS Pollack J. B. , Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996 , Icarus , 124 , 62 10.1006/icar.1996.0190 Crossref Search ADS Quintana E. V. , Adams F. C., Lissauer J. J., Chambers J. E., 2007 , ApJ , 660 , 807 10.1086/512542 Crossref Search ADS Riols A. , Lesur G., 2018 , A&A , 617 , A117 10.1051/0004-6361/201833212 Crossref Search ADS Safronov V. S. , 1969 , Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets . Keter Publishing House , Jerusalem Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Schräpler R. , Henning T., 2004 , ApJ , 614 , 960 10.1086/423831 Crossref Search ADS Sellek A. D. , Booth R. A., Clarke C. J., 2020 , MNRAS , 492 , 1279 10.1093/mnras/stz3528 Crossref Search ADS Shakura N. I. , Sunyaev R. A., 1973 , A&A , 24 , 337 Shampine L. F. , Reichelt M. W., 1997 , SIAM J. Sci. Comput. , 18 , 1 Crossref Search ADS Shampine L. F. , Reichelt M. W., Kierzenka J. A., 1999 , SIAM Review , 41 , 538 10.1137/S003614459933425X Crossref Search ADS Shi J.-M. , Chiang E., 2014 , ApJ , 789 , 34 10.1088/0004-637X/789/1/34 Crossref Search ADS Suzuki T. K. , Inutsuka S.-i., 2009 , ApJ , 691 , L49 10.1088/0004-637X/691/1/L49 Crossref Search ADS Takeuchi T. , Lin D. N. C., 2002 , ApJ , 581 , 1344 10.1086/344437 Crossref Search ADS Takeuchi T. , Clarke C. J., Lin D. N. C., 2005 , ApJ , 627 , 286 10.1086/430393 Crossref Search ADS Tanaka H. , Takeuchi T., Ward W. R., 2002 , ApJ , 565 , 1257 10.1086/324713 Crossref Search ADS Testi L. et al. , 2014 , Protostars and Planets VI . University of Arizona Press 10.2458/azu_uapress_9780816531240-ch015 Crossref Throop H. B. , Bally J., 2005 , ApJ , 623 , L149 10.1086/430272 Crossref Search ADS van Boekel R. , Min M., Waters L. B. F. M., de Koter A., Dominik C., van den Ancker M. E., Bouwman J., 2005 , A&A , 437 , 189 10.1051/0004-6361:20042339 Crossref Search ADS Ward W. R. , 1986 , Icarus , 67 , 164 10.1016/0019-1035(86)90182-X Crossref Search ADS Ward W. R. , 1997 , Icarus , 126 , 261 10.1006/icar.1996.5647 Crossref Search ADS Youdin A. N. , Goodman J., 2005 , ApJ , 620 , 459 10.1086/426895 Crossref Search ADS Youdin A. , Johansen A., 2007 , ApJ , 662 , 613 10.1086/516729 Crossref Search ADS Youdin A. N. , Lithwick Y., 2007 , Icarus , 192 , 588 10.1016/j.icarus.2007.07.012 Crossref Search ADS APPENDIX A: DERIVATION OF THE JUMP CONDITIONS In our model, the gas within the disc flows vertically until it reaches an inclined ionization front that separates the neutral disc from the ionized wind. Across the front, the perpendicular velocity components un⊥ and ui⊥ are constrained by the Rankine–Hugoniot equations in equations (5) and (6), while the parallel components remain unchanged (i.e. ui∥ = un∥ = u∥). The inclination angle iIF, the sound speed in the disc (cs) and wind (cs, i), and the ionized gas velocity (ui) and density (ρg, i) at the base of the wind are all known parameters of the model. This leaves a total of five unknown variables that need to be determined: un, un⊥, u∥, ui⊥, and ρg, n. Together with the jump conditions, we close the system of equations by using trigonometric relations between the velocity components: $$\begin{eqnarray*} u_{\mathrm{i}\perp }^2 & = u_{\mathrm{i}}^2 - u_\parallel ^2, \end{eqnarray*}$$(A1) $$\begin{eqnarray*} u_{\mathrm{n}\perp }& = u_\parallel \cot {i_{\rm IF}}, \end{eqnarray*}$$(A2) $$\begin{eqnarray*} u_{\mathrm{n}}& = u_\parallel \csc {i_{\rm IF}}. \end{eqnarray*}$$(A3) Using equation (5) to replace ρg, n in equation (6), we obtain a quadratic equation for un⊥, $$\begin{eqnarray*} u_{\mathrm{i}\perp }u_{\mathrm{n}\perp }^2 - \left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2\right) u_{\mathrm{n}\perp }+ c_{\mathrm{s}}^2 u_{\mathrm{i}\perp }= 0, \end{eqnarray*}$$(A4) which is satisfied by equation (7). Substituting this relation for un⊥ back into equation (5) and solving for ρg, n yields (after rationalizing the denominator) the second equality in equation (8). Meanwhile, by setting equation (7) equal to equation (A2) and isolating the radical, we find $$\begin{eqnarray*} c_{\mathrm{s,i}}^2+u_{\mathrm{i}\perp }^2 -2 u_{\mathrm{i}\perp }u_\parallel \cot {i_{\rm IF}} = \sqrt{\left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \right)^2 - 4 c_{\mathrm{s}}^2 u_{\mathrm{i}\perp }^2}. \end{eqnarray*}$$(A5) Squaring both sides, canceling like terms, and rearranging those remaining gives $$\begin{eqnarray*} \left(c_{\mathrm{s,i}}^2 + u_{\mathrm{i}\perp }^2 \right) u_\parallel \cot {i_{\rm IF}} = u_{\mathrm{i}\perp }\left(c_{\mathrm{s}}^2 + u_\parallel ^2 \cot ^2{i_{\rm IF}} \right). \end{eqnarray*}$$(A6) Squaring the equation again, we then use equation (A1) to eliminate ui⊥ and acquire an equation for u∥ in terms of known variables. Finally, by expanding bracketed terms, collecting powers of u∥, and simplifying the coefficients, we arrive at equations (9)–(13) from the main text. Once u∥ is obtained, the remaining quantities are easily computed from the constraint equations. © 2020 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Dust delivery and entrainment in photoevaporative winds
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/staa3608
DA - 2021-01-05
UR - https://www.deepdyve.com/lp/oxford-university-press/dust-delivery-and-entrainment-in-photoevaporative-winds-M0XiXqFX0b
SP - 1127
EP - 1142
VL - 501
IS - 1
DP - DeepDyve
ER -