Radiation-Hydrodynamic Models of X-Ray & EUV Photoevaporating Protoplanetary Discs

Radiation-Hydrodynamic Models of X-Ray & EUV Photoevaporating Protoplanetary Discs

J. E. Owen, B. Ercolano, C. J. Clarke & R. D. Alexander
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Department of Physics and Astronomy, University College London, Gower Place, London, WC1E 6BT
Sterrewacht Leiden, Universiteit Leiden, Niels Bohrweg 2, 2300 RA Leiden, the Netherlands
E-mail: jo276@ast.cam.ac.uk

We present the first radiation-hydrodynamic model of a protoplanetary disc irradiated with an X-EUV spectrum. In a model where the total ionizing luminosity is divided equally between X-ray and EUV luminosity, we find a photoevaporation rate of 1.4 Myr, which is two orders of magnitude greater than the case of EUV photoevaporation alone. Thus it is clear that the X-rays are the dominant driving mechanism for photoevaporation. This can be understood inasmuch as X-rays are capable of penetrating much larger columns (cm) and can thus effect heating in denser regions and at larger radius than the EUV can. The radial extent of the launching region of the X-ray heated wind is 1-70AU compared with the pure EUV case where the launch region is concentrated around a few AU. When we couple our wind mass-loss rates with models for the disc’s viscous evolution, we find that, as in the pure EUV case, there is a photoevaporative switch, such that an inner hole develops at 1 AU at the point that the accretion rate in the disc drops below the wind mass loss rate. At this point, the remaining disc material is quickly removed in the final 15-20% of the disc’s lifetime. This is consistent with the yr transitional timescale estimated from observations of T-Tauri stars. We however note several key difference to previous EUV driven photoevaporation models. The two orders of magnitude higher photoevaporation rate is now consistent with the average accretion rate observed in young stars and will cut the disc off in its prime. Moreover, the extended mass-loss profile subjects the disc to a significant period (20% of the disc’s lifetime) of ‘photoevaporation starved accretion’. We also caution that although our mass-loss rates are high compared to some accretion rates observed in young stars, our model has a rather large X-ray luminosity of erg s; further modeling is required in order to investigate the evolutionary implications of the large observed spread of X-ray luminosities in T-Tauri stars.

accretion, accretion discs - circumstellar matter - planetary systems: protoplanetary discs - stars: pre-main-sequence - X-rays: stars.
pagerange: Radiation-Hydrodynamic Models of X-Ray & EUV Photoevaporating Protoplanetary DiscsReferencespubyear: 2002

1 Introduction

The structure and evolution of protoplanetary discs is currently an important and much debated topic in astrophysics. Detailed knowledge of disc evolution is needed to constrain the conditions under which planets form. In particular, an important issue in disc evolution is the mechanism by which the disc is eventually dispersed and the timescales involved, which consequently sets the time available for planet formation to occur.

For solar type stars it is now well established observationally that at an age of a few Myr a large percentage of stars are still surrounded by discs that are optically thick at infrared wavelengths (e.g. Haisch et al., 2001). However after yr most stars are observed to no longer be surrounded by discs. The relatively small number of discs that are observed to be in transition (between being surrounded by an optically thick disc and having no detectable disc) implies that the transition occurs on a much shorter yr timescale (Kenyon & Hartmann, 1995; Duvert et al., 2000). This ‘two timescale’ phenomenology provides a strong constraint on models of disc evolution and dispersal and is inconsistent with standard viscous-accretion models which predict a power-law decline (Hartmann et al., 1998)111 It is a matter of current debate whether this two-timescale description also applies to low mass stars (see Sicilia-Aguilar et al., 2008; Ercolano et al., 2009a; Currie et al., 2009)..

One of the most successful models in reproducing this ‘two-timescale’ behavior is photoevaporation by radiation from the central star. Recent detection of [Ne ii] emission lines (Pascucci, 2007; Herczeg et al., 2007; Najita et al., 2009; Flaccomio et al., 2009; Pascucci et al., 2009) has been interpreted as evidence for irradiation of the disc’s atmosphere by energetic photons (extreme ultraviolet or X-ray) and possibly associated photoevaporation (Alexander, 2008b; Pascucci et al., 2009). However there is still no consensus over whether photoevaporation is best explained in terms of far ultraviolet (FUV, eV), extreme ultraviolet (EUV, eV) or X-ray irradiation (eV), or some combination of the above. EUV photoevaporation (of optically thick discs) is theoretically well developed (Hollenbach et al., 1994; Font et al., 2004) with a predicted mass-loss rate of 10Myr. The mass-loss profile is dominated by a region within several AU of the characteristic (gravitational) radius AU where the internal energy of the gas exceeds the escape potential. Coupling this mass-loss profile to viscous evolution models (Clarke et al., 2001), this ‘EUV-switch’ model shows that once the accretion rate through the disc drops to the photoevaporative mass-loss rate, a hole develops and the inner disc drains onto the central star on the short viscous accretion timescale at a few AU. The inner region is now optically thin to EUV photons and a higher ‘direct-wind’ mass-loss rate of 10Myr is obtained (Alexander et al., 2006a). The disc then photoevaporates from the inside out on a timescale of 10yr (Alexander et al., 2006b). The relevance of the EUV-switch model, however, depends, albeit weakly, on the EUV luminosity of the star; the canonical mass loss rates quoted above correspond to an EUV luminosity of 10 photons s. This cannot be verified by direct observation owing to the strong absorption of ionizing photons by intervening neutral hydrogen. While Alexander et al. (2005) derived EUV fluxes from young solar type stars which are consistent with those required by photoevaporation, it has been pointed out (e.g. Glassgold et al., 2004; Ercolano et al., 2009b; Glassgold et al., 2009) that the EUV photons emitted by the star do not necessarily reach the disc, given that EUV photons are extremely easily absorbed, by small columns of neutral material surrounding the young star. X-ray observations of young T-Tauri stars show that neutral hydrogen columns of up to 10cm are common around many objects, probably due to accretion funnels, implying that EUV photoevaporation is unlikely to be effective until late times.

FUV driven neutral flows have also been considered as a possible mechanism for disc dispersal. In rich cluster environments, Johnstone et al. (1998) showed that FUV radiation from nearby stars can dominate the photoevaporation rate. This, however doesn’t explain the rapid dispersal of discs around young stars with no nearby source of FUV flux. Gorti & Hollenbach (2009) have considered FUV driven flows from the central star and argue that mass-loss rates of order Myr can be obtained at large radius (AU).

In contrast to EUV radiation, X-rays are able to penetrate larger columns of neutral material. Also X-rays are largely unaffected by accretion funnels (since they are emmitted in the magnetosphere, above the accretion funnels); furthermore the X-ray luminosity function of young T-Tauri star is well known (Preibisch et al., 2005; Albacete Colombo et al., 2007). For these reasons several authors have considered X-ray irradiation as a possible driving mechanism for photoevaporation (e.g. Alexander et al., 2004; Gorti & Hollenbach, 2009; Ercolano et al., 2008b, 2009b). While Alexander et al. (2004) and Gorti & Hollenbach (2009) quantitatively disagree on the mass-loss rates due to X-ray irradiation, they both conclude that X-ray heating is considerably less important than EUV heating. Ercolano et al. (2008b, 2009b) obtained mass-loss rates in the range Myr and argue that the X-ray heating will dominate the mass-loss and hence the evolution of the disc. While the hydrodynamics of EUV photoevaporation have been solved in detail, there are no prior hydrodynamical studies of X-ray irradiated gas. Alexander et al. (2004); Gorti & Hollenbach (2009) both use a 1+1D description of the radiative transfer, and Ercolano et al. (2008b, 2009b) uses a 3D approach; however, in all previous studies the disc structures are assumed to be in hydrostatic equilibrium and the mass loss-rates are obtained either by assuming a sonic flow from the location where the local gas temperature equals the escape temperature (the method, e.g. Alexander et al., 2004; Ercolano et al., 2008b, 2009b), or by coupling a spherical Parker wind solution to the hydrostatic structure (Gorti & Hollenbach, 2009). While these methods may (or may not) give a reasonable order of magnitude estimate of the total mass-loss rate, they are unlikely to give a reliable radial profile of the mass-loss rate, as is needed to determine the evolution of viscously evolving photoevaporating discs.

In this paper we couple the radiative transfer models of Ercolano et al. (2009b) to two-dimensional hydrodynamical models to obtain a self consistent flow solution for a primordial disc (i.e. an optically thick disc extending to the dust destruction radius) and for inner hole discs (i.e. gas discs with cleared inner holes). In section 2 we present our approach to the radiation-hydrodynamics. In sections 3 and 4 we apply these methods to a primordial disc and discs with inner holes and present the basic results. In section 5 we show the results of various tests we conducted. In section 6 we present the results of coupling our wind profiles to a viscous evolution model. In section 7 we compare our work to previous results and discuss the model’s limitations in section 8. In section 9 we discuss the observational consequences of our model and summarise our main findings in section 10.

2 Methods

In what follows we undertake 2D hydrodynamic calculations in which - under the assumption that the gas is in radiative equilibrium - we parametrize the gas temperature as a function of local variables. This parametrization is based on the results of 3D Monte Carlo radiative transfer calculations using mocassin (Ercolano et al., 2003, 2005, 2008a). We then verify a posteriori that the cooling timescale is less than the flow timescale and hence justify our assumption of radiative equilibrium. Below we set out the details of the radiative and hydrodynamic elements of the calculation.

2.1 Radiative Transfer

Ercolano et al. (2008b, 2009b) employed the three-dimensional Monte Carlo photoionization and dust radiative transfer code, mocassin (Ercolano et al., 2003, 2005, 2008a), to calculate the temperature and ionization structure of a typical T-Tauri disc irradiated by chromospheric X-ray and EUV flux. mocassin uses a stochastic approach to the transfer of radiation and allows for arbitrary geometries and density distributions. The version of the code used for these models was modified by Ercolano et al. (2008b, 2009b) and we refer to these articles for a detailed description of the code and model setup. We use the model labeled FS0H2Lx1 by Ercolano et al. (2009b) and summarise here only the basic input parameters.

2.1.1 Intial Disc Configuration

The initial density structure of a protoplanetary disc surrounding a 0.7 M star (T = 4000 K, R = 2.5 R) was taken from the D’Alessio et al. (2001) set and was chosen as the model that best fit the median SED in Taurus. This model is in vertical hydrostatic equilibrium in the case that the dust temperatures are set by irradiation by the optical radiation from the central star and under the assumption of full thermal coupling between the gas and dust. The total mass in the disc is 0.026 M with an outer radius of 500 AU.

This disc structure, which provides our starting density distribution, employs a bimodal dust distribution, where atmospheric dust follows the standard MRN model but with a very low dust to gas mass ratio (2.5 for graphite and 4.0 for astronomical silicates) and interior dust consists of larger grains with a size distribution still described by a power law of index -3.5, but with = 0.005 m and = 1 mm and a dust-to-gas mass ratio of 3.1 for graphite and 5.0 for astronomical silicates. The transition between atmospheric and interior dust occurs at a height of 0.1 times the midplane gas scale height.

High-energy dust absorption and scattering coefficients are calculated from the dielectric constants for graphite and silicates of Laor & Draine (1993), which extend to the X-ray domain. Spherical grains are assumed and we use the standard Mie scattering series expansion for complex refractive indices, , where is the scattering parameter (see Laor & Draine, 1993). For and we use Rayleigh-Gans theory (Bohren & Huffman, 1983), and for and we use the treatment specified by Laor & Draine (1993), which is based on geometric optics.

2.1.2 Method

This initial disc was then irradiated by a synthetic coronal spectrum of X-ray luminosity L = 210erg s (see Figure 1

Figure 1: Input spectrum used in mocassin calculation. See text and Ercolano et al. (2009b) for details.

) and total ionizing luminosity of L = 410erg s. This synthetic spectrum was a thermal spectrum generated by the plasma code of Kashyap & Drake (2000); the emission measure distribution is based on that derived for RS CVn type binaries by Sanz-Forcada et al. (2002) (which peaks at K) and fits to Chandra spectra of T Tauri stars by Maggio et al. (2007) ( which latter peaks at around logK). Our use of the input spectrum shown in Figure 1 implies that we assume a zero column of neutral material close to the source. The experiments described in Ercolano et al. (2009b) however show that screening columns of up to cm, which would only absorb the EUV component of the spectrum, have little effect on the resulting photoevaporation estimates; we thus expect our results to generalise to the case of moderate screening ( cm) since this has little effect on the X-rays in the keV range which dominate the disc heating (see below). As described in Ercolano et al. (2009b), the disc temperatures are updated in response to irradiation by the X-EUV spectrum and the density structure updated so as to attain a situation of hydrostatic equilibrium for this updated temperature profile. These temperature and density distributions are then updated iteratively until a situation of both thermal and hydrostatic equilibrium is attained. It is these converged hydrostatic models that form the basis of the parametrization described below, even though the hydrodynamical calculation evolves away from this initial static structure towards a steady flow solution. We however check retrospectively that our parametrization also applies to the evolved structure (see Section 5).

Under the assumption that the thermal and hydrodynamical timescales are decoupled222The thermal timescale in the flow is typically less than a year and thus more than an order of magnitude less than the hydrodynamic time scale. a coupling of the mocassin code with a hydrodynamical code is trivial in principle. It simply consists of feeding the density fields produced by the hydrodynamical code into mocassin to calculate gas temperatures that are then in turn fed back to the hydrodynamical code. While conceptually easy, this type of coupling is computationally prohibitive.

In order to render the calculation tractable, we exploit the fact that the wind will be optically thin to high energy photons. In the density and temperature regimes expected in a photoevaporative wind, radiative heating and ionization compete with collision-initiated cooling and recombination and hence the ionization equilibrium strongly depends on the ionization parameter (see Tarter et al., 1969). The ionization parameter can be conceptually considered as the ionizing energy available to a gas particle at a distance from the source and will determine the ionization state of the gas and hence the gas temperature at that point (e.g. Igea & Glassgold, 1999). We show in Figure 2 the mean relation between ionization parameter and temperature in the X-ray heated region of the hydrostatic models of Ercolano et al. (2009b) (upper panel), where we define the X-ray heated region as comprising gas for which the absorbing column to the central star is less than cm (roughly the penetration depth of 1 keV photons). The scatter of the mocassin grid cells about the relation shown is small ( dex) thus showing that the ionization parameter is indeed a good determinant of the gas temperature. This confirms the central importance of X-ray heating in determining the thermodynamics of the heated atmosphere (note that a tight relationship between ionization parameter and temperature is not expected in general in the case of heating dominated by EUV radiation). In Section 4 we will consider the case of photoevaporation from discs containing inner holes for which we need to improve our parametrization in the lower temperature regime. The lower panel of Figure 2 depicts the result of a similar calculation of a disc which is truncated at an inner radius of AU; we have confirmed that this provides also a good fit to models with different inner hole sizes and hence use this relation in all the hydrodynamic models in Section 4. Note the similarity of the relations for the identical ionization parameter range. This implies that it is the direct field that is dominating in both cases. This should be contrasted with the pure EUV case, where the diffuse field dominates for the optically thick primordial disc and the direct field dominates for discs with inner holes.

(a) Relation for Primordial Disc
(b) Relation for Inner Hole Discs
Figure 2: Plot of Temperature as a function of ionization parameter for FS0H2Lx1 model (Ercolano et al., 2009b) with numerical parametrization shown in red. The black asterisks show the mean temperature in the ionization parameter bins obtained from mocassin, with a scatter about the mean of dex.

Since these relations show that the gas temperature varies monotonically with ionization parameter in the Xray heated regime, we can thus use the ionization parameter (which depends only on the incident flux and local number density) in order to prescribe the gas temperature for regions lying at columns of cm from the central source.

For columns larger than cm Ercolano et al. (2009b) (ECD09) found that the gas and dust were thermally coupled. Hence we use the dust temperatures of D’Alessio et al. (2001) to fix the gas temperatures for columns greater than cm. In order to check that our temperature parametrization is appropriate for the final converged flow solution, we feed back this density field into mocassin and compare the temperatures thus obtained; these results are presented in §5.

2.2 Radiation-Hydrodynamic Methods

In order to accurately determine the mass-loss rates and the kinematic and morphological structure of the evaporative flow, we need to compute the temperature distribution in the flow and feed this into a hydrodynamic code. We have modified the publicly available zeus-2d hydrodynamic code (Stone & Norman, 1992a, b; Stone et al., 1992) to include the heating from X-ray and EUV radiation, by writing a new module that calculates the gas temperature using the methods discussed in §2.1. We then use this temperature to update the internal energy density of gas () in the wind via the ideal equation of state:


where is the Boltzmann constant, is the hydrogen mass and is the ratio of the specific heats. The mean molecular weight is fixed in the flow, to a value of . The internal energy is only required to update the pressure needed for evaluation of the momentum equation. The pressure is thus determined again from the ideal gas law, via:


and hence the procedure is independent of choice of . This thermal update is done during both the source step and transport step in order to minimise the time to a converged solution.

We assume azimuthal and mid-plane symmetry and since we expect a photoevaporative flow to be approximately spherical at large radius, we use a spherical grid. We neglect magnetic fields and self- gravity. The rotation option in zeus-2d is enabled, which introduces the necessary pseudo-forces arising from rotating frames. The inner and outer radial boundaries are set to outflow, while the angular boundaries are given the appropriate symmetry boundary conditions. We choose to use the standard van-Leer (second-order) interpolation scheme and the von-Neumann & Richtmyer parametrization for the artificial viscosity with .

3 Primordial Disc

In order to determine the flow structure for the disc for the majority of its standard viscously accreting phase we run a single hydrodynamic model. We expect the structure of the flow to be set by the temperature of the X-ray heated region and, to a lesser extent, by the temperature of the disc atmosphere which is controlled by optical heating of dust in the disc’s surface layers. The structure of the underlying disc is insensitive to this provided that it is very optically thick, both in the X-ray and optical regions. This is amply fulfilled by primordial (untruncated) discs over many orders of magnitude in accretion rate, thus justifying the fact that we use a single hydrodynamic model to evaluate the wind mass loss profile for the entire period of the disc’s evolution prior to gap opening (see Sections 4 and 6). A similar argument has previously been applied also to the EUV case, where the mass loss profile is unaffected by the total column in the disc, provided it is optically thick in the EUV (Hollenbach et al., 1994).

3.1 Numerical Method

The grid is chosen to be in the range and AU333Note throughout this paper we use and to distinguish between spherical and cylindrical polar co-ordinates respectfully. with 100 grid cells in the angular direction and 250 in the radial. Grid cells are logarithmically spaced in the radial direction, so that we have adequate resolution at small radius: in particular we ensure that we can resolve the onset of any flow and can resolve the scale height of the X-ray ‘dark’ disc at all radii. As initial conditions, the vertically hydrostatic density and temperature structures generated by Ercolano et al. (2009b) were interpolated onto the zeus-2d grids and assigned Keplerian rotation (suitably adjusted for the small radial pressure gradient). The disc structure was then allowed to evolve hydrodynamically using the modified code, with the gas temperatures being updated as described in §2.1 at every timestep, until a ‘steady-state444Where the term ‘steady-state’ strictly means quasi-steady-state or more formally .’ disc/wind system was obtained. It should be noted that the density structure of the X-ray ‘dark’ regions is also allowed to evolve hydrodynamically, i.e. it is not reset to the original density structure at every time-step (compared with previous EUV hydrodynamic simulations e.g. Alexander 2008b; Font et al. 2004 which - since they did not model the EUV ‘dark’ disc - reset the base density of the flow to its original value at each time-step). We emphasise that the mass loss over the duration of the simulation is much less than the initial disc mass.

The evaporative flow was found to have converged after 4-6 orbital timescales at the outer boundary (convergence was also checked by determining the Jacobi potential was conserved along streamlines: See §5). The flow structure after 10 orbital timescales at the outer boundary is used in our analysis.

3.2 Results

(a) Initial Disc
(b) Flow Solution
Figure 3: The density structure of a primordial disc; with the initial disc shown in (a) and the ’steady-state’ flow solution in (b). The red line indicates the =2/3 surface from the central star in both plots. In (b); streamlines are plotted at 5% intervals of the integrated mass loss rate, with the sonic surface plotted in magenta. Note the smoothness of the sonic surface indicates that we are extremely close to a true ‘steady-state’. The snapshot shown is at a simulation time of 10 orbital periods at 100AU.

In Figure 3 we show the density structure of the wind and underlying primordial disc with streamlines overlain. The Figure shows that the flow originates near the point at which the X-ray flux can penetrate, i.e. a column of cm. This also indicates that the flow launching surface is significatly higher than the disc’s photosphere as seen by the star (red solid line in 3(b), where we have used the opacities of D’Alessio et al. 2001 to determine this surface). Also, comparing the intial disc structure (Figure 3(a)) and the final flow solution (Figure 3(b)) it is clear that the optically thick disc region is very similar and is thus observationally consistent with the SED of classical T-tauri stars by construction, with the surface of optical depth to stellar photons of lying at 0.1. We stress that although the X-ray heated wind launches from well above the disc mid-plane, both the wind and the upper regions of the underlying disc are optically thin to stellar photons and thus the fraction of young stars that are optically obscured by their discs remains at the level implied by the D’Alessio et al. (2001) models.

Figure 3 also shows that the mass loss is dominated from the region 5-20 AU, although there is still significant mass-loss at larger radii. We note that the streamlines are close to radial at large radius similar to those obtained from pure EUV models (Font et al., 2004; Alexander et al., 2006a; Alexander, 2008b). The observation that the flow is approximately radial at large radius allows us to be confident that the outer boundary has not affected our results, since radial supersonic outflow is treated exactly in zeus-2d (Stone & Norman, 1992a). The total mass loss rate is found by integrating the mass flux across a spherical surface situated at a radius of 85% of the outer grid boundary (sufficiently far from the outer radial boundary to avoid any effects from spurious reflections due to non-radial outflow); and X-ray ‘dark’ gas is explicitly excluded from our measurement. We then determine the cylindrical radius from within which the integrated mass loss rate occurs by following the streamlines down to the base of the flow. Figure 4 shows the cumulative mass-loss rate as a function of cylindrical radius, and shows that the mass-loss rate is dominated by the region between 5 and 20 AU, but significant mass loss occurs also at larger radii with 50% of the mass loss occuring between 18-70AU. This should be compared to the pure EUV case where only 10% of the total mass-loss occurs outside 18AU (Font et al., 2004). The X-rays cause a much broader wind profile by being able to heat gas at larger radii than pure EUV irradiation due to their greater penetrative power. The total integrated mass-loss rate over the entire disc was found to be 1.4yr.

Figure 4: Cumulative mass loss rate as a function of cylindrical radius, with the numerical fit used in the viscous evolution shown as the solid line. Note the smooth profile with only one ‘bump’, indicates the flow is extremely close to a true ‘steady-state’.

4 Inner Hole Flows

As will be shown in §6.2, photoevaporation combined with viscous evolution results in a gap opening in the primordial disc; in contrast to the case of planet clearing, where streams cross the gap at certain azimuths, in the present (axisymmetric) case the inner disc is starved of re-supply once the gap is opened. The inner disc then drains onto the star leaving the outer disc open to direct irradiation by the central star. The flow expected from such a system is of course different from a primordial disc’s photoevaporative flow, and in this section we set out the methods used to calculated mass-loss rates from discs with inner holes.

4.1 Numerical Method

Given we expect the mass-loss rates to vary with inner hole size as seen in the EUV case (Alexander et al., 2006a), we have run several models with inner holes at different radii, out to a radius beyond which disc holes become hard to detect observationally (i.e. at around AU). We used a regularly space spherical grid in order to maximise resolution while minimising computational effort. The grid was constructed with 400 cells in the radial direction and 200 cells in the angular direction, again sufficient to resolve the onset of any flow and the scale height of the ‘dark’ disc at all radii. All other zeus-2d parameters were kept the same as above. Table 1 summarises the input parameters used. As initial condition, the hydrodynamic primordial disc model was cut at the inner hole radius. Then the density profile was smoothed over several pressure scale lengths at inner edge, the we adjusted the angular velocity to match the new imposed radial pressure gradient. The model was then allowed to evolve hydrodynamically until it attained a steady flow solution.

Model (AU) Grid range (AU)
A 8.3 [1.33,90]
B 9.7 [1.67,90]
C 11.7 [3.35,100]
D 14.2 [3.35,100]
E 17.7 [3.35,100]
F 21.1 [3.35,100]
G 30.5 [3.35,100]
H 64.5 [10,450]
I 86.9 [15,500]
Table 1: Table listing the input parameters for simulations of inner hole discs, where is the inner edge of the X-ray ‘dark’ disc

4.2 Analytic Model

Before we perform detailed hydrodynamic simulations to determine the flow structure of transition discs it is useful to develop a theoretical framework in which to analyze the simulations. In order to estimate the mass-loss rate scaling we adopt a similar method to Alexander et al. (2006a) and assume the mass loss rate per unit area is where is the local launch velocity of the ionized gas. However unlike the calculations of Alexander et al. (2006a) where the gas was assumed to be isothermal, we have no such restriction in our calculations. We expect streamlines that are launched from the disc’s inner rim (which will turn out to dominate the mass loss) to start on a nearly radial inward trajectory (see Figure 4.3). Inevitably such a trajectory will result in a centrifugally driven wind (whose force term scales as ) resulting in a launch velocity that scales as the escape velocity . Furthermore, if we assume that the effective emitting area of the wind scales as the square of the inner disc radius () then we can write the wind mass-loss rate as scaling with . We will use the code to determine the number density at the base of the flow at the disc’s inner rim () and check out our estimate that in this case


The actual scaling of with however needs to be determined numerically since the ionization parameter at the flow base (and hence ) depends on the temperature at the inner rim and this needs to be calculated using the radiation-hydrodynamic algorithm for each hole size. We parametrize the dependence of mass loss rate on in Section 4.3.

4.3 Results

All inner hole flows share the same general morphology: the mass-loss is dominated by regions close to the inner edge, but significant mass-loss still occurs at larger radii. Figure 5

Figure 5: The ‘steady-state’ density structure of a transition disc and its photoevaporative wind with an inner hole of 14.2 AU (model D), the streamlines are plotted at 5% intervals of the integrated mass loss rate, with the sonic surface plotted in red. Note the sonic surface being close to smooth indicates that we are extremly close to a true ‘steady-state’. The snapshot shown is at a simulation time of 10 orbital periods at 100AU.

shows the flow structure of model D with an inner hole at 14.2AU. We extract mass-loss rates in an identical fashion to the primordial disc case and note again that the radial structure of the streamlines at large radii allows us to be confident that the outer boundary has not affected our mass-loss rates.

Figure 6 shows the ratio reproduces the predicted scaling discussed at the end of §4.2, i.e. it is consistent with our assumption that the effective launch area scales as and that the launch velocity scales as the escape velocity.

Figure 6: Plot showing analytic model () agrees well with the simulations. Error bars shown indicate errors in obtained by averaging over cells at the inner edge.

In Figure 7

Figure 7: Plot showing at different inner hole sizes; The blue point represent mass-loss rates from models A-I, while the black line shows the numerical fit used in the viscous evolution.

we present the mass-loss rates as a function of inner hole radius complete with the numerical fit used in the viscous evolution. This shows that for out to AU the mass-loss rates fall approximately as while at larger the mass-loss rates begin to increase. We can understand this behavior by noting that equation (3) and the definition of the ionization parameter implies that so that the scaling of photoevaporation rate with radius depends on the radial variation of and hence on the form of the temperature versus ionization parameter relation. From Figure 2(b) we see that this contains a point of inflection at , which is encountered when the is around AU. At smaller radii and hence higher temperatures the variation in ionization parameter is small for a given temperature variation and hence the modest decline in the flow temperature with inner hole radius is achieved at nearly constant . Beyond AU, a modest decline in temperature corresponds to a larger decline in and hence the fall off in becomes less steep and even starts to increase somewhat with increasing . It should however be stressed that the radial variation of with is very mild (see Figure (7)).

5 Tests

In order to be confident that the solutions obtained are accurate and are unaffected by the assumptions discussed in §2.2, we have conducted various tests which we discuss below:

  1. A numerical convergence test using double the resolution confirms the flow is converged to 5%.

  2. In order to be confident the parametrization of X-ray heating is correctly reproducing the temperature structure in the flow we have re-calculated the gas temperatures of the steady-state solution using mocassin. Figure 8 shows a comparison between the temperatures determined by the zeus-2d algorithm and those obtained using mocassin. This shows that we have errors 15% for 50% of all grid cells and errors 40% for 95% of the grid. In Figure 9 we show the comparison between the gas temperatures determined using the zeus-2d algorithm and those obtained using mocassin at various cuts in cylindrical radii. It shows that we have good agreement throughout the grid especially in the launch region of the flow (region of large temperature gradient) and in the regions of the flow which are subsonic (red regions) where the thermal pressure gradient can be a first order effect. We note that the disagreement between the temperatures in the flow/disc transition region is an artifact of mocassin which imposes a hard transition between the molecular and atomic zones (see discussion in ECD09), hence the errors quoted can be considered an overestimate.

  3. It has been assumed that the gas and dust are thermodynamically coupled in regions that are optically thick to X-ray irradiation. In order to test this assumption we have performed a run where we allow the gas in the X-ray dark regions to evolve adiabatically. We then compare the timescale for the gas and dust to return to thermal equilibrium through collisions (Whitworth & Clarke, 1997) to the local dynamical timescale and we find that the dynamical timescale is the thermal timescale at the minimum. Thus we can conclude that the gas and dust will be thermodynamically coupled when not penetrated by the X-ray flux.

  4. In the following sections, we assume that the mass-loss profile obtained for the primordial disc is valid for the majority of the discs lifetime (i.e. the viscously accreting phase). During this period the disc structure will change as material drains onto the central star. In order to check our mass-loss rate is is valid for a varying (optically thick) disc structure, we compare it to the mass-loss rates obtained for the adiabatic runs mention above, which produce markedly different (adiabatic) disc structures and find good agreement between the models.

  5. In our hydrodynamic modeling we have used the dust temperatures (and hence gas temperature of the X-ray ‘dark’ region) obtained by D’Alessio et al. (2001). Since the dust temperatures will be linked to the disc structure and vice-versa, and since the disc structure has evolved from the original hydrostatic equilibrium structure, is important to assess how accurately the dust temperatures of the D’Alessio et al. (2001) model can describe the dust temperatures in the hydrodynamic solution. We’ve calculated the dust temperatures of both the original model and the flow solution using mocassin and find we have agreement to within the Monte-Carlo errors.

  6. In order to determine whether our models are close to a true steady-state we compute the Jacobi555The Jacobi potential is the analogue to the Bernoulli potential in a rotating frame. potential along the streamlines, noting that this will be constant in a steady flow. For each streamline we determine the enthalpy term by paramatrizing the run of pressure and density as a baratropic relation (i.e. ) along the streamline, then integrating along the streamline. However we note that, the flow is not strictly baratropic or irrotational and thus don’t expect the Jacobi potential to be constant for all streamlines. Figure 10 shows that the Jacobi constant is conserved along the streamlines. We find that the Jacobi constant is conserved to within 3-4% indicating we are very close to a true steady-state flow.

  7. ECD09 models do not include molecular cooling in the thermal balance. While (as shown in Figure 3(b)) the flow is launched in a region higher than the optically thick disc, the lack of molecular cooling could in principle be a potential source of error on the gas temperature in the disc interiors, and as such the cause of a hydrostatically ‘puffed up’ disc. Although as shown by Glassgold et al. (2004, their Figure 4) molecular cooling only comes in at large column densities and, even there, it only has a very limited role when compared to dust-gas collisions, it is useful to test whether a cooler (flatter) disc interior would result in significant changes in the photoevaporative mass loss rate of the system (due to neglecting molecular cooling, or some other, as yet unidentified cooling mechanism). We can do this by considering the mass-loss rates obtained from the adiabatic models described above (point iii & iv). The use of an adiabatic equation of state results in a cooler disc structure and hence a lower launching height. While the exact quantitative result is sensitive to disc heating, the overall mass-loss rates in all models were found to be M yr. This effect is small since the change in gravitational energy between the midplane and the launching height is of order 10% of the total gravitational energy gained by the wind, as shown in Figure 14.

Figure 8: Histogram of cells with temperature errors less than a given value. Black solid line refers to errors in cells in entire grid, red dashed line to errors from cells that are locally sub-sonic (in the rotating frame).
Figure 9: Radial cuts showing temperature comparison for R=2,10,20AU. Solid line indicates temperatures calculated by the zeus-2d algorithm, where the red solid line indicates where those cells are locally subsonic (in the rotating frame) and in the photoevaporative flow. Blue stars represent the temperatures determined by mocassin.
Figure 10: Plot of the Jacobi potential along various streamlines in the range 2-50AU, plotted against the normalised streamline length. The variation in the mean value along each streamline is 3%

6 Viscous Evolution

The evolution and dispersal of an irradiated protoplanetary disc is controlled by the simultaneous actions of the photoevaporative wind and viscous torques in the disc. A numerical simulation of disc dispersal must therefore simultaneously account for photoevaporation and viscous evolution.

6.1 Method

The diffusion equation for the evolution of the surface density is (Lynden-Bell & Pringle, 1974; Pringle, 1981; Alexander et al., 2006b):


Given that the initial disc structure taken from D’Alessio et al. (1998) has a surface density that scales approximately as , an initial surface density structure of the form


is used, where is the initial disc mass and is the initial scale size of the disc. This initial structure is a similarity solution of the viscous diffusion equation (Lynden-Bell & Pringle, 1974) when the viscosity scales linearly with radius and can be parametrized as


and the mass-loss in the wind can be neglected. The viscous scaling constant can be determined from (see Hartmann et al., 1998; Clarke et al., 2001; Alexander et al., 2006b):


The integration of equation (4) is performed using a first-order explicit scheme, the grid consists of 1000 points equispaced in and spans the range [0.0025, 2500] AU with zero-torque boundary conditions (Pringle et al., 1986). The radius beyond which of the disc mass is initially located () is taken to be 10AU and the initial accretion rate is yr. While we adopt a viscosity scaling in order to be consistent with the surface density profile of our initial disc structure, Clarke et al. (2001) & Alexander et al. (2006b) considered the effect of using a different viscosity parametrization and concluded that, while the absolute timings would change, the overall qualitative nature of a photoevaporative ‘switch’ was unchanged.

The relative timescales between the three stages (viscously evolving, gap opening and disc clearing) are independent of the initial disc mass (for a fixed initial accretion rate) and can thus be scaled accordingly for any chosen value. In order to set a reasonable absolute value for the disc lifetime, we have chosen an initial disc mass of 0.2M, which (given the initial accretion rate) corresponds to an initial viscous timescale at the disc’s outer edge of yr; this in turn corresponds, for our chosen case of A.U. to a viscous alpha parameter of at .

6.1.1 Surface Mass-loss rates

The functional forms of used are those obtained for the primordial disc models and inner hole models. Since we expect the inner hole to drain on a short timescale, as a first approximation we either use the primordial disc wind profile, or the inner hole mass-loss profile and do not attempt to model the transition in any detail. The transition occurs once the radial column density to the inner edge is cm.

6.2 Results

Figure 11: Snapshots of the surface density plotted at 0%, 20%, 40%, 60%, 70%, 75%, 77%, 79%…99% of the disc’s total lifetime (see footnote 4), The surface density is in dimensionless form and in units of (see equation 5).

Figure 11 shows the evolution of the surface density of the disc. At early times the evolution of the disc is similar to a standard viscous accretion disc with a power-law decline in surface density and accretion rate. Once the accretion rate is comparable to the mass-loss rate due to the wind (which occurs at of the disc lifetime666At this point we define the disc lifetime as being the total lifetime when photoevaporative clearing is taken into account; i.e. the time to clear out to 87AU (beyond which is not accurately known); Due to the power law nature of viscous draining, it is hard to define a ‘lifetime’ of a control disc (with no photoevaporation). However, to be indicative: the model disc that we have run has a total lifetime when photoevaporation is included of Myr; if such a disc were simply viscously evolving without photoevaporation it would become optically thin in the region of the disc that dominates the K band emission after Myr.), the photoevaporative wind begins to affect the surface density profile in the range 5-30 AU; a gap then opens at 1AU after 78% of the disc lifetime. Once the gap opens, the inner disc quickly drains on a timescale of order the viscous timescale at 1 AU. During this draining time the inner edge of the outer disc is continually photoevaporated outward due to the significant mass-loss at larger radius.

Once the inner disc has drained, the outer disc is exposed to the direct photoevaporative wind, and the inner edge of the disc has now moved to 9AU. The total mass remaining in the disc after the inner disc has drained onto the star was 9% of the intial disc mass and the outer disc is then dispersed entirely in the remaining 15% of the disc’s lifetime. We find, in total 84% of the intial disc mass is accreted onto the central star, with the remaining being lost in the wind.

In Figure 12

Figure 12: Plot of accretion rate onto the star against time. The black solid line shows the model for the X-EUV wind, while the blue dashed line shows the pure EUV canonical model (see Alexander et al., 2006b, Figure 1).

we show the accretion rate onto the star throughout the disc’s lifetime and show that the disc evolves much faster compared to a standard viscously accreting disc. We also see that, even prior to gap opening, the disc spends a significant portion of its lifetime () with the accretion rate onto the star being less than the photoevaporative wind rate: the accretion rate onto the star is 1 yr at the point that the gap opens, which is about an order of magnitude less than the wind mass loss rate. This effect (which has been termed ‘photoevaporation starved accretion’ by Drake et al. 2009) is a consequence of the radially extended mass loss in X-ray heated discs so that, prior to gap opening, there is a significant region of the disc where the accretion rate decreases towards small radii. This behavior was not seen to be a prenounced effect in previous pure EUV models since, in the EUV case, the mass loss is more strongly localised and the gap opens in a disc that is hardly depleted by the wind at large radii.

7 Comparison with previous work

Several other authors have considered photoevaporative flow driven by either, EUV or X-ray radiation and Ercolano et al. (2009b) have considered the combined effect of X-EUV radiation. The predicted mass loss rates due to pure EUV photoevaporation have been accurately determined as Myr (Font et al., 2004) for a 0.7M star and an ionizing flux of s. Also Alexander et al. (2006a) found that when the inner hole first forms, the mass loss rate rises by about a factor ten and then increases mildly as the hole grows according to .

As discussed in Section 1, previous estimates of X-ray photoevaporation rates to be found in the literature are qualitatively different from the values presented in the present work since they are not based on hydrodynamic calculations. Alexander et al. (2004) estimated mass-loss rates per unit area at AU which were similar to those obtained in this work. ECD09 argued for similar surface mass-loss rates as those of Alexander et al. (2004), but estimated a much larger launching regions (i.e.over a region 5-50AU). They therefore concluded that X-ray photoevaporation would dominate over EUV photoevaporation at all but the lowest of X-ray luminosities. The estimates of Gorti & Hollenbach (2009) suggested much lower mass loss rates than the EUV case, implying that X-rays would only have an indirect effect on disc dispersal. ECD09 include a detailed discussion of the differences in the thermochemical calculations of Gorti and Hollenbach 2009, which may have led to their result, and concluded that the rather harder X-ray spectrum employed by Gorti and Hollenbach may have been an important contributory factor.

The mass loss rate of 1.4Myr that we obtain from our radiation hydrodynamical calculations can be compared to the estimate of ECD09 for the same input parameters (model FS0H2Lx1 of ECD09) of 6.7 Myr. This difference arises from the fact that ECD09 assumed that the gas escaped at the sound speed at the point at which the internal energy of the gas was equivalent to the escape potential. In reality the flow is launched sub-sonically from deeper in the disc where the density is higher. It is nevertheless reassuring that the two values are of a similar order of magnitude. We also find that the mass-loss occurs over a large radial range of 1-70AU, in contrast to the EUV case where it is steeply peaked about one characteristic radius (for a 0.7M star this would be 6.2AU). Given the mass loss rates obtained and the monotonic temperature : ionization parameter relations found (see Figure 2), it is clear in this case that X-ray photoevaporation will be a dominant effect in disc evolution.

Coupling the mass-loss profiles to secular viscous evolutions models has only been previously attempted for the pure EUV case by Alexander et al. (2006b) where the mass-loss profile had been accurately determined. In order to compare our model with what the EUV model of Alexander et al. (2006b) would predict with the same input parameters, we take the canonical model (see Alexander et al., 2006b, Figure 1) and multiply both the initial disc mass and the initial disc viscous timescale by a factor , thus ensuring that the rescaled model has the same initial accretion rate, disc mass and viscous timescale as the present model. 777Ignorance of the magnitude of the viscosity in discs means that the absolute values for disc lifetimes that we quote here in the EUV and X-ray heated case should not be used as grounds for favoring one model or another: we have scaled the model used here so that its total lifetime is an observationally reasonable 1Myr and the EUV lifetime is then rather long, but we could have just as easily scaled down the initial mass and initial viscous timescale in such a way that the EUV lifetime was observationally reasonable and the X-ray lifetime was unreasonably short. The purpose of our present comparison is therefore to look at the ratio of timescales in X-ray and EUV models with the same underlying viscous parameters. In the (rescaled) EUV case, it was found that the disc dispersed entirely after 24Myr with the transition stage occurring on a short yr timescale (approximately 3% of the total disc lifetime). Our results are consistent with this ‘switch’ nature of disc dispersal (albeit in somewhat softened form: see Figure 12). In the present (X-EUV) models we find that a gap in the disc opens at 1AU after 78% of the discs lifetime (i.e. nominally at an age of 1Myr). Once the inner edge of the outer disc is exposed to direct radiation we find the outer disc is quickly removed on a timescale yr with a mass loss rate scaling as .

While X-ray driven photoevaporation produces the same qualitative results when it comes to disc dispersal, there are several distinct differences. Prior to gap opening the surface density begins to flatten, indicating that the accretion rate is dropping below that of a steady disc (i.e. that the disc is in a state of ‘photoevaporation starved accretion’ as discussed above). This effect is much less obvious in the EUV case owing to the more restricted radial extent of the EUV wind. Once the gap opens at around a few AU (similarly in both cases), the inner edge of the outer disc evolves rapidly to AU in the X-ray case due to the fact that this region was already significantly depleted, due to the broad mass-loss profile. In the EUV model, the inner edge remains at 1AU during disc draining, due to the steep mass-loss profile and then quickly moves outwards to AU as the photoevaporation rate rises by an order of magnitude when the inner edge of the disc is directly irradiated. The difference in behavior of the edge of outer disc during inner disc draining, can be readily understood when the surface mass-loss profiles are compared for the EUV case and X-EUV case. Figure 13 clearly shows the X-EUV profile is much broader than the EUV profile; the EUV wind thus has very little effect on the disc beyond whereas the X-EUV wind subjects the disc to a period of depleted accretion inwards of 25AU prior to gap opening.

Figure 13: In this Figure we plot the surface mass-loss rate for the X-EUV wind (solid black line) obtained in this work and the EUV wind (dashed blue line) (Font et al., 2004; Alexander, 2008b). We have scaled the peaks of the profile to be the same value in order to compare the relative broadness of the different profiles. Clearly this shows the X-EUV profile is much broader than the EUV profile.

During the clearing of the outer disc, the mass-loss rate falls as (out to 30-40AU) in the X-ray case, compared to increasing as in the EUV case. This difference is due to the fact that in the EUV case, where the ionized region is approximately isothermal, the flow is launched at roughly constant ( sonic) velocity and a simple Strœmgren sphere argument implies that the density in the launch zone scales as ; in the X-ray case, where a range of temperatures are present, the effective launch speed from the exposed inner rim of the disc is the local escape speed which scales as . The density in the launch zone (which corresponds to a nearly constant value of the ionization parameter) falls off slightly less steeply than . Both these effects mean that EUV photoevaporation is relatively favored at large radius. We estimate that the radius at which the direct EUV photoevaporation rate becomes comparable to the X-ray photoevaporation rate is AU and hence the EUV will remain unimportant in the case of discs around solar type stars with X-ray luminosities 10erg s and greater. It may however be relevant for systems with lower X-ray luminosities. Given the large scatter of X-ray luminosities for stars of similar masses, a full discussion of the parameter space is beyond the scope of this paper, but will be considered in future work.

8 Model limitations

While this model is an important step forward in understanding disc dispersal, it is important to discuss the limitations of our models. We have neglected the effects of non-ionizing FUV radiation in our calculations. Gorti & Hollenbach (2009) have shown that FUV radiation is capable of heating gas to above the escape temperature at large radius (100-200AU), and predict mass loss rates of a similar magnitude to those obtained from our X-EUV models. Given that these mass-loss rates are far more uncertain than those obtained from a self-consistent hydrodynamic model we cannot make any definite predictions of the effect FUV inclusion would have on our model. Qualitatively, however, the existence of a further FUV wind at large radius would starve the inner disc of material and would thus accelerate the final inside-out clearing described here.

Throughout our modeling we assume that the gas is always in thermal equilibrium with the radiation based on the comparisons between the thermal and flow timescales discussed in §2.1. This assumption can be tested by comparing the gained energy flux of the wind to the ionizing luminosity. Figure 14 shows a plot of specific energy for the primordial disc. Considering the streamlines around the median mass-loss value (launching from 15-20AU) the total change in specific energy along these streamlines from launch (with a negative specific energy of erg g) to exit from the grid (with a positive specific energy of erg g) is erg g. Given the grid boundary is far from the sonic surface (see Figure 3(b)) this can be considered close to the maximum gained energy flux of the wind (c.f. Parker wind; Parker 1958). In order to compare with the input energy rate (), we can compute the mechanical luminosity of the wind using (where is the specific energy) which evaluates to erg s. Thus it is clear that the energy gained in the wind is much less than the input energy from the X-EUV spectrum ( erg s). For all flows (both the primordial and inner-hole flows) we explicitly calculate the gained energy flux in the wind and find that its (8%). Hence the assumption of thermal equilibrium is valid in this case. Hence the terms in the energy equation from advection and work can be neglected.

Figure 14: Figure showing the specific energy of gas particles throughout the grid. Streamlines are plotted at 5% of the total mass-loss rate. The dark blue lane indicates the transition between a overall negative specific energy (bound gas) and an overall positive specific energy (unbound gas). This clearly shows that the total energy gained by gas launched into the wind is less than with a total efficiency compared to the input X-ray energy.

During the transition between a primordial disc and a transition disc, we assume the mass-loss profile for the former, until the column to the inner edge of the outer disc is less than cm. In reality we would expect the transition to be a smooth process. Given the difference in mass-loss rates between the primordial and transition disc models is small, and the transition time is rapid, such an approximation will have little effect on the overall evolution. However a more detailed approach will be required when it comes to predicting the statistical distribution of inner hole sizes as a function of stellar parameters (Alexander et al., 2006b).

9 Observational Consequences

We have shown that the X-EUV photoevaporation model described here can reproduce the ‘two-timescale’ evolution observed in young solar type stars, in agreement with current observational constraints on relative timescales. In this sense it is qualitatively similar to the EUV switch model as set out in Alexander et al. (2006b): in particular the second (rapid clearing) timescale is similar and so is the radius ( AU) at which the gap first opens. The temperature distribution in the flow is however rather different in the two cases (the X-EUV wind is somewhat cooler: compare Figure 9 with the EUV case where the wind is close to K throughout). It is thus an immediate task to compare the predictions of the new model with observational diagnostics that have been studied in the EUV case (Ercolano et al in prep.). In particular, the strength and line profile of the [Ne II] emission at m (Pascucci, 2007; Pascucci et al., 2009; Herczeg et al., 2007; Najita et al., 2009; Flaccomio et al., 2009) has been suggested as an observational test for X-ray irradiated discs e.g.(Glassgold et al., 2007; Meijerink et al., 2008; Güdel, 2008; Schisano et al., 2009); our calculations open up the possibility of detailed modeling, analogous to that of Alexander (2008b) for the EUV only case. Likewise, the X-EUV model needs to be assessed with regard to the production of forbidden line emission (including O[I]) as in the study of Font et al. (2004) for the EUV case.

The important difference in terms of disc secular evolution between the X-EUV and the EUV only model is the magnitude of the photoevaporative mass loss rate which is larger in the X-EUV case by 2 orders of magnitude. The wind mass loss rate is now comparable with typical accretion rates measured in Classical T Tauri stars (Natta et al., 2006). This means that whereas the EUV switch provided a mechanism for the final shutdown of a protostellar accretion disc when its accretion rate had fallen to very low (nearly undetectable) values, the X-EUV wind (at least for the rather high X-ray luminosity employed here) should cut off a disc in its prime. This raises a number of issues:

9.1 Accretion rates in T Tauri stars and their relationship with X-ray luminosity

If all discs were subject to the vigorous X-EUV winds found in this model, then young stars should rarely exhibit accretion rates less than a few times yr, since subsequent draining is rapid (see Figure 12). This is clearly untrue observationally, as this accretion rate is close to the median for Classical T Tauri stars. However it needs to be borne in mind that T Tauri stars exhibit a orders of magnitude scatter in in X-ray luminosities at a given mass (Preibisch et al., 2005) and that the value that we have adopted here is modestly above the median value. We have not yet undertaken full hydrodynamical calculations with a range of X-ray luminosities but the hydrostatic estimates of ECD09 tentatively suggest that the photoevaporation rate may scale roughly linearly with X-ray luminosity (cf the EUV case which, being recombination limited, predicts a mass loss rate that scales as the square root of the EUV output of the star; Hollenbach et al. 1994). Therefore we would expect that in stars with a lower X-ray luminosity the plot of accretion rate versus time (equivalent to Figure 12) would steepen at a lower accretion rate.

What is the expected sign of the dependence of accretion rate upon X-ray luminosity within the subset of stars that have not yet cleared their discs (i.e. in Classical T Tauri stars)? Some care is needed in phrasing this question. If we have a coeval cohort of stars, then evidently the stars that are stronger X-ray emitters would be expected to exhibit lower accretion rates, as a consequence of ‘photoevaporation starved accretion’. On the other hand, if we have a mixture of (disc bearing) stars of all ages then in a statistical sense the strong X-ray emitters will be younger (because the predicted disc lifetime is less) and they will have higher accretion rates on average (because the low luminosity X-ray sources live for a longer time and spend most of this time with accretion rates modestly above the (low) wind rate).

A clear observational consequence of the lower disc lifetime in luminous X-ray sources is that such sources spend a higher fraction of the pre-main sequence state in a disc-less (Weak Line T Tauri) state: consequently one expects that Weak Line T Tauri stars should be more luminous X-ray sources on average than their disc bearing (Classical T Tauri) counterparts. This is indeed the case (e.g. Neuhäuser et al., 1995; Flaccomio et al., 2003; Preibisch et al., 2005) though this has also been variously explained in terms either of X-ray absorption in accretion columns (Gregory et al., 2007) or confinement of the X-ray producing corona in accreting systems (Preibisch et al., 2005) (in other words, it is conventionally assumed that accretion suppresses X-ray emission rather than the effect predicted here whereby X-ray emission suppresses accretion: see also Drake et al. 2009 for a discussion of this possibility). Naturally, this effect cannot be quantified until we have verified through hydrodynamic simulations how the wind mass loss rate scales with X-ray luminosity.

9.2 Disc lifetimes

The fact that the X-EUV photoevaporation rate is much higher than that for EUV only has obvious implications for disc lifetimes. For example, we find that for the same initial parameters, the disc lifetime is around 1.5 orders of magnitude lower in the present case than when compared with the EUV only calculations of Alexander et al. (2006b). However it should be borne in mind that in both cases the absolute lifetime depends on the the time for viscous evolution down to a critical accretion rate and this therefore depends both on the initial disc mass and the viscosity in the disc. The magnitude of the latter can in any case be adjusted so as to reproduce any reasonable disc lifetime. What will be interesting, however, is to examine the distribution of relative disc lifetimes that results for a plausible distribution of disc initial properties and X-ray luminosities (c.f. Armitage et al., 2003).

9.3 Transition discs (inner hole sources)

We have shown that X-ray luminous T Tauri stars should evolve through the following sequence of states: a) optically thick (untruncated) stage dominated by viscous evolution b) inner gap phase, with a low density inner disc and accretion on to the central star that is less than the wind mass loss rate ( yr) and c) inner hole phase with an evacuated inner disc and essentially no accretion onto the star and an outer disc mass of . This is qualitatively similar to the EUV case described in Alexander et al. (2006b); important quantitative differences are the () higher upper limit on the accretion rate in stage b) in the X-EUV case and also the somewhat higher upper limit (roughly an order of magnitude) on the outer disc mass during stage c). It has often been argued that the relatively high accretion rates and outer disc masses in observed inner hole sources is a counter-argument against the relevance of (EUV) photoevaporation (Najita et al., 2007; Alexander, 2008a; Kim et al., 2009). We thus see that X-EUV photoevaporation alleviates both of these problems.

It would however be unwarranted to necessarily conclude that all (or even the majority) of observed T Tauri stars with spectral evidence for cleared inner holes are created in this way. Indeed it should be stressed that even if we judge that the bulk of such sources are created by other mechanisms, this does not diminish the important role of photoevaporative clearing in the majority of discs. It needs to be recalled that the evolutionary sequence described above implies - in the absence of other intervening processes (such as the formation of planets or grain growth) - that all disc systems should evolve through a short-lived inner hole phase. The observed frequency of inner hole sources (of order of all disc-bearing sources) is broadly compatible either with a scenario in which all observed inner hole sources are a short-lived evolutionary phase undergone by all stars or one in which the majority of observed inner hole sources are instead a long-lived evolutionary phase undergone by a subset of all T Tauri stars. There is no shortage of other mechanisms proposed for inner hole creation, including for example clearing by planets or binary companions (Quillen et al., 2004; Varnière et al., 2006; Rice et al., 2003, 2006), grain growth in the inner disc (Dullemond & Dominik, 2005; Ciesla, 2007), photophoresis (Krauss et al., 2007) and inside out evacuation by the magnetoroational instability (Chiang & Murray-Clay, 2007) (see e.g. reviews by Alexander 2008a and discussion in Kim et al. 2009); what is notable is that these alternative mechanisms are not linked to a mechanism for clearing the residual outer disc and are thus not associated with a short-lived phase prior to complete disc dispersal. Therefore although they may apply to some of the observed inner hole sources they cannot represent an evolutionary norm for young stars without over-predicting the incidence of inner hole sources.

Given that we cannot distinguish on statistical grounds between a short-lived transistion phase experienced by all stars or a long-lived phase undergone by a minority, we need instead to look at the properties of individual sources, bearing in mind that the properties of these sources will partly depend on the observational criteria used for object selection. Therefore, for example, the sample of inner hole sources compiled by Kim et al. (2009) on the basis of a dip in the Spitzer IRS spectrum (in the range m) are largely a group of objects with high accretion rates on to the central star ( yr) and high outer disc masses (). Whereas the X-EUV models are close to being able to reproduce these conditions in a short period briefly after gap opening, they should not be representative of photoevaporative holes in general. In particular, X-EUV models are not compatible with large ( A.U.) holes that are vigorously accreting (such as SZ Cha, GM Auriga or UX TauA; Najita et al. 2007, Calvet et al. 2005, Espaillat et al. 2007). On the other hand, the inner hole sources compiled by Cieza et al. (2008) are qualitatively different, being associated with lower disc masses and low or undetectable levels of accretion: these may indeed include sources with photoevaporating holes in the evolutionary stage c) above. (Likewise the mm survey of Andrews Andrews & Williams 2005 contains objects that may be in this category, although none of these have the high disc masses that are predicted by X-EUV models at the onset of outer disc clearing). It should be stressed that EUV photoevaporation models are mildly biased towards larger hole sizes and that this effect is somewhat stronger in X-EUV models, due to the larger outer disc masses and the gentle decline in wind rate at large radii (Figure 13). We would therefore expect that photoevaporating holes would be more strongly represented amongst sources with large inner holes as one would expect to uncover with Herschel or ALMA.

A final point about the formation mechanism for observed inner hole sources is that it is a robust prediction of photoevaporation models that there should be no flow across the gap once it is opened. Thus following gap opening, the draining inner disc is not re-supplied and the ratio of its mass to the accretion rate onto the star should be a good measure of its lifetime in the draining phase. Observed inner hole sources are generally characterised by the dust content of their inner discs and in most cases the gas mass is unknown. However, measurements of CO emission in TW Hydra and GM Auriga (Salyk et al., 2007) indicate very low inner disc gas masses, which, when combined with the measureed accretion rates in these objects would imply improbably short lifetimes in the absence of re-supply. This argument therefore is a strong indication that the inner disc is re-supplied in these objects, as is a natural expectation in planet clearing models.

10 Conclusions

In this work we have presented the first self-consistent hydrodynamic model of an X-EUV driven photoevaporation wind. The mass-loss rate of 1.4Myr occurs over a large radial extent of 1-70AU. This should be contrasted with the much lower yr values obtained from EUV photoevaporation alone and with the fact that in the EUV case the mass loss rate is more concentrated at a scale of a few AU.

When combined with viscous evolution, X-EUV photoevaporation produces a qualitatively similar pattern of inner hole creation followed by inside out clearing as is found in EUV models. The initial size of the inner hole (about an AU) and the rapid clearing following gap opening is similar in X-EUV and EUV models. The important difference is that this occurs earlier, and at much higher accretion rate, in the X-EUV case. This implies that inner holes can be associated with more massive outer discs and with higher accretion rates on to the star than was possible in the EUV case. Our results thus weaken the objections made previously to EUV photoevaporation with regard to inner hole creation. For example, the X-EUV model can open a gap even in the case of massive outer discs such as that in TW Hydra.

Nevertheless there are generic predictions of photoevaporation models (whether EUV or X-EUV) that appear not to be compatible with features of some of the best studied inner hole sources. Namely: i) photoevaporation models predict that gas and dust in the inner hole, as well as diagnostics of accretion onto the star, should be observable only during inner disc clearing and that during this phase the ratio of inner disc gas mass to accretion rate should be of order the lifetime of this phase. In TW Hydra and GM Auriga, however, the observed ratio is very low and would imply an unacceptably short lifetime of this phase. This low value probably implies that, in contrast to what is predicted in photoevaporation models, the gas and dust in the inner hole is continuously replenished from the outer disc. ii) Photoevaporation models also predict that the phase of inner hole draining is associated with small holes (i.e. of AU scale). There is no mechanism to produce large holes that contain gas and undergo vigorous accretion onto the star (as in GM Auriga).

On the other hand, both EUV and X-EUV models can explain the detection of outer discs in around of stars which show no evidence for inner disc emission/accretion diagnostics (Andrews & Williams, 2005). Future measurements with Herschel or ALMA are likely to uncover a class of objects with large inner holes and we would expect this population to be dominated by photoevaporative holes.

Finally, we stress that the X-EUV case investigated here corresponds to the case of a rather large X-ray luminosity (erg s) and further work is required in order to explore the effect on disc evolution of the wide range of X-ray luminosities measured in T Tauri stars. A clear qualitative prediction is that disc lifetimes are expected to be shorter in X-ray luminous sources and that these sources should be be over-represented among disc-less (Weak Line) T Tauri stars. This prediction is in line with observations: though the larger X-ray luminosity in Weak Line T Tauri stars is conventionally ascribed to suppression of X-rays in accreting systems, we here follow Drake et al. (2009) in suggesting that another factor is that X-rays may instead suppress accretion.


We thank Paola D’Alessio for providing us with the electronic data for the gas density distribution in the disc. We would also like to thank the referee for their comments and suggested imporvements. JEO acknowledges support of a STFC PhD studentship. BE acknowledges funding from a STFC advanced fellowship. RDA acknowledges support from the Netherlands Organisation for Scientific Research (NWO), through VIDI grants 639.042.404 and 639.042.607. This work was partly performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England.


  • Alexander et al. (2004) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2004, MNRAS, 354, 7
  • Alexander (2008b) Alexander, R. D. 2008b, MNRAS, 391, L64
  • Alexander et al. (2005) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2005, MNRAS, 358, 283
  • Alexander et al. (2006a) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216
  • Alexander et al. (2006b) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229
  • Alexander (2008a) Alexander, R. 2008a, New Astronomy Review, 52, 60
  • Andrews & Williams (2005) Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134
  • Armitage et al. (2003) Armitage, P. J., Clarke, C. J., & Palla, F. 2003, MNRAS, 342, 1139
  • Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, New York: Wiley, 1983,
  • Calvet et al. (2002) Calvet, N., D’Alessio, P., Hartmann, L., Wilner, D., Walsh, A., & Sitko, M. 2002, ApJ, 568, 1008
  • Calvet et al. (2005) Calvet, N., et al. 2005, ApJL, 630, L185
  • Cieza et al. (2008) Cieza, L. A., Swift, J. J., Mathews, G. S., & Williams, J. P. 2008, ApJL, 686, L115
  • Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485
  • Albacete Colombo et al. (2007) Albacete Colombo, J. F., Flaccomio, E., Micela, G., Sciortino, S., & Damiani, F. 2007, A&A, 464, 211
  • Ciesla (2007) Ciesla, F. J. 2007, ApJL, 654, L159
  • Cieza et al. (2008) Cieza, L. A., Swift, J. J., Mathews, G. S., & Williams, J. P. 2008, ApJL, 686, L115
  • Chiang & Murray-Clay (2007) Chiang, E., & Murray-Clay, R. 2007, Nature Physics, 3, 604
  • Currie et al. (2009) Currie, T., Lada, C. J., Plavchan, P., Robitaille, T. P., Irwin, J., & Kenyon, S. J. 2009, ApJ, 698, 1
  • D’Alessio et al. (1998) D’Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
  • D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321
  • Drake et al. (2009) Drake, J. J., Ercolano, B., Flaccomio, E., & Micela, G. 2009, ApJL, 699, L35
  • Dullemond & Dominik (2005) Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971
  • Duvert et al. (2000) Duvert, G., Guilloteau, S., Ménard, F., Simon, M., & Dutrey, A. 2000, A&A, 35
  • Eisner et al. (2006) Eisner, J. A., Chiang, E. I., & Hillenbrand, L. A. 2006, ApJL, 637, L133
  • Espaillat et al. (2007) Espaillat, C., Calvet, N., D’Alessio, P., Hernández, J., Qi, C., Hartmann, L., Furlan, E., & Watson, D. M. 2007, ApJL, 670, L135
  • Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X.-W. 2003, MNRAS, 340, 1136
  • Ercolano et al. (2005) Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038
  • Ercolano et al. (2008a) Ercolano, B., Young, P. R., Drake, J. J., & Raymond, J. C. 2008a, ApJS, 175, 5345, 165
  • Ercolano et al. (2008b) Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008b, ApJ, 688, 398
  • Ercolano et al. (2009a) Ercolano, B., Clarke, C. J., & Robitaille, T. P. 2009a, MNRAS, 394, L141
  • Ercolano et al. (2009b) Ercolano, B., Clarke, C. J., & Drake, J. J. 2009b, ApJ, 699, 1639
  • Flaccomio et al. (2003) Flaccomio, E., Damiani, F., Micela, G., Sciortino, S., Harnden, F. R., Jr., Murray, S. S., & Wolk, S. J. 2003, ApJ, 582, 398
  • Flaccomio et al. (2009) Flaccomio, E., Stelzer, B., Sciortino, S., Micela, G., Pillitteri, I., & Testi, L. 2009, arXiv:0906.4700
  • Font et al. (2004) Font, A. S., McCarthy, I. G., Johnstone, D., & Ballantyne, D. R. 2004, ApJ, 607, 890
  • Glassgold et al. (2004) Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972
  • Glassgold et al. (2007) Glassgold, A. E., Najita, J. R., & Igea, J. 2007, ApJ, 656, 515
  • Glassgold et al. (2009) Glassgold, A. E., Meijerink, R., & Najita, J. R. 2009, arXiv:0905.4523
  • Gregory et al. (2007) Gregory, S. G., Wood, K., & Jardine, M. 2007, MNRAS, 379, L35
  • Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539
  • Gregory et L. (2007) Gregory, S., Wood, K., Jardine, M. 2007,MNRAS379,L35
  • Güdel (2008) Güdel, M. 2008, Astronomische Nachrichten, 329, 218
  • Haisch et al. (2001) Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJL, 553, L153
  • Hartmann (1995) Hartmann, L. W. 1995, APSS, 224, 3
  • Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
  • Herczeg et al. (2007) Herczeg, G. J., Najita, J. R., Hillenbrand, L. A., & Pascucci, I. 2007, ApJ, 670, 509
  • Hollenbach et al. (1994) Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654
  • Hughes et al. (2009) Hughes, A. M., et al. 2009, ApJ, 698, 131
  • Igea & Glassgold (1999) Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848
  • Johnstone et al. (1998) Johnstone, D., Hollenbach, D., & Bally, J. 1998, ApJ, 499, 758
  • Kashyap & Drake (2000) Kashyap, V., & Drake, J. J. 2000, Bulletin of the Astronomical Society of India, 28,475
  • Kenyon & Hartmann (1995) Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117
  • Kim et al. (2009) Kim, K. H., et al. 2009, ApJ, 700, 1017
  • Krauss et al. (2007) Krauss, O., Wurm, G., Mousis, O., Petit, J.-M., Horner, J., & Alibert, Y. 2007, A&A, 462, 977
  • Laor & Draine (1993) Laor, A., & Draine, B. T. 1993, ApJ, 402, 441
  • Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 60
  • Maggio et al. (2007) Maggio, A., Flaccomio, E., Favata, F., Micela, G., Sciortino, S., Feigelson, E. D., & Getman, K. V. 2007, ApJ, 660, 1462
  • Meijerink et al. (2008) Meijerink, R., Glassgold, A. E., & Najita, J. R. 2008, ApJ, 676, 518
  • Najita et al. (2007) Najita, J. R., Strom, S. E., & Muzerolle, J. 2007, MNRAS, 378, 369
  • Najita et al. (2009) Najita, J. R., et al. 2009, ApJ, 697, 957
  • Neuhäuser et al. (1995) Neuhäuser, R., Sterzik, M. F., & Schmitt, J. H. M. M. 1995, APSS, 223, 178
  • Natta et al. (2006) Natta, A., Testi, L., & Randich, S. 2006, A&A, 452, 245
  • Parker (1958) Parker, E. N. 1958, ApJ, 128, 664
  • Pascucci (2007) Pascucci, I. 2007, Bulletin of the American Astronomical Society, 38, 872
  • Pascucci et al. (2009) Pascucci, I., & Sterzik, M. 2009, ApJ, 702, 724
  • Preibisch et al. (2005) Preibisch, T., et al. 2005, ApJS, 160, 401
  • Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137
  • Pringle et al. (1986) Pringle, J. E., Verbunt, F., & Wade, R. A. 1986, MNRAS, 221, 169
  • Sicilia-Aguilar et al. (2008) Sicilia-Aguilar, A., Henning, T., Juhász, A., Bouwman, J., Garmire, G., & Garmire, A. 2008, ApJ, 687, 1145
  • Rettig et al. (2004) Rettig, T. W., Haywood, J., Simon, T., Brittain, S. D., & Gibb, E. 2004, ApJL, 616, L163
  • Rice et al. (2003) Rice, W. K. M., Wood, K., Armitage, P. J., Whitney, B. A., & Bjorkman, J. E. 2003, MNRAS, 342, 79
  • Rice et al. (2006) Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619
  • Salyk et al. (2007) Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2007, ApJL, 655, L105
  • Sanz-Forcada et al. (2002) Sanz-Forcada, J., Brickhouse, N. S., & Dupree, A. K. 2002, ApJ, 570, 799
  • Salyk et al. (2007) Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2007, ApJL, 655, L105
  • Schisano et al. (2009) Schisano et al. 2009 In Press
  • Stone & Norman (1992a) Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753
  • Stone & Norman (1992b) Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791
  • Stone et al. (1992) Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819
  • Tarter et al. (1969) Tarter, C. B., Tucker, W. H., & Salpeter, E. E. 1969, ApJ, 156, 943
  • Quillen et al. (2004) Quillen, A. C., Blackman, E. G., Frank, A., & Varnière, P. 2004, ApJL, 612, L137
  • Varnière et al. (2006) Varnière, P., Blackman, E. G., Frank, A., & Quillen, A. C. 2006, ApJ, 640, 1110
  • Whitworth & Clarke (1997) Whitworth, A. P., & Clarke, C. J. 1997, MNRAS, 291, 578
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description