The dispersal of protoplanetary discs I: A new generation of X-ray photoevaporation models

TL;DR

This paper introduces a new set of X-ray photoevaporation models for solar-type stars, improving the understanding of disc dispersal and its impact on planet formation through detailed simulations and population analysis.

Contribution

It presents a novel hydrodynamical modeling approach with an improved temperature parameterization, providing recipes for disc mass loss rates applicable to population synthesis.

Findings

Mass loss rates depend on stellar X-ray luminosity.

Models produce line profiles consistent with observations.

Approximately half of observed transition disc features are explained.

Abstract

Photoevaporation of planet forming discs by high energy radiation from the central star is potentially a crucial mechanism for disc evolution and it may play an important role in the formation and evolution of planetary system. We present here a new generation of X-ray photoevaporation models for solar-type stars, based on a new set of hydrodynamical simulations, which account for stellar irradiation via a new, significantly improved, parameterisation of gas temperatures, based on detailed photoionisation and radiation transfer calculations. This is the first of a series of papers aiming at providing a library of models which cover the observed parameter space in stellar and disc mass, metallicity and stellar X-ray properties. We focus here on solar-type stars (0.7 M$_\odot$) with relatively low-mass discs (1% of the stellar mass) and explore the dependence of the wind mass loss rates on stellar X-ray luminosity. We model primordial discs as well as transition discs at various stages of evolution. Our 2D hydrodynamical models are then used to derive simple recipes for the mass loss rates that are suitable for inclusion in one-dimensional disc evolution and/or planet formation models typically employed for population synthesis studies. Line profiles from typical wind diagnostics ([OI] 6300 $\overset{\lower.5em\circ}{\mathrm{A}}$ and [NeII] 12.8 $\mu$m) are also calculated for our models and found to be roughly in agreement with previous studies. Finally, we perform a population study of transition discs by means of one-dimensional viscous evolution models including our new photoevaporation prescription and find that roughly a half of observed transition discs cavities and accretion rates could be reproduced by our models.

Tables

Table 1. Table 1: Coefficients for the Temperature fit

$N_H$	b	c	d	m
$1\text{×}{10}^{20}\mathrm{pp}/{\mathrm{cm}}^2$
$0-25$	$-26.5735$	$-6.5924$	$4.0359$	$0.2086$
$25-50$	$-53.6260$	$-5.8336$	$3.9421$	$0.1220$
$50-75$	$-65.0650$	$-5.7396$	$3.9160$	$0.1027$
$75-100$	$-33.4958$	$-5.4189$	$3.9088$	$0.2039$
$100-125$	$-32.0490$	$-5.2748$	$3.8974$	$0.1972$
$125-150$	$-35.0241$	$-5.1722$	$3.8840$	$0.1699$
$150-175$	$-40.8293$	$-5.1922$	$3.8705$	$0.1422$
$175-200$	$-36.0245$	$-5.2231$	$3.8720$	$0.1527$
$200-225$	$-54.3903$	$-5.0610$	$3.8584$	$0.1053$
$225-250$	$-36.0854$	$-4.9851$	$3.8635$	$0.1579$

Table 2. Table 2: Parameter space

variable	value
disc extent
radial [au]	$0.33-1000$
polar [rad]	$0.005-\pi /2$
grid resolution
radial	$412$
polar	$160$
physical properties
${\log }_{10}\left(L_X\right)$ [erg/s]	$28.3,29.3,29.8,30.3,31.3$
$R_{\mathrm{hole}}$ [au]	$6.2,9.9,17.9,28.2$
$\mu$	$1.37125$

Table 3. Table 3: Mass-loss rates

simulation	Rhole	$\dot{𝐌}$
	au	${10}^{-8}$ M⊙ yr-1
PR	$0.33$	$2.644$
TR5	$6.2$	$2.940$
TR10	$9.9$	$3.036$
TR20	$17.9$	$2.923$
TR30	$28.2$	$2.887$

Table 4. Table 4: Line profiles from primordial discs irradiated by log 10 ⁡ ( L X ) subscript 10 subscript 𝐿 X \log_{10}{(L_{\mathrm{X}})} = 30.3   erg ​ s − 1 times 30.3 erg superscript s 1 30.3\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\ \mathrm{s}^{-1} . The table lists the velocity of the peak and the full width at half-maximum (FWHM) for 4 4 4 inclinations from 0 0 to 90   ∘ times 90 90\text{\,}{}^{\circ} . The profiles were degraded to an instrumental resolution of R = 100 , 000 𝑅 100 000 R=100,000 .

	Ne II $12.8\mu m$		[OI] $6300\AA$
i	vpeak	FWHM	vpeak	FWHM
∘	km s-1	km s-1	km s-1	km s-1
$0$	$-3.56\pm 0.01$	$11.50$	$-6.707\pm 0.006$	$12.427$
$30$	$-3.94\pm 0.03$	$15.94$	$-6.28\pm 0.02$	$14.22$
$60$	$-1.93\pm 0.06$	$18.15$	$-3.21\pm 0.05$	$15.97$
$90$	$0.10\pm 0.04$	$8.25$	$0.11\pm 0.01$	$16.36$

Equations

$${\log }_{10}(T(\xi ))=d+\frac{1.5-d}{{\left[1.0+({\log }_{10}(\xi )/c{)}^b\right]}^m}$$

$${\dot{\Sigma }}_w(R)$$

$${\dot{M}}_w(R)=\dot{M}(L_X)10^{a\log R^6+b\log R^5+c\log R^4+d\log R^3+e\log R^2+f\log R+g}$$

$${\dot{\Sigma }}_{w,T}(R)=a{b}^x{x}^{c-1}[x\ln (b)+c]\frac{1.12\dot{M}(L_X)}{2\pi R}[M_{\odot }{au}^{-2}{yr}^{-1}]$$

$${\log }_{10}(\dot{M}(L_X))=A_{\text{L}}\ast \exp \left[\frac{(\ln ({\log }_{10}(L_X))-B_{\text{L}}{)}^2}{C_{\text{L}}}\right]+D_{\text{L}},$$

$$\alpha =2\cdot 10^{-3}+7.5\cdot 10^{-5}\cdot R_1.$$

$${\tau }_{\mathrm{rec}}=\frac{1}{{\alpha }_{\mathrm{A}}(T_{\mathrm{e}})n_{\mathrm{e}}}\simeq 1.5\cdot 10^9{T}_{\mathrm{e}}^{0.8}{n}_{\mathrm{e}}^{-1},$$

$${\tau }_{\mathrm{adv}}>{\tau }_{\mathrm{rec}}.$$

Full text

††pubyear: 2019††volume: 487††pagerange: The dispersal of protoplanetary discs I: A new generation of X-ray photoevaporation models–14

1 Introduction

The modality of disc dispersal is thought to be of fundamental importance to planet formation, yet the responsible mechanism is still largely unconstrained. Photoevaporation from the central star is currently a promising avenue to investigate, but the models developed to date do not yet have enough predictive power for a detailed comparison with the observations. Photoevaporative profiles and rates play a very important role in understanding disc dispersal as well as in the feedback of this process on planet formation and migration (e.g. Ercolano & Rosotti, 2015; Ercolano & Pascucci, 2017; Ercolano et al., 2017; Carrera et al., 2017; Jennings et al., 2018).

Previous efforts trying to infer wind mass-loss rates and profiles from observed atomic and ionic emission line intensities and profiles have been inconclusive. These have included observational and theoretical studies of the Ne II fine structure line at $12.8\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (Herczeg et al., 2007; Glassgold et al., 2007; Pascucci & Sterzik, 2009; Schisano et al., 2010; Ercolano & Owen, 2010), the [OI] collisionally excited lines at $6300\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ and $5577\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (Rigliaco et al., 2014; Simon et al., 2016; Ercolano & Owen, 2016) as well as a number of other optical collisionally excited metal lines (Natta et al., 2014; Canovas et al., 2018). As discussed in Ercolano & Owen (2016) collisionally excited lines are unsuitable to infer mass-loss-rates or to constrain the wind driving mechanism due to their strong temperature dependence allowing them only to trace the hot layer of the wind heated by the EUV radiation and not the bulk of the wind. Ercolano & Pascucci (2017) further discuss this issue also in the context of new observational evidence of two distinct components in the emission lines of transition discs, possibly associated with MHD and photo-evaporative winds (see e.g. Simon et al., 2016).

Mid-infrared observations of molecular species (e.g. CO) may provide a viable alternative, and indeed recent observations (e.g. Pontoppidan et al., 2011; Brown et al., 2013; Klaassen et al., 2013; Klaassen et al., 2016) have pointed out that these lines may be tracing a slow and partially molecular wind. Current photoevaporative models however lack the ingredients necessary to quantitatively model these types of observations. Specifically, the models of Owen et al. (2010); Owen et al. (2011, 2012) and Alexander et al. (2006a, b) perform hydrodynamical calculations of the photoevaporative wind but do not include a chemical network beyond the atomic state. The models of Gorti et al. (2009), on the other hand, include a thermochemical calculation, but only consider hydrostatic disc atmospheres, thus precluding the calculation of line profiles. A first step in this direction has been made by Wang & Goodman (2017b) and Nakatani et al. (2018), where both the hydrodynamics and the chemistry are evolved together, employing a simplified radiative transfer scheme.

We have started a comprehensive effort, in the context of a coordinated program aiming at understanding the nature of transition discs 111DFG Research Unit FOR: 2634/1 ”Planet Formation Witnesses and Probes: Transition Discs”, www.transitiondiscs.com, that finally aims at performing quantitative spectroscopy of disc winds to allow the direct determination of their strength (mass-loss rates and profiles). To achieve these aims we are building new radiation-hydrodynamical calculations of irradiated discs, coupled with photoionisation, chemistry and radiative transfer calculations for a large parameter space, covering stars of different X-ray properties. The comparison of our model grids with existing and upcoming observations will allow us to constrain the mass-loss rates and the launching regions of the wind and thus pin down the underlying driving disc dispersal mechanism.

This paper is the first in a series, where a new, significantly improved, set of hydrodynamical models of X-ray photoevaporation, using an improved temperature scheme are presented. In a following paper, the wind solutions obtained here will be coupled with a chemical network to provide a first assessment of possible molecular diagnostics. As well as providing the basis for our future investigation, this paper specifically aims at addressing a number of questions left open by the widely used models of Owen et al. (2010); Owen et al. (2011, 2012); Owen (2016). Namely, the new temperature scheme allows to follow more carefully the evolution of the outer disc during the transition phase addressing the current discrepancy between the observations and the prediction of non-accreting large hole transition discs (Owen et al. 2012, but also see Ercolano et al. 2018, Wölfer et al. 2019, in prep.). The evolution of the gas surface density during the dispersal phase is an important ingredient to planet formation models based on population synthesis (e.g. Mordasini et al., 2012). To that aim, we provide here an accurate scheme that can be easily imported in population synthesis codes.

The paper is organised as follows. In Section 2 we describe the numerical methods used and the parameter space studied, focussing on the improvement with respect to previous studies. In Section 3 we describe the main results, providing analytical fits. Finally, in Section 4 we draw the main conclusions from this work and outline the future development.

2 Methods

Photoevaporative flows originate from gas that is thermally heated from the central star, launched by the pressure gradient and then centrifugally accelerated. A hydrodynamical calculation of such flows thus requires the calculation of gas temperatures at each hydrodynamical time-step. Three approaches have been proposed in the literature: (i) use temperature parameterisations, from detailed radiative transfer calculations that solve the heating and cooling equations, for a range of physical and irradiation properties (e.g. Owen et al., 2010; Owen et al., 2011, 2012; Haworth et al., 2016). (ii) perform a radiative transfer and thermochemical calculation at each hydrodynamical step (e.g. Wang & Goodman, 2017b; Nakatani et al., 2018). (iii) Alternatively, an analytic model of the self-similar solution of thermal winds topology has been derived for an isothermal disc (Clarke & Alexander, 2016).

In this work, we choose the first approach (see below) as this allows for a much more accurate estimation of the temperatures, which is crucial for our aims. This is because the computational costs of performing a radiative transfer and thermal calculation at each hydrodynamical time-step are so high that extreme simplifications must be adopted in order to make the task viable. This point is discussed in more detail throughout the paper, where our results are considered in the context of previous work, including the recent calculations by Wang & Goodman (2017b).

2.1 Thermal Calculations

We have used the gas photoionisation and dust radiative transfer code mocassin (Ercolano et al., 2003; Ercolano et al., 2005; Ercolano et al., 2008a) to obtain the gas temperature for typical protoplanetary disc conditions, irradiated by the X-ray+EUV spectrum presented by Ercolano et al. (2008b); Ercolano et al. (2009). This synthetic thermal spectrum was generated by the plasma code of Kashyap & Drake (2000), assuming an emission measure distribution based on that derived from RS CVn type binaries by Sanz-Forcada et al. (2002) (which peaks at ${10}^{4}\text{\,}\mathrm{K}$) and fits Chandra spectra of T Tauri stars by Maggio et al. (2007).

We assume solar abundances for the gas (Asplund et al., 2005), depleted for the amount of gas locked in dust grains, according to (Savage & Sembach, 1996). Possibly, significantly larger amounts of metal depletion are expected to occur and have indeed been reported in recent observations of protoplanetary discs (e.g. Hogerheijde et al., 2011; Favre et al., 2013; Ansdell et al., 2016; Kama et al., 2016; Du et al., 2017; Miotello et al., 2017). The effect of metal depletion on photoevaporation rates can be severe (Ercolano & Clarke, 2010) and can have dramatic effects on the surface density evolution of discs (Ercolano et al., 2018). A detailed investigation of this effect on the wind rates and profiles by means of hydrodynamical simulations similar to those presented here is currently underway and will be presented in a forthcoming paper (Wölfer et al., in preparation).

We considered $10$ different slabs of $2.5\text{\times}{10}^{21}\text{\,}\mathrm{p}\mathrm{p}\,\mathrm{c}\mathrm{m}^{-2}$ each, up to a column density of $2.5\text{\times}{10}^{22}\text{\,}\mathrm{p}\mathrm{p}\,\mathrm{c}\mathrm{m}^{-2}$, and obtained for each a gas temperature prescription. The column is resolved in 100 grid points and we checked that larger resolution does not change the result. Beyond the maximum column density, we assume thermal coupling of the gas and dust, and the dust temperatures from the models of D’Alessio et al. (2001) are mapped to our models.

The models consist of 3D slabs with reflecting boundary conditions on the xy and zy planes and outflowing boundary conditions on the xz planes. The energy packets are injected plane parallel, but they are allowed to be absorbed, re-emitted and scattered in 3D. The main heating channel is X-ray photoionisation, although photoelectric heating from dust grains is also included, but is mostly negligible. Cooling channels include continuum process (2-photon continuum, free-bound recombination, free-free), however the dominant cooling mechanisms in these models are collisionally excited (forbidden) lines and recombination lines, which are generally optically thin and escape without further interactions, except in the case of resonant emission lines, which could be absorbed by dust and are treated using and escape probability approach (see Ercolano et al. 2005). HeI and HeII recombination lines (the latter not relevant in these models) emitted at energies $>$ 13.6eV are followed as they could go on ionising H or HeI respectively.

Based on our model results we propose a new parameterisation of the gas temperature as a function of ionisation parameter $\xi=L_{X}/nr^{2}$ (Tarter et al., 1969) and gas column density.

The following relation produces a good fit to the model data:

[TABLE]

where the coefficients, as a function of the hydrogen column densities $N_{H}$, are given in Table 1.

Figure 2 shows a comparison of our parameterisation scheme to the one used by Owen et al. (2010); Owen et al. (2011, 2012), which was independent of the gas column density, and is likely a column density biased average. Our new radiative transfer and thermal calculations extend to much lower values of $\xi$ (down to $\log{\xi}=-8$, instead of $\log{\xi}=-6$), allowing us to better model the outer regions of the disc, which are particularly relevant to follow the evolution of transition discs. The higher maximum temperature reached by the previous prescription (as seen in the upper right part of Fig. 2) is due to an integration over a finer grid, which allowed them to solve the low-density region heated by the EUV radiation. We chose not to over-resolve this region since it does not contribute to the total mass-loss rate. Including the column density as a parameter, allows us to better reproduce the temperatures at different locations of the disc. While the single-slab parameterisation of Owen et al. (2010); Owen et al. (2011, 2012) can produce errors of order $30\text{\,}\mathrm{\char 37\relax}$ for column densities below ${10}^{22}\text{\,}\mathrm{p}\mathrm{p}\,\mathrm{c}\mathrm{m}^{-2}$, we reduce this value to less than $1\text{\,}\mathrm{\char 37\relax}$ as shown in Figure 1

The consequences of the new temperature prescription on the photoevaporation rates and profiles are discussed in more details in Section 3.

Finally, Haworth et al. (2016) presented a parameterisation scheme of temperature as a function of ionisation parameter to model the very late evolution of transition disc and the onset of the thermal sweeping instability (Owen et al., 2013). We draw attention to the fact, also discussed in Haworth et al. (2016), that the temperature ionisation prescription presented in their work (see their Figure 1) is only valid under the assumptions made there. Haworth et al. (2016) used gas that was optically thin to low energy X-ray photons to compute the temperatures and assumed that the profile is independent of column density. These assumptions are not justified in our disc models.

2.2 Hydrodynamics

We used the open source, modern hydro-code PLUTO (Mignone et al., 2007), adopting a spherical coordinate system centred onto the star. This choice of the reference frame allows us to determine, without increasing the computational load, the column density for each cell in the domain. We then modified the temperature calculation in the code in order to accurately model the stellar irradiation in the upper layers of the disc, using the parameterisation presented in the previous section, without impacting the code performance.

The initial density and temperature distributions for the hydrodynamical simulation of the primordial disc were taken from the hydrostatic equilibrium models of Ercolano et al. (2008b); Ercolano et al. (2009). We note that the choice of initial distribution does not affect the final result, but has an impact on the computational time required to achieve a steady state solution of the wind. In order to avoid numerical issues due to the very low density in the region near the pole and for large radii, the grid has been defined with an inverse logarithmic scaling in the polar direction, and a logarithmic scaling in the radial one, taking advantage of the polar coordinates.

Particular attention has been made to correctly model the outer boundary of the computational domain that may cause spurious oscillations (observed also in Wang & Goodman (2017b)), and that may affect the final result. In order to damp spurious numerical reflection, we adopted an outer boundary inside the computational domain. Practically, we set the right-hand side of the conservative equations in the region outside $980\text{\,}\mathrm{au}$ to zero. These cells are thus not evolved in time. A more detailed analysis of this issue is included in Appendix C.

The transition discs were derived from the primordial disc distribution by adding an exponential decay of the density close to the gap edge, from the hydrostatic models of Ercolano et al. (2008b); Ercolano et al. (2009). The disc is then readjusting its shape close to the gap as the system reaches a stable state where the gas inflow is balanced by the photoevaporative wind.

We checked the consistency of our results by increasing the grid resolution for our fiducial model (PR). The list of the systems modelled is reported in Table 3.

The parameter space analysis performed is summarised in Tab. 2.

2.3 Mass-loss rates calculations

We derive the mass-loss rates as a function of the cylindrical radius in the disc, adopting a similar approach as in Owen et al. (2010). First, we remapped the grid onto a Cartesian grid $2000\text{\times}2000$ in extent. Then, we defined an outer radius and followed the gas streamlines from it to the base of the flow, obtaining the radius at which the flow was launched. The outer radius was chosen to be $200\text{\,}\mathrm{a}\mathrm{u}$ in order to avoid the outermost regions, affected by numerical oscillations caused by the reflection of sound waves onto the outer boundary, and maximise the number of orbits performed at the given location within the simulation (see Appendix C for more details).

The flow base is defined at the location where a local maximum is observed in the temperature profile for each cylindrical radius. We checked that this definition was consistent also with the Bernoulli parameter (e.g. Wang & Goodman, 2017b).

3 Results

In this section, the results from our calculations are described in the context of previous analysis and the consequences of our improved temperature scheme on mass-loss rates and wind profiles are investigated. We provide, where possible, simple analytical formulae for implementation into disc evolution and planet formation models, which also include dependences on X-ray luminosities, as well as on the inner hole size in the case of transition discs. Transition disc demographics deriving from our new models are also calculated and compared with previous calculations and with current observations. Classical wind tracers, including the Ne II fine structure line at $12.8\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and the collisionally excited neutral oxygen line at $6300\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ are calculated for our fiducial model and compared to previous calculations by Ercolano & Owen (2016) and to observations. The raw model grids for each individual case are available from the authors upon request.

3.1 Mass-loss rates and profiles

Fig. 3 shows the density, temperature and radial velocity plots at the end of the simulation for the primordial disc in our reference model (PR). The streamlines of the wind flow are overlaid and plotted every $5\%$ intervals of the integrated mass-loss rate, from the outer radius, placed at $200\text{\,}\mathrm{au}$ (well inside the computational domain), down to the base of the flow. We can see that they originate mostly inside $50\text{\,}\mathrm{au}$ and they are radially spreading at larger radii as also seen by Owen et al. (2010); Wang & Goodman (2017b) as well as in pure EUV models (Font et al., 2004; Alexander et al., 2006a; Alexander, 2008). The red dashed line marks the sonic surface, the smoothness of which indicates that we are in a stable state condition.

The differences with the parameterisation by Owen et al. (2010) are subtle but important and are highlighted in Fig. 4. The new cumulative mass-loss rate is consistently higher between $5$ and $50\text{\,}\mathrm{au}$, boosting the total mass-loss rate which is reported in Tab. 3 and is $2$ times bigger than the one previously obtained. This is expected considering the arguments put forward in Owen et al. (2012), which is that the temperature at the sonic surface and hence the sound speed is fixed by the depth of the potential. At lower column densities a fixed temperature requires a lower ionization parameter for our new relations, thus yielding a higher mass loss rate as a consequence of the higher density.

This difference results also in a modified surface density mass-loss profile (as shown in Fig. 5), where two peaks are visible around $1$ au and $20$ au.

The fit to the surface density mass-loss profile for the primordial disc is:

[TABLE]

with

[TABLE]

and a = -0.5885, b = 4.3130, c: -12.1214, d = 16.3587, e = -11.4721, f = 5.7248, g = -2.8562. In the previous equations we use the standard convention and we refer to the natural logarithm with $\ln{}$ and the logarithm in base $10$ with $\log{}$. The radius is expressed in au, and $\dot{M}(L_{X})$ is the normalization factor given by eq. 5.

3.2 Hole Radius dependence

As described in Section 2.2, we created a set of transition discs from the primordial disc and for the reference X-ray luminosity of $2\text{\times}{10}^{30}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\ \mathrm{s}^{-1}$, in order to test the dependence of the mass-loss rate due to photoevaporation as a function of the inner hole radii. The result of our parameter study is reported in Tab. 3, as well as in Figure 7. The total mass-loss rate is increased with respect to the primordial disc, by a factor $\sim 1.12$, but it results independent of the inner radii of the disc. The surface density mass-loss profile for the transition discs are given by

[TABLE]

where $x=(R-R_{\mathrm{gap}})$, a = $1.1843\text{\times}{10}^{-1}\text{\,}\mathrm{,}$ b = $9.9695\text{\times}{10}^{-1}\text{\,}\mathrm{,}$ c = $4.8835\text{\times}{10}^{-1}\text{\,}\mathrm{,}$ and $\dot{M}(L_{X})$ is given by eq. 5.

3.3 Stellar X-ray luminosity dependence

The dependence of the integrated mass loss rate on the stellar X-ray luminosity was also investigated.

The result is shown in Fig. 8 and it can be described by a sigmoidal function:

[TABLE]

in $M_{\odot}\,\mathrm{yr}^{-1}$, with $A_{\textrm{L}}=-2.7326$, $B_{\textrm{L}}=3.3307$, $C_{\textrm{L}}=-2.9868\cdot 10^{-3}$, $D_{\textrm{L}}=-7.2580$. An interesting behaviour is observed, where for high luminosities, the integrated mass-loss rate tends to saturate around $1\text{\times}{10}^{-7}\text{\,}\mathrm{M}_{\odot}\mathrm{y}\mathrm{r}^{-1}$. At high X-ray luminosities, only the flat region of the $T-\xi$ relation (see Fig. 2) is accessible to the flow, and the theory from Owen & Clarke (2012) breaks down. In this regime, the disc jumps straight from the disc temperature to $1\text{\times}{10}^{4}\text{\,}\mathrm{K}$, and the depth of the potential is irrelevant, but the mass-loss rate is set solely by the disc structure. Thus the surface-mass loss profile becomes dependent on the X-ray luminosity as shown in Fig. 9, in contrast with Fig. 4 of Owen et al. (2011). Moreover, the disc structure is affected by the higher X-ray luminosity and the area accessible by the X-ray radiation is increased with higher X-ray luminosities since the disc is more puffed up.

3.4 Transition Disc demographics

Photoevaporation models, naturally produce transition discs as a result of inside-out dispersal. Other processes, including giant planet formation, could also produce dust-depleted cavities in discs. It is thus important to differentiate amongst the possible origins of observed transition discs, if the latter are to be used as laboratories to study disc dispersal and/or planet formation. The observed transition disc demographics span a large parameter range in terms of the size of the central cavity (from sub-au to hundreds of au) and the measured accretion rate onto the central star (from a few non-accreting sources to accretion rates similar to those of full discs). This points to the notion of transition discs being a very diverse class of objects (e.g. (e.g. Owen & Clarke, 2012; Owen, 2016; Ercolano & Pascucci, 2017), where the X-ray photoevaporation models of Owen et al. (2010); Owen et al. (2011), can only explain a fraction of the observed objects, namely discs with small (a few $\mathrm{au}$) cavities and low accretion rates (less than 1\text{\times}{10}^{-8}\text{\,}\mathrm{M_{\odot}}\,\mathrm{yr}^{-1}).

The wind profile has a strong influence on the evolution of the surface density of the disc (Ercolano & Rosotti, 2015; Ercolano et al., 2018; Jennings et al., 2018), for that reason, we have explored here the consequences of our improved X-ray photoevaporation model on the transition disc demographics. In Figure 10 we show the results of our transition disc population syntheses (panel a) and compare it to the calculations of Owen et al. (2011), shown in panel c of the same figure. The observations are represented by grey points and exclude stars of spectral type A, G and earlier (see Ercolano & Pascucci, 2017, for the individual references to the observational data).

We refer the reader to Owen et al. (2011) for details of the methods employed, which is based on one-dimensional viscous evolution models of discs under the effect of photoevaporation (see also Ercolano & Rosotti, 2015, for details of the code used for the calculations). We sample the X-ray luminosity function of the central stars stochastically using observational data for the Taurus region. We use a cumulative X-Ray luminosity function built from observed luminosities in the range $0.310\text{\,}\mathrm{keV}$ of pre-main sequence stars with masses $0.5\leq M\leq$1.0\text{,}\mathrm{M}_{\sun} in the Taurus cluster (Güdel et al., [2007](#bib.bib34)). In order to compare with Owen et al. ([2011](#bib.bib53)) we also initially assume a single disc set up consisting of a viscosity parameter, $\alpha=$2.5\text{\times}{10}^{-3} and initial disc scaling radius, $R_{1}=$18\text{,}\mathrm{au}$$.

The new calculations show a better overlap with the observations. Accreting transition discs with cavities up to $\sim 30\text{\,}\mathrm{au}$ are predicted by the new models, compared to $<10\text{\,}\mathrm{au}$ for the Owen et al. (2010); Owen et al. (2012) models. The reason for that is that the new profiles are more efficient at removing material at larger disc radii, resulting in the more efficient outward dispersal of the outer disc, whose inner edge can then reach larger radii before the inner disc is viscously accreted onto the star and the accretion signature is lost.

Our new models however still fall short of being able to explain the several observed objects with larger inner cavities and more vigorous accretion rates. The problem persists also if we relax the assumption of all sources having the same $\alpha$ and $R_{1}$ and follow instead the approach of Ercolano et al. (2018) by randomly sampling the initial scaling radius, $R_{1}$ from a uniform distribution allowing $R_{1}$ to vary between $10100\text{\,}\mathrm{au}$. In these models, $\alpha$ is then chosen accordingly in order to ensure that the dispersal timescale of the populations match observations (Güdel et al., 2007). Assuming a simple inverse dependence of $\alpha$ on $R_{1}$ one finds:

[TABLE]

The resulting population synthesis is shown in panel b of Figure 10, where we see that the relaxation of initial disc conditions leads to an increase in the inner cavity sizes of accreting transition discs that can be covered by the models up to $30\text{\,}\mathrm{a}\mathrm{u}$. However, as mentioned above, these models can still only explain half of the observed sources.

Large cavity, strongly accreting transition discs are often taken as possible signposts for (multiple) giant planet formation (e.g Zhu et al., 2011; Dodson-Robinson & Salyk, 2011) although this interpretation is not supported by theoretical models for many individual sources (e.g. Zhu et al., 2012; Owen, 2014, 2016) and also on statistical grounds (e.g. Dong & Dawson, 2016). Ercolano et al. (2018) provide an alternative explanation, where large cavity, strongly accreting transition discs are a result of X-ray photoevaporation in metal (C and O) depleted discs. This interpretation, however, needs to be confirmed by more rigorous hydrodynamic modelling at variable metallicities, which is the focus of a forthcoming paper (Wölfer et al. 2019, in prep.). An alternative explanation for the large accretion rates measured in large cavity TDs was also presented by Wang & Goodman (2017a), involving a magnetically driven supersonic accretion of low surface density material. Ercolano et al. (2018) however argue that this model is at odds with recent observations from the Atacama Large Millimiter Array (ALMA) of the inner cavities of some of these objects, showing that a large quantity of gas is present inside the dust cavity (Bruderer et al., 2014; van der Marel et al., 2015; van der Marel et al., 2016; Dong et al., 2017).

Finally, our models predict a large population of non-accreting transition discs with large cavities, which were named ’relics’ by Owen et al. (2011) and have not yet been observed. Owen et al. (2011) predict that $\sim 97\%$ of all transition discs should be relics. Our new models predict $\sim 88.8\%$ and $\sim 96.4\%$ of relics when fixed and variable initial disc parameters are assumed, respectively. As final thermal sweeping seems to be ruled out by recent calculations (Haworth et al., 2016), this discrepancy likely points to additional processes taking over disc dispersal at large radii (e.g. FUV-photoevaporation or enhanced X-ray photoevaporation at large radii in C-depleted sources, Wölfer et al., in prep., Owen & Kollemeier, in prep.).

3.5 The thermal and ionisation state of the wind and its observables

Ercolano & Owen (2010, 2016) presented a study of possible wind tracers, including the [Ne II] fine structure line at $12.8\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and a number of collisionally excited lines (CELs) from atomic and low ion species like C, O and S. These lines were targeted in those studies as they are often detected with few km/s blue-shifts in observational surveys (e.g. Rigliaco et al., 2013; Natta et al., 2014; Simon et al., 2016), thus bearing the signature of a slow-moving disc wind which may be photoevaporative in nature. While Ercolano & Owen (2010, 2016) show that X-ray photoevaporation models can indeed reproduce many of the reported observations (but see also Simon et al., 2016), Ercolano & Owen (2016) also demonstrate that CELs, with their strong temperature dependence, are not suitable as measures of mass-loss rates as they just trace a small heated volume of gas high in the wind rather than its bulk at its base (see also discussion in Ercolano & Pascucci, 2017). Nonetheless, the ability to reproduce the observed line intensities and their approximate profiles is an important requirement for a viable photoevaporation model.

One example is the observation of a blue-shift in the [OI] line at $6300\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, which immediately points to a warm quasi-neutral wind, hence ruling out EUV-driven photoevaporation, which produces a fully ionised wind, as in the case of the recent simulations by Wang & Goodman (2017b). In this context, we have calculated here line intensities and line profiles from our new models, via a post-processing of the hydrodynamical grids through the 3D photoionisation and dust radiative transfer code MOCASSIN (Ercolano et al., 2003; Ercolano et al., 2005; Ercolano et al., 2008a), and showed that these are in agreement with current observational data. The results are summarised in Fig. 11 and in Tab. 4. On the upper-left panel the surface density of the primordial disc (PR) is shown, with superimposed the emission regions of [OI] $6300\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ in yellow, and [NeII] $12.8\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ in purple. In the middle and low-left panels the relative emission lines are then displayed for $4$ different inclinations ($0^{\circ}$, $30^{\circ}$, $60^{\circ}$, and $90^{\circ}$). A zoom-in onto the innermost disc is presented in upper plots, in order to show that the emission regions becomes very close to the bound material in the disc, reducing the blue-shift of the produced lines. As a result, the peak for the [OI] in the edge-on configuration is significantly blue-shifted while, for the [NeII], the more extended bound component is reducing the blue-shift to only a few km/s (see also Tab. 4). On the right-hand side of Fig. 11 the same plots are shown for a transition disc (TR10), where the profiles are much more symmetric since both sides of the disc are visible through the central hole, thus allowing the detection of the red-shifted as well as the blue-shifted portion of the wind. A more detailed study of spectroscopic wind diagnostics predicted by our new models, compared with observations, will be presented in an forthcoming paper (Ercolano & Picogna, in prep.).

4 Conclusions

In this work, we performed hydrodynamical models of an X-EUV-driven photoevaporative wind for a primordial disc and transition discs with different inner holes. We improved on the models of Owen et al. (2010); Owen et al. (2011, 2012) by parameterising the temperature as a function not only of local properties of the gas (via the ionization parameter) but also on the column density to the star. This new approach resulted in new mass-loss profiles and total mass-loss rate of $2.644\text{\times}{10}^{-8}\text{\,}\mathrm{M}_{\odot}\mathrm{y}\mathrm{r}^{-1}$ (for a star X-ray luminosity of $2\text{\times}{10}^{30}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{s}^{-1}$). The main results from our models can be summarised as follows:

• •

the mass-loss rates generated by photoevaporative winds flattens for high stellar X-ray luminosities (Figure 8), reflecting the flattening of the temperature profile for high ionisation parameters (see Figure 2);

• •

the mass-loss rate does not depend on the inner-hole radii for transition discs (Figure 7);

• •

the demographics of transition discs obtained using the new scaling laws can explain a larger fraction of the observed transition discs (Figure 10), and in general it can explain accreting transition discs with cavities up to $30\text{\,}\mathrm{au}$;

• •

the line profiles generated from the hydrodynamical models can reproduce the two-components profile found by observational surveys (Figure 11), where the peak at $\theta=0$ is given by bound material in the inner part of the protoplanetary disc. These models, however, cannot reproduce the blueshift in the broad component that is often reported in the observations (Simon et al., 2016; Banzatti et al., 2018).

• •

the interaction between the photoevaporative wind and the outer boundary of the computational domain needs to be taken with caution, since it can generate spurious numerical oscillation, affecting the final mass-loss rate calculation (as observed in Wang & Goodman (2017b)).

• •

as observed in Appendix A, the differences in the cumulative mass-loss rate and surface density profiles are limited when comparing our numerical set-up with the one of Owen et al. (2010), pointing out that the results are quite robust. Nevertheless, the dependence of the temperature profile on the column density, leads to an increased mass-loss rate at large radii (see Fig. 4) and to the observed higher fraction of transition discs with large holes in the population synthesis models (see Fig. 10).

• •

even though we include EUV in the spectrum irradiating the planet-forming disc, this cannot dominate the total mass-loss rate (as proposed by Wang & Goodman (2017b)) because, as observed also in Owen et al. (2010), the EUV radiation field is absorbed for small column densities and it is not able to reach the high-density parts of the disc, thus its contibution to the total mass loss rate is negligible

• •

Wang & Goodman (2017b) found adiabatic cooling to be an important channel to the thermal balance in their simulations. We find, on the contrary (and in agreement with previous calculations by Owen et al. (2010) that this can be neglected, given that the energy flux gained in the wind is always much smaller ($<8\%$) that from the XEUV spectrum. Thus, we conclude that the gas is approximately in thermal equilibrium. Further to this point, we performed an explicit comparison of the advection and recombination timescales throughout our grids, finding the former to significantly exceed the latter (see App. B). Rather than advection, we suspect that a more likely source of discrepancy between the temperature structure of the models of Wang & Goodman (2017b) and ours (and Owen et al., 2010), is the limited sample of frequency points ($6$) used in their models (we use circa $1000$), which is unable to describe the ionisation structure of the gas.

The models presented here form the basis for future investigations aiming at finally being able to perform quantitative spectroscopy of disc winds from observations. Work is already in progress on the next steps, which include a study of the metallicity dependance of on the mass-loss rates (Wölfer et al., in preparation), detailed calculation of synthetic spectra from our calculations (Ercolano & Picogna, in preparation), and the inclusion of a realistic chemical network to determine the importance of FUV heating compared to the XEUV component in the different regions of the disc.

Acknowledgements

We thank the anonymous referees for the useful comments and suggestions which substantially helped improving the quality of the paper. GP acknowledges support from the DFG Research Unit “Transition Disks” (FOR 2634/1, ER 685/8-1). JEO is supported by a Royal Society University Research Fellowship.

We thank T. Haworth for the thermal sweeping models, T. Grassi for the insightful discussions, and L. Wölfer for the help with the MOCASSIN fits.

This work was performed partly on the computational resource ForHLR I funded by the Ministry of Science, Research and the Arts of Baden-Württemberg and the DFG, and partly on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP).

Appendix A Comparison test with previous model

This study is an improvement over the series of models made by Owen et al. (2010). As a first step we recreate the set-up of Owen et al. (2010) for their standard case (a primordial disc with $\log_{10}{(L_{x})}=30.3$) and a column-density independent temperature prescription. In Fig. 12 we show the direct comparison of the density, temperature and radial velocity structure, where no major difference is observed both in the cumulative value of the mass-loss rate, and in the surface density mass-loss rate.

Appendix B Radiative equilibrium

The approach used in this paper is based on the assumption of radiative equilibrium within the disc. This condition is met when the micro-physical time-scale, affecting the temperature equilibrium, is shorter than the hydrodynamical time-scale. Usually, hydrogen recombination is the microphysical process with the longest time-scale (Ferland, 1979; Salz et al., 2015)

[TABLE]

where $T_{\mathrm{e}}$ is the electron temperature, $n_{\mathrm{e}}$ is the electron density, and $\alpha_{\mathrm{A}}(T_{\mathrm{e}})$ is the temperature-dependant recombination rate. Thus, we did an a-posteriori check of the whole computational domain in order to check whether the hydrodynamical time-scale was longer than the hydrogen recombination time-scale in the regions important for the wind dynamics

[TABLE]

The result is shown in Fig. 13, where the fraction of the advection and recombination timescales is plotted. In the entire computational domain the hydrodynamical timescale is several order of magnitudes greater than the microphysical one, showing that our assumption is valid.

Appendix C Choice of the disc outer radius

The choice of the disc outer radius is important to discriminate between the material that is removed from the disc and the one that is simply redistributed within it. In principle, we can use our hydro model self-consistently to determine if the material becomes unbound or not. However, in order to have stable streamlines at hundreds of au from the central star, we need many more orbits than the ones modelled. We present in Fig. 14 the evolution of the total mass loss rate for the primordial disc considering the outer radius at the different locations. We see that in the first hundred orbits the mass-loss rate is increased for lower outer radii since we are considering streamlines that are falling back towards the disc mid-plane at further locations. Nevertheless, as the disc reaches a quasi-stable equilibrium, the differences become smaller. We decided to use the smallest radii which was not cutting out regions where the photoevaporation was effective, in order to have streamlines as stable as possible, thus we choose 200 au as our outer radius.