© ESO 2020

1 Introduction

Planets form from the gas and dust surrounding newly born stars, whose physical properties and final dispersal are strongly influenced by the stellar irradiation from their host star. In particular, high-energy radiation may warm up the disk atmosphere, launching a thermal wind (see e.g. Hollenbach et al. 1994; Gorti & Hollenbach 2009; Alexander et al. 2014). Models predict that this photoevaporative wind can ultimately disperse the disk and may have important consequences for the formation and evolution of planetary systems (Alexander & Pascucci 2012; Ercolano & Rosotti 2015; Ercolano et al. 2017; Carrera et al. 2017; Jennings et al. 2018; Monsch et al. 2019).

Despite the potential influence of this process on the formation of planets, the magnitudes of photoevaporative winds are still largely uncertain, with model predictions diverging by several orders of magnitude (Armitage 2011; Alexander et al. 2014; Ercolano & Pascucci 2017). One problem is that to date, the only direct evidence of these winds is blue-shifted forbidden-line emission towards T Tauri stars, including [Ne II] 12.8 μm and [O I] 6300 Å (e.g. Hartigan et al. 1995; Pascucci et al. 2008; Rigliaco et al. 2013; Natta et al. 2014; Simon et al. 2016; Banzatti et al. 2019). While the intensity and low-resolution profiles of these lines can be matched very well by X-ray photoevaporation models (Alexander 2008; Ercolano & Owen 2010, 2016), it has been demonstrated that these lines do not trace the base of the wind. Furthermore, their extreme temperature dependence makes them a tracer of the heating mechanism of an already unbound wind rather than tracing the wind-driving mechanism itself (Ercolano & Owen 2016). An additional problem is that high resolution data has revealed very complex line profiles which may include components emitted in a magnetically-driven wind (Banzatti et al. 2019).

Different types of wind diagnostics would be desirable for constraining disk dispersal models; small dust grains – which can be entrained by the wind – may provide an interesting avenue towards this end (Giacalone et al. 2019). A previous work by Owen et al. (2011a) has shown that grains up to about 2 μm in size (i.e. radius) can be lifted up and blown out by an EUV-driven wind around a Herbig Ae/Be star. More recently, Hutchison et al. (2016a) have investigated EUV-driven dust outflow by means of a two-fluid smoothed particle hydrodynamics (SPH) code (Hutchison et al. 2016b); they find entrainment of grains of up to 4 μm around a 0.75 M_⊙ T Tauri star (although they note that this value may drop to about 1 μm due to grains settling towards the disk midplane). Both Owen et al. (2011a) and Hutchison et al. (2016a) show that the wind selectively entrains grains of different sizes from different radii. This results in a dust population which spatially varies in the wind, due to the topology of the gas streamlines which propagate almost radially outwards. At NIR wavelengths, this variable grain population produces a “wingnut” morphology which may already have been observed in the case of PDS 144N (Perrin et al. 2006). Yet Owen et al. (2011a) could not reproduce the color gradient of the observations, which show redder emission at larger heights above the disk; they suggest that this could be because the synthetic observations might be dominated by emission from the smallest grains entrained inthe flow. Grain growth in the underlying disk (see Testi et al. 2014), which they neglected in their calculations for simplicity, could reduce the population of small grains, and may hence provide a solution to this color problem. While it is unclear whether the observations of PDS 144N can be explained by dust entrainment in a photoevaporative wind, Owen et al. (2011a) have demonstrated that a significant amount of small grains – which dominate the opacity in the FUV – do populate disk winds, and hence play an important role in their chemistry.

In this work, we study the entrainment of dust grains in an X-ray driven wind around a T Tauri star. For this we use a particle approach (Picogna et al. 2018), bootstrapped onto a steady-statehydrodynamical simulation of a photoevaporating disk (Picogna et al. 2019). Our results aim to facilitate more detailed studies of the detectability of winds in scattered light, as well as wind opacity models. The latter should allow for more realistic chemical modelling of the gas in the wind, enabling us to search for new wind diagnostics.

This paper is organised as follows: We present the numerical setup of the gas disk and photoevaporative disk wind, and the dust grain evolution in Sect. 2. In Sect. 3, we take a detailed look at what we can extract from the dust grain trajectories we obtain. We discuss our findings in Sect. 4 and summarize them in Sect. 5.

2 Methods

The dynamics of dust grains in protoplanetary disks can be studied either by directly integrating the orbits of a large number of dust “super-particles”, which sample the local properties of the dust population, or by solving the collisionless Boltzmann equation forthe particle distribution function. For a population of very small (i.e. tightly coupled to the gas) dust particles, the Boltzmann equation can be reduced to the zero-pressure fluid equation (Cuzzi et al. 1993; Garaud et al. 2004); this “two-fluid” approach has already been used to study planet-disk interactions (e.g. Paardekooper & Mellema 2004, 2006; Zhu et al. 2012). However, it is limited to a single population of small particles as it cannot account for the full velocity distribution of the grains at a single location, and it is not able to capture strong density gradients.

In contrast, a particle approach as implemented by Picogna et al. (2018) has the notable advantage of following the evolution of solid particles with different physical properties, recovering the dust dynamics very well also in the limit where the grains are decoupled from the gas (Youdin & Johansen 2007; Miniati 2010; Bai & Stone 2010). This method has been successfully applied to the study of planet-disk interaction with both SPH and grid-based codes (Fouchet et al. 2007, 2010; Lyra et al. 2009; Ayliffe et al. 2012; Zhu et al. 2014), and to modelling the draining of dust grains from the inner region of a photoevaporating transition disk (Ercolano et al. 2017).

2.1 Gas disk with XEUV wind

The set-up of the hydrodynamical calculations for the gas disk has been described in detail in Picogna et al. (2019), so here we limit ourselves to summarising the basic parameters of the specific run employed in this work. We studied a protoplanetary disk of M_disk ≃ 0.01 M_* around a M_* = 0.7 M_⊙ T Tauri star. This star was set to emit X-ray and EUV radiation according to the emission line spectrum presented by Ercolano et al. (2008a, 2009); the X-ray luminosity of the star was L_X = 2 × 10³⁰ erg s⁻¹, which is close to the median of the X-ray luminosity distribution for this stellar mass (Preibisch et al. 2005).

The hydrodynamics (HD) simulations were performed via a modified version of the Pluto code. In this version, at each HD step, temperatures are updated according to the local ionization parameter¹ and column density to the central source (for further details, see Picogna et al. 2019). The temperature parametrization was obtained via detailed radiative-transfer calculations using the Mocassin code (Ercolano et al. 2003, 2005, 2008b). Within the Pluto code used for the gas evolution (Mignone et al. 2007), we employed a 2.5D Eulerian grid in spherical coordinates². To avoid any boundary effects, a large radial range of 0.33 ≤ r [AU] ≤ 1000 in 412 logarithmically-spaced steps was modelled, with 0.005 ≤ ϑ ≤ π∕2 in 320 uniform steps.

The v_r- and ϱ-profiles of the gas disk employed for the simulations are showcased in Fig. 1; we used only the inner ≈ 300 AU (352 cells) of the hydrodynamical grid in order to follow the dust evolution because the mass loss due to the photoevaporative wind becomes negligible at larger radii (see Picogna et al. 2019).

The (sonic) disk surface is defined as the plane where the gas velocities change from locally sub- to super-sonic. The base of the X-ray and extreme ultraviolet (XEUV) flow, by contrast, is given by the location of the largest temperature gradient (see e.g. Ercolano et al. 2009), and lies slightly below the sonic surface. It coincides with a strong drop in gas density (see Fig. 1).

Due to the high grid resolution employed and the disk having settled into a stable (quasi-equilibrium) state, the (numerically computed) base of the wind is quite smooth (as would be expected; see also Bai 2017); hence, we forgo additional artificial smoothing. In Fig. 1, we have added lines indicating one, 3.5, and 5.5 scale heights H, the latter two framing the base of the flow; here we use H = h ⋅ R, and H = c_s∕Ω_K, with c_s the local sound speed and ${\Omega }_{\text{K}}=\sqrt{G{M}_{*}/r^3}$ the Keplerian orbital velocity.

The flow base being located at a scale height of h ≃ 0.08 (0.14; 0.23) at R = 10 AU (100 AU; 300 AU) implies that the disk is still rather hot and puffed up and that the dust grains have to travel rather far above the midplane if they are to enter the wind region (see comparable simulations by Owen et al. 2012; Bai 2016).

Fig. 1 Density (green) and velocity map (yellow arrows: v) of the gas disk model from Picogna et al. (2019). The base of the XEUV flow (i.e. the location of the largest gradient in temperature) is shown in black, and coincides with a strong drop in density. vr is pointing radially outwards everywhere in the wind but very close to the sonic surface where v points outwards, but away from the disk (for R < 200 AU). The cyan line represents z(R) = H; additional lines for 3.5 H (dotted blue) and 5.5 H (dashed blue) show the range of scale heights the disk surface covers.

Fig. 2 Left: initial placement of the dust grains at Δt = 0 at the base of the photoevaporative flow (black), slightly below the sonic surface. The grains are colored according to their size and the color scale (scaling with $\sqrt{\left|v_r\right|}$) of the gas disk represents the extent of its local radial velocity (with − 2 ≲ vr [km s−1] ≲ 30). The gas map is mostly smooth in the region of interest, that is everywhere but close to the midplane at high R. Right: snapshot of the grain positions at Δt ≃ 100 yr, all else equal. The very low spread of the lines of individual grain sizes is due to the initial setup, placing particles directly on the disk surface without a spread in their initial velocities. For comments on this, see also Sect. 3 and Appendix B.

2.2 Dust grains

The dust grains were modelled as passive Lagrangian particles inserted in the steady-state gas solution, as originally implemented by Picogna et al. (2018), to whom we refer for the details of the implementation³. The motion of these particles is driven by their gravitational attraction towards the central star, the drag force from the surrounding gas, and turbulent diffusion below the disk surface as prescribed by Charnoz et al. (2011).

Ormel & Liu (2018) give a concise comparison of their stochastic equation of motion to the strong-coupling approximation of Charnoz et al. (2011). Although the former may be preferable for modelling grain motion within the disk, we are mainly interested in what happens once a grain enters the wind region where the gas density and thus turbulence are low; hence, just like Giacalone et al. (2019) proceed for their magnetohydrodynamic (MHD) wind model of dust motion, we do not optimize our model for the disk interior.

Above the disk surface, gas densities are too low to induce kicks (see Fig. 1), allowing us to neglect an otherwise necessary (Flock et al. 2017), more intricate distinction between MRI and VSI⁴.

2.2.1 Grain sizes

In their EUV-only simulations, both Owen et al. (2011a) and Hutchison et al. (2016a) find that grains around μm size are entrained. In the MRI computations of Miyake et al. (2016), the grain distribution considered is 0.1 ≤ a₀ [μm] ≤ 100, and the maximum entrainable grain size is found to decrease very steeply with the (cylindrical) midplane radius $R=\sqrt{x^2+y^2}.$ Furthermore, Giacalone et al. (2019) investigate 5 × 10⁻³ ≲ a₀ [μm] ≲ 5.

Therefore, we ran an initial set of simulations with grain sizes 10⁻³ ≤ a₀ [μm] ≤ 10², which established that for our model, the size barrier for wind blow-out lies between 5 and 15 μm. On the basis of this initial experiment, we restricted our size range to 0.01 ≤ a₀ [μm] ≤ 20, with steps of Δa₀ = 1 μm for 1 ≤ a₀ [μm] ≤ 15. We forwent a higher size resolution in favor of increasing spatial resolution, that is tracing more dust grains per each size. All different a₀ modelled are listed in Table 1 (and Fig. 2); per size, we simulated the trajectories of at least 5000 dust grains, yielding at least 25 grains per 1 AU of launching radius along the disk surface.

2.2.2 Internal grain density

Following Owen et al. (2011a) and in order to facilitate a direct comparison to their results, we assume a uniform internal density of ϱ_grain = 1 g cm⁻³ for the dust particles. This value is on the lower end of the 0.3 ≲ ϱ_grain [g cm⁻³] ≲ 6.2 interval established by Love et al. (1994), and agrees best with the values Joswiak et al. (2007) find for material of cometary origin (0.6≲ ϱ_grain [g cm⁻³] ≲ 1.7). Similar values are used in other works, too (e.g. Tamfal et al. 2018; Owen & Kollmeier 2019).

Other models employ somewhat different values for ϱ_grain. For instance, Li & Greenberg (1997) and Miyake et al. (2016) use the average value given by Love et al. (1994), ⟨ϱ_grain⟩≃ 2 g cm⁻³; Hutchison et al. (2016a) and Flock et al. (2017) employ ϱ_grain = 3 g cm⁻³, and Weingartner & Draine (2001) and Giacalone et al. (2019) use ϱ_grain = 3.5 g cm⁻³. These values are closer to the ones Joswiak et al. (2007) find for asteroidal material, which may have been heated slightly less than the cometary grains. Future on-site analysis of interplanetary and interstellar dust grains will provide further constraints on these intervals (e.g. Arai et al. 2018, Destiny+).

2.2.3 Initial positioning

We position our grains directly on the base of the flow which is located slightly below the disk surface. This allows us to study their trajectories from when they enter the wind-dominated region above the disk. Within 0.33 ≤ r [AU] ≤ 200, we use a random distribution uniform in r for the initial placement.

The left panel of Fig. 2 shows the initial grain positioning along the base of the flow (in black), with the dust grains colored according to their size.

Our model is intended to be combined with a vertical mixing prescription later on in order to extract a realistic dust density distribution in the wind. So we may very well, just as Hutchison et al. (2016a) note they did, model grain sizes that will not migrate far enough vertically to actually reach the wind region (see also Youdin & Lithwick 2007; Krijt & Ciesla 2016).

The dust grains were placed – and remain – well outside of the sublimation radius applicable for a M_* = 0.7 M_⊙ star (Giacalone et al. 2019). However, we note that the intense X-ray radiation from the young stellar object may destroy PAH-size grains in the disk atmosphere before they can be entrained in the wind (Siebenmorgen & Krügel 2010; Siebenmorgen & Heymann 2012); we did not include such events in our simulation.

2.2.4 Initial velocities

We initialised our dust particles to start from a quasi-equilibrium; so we set both v_r,0 and v_ϑ,0 to 0 because the local gas velocities are quite low anyways and will hence not cause a strong upwards motion ($\left|v_{\vartheta }\right|\lesssim 50$m s⁻¹ along the disk surface compared to $\left|v_r\right|\lesssim 200$m s⁻¹).

For v_φ,0, we assumed a Keplerian velocity of v_φ,0 = r Ω_K; since the starting positions are at z ≳ 3.5 H (see Fig. 1), we computed the Keplerian speed for the spherical radius r, and not for the midplane radius R.

2.2.5 Further limitations

To cut computational costs and allow for a reasonable amount of particles to be modelled, we made a series of simplifying assumptions:

Firstly, we neglected MHD effects⁵.

Secondly, we did not include self-gravity from the disk. We show in Appendix A that this simplification should not significantly affect our results.

Thirdly,we did not include dust-gas back reactions. As shown by Dipierro et al. (2018) and Tamfal et al. (2018), these are important in the disk midplane; but we focus our modelling efforts on the wind regions above the disk, where dust-to-gas ratios are not expected to be enhanced (Krijt & Ciesla 2016).

Fourthly, dust-dust interactions were neglected. The gas drag accelerates the dust grains to at least a few km s⁻¹ (i.e. $v_{r,\text{esc}}=\sqrt{2GM/r}\approx 11\text{\hspace{0.17em}km}{\text{s}}^{\text{-1}}$ at r = 10 AU, or v_r,esc ≈ 2 km s⁻¹ at r = 300 AU), but the dust densities in the wind are much lower than around the disk midplane. At the latter, the growth time scale is already around 10² –10³ yr (Birnstiel et al. 2016). We shall see below that therefore, assuming no interactions provides a reasonable approximation.

3 Results

We traced the trajectories of the dust grains over time until they either leave the computational domain or until the simulation time frame of Δt_sim ≃ 2.2 kyr ends.

If they left the domain, they were replaced by a new grain of the same size, placed as described in Sect. 2; the actual amount of dust particles modelled per a₀ is listed in Table 1⁶. Merely for visualizing the actual simulation, the right panel of Fig. 2 shows a snapshot of the simulation after around 100 yr.

A selection of trajectories obtained from the simulation is shown in Fig. 3; for clarity, we limit ourselves to plots for three distinguished grain sizes (0.1, 4, and 10 μm) below. These were chosen because they represent the three major varieties of grains encountered (see Appendix C). Panels containing the complete set of 20 different a₀ are included in Appendix C.

In general terms, the grains analyzed are either fully entrained (blown out by the XEUV wind, leaving the computational domain above the disk surface), fall back below the base of the flow at R ≳ 160 AU, or are not even picked up by the wind despite the turbulent kicks allowing for upwards motion. Additionally, we find that trajectories for a given grain size do (almost) never intersect, wherefore different starting positions will lead to different paths in the wind. Thus, the initial positioning of a grain of size a₀ along the launching region pre-determines which wind regions it can populate. This matches with – and is a direct result of – the gas wind velocity map seen in Fig. 1 pointing radially outwards almost everywhere.

3.1 Robustness of the initial velocity setup

In Sect. 2, we have described a rather simplistic initial velocity setup, with v₀ depending only on Ω_K. However, in order for our grains to even reach the 3.5...5.5 H which we launch them from (see Fig. 1), the gas must have some degree of turbulence, implying some variety in the initial velocities.

To account for this, we have run a series of tests with a Gaussian spread of σ(v_i) = 100 m s⁻¹ in all three directions i ∈{r;ϑ;φ} of the initial velocity vector v₀. This value was chosen because it is an overestimate of the fragmentation speeds given by Birnstiel et al. (2009, 10 m s⁻¹) and Wada et al. (2013, ≲ 8 m s⁻¹ for silicates, ≲ 80 m s⁻¹ for icy aggregates), and because it is slightly higher than the upwards speed of the gas which we find along the base of the flow, |v_ϑ | ≲ 50 m s⁻¹. So it should make for a suitable approximation of a velocity spread introduced by turbulent vertical mixing.

The results obtained do not deviate significantly from those retrieved without this spread; therefore we will proceed to show only the latter. For a more extensive elaboration on the similarities and differences identified, see Appendix B.

Fig. 3 Randomly selected dust trajectories for a0 = 0.1, 4, and10 μm (left, middle, and right panels, respectively). The trajectory color represents the local value of the Stokes number. Entrained dust grains, launching from the launching region (black), migrate upwards on the colorbar the as they move to regions of lower gas density. Gas streamlines are shown in dash-dotted grey. (For all simulated grain sizes, see Fig. C.1.)

3.2 Dust coupling to the gas

In Fig. 3, the trajectories are colored by their Stokes number St = t_stop ⋅ Ω_K, with t_stop = m_dust v_dust∕F_drag their local stopping time⁷. Besides, gas streamlines spaced by 5% of the total mass-loss rate of the gas in the wind region (i.e. Ṁ_w) have been added in (dash-dotted) grey for direct comparison.

St ≪ 1 indicates that the dust motion is well-coupled to the gas flow; hence St ≪ 1 is needed for a particle to be lifted up by the wind, since it is only affected by gas drag and stellar gravity. Indeed, we find all entrained grains to have St ≲ 0.4 when they are picked up by the wind (i.e. colors from blue to green). While in the wind, they are sped up by the rather fast photoevaporative flow (of up to v_r ≲ 30 km s⁻¹, see Fig. 2), which leads to a steady increase in their speed, and thus also St; the latter may grow by up to an order of magnitude.

At low St, the dust grains follow the gas flow (see especially the left panel of Fig. 3); at St → 1 however, they decouple from the gas flow (see the middle and right panels of Fig. 3, especially for higher R). For R ≳ 160 AU, the gas streamlines – in particular those close to the disk – start to bend towards it because the stellar irradiation is starting to decline this far out. As a result, the dust grains that have already reached a relatively high radial velocity at this point overshoot the gas streamlines. Further inwards and at higher z, the dust trajectories fall below the gas streamlines if they become decoupled.

For an in-depth analysis of the dust motion, Fig. 4 shows two randomly selected dust particles launched from R ≃ 20 AU. These are representative of the dust grains picked up by the wind from this R, with other grains entrained from around this R showing very similar trajectories; we opted for a launching point rather close to the star in order to showcase fully-entrained grains⁸.

The 0.1 μm particle (left column) remains entrained in the photoevaporative flow and follows the gas motion almost perfectly; its Stokes number remains small (St < 0.1) throughout its trajectory. Its t_stop also stays small, even after the strong increase (of a factor of about 4) it experiences when being picked up by the wind, simultaneous to the strong decrease in the density of the surrounding gas (of a factor of almost 100). While the grain is at r < 300 AU, t_stop is always smaller thanthe time needed for blow-out to 300 AU (dashed grey line in the second panel), which may serve as a further indication that the particle stays coupled to the gas. The fourth through sixth panels show a comparison of gas (blue) and dust (red) velocities for v_R and v_z, where we can also observe a strong coupling. It is only at r → 300 AU that a slight deviation of v_dust from v_gas occurs; the decoupling does not necessarily coincide with the particles reaching escape velocities (dashed grey lines in the fourth panel). This means that even sub-micron particles start decoupling from the gas flow at high r, that is after picking up enough momentum from the wind. Also, with the curve for ϑ flattening down and v_ϑ being rather small in comparison to v_R, the motion of the showcased particle is almost fully radial at larger R. As v_ϑ < 0, a small additional upwards component remains.

On the other hand, the 10 μm particle (right column) decouples from the gas flow within the first few AU of entering the wind region; simultaneously, its Stokes number quickly becomes St > 0.5. Its t_stop grows even larger than the time passing between wind pick-up and leaving the domain (again, dashed blue line). The XEUV wind drags the grain along, increasing its v_R which remains comparable to v_R,gas for the first part of its trajectory. By contrast, v_z quickly diverges from the gas flow, which leads to the direction of the dust motion clearly (visually) differing from the gas streamlines, intersecting multiple ones. At R ≳ 160 AU, the gas flow starts pointing back down towards the disk, but since the grain is already decoupled, it does not seem to be affected by this; in the sixth panel, we see that its direction of motion remains almost constant after the initial acceleration, as would be expected for high t_stop.

Despite the gravitational pull acting on the 10 μm grain, its v_R and v_z are slightly increasing for R ≳ 100 AU. This is caused by the relatively high difference in gas and dust velocity $\left|v_{\text{gas}}-v_{\text{dust}}\right|$; even if this additional speed-up were missing, the grain would still reach v_r,esc well within r < 300 AU, as the grey dashed line in the fifth panel shows.

Fig. 4 Analysis of two dust trajectories (left: a0 = 0.1 μm, right: a0 = 10 μm) entrained in the photoevaporative outflow from R ≃ 20 AU. Top panels: actual motion in the (R, z)-plane colored by St. Lower panels, from top to bottom: tstop and ϱgas, St and ∂t St, the horizontal speed vR of gas and dust, the vertical speed vz of gas and dust, a comparison of the direction of the motion vz∕vx for gas and dust, and ϑ and vϑ. See the text for an in-depth commentary.

3.3 Dust timescales in the wind

In Fig. 4, we have included a timescale Δt_bnd for the motion of the dust grains. This timescale was computed as the difference between the time at which the particle is picked up by the wind (t_wind,0) and the time at which it crosses the domain at r > 300 AU while being entrained (t_bnd), thus Δt_bnd = t_bnd − t_wind,0. The full distribution of timescale data points for all trajectories is shown in Fig. 5 (cyan and blue).

In general, wind entrainment timescales appear to span 10² –10³ yr. This is similar to the range given by Birnstiel et al. (2016) for midplane dust growth; since dust densities at z ≳ 3.5 H (see Fig. 1) are considerably lower, dust-dust interactions should indeed be negligible, as claimed in Sect. 2. More concisely, Kornet et al. (2001) have established that the dust growth timescale goes as t_grow ∝ ϱ_dust ⋅ v_dust, with ϱ_dust the local dust density⁹. As Table 2 shows, ϱ_gas drops off heavily towards the disk surface while $v_{\text{gas}}\equiv \left|v_{\text{gas}}\right|$ remains very comparable; Figs. 1 and 2 illustrate that the wind speed picks up slightly above – and not at – the base of the flow (as would be expected considering the latter is located slightly below the disk surface). Hence the dust-dust interaction timescale is much longer than the wind blowout time, if we assume a constant dust-to-gas ratio. In case of a more realistic relation, dust would be even more scarce than gas for similar z (Krijt & Ciesla 2016).

The 0.1 μm grains of Fig. 5 have Δt_bnd as low as 70 yr if they are picked up at small R; this means that the longer distance to the domain boundary is outweighed by the higher acceleration the dust experiences close to the star. For R ≳ 80 AU, Δt_bnd decreases; if we assume the general trend of slower speed-up from larger R to persist, this would mean that the grains are picked up close enough to the computational boundary (at R ≃ 300 AU) to be blown out faster than those from slightly further in. So this drop-off is caused by the numerical setup, not by the actual physics involved.

For a₀ ≥ 0.5 μm, we do not see this fall-off anymore. Disregarding the various local peaks in Δt_bnd which are caused by the base of the wind notbeing perfectly smooth¹⁰, a clear trendof the blow-out time Δt increasing with the launching position R emerges. For10 μm grains, we find max (Δt) ≈ 10³ yr, which is still well below the simulation time frame Δt_sim, validating the latter a posteriori.

Since the cutoff at r ≃ 300 AU is somewhat arbitrary, the orange and red parts of Fig. 5 show a different approach to defining a timescale: these represent the time between wind pick-up (t_wind,0, as above) and reaching v_r,esc at t_esc, that is Δt_esc = t_esc − t_wind,0. Because the velocity field of the gas flow in the wind is pointing outwards (see Fig. 1), it is highly unlikely that a grain will not be fully blown out by the XEUV wind once it has reached v_r,esc.

The overall appearance of the average values for ⟨Δt_bnd⟩ (blue) and ⟨Δt_esc⟩ (red) is quite similar; both show a rather distinct upwards trend, mitigated only for R ≳ 80 AU and a₀ < 0.5 μm. Interestingly, this feature also holds for Δt_esc, which in contrast to Δt_bnd does not depend on the choice of the computational boundary.

The values retrieved for Δt_esc are different – and always smaller. This difference is most pronounced for small grains, which reach v_r,esc within 5–60 yr (for a₀ = 0.1 μm), while needing 70–170 yr to leave the simulation domain. Dust particles with a₀ ≤ 8 μm launched close to the star are accelerated to v_r,esc within merely a few years; bigger grains take longer to pick up speed, or are too heavy to be picked up by the wind at all (a₀ ≥ 12 μm).

The dashed green lines in Fig. 5 indicate the steady-state dynamical timescale $t_{\text{dyn}}=2\pi /{\Omega }_{\text{K}}=2\pi /\sqrt{G{M}_{*}/r^3}$ for a Keplerian orbit at the base of the XEUV-driven flow. For a₀ = 0.1 μm, t_dyn dominates the dust motion only for starting points R ≲ 20 AU; this increases to R ≲ 30 AU for a₀ = 4 μm. By contrast, the photoevaporation of the 10 μm grains is largely dominated by t_dyn, which for those is mostly larger than Δt_bnd.

Because Δt_esc < Δt_bnd, we find t_dyn > Δt_esc for both 0.1 μm and 4 μm for R > 10 AU, meaning that for almost all grains the blow-out happens (considerably) faster than their “usual” timescale. It is only for the large particles (a₀ = 10 μm) that the time needed for acceleration to the escape velocity becomes comparable to the Keplerian timescale.

Fig. 5 Time needed to fully blow out dust particles from when they first enter the wind at twind,0 to the domain boundary, which they reach at tbnd (data points in cyan, mean in blue), and time needed to accelerate dust grains to vr,esc (data points in orange, mean in red); both for a selection of three a0. Keplerian orbital times (tdyn ≡ ΔtK) at the disk surface are included as dashed green lines. The rasterization of the data points results from a time-discrete particle tracking and a binning in R-direction; thus, one raster point may represent multiple data points. (For all a0, see Fig. C.2.)

3.4 Maximum entrained grain size

As has already been suggested in Fig. 2, dust grains may be too heavy to be blown out by the XEUV wind. This is investigated further in Fig. 6 where the blue line shows the maximum grain size ${\text{max}(a_0)|}_R$ that can be fully entrained in the wind from a starting position R along the base of the flow. As noted in Sect. 2, the latter is not entirely smooth; this, in turn, causes the craggy appearance of the graph for max (a₀).

To facilitate direct comparisons, we have included the corresponding EUV-only curve of Hutchison et al. (2016a) in orange and the MHD-wind results of Miyake et al. (2016) in green. The figure shows that overall, an X-ray driven wind is able to entrain larger grains over a larger radial range than its EUV-only or MHD-driven counterparts; this may have a noticeable impact on deduced opacity maps, and also on the detectability of the wind in scattered light.

The biggest grains that are blown out, that is a₀ = 11 μm, are entrained from 15 ≲ R [AU] ≲ 30. Closer to the star, stellar gravity counteracts the gas drag force; but since the gravitational pull drops off with r², the drag force dominates particle motion further out¹¹. With increasing R, max (a₀) slowly decreases out to R_max ≃ 140 AU, where it quickly drops to max (a₀) = 0. This coincides with the maximum radius at which XEUV photoevaporation is effective for the gas component of the disk; Picogna et al. (2019) show that the surface mass-loss rate (${\dot{\Sigma }}_{\text{gas}}$) drops to negligible values at R ≈ 140 AU. So both gas and dust residing at the disk surface at R ≳ 140 AU are very unlikely to be thermally unbound from there. Thus, for both their and our simulations this marks the outer boundary of the XEUV-wind-dominated region of the protoplanetary disk. Furthermore, it validates a posteriori our choice to limit the computational domain to r ≲ 300 AU and the initial particle placement to r ≤ 200 AU¹².

As mentioned in Sect. 2, ϱ_grain is not well-constrained. A variation of the dust density is found to directly anti-correlate with max (a₀), that is ϱ_grain ∝ 1∕max (a₀); this agrees with the analytical findings of Hutchison et al. (2016a, their Eq. (15)). Hence, for a threefold internal grain density of ${\varrho }_{{\text{grain}}^{\prime }}=3\text{\hspace{0.17em}g}{\text{cm}}^{-3}=3{\varrho }_{\text{grain}}$, our global max (a₀) drops to $\text{max}{(a_0)}^{\prime }=\text{max}(a_0)/3\approx 3.5\mu$m. An accordingly scaled version of our results is included in Fig. 6 as the blue dashed line in order to allow for simpler comparison to the results of Hutchison et al. (2016a,b); but it should be kept in mind that their and our mass-loss rates are not entirely similar. While they quote a surface mass-loss rate of ${\dot{\Sigma }}_{\text{gas}}=3.6\times {10}^{-12}(\text{g}{\text{cm}}^{\text{-2}}\text{\hspace{0.17em}s})$ at R = 5 AU, we have ${\dot{\Sigma }}_{\text{gas}}\approx 2\times {10}^{-13}(\text{g}{\text{cm}}^{\text{-2}}\text{\hspace{0.17em}s})$ (Picogna et al. 2019, their Fig. 5).

A 2D map of the grain sizes that can populate different regions of the wind is shown in Fig. 7¹³. The larger grains remain rather close to the disk surface; by contrast, smaller ones are lifted up to larger scale heights. For the same launching position R, smaller grains reach higher z. The abrupt decline of the 20 μm grains at R ≳ 160 AU results from the initial placement of the particles within 0.33 ≤ r [AU] ≤ 200, or equivalently 0.3 ≲ R [AU] ≲ 160.

In Fig. 8, we see the regions which are populated by grains of a given a₀ in green. We find that there are no deserts of smaller dust particles in regions populated by larger ones; for instance, wherever we find grains with a₀ = 10 μm, we also find grains with a₀ = 0.1 μm and a₀ = 4 μm. From Fig. 7, we have learnt that smaller grains will reach higher scale heights; to quantify this behaviour, we have included the maximum height ${\text{max}(z)|}_R$ for a certaina₀ at R in orange in Fig. 8. Fits with a simple second-order polynomial, $${\text{max}(z)|}_R=c_{1}R+c_2{R}^2,$$(1)

are shown asblue dotted lines; they match ${\text{max}(z)|}_R$ quite well except for the very inner region where the dust distribution is slightly more flared¹⁴. Higher-order polynomials fit the inner region better, however they are not included here since these fits are not a physically-derived, but merely a numerical prescription which we intended to keep quite simple.

The fit parameters c₁ and c₂ for all grain sizes are shown in Fig. 9 together with merely phenomenological prescriptions for the scaling of c₁ and c₂ with a₀ which are intended for comparisons to observational data from edge-on disks in future work. The fit formulas used are meant solely for a simplified reproduction of the values of c₁ and c₂; they are not based on physical considerations.

As may have been expected, both c₁ and c₂ decrease with increasing a₀, indicating the decline of the slope of ${\text{max}(z)|}_R$ with a₀ already seen in Fig. 8. The higher errors at lower a₀ result from the larger inclination of ${\text{max}(z)|}_R$ for smaller grains; in addition, the population map shown in Fig. 8 has a finite resolution as outlined above, which further contributes to a larger uncertainty especially for very steep lines.

Furthermore, we find c₁ ≫ c₂ for all a₀ modelled; this corresponds to the mostly linear appearance of ${\text{max}(z)|}_R$ especially for larger R. Yet c₂ does not drop to zero, so a certain amount of flaring is preserved for all particle sizes.

Fig. 6 Size of the largest grains entrained from a point R along the base of the wind (blue, peak at about 20 AU); the saw-tooth appearance of the curve at larger R is caused by the finite resolution of the underlying gas grid. When comparing to Hutchison et al. (2016a, their Fig. 7 for M* = 0.75 M⊙) (orange, peak at around 40 AU), we can see the size enhancement – especially at smaller R – caused by the inclusion of X-rays in our photoevaporative wind model. The blue dashed line represents our results scaled down by a factor of 3, to compensate for the differing internal grain densities of Hutchison et al. (2016a) and this work; yet, Ṁw still differsbetween the models, making a direct comparison difficult. The MHD wind model investigated by Miyake et al. (2016, their Fig. 4) (green) shows a distinctly different entrainment curve, starting off at very high a0 in the jet region but dropping towards max (a0) = 0 very quickly.Around R ≃ 140 AU, our max (a0) plummets to 0.

Fig. 7 Maximum size max (a0) of thedust grains in the wind in analogy to Owen et al. (2011a, their Fig. 2). The brightest color found above the sonic surface indicates a global maximum of max (a0) = 11 μm; the region below the base of the wind (black) is included merely for completeness. The visible correlation between grain size and maximum scale height is further investigated in Fig. 8.

Fig. 8 Areas populated by the dust grains (green); since we do not start from a realistic distribution along the base of the wind, we do not portray a density map. Smaller grains reach higher z at similar R. Wind base in black, numerical ${\text{max}(z)|}_R$ in orange, and corresponding fit in (dotted) blue; fits according to Eq. (1) annotated. (For all a0, see Fig. C.3.)

Fig. 9 Fit parameters for all fits for ${\text{max}(a_0)|}_R$ as a functionof a0 in blue, with 3σ-errors included as shaded regions. Non-physical fits to these fit parameter curves in orange, with the parameterizations given in the text boxes.

4 Discussion

We have numerically simulated the trajectories of dust grains in the XEUV-irradiated wind regions of a gaseous protoplanetary disk. As was to be expected (see e.g. Armitage 2015), we found small dust particles to show very good agreement with the gas streamlines and, in contrast, bigger grains to noticeably deviate from them. Thus, analytical models of gas motion as provided by Clarke & Alexander (2016) cannot be used to accurately model the trajectories of dust grains attaining St ≳ 0.5 in the wind.

4.1 Photoevaporative winds and radiation pressure

Owen & Kollmeier (2019) use the X-ray photoevaporation model of Owen et al. (2011b) to quantify grain blow-out due to direct radiation pressure alone; since their M_*, M_disk, and L_X are the same as used in this work and their internal grain density of 1.25 g cm⁻³ only slightly differs from our ϱ_grain = 1 g cm⁻³, this allows for an almost perfect comparison of their results to grain entrainment by XEUV winds¹⁵. Using an effective stellar surface temperature of 4500 K, they find max_rad(a₀) ≃ 0.6 μm for the largest grains for which the radiation pressure from an equivalent black body still outweighs stellar gravity. Thus, a photoevaporative XEUV wind enhances the size of dust grains blown out by a factor of almost 20 over radiation pressure alone.

Hence we can conclude that a photoevaporative XEUV-driven wind is much more effective at removing larger dust grains than direct radiation pressure.

4.2 Comparison to MHD wind models

Miyake et al. (2016), building on an MHD model established by Suzuki & Inutsuka (2009), performed 1D simulations of dust motion in MHD-driven winds. While they mainly focused on floating grains¹⁶, they also give a maximum entrainable grain size for their model along R; we have included it in Fig. 6 in orange. This shows that whereas MHD winds excel at removing large dust grains from regions very close to the star, photoevaporative winds start to dominate – in terms of entrained grain size – at R ≳ 2 AU. Thus, MHD winds would seem to be limited to the jet region and its immediate surroundings.

Recently, Giacalone et al. (2019) have investigated dust transport in (cold) magneto-centrifugally-driven disk winds in 2D; they conclude that the region of interest for dust pick-up from and re-deposition on the disk surface covers their full modelled range of R, which they chose to set up as 0.1 ≤ R [AU] ≤ 100. For their T Tauri model, they find entrainment of grains with a₀ ≲ 2 μm (their Fig. 3). This value was obtained for a disk surface temperature of T_surf = 600 K, which seems very plausible from an MHD point of view, but somewhat low from a photoevaporative one. They show as well that max (a₀) clearly depends on said T_surf, with higher T_surf increasing their max (a₀).

In contrast to our numerical setup, Giacalone et al. (2019) opted for a semi-analytical approach to trajectory modelling (see their Eqs. (1) through (3)); Fig. 10 compares the grain velocities we extract for the 10 μm grain of Fig. 4 to their prescription. We are only showing a comparison for a large dust particle here because for small a₀, gas and dust velocities are very similar. It should be noted, in this regard, that Giacalone et al. (2019) use their semi-analytical equations just for relatively small dust grains, considering their max (a₀) as quoted above; we present Fig. 10 merely to show that numerical simulations are still necessary for modelling the trajectories of the larger grains whose velocity distinctly decouples from the gas flow. Thus, until more intricate analytical prescriptions are introduced for tracing dust motion in a photoevaporative wind, a numerical approach must be employed for accurate results.

Fig. 10 Re-plot of vR and vz for the same dust grain as in the right column of Fig. 4 (a0 = 10 μm). The dust parameters are colored red, the surrounding gas in blue; the motion according to Eqs. (1) and (3) of Giacalone et al. (2019), with tstop as retrieved from our model, is added in green. Especially towards larger R, the semi-analytical prescription does not reproduce the dust motion very well. It should be noted, however, that Giacalone et al. (2019) do not use it for grains this big.

4.3 Comparison to EUV wind models

Previous investigations of dust entrainment in photoevaporative winds have limited themselves to EUV-only photoevaporation (Owen et al. 2011a; Hutchison et al. 2016a).

Owen et al. (2011a) base their work on an EUV-luminosity optimized model of a 2.5 M_⊙ Herbig Ae/Be star; while this means that our results are not directly comparable in terms of stellar parameters, both their and our work test the respective highest-luminosity (i.e. best-case) scenario, and hence aim to provide an upper limit to max (a₀). Comparing their max (a₀) ≃ 2.2 μm (see their Fig. 2) to our value of 11 μm, both retrieved for ϱ_grain = 1 g cm⁻³, we find a clear enhancement of particle sizes blown out by the wind. Thus, depending on the mass-loss rates caused by XEUV photoevaporation, the inclusion of X-ray photons in the disk irradiation model may be a crucial component for accurately predicting dust entrainment.

The basic mechanics do not change, though; when comparing our Fig. 7 to Owen et al. (2011a, their Fig. 2), we see the same qualitative behaviour. In both cases, the maximum height which a grain of a given a₀ at a given R can be lifted to decreases with increasing a₀. As expected, the largest entrainable grains are not lifted up very high above the disk by the photoevaporative wind.

Hutchison et al. (2016a) explore a wide range of stellar parameters; apart from choosing a penetration depth typical for T Tauri stars (provided by Woitke et al. 2016), they also ran their models for a stellar mass of 0.75 M_⊙, very similar to our M_* = 0.7 M_⊙. Just like us, they find an anti-correlation M_*∝ 1∕max (a₀) which explains their value of max (a₀) ≃ 4 μm in comparison to Owen et al. (2011a).

This value is still distinctly smaller than the 11 μm we have found to be picked up; yet as already noted in Sect. 3, we get very close to their EUV-only results with our model when using their ϱ_grain = 3 g cm⁻³. However, as pointed out above as well, the entrainable grain sizes are still not comparable since gas mass-loss rates strongly differ between Hutchison et al. (2016a,b) and our model. This is due merely to the different numerical setup of their and our models; Owen et al. (2012) and Owen & Jackson (2012) show that for viable L_X, X-ray driven winds clearly dominate over EUV-driven ones in terms of Ṁ_w. In addition to this, we have now seen that even at lower ${\dot{\Sigma }}_{\text{gas}}$, XEUV winds may entrain larger grains than EUV-only ones.

Furthermore, Fig. 6 shows that independent of the exact grain size, X-ray irradiation shifts the peak in max (a₀) towards lower R. The T Tauri star of Hutchison et al. (2016a) entrains its max (a₀) from 40 ≲ R [AU] ≲ 50, and the more luminous HerbigAe/Be star of Owen et al. (2011a) from 30 ≲ R [AU] ≲ 40. Interestingly, albeit exerting a stronger gravitational pull, the higher-mass star blows out its largest grains from further in; thus it seems that the higher gravity is outweighed by an even stronger wind launched.

In comparison, our XEUV wind picks up its largest grains from 20 ≲ R [AU] ≲ 30, indicating that dust grains are entrained more efficiently closer to the star. In contrast to Giacalone et al. (2019) and their MHD wind model, we only find grain fallback for R > 200 AU (a₀ = 0.1 μm), or R ≳ 150 AU (a₀ = 10 μm; see Fig. 8) for launching positions within a similar range of R – as noted above, the trajectories do not intersect (see Fig. 3). Hence, the photoevaporative wind re-deposits only little material on the disk surface. This agrees with the conclusions of Owen et al. (2011a) that once entrained, a dust grain will almost always remain in the wind, and be carried out to large radii. Furthermore, as noted above, all grains leaving the computational domain above the disk surface have v_r > v_r,esc, and are hence very likely to entirely leave the protostellar environment.

5 Summary

We have modelled dust trajectories for grain sizes 10⁻³ ≤ a₀ [μm] ≤ 10² in the wind region of a M_disk ≃ 10⁻² M_* gas disk around a M_* = 0.7 M_⊙ T Tauri star, irradiating its surroundings with L_X = 2 × 10³⁰ erg s⁻¹ on top of anEUV spectrum. Our main findings are as follows:

• X-ray driven winds are able to entrain grains up to a size of a₀ ≲ 11 μm; this is larger than the maximum entrained grain size from EUV-only models.

• XEUV winds pick up the largest particles from R ≃ 20 AU. By contrast, EUV-only winds entrain their largest dust from further out (i.e. R ≃ 40 AU). MHD winds show a very different profile, picking up very large grains from regions very close to the star (R ≪ 3 AU), but then rapidly loosing momentum farther from the star.

• Dust grains are launched with Stokes numbers St < 0.4 (i.e. St ≪ 1). Once entrained, the large grains decouple from the gas flow; smaller dust particles decouple at later times.

• μm-sized dust grains of are blown out of the inner 300 AU of a protoplanetary disk on a timescale of 10²–10³ yr.

• For a given grain size, the launching point of a grain determines its further trajectory.

• Smaller dust grains may be lifted up higher by the wind, with the maximum height max (z) at a given R decreasing with grain size a₀.

• An anti-correlation between ${\text{max}(z)|}_R$ and a₀ may be a typical signature of dusty photoevaporative winds.

Much like Owen et al. (2011a), we have found a strong dependence of the max (a₀) in the wind on R (see Figs. 7 and 8). This signature structure may be detectable in observations of edge-on disks if vertical mixing is strong enough to transport large grains to the disk surface; however, as has been shown by Hutchison et al. (2016a), this may not be the case, or may need additional MHD effects (Miyake et al. 2016).

In a future work we aim to present dust opacity maps and synthetic observations of typical protoplanetary dusty XEUV winds in order to investigate their detectability with current and future instrumentation. As the mass-loss profiles of winds due to the different mechanisms – that is, EUV, X-ray, MHD, etc. – are shown (or expected) to be different, this will reflect in the size distribution and density of dust particles entrained at different locations in the wind. While a quantitative discussion requires a calculation of detailed emission maps via radiative transfer modelling, we can speculate from the wind profiles of the gas that an X-ray-driven wind might produce a more extended launching region for larger grains than an EUV-only wind for which most of the entrainment is expected to occur near the gravitational radius of the disk. The dust emission is expected to be concentrated even closer to the star in the case of an MHD wind. Whether these differences are detectable with current instrumentation remains a matter of future investigation.

Acknowledgements

We would like to thank Cathie Clarke and Mark Hutchison for very helpful discussions, and the (anonymous) referee for a constructive report that improved the manuscript. This research was supported by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft), grants FOR 2634/1, ER 685/8-1, and ER 685/9-1, the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence Origin and Structure of the Universe. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769. The simulations have been carried out on the computing facilities of the Computational Center for Particle- and Astrophysics (C2PAP).

Appendix A Omitting disk gravity

As noted in Sect. 2, we decided not to include the gravitational pull from the gas disk in our model in order to cut computational costs. We verified a posteriori that this does not strongly impact our results.

Fig. A.1 Maximum of the per-a0 mean values for |Fdisk,z∕Fdrag,z| (see colorbar), mapped to a 2 AU × 2 AU grid. The gas drag clearly dominates over the (overestimated) disk gravity, especially in the wind region (base of the wind in black). Empty (white) regions above the sonic surface indicate that no particle has been recorded while in this cell, owing to the high grain speeds there (in contrast to Sect. 3, we did not interpolate between recorded particle positions).

To this end we compared, at each recorded position of all of our dust particles, the gas drag force $$F_{\text{drag},z}=\frac{\Delta v_{\text{gas},z}m}{t_{\text{stop}}},$$(A.1)

with a grain mass of $m=\frac{4\pi }{3}{a}_0^3{\varrho }_{\text{grain}}$, to an overestimate of the gravitational pull the disk would produce, $$F_{\text{disk},z}=-\frac{G{M}_{\text{disk}}m}{z^2}.$$(A.2)

The latter formula leads to a slight overestimate of the z-component of F_disk because the disk is not centered exactly below a given particle. If F_disk,z ≪ F_drag,z, or rather F_disk,z∕F_drag,z ≪ 1, neglecting disk gravity should not noticeably affect our results. To determine whether this is the case, we compute the mean value $\langle \left|F_{\text{disk},z}/F_{\text{drag},z}\right|\rangle$ of all grains of size a₀ on a 2 AU × 2 AU grid. The maximum values (over all a₀) for all these means are shown in Fig. A.1. As the colors demonstrate, the gas drag clearly dominates over the overestimated gravitational pull from the protoplanetary disk. We find only a few exceptions: firstly, well within the disk, from where it would be rather unlikely to see wind entrainment (as noted in Sect. 2, we have concentrated on realistically modelling the wind region, not the disk interior); secondly, very close to the host star, that is at R ≪ 10 AU, from where we do not see substantial entrainment (see Fig. 6); and thirdly, around the base of the flow for R ≳ 180 AU, which lies beyond the region from which wind pick-up happens in the first place, as noted in Sect. 3.

In other words, we find that the photoevaporative wind is – at low r – strong enough to compensate for both the closeness of the grain launching area to the disk midplane (and hence the center of mass) and also the concentration of the disk mass within r ≲ 100 AU (see Fig. 1).

At high r, the wind will not lose momentum (see the green velocity maps in Fig. 2), but the center of mass of the disk will have moved farther away from the base of the photoevaporative flow, meaning that gravity is even less likely to play a major role for wind entrainment.

Appendix B Including a spread in the initial velocities

As outlined in Sects. 2 and 3, our model was set up with no initial velocity spread, that is for all our grains, we set $v_0=r{\Omega }_{\text{K}}{\widehat{v}}_{\phi }$. Fig. B.1 shows that even large entrained grains reach v_r,esc well inside the computational domain (i.e. at r ≪ 300 AU); so we may assume that a variation of the starting velocities will not strongly affect max (a₀).

However, some turbulence is needed to vertically transport the dust particles to base of the wind; hence, we would realistically assume some spread in v₀. The main contribution to vertical mixing stems from v_ϑ which we found to be rather low in the vicinity of the base of the wind (|v_ϑ |≲ 50 m s⁻¹). So in order to check the effects of varying starting velocities, we introduced a Gaussian spread of σ(v_i) = 100 m s⁻¹ in all three coordinate directions i ∈{r, ϑ, φ}.

As a first step, we compared entrainment ratios from along the base of the flow. We expect them to differ due to grains with a reduced (enhanced) upwards velocity being less (more) likely to be picked up by the photoevaporative wind. In Fig. B.2, we see a comparison of said entrainment ratios for σ(v) = 0 (labelled η₀) and σ(v) = 100 m s⁻¹ (labelled η_σ(v)).

While we encounter some statistical variation between the two, there is no clear systematic distinction. Grains with a₀ ≲ 5 μm show little variation, with 0.95≲ η₀∕η_σ(v) ≲ 1.2; so overall, they are slightly less likely to be picked up if σ(v) > 0. This may seem counterintuitive at first; it is probably a consequence of the combination of three individual directional, randomly positive or negative offsets to v₀. A reduced speed in either of the three coordinate directions thus may be difficult to compensate via possibly positive changes along the other two axes of motion.

Dips to 0 indicate R-bins where η₀ = 0 and η_σ(v) > 0; in order to avoid a series of minima produced by only very few stray grains, we have introduced a threshold of at least twenty grains to be entrained for either η > 0.

For 6 ≲ a₀ [μm] ≤ 11, entrainment fractions may vary by a factor of up to 2; in addition, the peaks are more pronounced due to the lower-number statistics forlarger grains (see Table 1).

So while we may expect the dust content of the wind to vary according to the strength of the turbulent mixing, this should not be detrimental to the rest of our findings.

Furthermore, the maximum entrained grain size is affected merely marginally; this can be checked when comparing Figs. B.3 and B.4 to Figs. 6 and 7, respectively. max (a₀) does vary slightly along the base of the wind; this is to be expected due to the randomized initial placement of the dust grains. Yet the position of the peak is well-preserved – just as for σ(v) = 0, we find max (a₀) = 11 μm at R ≃ 20 AU. The additional local maxima stem from the numerical unevenness of the launching area, demonstrated in Fig. B.5.

Fig. B.1 Acceleration of entrained grains of a0 = 10 μm: wherever there is at least one particle with vr < vr,esc, the area is colored in dark blue; if all grains have reached vr > vr,esc, it is yellow; cyan areas are not traversed by any grains of this size. Even for these comparatively large grains, vr,esc is reached well within the computational domain (if at all).

Fig. B.2 Ratio of the fraction of grains entrained for σ(v) = 0 (i.e. η0) and the fraction of grains entrained for σ(v) = 100 m s−1 (i.e. ησ(v)) for bins of 10 AU along the launching area. For clarity, we required at least 20 grains to be entrained for either η > 0. While there are distinct deviations between η0 and ησ(v), they do not exhibit a clear pattern. Dips to 0 indicate regions where σ(v) > 0 allows for additional wind pick-up of some dust particles that could not have been entrained with σ(v) = 0.

Fig. B.3 Maximum entrained grain size max (a0) along the disk surface when σ(vi) = 100 m s−1 is included in the grain setup. The differences to Fig. 6, which forgoes the velocity spread, are minor, and mostly due to numerical variations of the initial positioning.

Fig. B.4 Maximum entrained grain size max (a0) in 2D, to becompared to Fig. 7. As with Figs. 6 and B.3, the differences caused by σ(vi) are minor at best.

Fig. B.5 ϑ-coordinate of the base of the photoevaporative flow. Its (slightly) craggy shape is especially pronounced around 50 ≲ R [AU] ≲ 90, which leads to small numerical artifacts in our results.

If a grain with a velocity vector that is especially enhanced in direction of the wind motion is launched from an edge point along the base line, this may allow it to enter the wind region in contrast to another particle starting from the same location, but with a more downwards-biased velocity vector. The 2D maps for max (a₀) in Figs. 7 and B.4 are almost identical, apart from few very narrow lines from individual grains of a₀ ≥ 9 μm; it stands to reason that these are caused by the randomness in the initial conditions just described.

The fits to the parameters c₁ and c₂ of Eq. (1), shown in Fig. 9, are almost entirely unaffected by the additional σ(v_i); for brevity, suffice to state that both c₁ and c₂ deviate but in the third significant figure.

To summarize, while we must assume that there is a sensible amount of gas turbulence around the disk surface, this does not strongly affect our findings.

Appendix C Plots and parameters for all dust grain sizes

We have modelled 20 distinct grain sizes, listed in Fig. 2 and Table 1; for clarity we have chosen to only show plots for three distinct sizes above. In the following, we present the corresponding plots for all a₀, and argue why we have chosen exactly these to represent the full sample.

An arbitrary selection of dust grain trajectories, spaced by roughly 5 AU intervals, and entering the wind-dominated region at various points along its base is shown in Fig. C.1, the full version of Fig. 3. The individual trajectories are colored according to their local Stokes number St.

There are four main scenarios to be found here:

Firstly, full wind entrainment. This applies if St ≪ 1 throughout the trajectories, that is for a₀ ≲ 0.1 μm. Visually, the trajectories shown for these grain sizes are mostly blue, indicating St < 0.1. As discussed above, these dust particles follow the gas motion very closely; this causes grains to fall back below the base of the wind along the gas streamlines, that is for R ≳ 180 AU. The most massive grains of this grouphave sizes a₀ = 0.1 μm, which we have chosen for the plots in Sect. 3.

Secondly, slow decoupling from the gas flow. For 0.5 ≲ a₀ [μm] ≲ 5, the grains are picked up with St ≪ 1. Yet while the particles are blown out and hence within r < 300 AU, their St approaches 1. Thus the trajectories decouple from the gas; high above the disk surface, they fall below the gas streamlines, whereas close to it, they deviate upwardly, leading to a more radial outflow in all cases. Grain fallback occurs for R ≳ 170 AU. a₀ = 4 μm represents a grain size from this interval for which St → 1 is readily apparent.

Thirdly, quick decoupling. For 6 ≲ a₀ [μm] ≤ 11, the interval during which the dust grains follow the gas stream is shorter, and the grains reach higher St while being blown out. For R ≳ 160 AU, we observe particles falling back below the base of the flow because the wind cannot provide enough momentum. The almost-largest grain size in this group is a₀ = 10 μm (11 μm grains are more sparse in the wind region).

Fig. C.1 Full version of Fig. 3: arbitrary selection of dust grain trajectories for all 20 a0, colored by their local St. Wind base in black, gas streamlines in 5% steps of Ṁw in dash-dotted grey.

Fig. C.2 Times needed to fully blow out dust particles from their starting position along the base of the wind to the domain boundary (i.e. Δtbnd) and to accelerate them to vr,esc (i.e. Δtesc). Individual data points in cyan and orange and mean values in blue and red, respectively. The rasterization of the former points results from a time-discrete particle tracking and a binning in R-direction; a raster point may therefore represent multiple data points. Keplerian orbital times at the base of the wind are included as dashed green lines.

Fourthly, no wind pick-up. For a₀ > 11 μm, the grains are too heavy to be lifted up by the wind; even if they reach the disk surface, they fall back below it and successively follow the gas streams in the disk because of the comparably high densities there. Since we find no entrainment for these grains, we have omitted this group from the plots in Sect. 3.

Expanding on Fig. 5, Fig. C.2 shows the times it takes the dust grains to reach r ≳ 300 AU after wind pick-up – labelled Δt_bnd (data points in cyan, mean in blue) – and the (much lower) times for acceleration to the escape velocity v_r,esc – labelled Δt_esc (data points in orange, mean in red)¹⁷.

Choosing the same size categories as above, with the panels for a₀ > 11 μm omitted since such grains are not entrained by the photoevaporative flow, we find the following:

Firstly, very small dust grains. These are accelerated very strongly when entering the wind close to the star, resulting in low Δt_bnd and even lower Δt_esc; at higher r, the timescales increase. At R ≈ 80 AU, Δt_bnd peaks and starts falling off again because grains launching from further out have less distance to cover to the computational boundary at r ≃ 300 AU. Entrainment occurs for R ≲ 140...160 AU, depending on a₀; the according time frames span 70 ≲ Δt_bnd [yr] ≲ 170 and 5 ≲ Δt_esc [yr] ≲ 50. The timescale of the particle motion in the wind is smaller than the Keplerian orbital timescale for grains launched from R ≳ 20 AU (Δt_bnd) or R ≳ 1 AU (Δt_esc), meaning that the motion in the wind dominates the dynamic evolution.

Secondly, small grains. They show a monotonic increase of Δt_bnd and Δt_esc with R; the peak we have seen for smaller a₀ disappears. Hence, the farther out a dust grain is picked up by the XEUV wind, the longer it takes to be blown out of the computational domain and to be accelerated to v_r,esc; while this may appear counter-intuitive at first, it merely illustrates that the gas flow is much stronger closer to the star. The grain transport happens on a timescale of a few yr to a few 10² yr; the Keplerian motion dominates the dynamic timescale for R ≳ 30 AU (Δt_bnd) or R ≳ 5 AU (Δt_esc).

Thirdly, medium-sized grains. As for the small grains, we see a mostly monotonic relation between R and Δt_bnd. However, the graph exhibits very distinct peaks for a₀ [μm] ∈{8;9}, which we have not commented on in Sects. 2 and 4 because they are numerical artifacts caused by the craggy launching area (see Sect. 3 and Appendix B), and interestingly disappear in the setup including an initial velocity spread¹⁸. Apart from that, we find t_dyn < Δt_bnd for R ≲ 60 AU and an almost identical value for Δt_esc, indicating that speed-up is slower for heavier grains, as would be expected.

For R → 0, we find Δt_esc → 10 yr for a₀ ≲ 9 μm. This is another indicator that the photoevaporative flow is strongest close to the star, where it is most effective at accelerating the dust grains.

Last, Fig. C.3 shows the areas occupied by the dust grains, and the according fits for ${\text{max}(z)|}_R$. The fit parameters are also listed in Table C.1, with c₁ and c₂ as defined in Eq. (1).

Firstly, as claimed in Sect. 3, these plots illustrate that all a₀ < max (a₀) populate the regions occupied by the max (a₀) shown in Fig. 7. Secondly, the unpopulated regions found in Fig. C.3 are once again caused by the slight bumpiness of the base of the flow; additional test runs with a much increased spatial resolution (i.e. 1000 particles per 1 AU along the launching plane) were used to confirm this – the plots for them look very similar (and have therefore not been included here).

Fig. C.3 Full version of Fig. 8: dust population in the wind for all a0. We see a clear correlation between R and max (z) in the individual plots, fitted with dashed blue lines. The parameters given in the plots also listed in Table C.1. Base of the wind in red, non-fitted population boundaries in orange.

References

All Tables

All Figures

Fig. 1 Density (green) and velocity map (yellow arrows: v) of the gas disk model from Picogna et al. (2019). The base of the XEUV flow (i.e. the location of the largest gradient in temperature) is shown in black, and coincides with a strong drop in density. vr is pointing radially outwards everywhere in the wind but very close to the sonic surface where v points outwards, but away from the disk (for R < 200 AU). The cyan line represents z(R) = H; additional lines for 3.5 H (dotted blue) and 5.5 H (dashed blue) show the range of scale heights the disk surface covers.

In the text

Fig. 2 Left: initial placement of the dust grains at Δt = 0 at the base of the photoevaporative flow (black), slightly below the sonic surface. The grains are colored according to their size and the color scale (scaling with $\sqrt{\left|v_r\right|}$) of the gas disk represents the extent of its local radial velocity (with − 2 ≲ vr [km s−1] ≲ 30). The gas map is mostly smooth in the region of interest, that is everywhere but close to the midplane at high R. Right: snapshot of the grain positions at Δt ≃ 100 yr, all else equal. The very low spread of the lines of individual grain sizes is due to the initial setup, placing particles directly on the disk surface without a spread in their initial velocities. For comments on this, see also Sect. 3 and Appendix B.

In the text

Fig. 3 Randomly selected dust trajectories for a0 = 0.1, 4, and10 μm (left, middle, and right panels, respectively). The trajectory color represents the local value of the Stokes number. Entrained dust grains, launching from the launching region (black), migrate upwards on the colorbar the as they move to regions of lower gas density. Gas streamlines are shown in dash-dotted grey. (For all simulated grain sizes, see Fig. C.1.)

In the text

Fig. 4 Analysis of two dust trajectories (left: a0 = 0.1 μm, right: a0 = 10 μm) entrained in the photoevaporative outflow from R ≃ 20 AU. Top panels: actual motion in the (R, z)-plane colored by St. Lower panels, from top to bottom: tstop and ϱgas, St and ∂t St, the horizontal speed vR of gas and dust, the vertical speed vz of gas and dust, a comparison of the direction of the motion vz∕vx for gas and dust, and ϑ and vϑ. See the text for an in-depth commentary.

In the text

Fig. 5 Time needed to fully blow out dust particles from when they first enter the wind at twind,0 to the domain boundary, which they reach at tbnd (data points in cyan, mean in blue), and time needed to accelerate dust grains to vr,esc (data points in orange, mean in red); both for a selection of three a0. Keplerian orbital times (tdyn ≡ ΔtK) at the disk surface are included as dashed green lines. The rasterization of the data points results from a time-discrete particle tracking and a binning in R-direction; thus, one raster point may represent multiple data points. (For all a0, see Fig. C.2.)

In the text

Fig. 6 Size of the largest grains entrained from a point R along the base of the wind (blue, peak at about 20 AU); the saw-tooth appearance of the curve at larger R is caused by the finite resolution of the underlying gas grid. When comparing to Hutchison et al. (2016a, their Fig. 7 for M* = 0.75 M⊙) (orange, peak at around 40 AU), we can see the size enhancement – especially at smaller R – caused by the inclusion of X-rays in our photoevaporative wind model. The blue dashed line represents our results scaled down by a factor of 3, to compensate for the differing internal grain densities of Hutchison et al. (2016a) and this work; yet, Ṁw still differsbetween the models, making a direct comparison difficult. The MHD wind model investigated by Miyake et al. (2016, their Fig. 4) (green) shows a distinctly different entrainment curve, starting off at very high a0 in the jet region but dropping towards max (a0) = 0 very quickly.Around R ≃ 140 AU, our max (a0) plummets to 0.

In the text

	Fig. 7 Maximum size max (a0) of thedust grains in the wind in analogy to Owen et al. (2011a, their Fig. 2). The brightest color found above the sonic surface indicates a global maximum of max (a0) = 11 μm; the region below the base of the wind (black) is included merely for completeness. The visible correlation between grain size and maximum scale height is further investigated in Fig. 8.
In the text

	Fig. 8 Areas populated by the dust grains (green); since we do not start from a realistic distribution along the base of the wind, we do not portray a density map. Smaller grains reach higher z at similar R. Wind base in black, numerical ${\text{max}(z)|}_R$ in orange, and corresponding fit in (dotted) blue; fits according to Eq. (1) annotated. (For all a0, see Fig. C.3.)
In the text

	Fig. 9 Fit parameters for all fits for ${\text{max}(a_0)|}_R$ as a functionof a0 in blue, with 3σ-errors included as shaded regions. Non-physical fits to these fit parameter curves in orange, with the parameterizations given in the text boxes.
In the text

Fig. 10 Re-plot of vR and vz for the same dust grain as in the right column of Fig. 4 (a0 = 10 μm). The dust parameters are colored red, the surrounding gas in blue; the motion according to Eqs. (1) and (3) of Giacalone et al. (2019), with tstop as retrieved from our model, is added in green. Especially towards larger R, the semi-analytical prescription does not reproduce the dust motion very well. It should be noted, however, that Giacalone et al. (2019) do not use it for grains this big.

In the text

Fig. A.1 Maximum of the per-a0 mean values for |Fdisk,z∕Fdrag,z| (see colorbar), mapped to a 2 AU × 2 AU grid. The gas drag clearly dominates over the (overestimated) disk gravity, especially in the wind region (base of the wind in black). Empty (white) regions above the sonic surface indicate that no particle has been recorded while in this cell, owing to the high grain speeds there (in contrast to Sect. 3, we did not interpolate between recorded particle positions).

In the text

	Fig. B.1 Acceleration of entrained grains of a0 = 10 μm: wherever there is at least one particle with vr < vr,esc, the area is colored in dark blue; if all grains have reached vr > vr,esc, it is yellow; cyan areas are not traversed by any grains of this size. Even for these comparatively large grains, vr,esc is reached well within the computational domain (if at all).
In the text

Fig. B.2 Ratio of the fraction of grains entrained for σ(v) = 0 (i.e. η0) and the fraction of grains entrained for σ(v) = 100 m s−1 (i.e. ησ(v)) for bins of 10 AU along the launching area. For clarity, we required at least 20 grains to be entrained for either η > 0. While there are distinct deviations between η0 and ησ(v), they do not exhibit a clear pattern. Dips to 0 indicate regions where σ(v) > 0 allows for additional wind pick-up of some dust particles that could not have been entrained with σ(v) = 0.

In the text

	Fig. B.3 Maximum entrained grain size max (a0) along the disk surface when σ(vi) = 100 m s−1 is included in the grain setup. The differences to Fig. 6, which forgoes the velocity spread, are minor, and mostly due to numerical variations of the initial positioning.
In the text

	Fig. B.4 Maximum entrained grain size max (a0) in 2D, to becompared to Fig. 7. As with Figs. 6 and B.3, the differences caused by σ(vi) are minor at best.
In the text

	Fig. B.5 ϑ-coordinate of the base of the photoevaporative flow. Its (slightly) craggy shape is especially pronounced around 50 ≲ R [AU] ≲ 90, which leads to small numerical artifacts in our results.
In the text

	Fig. C.1 Full version of Fig. 3: arbitrary selection of dust grain trajectories for all 20 a0, colored by their local St. Wind base in black, gas streamlines in 5% steps of Ṁw in dash-dotted grey.
In the text

Fig. C.2 Times needed to fully blow out dust particles from their starting position along the base of the wind to the domain boundary (i.e. Δtbnd) and to accelerate them to vr,esc (i.e. Δtesc). Individual data points in cyan and orange and mean values in blue and red, respectively. The rasterization of the former points results from a time-discrete particle tracking and a binning in R-direction; a raster point may therefore represent multiple data points. Keplerian orbital times at the base of the wind are included as dashed green lines.

In the text

	Fig. C.3 Full version of Fig. 8: dust population in the wind for all a0. We see a clear correlation between R and max (z) in the individual plots, fitted with dashed blue lines. The parameters given in the plots also listed in Table C.1. Base of the wind in red, non-fitted population boundaries in orange.
In the text