Modeling the delivery of dust from discs to ionized winds

A necessary first step for dust removal in protoplanetary disc winds is the delivery of dust from the disc to the wind. In the case of ionized winds, the disc and wind are sharply delineated by a narrow ionization front where the gas density and temperature vary by more than an order of magnitude. Using a novel method that is able to model the transport of dust across the ionization front in the presence of disc turbulence, we revisit the problem of dust delivery. Our results show that the delivery of dust to the wind is determined by the vertical gas flow through the disc induced by the mass loss, rather than turbulent diffusion (unless the turbulence is strong, i.e. $\alpha \gtrsim 0.01$). Using these results we provide a simple relation between the maximum size of particle that can be delivered to the wind and the local mass-loss rate per unit area from the wind. This relation is independent of the physical origin of the wind and predicts typical sizes in the 0.01 -- $1\,\mu m$ range for EUV or X-ray driven winds. These values are a factor $\sim 10$ smaller than those obtained when considering only whether the wind is able to carry away the grains.


INTRODUCTION
Mass-loss in protoplanetary disc through winds is important for understanding their evolution. In particular, photoevaporative winds driven by either X-rays, extreme, or far ultra-violet (EUV, FUV) radiation are thought to be responsible for the final rapid clearing of protoplanetary discs (e.g. Clarke et al. 2001;Owen et al. 2011b;Ercolano et al. 2015;Gorti et al. 2016). More recently, magnetohydrodyanmic (MHD) winds have replaced turbulence as being the most promising processes responsible for driving accretion in protoplanetary discs (Salmeron et al. 2007;Suzuki & Inutsuka 2009;Bai 2017;Béthune et al. 2017). As a result, the entrainment of dust in these winds has also become an important issue. Throop & Bally (2005) suggested that the preferential removal of gas by winds might aid planet formation although more recent studies (Ercolano et al. 2017;Sellek et al. 2020) suggest that this is unlikely unless radial drift of dust can be suppressed. In addition, the dust entrained in winds may provide a way to probe them observationally (e.g. Owen et al. 2011a;Miotello et al. 2012;Franz et al. 2020).
The problem of dust entrainment may be thought of in two parts. Firstly, there is the question of whether dust particles entering the wind region are sufficiently well coupled to the gas so as to be carried away by the wind. However, more importantly, such escaping dust grains also need to be delivered to the wind from the underlying disc. In the case of winds driven by EUV radiation, the wind and disc are sharply delineated by a narrow ionization front where the gas density and temperature vary by many orders of magnitude. Previously, Hutchison et al. (2016); Hutchison & Clarke (2020) have E-mail: r.booth@imperial.ac.uk found that delivery is the limiting step in controlling the range of dust sizes that are lost in such winds.
In some previous studies, (e.g. Hutchison et al. 2016) it has been assumed that dust delivery to the wind occurs diffusively, with turbulence in the disc competing against settling due to gravity to loft the grains into the wind. However, it has recently become clear that dust may instead be delivered to the wind advectively, through coupling of the motion of dust grains to the upward motion of gas towards the ionization front. Although this gas motion is strongly subsonic (likely even below the typical turbulent speeds), Hutchison & Clarke (2020) argued that advection is an important component in the delivery of dust to the wind. There is precedence for this, since in their MHD simulations of disc winds Riols & Lesur (2018) showed that advection resulted in an increase of the dust scale height over the height expected from purely turbulent transport.
However, Hutchison & Clarke (2020) encountered a problem in quantifying how efficiently dust is delivered to the wind. This problem was associated with the non-convergence of their results as the width of the ionization front was reduced. This non-convergence, as discussed in Hutchison & Clarke (2020), results from the fact that there is a steep gradient in gas density at the ionization front, where the gas goes from being cold in the disc to hot in the ionized wind. Although the gas velocity changes rapidly across the ionization front, changes in the dust velocity occur on a 'stopping time', ts, the time over which drag forces act. This means that the dust density varies over the 'stopping length', vts. Since the stopping length can be much larger than the ionization front width, this leads to a steep increase of the dust-to-gas ratio across the ionization front. Hutchison & Clarke (2020) modelled the effects of turbulence as a diffusion equation (following Dubrulle et al. 1995) so that this large gradient in dust-to-gas ratio produced a large negative diffusive flux. This flux then suppresses the delivery of dust to the wind, by an amount that depends on the width of the ionization front.
However, this approach is not fully consistent because the diffusion is ultimately driven by the coupling of the dust dynamics to turbulent gas motions via drag forces. Diffusive motions are therefore subject to the same constraints as the mean flow in being limited by the finite coupling between dust and gas. A reduction of the effective difffusion coefficient, in cases where the gas flow changes on scales less than the stopping distance, is however not captured by the formulation of Dubrulle et al. (1995), which therefore gives erroneous results in this limit.
The goal of this paper is to rectify this deficiency in modeling dust transport across narrow fronts, thus determining how efficiently dust is delivered to the wind. Instead of solving an advection-diffusion equation, we use a Monte-Carlo model to trace the dynamics of individual dust grains, explicitly treating the coupling of the dust to turbulent velocity fluctuations in the disc gas. This approach is similar to the one used by Youdin & Lithwick (2007) to model the diffusion of large dust particles in discs, for example. We present our model in section 2, and in section 3 demonstrate that, in contrast to other formulations for modeling dust in turbulent flows in the literature, we are able to correctly recover the structure of the gas and dust across steep transitions in the gas density, an important prerequisite for tackling problems involving ionization fronts. In section 4 we consider the two criteria suggested by Hutchison & Clarke (2020) as limiting the maximum size of grains that are a) deliverable to the ionization front and b) entrainable by the wind above the ionization front, for which the corresponding Stokes numbers (evaluated just below the ionization front) are denoted by Stcrit and Stmax respectively). We then use these limits to estimate the level of turbulence at which a transition between diffusive and advective feeding of dust into the wind base is expected. In section 5 we demonstrate that, as anticipated by Hutchison & Clarke (2020), Stcrit (which turns out to be ∼ 0.01 for a wide range of input parameters) indeed represents a good limit for setting the maximum size of dust delivered into the wind: although somewhat larger dust grains may enter the wind in the limit of strong turbulence, the grains entering the wind have Stokes number significantly less than Stmax and hence are all capable of being fully entrained in the ionized flow. A discussion of our results and conclusions are presented in section 6 and section 7.

Modeling the disc/wind base: gas
We model the entrainment of dust in a background disc undergoing photoevaporation. The vertical structure of the disc is computed by solving the momentum equation of hydrodynamics in one dimension, assuming steady-state such that ρuz is constant. Here uz is the gas velocity, and ρ is the gas density.
The sound-speed profile, cs(z), is chosen to model the transition from a cold disc to a hot photoionized wind at the ionization front, Here disc is assumed to be vertically isothermal, where c s,disc and c s,wind denote the sound speed in the disc and wind, with zIF and W specifying the location and width of the transition. By default we take M = 1 M and R = 10 au. We assume the disc aspect ratio is given by H/R = 0.05(R/au) 0.25 , thus c s,mid ≈ 0.84 km s −1 . The sound speed in the wind, c s,wind = 12.85 km s −1 , appropriate for a fully ionized hydrogen gas at 10 4 K. Where required, the mid-plane density is taken to be ρ0 = Σ(R)/ (2π)H, assuming the gas surface density Σ(R) = 30(R/au) −1 g cm −3 .
To set the mid-plane velocity, we assume that photoionization drives an outflow with a velocity of v wind = 0.5c s,wind at the ionization front, this being motivated by typical launch velocities for self-similar solutions for isothermal winds (Clarke & Alexander 2016). Explicitly, we find the velocity at z = 0 iteratively by integrating Equation 1 to z = zIF + 9W and requiring that vz at this point is 0.5c s,wind . The parameter zIF controls the density at the base of the ionized wind ρion, once the mid-plane density and temperature of the disc are assigned. Since the disc below the ionization front is very close to a state of hydrostatic equilibrium, ρion ∼ ρ0 exp(−z 2 IF /2H 2 )fi, where fi is the factor by which the density drops across the ionization front: fi ∼ c 2 s,mid /(c 2 s,wind +v 2 wind ). The canonical parameters detailed above and ρion derived from Equation 23 for an ionizing flux of 10 42 s −1 corrsponds to zIF ∼ 4H. Assuming the standard profiles for photoevaporative winds driven by EUV radiation (Hollenbach et al. 1994), these values would correspond to integrated mass loss rates from the disc of order 10 −10 − 10 −9 M yr −1 .
Equation 1 is solved numerically using the 4th-order Runge-Kutta method of Dormand & Prince (e.g. Press et al. 2007) as implemented in the package in the library 1 . The solution at intermediate points is then obtained via piecewise-cubic Hermite interpolation (Fritsch & Carlson 1980).

Modeling the disc/wind base: dust
The dust component is treated using a Stochastic Lagrangian Model. We compute the trajectory of a large number of tracer particles under the action of gravity, coupled to the gas via drag forces: where vz is the vertical velocity of the dust. Here we have decomposed the gas velocity into its background component and a fluctuating part, u z (z, t), which represents the motions due to turbulence in the disc that are responsible for diffusion. We will assume that the turbulent fluctuations are Gaussian in nature with a correlation time, te. We assume linear Epstein drag such that, where ρgrain is the internal density of a dust grain and s is its size. Typically, we label grain size by their Stokes number, St = tsΩ, but where relevant we will assume ρgrain = 1 g cm −3 . The background gas velocity, uz(z), sound speed, cs, and density, ρ, are taken from the model described in subsection 2.1. The turbulent fluctuations are treated using a Langevin model based on Thomson (1984) (see also Wilson et al. 1983). Thomson (1984) derived a Stochastic Lagrangian model for the motion of tracer particles in the atmosphere by requiring that the statistical distribution of the particles must be the same as that of the underlying atmosphere. Thomson (1987) showed that this requirement -that the particles must remain 'well mixed' with the gas -is rather general, with Stochastic Lagrangian models that satisfy this criterion being consistent with the Euler equations and able to reproduce both the short and long term behaviour of the gas. Under the assumption of Gaussian velocity fluctuations, the model for the gas is where The w hs term corrects for the fact that that a Gaussian distribution of turbulent velocities with mean zero will drive a non-zero net flux when there are gradients in ρ(z) or σ(z) (Thomson 1984;Ciesla 2010). Here te is the Lagrangian correlation time, σ(z, t) 2 is the variance of the velocity fluctuations and W is Wiener Process, i.e. δW is a random number distributed as δW ∼ N (0, δt) (where N (0, δt) is a normal distribution with mean zero and variance δt). The diffusion coefficient D is linked to te and σ via D = σ 2 te (e.g. where Ω is the Keplerian frequency and D = αc 2 s,disc Ω −1 , giving σ = √ αc s,disc so that the dimensionless parameter α relates the sound speed to the turbulent velocity as in the viscous α parametrisation of Shakura & Sunyaev (1973). We extend this model to treat dust grains in the simplest way possible, which is to use take u z (z) = σ(z, t)wt and use the particle's position to define z in Equation 7. Explicitly, we use which reduces to Thomson (1984)'s model in the limit tstop → 0. Our model is similar to, but differs from, existing Stochastic Lagrangian Models for dust in the literature. For uz = 0 and constant ρ, our model reduces to that of Youdin & Lithwick (2007). The model of Ormel & Liu (2018) is the most similar to ours, differing by the way in which the correction term w hs is implemented. Ormel & Liu (2018) add a term σw hs to uz in Equation 10 while neglecting w hs in Equation 11. When w hs is slowly varying the effect of this on the dynamics is small; however, in the presence of a sharp transition in the density (or turbulence), as is the case at an ionization front, the difference becomes significant. Since we apply the correction in Equation 11 the effects of steep transition in density are averaged over te, whereas in the case of Ormel & Liu (2018) they are applied locally. The model proposed by Laibe et al. (2020) is equivalent to assuming w hs = 0 in Equation 11. These differences are highlighted in section 3 2 .
2 Ciesla (2010) also provided a Stochastic Lagrangian Model that satisfies Thomson (1987)'s well-mixed condition. However, Ciesla (2010) used the terminal velocity approximation for the mean flow and, by imposing fixed (i.e. tstop independent) velocity impulses, neglected the finite Lagrangian correlation time of the turbulence, making it inappropriate for the ionization front problem. Thus we do not consider it further.

Implementation
We have adopted a semi-implicit approach to solve equations 9-11 efficiently for particles with short stopping times. First, we move the particle from z to z + vzδt. Next we update wt, evaluating the righthand side of Equation 11 at the new position. Finally, we update the dust velocity, vz, using the new position and wt. This update is done implicitly in vz to avoid limits on the time-step due to small tstop, such that It is straightforward to verify that this expression is correct in the limits tstop → 0 and tstop → ∞.
The time-step, δt, is chosen to satisfy a number of constraints: The first three constraints are designed to ensure that wt and w hs do not change significantly in one time-step, while the last one is included to make sure that the particles do not jump across the ionization front in a single time-step. Rather than setting δt = δtmax, we instead set δt = 2 − max Ω −1 , where max = 63 and is the largest integer such that δt ≤ δtmax. Particle time steps are always allowed to decrease; however, increases are only allowed if the particle would remain synchronised (i.e. t is exactly divisible by δt). This decision is made to ensure that all particles are synchronised every te, so that their positions and velocities can be sampled at the same time.

TESTS
To test the code, we compute the distribution of particles in a Gaussian disc with a constant sound speed and uz = 0, comparing the results to the methods proposed by Ormel & Liu (2018) and Laibe et al. (2020). For each test, 10 3 particles were injected at z=0. After a burn-in period of 10 4 Ω −1 the positions of each particle in the range [−5, 5] (for St = 0 and [−2, 2] for St = 0.05) were recorded every 10Ω −1 for the next 10 5 Ω −1 . The density was computed by binning the particles into 100 bins, normalised such that the total mass is 1.
First we consider particles with St = 0, which should be distributed with the same density as the gas. Figure 1 (left panel) shows this for the case of a constant α = 0.01. Here both our method and Ormel & Liu (2018) 3 produce similar results, with the particles well-mixed with the gas. However, the method of Laibe et al. (2020) produces a constant dust density rather than dust-to-gas ratio, which is a consequence of neglecting the w hs term.
Next, we consider the same model but with α varying from 0.01 to 0.1, using the functional form of Equation 2 with zIF = H and W = H/60. Note that cs is constant in this test. Under these conditions particles with St = 0 should still have the same density distribution as the gas. Now we see a difference between the method of Ormel & Liu (2018) and the one presented in this work, with only our method producing a constant dust-to-gas ratio. The 'blip' in the density around z = H produced by Ormel & Liu (2018)'s method arises because w hs is large at this location. When applied in Equation 10, as in Ormel & Liu (2018), this leads to large velocities for the dust particles, which are responsible for the 'blip'. In our method, these large velocities are not produced because w hs gets averaged over te. We note that these differences are only significant if w hs varies over a small length scale -this was not the case in the tests presented in Ormel & Liu (2018), but such variations do occur close to the ionization front in our model. Again, the method of Laibe et al. (2020) produces constant dust density in regions where α is constant.
The right panel of Figure 1 shows a repeat of the test with α varying with height for particles with St = 0.05 at z = 0. In this case all of the methods produce similar results close to z = 0, which are in good agreement with the analytical solution of Youdin & Lithwick (2007, away from the mid-plane the analytical solution is no longer valid). Again the Ormel & Liu (2018) method shows an artefact at the transition, while in this case settling reduces the difference between our method and that of Laibe et al. (2020).
As a final confirmation of the ability of our code to deal with sharp gradients in the density we show in Figure 2 the density of St = 0 particles in a ionization front test with a width, W = 10 −5 au (and H ≈ 0.89 au). This shows excellent agreement with the background profile, as expected.

ANALYTIC ESTIMATES
Here we provide some simple estimates of the maximum size of dust grains that can be entrained in the wind. For sufficiently weak turbulence, the delivery of dust to the ionization front can be estimated by neglecting the u z term in Equation 4: in this case the passage of dust from the disc to the wind is simply set by the extent to which grains can couple to the advective flow of gas in the disc induced in response to mass loss at the ionization front. The maximum size of particles delivered to the wind may be estimated as the biggest particle for which vz > 0 at the ionization front, which, in the terminal velocity limit may be written as as suggested by Hutchison & Clarke (2020). Note that this refers to the Stokes number measured just below the ionization front. Since uz cs in the disc, Stcrit can be approximated as: where we have used the Rankine-Hugoniot relations to relate the velocity in the disc to the Mach number at the base of the wind, Mw = 0.5. For typical values of c s,disc and c s,wind , we find Stcrit ≈ 0.01. For comparison, the maximum particle size that, once in the wind, can escape from the disc is approximately given by the particle size for which the terminal velocity is zero at the base of the wind, i.e. downwind of the ionization front 4 . Following, Hutchison & Clarke (2020), we refer to this as Stmax, which is given by We emphasise that although Stcrit and Stmax relate to situations of force balance applied on either side of the ionization front, they relate to stopping times that are evaluated at the same location (i.e. just below the ionization front) and their ratio therefore directly relates to the ratio of dust sizes that achieve this condition on each side of the front. Comparison of equations (14) and (15) immediately shows that in the case of dust that is advected through an ionization front, the grains that are just able to reach the ionization front are a factor ∼ c s,disc c s,wind in size below the maximum size that can be entrained in the ionized wind. Thus delivery of grains to the ionization front is the limiting step in removing dust from the disc rather than the subsequent ability of the ionized wind to carry it away (Hutchison & Clarke 2020).
When turbulence is strong, we expect that dust may be delivered to the ionization front diffusively instead of being delivered by advection. Neglecting the contribution from advection (i.e. setting u(z) = 0), the dust density is given by (Dubrulle et al. 1995;Takeuchi & Lin 2002). Therefore, in a purely diffusive disc the delivery of dust to the ionization front should drop once St(zIF) > α (note we have neglected the influence of the ionization front itself). Since advection can efficiently supply dust to the ionization front for sizes below Stcrit, we therefore expect the transition to the diffusive regime to occur at α ∼ Stcrit ≈ 0.01.

Advective and diffusive dust delivery
For the results presented in this section we set up the gas profile according to the description in subsection 2.1. Dust particles are then injected continuously at z = 0 at a rate of 4Ω and the simulation is run for 10 5 Ω −1 . For the boundary conditions, particles are removed once they cross z = 5H while at z = 0 we use reflecting boundaries (i.e. particles that cross z = 0 have the sign of z, vz and wt flipped). The mass-loss time-scale of dust particles was computed by comparing the rate of particle injection to the total number of particles in the domain once the simulation has reached a steady state. Since for very low mass-loss rates (i.e. for St > Stcrit) steady state is not reached within 10 5 Ω −1 we instead fit a model to the total number of particles in the domain over time using least-squares (via 's _ routine 5 ). The model we use is whereṄ0 = 4Ω. We then compare τ (i.e. N/Ṅ0 in steady-state) to the mass-loss time-scale of the gas Σ/Σ to determine the efficiency of dust entrainment in the wind, With this definition, the mass-loss rate of dust is simply the product of F, the mass loss rate of gas and the dust to gas ratio. We note that this definition of the entrainment efficiency is slightly different to the definition used by Hutchison & Clarke (2020), who useḋ N0/(ρ d (z = 0)uz(z = 0)) (where ρ d is normalised such that the total mass is 1). This choice is dictated by practicality: the definition used here is easier to determine accurately. However, the estimates only differ by a factor (Σg/Σ d )(ρ d (z = 0)/ρg(z = 0)) ≈ 1.
For the simulations presented in this section, we choose zIF in the approximate range of 3H -4H as a compromise between realistic mass-loss rates and computational expediency. Note the mass-loss time-scale, Σ/Σ, is (approximately) given bẏ Therefore zIF = 4H corresponds to a reasonable mass-loss timescale of 10 6 yr at R = 10 au. The efficiency of dust entraiment for models with R = 10 au and zIF = 3.5 H are shown in Figure 3 for α = 5 × 10 −4 to 5 × 10 −2 and a range of ionization front widths, W .
For small α, we find a flux efficiency F ≈ 1 for small St, transitioning to F ≈ 0 at St ≈ Stcrit for all but the largest of ionization front widths. This is the expected result for advection-dominated delivery of dust to the wind. Comparing the α = 5 × 10 −4 results to those α = 5 × 10 −3 shows that increasing α mildly increases F for St close to Stcrit, but in both cases Stcrit remains a good estimator of the maximum dust size that can be entrained.
For α = 0.05 the delivery of dust to the ionization front is now diffusion dominated according to our estimate in section 4 (since α > Stcrit), resulting in different behaviour. For St Stcrit, the mass-loss rate of dust is now higher than that of the gas and also dependent on the width of the ionization front. We also find that for sufficiently narrow ionization fronts the dust flux eventually converges. Convergence occurs at progressively smaller ionization front widths as the Stokes number decreases. This convergence can be explained by considering the stopping distance at the two different sides of the ionization front: where in both cases the Stokes number is measured in the disc immediately before the ionization front. As the width of the ionization front is successively decreased, particles with a given Stokes number will first decouple on the down-wind side of the ionization front. Finally once W l disc ≈ 0.03StH, the particles will cross the ionization front before being able to react, and thus the width of the ionization can no longer affect the flux. Figure 3 confirms this, with convergence by W ≈ 0.1l disc .
The dependence of the flux on α can be understood by looking at the dust-to-gas ratio profiles, shown for two examples in Figure 2. At α = 5×10 −4 and St mid = 10 −5 , the dust-to-gas ratio increases with z between the mid-plane and the ionization front. This follows from mass conservation and the fact that dust velocity is lower than the gas velocity since the gravitational acceleration is not nearly balanced by pressure, as in the case of the gas, and finite tstop prevents the dust from keeping up with the gas flow. However, already for this low α the dust-to-gas ratio gradient is lower than predicted by the purely advective regime. This is the result of turbulent diffusion and acts to reduce the mass flux (as can be seen from Figure 3 and also the dust-to-gas ratio being below 1 at z = 5H).
Increasing α, one would expect the dust-to-gas ratio to be driven towards a constant value even more strongly. However, in the simulation with α = 0.05, we see a negative gradient in dust-to-gas ratio, evidence that diffusion is now driving an outward flux of dust (which is further supported by the average dust velocity in Figure 2 being larger than the mean gas velocity). This gradient is particularly strong close to the ionization front. Our proposed explanation for this is that the diffusive supply of dust to the ionization front from the disc is not matched by a return flux from the wind since the low density downwind of the front means that such particles are not effectively coupled to turbulent motions driving them back through the front. Conversely, for particles small enough that they remain coupled to the gas through the ionization front, the diffusive flux is cancelled by the w hs correction term resulting in F ≈ 1.
Even in the diffusive regime the mass-loss rate of dust becomes  In each case the Stokes number shown is the Stokes number in the disc as measured at the ionization front.
negligible once the Stokes number exceeds Stcrit by more than a factor of ∼ 10. For grains of this size, St > Stmax, and the drag force on particles that pass through the ionization front is no longer sufficient to overcome gravity in the wind. Therefore, even if dust can be supplied to the wind, ultimately it can not escape the disc. The results are not sensitive to the underlying parameters of the disc. This is demonstrated in Appendix A and Figure 4, where we compare the Stokes number at which F drops to 0.5 to the values of Stcrit and Stmax. We do this varying zIF at R = 10 au (left panel) and also for a model where zIF is computed according to Equation 23 for Φ = 10 42 s −1 . The EUV model follows Hutchison & Clarke (2020), who assume the density at the wind base is controlled by recombination: where the velocity at the wind base is 0.5c s,wind as before. Here mH is the mass of a hydrogen atom, α2 = 2.6 × 10 −13 cm 3 s −1 is the Case B recombination coefficient and Φ is the stellar EUV luminosity.
In Figure 4 we see that the maximum size entrained is close to Stcrit for α < Stcrit independent of the height of the ionization front or the location in the disc. Furthermore, although the size increases above Stcrit for α > Stcrit, it always remains smaller than Stmax.

Typical grain sizes entrained
Now that we have ascertained that the maximum size of dust grains delivered to the wind is determined by Stcrit in the advective regime while Stmax limits the size of grains removed in the diffusive regime, we consider what these Stokes numbers mean in terms of grain size. From the definition of these limits (i.e. zero acceleration of a stationary grain just below and just above the ionization front respectively) we can write scrit = 8 πΣ ρgrainΩ This shows that maximum size of dust particle that can be entrained is insensitive to the disc mass, which only enters through the dependence of zIF on disc mass, which is weak. Note that this equation is valid even if the disc is not vertically isothermal (as assumed in this paper) as long as HIF is determined from c s,disc measured at the ionization front. Similarly the definition for smax follows by replacing c s,disc with c s,wind (in the definition for HIF).
We show scrit for representative values ofΣ and R in Figure 5, over which we plot the mass-loss profiles from representative EUV (Hollenbach et al. 1994) and X-ray (Picogna et al. 2019) driven wind models assuming zIF = 4HIF. Typical grain sizes vary between 0.01 and 1 µm.
These values of scrit can be estimated analytically from the massloss rates, i.e. in the case of an EUV driven wind with density profile given by (23) (as is appropriate to disc radii interior to Rg = GM * /c 2 s,wind ≈ 5 au): Note that in the original model of Hollenbach et al. (1994), the density at the ionization front falls off more steeply with radius beyond Rg, scaling as R −2.5 , outside the gravitational radius: this effect has been included in the estimate forΣ(R) in the EUV case shown in Figure 5.
A simple estimate for scrit in X-ray driven winds may be esti- Figure 5. Maximum size of dust particles that can be delivered to a photoevaporative wind for a given mass-loss rate per unit area and radius (colour map and white contours). The black lines show the mass-loss profiles for the EUV driven wind model with Φ = 10 42 s −1 (Hollenbach et al. 1994) and X-ray driven wind with L X = 2 × 10 30 erg s −1 (Picogna et al. 2019), for which the integrated mass loss rates are 7.4 × 10 −10 and 2.7 × 10 −8 M yr −1 , respectively. Note that s crit is given here for the case ρ grain = 1 g cm −3 and an ionization front height of z IF = 4H and scales with H/z IF ρ grain (Equation 25). mated fromΣ ∼ 2 × 10 −12 (R/10 au) −3/2 g cm −2 s −1 for an Xray luminosity, LX = 2 × 10 30 erg s −1 (Picogna et al. 2019). The corresponding estimate for the maximum grain size entrained is then More precise numbers can obtained by directly using the fits forΣ provided by Picogna et al. (2019), which were used in Figure 5. Finally, the distribution of dust entrained in the wind can be computed from the flux efficiency, F (Figure 3). Since F falls off rapidly for St > Stcrit the contribution from sizes much beyond Stcrit can be neglected.

DISCUSSION
In this paper we have demonstrated that the removal of dust from protoplanetary discs by winds is driven by advection unless turbulence in the disc is strong, i.e. α Stcrit ≈ 0.01 (Equation 14). Note that Stcrit is measured in the disc immediately below the ionization front and thus varies only weakly with system parameters (see Figure 4), being mainly set by the difference in temperature between the disc and the ionized wind. We find that strong turbulence acts to increase the amount of dust supplied to the wind by a factor of a few in certain size ranges (see left hand panel of Figure 3), rather than decreasing it, as was found by Hutchison & Clarke (2020). This difference can be attributed to the way in which diffusion was treated in the two studies. Hutchison & Clarke (2020) treated diffusion by adding a diffusive flux to the mass-conservation equation using the model of Dubrulle et al. (1995). Here we have used a Monte-Carlo model for the dust in which diffusion is treated through the direct coupling of dust to turbulent motions in the disc gas via drag forces, modelling the turbulence assuming isotropic Gaussian turbulence with a constant velocity dispersion. The explanation for this difference is that if one simply adds a diffusive flux using the Dubrulle et al. (1995) model, it implies an increase in the dust to gas ratio across the front which can drive a strong negative diffusive flux. The reason why this does not happen in reality is that, as particles cross the ionization front, they decouple from the gas flow and do not participate in diffusive motions. This behaviour can only be captured by a treatment that explicitly models the ability of particles with finite stopping time to decouple from the diffusive motions over a region where there is a steep gradient in background gas properties. We therefore caution against applying the Dubrulle et al. (1995) model in situations involving ionization fronts.
Although isotropic turbulence is likely a poor approximation at the ionization front, our results are unlikely to be substantially affected by this. This is obviously the case when the delivery of dust is dominated by advection, which is the case for both weak turbulence and sufficiently small particles in the regime of strong turbulence. Since large particles cross the ionization front within a stopping time, they cannot couple to the gas within the ionization front, and therefore the details of the turbulence at the ionization are not important. For intermediate grain sizes in conditions of strong turbulence, particles begin to decouple within the ionization front. The mass-loss rate of these particles is sensitive to the width of the ionization front, and therefore possibly also sensitive to details of the turbulence there. However, the dependence of the flux on the properties of the ionization front is weak, and the phenomenological behaviour is unlikely to be affected.
Our results suggest that vertical advective transport in discs could play an important role in determining the vertical height of discs measured in scattered light. Recent non-ideal MHD simulations suggest that discs may have weak turbulence (α 10 −4 ), with observational studies providing supporting evidence (see, e.g. Mulders & Dominik 2012;Flaherty et al. 2015;Simon et al. 2018;Flaherty et al. 2020). Under such conditions our models show that advective transport due to the wind should dominate the lofting of small grains -this could be tested by comparing resolved observations of disc thickness in the scattered light to the thickness derived for millimetre grains (e.g. Pinte et al. 2016;Avenhaus et al. 2018;Villenave et al. 2020). If turbulence is stratified, i.e. α increases with height, then turbulence might still play an important role. However, Riols & Lesur (2018) found that the vertical variation of α could not explain the lofting of grains seen in their MHD simulations unless advective transport was also included. We suggest that this is likely to be a general feature of discs undergoing mass loss due to winds, independent of the winds' origin.
In our calculations we have neglected the influence of radiation pressure on the dust grains. Since the optical photosphere is at lower altitudes than the EUV (or X-ray) photosphere, Owen & Kollmeier (2019) argued that radiation pressure could remove grains efficiently. We now show that when including advective transport, radiation pressure does not greatly change the picture. Neglecting turbulence, but including radiation pressure, the vertical and radial velocities of dust grains are given by where β is the ratio of the radiation pressure force to the gravitational force. Here we have assumed that β is large enough that the radiation pressure term dominates over all other components (such as the radial gas pressure gradient) in the the equation for vR from Owen & Kollmeier (2019). Re-writing the gas velocity in terms of Stcrit, If the height of the ionization front scales as zIF ∼ R 1+δ , then particles with St Stcrit/(1 + βδ) will be delivered to the wind.
For typical values of β, the maximum grain size delivered to the wind is not much affected. Note that although radiation pressure can increase the maximum size of particles that can be entrained when δ < 0, this is not the case for our EUV model with Σ ∝ 1/R. Another factor that could reduce the size entrained would be if there was a steep dependence of scrit with radius, such that radiation pressure drives the grains to larger radii where they can no longer be entrained. However, given that the dependence of scrit is not particularly strong ( Figure 5), this will also not dramatically affect the maximum size of grains delivered to the winds.

CONCLUSIONS
We have investigated the entrainment of dust grains in photoevaporative winds using a novel Monte-Carlo dust dynamics model which correctly models dust transport across the ionization front separating the disc and ionized wind. This treatment avoids spurious effects previously found when solving the advection-diffusion equation in the limit that the width of the ionization front is less than the dust stopping distance. Our calculations yield dust transport efficiencies that converge in the limit of narrow ionization fronts as expected. We highlight that special care needs to be taken in the choice of Monte-Carlo dust modeling algorithm in the demanding case of a steep density feature such as an ionization front and that algorithms in the literature produce numerical artefacts under these conditions. Our simulations show that the delivery of dust to the wind base is dominated by the advection of small dust grains by the vertical gas flow that appears as a consequence of the photoevaporative mass loss. This is contrary to the usual assumption that turbulent diffusion is responsible for lofting grains to the ionization front, which we show only occurs if disc turbulence is strong (i.e. for values of the Shakura & Sunyaev 1973 α-parameter 0.01).
Our results confirm the hypothesis of Hutchison & Clarke (2020) that the maximum size of dust grains entering the wind is set by the condition of zero force on a stationary dust grain immediately below the ionization front (a limit that we denote as scrit). This is not the same as the commonly assumed limit (which, following Hutchison & Clarke 2020 we designate smax) which corresponds to the condition of zero force on a stationary dust grain immediately above the ionization front. The drag force scales as the product of the gas flux and the local sound speed: since the flux is conserved across the front, this means that the ratio of scrit to smax is given by the ratio of local sound speeds, and is thus typically around ∼ 0.1 in the case of ionized winds from protostellar discs. Equation (25) allows the value of scrit to be estimated for any wind where the local mass flux and height of the base of the heated region is known; Figure  5 illustrates the typical values that apply in the case of mass loss profiles for canonical EUV and X-ray driven winds. These values are lower by around a factor 10, for equivalent parameters, than those previously proposed (Takeuchi et al. 2005;Owen et al. 2011a;Franz et al. 2020), a result that we ascribe to the aforementioned difference between scrit and smax.   , 1983, Boundary-Layer Meteorology, 27, 163 Youdin A. N., Lithwick Y., 2007 APPENDIX A: EXTRA FLUX-EFFICIENCY PLOTS In Figure 6 and Figure 7 show the flux efficiency, F, for models with ionization front heights, zIF of 3H and 4H respectively. The results are nearly identical to case with zIF = 3.5H presented in section 5. This paper has been typeset from a T E X/L A T E X file prepared by the author.