## Photoevaporation of protoplanetary discs with PLUTO+PRIZMO

## I. Lower X-ray-driven mass-loss rates due to enhanced cooling

A. D. Sellek 1 2 , , T. Grassi 3 , G. Picogna 4 , Ch. Rab 3 4 , , C. J. Clarke 2 , and B. Ercolano 4 5 ,

- 1 Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands
- e-mail: sellek@strw.leidenuniv.nl
- 2 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
- 3 Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, D-85748 Garching, Germany
- 4 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 Munich, Germany
- 5 Exzellenzcluster 'Origins,' Boltzmannstr. 2, D-85748 Garching, Germany

Received March 28, 2024; accepted July 31, 2024

## ABSTRACT

Context. Photoevaporation is an important process for protoplanetary disc dispersal but there has so far been a lack of consensus from simulations over the mass-loss rates and the most important part of the high-energy spectrum for driving the wind.

Aims. We aim to isolate the origins of these discrepancies through carefully-benchmarked hydrodynamic simulations of X-ray photoevaporation with time-dependent thermochemistry calculated on the fly.

Methods. We conduct hydrodynamic simulations with pluto where the thermochemistry is calculated using prizmo . We explore the contribution of certain key microphysical processes and the impact of using di ff erent spectra used previously in literature studies.

Results. We find that additional cooling results from the excitation of O by neutral H, which leads to dramatically reduced mass-loss across the disc compared to previous X-ray photoevaporation models, with an integrated rate ∼ 10 -9 M yr ⊙ -1 . Such rates would allow for longer-lived discs than previously expected from population synthesis. An alternative spectrum with less soft X-ray produces mass-loss rates around a factor of 2-3 times lower. The chemistry is significantly out of equilibrium, with the survival of H2 into the wind aided by advection. This leads to its role as the dominant coolant at 10s au - thus stabilising a larger radial temperature gradient across the wind - as well as providing a possible wind tracer.

Key words. Astrochemistry - Hydrodynamics - Methods: numerical - Protoplanetary disks - Stars: winds, outflows - X-rays: stars

## 1. Introduction

Protoplanetary discs are the site of planet formation, as evidenced by direct observations of protoplanets around the T Tauri star PDS70 (Keppler et al. 2018; Müller et al. 2018) as well as candidates around the more massive Herbig stars AB Aurigae b (Currie et al. 2022) and HD 169142 (Gratton et al. 2019; Hammond et al. 2023). The formation, growth, and migration of planets in protoplanetary discs is ultimately limited by disc dispersal, thought to happen rapidly (Skrutskie et al. 1990; Simon &amp; Prato 1995; Wolk &amp; Walter 1996) after around 2-8 Myr of evolution (Pfalzner et al. 2022) depending on the star-forming environment (Michel et al. 2021) and stellar mass (Ribas et al. 2015).

cretion become less significant as the disc accretion rate drops. Therefore, these winds become more tenuous and hence more optically thin to high-energy radiation, allowing it to more easily reach the disc and drive photoevaporation (Kunitomo et al. 2020; Weder et al. 2023). On the other hand, photoevaporation rates are largely insensitive to the properties of the underlying disc (at least until any transparent gas cavity opens): the disc adjusts to deliver the mass flux expected from the thermochemistry of the wind (Owen et al. 2012) and the mass-loss rate remains high until the disc becomes optically thin to the dominant radiation field driving photoevaporation.

Disc winds are widely evidenced throughout the disc lifetime (Pascucci et al. 2023) and are one important factor contributing to the evolution and dispersal of these discs. Typically two paradigms are considered: photoevaporation, where thermal pressure gradients resulting from heating by high-energy stellar radiation drive the wind, and magnetically-driven winds, where gradients in the magnetic pressure and centrifugal forces (resulting from material being forced to co-rotate with the magnetic field) launch gas outwards (although Bai 2017, define the intermediate case of magnetothermal winds as those where both magnetic and thermal e ff ects play a role).

Photoevaporation is a particularly attractive paradigm for the late-stage dispersal of discs. The mass-loss rates that can be produced by cold magnetically-driven wind models powered by ac-

Photoevaporation also explains a number of observational constraints which suggest dispersal progresses from the inside out. The 'UV switch' model of Clarke et al. (2001) shows how photoevaporation should carve a gap, and then a cavity, inside the inner few au. Such a mode of dispersal is supported by the relative disc lifetimes at di ff erent wavelengths (Ribas et al. 2014) and the evolutionary loci of discs in colour-colour space (Koepferl et al. 2013), which imply the loss of the warmer inner material - which emits at shorter wavelengths - first.

However, theoretical predictions for the strength of photoevaporative winds lack consensus. Mass-loss rates in the literature have varied from 10 -10 -10 -9 M yr ⊙ -1 for models driven by Extreme Ultraviolet (EUV) (Hollenbach et al. 1994) to 10 -8 -10 -7 M yr ⊙ -1 for X-ray-driven models (Owen et al. 2010, 2011a, 2012; Picogna et al. 2019; Ercolano et al. 2021; Picogna et al.

2021) or models driven by Far Ultraviolet (FUV) (Gorti et al. 2015; Nakatani et al. 2018a,b).

In part, these di ff erences may be attributed to di ff erent methodologies - since the problem of thermochemistry on the fly with hydrodynamics is very expensive, di ff erent works have compromised on di ff erent factors. For example, Hollenbach et al. (1994) or Gorti &amp; Hollenbach (2009); Gorti et al. (2015) do not explicitly solve the hydrodynamics, but estimate what solution should result given their calculated temperature structure. On the other hand, X-ray photoevaporation models since Owen et al. (2010) have typically solved the hydrodynamics with temperatures prescribed as a function of the ionisation parameter ξ = LX nr 2 (Tarter et al. 1969), which roughly speaking represents the ratio of the X-ray photoionisation heating to the two-body collisional cooling. Such an approach allows a large parameter space to be explored at the expense of not being able to consider e ff ects such as out-of-equilibrium thermochemistry or adiabatic cooling resulting from the PdV work done by the expanding gas. The ξ -T relationships have also been calculated assuming purely atomic gas using the mocassin code (Ercolano et al. 2003, 2005, 2008) and only EUV X-ray radiation. In recent years, there have / been moves towards on-the-fly thermochemistry and hydrodynamics both for photoevaporation (e.g. Wang &amp; Goodman 2017; Nakatani et al. 2018a,b; Komaki et al. 2021) and also for magnetohydodynamic (MHD) winds (Wang et al. 2019; Gressel et al. 2020), but such e ff orts naturally come at the expense of solving the thermochemistry in detail and hence are very sensitive to the simplifying assumptions made.

The fact that, typically, several things - such as the irradiating spectrum, possible coolants, radiative transfer method, grid extent, explicit calculation of hydrodynamics, and assumption of thermochemical equilibrium - di ff er between each methodology, has made it very di ffi cult to make direct comparisons between the above works and clearly assess the dominant factors responsible for the divergent results. For example, Wang &amp; Goodman (2017) - who favoured a low-mass-loss-rate EUV-driven wind attributed the higher mass-loss rates of X-ray models to the lack of molecular cooling. On the other hand, Sellek et al. (2022) demonstrated that a lack of e ff ective optical line cooling in the work by Wang &amp; Goodman (2017) meant that their predicted wind temperatures were too high so despite the additional molecular cooling, this thermochemical treatment was still not su ffi -ciently complete. More importantly, Sellek et al. (2022) showed that the results were very sensitive to the shape of the X-ray spectrum, finding that the simplified spectrum used by Wang &amp; Goodman (2017) where they assumed all the X-rays to have an energy of 1 keV meant that the X-rays interacted too weakly with the gas to achieve significant heating; that is to say Wang &amp; Goodman (2017) were prevented from finding an X-ray wind primarily by their lack of soft X-rays.

In this work, we seek to build the most comprehensive photoevaporation model to date, using prizmo (Grassi et al. 2020) to combine hydrodynamics with thermochemistry calculated on the fly; for this first application, we focus on X-ray-driven winds. We stress the importance of making sure that our approach is su ffi ciently complete to capture thermochemistry in all possible wind physical environments, including ionised, neutral atomic, and molecular gas. As such, we emphasise benchmarking the temperatures calculated by prizmo against previous models and other codes. Important aspects of this include explicitly exploring the role of what we find to be the most critical thermochemical processes, and enabling direct comparison to other models at each stage with appropriate spectra.

Section 2 sets out the key features of prizmo and any choices we make for the thermochemistry, while Sect. 3 explains how we use prizmo in conjunction with pluto for the hydrodynamics. Section 4 presents benchmarks of the temperature structure of winds in di ff erent regimes to illustrate the key successes of prizmo and elucidate important di ff erences to previous works. We then demonstrate that hydrodynamic simulations using prizmo substantially alter the temperature profile and reduce the mass-loss rates of X-ray-driven winds in Sect. 5. We also derive the corresponding mass-loss profile, establish the importance of non-equilibrium thermochemistry, and explore how the mass-loss rate is a ff ected by the microphysics and irradiating spectrum. In Sect. 6, we briefly discuss the implications of these results for simulations of photoevaporation, disc evolution and observations of disc winds. We summarise our conclusions in Sect. 7.

## 2. Methods: thermochemistry

Previous works have highlighted the need for an approach to photoevaporation modelling which a) uses a broad, detailed, spectrum for the radiation, b) includes su ffi cient atomic transitions to provide cooling regardless of the temperature (and ionisation), c) includes cooling of molecules with abundances calculated self-consistently with the local radiation field, d) is done on the fly with hydrodynamics to avoid assuming radiative thermal equilibrium in the regions where adiabatic cooling may be prevalent.

Such a comprehensive approach demands a fast code that can solve the thermochemistry - a term here used to describe collectively the heating and cooling that results from the chemical composition of the gas and interactions between its constituents - across molecular and atomic regimes and can be called as a library by hydrodynamics codes using the operator splitting technique (e.g. Teyssier 2015). prizmo (Grassi et al. 2020) was developed specifically for this purpose. prizmo has been significantly updated since its first release, and the main changes are described in Sect. 2.1 1 .

prizmo is configured by means of a python preprocessor to prepare the appropriate thermochemistry fortran modules and corresponding required data files. These modules include C wrappers, allowing them to be compiled with codes such as pluto (Sect. 3).

The key inputs to the preprocessor are:

- -Atomic data (level energies, degeneracies, Einstein coe ffi -cients, and fits to collisional de-excitation rates)
- -A list of (photo)chemical reactions with rate coe ffi cient expressions
- -A user-defined spectrum and number of energy bins
- -Frequency-dependent dust optical constants

Our preprocessor inputs are detailed in Sections 2.2-2.5.

A key feature of prizmo is its multi-frequency approach: all photochemical cross-sections and dust opacities (which are integrated over a grain size distribution) resulting from the selected reaction network and dust optical constants are tabulated as a function of photon energy. All the radiation-related processes (e.g., photoheating) therefore use the same multi-frequency approach where possible. This significantly improves the consistency of the results, by allowing prizmo to calculate photochemistry using the local radiation field (rather than assuming rates based on a standard field).

1 This version of prizmo is available from GitHub at: https:// github.com/tgrassi/prizmo , commit b2158ae .

## 2.1. Updates to prizmo

## 2.1.1. Structure/architecture

Since Grassi et al. (2020), prizmo has been modified to reduce the global computational time and improve the code readability, and further modifications were made to the physics and chemistry involved in the present application. Nevertheless, despite these improvements a ff ording an approximately 50% speed-up, the thermochemistry still dominates the execution time and is expected to typically be around 100 times longer than would be required to have no additional impact compared to pure hydrodynamics. This section describes the improvements from Grassi et al. (2020).

## 2.1.2. Rate coefficients calculation

We increased the vectorisation of rate coe ffi cient calculations in the chemical network by grouping chemical rate coe ffi cients by the same arithmetic expression, for example α ( T / 300 K) β exp ( -γ/ T ) and α ( T / 300 K) β belong to di ff erent groups. Since most of the rates have relatively simple algebraic forms, with this approach we obtain a speed-up similar to interpolating the rate coe ffi cients in temperature in the log-log space.

We have also approximated some rates in order to precompute common algebraic quantities; for example, ( T / 300 K) -0 49 . has been approximated to ( T / 300 K) -0 5 . , which is then a common term to several reaction rates, and hence can be precomputed. These small rate variations do not impact the final results, and have a ∼ 1 -2% reduction on the computational time. However, it is important to notice that small speedups have a relatively large impact on simulations that require thousands of core hours.

## 2.1.3. Atomic cooling solver

The linear solver in the atomic cooling routine has been extended to solve linear systems analytically up to 5 × 5 matrices without using L apack . This considerably improved the speed of the atomic cooling linear solver for the larger atomic line systems needed to accurately capture atomic cooling across the range of temperatures found in protoplanetary discs and winds (see Sect. 2.3, also Sellek et al. 2022). However, the atomic cooling routine occupies a large fraction of the computational time, mainly due to fitting functions in the log-log space employed to evaluate the rate coe ffi cient with collisional partners. Fits have a relatively large computational impact because they require 10 x -type operations, where x is a double precision number. We reduce the impact of their original algebraic functions by vectorising fitting operations over temperature. In addition to this, the combination of rate interpolations and solving the N × N linear system analytically provides a good computational speed with respect to, for example, multidimensional linear fits that require a dimension for each collisional partner and for temperature and, analogously to the other method, employ 10 -type operations. x

## 2.1.4. Photoionisation

Photoionisation is calculated using the cross-sections of Verner et al. (1996) for each of the neutral atomic species in our network (H, He, C, O). Modelling each species separately - rather than using a single approximated total photoionisation cross-section -is necessary to follow the ionisation balance of each species individually and self-consistently with the heating. It also al- lows flexiblity in treating variations in the gas compositions, for example due to altered elemental abundances or regions of enhanced ionisation. However, note that as well as due to the direct radiation field of the star, ionisation can occur as a result of the di ff use radiation field created by recombinations. Since we only conduct radial ray tracing in this work, we cannot follow these photons accurately and hence assume they are reabsorbed on the spot. Although the di ff use field of one element can potentially ionise those with a lower ionisation potential (e.g., ground-state recombinations of He may ionise H), for simplicity we assume that the reabsorption occurs by the same element, noting that this may somewhat underestimate the energy liberated into the gas by photoelectrons but this will only make any significant di ff erence in regions dominated by EUV heating while it is a much smaller correction if the original ionisation was by an X-ray photon.

Moreover, since doubly-ionised and higher species were not included in the network we use in this work (see Sect. 2.2), we cannot include here photoionisation reactions of the singly ionised species. Similarly, neither can the code account for multiple ionisation by the Auger-Meitner e ff ect, although the crosssection due to ionisation of inner (K) shell electrons is included (Verner &amp; Yakovlev 1995) and is assumed to always result in a single electron being ejected. Consequently, the entire yield over the (outer shell) ionisation energy was available to contribute to heating (assuming that any excited products from inner shell ionisation de-excite rapidly through collisions). We do not produce secondary electrons, either by the Auger e ff ect or by rapid collisional ionisation, save for what may already be produced by the thermalised electrons in the network. This may overestimate the heating by a small factor (Maloney et al. 1996, e.g.). The e ffi ciency of the heating is thus 100% for photoionisation of completely neutral gas and 0% for photoionisation of fully (singly-)ionised gas. We intend to address these limitations in future works.

Neither X-ray photoionisation of molecules nor of dust is included in the present code (while the cross-section of molecules can be well-approximated by the sum of that of the atoms, the branching ratios of the outcomes are poorly known).

## 2.1.5. Molecular cooling and heating

It has been suggested that a key process missing in previous Xray photoevaporation calculations that relied on mocassin was molecular cooling (Wang &amp; Goodman 2017). Hence, we include it here according to tables computed from Omukai et al. (2010) for CO cooling and from the piece-wise functions of Glover &amp; Abel (2008); Glover (2015) for H2 cooling. While we do not include cooling from H2O, we do not see it at su ffi cient levels in the wind for it to be an important coolant.

Molecules can also contribute to heating: photodissociation and pumping of H2 are included, each proportionally to the photodissociation rate calculated by the chemical network.

## 2.1.6. Photoelectric heating

We have changed the photoelectric routine following Kamp &amp; Bertoldi (2000) as reported in Woitke et al. (2009), their Eq. (93), that is related to the strength of the FUV radiation field χ . This relation is less accurate and tunable than the full formalism of Weingartner &amp; Draine (2001) and Weingartner et al. (2006), but we found it to be less prone to show unphysical heating when coupled with the radiative transfer, mainly due to implementation details. Note that this formulation gives only the heating

from silicate grains; in the current code no heating from Polyaromatic Hydrocarbons (PAHs) is included. PAH heating can be of consequence in the upper layers of the disc at 10s au and in the outer disc (Woitke et al. 2009) and while this shouldn't a ff ect the models presented in this paper (which include little FUV in the spectra), it will need improving in future work where we aim to apply the code to FUV photoevaporation.

## 2.1.7. Dust temperature and dust cooling

The previous code version computed the dust temperature using a bisection method applied to Eq. (41) in Grassi et al. (2020). Weobtain faster integration times by linearly interpolating a pregenerated table of the form T d = f ( Γ a , T n , ), where Γ a is the absorption integral using the impinging radiation, as described in Eq. (43) of Grassi et al. (2020), T is the gas temperature, and n is the total gas density. Analogously, we compute the dust cooling and heating with the same approach using two pre-computed separated tables. Using two tables allows dividing heating (positive) and cooling (negative) in the log-log space.

## 2.2. Species &amp; reactions in our network

Our network is based around the elements H, He, C, and O with interstellar medium (ISM) gas-phase abundances consistent with previous photoevaporation models (e.g. Ercolano et al. 2009; Wang &amp; Goodman 2017): He H / = 0 1, C . / H = 1 4 . × 10 -4 and O H / = 3 2 . × 10 -4 (Savage &amp; Sembach 1996).

Weuse the chemical network described in Grassi et al. (2020, Appendix G), which is designed to obtain the correct temperature structure, rather than the correct chemical structure of the less abundant chemical species. It is loosely based on Photodissociation Region (PDR) chemistry and includes the neutral and singly-ionised versions of each atom, electrons, 22 neutral ionised gaseous molecules, and finally CO and H2O ices, all / for a total of 33 species. A total of 290 reactions are included, covering two-body gas-phase reactions, photoionisation, photodissociation, cosmic-ray-induced reactions, formation of H2 on dust grains, and freeze out thermal desorption of CO and H2O / onto from dust grains. /

Note that we exclude S from our model as its abundance and main carriers in protoplanetary discs are somewhat debated. Although previous models used the ISM measurements of Savage &amp; Sembach (1996) to determine gas-phase abundances - in which S was not measured to be depleted with respect to solar S does appear to become progressively depleted as material becomes increasingly dense during the molecular cloud and core phases of star formation (Hily-Blant et al. 2022; Fuente et al. 2023). Similarly, protoplanetary discs exhibit gas-phase S depletion (Kama et al. 2019) of around 90 percent, implying a highly abundant refractory reservoir; meteoritic studies show abundant FeS minerals in chondrites (Kallemeyn et al. 1989), which may therefore be a likely candidate. As a consequence, previous models of wind emission overestimate how bright and frequentlydetected S emission lines should be (Pascucci et al. 2020) and how much S can contribute to cooling.

## 2.3. Atomic cooling

Atomic line cooling happens via collisional excitation of an atom or ion to an electronically excited state, followed by radiative de-excitation. To find the total cooling emission, we therefore need to compute the population of each atomic level. To do so, we assume equilibrium between collisional and radiative transitions and solve the linear system formed by the collisional (de)excitation rates and the spontaneous emission rates. In particular, the relative abundance of the i th level xi is determined by

$$\ e \, \text{we aim} \quad - x _ { i } \sum _ { \ell } \left ( n _ { \ell } \sum _ { i \neq j } k _ { i j } ^ { \ell } \right ) - x _ { i } \sum _ { j < i } A _ { i j } + \sum _ { \ell } \left ( n _ { \ell } \sum _ { i \neq j } k _ { i j } ^ { \ell } x _ { j } \right ) + \sum _ { j > i } A _ { j i } x _ { j } = 0 \,, \\ \quad \,. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

where n ℓ is the abundance of the ℓ th collider ( ℓ ∈ { e -, H H , + , ortho -H2 , para -H2 ), } k ℓ i j is the temperaturedependent collisional rate coe ffi cient that takes into account the transition between the i th and j th levels when (de-)excited by the ℓ th collisional partner, and Aij is the Einstein coe ffi cient of the i → j radiative de-excitation.

By solving the linear system using Eq. (1) for each level and replacing the last level equation with the mass conservation (i.e. P i xi = 1), it is possible to compute the level population xi . This system of equations can be solved with a matrix approach as discussed in 2.1.3.

After xi are computed, the cooling emission of each transition becomes Λ i j = n Aij xi ∆ Eij , where n is the total abundance of the specific atomic coolant, and ∆ Eij is the energy di ff erence between the i th and j th levels. All of the transitions included here (with the exception of Lyman α , discussed below), are forbidden transitions, to which the wind is generally optically thin and thus we assume an escape probability of unity. To confirm this, we checked the line-centre cross-section of the lines assuming a thermal broadening vrmth ∼ 1 km s -1 according to σ = λ 3 A 21 8 √ 2 π / 3 2 v th (Osterbrock &amp; Ferland 2006, Equation 4.45). For the shorter wavelength optical lines this is typically ≲ 10 -21 cm 2 , while those of the FIR lines are somewhat higher at ≲ 10 -17 cm 2 (the highest being O i 145 µ m with 8 6 . × 10 -18 cm 2 . However, given the atomic abundances of ∼ 10 -4 , the cross-sections per hydrogen atom are ≲ 10 -21 cm 2 . Therefore, given that soft Xrays typically penetrate to ∼ 10 21 cm -2 , we expect the lines to be at most marginally optically thick in the heated wind gas. The total cooling is the sum of all the allowed (or relevant) atomic transitions.

In Table 1, we summarise the number of levels used for each of our atomic coolants (O, O + , C, C + ) and the provenance of the data for the collisional (de-)excitation strengths rates. As de-/ tailed above, we include e -, H H , + , ortho -H2 , and para -H2 as colliders wherever data is available. This ensures that we capture collisional excitement and cooling as far as possible in ionised, neutral atomic and molecular environments. In most cases, suitable numerical fits to the appropriate rates calculated at di ff erent temperatures were available in Draine (2011, Appendix F). For excitation of O and C by H + we used piecewise fits from Glover &amp; Jappsen (2007). For collisional (de)excitation of O by H, ortho-H2 and para-H2, we use the latest calculations by Lique et al. (2018) 2 to which we perform our own fits, presented in Appendix A. Level energies, degeneracies, and Einstein coe ffi cients of spontaneous emission Aij , were taken from the National Institute of Standards and Technology database (Kramida et al. 2022).

At the temperature range of interest ≲ 10 4 K, we expect emission from lines with excitation temperatures up to a few

2 Note that the values quoted in the paper have been corrected (Wolfire et al. 2022), and are available in the Leiden Atomic and Molecular Database (LAMDA, van der Tak et al. 2020) at https://home.strw. leidenuniv.nl/~moldata/datafiles/oatom@lique.dat .

Table 1: References for the atomic collisional rate data used by our version of prizmo .

| Species   | n levels                                                             | Source of Collisional Rate Data by Collider                                    | Source of Collisional Rate Data by Collider   | Source of Collisional Rate Data by Collider   |
|-----------|----------------------------------------------------------------------|--------------------------------------------------------------------------------|-----------------------------------------------|-----------------------------------------------|
|           | e -                                                                  | H                                                                              | H +                                           | ortho-H 2 &para-H 2                           |
| O         | 5 3 P J j - 3 P J i : Bell et al. (1998); otherwise Pequignot (1990) | 3 P J j - 3 P J i : Lique et al. (2018); 1 D 2 - 3 P J i : Krems et al. (2006) | 3 P J j - 3 P J i : Glover &Jappsen (2007)    | 3 P J j - 3 P J i : Lique et al. (2018)       |
| O +       | 5 Tayal (2007)                                                       | ...                                                                            | ...                                           | ...                                           |
| C         | 3 Glover &Jappsen (2007)                                             | Abrahamsson et al. (2007)                                                      | Glover &Jappsen (2007)                        | Schroder et al. (1991)                        |
| C +       | 2 Tayal (2008)                                                       | Barinovs et al. (2005)                                                         | ...                                           | Flower &Launay (1977)                         |

times ≲ 10 4 K, that is, with wavelengths into the optical. If these optical coolants are missing, wind temperatures may be overestimated (Sellek et al. 2022). Consequently, we include all of the fine-structure terms in the ground configuration of O and O + (resulting in 5-level systems) such that in both neutral and ionised environments, we can have optical line cooling. However, we keep C and C + to a 3- and 2-level system, respectively, due to the limited availability of collision strengths and the increased computational overhead of calculating the level population for higher-level systems; O has typically been found to be a more important coolant than C.

By comparing to a list of the most important coolants in the wind region in the fiducial mocassin model of Sellek et al. (2022), we confirmed that modulo the above caveats, our approach captures most of the lines that contribute to cooling at more than the percentage level, and that beyond this the model is limited by the absence of additional coolants (e.g. S, Si, Fe) rather than the number of levels.

Lyman α cooling is included according to (Spitzer 1978) and assumed to be e ff ectively thin as at the densities of interest, any reabsorbed photons are likely to be quickly re-emitted or reprocessed to optically thin wavelengths by dust (Sellek et al. 2022).

## 2.4. Spectra

The choice of the irradiating spectrum can be highly important for the outcome of photoevaporation simulations (Sellek et al. 2022). To enable comparison with the most recent generation of X-ray photoevaporation models, our fiducial spectrum is the 'Spec30' X-ray spectrum of Ercolano et al. (2021). We normalise this spectrum to a typical LX = 2 04 . × 10 30 erg s -1 (as measured between 0 5 . -5 keV) for a solar-mass star as per Picogna et al. (2021). To this, we also add a 5000 K blackbody assuming a stellar radius of 1 M ⊙ ; this is to ensure reasonable dust temperatures through heating by optical UV photons. The / spectrum is shown in Fig. 1.

Our fiducial spectrum is used throughout our benchmarks, our main simulation and our explorations of cooling and microphysics. In addition, we run a single simulation (see Sect. 5.2.1) with the X-ray spectrum used by Nakatani et al. (2018b); Komaki et al. (2021) - which is based on a model of the observed X-ray spectrum of TW Hya (Nomura et al. 2007) - which we normalise and add to the blackbody in the same way as the fiducial spectrum.

## 2.5. Dust

The dust number density n d was assumed to follow a grain size distribution of φ ( a ) = d n d( a ) d a ∝ a -3 5 . power law between 5 × 10 -7 and 2 5 . × 10 -5 cm (Mathis et al. 1977), with bulk density of 3 g cm -3 (e.g., see Zhukovska et al. 2008; Grassi et al. 2017). The optical properties are calculated with the Mie scattering theory

Fig. 1: Comparison of the two spectra (E21: Ercolano et al. 2021, blue); (N18: Nakatani et al. 2018b, green) investigated in this work with di ff erent regions of the spectrum highlighted and the respective luminosities given in the legend.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-6ed4be1a-ee9c-4e8a-a7aa-7eb95f4da192.jpg)

using the dielectric functions of Draine (2003a,b) which assume a mix of a graphitic carbonaceous component and amorphous astronomical silicates.

These properties are all suggestive of ISM-like grains which have not undergone grain growth. As such, all these grains would be expected to be entrained by a wind (e.g. Hutchison et al. 2016; Franz et al. 2020; Booth &amp; Clarke 2021) so we assume a uniform ISM dust-to-gas mass ratio of 10 -2 everywhere. However, note that grain growth is expected to occur within 10 5 yr of disc evolution and thereafter, maybe as little as ≲ 10% of the dust mass can enter the wind (Franz et al. 2022a). Our results, therefore, likely overestimate the abundance of dust in the wind and hence provide an upper limit on its e ff ects on the thermal structure. While we test our sensitivity to this in Sect. 4.2 and discuss what a more reasonable value might be in light of our results in Sect. 6.2, a more consistent approach with variable dust entrainment is beyond the scope of this exploratory study.

## 3. Methods: hydrodynamics

## 3.1. Overall workflow

Our hydrodynamic simulations are conducted in a three-stage process using pluto (Mignone et al. 2007).

We start with an initial disc model (originally presented in Picogna et al. 2021) generated using the D'Alessio Irradiated Accretion Disc (DIAD) models (D'Alessio et al. 1998, 1999, 2001, 2006) for a gas disc mass of 0 045 M . ⊙ (extending to 400 au) orbiting a K6 star (with M ⋆ = 1 M , ⊙ L ⋆ = 2 335 L . ⊙ , R ⋆ = 2 615 R . ⊙ , T ⋆ = 4278 K). The dust is assumed to be well-mixed, with a maximum grain size in the atmosphere of 2 5 . × 10 -5 cm and a minimum grain-size of 5 × 10 -7 cm(the same as we assume in prizmo ). We use this to calculate initial conditions for the wind by evolving the hydrodynamics for 50 orbits at 10 au, with the wind temperatures calculating temperatures using the 'Spec30' ξ -T relationship of Ercolano et al. (2021). Aside from details of the extent resolution of the numerical grid, this is e / ff ectively identical to the 1 M ⊙ case from Picogna et al. (2021) and so we refer to this as our Picogna et al. (2021) analogue simulation.

We use the Picogna et al. (2021) analogue simulation to provide a density input to pluto+prizmo that already has a wind such that the initial conditions are likely to be closer to the solution; this results in a smoother and more rapid approach towards the steady state without producing transient features due to the sudden expansion of the heated disc material into the wind region.

Secondly, we then run pluto+prizmo with all the hydrodynamic processes disabled in order to perform a 'burn-in' in which only the chemical network and temperature are evolved for 0.1 Myr. This allows us to reach an equilibrium where these quantities are consistent with the density profile. This also allows us to benchmark the temperatures calculated by prizmo against the temperatures prescribed from the ξ -T relationships and identify areas of major di ff erence (see Sect. 4.2).

Finally, the hydrodynamic processes are re-enabled. The simulation is first run for 50 orbits at 10 au (approximately one gravitational radius for 10 4 K gas). Since the structure and composition of the wind settle from the inside out, we are then able to freeze the inner half of the grid ( r &lt; 17 3 au) in order to speed . up the calculation during a further 150 orbits.

## 3.2. Hydrodynamics-thermochemistry interface

Weperform several necessary modifications to PLUTO to couple it with prizmo 3 .

The hydrodynamics and thermochemistry are solved alternately within pluto 's Strang splitting scheme. On evennumbered steps, hydrodynamics, including advection of species, first takes place over the hydrodynamic time step δ t ; assuming no errors are encountered (see below), prizmo is then called as part of the 'SourceStep'. On odd-numbered steps, the order of the two processes is reversed with thermochemistry taking place first followed by the hydrodyanmics.

The mass fraction of each species in the network is tracked during the hydrodynamics using pluto 's tracer variables. By default, pluto does not know that these should - like the density - be non-negative and therefore in principle prizmo can receive negative abundances leading to a fatal error. To avoid this, we added checks to pluto 's 'ConsToPrim' function where similar

3 This version of the code is available from GitHub at: https: //github.com/AndrewSellek/PLUTO\_PRIZMO/tree/chemistry , commit 8e525b9

checks are carried out on the energy, pressure and density. If the tracer mass-fraction ϵ i &gt; -1 × 10 -6 and the tracer density ρ i = ϵ ρ &gt; i -ρ floor the mass-fraction is simply set to 0, otherwise the cell is flagged for recalculation, first using di ff erent solvers ( pluto 's default FAILSAFE method), or, if that fails, by reducing the timestep by a factor 10 until significant negative values are no longer produced. If no solution is found after ten retries, we consider that the solution does not converge so the code returns an error and aborts; once in our fiducial run we required four retries otherwise no more than two were ever needed.

Within the thermochemical step, the first task is to update the radiation field. For this we conduct radial tracing, considering both geometric dilution and the attenuation:

$$\L i h e \, \text{wind} \\ \text{specific} \, \text{$j_{\nu}^{i+1}=j_{\nu}^{i+1}e^{-\tau_{\nu}}\frac{r_{i}^{2}}{r_{i+1}^{2}}$.}$$

The attenuation e -τν calculated by prizmo includes absorption of radiation by the dust as well as all relevant gas-phase species. Note that by conducting radial ray tracing we neglect the scattering of radiation which can be important for UV photons but is negligible for the X-ray photons that ultimately drive the wind we see. Subsequently, prizmo then advances the chemical network in each cell of the simulation for a total time δ t . In so doing, the gas temperature is evolved along with the number densities of each species; the updated values are passed back to pluto in the form of the gas pressure and the tracer variables.

However, within the underlying disc, the dust radiative transfer is an inherently di ff usive problem. Since we only conduct 1D radiative transfer with no scattering of radiation - rather than adopting an expensive Monte Carlo (Pinte et al. 2006; Dullemond et al. 2012, e.g.) or Flux-Limited Di ff usion (Levermore &amp; Pomraning 1981; Kuiper et al. 2010) approach prizmo alone significantly underpredicts dust temperatures below the optical photosphere. This is a problem since, due to the tight dust-gas coupling expected here, the gas temperature would be similarly low and hence the disc would lose a lot of vertical pressure support and shrink. We thus follow the typical approach, which (Owen et al. 2010; Wang &amp; Goodman 2017) is to pin the temperature structure to a pre-calculated dust temperature structure (e.g., Chiang &amp; Goldreich 1997; D'Alessio et al. 2001) beyond a certain column density. Therefore at high column densities we smoothly interpolate between the prizmo temperatures T PRIZMO and the initial conditions (which given the equivalent pinning done in the ξ -T step will be the D'Alessio et al. 1998, 1999, 2001, 2006, temperature, T DIAD) according to the following equation (where N 23 = N tot / 10 23 cm -2 ):

$$\overset { \text{$n$ even-} } { \text{$i$ species} }, \quad \underset { \text{$issuming} } { T } _ { \text{gas} } = \begin{cases} T _ { \text{PRIZMO} } & N _ { \text{tot} } / \text{cm} ^ { - 2 } \leq 1 0 ^ { 2 1 } \\ \frac { T _ { \text{PRIZMO} + N _ { 2 3 } ^ { 2 } T _ { \text{DIAL} } } } { 1 + N _ { 2 3 } ^ { 2 } } & 1 0 ^ { 2 1 } \leq N _ { \text{tot} } / \text{cm} ^ { - 2 } \leq 1 0 ^ { 2 5 } \\ T _ { \text{DIAL} } & 1 0 ^ { 2 5 } \leq N _ { \text{tot} } / \text{cm} ^ { - 2 } \,. \end{cases}$$

We show and discuss an example of how this a ff ects the dust gas / temperature in Sect. 5.

This method also significantly benefits the computational cost by speeding up the calculation in high-density cells, since the temperature and chemistry are no longer coupled.

## 3.3. Hydrodynamics setup

We conduct the full hydrodynamics + thermochemistry simulations on a 2D spherical grid with 400 cells in the radial direction extending from 0 75 . -400 au and 70 in the latitudinal direction

between the midplane and polar axis (spaced to maximise the resolution in the wind-launching region). For the burn-in with thermochemistry only (and Picogna et al. (2021) analogue), we extend the grid inwards to r = 0 25 au using an extra 70 cells. . This is to ensure a steady solution at the inner radial boundary, as testing showed that if the inner cells received completely unattenuated radiation, the results were quite unstable as the disc would develop a very thin, hot region in the inner row of cells that should not be present (assuming a full disc has no cavity). The region 0 25 . -0 75 au therefore acts to attenuate the radiation . reaching the disc and the subsequent hydrodynamic solution is obtained using the radiation field measured in the burn-in runs at 0 75 au. .

In the latitudinal direction we used pluto 's 'polaraxis' and 'eqtsymmetric' boundary conditions at the rotation axis and disc midplane. At the outer radial boundary, we used outflowing boundary conditions, while at the inner radial boundary we adopted user-defined conditions to ensure stability; these are identical to outflowing conditions except that in order to prevent material entering the grid from radii where we do not expect any wind, the radial component of velocity is set in the manner of reflective boundary conditions if it has a positive sign.

Finally, we used the following pluto settings: linear reconstruction, RK2 timestepping, the MINMOD\_LIM limiter and the HLLC solver. Moreover, since early tests with our code were susceptible to negative pressures (a common problem which can result from the fact that the usual equations of hydrodynamics aim to conserve total and kinetic energy and obtain the thermal pressure as the di ff erence, which can become negative due to discretisation errors) we used pluto 's option to instead conserve entropy, which conserves thermal energy pressure more directly. / We found that using this option adequately avoided these errors.

## 4. Static benchmarks of prizmo

To validate the working of prizmo we first perform two benchmarks in the limits corresponding to the extremes of literature photoevaporation models: a low-density high ionisation param-/ eter EUV model from (Wang &amp; Goodman 2017) and a higherdensity X-ray model from our Picogna et al. (2021) analogue.

## 4.1. High ionisation parameter regime

Sellek et al. (2022) previously benchmarked the temperature structure of the EUV wind of Wang &amp; Goodman (2017) against mocassin . This work demonstrated that the wind temperatures were strongly overestimated by Wang &amp; Goodman (2017) compared to mocassin due to a lack of e ff ective optical lines available for cooling in their model: this could a ff ect the radii at which a wind can launch, as well as the dominant role for adiabatic cooling that they concluded, and the chemical composition of the wind. On the other hand, mocassin found higher wind temperature just below the wind base, likely as a result of it not including molecular cooling.

In Fig. 2, we compare the gas temperatures calculated by prizmo to the original temperatures of Wang &amp; Goodman (2017) and those calculated by Sellek et al. (2022). Since prizmo is designed to be the most comprehensive of the three codes - covering atomic cooling at low and high temperatures in neutral and ionised regions and also molecular cooling in warm dense gas - it is no surprise that everywhere it generally agrees with the cooler of the other two codes. Thus, in the atomic wind, optical lines help prizmo successfully avoids the overly high 2 - × 3 10 4 Ktemperatures predicted by Wang &amp; Goodman (2017) and thus agree with mocassin . Conversely, just below the wind base, molecular cooling in prizmo shrinks the large warm region at 10s au predicted by mocassin , instead limiting it to the inner ∼ 25 au as seen by Wang &amp; Goodman (2017).

Note that deeper into the disc where the dust and gas temperatures are coupled, prizmo and Wang &amp; Goodman (2017) continue to find even lower temperatures than mocassin , which cannot result from the additional cooling. This is since, unlike mocassin , prizmo and Wang &amp; Goodman (2017) do not calculate the di ff use radiation field that heats the dust, and hence underestimate the dust (and coupled gas) temperatures below the optical photosphere. Consequently, to avoid this problem, the temperatures Wang &amp; Goodman (2017) report - and which we plot in Fig. 2 - are set to be the maximum of that predicted from radial ray tracing alone and the temperature profile from Chiang &amp; Goldreich (1997), which reduces the extent of the underestimate. Therefore, prizmo shows the coolest temperatures; this would a ff ect the hydrostatic support of the disc and hence when using the temperatures for hydrodynamics we similarly have to pin the temperature as described by Eq. 3.

Overall, we conclude that prizmo performs well in the two extremes of the low-density winds with high ionisation parameter and in the dense underlying disc.

## 4.2. Lower ionisation parameter regime

Previous generations of X-ray photoevaporation models have generally relied on the ionisation parameter to control the relationship through a prescribed relationship. The temperatures after 50 orbits of our Picogna et al. (2021) analogue, which are calculated in this way, are shown in the leftmost panel of Fig. 3 - they show the typical appearance of an X-ray-driven wind: a nearly isothermal temperature of several 1000 K.

In the central panel of Fig. 3, we compare these to the temperatures calculated by post-processing this model using prizmo during the 'burn-in' step. These are vastly di ff erent, showing generally cooler temperatures and a strong radial dependence. In the right-hand panel, we confirm that this is not a particular issue with prizmo by showing that post-processing with another thermochemical code P ro D M i o (PROtoplanetary DIsk MOdel 4 , Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et al. 2016; Rab et al. 2018) shows the same qualitative features. We note that the P ro D M i o model uses a larger chemical network (100 species, see Kamp et al. 2017 for details) and includes heavy elements such as Fe or S, and more complete atomic line lists. Hence, the temperatures from P ro D M i o tend to be a bit cooler. Nevertheless, the overall agreement of PRIZMO with P ro D M i o is very good, indicating that the choices for the reduced chemical network and simplified heating cooling are well / justified.

The cooler temperatures in this regime suggest that there are cooling channels included in prizmo and P ro D M i o but not in mocassin , hence we now experiment with turning o ff di ff erent (combinations of) processes in prizmo to see what drives the difference. In this, we must keep in mind two questions with potentially separate answers: a) which cooling destabilises the temperatures from the 10 4 K temperatures predicted by the ξ -T relationship? and b) and which ends up dominating the cooling in the final equilibrium?

The most obvious candidate to make a di ff erence might be molecular cooling. However, since H2 and CO are readily photodissociated (and, unlike in later sections, are not replenished

4 https://prodimo.iwf.oeaw.ac.at rev.: 4fee3902 2023 06 19 / /

Fig. 2: Comparison of the gas temperatures predicted by (L-R) Wang &amp; Goodman (2017), mocassin (Sellek et al. 2022) and prizmo (this work) for the EUV-wind density model of Wang &amp; Goodman (2017). Inset panels show the inner 10 au.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-6a03916b-3771-4377-98c0-e5d30e69e4ec.jpg)

Fig. 3: Comparison of the gas temperatures predicted by (L-R) the ξ -T relationship from Ercolano et al. (2021) (c.f. Picogna et al. 2021), prizmo (this work) and P ro D M i o for the Picogna et al. (2021) analogue model. Inset panels show the inner 10 au.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-dfc21094-8a6c-4351-ac50-d626b9089bc0.jpg)

here by advection as this is a static benchmark), their abundances are not high enough to dominate the cooling. Thus, in panel b) of Fig. 4, turning o ff the molecular cooling does not raise the temperatures significantly save for in a limited region around 20 -30 au. Moreover, especially since thermal dissociation makes molecular survival at high temperatures even less likely, molecular cooling cannot be the process that destabilises the hot wind temperatures.

Another di ff erence is that the version of mocassin that was used in the works described above only includes collisional excitement by electrons, whereas in this work, we include electrons, H, H , + and both ortho- para-H2. / Oxygen is generally

Fig. 4: Comparison of the temperatures predicted prizmo when various combinations of cooling processes are switched o ff or reduced: "No Mol" means that molecular cooling due to H2 and CO are not included; "No O + H" means that collisional excitation of O by H is not included; "Low dust" means that the dust-to-gas mass ratio was lowered to 10 -10 . Inset panels show the inner 10 au.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-d3ba9e7a-e250-4518-8006-0c4f9b068349.jpg)

found to be a bright, dominant coolant and H will be the most abundant collider in a largely neutral atomic wind; indeed Ercolano &amp; Owen (2016) found that these collisions were necessary to match observed [O ] line luminosities. Therefore, in i panel c) of Fig. 4, we consider the e ff ect of switching o ff the collisions between O and H. We see that this has something more of an e ff ect in raising the temperature on scales of 10s-100s au, though still note much above 1000 K. This suggests that this is the process dominating the cooling in the final equilibrium, but alone is not destabilising the hottest temperatures. Switching o ff both molecular cooling and collisions between O and H is somewhat more powerful than either alone (Fig. 4, panel d). This is in line with the fact that H2 cooling is highly sensitive to the temperature at 100s K and will be boosted by the small rises; moreover, the decrease in the excitation of the oxygen lines due to fewer colliders limits the role of atomic cooling. The overall result is that H2 cooling is therefore what dominates at the 'new' equilibrium shown in panel c), and can cool the wind quite considerably but not quite so much as the collisional excitation of O by neutral H (and still cannot be destabilising the hottest temperatures). However, even without both processes, the temperatures are still lower than predicted from the ξ -T relation (Fig. 3).

Furthermore, our model includes a large amount of small dust; in the wind T dust &lt; T gas and hence it can have a cooling e ff ect. To see the e ff ect of this, we reduce the dust-to-gas mass ratio to a negligible 10 -10 . This has relatively little e ff ect on its own (Fig. 4, panel e) or in combination with switching o ff molecular cooling (Fig. 4, panel f). This further underscores the fact that collisional excitation of O by neutral H is dominating the cooling at the equilibrium achieved by the full prizmo model. Moreover, as argued above for O + H collisions, thermal accommodation with dust cannot be the sole process destabilising the hot temperatures. However, when depleted dust is used in combination with switching o ff collisions between O and H (Fig. 4, panel g), we find approximately isothermal several 1000 K temperatures comparable to those predicted by the Ercolano et al. (2021) ξ -T relation (very little additional di ff erence is made by also removing molecular cooling; Fig. 4, panel h). This confirms that both thermal accommodation with dust and collisional excitation of O by neutral H can destabilise the hot temperatures and both must be removed to prevent this, while cooling by H2 does not contribute to this picture and only dominates the cooling under conditions where it is cool enough to be abundant enough but warm enough to excite its emission lines.

To summarise,

- -either thermal accommodation with dust or collisional excitation of O by neutral H can destabilise the hot 10 4 K temperatures and both processes must be removed from prizmo to recover them; H2 cooling is irrelevant here as it would not survive at equilibrium at such high temperatures.
- -Collisional excitation of O by neutral H dominates the cooling at the thermal equilibrium reached in the full model; thus removing the other processes alone has little e ff ect on the temperatures.
- -Once the temperatures are cool enough for some H2 survival, it is the next most important coolant after collisional excita-

tion of O by neutral H and prevents the temperatures rising too significantly if only O + H collisions are removed.

Therefore, in conclusion, the cooler temperatures result for two combined reasons, one more fundamental, one more related to the modelling choices we make: while the dust abundance and grain size distribution could reasonably be such that dust is not a significant coolant (in this case panel e) would be most realistic), we cannot escape the fact that mocassin underestimated the atomic cooling in largely neutral media (where H is much more abundant than electrons) due to its omission of O + H collisions. O + Hcollisions are su ffi cient to lower the temperature in the intermediate density ionisation parameter conditions seen at / 10s au in previous X-ray photoevaporation models, and dominate the cooling in these static models so long as the process is included. Conversely, wherever the degree of ionisation is high, such as inside 10 -20 au where the X-ray flux is large (or in the lower-density EUV photoevaporation scenario; Sect. 4.1), the relative contribution of O + H collisions to the cooling rate is small, and the optical lines (including the Lyman lines of H) lead to the usual 10 4 K thermostat. The result is that lower densities would be needed to achieve the same temperatures seen in previous photoevaporation models, this could be visualised as translating the ξ -T relation to higher ionisation parameters (see Sect. 5.1.4).

## 5. Results

## 5.1. Structure of the wind

## 5.1.1. Temperature, density and velocity

Figure 5 shows the gas density (left) and temperature (centre) of our fiducial model after 200 orbits at 10 au; the white dotted line in each panel indicates the sonic surface. The top panels in Fig. 6 then show the temperature and density extracted along this sonic surface; we extract these values for r &gt; 6 2 au since inside this . radius there is no coherent transsonic flow. The inner streamline that crosses the sonic surface at 6 2 au launched from the disc at . 3 2 au at an elevation of 28 . ◦ . These values are consistent with the idea that a spherical, isothermal, 10 4 Kwind becomes transsonic at r ≳ 0 5 . rG ∼ 4 5 au (Parker 1958; Owen et al. 2012). While . some mass-loss may be driven from within by the enthalpy (Li ff -man 2003; Dullemond et al. 2007; Alexander et al. 2014), it is generally subsonic and here becomes supersonic at ∼ 6 au.

Qualitatively, the temperatures appear similar to those predicted by prizmo and P ro D M i o in Fig. 3 from post-processing the density field of our Picogna et al. (2021) analogue: there is a warm, vertically extended, inner region reaching 10 4 K, and declining temperatures outwards. Tracing along the sonic surface, Fig. 6 shows that the temperature is to a very good approximation isothermal in the inner 10-15 au before dropping suddenly. In the bottom left panel we see that this sudden drop occurs with the onset of the presence of trace amounts of molecular hydrogen in the wind, something we explore further later. Beyond this point - at around 17 au - out to around 120 au, the temperature then becomes a smoothly decreasing function of radius, to which we find a power law fit of T ∝ r -0 97 . . This is much steeper than typically found in previous X-ray photoevaporation simulations which showed essentially isothermal winds (Picogna et al. 2021), which we explore further in Sect. 5.1.4. It is also much steeper than the values found by Nakatani et al. (2018b). However, it is very close to the T esc ∝ r -1 escape temperature predicted by Owen et al. (2012) by assuming that nozzle flow e ff ects are negligible in the wind, and thus balancing the gravitational force and divergence terms in a manner similar to a spherical wind (Parker 1958). We overplot the values from Owen et al. (2012) and while they slightly overpredict the temperature at the sonic surface compared to the simulations by around 20 per cent, they are generally in good agreement. The di ff erence could be explained by the assumptions of a spherically diverging flow or the lack of centrifugal force.

The right-hand panel of Fig. 5 allows us to highlight the different regimes in dust and gas temperature and how these correspond to the pinning approach of Eq. 3. At columns ≲ 10 21 cm -2 , we are in the wind and T gas &gt; T dust . For a column density &gt; 10 21 cm -2 , the two temperatures are coupled and T gas ≈ T dust . However, since prizmo underpredicts the dust temperature below the photosphere, to avoid the loss of hydrostatic support to the gas, we smooth between the prizmo temperature and the DIAD temperature; the two diverge su ffi ciently little that the plotted gas and dust temperatures are the same down to a column density 10 23 -10 24 cm -2 , but start to diverge again beyond that, with T gas &gt; T dust since prizmo underestimates T dust .

We can also analyse the run of density in these regions. In the isothermal part of the wind, it scales as ρ ∝ r -1 01 . , similar to previous isothermal X-ray photoevaporation models (Picogna et al. 2021). When the temperature drops at ∼ 17 au, the density jumps in order that the pressure remains essentially smooth. It then declines again, though more shallowly, on average as ρ ∝ r -0 51 . . The peak number densities are therefore in the vicinity of 10 5 cm -3 , more than an order of magnitude lower than previous X-ray photoevaporation models where they peaked at ∼ 3 × 10 6 cm -3 . Even though similar 10 4 K temperatures are reached in these inner parts of the wind to previously, the additional cooling provided by collisional excitation of O by H means the wind can only reach them by having a much lower density (corresponding to ξ ≳ 10 -3 , compared to ξ ≳ 10 -4 previously; see Sect. 5.1.4)

The panels in Fig. 5 are also overlaid with streamlines at regular intervals. We start the streamlines at the location where the Bernoulli constant defined by

$$\underset { \substack { w h i a } { \text{infty} \\ \text{with the} } } { \text{infty} } \ B = \frac { u ^ { 2 } } { 2 } + \frac { \gamma } { \gamma - 1 } \frac { P } { \rho } - \frac { G M _ { * } } { r },$$

(for gas velocity u , adiabatic constant γ , pressure P , density ρ , stellar mass M ∗ and distance to the star r ), becomes positive. This is an approximation - assuming an adiabatic case - for the location that the gas becomes unbound, but generally matches well the point that the streamlines start to flow away from the disc rather than trailing along its surface.

Notably, the streamlines subtend an angle of no more than 30 ◦ with respect to the inclined wind base. This is much less than the 90 ◦ often assumed (e.g. Clarke &amp; Alexander 2016): when the temperature jump across the base is very strong, the gas is only accelerated perpendicular to the base (but not parallel to it) resulting in velocity vectors that are approximately normal to it. In this case however, the temperatures vary more smoothly from disc to wind so the perpendicular acceleration is less strong and the streamlines are more radial from the start. This means that the opening angle of the wind (the angle between the streamline and rotation axis) at the base is also larger, and ranges from 20 -35 . Notably these values are centred around 30 , the min-◦ ◦ imum opening angle for a magnetocentrifugal wind (Blandford &amp; Payne 1982). Therefore if such an angle is inferred from observations (e.g. Banzatti et al. 2019; Arulanantham et al. 2024), caution should be exercised in immediately interpreting it as a signature of an MHD wind.

$$t = 6 3 2 5 \ y r$$

Fig. 5: The gas density (left), gas temperature (centre) and ratio of gas and dust temperatures (right) after 200 orbits of the fiducial simulation. The latter is defined as min( T gas PRIZMO , , T DIAD) / T dust PRIZMO. , The solid purple line is the sonic surface, while the dashed lines show streamlines at regular intervals. The grey lines in the right-hand panel show total column densities of N = 10 21 , 10 23 , 10 25 cm -2 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-503836b8-56d9-40cd-95b4-9f4a18292ff1.jpg)

## 5.1.2. Composition of the wind

The left-hand panel of Fig. 7 shows a composite colour image representing the major H-carrying species where red represents H , green neutral H and blue H2. As one would expect, the un-+ derlying disc is strongly molecular, while the innermost parts of the wind can become highly ionised due to EUV photons. What is perhaps more surprising is the amount of H2 present in the outer parts of the wind. Figure 6 shows that its abundance rapidly grows to around 10 per cent around ∼ 17 au and then steadily grows at the expense of neutral H until it becomes dominant outside of 50 -60 au. Indeed, there is more H2 than at the end of the burn-in step, despite the lower density wind conditions. This is a surprise since H2 can be photodissociated by FUV &amp; EUV photons (while we do not have an explicit significant FUV component in our spectrum, there is still some contribution from the photosphere) . Moreover, our network includes 5 the following reactions resulting in the thermal dissociation of H2 into atoms ions: /

at the outer edge of the wind. However, the self-shielding of H2 can help protect it once it accumulates a column of ≳ 10 14 cm -2 (Draine &amp; Bertoldi 1996; Richings et al. 2014). Given the approximate density of 10 5 cm -3 and a molecular fraction of 10 -6 , inside the H-H2 transition, the column density can only be N ≈ 10 -6 × 10 5 cm -3 × 10 au = 10 13 cm -2 . Only once the wind reaches a molecular fraction of 10 -1 outside the H-H2 transition, can a column density of N ≈ 10 14 cm -2 be reached in the space of only 10 -3 au. Hence, self-shielding can contribute to the steady increase in H2 in the wind at 10s au, but not the initial onset.

$$H _ { 2 } + H \to 3 H & & ( 5 ) ^ { \text{ would } n } _ { \text{c:d:d:d} }$$

$$H _ { 2 } + e ^ { - } & \to 2 H + e ^ { - } & ( 6 ) \begin{matrix} \sin \\ \ e a u i l i b r i \end{matrix}$$

$$H _ { 2 } + H _ { 2 } & \to 2 H + H _ { 2 } & ( 7 ) \quad \text{process}$$

$$H _ { 2 } + H e ^ { + } & \to H + H ^ { + } + H e. & ( 8 ) \stackrel { \text{years with} } { \underset { \text{$\text{$\dots$} } { \text{$\dots$} } }$$

Given the wind densities suggested by our simulations, the wind will be optically thin to the photodissociating radiation. By integrating the photodissociation cross-section across our input spectrum, we find that the photodissociation rate per H2 molecule is k PD = 3 7 . × 10 -7 s -1 L FUV 10 30 erg s -1 GLYPH&lt;16&gt; r au GLYPH&lt;17&gt; -2 , corresponding to an expected lifetime of around 25 years at 17 au or 2800 years

5 H2 can also be destroyed by X-ray photodissociation, but this is not included in the current version of PRIZMO.

Moreover, for thermal dissociation, given the abundances in the inner disc and the 10 4 K temperatures involved, we expect neutral hydrogen to be the most important destruction route. This provides a thermal dissociation rate per H2 molecule of k th = 5 2 . × 10 -11 cm 3 s -1 n H; given the gas densities ∼ 10 5 cm -3 , which will be dominated by neutral H, a thermal dissociation rate of k th ≲ 6 × 10 -6 could be produced. This is a factor of a few higher than the photodissociation rate, so should shift the transition even further outwards, and further means that self-shielding would not protect the H2.

Since self-shielding cannot explain the onset of H2, nonequilibrium chemistry must therefore play a role. We postprocess the simulation with prizmo for an additional ∼ 6500 years with hydrodynamics switched o ff (Fig. 8, left); we then expect that multiple photodissociation timescales will have elapsed everywhere and the H2 fraction will have more or less reached equilibrium. Comparing to the left-hand panel of Fig. 7, we can see that the colours in the wind shift away from dark and light blues towards light blues and greens - and the border between atomic and molecular winds moves outwards - indeed demonstrating an increase in neutral H as H2 is photodissociated and confirming that the molecular abundances are out of equilibrium.

We now estimate the hydrodynamic timescale to show that advection of H2 into the wind is responsible for its survival. Ma-

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-0c55604b-256f-49e3-98a7-25900286d012.jpg)

r

/ au

Fig. 6: Quantities extracted along the sonic surface. We only show r &gt; 6 au (where r is the spherical radius) since there is no coherent transsonic flow interior to 6 au. Top left: the gas (blue) and dust (orange) temperatures. For reference the dashed line shows the escape temperature of Owen et al. (2012) and the dotted lines are fits to the hot isothermal section, molecular section and outer section. Top right: the gas density, with fits to the same sections. Bottom left: the proportion of H in H + (red), H (green) and H2 (blue). Bottom right: the fraction of the thermochemical cooling given by each major category of cooling.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-0f05afc5-b00a-4685-a6b4-b5f6977e398f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-8429e44c-ed06-4a63-aaca-a7c63f2355ca.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-b19ac65a-4c48-487c-b890-638f5ac5b5fc.jpg)

terial at each point on the sonic surface has survived to a distance r and is moving at speed c S( r ), where, as shown above, this can be approximated as the sound-speed at the escape temperature c S( r ) = q GM ∗ 2 r (Owen et al. 2012). This corresponds to a rate k ad ∼ c S / r ≈ 1 4 . × 10 -7 GLYPH&lt;16&gt; r au GLYPH&lt;17&gt; -1 5 . . Realistically, this is an upper limit, since the gas is initially moving subsonically. In the right-hand panel of Fig. 8, we also compare the rate at which H2 is advected into the wind in our models to the rate at which it forms on grain surfaces from purely atomic gas (the upper limit on its formation rate). This further supports our above estimate

that advection is the dominant source of H2 in the wind on scales of 10s au. Note that we also have a rather large dust-to-gas ratio in the wind and lowering this may further reduce the in-situ production of H2 relative to its advection. Overall we conclude that non-equilibrium chemistry is important for the composition of the wind.

Having concluded that advection is responsible for the presence of H2 in the wind, in Fig. 9, we now more carefully compare our estimate for the advective resupply timescale to the above estimate for the photodissociation rate and the ex-

$$t = 6 3 2 5 \ y r$$

Fig. 7: Three-colour images of the major carriers of hydrogen (left), carbon (centre) and oxygen (right). In each case, red colours represent ionised species (H + / C + / O ), green neutral atomic species (H C O) and blue molecular species (H2 + / / / CO vapour + ice H2O / vapour + ice). Therefore darker blues indicate highly molecular gas while light blues indicate gas with a significant mix of atomic and molecular species. The solid purple line in the sonic surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-93486e9c-3fd9-44e6-b3d0-b87f1609e011.jpg)

Fig. 8: E ff ects of non-equilibrium (thermo)chemistry. Left: Composite colour image as per the right hand panel of Fig. 5, after 6507 yr of additional evolution without hydrodynamics. Centre: Estimate of where the rate of advection of H2 into the wind dominates over its in situ formation on dust grains. Right: The adiabatic cooling as a fraction of the total thermochemical cooling (N.B. white regions are where PdV work instead heats the gas). In each panel, the solid purple line in the sonic surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-8c8e5aa5-a465-4cb7-bc84-51de55c48f24.jpg)

pected thermal dissociation rate assuming T r ( ) = T esc ( r ) 6 and n H( ) r ≈ 10 5 5 . ( r / au) -1 , in order to understand where the location of the H-H2 transition. Equating the advection and photodissoci- ation rates, we can demonstrate analytically that H2 would survive against photodissociation alone at

6 While inside the onset of the molecular survival, the gas is actually fairly isothermal at T ∼ 10 4 K, if the H 2 did survive it would be expected to stabilise T r ( ) ∼ T esc ( r ) (see later discussion), hence this is a reasonable scenario for explaining why, under those conditions H2 wouldn't survive if it were present (a situation which is unstable since if destroying H2 raises the temperature, it further increases the destruction rate).

$$\text{vacuum} \quad vivec \, \text{agam\,p\nv\nucuss\nucuss\nolum\,a\nomega\,a} \\ \text{dissoci} \quad \frac { r _ { H 2 } } { \text{au} } & = \left ( \frac { 3. 7 \times 1 0 ^ { - 7 } } { 1. 4 \times 1 0 ^ { - 7 } } \right ) ^ { 2 } \left ( \frac { L _ { F U V } } { 1 0 ^ { 3 0 } \, \text{erg} \, s ^ { - 1 } } \right ) ^ { 2 } \\ \text{s actually} \quad \text{expected} \quad \text{$=7\left(\frac{L_{F U V }}{1 0 ^ { 3 0 } \, \text{erg} \, s ^ { - 1 } } \right ) ^ { 2 }, \\ \text{would} ^ { t } \, \text{$\dots$} \quad \text{ac\,can\,in\,Eq\,0\,\text{Allow}$ at this location} \text{ } \text{thormal\,dicecio} \quad.$$

as seen in Fig. 9. However, at this location, thermal dissociation is dominant and so the transition does not occur here but

Fig. 9: The rates (per H2 molecule) at which H2 is supplied to the wind by advection (black), H2 is photodissociated (red, dotted) and H2 is thermally dissociated by H (red, dot-dashed). The solid red line represents the total dissociation rate given by the sum of the other two, while the vertical grey lines mark the position where H2 survives against each dissociation rate.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-c4723e83-e31b-447d-b062-e63002588712.jpg)

it pushed out somewhat further. The estimated thermal dissociation rate is exponentially sensitive to temperature and therefore a much stronger function of radius than the other two rates. It falls below the photodissociation rate at ∼ 16 au; shortly thereafter at ∼ 17 au, the total dissociation rate falls below the advection rate, which is in excellent agreement with where we see the onset of H2 in the simulation (Fig. 6) While photodissociation here dominates the destruction of H2 at its onset, albeit mildly, we note that this model has a fairly low FUV flux, and that raising this could change substantially shift the balance of destruction further towards photodissociation and push the onset of H2 even further out (Eq. 9).

Finally, in Fig. 7, we also show similar colour maps for C and O as for H. Here the blue regions are overall more confined to the disc, showing that molecules are less important carriers of these atoms. For example, CO is the most important molecular carrier of C, but since it is photodissociated by FUV (at a rate around 13 times faster per molecule than H2) - and less abundant than H2 and so less able to self-shield at large radii - we only see trace quantities in the wind. This is consistent with the results of Panoglou et al. (2012) whose model of MHD winds suggested CO could be abundant in winds during the Class 0 phase, but highly depleted relative to H2 by the Class II phase. For O, H2O ice becomes a surprisingly larger contributor at large radii where the dust temperatures fall below the 150 K freeze-out temperature. In our models this happens because any water that forms freezes out and acts as a sink of O although in reality photodesorption will likely liberate the molecules which will then rapidly photodissociate (if they do not already do so in the solid state); we aim to include such processes in future work to ensure we can accurately study the outer disc. Otherwise, O is mainly neutral and atomic, and shows the same ionisation pattern as H as expected due to their near identical ionisation potentials. In contrast, C is overall somewhat more ionised since its ionisation energy is lower and thus it can be e ff ectively ionised by highly penetrating FUV photons

## 5.1.3. Main heating &amp; cooling processes

The final panel of Fig. 6 shows the fractional contribution of each process to the total thermochemical cooling. In the inner parts of the wind, there are few molecules, and the atomic lines dominate the cooling. Once significant ≳ 10 per cent, H2 enters the wind, at 20-30 au, it dominates the cooling with no more than 1 per cent coming from other processes. However, we see that despite H2 becoming the increasingly dominant carrier of H as we move outwards, the atomic cooling slowly rises, taking over again beyond ∼ 120 au. This occurs because the H2 cooling is strongly temperature dependent and beyond a certain radius where the temperature falls below ∼ 600 K, the mid-IR lines are no longer excited. On the other hand, several atomic emission lines are present at far-IR wavelengths and so are still excited at lower temperatures. However, the lower temperatures result in lower overall cooling. There is also a narrow region around the onset of the molecular wind where 'chemical' cooling - consisting of the energy lost to recombination lines and collisional dissociation of H2 - dominates the cooling. This is because this is the region where the wind crosses from being cooler and significantly molecular, to being hotter and fully atomic - streamlines of material do cross this border and thus H2 suddenly experiences hotter gas where it it rapidly dissociated.

The ratio of the adiabatic cooling to this total is also plotted in Fig. 6. Apart from in the narrow transition region around the onset of the molecular wind, adiabatic cooling generally contributes no more than 10 per cent. This is less than found by, for example, Wang &amp; Goodman (2017) because lower temperatures and higher densities typically promote two-body, line cooling processes (which are quadratic in density) at the expense of adiabatic cooling (which is only linear in density) (Sellek et al. 2022). Moreover, Wang &amp; Goodman (2017) excluded optical emission lines and so the only active thermochemical cooling in the wind was recombination, which we do find to be everywhere less significant than adiabatic cooling in line with their results.

The strong contribution of adiabatic cooling in the transition region occurs because gas crosses from a cooler, more molecular region where significant H2 cooling is possible, to where the wind is completely atomic, H2 becomes too rare to contribute to the cooling, and the temperatures rise. The gas dilutes strongly as the temperature rises - to conserve mass it must therefore also accelerate strongly, which leads to a localised increase in adiabatic cooling.

At very large radii, where line cooling becomes less e ff ective at temperatures T &lt; hc k λ , adiabatic cooling becomes relatively more important, rising to around 30 per cent of the thermochemical cooling rate. Again, this must imply stronger acceleration, which brings the sonic surface increasingly close to the base. This can be seen in the panels of Fig. 5 and means that going to large radius, the sonic surface is not a line of constant latitude but probes denser, cooler, gas, producing the increase in ρ and sharper decrease in T seen beyond 120 au in Fig. 5.

The role of adiabatic cooling is also shown across the whole of the wind in the right-hand panel of Fig. 8. The same features can be seen: the steady increase towards larger radii accompanying the drop in the sonic surface towards the base amd the narrow strip of significance across the molecular-atomic transition running diagonally across the the inset. The latter appears to wider out at higher altitudes to cover most of the more highly-ionised parts of the flow (see for example the orange region in the carbon carriers in Fig. 7) suggesting that this region is characterised by a relatively strong acceleration and poor radiative cooling.

## 5.1.4. Explaining differences to previous models

Previous X-ray photoevaporation models calculated the gas temperatures using pre-calculated functions of the ionisation parameter ξ = LX nr 2 . In order to understand the salient di ff erences in our models from such works - the lower densities and the stronger temperature gradient in Fig. 10 we plot a 2D histogram of the temperature and ionisation parameter for our Picogna et al. (2021) analogue (which uses the ξ -T relationship of Ercolano et al. 2021) and our fiducial model. We include only cells where T gas . T DIAD: for the fiducial model, this is everything with N &lt; 10 25 cm -2 (Eq. 3) and for the analogue model those with N &lt; 2 × 10 22 cm -2 and ξ &gt; 10 -8 . In both cases, we observe a tight correlation between ξ and T for ξ ≳ 10 -6 ; for the fiducial model somewhat more scatter occurs above the main trend for ξ ≲ 10 -6 due to cells with 10 21 &lt; N / cm -2 &lt; 10 25 which are interpolated between T PRIZMO and T DIAD &gt; T PRIZMO.

For the Picogna et al. (2021) analogue model, there are very few points in the regime 500 -2000 K, which separates the cold disc from the hot wind. Since the ξ -T relationship is quite flat in the hot wind, it occupies a narrow range of temperatures and cannot achieve a significant temperature gradient.

We can understand this behaviour by considering that in order to accelerate a wind (and also in the underlying hydrostatic disc), we expect pressure to be decreasing upwards,

$$\begin{pmatrix} \partial \ln P \\ \partial \theta \end{pmatrix} _ { r } > 0. & \frac { \text{achieve} \frac { 1 } { 2 } } { ( 1 0 ) } & \text{that the t} \\ & \text{depends}$$

Moreover, given that the wind material accelerates upwards, we also expect an upwardly decreasing density

$$\left ( \frac { \partial \ln \rho } { \partial \theta } \right ) _ { r } > 0. & & ( 1 1 ) \quad \text{theated $\xi$} \\ & \quad \text{fribes}.$$

Combining these conditions with the ideal gas law leads to the following requirement

$$\left ( \frac { \partial \ln T } { \partial \ln n } \right ) _ { r } > - 1, & & ( 1 2 ) & \xi \sim 1 0 ^ { - 1 } \\ & \quad \text{be a pos} \\ \text{More}$$

which, given the definition of ξ , we can also express as

$$\left ( \frac { \partial \ln T } { \partial \ln \xi } \right ) _ { r } & < 1. & ( 1 3 ) \quad & \text{disc win} \\ & \quad & \text{in} \quad \text{the} \\ \quad & \text{under} \quad \text{$\dots\mu$}.$$

Thus, we see that if the temperature increases too steeply with increasing ξ then it leads to an increasing pressure, which would oppose the flow. For the ξ -T relationship at each column from Ercolano et al. (2021), we calculate the upper and lower points where this condition is met, which we indicate with the white dash-dotted lines on Fig. 10 - we can see the upper set is in good agreement with the deficit of points separating the locations of the cold disc and hot wind (while there are no points below the lower set since the cold disc is generally found at lower ξ ).

For our fiducial model, the ξ -T relationship in Fig. 10 is shifted to the right with respect to those calculated by Ercolano et al. (2021) - this implies that for the same temperature to be reached, higher ξ - or equivalently lower densities - are needed. This happens, in part, due to the more e ffi cient atomic cooling resulting from the excitation of O by H, and is seen in models even without molecular cooling. As discussed, this is important for explaining the lower density achieved in the inner 10 4 K part of the wind (where H2 does not survive and cannot contribute to the cooling).

Moreover, the relationship is much shallower than before. The shape of the ξ -T relationship is determined by the shape of the cooling function. For calculations with mocassin there is a lot of cooling from FIR atomic lines at low temperatures, and from optical atomic lines at high temperatures. This means that there is an intermediate temperature range where the excitation of the FIR lines has saturated, while that of the optical lines has not yet become e ffi cient. In this regime, even a small amount of additional heating cannot be o ff set by additional cooling, without a large rise in the temperature in order to activate the optical lines, leading to a steep relationship between ξ and T. Once the optical lines are e ffi ciently excited, a small perturbation in temperature results in a large change in cooling and many di ff erent heating rates can be balanced in a narrow range of temperatures, leading to the roughly isothermal upper branch.

The introduction of molecular cooling, along with the presence of H2 in the wind due to its advection, results in a large number of lines at mid-IR wavelengths. These are able to provide cooling in that previously unstable branch, while not being so strong as to prevent additional rises in temperature entirely, thus leading to a much shallower relationship between ξ and T, with maximum gradient ∂ ln( T ) ∂ ln( ξ ) = 0 79. This means that inequal-. ity 13 is now satisfied at these temperatures, and allows them to be populated by wind material. Therefore, the full 1 / r range of escape temperatures predicted by Owen et al. (2012) can now be achieved. In turn, at the radii 20 -120 au where this occurs, given that the temperature profile is fixed, the radial run of density then depends only on the slope of the ξ -T profile.

We note that there is a similarity to the concept of thermal sweeping (Owen et al. 2013), which Haworth et al. (2016) showed happens if there is a maximum pressure along the X-ray heated ξ -T curve that exceeds the pressure in the disc midplane. In this case, such a scenario is unlikely as the shallow ξ -T profiles suggest that the X-ray heated gas can never exceed the pressure of the disc gas. At the very least, a very large ( ∼ 100 × ) depletion in the disc gas would be necessary by the time the wind opens a cavity in order to move the disc points from ξ ∼ 10 -7 to ξ ∼ 10 -5 for thermal sweeping by hot X-ray-heated gas to then be a possible solution to the relic disc problem.

Moreover, while the microphysics of the heating and cooling are substantially di ff erent, a similar e ff ect is known for thermal disc winds from X-rays binaries, where thermal equilibrium stability curves may be employed to find the maximum and minimum temperatures of the cool and hot branches respectively, with the stable wind conditions being those on the hot branch (e.g. Higginbottom &amp; Proga 2015; Higginbottom et al. 2017).

The retention of a tight correlation between ξ and T is somewhat to be expected, since like, the case the atomic cooling previously assumed in the ξ -T relationships, the H2 cooling happens via two-body processes. Balancing a two-body cooling rate Γ = n H2 n coll f ( T ) with the optically-thin X-ray photoionisation heating rate Λ = n ¯ σ L X 4 π r 2 for some spectrally-averaged cross-section ¯ σ allows one to find the temperature as T = f -1 ( ¯ σ 4 π L X nr 2 ) = f -1 ( ¯ σ 4 π ξ ) assuming that n H2 , n coll ∼ n . For a fixed spectrum and gas composition, the temperature then depends only on ξ .

However, since the gas composition realistically varies throughout the wind, achieving a ξ -T relationship also relies on the composition (i.e., in this case, n H2 / n ) varying with ξ . This is generally true for the photoionisation level (e.g. Glassgold et al. 2007), but is less obvious for the H2 abundance as a result of the balance of photodissociation and advection for two reasons. Firstly, the photodissociation is driven mainly by UV, not X-ray, photons i.e. the overall dependence is on both L FUV and L X; in

0

Fig. 10: The distribution of cells in the plane of gas temperature T and ionisation parameter ξ for the Picogna et al. (2021) analogue model (left) and our fiducial model (right). We only show cells which are not fixed to the initial T DIAD. Overlaid are various parametrizations of the relationships between the variables: on the left, the yellow dash-dotted lines are the parametrizations for three di ff erent column densities from Ercolano et al. (2021); on the right, the yellow dash-dotted line is the lowest column density Ercolano et al. (2021) parametrization, the orange dashed line is from Haworth et al. (2016) and the white pink solid lines are from / this work. The other white lines delineate the unstable region in the Ercolano et al. (2021) prescription (dash-dotted) and lines of constant pressure at a given radius (dotted).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-aa31a588-4233-4801-a89a-38ebd9044f3f.jpg)

this case, since the radiation field is fixed the two are directly correlated, but some caution should be exercised extrapolating to a di ff erent spectrum. Secondly, both photodissociation and advection are one-body processes, depending only on radiation field and r as per Sect. 5.1.2 which would suggest a relationship independent of n . However, at larger columns the self-shielding of H2 introduces a dependence of N k -(where k = 3 4 is often / assumed, e.g. Draine &amp; Bertoldi 1996) and thus, given the shallow density gradients seen here, a dependence on the local value of n . Balancing this with the advection rate described earlier, a dependence of L FUV ( nr 5 3 / ) 3 4 / ∼ L FUV L X ξ 3 4 / L 1 4 / X r 1 4 / may result. While not solely dependent on ξ , given the fixed L X and L FUV / L X and weak dependence on r , we expect the composition to track ξ closely enough to support the tight correlation seen.

Given that a tight correlation between ξ and T is seen despite these caveats, we can produce a ξ -T fit to our simulations. We explore two fits: the first on all cells with a column density N ≤ 10 25 cm -2 , downweighting those with N ≥ 10 21 cm -2 according to the weighting scheme used to interpolate between the prizmo and DIAD temperatures (Eq. 3); the second where we only fit the subset of those cells with T gas &gt; T dust , so that we only include cells whose thermochemistry is dominated by the processes discussed above that lead to a ξ -T relationship, rather than dust heating. On the other hand, for practical usage, the first fit may more smoothly transition to the dust temperatures. We

Table 2: Parameters derived for the fits to the ξ -T relationship

| Fit Parameter   |   Value (All cells) |   Value ( T gas > T dust ) |
|-----------------|---------------------|----------------------------|
| a               |            1.04645  |                   1.19601  |
| b               |           -2.802    |                  -2.89746  |
| c               |           -0.145929 |                  -0.165763 |
| d               |            4.03277  |                   4.04502  |
| m               |        11654.6      |               10373.5      |

Notes. See Eq. 14 for functional form used in fits use the following fitting function (Ercolano et al. 2021):

$$\cdot \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \quad \quad \quad$$

The best fitting values of the parameters are given in Table 2 and show the fits on Fig. 10. Practically, little di ff erence is found between the two fits except at the coolest end. We briefly discuss the usefulness of such a fit in Sect. 6.1 but reserve testing whether it enables a return to a less-expensive method of determining wind temperatures for future work.

## 5.2. Mass-loss rates

The single key quantity from the simulations controlling the importance of photoevaporation and the lifetime of the disc is the

total mass-loss rate. In each snapshot, we calculate the cumulative mass-loss rates as a function of radius by summing the mass flux through the spherical surfaces at di ff erent radii, counting only cells where the Bernoulli constant (Wang &amp; Goodman 2017) is positive and hence the material is unbound. Note that if the chosen radius is too close to the simulation outer boundary (400 au) then there may be oscillations caused by sound waves reflected o ff of the boundary (Wang &amp; Goodman 2017; Picogna et al. 2019). To provide a controlled measure, and for ease of comparison to the wider literature (Komaki et al. 2021), in Fig. 11a we focus on values integrated out to 80 au.

The mass-loss rate of our fiducial model settles down after just under 4000 years (or around 120 orbits at 10 au) of evolution to a value of M ( &lt; 80 au) = 1 06 . × 10 -9 M yr ⊙ -1 . This is just over an order of magnitude lower than the Picogna et al. (2021) analogue simulation which had settled to a mass-loss rate of M ( &lt; 80 au) = 1 86 . × 10 -8 M yr ⊙ -1 after 50 orbits. This suggests that the previous typical X-ray mass-loss rates of ≳ 10 -8 M yr ⊙ -1 (Owen et al. 2010, 2011a, 2012; Picogna et al. 2019; Ercolano et al. 2021; Picogna et al. 2021) have been overestimates, consistent with evidence from population synthesis. However, despite several di ff erences in methodology - particularly regarding the illuminating spectrum - the mass-loss rates are now within a factor 2 of those calculated for a solar-mass star by Komaki et al. (2021), who also performed on-the-fly thermochemistry with hydrodynamics. This is an encouraging step towards a convergence of mass-loss rate estimates between works and in stark contrast to disagreement by several orders of magnitude seen previously.

We note however, that in our model there is still significant mass-loss beyond 80 au: fitting out to the observed drop-o ff in the surface density profile at ∼ 160 au suggests a total value 3-4 times higher at ˙ M tot = 4 3178 . × 10 -9 M yr ⊙ -1 . This is still ∼ 9 times lower than the total rate found by Picogna et al. (2021) for the same stellar mass and spectrum.

## 5.2.1. Dependence on spectrum

Due to the strong dependence of the cross-section for photoionisation on photon energy, X-ray photoionisation heating is very sensitive to the shape of the X-ray spectrum (Sellek et al. 2022). The results shown here - for example, the order of magnitude of the mass-loss rate and the presence of a molecular photoevaporative wind - are qualitatively similar to those presented in Nakatani et al. (2018b); Komaki et al. (2021). However, these works argued that the contribution of X-rays to driving the wind was negligible, unlike the results of our simulations. Therefore we consider if the di ff erent shapes of the X-ray spectra employed could be driving the di ff erence.

The TW Hya spectrum (Nomura et al. 2007) used by Nakatani et al. (2018b); Komaki et al. (2021) (see Fig. 1) clearly lacks significant amounts of flux in the soft X-rays compared to that used by Ercolano et al. (2021), and the EUV is almost entirely absent. This is despite us normalising the two spectra to the same total X-ray luminosity (2 × 10 30 erg s -1 between 0 5 . -5 keV) and the fact that TW Hya has a notably 'soft' Xray spectrum (Kastner et al. 2002) compared to other T Tauri stars (which here manifests as the faster drop-o ff at E &gt; 1 keV due to the lack of the highest-temperature plasma component; Nomura et al. 2007). The di ff erence at the soft end arises because in X-ray spectra as observed at the Earth, these energies are attenuated due to absorption (primarily by neutral hydrogen)

in the foreground interstellar medium along the line of sight . 7 The Nomura et al. (2007) spectrum is a model for this observed flux - and thus includes attenuation by their estimated column of N H = 2 7 . × 10 20 cm -2 - while the Ercolano et al. (2021) spectrum is a model for the intrinsic emission of the star - and is corrected for this attenuation. The latter is what the wind receives and is thus the more physical model.

Figure 11a shows that the impact of using the observed spectrum instead of an intrinsic spectrum is to lower the mass-loss rate by a factor 2-3. This indeed suggests that the mass loss is driven by soft X-rays as usually argued (Ercolano et al. 2009). Calculating the ratio of the luminosites of the two spectra integrated from 100 eV to di ff erent upper limits E up (Fig. 12), we find that this ratio equals the ratio of mass-loss rates for E up ≈ 400 eV. Under the typical assumption that the mass-loss rates are approximately linear with the appropriate X-ray luminosity (Owen et al. 2012), this suggests that the photoevaporation in our simulations is being driven by X-rays with energies 100 -400 eV.

## 5.2.2. Dependence on microphysics

Since our fiducial spectrum is identical to that used by Ercolano et al. (2021); Picogna et al. (2021), this cannot be a reason for the di ff erence in mass-loss rates between this work and the previous simulations.

Instead, the same microphysics that results in the drastically di ff erent temperature profile can also explain the lower mass-loss rates. Figure 11b shows the mass-loss rate for models where we run the hydrodynamics for a short time with either the dust-togas ratio reduced to a negligible value and the collisions between neutral H and O disabled - since this was a particularly e ff ective combination for changing the temperature in our benchmarking -or the cooling from H2 excluded - since this dominates over most of the wind in our fiducial model.

In the latter case we see that neglecting H2 cooling boosts the mass-loss rate by only a factor ∼ 2. While important, this cannot fully explain the di ff erence to previous models. Instead, in the former case, we find a much more dramatic di ff erence, with the mass-loss rate being more than an order of magnitude higher. As shown in Fig. 3, this combination of reduced dust and no excitation of O by neutral H significantly reduces the cooling, which there resulted in a hot, approximately isothermal wind to persist. Compared to our fiducial model then, gas of a given density can reach a high temperature, meaning that the escape temperature can be reached for higher density gas, leading to a higher density wind. The fact that these high mass-loss rates are even higher than those seen in previous X-ray simulations likely results from weaker atomic cooling due to lower collisional excitation rates by electrons (Bell et al. 1998), as well as the smaller collections of lines we included (for example by neglecting S due to its likely strong depletion).

Therefore, the lack of dust and neutral colliders in mocassin are the most likely reason for the high-mass loss rates seen in previous X-ray photoevaporation simulations. Since our dust-togas ratio is likely too high (see Sect. 6.2), then the latter of these is the more important reason. However, the advective survival of H2 then allows it to become the most important coolant limiting the mass-loss rates, something that could not be captured without the on-the-fly thermochemistry.

7 Although in principle, there could be some contribution from the circumstellar material - including any wind - TW Hya is close to pole-on so this contribution should be negligible.

Fig. 11: Panel a: The evolution of the mass-loss rate inside 80 au calculated using pluto+prizmo with the X-ray spectrum of Ercolano et al. (2021, E21) (blue; 200 orbits at 10 au) and Nakatani et al. (2018b, N18) (green; 150 orbits at 10 au). For comparison, we also include 50 orbits of an analogue of the Picogna et al. (2021, P21) model (which uses a ξ -T relationship derived from mocassin ) as the black dashed line, and the value at 1 M ⊙ from the fit to ˙ ( M M ∗ ) given in Komaki et al. (2021) as the grey dotted line. Panel b: A zoom in on the evolution of the mass-loss rate inside 80 au during the first 50 orbits (at 10 au) of our fiducial simulation with the X-ray spectrum of Ercolano et al. (2021, E21) (dark blue) compared with models where we turn o ff the cooling due to collisions between O and H and include negligible dust (turquoise) or turn o ff the cooling due to H2 (light blue). For comparison, we again include 50 orbits of the Picogna et al. (2021, P21) analogue simulation as the black dashed line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-5b5d530d-0b64-414a-a23a-7832045abf04.jpg)

## 5.2.3. Mass-loss profile

Table 3: Parameters derived for the fit to the mass-loss profile

In Fig. 13 we show fits to the mass-loss profile in our fiducial simulations taken every 50 orbits at 10 au. We determine these in the same way as before: by calculating the cumulative massloss profile, ˙ ( M &lt; r ), fitting a sixth-order polynomial in log-log space and di ff erentiating:

$$\cdot & \quad \succeq \\ \log \left ( \frac { \dot { M } ( < R ) } { \dot { M } _ { \text{tot} } } \right ) = a \log ( R ) ^ { 6 } + b \log ( R ) ^ { 5 } + c \log ( R ) ^ { 4 } + d \log ( R ) ^ { 3 } \\ & + e \log ( R ) ^ { 2 } + f \log ( R ) + g,$$

$$\underbrace { \dot { \Sigma } } = \frac { 1 } { 2 \pi R } \frac { \partial \dot { M } ( < R ) } { \partial R }. \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The parameters of the fit for the final snapshot are given in Table 3; ˙ M tot is measured out to the point where the fit predicts a negative mass-loss rate per unit area (which we use as the outer edge of the wind) - and is the value that should be used in the fit - while ˙ M true ( &lt; 80 au) is the equivalent measurement out to 80 au, and ˙ M fit ( &lt; 80 au) is the value at 80 au predicted by Eq. 15.

As the simulation evolves, most of the change can be seen in the outer regions of the wind, while the innermost parts are unchanged (due to the freezing of the simulation inside r ≈ 17 au).

Notes. The form of the profile is described by Eq. 15. The bottom two values give, for reference, the cumulative mass loss rate inside 80 au predicted from the fit and that measured directly in the simulation (as used in Sect. 5.2). We recommend using the full precision given here since the parameters are fits to exponents which can result in truncation errors compounding to give a poor fit at large radii.

| Fit Parameter                        | Value                                   |
|--------------------------------------|-----------------------------------------|
| ˙ M tot a                            | 4 . 3178 × 10 - 9 M ⊙ yr - 1 - 1 . 2108 |
| M fit ( < 80 au) ˙ M true ( < 80 au) | 1 . 50 × 10 - 9 M ⊙ yr - 1 - 9 - 1      |
| b                                    | 9 . 6815                                |
| c                                    | - 31 . 226                              |
| d                                    | 52 . 293                                |
| e                                    | - 48 . 475                              |
| f                                    | 25 . 093                                |
| g                                    | - 7 . 6443                              |
| ˙                                    |                                         |
|                                      | 1 . 06 × 10 M ⊙ yr                      |

At the earliest time shown the outer parts of the simulation have not settled, and so an unphysical increase towards the outer boundary is seen. As the outer regions settle, the outer radius decreases to ∼ 160 au up to around 150 orbits, before increas-

Fig. 12: The ratio of the soft X-ray luminosity of the two X-ray spectra investigated in this work, as a function of the maximum energy deemed to be soft. The black dashed line shows the ratio of the two mass-loss rates for comparison.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-b9bc30f2-5ef4-49e6-9301-2900f8362998.jpg)

ing slightly again to 180 au by the final plotted snapshot. This suggests that the outer regions of the simulation have mostly stabilised, although running for a bit longer may demonstrate additional mass-loss from the outermost parts as the disc settles.

Aside from maybe at its very innermost edge, the profile is everywhere much lower than seen in the previous models by Picogna et al. (2021) - the additional microphysics included here compared to previously lowers the mass-loss rate from the disc everywhere. The slope of the profile is also a little shallower than previously due to the much less steep dependence of ρ on r. The inward shift of the peak and inner edge is likely due to slightly higher maximum temperatures in the wind, resulting in unbound material closer to the star. These temperatures are reached since the lower column densities are lower than the lowest slabs considered by Ercolano et al. (2021); Picogna et al. (2021).

When the dust-to-gas ratio is turned down to a negligible value and the collisions between neutral H and O are disabled, the mass-loss rate is increased by roughly the same factor ≳ 10 everywhere (though is slightly higher close to the star where temperatures are hotter); the reduction in mass-loss caused by this process is not isolated to any one location in the wind and is in line with our statement that the reduced cooling allows denser gas to reach the escape temperature. On the other hand, when cooling due to H2 is disabled, the mass-loss rates are fairly similar around R ≲ 3 au and at R ≳ 50 au. Note that smaller radii the mass-loss rates are still much smaller than previous works because the density of the wind has to be lower to overcome the additional heating from collisional excitation of O by neutral H and still reach the ∼ 10 4 K escape temperature. The greatest difference is seen around R = 10au where H2 starts to survive into the wind and where the wind is warm enough to excite the H2 lines leading to additional cooling. At smaller radii there is no H2 in the wind to provide cooling, while outside 50 au the wind is too cold to excite it.

## 6. Discussion

## 6.1. Necessity of this method

There remains the question of whether we needed to adopt such as expensive alternative method of calculating the thermochemistry. It has previously been argued that a ξ -T approach neglects the adiabatic cooling due to expanding gas and thus risks double counting the input of energy (Wang &amp; Goodman 2017). However, we find that this term is nearly everywhere subdominant to the atomic and or molecular line cooling and radiative ther-/ mal equilibrium is a good approximation . This is supported by 8 the close correlation between the ionisation parameter ξ and gas temperature T (as shown in Sect. 5.1.4), which demonstrates that this can remain an appropriate parametrization of the underlying microphysics.

However, we note that over much of the wind, we found cooling to be dominated by H2, the abundance of which is out of chemical equilibrium due to advection. Komaki et al. (2021), who derive similar photoevaporation rates also see partially molecular winds, and the importance of advection for the survival of molecules into the wind has also previously been argued for photoevaporation (Wang &amp; Goodman 2017; Nakatani et al. 2018a,b) and MHD winds (Panoglou et al. 2012; Wang et al. 2019). This is the most important impact of the hydrodynamics on the temperature structure in our work: such an e ff ect could not be captured by the sort of static slab models used to calibrate the ξ -T relationships before, thus suggesting it is necessary to perform the chemistry on the fly. However, since it does seem that despite this, a close ξ -T relationship can be found, and since it has previously been demonstrated that the ξ -T approach can be adapted to varying chemical abundances (Wölfer et al. 2019) - as well as column densities (Picogna et al. 2019) and di ff erent spectra (Ercolano et al. 2021) - a future intermediate possibility could be to produce a set of ξ -T relationships at di ff ering H2 abundances to simplify at least the temperature calculation.

## 6.2. Dust entrainment

Dust entrainment in winds is an important question both as a sink of mass, and as a potentially observable tracer of the wind (Owen et al. 2011b; Franz et al. 2022a,b). Both the delivery of dust to the base of the wind (Hutchison &amp; Clarke 2021; Booth &amp; Clarke 2021), and the ultimate trajectory in the wind (e.g. Franz et al. 2020) can play a role in this picture. Moreover, the level of dust may a ff ect the thermochemistry in the wind as discussed in Sect. 4.2; we therefore consider how self-consistent our simulations are.

Using our fitted mass-loss profile to the fiducial model in Eq. 25 of Booth &amp; Clarke (2021), we show the maximum entrained grain size as a function of radius in Fig. 14. The values lie in the range of 0 04 . -0 4 . µ m and are typically 1-2 orders of magnitude below those that Booth &amp; Clarke (2021) calculated using a previous X-ray photoevaporation model (Picogna et al. 2019), in line with our mass-loss rate which is approximately an order of magnitude smaller. This range straddles our assumed maximum grain size of 0 25 . µ m, suggesting that this is a reasonable value, though we also note that the entrained size is smaller in the inner disc compared to the outer.

8 The e ff ects of adiabatic cooling can be seen in Fig. 10 as the small number of points with ξ &gt; 10 -3 and T &lt; 10 4 Kand correspond to bound cells at small radii near the rotation axis.

Fig. 13: The mass-loss rate per unit area as a function of distance to the star. Left: di ff erent times (from light to dark blue with increasing time) for our fiducial model compared to those from Picogna et al. (2021) for the same stellar mass and X-ray spectrum (black). Right: comparing the fiducial model (dark blue) with those where we switch o ff certain areas of microphysics in order to test their impact.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-6ad5cb30-ef35-4efe-adce-4094c6e951ab.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-74de586d-0480-4ab0-b389-e24a33d79f2a.jpg)

By assuming that the maximum dust size in the disc is set by fragmentation according to Eq. 3 of Birnstiel et al. (2012) with α = 10 -3 and u frag = 10 2 cms -1 , and assuming a Lynden-Bell &amp; Pringle (1974) disc profile with M D = 0 1 M . ⊙ , R C = 100 au and T 1au = 270 K, we can estimate the maximum grain size at each radius in the disc a frag . For a Mathis et al. (1977) grain size that when dust is included, high photoevaporation rates for all stars would not allow discs to reach such low dust masses as are observed while remaining accreting.

distribution, then the entrained fraction is f ent = GLYPH&lt;18&gt; a max a frag GLYPH&lt;19&gt; 1 2 / . For the Picogna et al. (2019) mass-loss profile, this estimate suggests up to ∼ 10 per cent entrained mass (or a dust-to-gas ratio of 10 -3 ), consistent with the modelling by Franz et al. (2022a). Meanwhile, our mass-loss profile suggests values up to 5 times lower, largely in the range ϵ = 10 -4 -10 -3 . This suggests that the dust-to-gas ratio we assume is 10 -100 times too high, and we should explore lower values in future models.

## 6.3. Consequences for disc evolution &amp; population synthesis

One of the major usages of photoevaporation prescriptions such as we derive in Sect. 5.2.3 is to include in disc evolution models, often for the purposes of population synthesis. One limitation of the previous X-ray rates is that they needed fairly careful selection of the parameter space in order to get long-enough lived discs. For example, several works find it di ffi cult to get discs that live for 3 Myr or more (and therefore survive to the ages of older regions such as Upper Sco) when using the median X-ray luminosities (Emsenhuber et al. 2023; Appelgren et al. 2023). Similarly, Somigliana et al. (2020) predict a 'knee' in the plane of accretion rate vs disc mass that is not observed, which limits photoevaporation to low values, and Sellek et al. (2020) argue

One way to understand this is with the so-called 'UV-switch' model (Clarke et al. 2001): once the photoevaporation rate exceeds the accretion rate, accretion can no longer resupply the inner disc as all of the material it brings in is immediately ejected in the wind. This results in the opening of a gas cavity. When the accretion rates are higher than 10 -8 M yr ⊙ -1 , they are comparable to the median accretion rate (Manara et al. 2023) and consequently a disc may have very little chance to evolve before reaching this criterion. The lower values that we find thus delay the onset of the rapid disc dispersal.

This has motivated people to question which rates best fit the distribution of accretion rates. On the one hand, Alexander et al. (2023) showed that for 0 3 . ≤ M ∗ / M ⊙ ≤ 1 2 if . the mass-loss rate exceeds ∼ 10 -9 M yr ⊙ -1 then the process of gap opening, inner disc dispersal and cessation of accretion happens so quickly that there would be too few discs with low accretion rates to match the observed distribution. On the other hand, Ercolano et al. (2023) argue that if one constructs a population of discs to fit the disc fraction with age (i.e. the lifetime distribution) under the scenarios of EUV-, X-ray- and FUV- driven photoevaporation, the best fit to the distribution of low accretion rates is with X-ray photoevaporation (for which they used the higher rates of Picogna et al. 2021). This happens because EUV photoevaporation requires more tenuous discs in order for them to be dispersed in the same time, and therefore lower accretion rates.

Nevertheless, our lower mass-loss rates are a step in the right direction in terms of allowing photoevaporation to be consistent

Fig. 14: Left: the maximum entrained grain size in the wind according to the prescription of Booth &amp; Clarke (2021) for our fiducial model (blue), compared to that expected previous model of Picogna et al. (2019). Right: an estimate for the resulting fraction of dust mass dust-to-gas ratio in the wind. The grey lines indicate the values used in this model. /

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-eb3d372d-8ae8-48b1-890a-9b537be9fd07.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-5c8ed7fc-a6ff-4362-aac5-e43164b2422e.jpg)

with the distributions of disc accretion rates and lifetimes found most commonly in the literature. We illustrate this in Fig. 15, where we compare the evolution of identical 100 MJ, 100 au discs around solar-mass stars with undergoing viscous evolution ( α = 10 -3 ) and photoevaporation according to the prescriptions of Picogna et al. (2021) and Sect. 5.2.3 of this work. Since less mass is removed, our prescription results in larger gas masses at all times, and allows the disc to keep accreting at detectable levels for the full 10 Myr that we simulated. Moreover, around 10 per cent of the mass is left by the time of the cavity opening, suggesting that thermal sweeping would not kick in to prevent the existence of massive but non-accreting 'relic discs'.

The lower densities - and shallower radial run of density may also have important consequences. Sellek et al. (2024) showed that for appropriate column densities, the absorption of soft Xrays in the inner wind can lead to [Ar ii ] emission becoming more compact than the [Ne ii ] emission. In our models, more absorption would occur at larger radii, making such a configuration harder to achieve. Moreover, having shallower density profiles will produce narrower line profiles (Ballabio et al. 2020).

A final consequence of lowering the internal photoevaporation rate is that the relative importance of other mass-loss processes, such as external photoevaporation may be increased in a given environment. This may, for example, cap the possible disc lifetimes and lead to fewer cases of inside-out disc clearing (Coleman &amp; Haworth 2022).

## 6.4. Observables

The line emission from winds is collisionally excited and, consequently, sensitive to temperature. Depending on the excitation temperature, the reduction in the gas temperature at large radii in the winds is therefore also likely to reduce the population of the upper levels of the lines. For example, the commonly studied [Ne ii ] 12 81 . µ mline has an excitation temperature of 1123 K and therefore would not be emitted beyond about 60 au. However, in previous models, this line mainly originated inside 10 -20 au anyway (Ercolano &amp; Owen 2010; Ercolano et al. 2021), suggesting that any e ff ect would be small. However, given that the velocities in the wind are of order c S, the cooler wind will also be slower, and thus present high-spectral-resolution line profiles that are less blue-shifted and narrower (Ballabio et al. 2020).

Finally, the fact that molecules - particularly H2 - survive further into the wind (and closer to the star) than otherwise expected due to advection suggests that photoevaporation could also be traced by molecular lines, such as the 2 12 . µ m line (Gangi et al. 2020). This was already shown in Rab et al. (2022) for photoevaporative winds using a simpler post-processing approach for the thermochemistry. However, as noted in Rab et al. (2022), and as shown in this work, such an approach has to be used with care. Since the material is seen further along the streamlines, it will have higher velocities than static post-processing would suggest and, therefore, produce broader and more highly blue-shifted lines. We intend to use our model to produce synthetic observables in a future study in order to quantify more precisely the importance of advection for interpreting molecular wind tracers.

## 7. Conclusions

In this work we have presented improvements to the thermochemistry code prizmo (Grassi et al. 2020) - designed for use onthe-fly with hydrodynamics - and benchmarked its performance for describing the temperatures of photoevaporative winds in different regimes against previous works and other codes. We have then demonstrated the first science application to the problem of the thus-far divergent mass-loss rates for photoevaporation in the literature, focusing on the X-ray regime. By a) designing our ap-

## Picogna et al. 2021 This work

/circledot

Fig. 15: Disc gas mass (left) and accretion rate (right) over the course of 10 Myr for an example 100 MJ, 100 au disc around a solar-mass star undergoing viscous evolution ( α = 10 -3 ) and photoevaporation using the mass-loss prescription of Picogna et al. (2021) (black) or from this work (blue).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-30edc20a-55d6-4a4e-aa63-6bc92d5ed9f5.jpg)

proach to be su ffi ciently flexible and comprehensive to capture thermochemistry in ionised, neutral atomic and molecular environments at a range of temperatures and b) testing spectra used in di ff erent previous works in a controlled manner, we have ensured that we can confidently isolate the root of the di ff erences we find and probe the causes of the di ff erent claims in the literature. In particular, we find that dense gas; this lowers the mass-loss rates of X-ray photoevaporation for a solar mass star by over an order of magnitude to ∼ 10 -9 M yr ⊙ -1 (Fig. 11a). The mass-loss is reduced everywhere by a similar amount compared to previous works (Fig. 13).

- 1. The temperatures of X-ray driven photoevaporative winds, which are very weakly ionised, can be greatly overestimated if neutral H is not included as a collider to excite atomic forbidden line emission (Fig. 3).
- 2. The near-isothermal temperatures seen in previous X-ray photoevaporation models can be understood as a result of the very steep relationships between the temperature and ionisation parameter ξ used therein leading to a dearth of cells at unstable temperatures ∼ 1000 K (Fig. 10). In turn this occurs because of the bimodal distribution of atomic emission lines between the optical and far-IR regimes; molecular cooling in the mid-IR stabilises the intermediate temperatures leading to shallower relationship. As a result, the new simulations self-regulate to achieve the 1 / r escape temperature predicted by Owen et al. (2012) (Fig. 6).
- 3. Adiabatic cooling of the gas is only dominant in a very narrow region along the H2 dissociation front (Fig. 6), and very low-density gas along the rotation axis. However, it does become relatively more important at larger radii where the gas is cooler as the excitation of line cooling becomes harder.
- 4. As a result of the additional cooling mechanisms discussed above, the escape temperature can only be reached in less-
- 5. The winds are driven by soft X-rays; use of an observationally-derived spectrum without correcting for the foreground attenuation of these energies can lead to underestimates of the X-ray contribution to photoevaporation (Fig. 11a).
- 6. The winds contain a significant amount of molecular hydrogen; advection of material into the wind -i.e. nonequilibrium chemistry - is critical for its existence (Fig. 8). This will be important for matching the extent and velocity of line diagnostics probing molecular winds.

The lower mass-loss rates for X-ray photoevaporation, which we have understood to result from the additional cooling provided, in particular, by excitation of O by neutral H, help to resolve long-standing tensions in the field of photoevaporation simulations and move us towards a consensus on the expected photoevaporation rates.

These results have significant implications for the evolution and demographics of protoplanetary discs. As a result of the lower mass-loss rates, the lifetime of discs can be considerably longer than expected previously and they should be able to reach lower accretion rates before undergoing rapid cavity opening and dispersal. This helps resolve tensions with works on disc demographics which generally imply lower photoevaporation rates than the previous X-ray models suggested. Future work will need

to explore the parameter space, to see how these rates depend on the amount of FUV irradiation as well as the stellar mass, which will enable the production of a general prescription to be used in disc and planet population synthesis models.

## Data availability

The files needed to setup and run the fiducial model, as well as its final output, are available from Zenodo: https://doi.org/10.5281/zenodo.10891357 . The version of pluto with the prizmo interface is available from GitHub: https://github.com/AndrewSellek/PLUTO\_ PRIZMO/tree/chemistry , while prizmo itself is available from GitHub: https://github.com/tgrassi/prizmo .

Acknowledgements. We thank the reviewer, James Owen, for his detailed feedback which helped us to clarify various aspects of the microphysics and technical approach. We are grateful to Lile Wang for sharing the density and temperature grids from Wang &amp; Goodman (2017) used in the benchmarking, Catherine Espaillat for providing the DIAD models and Riouhei Nakatani for useful discussions about the comparison between our model and his and for sharing details of his X-ray spectrum. ADS also thanks James Matthews for useful conversations about the similarities to thermal winds from X-ray binaries. ADS thanks the Science and Technology Facilities Council (STFC) for a Ph.D. studentship for project 2277492 as part of Training Grant ST T505985 1 and was also supported / / by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 1010197S1 MOLDISK). CJC acknowledges support from the Science &amp; Technology Facilities Council (STFC) Consolidated Grant ST W000997 1 and Eu-/ / ropean Union's Horizon 2020 research and innovation programme No 823823 (DUSTBUSTERS) under the Marie Sklodowska-Curie grant agreement. We acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Research Unit 'Transition discs' - 325594231. This research was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311. CHR is grateful for support from the Max Planck Society. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP T022159 1), and DiRAC funding from the Science and Technology Facili-/ / ties Council (www.dirac.ac.uk). This work also benefited from the Core2disk-III residential program of Institut Pascal at Université Paris-Saclay, with the support of the program 'Investissements d'avenir' ANR-11-IDEX-0003-01.

## References

Abrahamsson, E., Krems, R. V., &amp; Dalgarno, A. 2007, ApJ, 654, 1171

- Alexander, R., Pascucci, I., Andrews, S., Armitage, P., &amp; Cieza, L. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, &amp; T. Henning, 475

Alexander, R., Rosotti, G., Armitage, P. J., et al. 2023, MNRAS, 524, 3948

Appelgren, J., Lambrechts, M., &amp; van der Marel, N. 2023, A&amp;A, 673, A139

Arulanantham, N., McClure, M. K., Pontoppidan, K., et al. 2024, ApJ, 965, L13

Bai, X.-N. 2017, ApJ, 845, 75

- Ballabio, G., Alexander, R. D., &amp; Clarke, C. J. 2020, MNRAS, 496, 2932
- Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76
- Barinovs, G., van Hemert, M. C., Krems, R., &amp; Dalgarno, A. 2005, ApJ, 620, ˘ 537
- Bell, K. L., Berrington, K. A., &amp; Thomas, M. R. J. 1998, MNRAS, 293, L83
- Birnstiel, T., Klahr, H., &amp; Ercolano, B. 2012, A&amp;A, 539, A148
- Blandford, R. D. &amp; Payne, D. G. 1982, MNRAS, 199, 883
- Booth, R. A. &amp; Clarke, C. J. 2021, MNRAS, 502, 1569
- Chiang, E. I. &amp; Goldreich, P. 1997, ApJ, 490, 368
- Clarke, C. J. &amp; Alexander, R. D. 2016, MNRAS, 460, 3044
- Clarke, C. J., Gendrin, A., &amp; Sotomayor, M. 2001, MNRAS, 328, 485
- Coleman, G. A. L. &amp; Haworth, T. J. 2022, MNRAS, 514, 2315
- Currie, T., Lawson, K., Schneider, G., et al. 2022, Nature Astronomy, 6, 751 D'Alessio, P., Calvet, N., &amp; Hartmann, L. 2001, ApJ, 553, 321
- D'Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., &amp; Servín, H. 2006, ApJ, 638, 314
- D'Alessio, P., Calvet, N., Hartmann, L., Lizano, S., &amp; Cantó, J. 1999, ApJ, 527, 893
- D'Alessio, P., Cantö, J., Calvet, N., &amp; Lizano, S. 1998, ApJ, 500, 411
- Draine, B. T. 2003a, ApJ, 598, 1017
- Draine, B. T. 2003b, ApJ, 598, 1026
- Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
- Draine, B. T. &amp; Bertoldi, F. 1996, ApJ, 468, 269
- Dullemond, C. P., Hollenbach, D., Kamp, I., &amp; D'Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, &amp; K. Keil, 555
- Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multipurpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015
- Emsenhuber, A., Burn, R., Weder, J., et al. 2023, A&amp;A, 673, A78
- Ercolano, B., Barlow, M. J., &amp; Storey, P. J. 2005, MNRAS, 362, 1038
- Ercolano, B., Barlow, M. J., Storey, P. J., &amp; Liu, X. W. 2003, MNRAS, 340, 1136
- Ercolano, B., Clarke, C. J., &amp; Drake, J. J. 2009, ApJ, 699, 1639
- Ercolano, B. &amp; Owen, J. E. 2010, MNRAS, 406, 1553
- Ercolano, B. &amp; Owen, J. E. 2016, MNRAS, 460, 3472
- Ercolano, B., Picogna, G., &amp; Monsch, K. 2023, MNRAS, 526, L105
- Ercolano, B., Picogna, G., Monsch, K., Drake, J. J., &amp; Preibisch, T. 2021, MNRAS, 508, 1675
- Ercolano, B., Young, P. R., Drake, J. J., &amp; Raymond, J. C. 2008, ApJS, 175, 534 Flower, D. R. &amp; Launay, J. M. 1977, Journal of Physics B Atomic Molecular Physics, 10, 3673
- Franz, R., Ercolano, B., Casassus, S., et al. 2022a, A&amp;A, 657, A69
- Franz, R., Picogna, G., Ercolano, B., &amp; Birnstiel, T. 2020, A&amp;A, 635, A53
- Franz, R., Picogna, G., Ercolano, B., et al. 2022b, A&amp;A, 659, A90
- Fuente, A., Rivière-Marichalar, P., Beitia-Antero, L., et al. 2023, A&amp;A, 670, A114
- Gangi, M., Nisini, B., Antoniucci, S., et al. 2020, A&amp;A, 643, A32

Glassgold, A. E., Najita, J. R., &amp; Igea, J. 2007, ApJ, 656, 515

- Glover, S. C. O. 2015, MNRAS, 451, 2082
- Glover, S. C. O. &amp; Abel, T. 2008, MNRAS, 388, 1627
- Glover, S. C. O. &amp; Jappsen, A. K. 2007, ApJ, 666, 1
- Gorti, U. &amp; Hollenbach, D. 2009, ApJ, 690, 1539
- Gorti, U., Hollenbach, D., &amp; Dullemond, C. P. 2015, ApJ, 804, 29
- Grassi, T., Bovino, S., Haugbølle, T., &amp; Schleicher, D. R. G. 2017, MNRAS, 466, 1259
- Grassi, T., Ercolano, B., Sz˝cs, L., Jennings, J., &amp; Picogna, G. 2020, MNRAS, u 494, 4471
- Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&amp;A, 623, A140
- Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126

Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51

- Haworth, T. J., Clarke, C. J., &amp; Owen, J. E. 2016, MNRAS, 457, 1905
- Higginbottom, N. &amp; Proga, D. 2015, ApJ, 807, 107
- Higginbottom, N., Proga, D., Knigge, C., &amp; Long, K. S. 2017, ApJ, 836, 42
- Hily-Blant, P., Pineau des Forêts, G., Faure, A., &amp; Lique, F. 2022, A&amp;A, 658, A168
- Hollenbach, D., Johnstone, D., Lizano, S., &amp; Shu, F. 1994, ApJ, 428, 654 Hutchison, M. A. &amp; Clarke, C. J. 2021, MNRAS, 501, 1127
- Hutchison, M. A., Price, D. J., Laibe, G., &amp; Maddison, S. T. 2016, MNRAS, 461, 742
- Kallemeyn, G. W., Rubin, A. E., Wang, D., &amp; Wasson, J. T. 1989, Geochim. Cosmochim. Acta, 53, 2747
- Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114
- Kamp, I. &amp; Bertoldi, F. 2000, A&amp;A, 353, 276
- Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&amp;A, 607, A41
- Kamp, I., Tilling, I., Woitke, P., Thi, W. F., &amp; Hogerheijde, M. 2010, A&amp;A, 510, A18
- Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., &amp; Weintraub, D. A. 2002, ApJ, 567, 434
- Keppler, M., Benisty, M., Müller, A., et al. 2018, A&amp;A, 617, A44
- Koepferl, C. M., Ercolano, B., Dale, J., et al. 2013, MNRAS, 428, 3327
- Komaki, A., Nakatani, R., &amp; Yoshida, N. 2021, ApJ, 910, 51
- Kramida, A., Ralchenko, Y., Reader, J., &amp; NIST ASD Team. 2022, NIST Atomic Spectra Database, accessed Sun 29th Jan 2023
- Krems, R. V., Jamieson, M. J., &amp; Dalgarno, A. 2006, ApJ, 647, 1531
- Kuiper, R., Klahr, H., Dullemond, C., Kley, W., &amp; Henning, T. 2010, A&amp;A, 511, A81
- Kunitomo, M., Suzuki, T. K., &amp; Inutsuka, S.-i. 2020, MNRAS, 492, 3849
- Levermore, C. D. &amp; Pomraning, G. C. 1981, ApJ, 248, 321 Li ff man, K. 2003, PASA, 20, 337
- Lique, F., Kłos, J., Alexander, M. H., Le Picard, S. D., &amp; Dagdigian, P. J. 2018, MNRAS, 474, 2313
- Lynden-Bell, D. &amp; Pringle, J. E. 1974, MNRAS, 168, 603
- Maloney, P. R., Hollenbach, D. J., &amp; Tielens, A. G. G. M. 1996, ApJ, 466, 561 Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed.
- S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, &amp; M. Tamura, 539 Mathis, J. S., Rumpl, W., &amp; Nordsieck, K. H. 1977, ApJ, 217, 425
- Michel, A., van der Marel, N., &amp; Matthews, B. C. 2021, ApJ, 921, 72

- Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
- Müller, A., Keppler, M., Henning, T., et al. 2018, A&amp;A, 617, L2
- Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., &amp; Kuiper, R. 2018a, ApJ, 857, 57
- Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., &amp; Kuiper, R. 2018b, ApJ, 865, 75
- Nomura, H., Aikawa, Y., Tsujimoto, M., Nakagawa, Y., &amp; Millar, T. J. 2007, ApJ, 661, 334
- Omukai, K., Hosokawa, T., &amp; Yoshida, N. 2010, ApJ, 722, 1793
- Osterbrock, D. E. &amp; Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (University Science Books)
- Owen, J. E., Clarke, C. J., &amp; Ercolano, B. 2012, MNRAS, 422, 1880
- Owen, J. E., Ercolano, B., &amp; Clarke, C. J. 2011a, MNRAS, 412, 13
- Owen, J. E., Ercolano, B., &amp; Clarke, C. J. 2011b, MNRAS, 411, 1104
- Owen, J. E., Ercolano, B., Clarke, C. J., &amp; Alexander, R. D. 2010, MNRAS, 401, 1415
- Owen, J. E., Hudoba de Badyn, M., Clarke, C. J., &amp; Robins, L. 2013, MNRAS, 436, 1430
- Panoglou, D., Cabrit, S., Pineau Des Forêts, G., et al. 2012, A&amp;A, 538, A2 Parker, E. N. 1958, ApJ, 128, 664
- Pascucci, I., Banzatti, A., Gorti, U., et al. 2020, ApJ, 903, 78
- Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, &amp; M. Tamura, 567
- Pequignot, D. 1990, A&amp;A, 231, 499
- Pfalzner, S., Dehghani, S., &amp; Michel, A. 2022, ApJ, 939, L10
- Picogna, G., Ercolano, B., &amp; Espaillat, C. C. 2021, MNRAS, 508, 3611
- Picogna, G., Ercolano, B., Owen, J. E., &amp; Weber, M. L. 2019, MNRAS, 487, 691
- Pinte, C., Ménard, F., Duchêne, G., &amp; Bastien, P. 2006, A&amp;A, 459, 797
- Rab, C., Güdel, M., Woitke, P., et al. 2018, A&amp;A, 609, A91
- Rab, C., Weber, M., Grassi, T., et al. 2022, A&amp;A, 668, A154
- Ribas, Á., Bouy, H., &amp; Merín, B. 2015, A&amp;A, 576, A52
- Ribas, Á., Merín, B., Bouy, H., &amp; Maud, L. T. 2014, A&amp;A, 561, A54
- Richings, A. J., Schaye, J., &amp; Oppenheimer, B. D. 2014, MNRAS, 442, 2780
- Savage, B. D. &amp; Sembach, K. R. 1996, ARA&amp;A, 34, 279
- Schroder, K., Staemmler, V., Smith, M. D., Flower, D. R., &amp; Jaquet, R. 1991, Journal of Physics B Atomic Molecular Physics, 24, 2487
- Sellek, A. D., Bajaj, N. S., Pascucci, I., et al. 2024, AJ, 167, 223
- Sellek, A. D., Booth, R. A., &amp; Clarke, C. J. 2020, MNRAS, 498, 2845
- Sellek, A. D., Clarke, C. J., &amp; Ercolano, B. 2022, MNRAS, 514, 535
- Simon, M. &amp; Prato, L. 1995, ApJ, 450, 824
- Skrutskie, M. F., Dutkevitch, D., Strom, S. E., et al. 1990, AJ, 99, 1187
- Somigliana, A., Toci, C., Lodato, G., Rosotti, G., &amp; Manara, C. F. 2020, MNRAS, 492, 1120
- Spitzer, L. 1978, Physical processes in the interstellar medium
- Tarter, C. B., Tucker, W. H., &amp; Salpeter, E. E. 1969, ApJ, 156, 943
- Tayal, S. S. 2007, ApJS, 171, 331
- Tayal, S. S. 2008, A&amp;A, 486, 629
- Teyssier, R. 2015, Annual Review of Astronomy and Astrophysics, 53, 325
- Thi, W. F., Woitke, P., &amp; Kamp, I. 2011, MNRAS, 412, 711
- van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., &amp; van Dishoeck, E. F. 2020, Atoms, 8, 15
- Verner, D. A., Ferland, G. J., Korista, K. T., &amp; Yakovlev, D. G. 1996, ApJ, 465, 487
- Verner, D. A. &amp; Yakovlev, D. G. 1995, A&amp;AS, 109, 125
- Wang, L., Bai, X.-N., &amp; Goodman, J. 2019, ApJ, 874, 90
- Wang, L. &amp; Goodman, J. 2017, ApJ, 847, 11
- Weder, J., Mordasini, C., &amp; Emsenhuber, A. 2023, A&amp;A, 674, A165
- Weingartner, J. C. &amp; Draine, B. T. 2001, ApJS, 134, 263
- Weingartner, J. C., Draine, B. T., &amp; Barr, D. K. 2006, ApJ, 645, 1188
- Woitke, P., Kamp, I., &amp; Thi, W. F. 2009, A&amp;A, 501, 383
- Woitke, P., Min, M., Pinte, C., et al. 2016, A&amp;A, 586, A103
- Wölfer, L., Picogna, G., Ercolano, B., &amp; van Dishoeck, E. F. 2019, MNRAS, 490, 5596

Wolfire, M. G., Vallini, L., &amp; Chevance, M. 2022, ARA&amp;A, 60, 247

- Wolk, S. J. &amp; Walter, F. M. 1996, AJ, 111, 2066
- Zhukovska, S., Gail, H. P., &amp; Trielo ff , M. 2008, A&amp;A, 479, 453

## Appendix A: Fits to OI collisional de-excitation data

We perform our own fits to the tabulation of the de-excitation rates of O i by H and H2 from Lique et al. (2018), in the form

$$k _ { \text{ul} } = ( a \times T _ { 2 } ^ { b + c \ln T _ { 2 } } + d \exp ( - e / T _ { 2 } ) ) \times 1 0 ^ { - 1 0 } \, \text{cm} ^ { 3 } \, s ^ { - 1 }, \quad ( \text{A.1} )$$

where T 2 = T 10 2 K . The values of the parameters are reported in Table A.1.

## Appendix B: Timing of prizmo

In Fig. B.1 we show information regarding the time taken by prizmo to calculate the thermochemistry in each cell in the simulation after 50 orbits at 10 au (the last snapshot before we freeze the simulation inside 17 au).

Generally, the most computationally expensive cells are found at heights corresponding to the disc - rather than the wind. However, the cells near the midplane are not the slowest since here we fix the temperature to the dust temperature from the DIAD models and hence the chemistry and temperature determination are decoupled. Moreover, the longest times are seen for cells inside 10 au (as shown in the inset plot), illustrating the benefit of freezing the simulation at smaller radii.

Nevertheless, the histogram of the times shows that the distribution of the time per cell is fairly smooth and the overall time is not dominated by expensive outliers. Indeed since we plot d n d log t , the y-axis represents the contribution of each timescale to the overall execution time. This suggests that the large number of cells with a prizmo execution time t prizmo ∼ 10 -3 s dominate the total time rather than the relatively smaller number of more expensive cells. Thus, there is no easy way to decrease the total computational timescale without reducing the number of cells considerably.

Of course, since the computation is parallelised, this does not ensure that the processors are load balanced. pluto allows the user to choose the arrangement of the processors but not to vary the number of cells between each. Therefore, if a single processor were to contain all of the most expensive cells, these would limit the computational timescale as the other processors may have to wait a long time for this single processor to finish its calculation. To minimise the risk of this, since more expensive cells are typically found in lower latitudes or along the disc-wind interface, we try to minimise the number of cells at a particular latitude by arranging the processors to have many in the r direction since that their radial extent is large and all spend some time working on the underlying disc. While refinements may be possible, we find this gives satisfactorily little variation between the total execution time for the chemistry per processor.

Table A.1: Parameters derived for the fits to collisional de-excitation rates of O i by H and H2

| Collider   | Transition    | Energy Levels   |   Wavelength ( µ m) |       a |     b |      c |      d |       e |
|------------|---------------|-----------------|---------------------|---------|-------|--------|--------|---------|
| H          | 3 P 0 → 3 P 1 | E 2 → E 1       |                 145 | 1.19    | 0.16  |  0.072 |  5.03  | 18.9    |
|            | 3 P 1 → 3 P 2 | E 1 → E 0       |                  63 | 1.37    | 0.348 |  0.018 |  2.96  |  3.77   |
|            | 3 P 0 → 3 P 2 | E 2 → E 0       |                  44 | 1.03    | 0.304 |  0.033 |  3.08  |  5.18   |
| para-H 2   | 3 P 0 → 3 P 1 | E 2 → E 1       |                 145 | 0.00719 | 0.624 |  0.051 |  0.338 | 36.2    |
|            | 3 P 1 → 3 P 2 | E 1 → E 0       |                  63 | 1.7     | 0.221 |  0.022 | -0.429 |  1.63   |
|            | 3 P 0 → 3 P 2 | E 2 → E 0       |                  44 | 0.00446 | 2.596 | -0.235 |  1.54  |  0.0757 |
| ortho-H 2  | 3 P 0 → 3 P 1 | E 2 → E 1       |                 145 | 0.0472  | 0.393 |  0.038 |  0.23  | 29.3    |
|            | 3 P 1 → 3 P 2 | E 1 → E 0       |                  63 | 1.84    | 0.215 |  0.022 |  1.54  |  8.13   |
|            | 3 P 0 → 3 P 2 | E 2 → E 0       |                  44 | 1.23    | 0.355 |  0.012 |  1.88  | 19.1    |

Notes. The form of the fits is given by Eq. A.1 and the data were taken from from Lique et al. (2018).

Fig. B.1: The distribution of the execution time per cell of prizmo after 50 orbits at 10 au.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-d2f7632f-41f4-485a-b642-941497876942.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00848v1-figure-55f47188-8506-421d-930e-1694560bc076.jpg)
Andrew D. Sellek Tommaso Grassi Giovanni Picogna Christian Rab Cathie J. Clarke Barbara Ercolano
Photoevaporation of protoplanetary discs with PLUTO+PRIZMO I. Lower X-ray-driven mass-loss rates due to enhanced cooling
Context: Photoevaporation is an important process for protoplanetary disc
dispersal but there has so far been a lack of consensus from simulations over
the mass-loss rates and the most important part of the high-energy spectrum for
driving the wind. Aims: We aim to isolate the origins of these discrepancies
through carefully-benchmarked hydrodynamic simulations of X-ray
photoevaporation with time-dependent thermochemistry calculated on the fly.
Methods: We conduct hydrodynamic simulations with pluto where the
thermochemistry is calculated using prizmo. We explore the contribution of
certain key microphysical processes and the impact of using different spectra
used previously in literature studies. Results: We find that additional cooling
results from the excitation of O by neutral H, which leads to dramatically
reduced mass-loss across the disc compared to previous X-ray photoevaporation
models, with an integrated rate of 10^-9 Msun/yr. Such rates would allow for
longer-lived discs than previously expected from population synthesis. An
alternative spectrum with less soft X-ray produces mass-loss rates around a
factor of 2-3 times lower. The chemistry is significantly out of equilibrium,
with the survival of H2 into the wind aided by advection. This leads to its
role as the dominant coolant at 10s au - thus stabilising a larger radial
temperature gradient across the wind - as well as providing a possible wind
tracer.