© The Authors 2024

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, which is thought to happen rapidly (Skrutskie et al. 1990; Simon & Prato 1995; Wolk & 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).

Disc winds are widely documented throughout the lifetime of a disc (Pascucci et al. 2023) and are one important factor contributing to the evolution and dispersal of 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 effects 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 accretion become less significant as the disc accretion rate drops. Therefore, the 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.

Photoevaporation also explains a number of observational constraints that 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, within the inner few au of the disc. Such a mode of dispersal is supported by the relative disc lifetimes at different wavelengths (Ribas et al. 2014) and the evolutionary loci of discs in colour-colour space (Koepferl et al. 2013), that 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⁻⁹ M_⊙ yr⁻¹ for models driven by extreme ultraviolet (EUV) (Hollenbach et al. 1994) to 10⁻⁸−10⁻⁷ M_⊙ yr⁻¹ 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 differences may be attributed to different methodologies since the problem of calculating thermochemistry on the fly with hydrodynamics is very expensive and thus different works have compromised on different factors. For example, Hollenbach et al. (1994) and Gorti & 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 $\xi =\frac{L_X}{n{r}^2}$ Tarter et al. 1969), that 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 effects 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 and X-ray radiation. In recent years, there have been moves towards on-the-fly thermochemistry and hydrodynamics both for photoevaporation (e.g. Wang & Goodman 2017; Nakatani et al. 2018a,b; Komaki et al. 2021) and for magnetohydodynamic (MHD) winds (Wang et al. 2019; Gressel et al. 2020), but such efforts naturally come at the expense of solving the thermochemistry in detail and hence are very sensitive to the simplifying assumptions that are made.

The fact that several things – such as the irradiating spectrum, possible coolants, radiative transfer method, grid extent, explicit calculation of hydrodynamics, and assumption of thermochemical equilibrium – typically differ between each methodology has made it very difficult to make direct comparisons between the above works and clearly assess the dominant factors responsible for the divergent results. For example, Wang & Goodman (2017) – who favoured a low mass-loss rate EUV-driven wind – attributed the higher mass-loss rates of X-ray models to their lack of molecular cooling. On the other hand, Sellek et al. (2022) demonstrate that a lack of effective optical line cooling in the work by Wang & Goodman (2017) meant that their predicted wind temperatures are too high, so despite the additional molecular cooling, this thermochemical treatment was still not sufficiently complete. More importantly, Sellek et al. (2022) show that the results were very sensitive to the shape of the X-ray spectrum, finding that the simplified spectrum used by Wang & Goodman (2017), where they assumed all the X-rays 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 & 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 sufficiently 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.

In this work, Section 2 sets out the key features of PRIZMO and any choices we made for the thermochemistry, while Sec. 3 explains how we used PRIZMO in conjunction with PLUTO for the hydrodynamics. Section 4 presents benchmarks of the temperature structure of winds in different regimes to illustrate the key successes of PRIZMO and elucidates important differences 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 Sec. 5. We also derive the corresponding mass-loss profile, establish the importance of non-equilibrium thermochemistry, and explore how the mass-loss rate is affected 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 Sec. 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 sufficient 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 Sec. 2.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 (Sec. 3).

The key inputs to the preprocessor are:

• Atomic data (level energies, degeneracies, Einstein coefficients, and fits to collisional de-excitation rates).

• A list of (photo)chemical reactions with rate coefficient expressions.

• A user-defined spectrum and number of energy bins.

• Frequency-dependent dust optical constants.

Our preprocessor inputs are detailed in Sects. 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).

2.1 Updates to PRIZMO

2.1.1 Structure and 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 affording 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 coefficient calculations in the chemical network by grouping chemical rate coefficients by the same arithmetic expression, for example α(T/300 K)^β exp(−γ/T) and α(T/300 K)^β belong to different groups. Since most of the rates have relatively simple algebraic forms, with this approach we obtain a speed-up similar to interpolating the rate coefficients in temperature in the log-log space.

We have also approximated some rates in order to pre-compute common algebraic quantities; for example, (T/300 K)^−0.49 has been approximated to (T/300 K)^−0.5, that is then a common term to several reaction rates and can hence be pre-computed. 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 LAPACK. 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 Sec. 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 coefficient 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 reduced 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^x-type operations.

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 allows for flexiblity in treating variations in the gas compositions, for example due to altered elemental abundances or regions of enhanced ionisation. However, we note that as well as due to the direct radiation field of the star, ionisation can occur as a result of the diffuse radiation field created by recombinations. Since we only conduct radial ray tracing in this work, we cannot follow these photons accurately, and hence we assumed they are reabsorbed on the spot. Although the diffuse 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 difference 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 Sec. 2.2), we could not include here photoionisation reactions of the singly ionised species. Similarly, neither could the code account for multiple ionisation by the Auger-Meitner effect, although the cross-section due to ionisation of inner (K) shell electrons was included (Verner & Yakovlev 1995) and was 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 did not include secondary electron production, either by the Auger effect 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 efficiency of the heating was 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 was 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 X-ray photoevaporation calculations that relied on MOCASSIN was molecular cooling (Wang & Goodman 2017). Hence, we included it here according to tables computed from Omukai et al. (2010) for CO cooling and from the piece-wise functions of Glover & Abel (2008); Glover (2015) for H₂ cooling. While we did not include cooling from H₂O, we do not see it at sufficient levels in the wind for it to be an important coolant.

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

2.1.6 Photoelectric heating

We changed the photoelectric routine following Kamp & 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 & 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. We note that this formulation gives only the heating from silicate grains; in the current code no heating from Polyaromatic Hydrocarbons (PAHs) was 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 should not affect the models presented in this paper (which include little FUV in the spectra), it will need improvement 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). We obtained faster integration times by linearly interpolating a pre-generated 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 computed the dust cooling and heating with the same approach using two precomputed separated tables. Using two tables allows for dividing the heating (positive) and cooling (negative) in the log-log space.

2.2 Species and reactions in our network

Our network was 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 & Goodman 2017): He/H = 0.1, C/H = 1.4 × 10⁻⁴ and O/H = 3.2 × 10⁻⁴ (Savage & Sembach 1996).

We used 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 or ionised gaseous molecules, and finally CO and H₂O 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 H₂ on dust grains, and freeze out/thermal desorption of CO and H₂O onto/from dust grains.

We note that we excluded 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 & 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 frequently detected 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, one therefore needs to compute the population of each atomic level. To do so, we assumed equilibrium between collisional and radiative transitions and solved the linear system formed by the collisional (de-)excitation rates and the spontaneous emission rates. In particular, the relative abundance of the ith level x_i can be determined by $$-x_i\sum _{\ell }\left(n_{\ell }\sum _{i\ne j}{k}_{ij}^{\ell }\right)-x_i\sum _{j<i}{A}_{ij}+\sum _{\ell }\left(n_{\ell }\sum _{i\ne j}{k}_{ij}^{\ell }{x}_j\right)+\sum _{j>i}{A}_{ji}{x}_j=0,$$(1)

where n_ℓ is the abundance of the ℓth collider (ℓ ∈ {e⁻, H, H⁺, ortho − H₂, para − H₂}), $k_{ij}^{\ell }$ is the temperature-dependent collisional rate coefficient that takes into account the transition between the ith and jth levels when (de-)excited by the ℓth collisional partner, and A_ij is the Einstein coefficient 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. ∑_i x_i = 1), it is possible to compute the level population x_i. This system of equations can be solved with a matrix approach as discussed in Sect. 2.1.3.

After x_i are computed, the cooling emission of each transition becomes Λ_ij = n A_ij x_i ΔE_ij, where n is the total abundance of the specific atomic coolant, and ΔE_ij is the energy difference between the ith and jth 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 v_th ~ 1 km s⁻¹ according to $$\sigma =\frac{{\lambda }^3{A}_{21}}{8\sqrt{2}{\pi }^{3/2}{v}_{\text{th\hspace{0.17em}}}}$$ (Osterbrock & Ferland 2006, Eq. (4.45)). For the shorter wavelength optical lines this is typically ≲10⁻²¹ cm², while those of the far-infrared (FIR) lines are somewhat higher at ≲10⁻¹⁷ cm² (the highest being O I 145 μm with 8.6 × 10⁻¹⁸ cm²). However, given the atomic abundances of ~10⁻⁴, the cross-sections per hydrogen atom are ≲10⁻²¹ cm². Therefore, given that soft X-rays typically penetrate to ~10²¹ cm⁻², 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 or rates. As detailed above, we included e⁻, H, H⁺, ortho − H₂, and para − H₂ as colliders wherever data is available. This ensured that we captured 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 different temperatures were available in Draine (2011, Appendix F). For excitation of O and C by H⁺ we used piece-wise fits from Glover & Jappsen (2007). For collisional (de-)excitation of O by H, ortho-H₂, and para-H₂, we used the latest calculations by Lique et al. (2018)² to which we performed our own fits, presented in Appendix A. Level energies, degeneracies, and Einstein coefficients of spontaneous emission A_ij, were taken from the National Institute of Standards and Technology database (Kramida et al. 2022)

At the temperature range of interest ≲10⁴ K, we expect emission from lines with excitation temperatures up to a few times ≲10⁴ 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 included 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 kept 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 effectively 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 normalised this spectrum to a typical L_X = 2.04 × 10³⁰ erg s⁻¹ (as measured between 0.5 − 5 keV) for a solar-mass star as per Picogna et al. (2021). To this, we also added a 5000 K blackbody assuming a stellar radius of 1 M_⊙; this was to ensure reasonable dust temperatures through heating by optical and UV photons. The spectrum is shown in Fig. 1.

Our fiducial spectrum was used throughout our benchmarks, our main simulation, and our explorations of cooling and microphysics. In addition, we ran 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 normalised and added to the blackbody in the same way as the fiducial spectrum.

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

2.5 Dust

The dust number density n_d was assumed to follow a grain size distribution of $\phi (a)=\frac{\mathrm{d}{n}_d(a)}{\mathrm{d}a}\propto a^{-3.5}$ power law between 5 × 10⁻⁷ and 2.5 × 10⁻⁵ cm (Mathis et al. 1977), with bulk density of 3 g cm⁻³ (e.g., see Zhukovska et al. 2008; Grassi et al. 2017). The optical properties were calculated with the Mie scattering theory 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 that 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 & Clarke 2021) so we assume a uniform ISM dust-to-gas mass ratio of 10⁻² everywhere. However, we note that grain growth is expected to occur within 10⁵ 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 effects on the thermal structure. While we test our sensitivity to this in Sec. 4.2 and discuss what a more reasonable value might be in light of our results in Sec. 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 were conducted in a three-stage process using PLUTO (Mignone et al. 2007).

We started 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 was assumed to be well-mixed, with a maximum grain size in the atmosphere of 2.5 × 10⁻⁵ cm and a minimum grain-size of 5 × 10⁻⁷ cm (the same as we assume in PRIZMO). We used this to calculate initial conditions for the wind by evolving the hydrodynamics for 50 orbits at 10 au, with the wind temperatures calculated using the ‘Spec30’ ξ − T relationship of Ercolano et al. (2021). Aside from details of the extent and resolution of the numerical grid, this is effectively 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 used 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 resulted 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 ran PLUTO+PRIZMO with all the hydrodynamic processes disabled in order to perform a ‘burn-in’ for 0.1 Myr in which only the chemical network and temperature were evolved. This allowed us to reach an equilibrium where these quantities were consistent with the density profile. This also allowed us to benchmark the temperatures calculated by PRIZMO against the temperatures prescribed from the ξ − T relationships and identify areas of major difference (see Sec. 4.2).

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

3.2 Hydrodynamics-thermochemistry interface

We performed several necessary modifications to PLUTO to couple it with PRIZMO³.

The hydrodynamics and thermochemistry were solved alternately within PLUTO’s Strang splitting scheme. On even-numbered 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 that lead to a fatal error. To avoid this, we added checks to PLUTO’s ‘ConsToPrim’ function where similar checks are carried out on the energy, pressure, and density. If the tracer mass-fraction ϵ_i > −1 × 10⁻⁶ and the tracer density ρ_i = ϵ_iρ > −ρ_floor the mass-fraction was simply set to 0, otherwise the cell was flagged for recalculation, first using different 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 conducted radial tracing, considering both geometric dilution and the attenuation: $$j_{\nu }^{i+1}=j_{\nu }^{i+1}{e}^{-{\tau }_{\nu }}\frac{r_i^2}{r_{i+1}^2}.$$(2)

The attenuation $e^{-{\tau }_{\nu }}$ calculated by PRIZMO includes absorption of radiation by the dust as well as all relevant gas-phase species. We note that by conducting radial ray tracing, we neglected the scattering of radiation, which can be important for UV photons but is negligible for the X-ray photons that ultimately drive the wind that results. Subsequently, PRIZMO then advanced the chemical network in each cell of the simulation for a total time δt. In so doing, the gas temperature was evolved along with the number densities of each species; the updated values were 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 diffusive problem. Since we only conducted 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 Diffusion (Levermore & 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 followed the typical approach, which (Owen et al. 2010; Wang & Goodman 2017) is to pin the temperature structure to a pre-calculated dust temperature structure (e.g. Chiang & Goldreich 1997; D’Alessio et al. 2001) beyond a certain column density. Therefore at high column densities, we smoothly interpolated 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₂₃ = N_tot/10²³ cm⁻²): $$T_{\text{gas\hspace{0.17em}}}=\{\begin{array}{ll}{T}_{\text{PRIZMO\hspace{0.17em}}}& N_{\text{tot\hspace{0.17em}}}/{\mathrm{c}\mathrm{m}}^{-2}\le {10}^{21}\\ \frac{T_{\text{PRIMO\hspace{0.17em}}}+N_{23}^2{T}_{\text{DIAD\hspace{0.17em}}}}{1+N_{23}^2}& {10}^{21}\le N_{\text{tot\hspace{0.17em}}}/{\mathrm{c}\mathrm{m}}^{-2}\le {10}^{25}\\ T_{\mathrm{D}\mathrm{I}\mathrm{A}\mathrm{D}}& {10}^{25}\le N_{\text{tot\hspace{0.17em}}}/{\mathrm{c}\mathrm{m}}^{-2}.\end{array}$$(3)

We show and discuss an example of how this affects the dust and gas temperatures in Sect. 5.

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

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

3.3 Hydrodynamics setup

We conducted 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 extended 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 was 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 that 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 difference, which can become negative due to discretisation errors) we used PLUTO’s option to instead conserve entropy, which conserves thermal energy (or 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 performed two benchmarks in the limits corresponding to the extremes of literature photoevaporation models: a low-density/high ionisation parameter EUV model from (Wang & Goodman 2017) and a higher-density 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 & Goodman (2017) against MOCASSIN. They demonstrate that the wind temperatures were strongly overestimated by Wang & Goodman (2017) compared to MOCASSIN due to a lack of effective optical lines available for cooling in their model: this could affect 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 & 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⁴ K temperatures predicted by Wang & 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 & Goodman (2017).

We note that deeper into the disc where the dust and gas temperatures are coupled, PRIZMO and Wang & Goodman (2017) continue to find even lower temperatures than MOCASSIN, that cannot result from the additional cooling. This is since, unlike MOCASSIN, PRIZMO and Wang & Goodman (2017) do not calculate the diffuse 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 & 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 & Goldreich (1997), which reduces the extent of the underestimate. Therefore, PRIZMO shows the coolest temperatures; this would affect 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.

Fig. 3 Comparison of the gas temperatures predicted by (L-R) the ξ − T relationship from Ercolano et al. (2021) (cf. Picogna et al. 2021), PRIZMO (this work), and PRODIMO for the Picogna et al. (2021) analogue model. Inset panels show the inner 10 au.

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 different, 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 ther-mochemical code PRODIMO (PROtoplanetary DIsk MOdel⁴, 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 PRODIMO 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 PRODIMO tend to be a bit cooler. Nevertheless, the overall agreement of PRIZMO with PRODIMO is very good, indicating that the choices for the reduced chemical network and simplified heating and cooling are well justified.

The cooler temperatures in this regime suggest that there are cooling channels included in PRIZMO and PRODIMO but not in MOCASSIN. Hence, we then experimented with turning off different (combinations of) processes in PRIZMO to see what drives the difference. In this pursuit, we kept in mind two questions with potentially separate answers: (a) Which is the cooling that destabilises the temperatures from the 10⁴ K temperatures predicted by the ξ − T relationship and (b) which process ends up dominating the cooling in the final equilibrium?

The most obvious candidate to make a difference might be molecular cooling. However, since H₂ and CO are readily photodissociated (and, unlike in later sections, are not replenished 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 off 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 difference 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-H₂. Oxygen is generally found to be a bright, dominant coolant and H will be the most abundant collider in a largely neutral atomic wind; indeed Ercolano & Owen (2016) found that these collisions were necessary to match observed [O I] line luminosities. Therefore, in panel c) of Fig. 4, we consider the effect of switching off the collisions between O and H. We observed that this has something more of an effect on raising the temperature on scales of 10s 100s au, though still not 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 off 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 H₂ 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 H₂ 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 < T_gas, and hence the dust can have a cooling effect. To see the effect of this, we reduced the dust-to-gas mass ratio to a negligible 10⁻¹⁰. This has relatively little effect on its own (Fig. 4, panel e) or in combination with switching off molecular cooling (Fig. 4, panel f). This further underscores the fact that collisional excitation of O by neutral H dominates 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 off collisions between O and H (Fig. 4, panel g), we find approximately isothermal temperatures of several 1000 K comparable to those predicted by the Ercolano et al. (2021) ξ − T relation (very little additional difference 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 H₂ 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⁴ K temperatures, and both processes must be removed from PRIZMO to recover them. The H₂ 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 effect on the temperatures. Once the temperatures are cool enough for some H₂ survival, it is the next most important coolant after collisional excitation 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+H collisions are sufficient to lower the temperature in the intermediate density and 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; Sec. 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⁴ 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).

Fig. 4 Comparison of the temperatures predicted by PRIZMO when various combinations of cooling processes are switched off 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.

Fig. 5 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(Tgas,PRIZMO, TDIAD)/Tdust,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 = 1021, 1023, 1025 cm−2.

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 > 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⁴ K wind becomes transsonic at r ≳ 0.5 r_G ~ 4.5 au (Parker 1958; Owen et al. 2012). While some mass-loss may be driven from within by the enthalpy (Liffman 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 PRODIMO 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⁴ 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, the power law that we find is very close to the T_esc ∝ r⁻¹ escape temperature predicted by Owen et al. (2012), who assumed that nozzle flow effects are negligible in the wind, and thus balanced 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%, they are generally in good agreement. The difference 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²¹ cm⁻², we are in the wind and T_gas > T_dust. For a column density > 10²¹ cm², 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 sufficiently little that the plotted gas and dust temperatures are the same down to a column density 10²³−10²⁴ cm⁻², but start to diverge again beyond that, with T_gas > 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⁵ cm⁻³, more than an order of magnitude lower than previous X-ray photoevaporation models where they peaked at ~3 × 10⁶ cm⁻³. Even though similar 10⁴ 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⁻³, compared to ξ ≳ 10⁻⁴ 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 $$B=\frac{u^2}{2}+\frac{\gamma }{\gamma -1}\frac{P}{\rho }-\frac{G{M}_{\ast }}{r},$$(4)

(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 & 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 minimum opening angle for a magnetocentrifugal wind (Blandford & 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 when immediately interpreting it as a signature of an MHD wind.

Fig. 6 Quantities extracted along the sonic surface. We only show r > 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.

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 is the sonic surface.

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 H₂. As one would expect, the underlying 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 H₂ present in the outer parts of the wind. Figure 6 shows that its abundance rapidly grows to around 10% 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 H₂ than at the end of the burn-in step, despite the lower density wind conditions. This is a surprise, since H₂ can be photodissociated by FUV & 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 the following reactions resulting in the thermal dissociation of H₂ into atoms or ions: $${\mathrm{H}}_2+\mathrm{H}\to 3\mathrm{H}$$(5) $${\mathrm{H}}_2+{\mathrm{e}}^{-}\to 2\mathrm{H}+{\mathrm{e}}^{-}$$(6) $${\mathrm{H}}_2+{\mathrm{H}}_2\to 2\mathrm{H}+{\mathrm{H}}_2$$(7) $${\mathrm{H}}_2+{\mathrm{H}\mathrm{e}}^{+}\to \mathrm{H}+{\mathrm{H}}^{+}+\mathrm{H}\mathrm{e}.$$(8)

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 H₂ molecule is $k_{\mathrm{P}\mathrm{D}}=3.7\times {10}^{-7}{\mathrm{s}}^{-1}\frac{L_{\mathrm{F}\mathrm{U}\mathrm{V}}}{{10}^{30}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-2}$, corresponding to an expected lifetime of around 25 years at 17 au or 2800 years at the outer edge of the wind. However, the self-shielding of H₂ can help protect it once it accumulates a column of ≳10¹⁴ cm⁻² (Draine & Bertoldi 1996; Richings et al. 2014). Given the approximate density of 10⁵ cm⁻³ and a molecular fraction of 10⁻⁶, inside the H-H₂ transition, the column density can only be N ≈ 10⁻⁶ × 10⁵ cm⁻³ × 10 au = 10¹³cm⁻². Only once the wind reaches a molecular fraction of 10⁻¹ outside the H-H₂ transition, can a column density of N ≈ 10¹⁴ cm⁻² be reached in the space of only 10⁻³ au. Hence, self-shielding can contribute to the steady increase in H₂ in the wind at 10s au but not the initial onset.

Moreover, for thermal dissociation, given the abundances in the inner disc and the 10⁴ K temperatures involved, we expect neutral hydrogen to be the most important destruction route. This provides a thermal dissociation rate per H₂ molecule of k_th = 5.2 × 10⁻¹¹ cm³s⁻¹ n_H; given the gas densities ~10⁵ cm⁻³, which will be dominated by neutral H, a thermal dissociation rate of k_th ≲ 6 × 10⁻⁶ s⁻¹ 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 H₂.

Since self-shielding cannot explain the onset of H₂, non-equilibrium chemistry must therefore play a role. We post-process the simulation with PRIZMO for an additional ~6500 years with hydrodynamics switched off (Fig. 8, left); we then expect that multiple photodissociation timescales will have elapsed everywhere and the H₂ 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 H₂ is photodissociated and confirming that the molecular abundances are out of equilibrium.

We now estimate the hydrodynamic timescale to show that advection of H₂ into the wind is responsible for its survival. Material 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_{\mathrm{S}}(r)=\sqrt{\frac{G{M}_{\ast }}{2r}}$ (Owen et al. 2012). This corresponds to a rate k_ad ~ c_S/r ≈ 1.4 × 10⁻⁷ s⁻¹ (r/au)^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 H₂ 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 H₂ in the wind on scales of 10s au. We 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 H₂ 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 H₂ 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 expected thermal dissociation rate assuming T(r) = T_esc(r)⁶ and n_H(r) ≈ 10^5.5 (r/au)⁻¹, in order to find the location of the H-H₂ transition. By equating the advection and photodissociation rates, we can demonstrate analytically that H₂ would survive against photodissociation alone at $$\begin{array}{rl}\frac{r_{\mathrm{H}2}}{\mathrm{a}\mathrm{u}}& ={\left(\frac{3.7\times {10}^{-7}}{1.4\times {10}^{-7}}\right)}^2{\left(\frac{L_{\mathrm{F}\mathrm{U}\mathrm{V}}}{{10}^{30}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}}\right)}^2\\ & =7{\left(\frac{L_{\mathrm{F}\mathrm{U}\mathrm{V}}}{{10}^{30}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}}\right)}^2,\end{array}$$(9)

as seen in Fig. 9. However, at this location, thermal dissociation is dominant and so the transition does not occur here but is 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 H₂ in the simulation (Fig. 6). While photodissociation here dominates the destruction of H₂ 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 H₂ 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 H₂) – and less abundant than H₂ 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 H₂ by the Class II phase. For O, H₂O ice becomes surprisingly abundant 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 effectively ionised by highly penetrating FUV photons.

Fig. 8 Effects 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 is the sonic surface.

Fig. 9 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, 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.

5.1.3 Main heating and 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%) H₂ enters the wind, at 20–30 au, it dominates the cooling with no more than 1% coming from other processes. However, we observe that despite H₂ 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 H₂ 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 FIR 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 H₂ – 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 H₂ suddenly experiences hotter gas where it is 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%. This is less than what has been found by, for example, Wang & 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 & 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 H₂ cooling is possible, to where the wind is completely atomic, H₂ 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 effective (at temperatures $T<\frac{hc}{k\lambda }$ where λ is the line wavelength), adiabatic cooling becomes relatively more important, rising to around 30% 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 and the narrow strip of significance across the molecular-atomic transition running diagonally across the inset. The latter appears to widen 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 $\xi =\frac{L_X}{n{r}^2}$. In order to understand the salient differences 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 < 10²⁵ cm⁻² (Eq. (3)) and for the analogue model those with N < 2 × 10²² cm⁻² and ξ > 10⁻⁸. In both cases, we observe a tight correlation between ξ and T for ξ ≳ 10⁻⁶; for the fiducial model somewhat more scatter occurs above the main trend for ξ ≲ 10⁻⁶ due to cells with 10²¹ < N / cm⁻² < 10²⁵ which are interpolated between T_PRIZMO and T_DIAD > T_PRIZMO.

For the Picogna et al. (2021) analogue model, there are very few points in the range 500–2000 K that 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, $${\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\theta }\right)}_r>0.$$(10)

Moreover, given that the wind material accelerates upwards, we also expect an upwardly decreasing density $${\left(\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }\theta }\right)}_r>0.$$(11)

Combining these conditions with the ideal gas law leads to the following requirement $${\left(\frac{\mathrm{\partial }\ln T}{\mathrm{\partial }\ln n}\right)}_r>-1,$$(12)

which, given the definition of ξ, we can also express as $${\left(\frac{\mathrm{\partial }\ln T}{\mathrm{\partial }\ln \xi }\right)}_r<1.$$(13)

Thus, we see that if the temperature increases too steeply with increasing ξ then it leads to an increasing pressure that 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 efficient 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⁴ K part of the wind (where H₂ 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 efficient. In this regime, even a small amount of additional heating cannot be offset 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 efficiently excited, a small perturbation in temperature results in a large change in cooling and many different 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 H₂ 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 $\frac{\mathrm{\partial }\ln (T)}{\mathrm{\partial }\ln (\xi )}=0.79$. This means that inequality Eq. (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⁻⁷ to ξ ~ 10⁻⁵ 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 different, a similar effect is known for thermal disc winds from X-ray 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 & Proga 2015; Higginbottom et al. 2017).

The retention of a tight correlation between ξ and T is somewhat to be expected, since, as in the case of atomic cooling previously assumed in the ξ − T relationships, the H₂ cooling happens via two-body processes. Balancing a two-body cooling rate $\mathrm{\Gamma }=n_{{\mathrm{H}}_2}{n}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}f(T)$ with the optically thin X-ray photoionisation heating rate $\mathrm{\Lambda }=\frac{n\overline{\sigma }{L}_X}{4\pi r^2}$ for some spectrally averaged cross-section $\overline{\sigma }$ allows one to find the temperature as T = $f^{-1}\left(\frac{\overline{\sigma }}{4\pi }\frac{L_X}{n{r}^2}\right)=f^{-1}\left(\frac{\overline{\sigma }}{4\pi }\xi \right)$ assuming that $n_{{\mathrm{H}}_2},n_{\text{coll\hspace{0.17em}}}\sim 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_{{\mathrm{H}}_2}/n$) varying with ξ. This is generally true for the photoionisation level (e.g. Glassgold et al. 2007), but is less obvious for the H₂ 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 this case, since the radiation field is fixed, the two are directly correlated, but some caution should be exercised when extrapolating to a different spectrum. Secondly, both photodis-sociation 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 H₂ introduces a dependence of N^−k (where k = 3/4 is often assumed, e.g. Draine & 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 $\frac{L_{\mathrm{F}\mathrm{U}\mathrm{V}}}{{\left(n{r}^{5/3}\right)}^{3/4}}\sim \frac{L_{\mathrm{F}\mathrm{U}\mathrm{V}}}{L_{\mathrm{X}}}{\xi }^{3/4}{L}_{\mathrm{X}}^{1/4}{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²⁵ cm⁻², downweighting those with N ≥ 10²¹ cm⁻² 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 > 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 use the following fitting function (Ercolano et al. 2021): $$\log T=d+\frac{a-d}{{[1+(\log \xi /c{)}^b]}^m}.$$(14)

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

Fig. 10 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 that are not fixed to the initial TDIAD. Overlaid are various parametrizations of the relationships between the variables: on the left, the yellow dash-dotted lines are the parametrizations for three different 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 and 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).

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 different radii, counting only cells where the Bernoulli constant (Wang & Goodman 2017) is positive (and hence the material is unbound). We 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 off of the boundary (Wang & 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. 11 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(≤80 au) = 1.06 × 10⁻⁹ M_⊙ yr⁻¹. 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(≤80 au) = 1.86 × 10⁻⁸ M_⊙ yr⁻¹ after 50 orbits. This suggests that the previous typical X-ray mass-loss rates of ≳ 10⁻⁸ M_⊙ yr⁻¹ (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 differences in methodology particularly regarding the illuminating spectrum the mass-loss rates are now within a factor two 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, and fitting out to the observed drop-off in the surface density profile at ~160 au suggests a total value three to four times higher at ${\dot{M}}_{\text{tot\hspace{0.17em}}}=4.3178\times {10}^{-9}{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$. This is still approximately nine times lower than the total rate found by Picogna et al. (2021) for the same stellar mass and spectrum.

Fig. 11 Evolution of the mass-loss rate inside 80 au calculated using PLUTO+PRIZMO. Panel a: Comparison between models using the X-ray spectrum of Ercolano et al. (2021, E21) (blue; 200 orbits at 10 au) and that of 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 $\dot{M}(M_{\ast }$) given in Komaki et al. (2021) as the grey dotted line. Panel b: Zoom-in on 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 off the cooling due to collisions between O and H and include negligible dust (turquoise) or turn off 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.

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 argue that the contribution of X-rays to driving the wind was negligible, unlike the results of our simulations. Therefore, we considered if the different shapes of the X-ray spectra employed could be the source of the difference.

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³⁰ erg s⁻¹ between 0.5–5 keV), and the fact that TW Hya has a notably ‘soft’ X-ray spectrum (Kastner et al. 2002) compared to other T Tauri stars (which here manifests as the faster drop-off at E > 1 keV due to the lack of the highest-temperature plasma component; Nomura et al. 2007). The difference 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⁷. 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²⁰ cm⁻² – 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 11 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 different 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.

Fig. 12 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.

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 difference in mass-loss rates between this work and the previous simulations.

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

In the latter case we see that neglecting H₂ cooling boosts the mass-loss rate by only a factor ~2. While important, this cannot fully explain the difference to previous models. Instead, in the former case, we find a much more dramatic difference, 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 allowed 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-to-gas ratio is likely too high (see Sec. 6.2), then the latter of these is the more important reason. However, the advective survival of H₂ 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.

5.2.3 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 determined these in the same way as before, namely by calculating the cumulative mass-loss profile, $\dot{M}(<r)$, fitting a sixth-order polynomial in loglog space, and differentiating $$\begin{array}{rl}\log \left(\frac{\dot{M}(<R)}{{\dot{M}}_{\text{tot\hspace{0.17em}}}}\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,\end{array}$$(15) $$\dot{\mathrm{\Sigma }}=\frac{1}{2\pi R}\frac{\mathrm{\partial }\dot{M}(<R)}{\mathrm{\partial }R}.$$(16)

The parameters of the fit for the final snapshot are given in Table 3; ${\dot{M}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ 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 ${\dot{M}}_{\text{true\hspace{0.17em}}}(\le 80\mathrm{a}\mathrm{u})$ is the equivalent measurement out to 80 au, and ${\dot{M}}_{\mathrm{f}\mathrm{i}\mathrm{t}}(\le 80\mathrm{a}\mathrm{u})$ 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). 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 increasing 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 H₂ is disabled, the mass-loss rates are fairly similar around R ≲ 3 au and at R ≳ 50 au. We note that at 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⁴ K escape temperature. The greatest difference is seen around R = 10au, where H₂ starts to survive into the wind and where the wind is warm enough to excite the H₂ lines leading to additional cooling. At smaller radii there is no H₂ in the wind to provide cooling, while outside 50 au, the wind is too cold to excite it.

Fig. 13 Mass-loss rate per unit area as a function of distance to the star. Left: different 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 off certain areas of microphysics in order to test their impact.

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 & Goodman 2017). However, we find that this term is nearly everywhere subdominant to the atomic and/or molecular line cooling, and radiative thermal equilibrium is a good approximation⁸. This is supported by 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 H₂, 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 & 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 effect 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 different spectra (Ercolano et al. 2021) – a future intermediate possibility could be to produce a set of ξ − T relationships at differing H₂ abundances to simplify at least the temperature calculation.

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

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 & Clarke 2021; Booth & 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 affect the thermochemistry in the wind as discussed in Sec. 4.2; we therefore considered how self-consistent our simulations are.

Using our fitted mass-loss profile to the fiducial model in Eq. (25) of Booth & 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 & Clarke (2021) calculate 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.

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⁻³ and u_frag = 10² cm s⁻¹, and assuming a Lynden-Bell & Pringle (1974) disc profile with M_D = 0.1 M_⊙, R_C = 100 au and T1au = 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 distribution, then the entrained fraction is $f_{\mathrm{e}\mathrm{n}\mathrm{t}}={\left(\frac{a_{\mathrm{m}\mathrm{a}\mathrm{x}}}{a_{\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{g}}}\right)}^{1/2}$. For the Picogna et al. (2019) mass-loss profile, this estimate suggests up to ~10% entrained mass (or a dust-to-gas ratio of 10⁻³), consistent with the modelling by Franz et al. (2022a). Meanwhile, our mass-loss profile suggests values up to five times lower, largely in the range ϵ = 10⁻⁴−10⁻³. 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 and population synthesis

One of the major uses of photoevaporation prescriptions, such as those 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 difficult 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 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.

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⁻⁸ M_⊙ yr⁻¹, 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) show that for 0.3 ≤ M_*/M_⊙ ≤ 1.2 if the mass-loss rate exceeds ~10⁻⁹ M_⊙ yr⁻¹ 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 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 M_J, 100 au discs around solar-mass stars with undergoing viscous evolution (α = 10⁻³) 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% 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’.

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 & 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 μm line 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 & Owen 2010; Ercolano et al. 2021), suggesting that any effect 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).

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 X-rays 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).

Finally, the fact that molecules particularly H₂ 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 on-the-fly use 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 approach to be sufficiently flexible and comprehensive to capture thermochemistry in ionised, neutral atomic, and molecular environments at a range of temperatures and (b) testing spectra used in different previous works in a controlled manner, we have ensured that we could confidently isolate the root of the differences we found and probe the causes of the different claims in the literature. In particular, we find that

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 being a result of the very steep relationships between the temperature and ionisation parameter ξ used therein, which leads 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 FIR regimes; molecular cooling in the mid-IR stabilises the intermediate temperatures, leading to a 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 H₂ dissociation front (Fig. 6) and in 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 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⁻⁹ M_⊙ yr⁻¹ (Fig. 11a). The mass-loss is reduced everywhere by a similar amount compared to previous works (Fig. 13).

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 underestimation of the X-ray contribution to photoevaporation (Fig. 11b).

6. The winds contain a significant amount of molecular hydrogen. Advection of material into the wind (i.e. non-equilibrium 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 understand as resulting from the additional cooling provided (in particular by excitation of O by neutral H), help resolve longstanding 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 previously expected, 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 suggested by previous X-ray models. 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 that can be used in disc and planet population synthesis models.

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).

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 & 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. A.D.S. also thanks James Matthews for useful conversations about the similarities to thermal winds from X-ray binaries. A.D.S. 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. 101019751 MOLDISK). C.J.C. acknowledges support from the Science & Technology Facilities Council (STFC) Consolidated Grant ST/W000997/1 and European 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 Facilities 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.

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 H₂ from Lique et al. (2018), in the form $$k_{\mathrm{u}\mathrm{l}}=(a\times T_2^{b+c\ln T_2}+d\exp (-e/T_2))\times {10}^{-10}{\mathrm{c}\mathrm{m}}^3{\mathrm{s}}^{-1},$$(A.1)

where $T_2=\frac{T}{{10}^2\mathrm{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 $\frac{\mathrm{d}n}{\mathrm{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⁻³ 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.

Fig. B.1 Distribution of the execution time per cell of PRIZMO after 50 orbits at 10 au.

References

All Tables

All Figures

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

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

	Fig. 3 Comparison of the gas temperatures predicted by (L-R) the ξ − T relationship from Ercolano et al. (2021) (cf. Picogna et al. 2021), PRIZMO (this work), and PRODIMO for the Picogna et al. (2021) analogue model. Inset panels show the inner 10 au.
In the text

	Fig. 4 Comparison of the temperatures predicted by PRIZMO when various combinations of cooling processes are switched off 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.
In the text

Fig. 5 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(Tgas,PRIZMO, TDIAD)/Tdust,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 = 1021, 1023, 1025 cm−2.

In the text

Fig. 6 Quantities extracted along the sonic surface. We only show r > 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.

In the text

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 is the sonic surface.

In the text

Fig. 8 Effects 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 is the sonic surface.

In the text

	Fig. 9 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, 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.
In the text

Fig. 10 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 that are not fixed to the initial TDIAD. Overlaid are various parametrizations of the relationships between the variables: on the left, the yellow dash-dotted lines are the parametrizations for three different 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 and 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).

In the text

Fig. 11 Evolution of the mass-loss rate inside 80 au calculated using PLUTO+PRIZMO. Panel a: Comparison between models using the X-ray spectrum of Ercolano et al. (2021, E21) (blue; 200 orbits at 10 au) and that of 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 $\dot{M}(M_{\ast }$) given in Komaki et al. (2021) as the grey dotted line. Panel b: Zoom-in on 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 off the cooling due to collisions between O and H and include negligible dust (turquoise) or turn off 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.

In the text

	Fig. 12 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.
In the text

Fig. 13 Mass-loss rate per unit area as a function of distance to the star. Left: different 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 off certain areas of microphysics in order to test their impact.

In the text

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

	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).
In the text

	Fig. B.1 Distribution of the execution time per cell of PRIZMO after 50 orbits at 10 au.
In the text