© The Authors 2025

1 Introduction

In recent years, significant observational progress has been made in detecting and characterising exoplanets. It is well known that planetary systems exhibit a large diversity in their architectures and planetary properties. However, the exact pathways leading to the observed final states are still under investigation (for a review, see Drążkowska et al. 2023). Since only a few observations of systems with forming planets exist (Keppler et al. 2018; Müller et al. 2018), it is crucial to study their formation environments through simulations. This approach allows us to gain greater insight into the structure and composition of protoplan-etary discs and therefore helps us constrain the properties and formation pathways of planets.

Planetary cores are believed to form via solid accretion in the so-called core accretion scenario. The solids can either be accreted in the form of planetesimals (Ida & Lin 2004; Guilera et al. 2014; Mordasini 2018; Miguel et al. 2019; Emsenhuber et al. 2021) or pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015, 2019; Ndugu et al. 2018; Lambrechts et al. 2019; Liu et al. 2019; Izidoro et al. 2021), with pebble accretion being more efficient than planetesimal accretion in the outer disc regions (e.g. Johansen & Bitsch 2019). Many population synthesis studies have used the core accretion scenario in order to reproduce observed masses and orbital distances of exoplanets. To further constrain planet formation models, we can now also use newly measured atmospheric abundances (e.g. Mollière et al. 2020; Line et al. 2021; Pelletier et al. 2021; August et al. 2023). This has been done in recent theoretical studies by Mordasini et al. (2016), Booth et al. (2017), Schneider & Bitsch (2021a,b), Bitsch et al. (2022), Mollière et al. (2022), and Penzlin et al. (2024).

The final composition of a planet depends crucially on the composition of the disc. Therefore, simulations of protoplane-tary discs need to include the full evolution of the disc’s structure and chemical composition. Many factors play a role in determining the latter, such as disc temperature and density, radiation fields, and the position of evaporation lines of different species (see e.g. Cuzzi & Zahnle 2004; Öberg et al. 2011; Henning & Semenov 2013; Schneider & Bitsch 2021a; Eistrup & Henning 2022; Mollière et al. 2022). The structure of the disc is influenced by pressure bumps that can block inward drifting pebbles (e.g. Pinilla et al. 2012), where the latter again plays a role in the time evolution of the disc’s composition. Indeed, James Webb Space Telescope (JWST) observations by Banzatti et al. (2023), Grant et al. (2023), Tabone et al. (2023), and Gasman et al. (2025) have shown that pressure bumps influence the chemical composition of inner protoplanetary discs. This effect is also seen in simulations, as demonstrated in works by Bitsch et al. (2021), Kalyaan et al. (2023), and Mah et al. (2024).

One possible mechanism for causing these pressure bumps is a planet reaching its pebble isolation mass and creating a gap in the protoplanetary disc (e.g. Paardekooper & Mellema 2006; Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018). Two other mechanisms that significantly influence the composition and evolution of protoplanetary discs are magnetohydrodynamic (MHD) winds (for a review, see Lesur et al. 2023) or photoevaporation (for a review, see Pascucci et al. 2023). Both of these processes carry away disc material, with photoevaporation carving deep gaps into the disc. Since photoevaporation determines the final stages of disc evolution (e.g. Ercolano et al. 2009; Pascucci & Sterzik 2009; Ercolano & Clarke 2010; Owen et al. 2012, 2013), it is especially important, given that observations indicate that discs on average only live for a few million years (Mamajek et al. 2009; Fedele et al. 2010) and therefore require a disc dissolving mechanism. The specific lifetime of discs, however, can vary significantly depending on the properties of the central star of the system (e.g. Michel et al. 2021; Pfalzner et al. 2022).

The process of photoevaporation describes the transfer of energy from high-energy photons to gas particles in the disc, causing the latter to increase their velocities to the point that they escape the system. Depending on the source of the energetic radiation, we either refer to the evaporation process as external photoevaporation (for a review, see Winter & Haworth 2022), when the radiation comes from nearby stars outside the studied system, or internal photoevaporation (for a review, see Pascucci et al. 2023), when the radiation’s origin is the host star of the studied system. After enough gas has escaped the protoplanetary disc, a gap opens, resulting in a pressure bump in the disc. The gap blocks inward-moving pebbles and, at the same time, cuts off the inner regions from their gas supply from the outer disc as the photoevaporative winds continuously carry the gas away.

This work is a follow-up on our first paper, see Lienert et al. (2024) (hereafter, referred to as Paper I), where we combined a viscous disc evolution model, including pebble drift and evaporation, with a model of internal photoevaporation and studied protoplanetary discs around solar-mass stars. Here, we extend this study to discs around lower-mass stars (M_★ = 0.1-0.5 M_⊙), especially focusing on the influence of the photoevaporative mass loss rate on the evolution of the elemental composition of the inner disc, which occurs on faster timescales compared to Sun-like stars (e.g. Mah et al. 2023). Here, and for the rest of the paper, the ‘inner disc’ is defined as the part of the disc close to the star, separated from the disc’s outer parts by the photo-evaporative gap. The gap’s inner edge, which defines the outer boundary of the inner disc, is generally located between 0.73 AU, depending on the stellar mass and the photoevaporation rate.

This paper is structured as follows: Section 2 includes a summary of our model. The obtained results are shown in Section 3, with their implications discussed in Section 4, before we conclude with Section 5.

2 Methods

Our numerical simulations were carried out with the code chemcomp, as described in Schneider & Bitsch (2021a), with added internal photoevaporation using the model of Picogna et al. (2021). Chemcomp is a 1D, semi-analytical model of protoplanetary discs combined with a planetary growth model. For more details on the planet formation model, we refer to Schneider & Bitsch (2021a). In this study, we mostly used the disc module of the code, where the disc evolution follows a classic viscous evolution model (see e.g. Lynden-Bell & Pringle 1974), using the α-viscosity description of Shakura & Sunyaev (1973).

Dust grains in the disc can grow into pebbles, with their evolution being modelled by the two-population approach developed by Birnstiel et al. (2012). This approach was originally calibrated for solar-mass stars (Birnstiel et al. 2012, see also the discussion in Pfeil et al. 2024). Without the calibration, our model might harbour small uncertainties in the grain growth and drift routine for lower-mass stars. The drift timescale of particles is directly dependent on their size. The particle’s size, in turn, depends on the fragmentation velocity and the disc’s viscosity. Both these parameters have large parameter ranges, with the fragmentation velocity being constrained to 1-10 m/s by experiments (Gundlach & Blum 2014), while the turbulence parameter α, which directly influences the viscosity, has a parameter range of α = 10⁻⁴−10⁻² in usual disc settings (see e.g. Dullemond et al. 2018 for the observational side, or Flock et al. 2015 for the theory perspective). Therefore, the uncertainty in the combination of those two parameters easily exceeds the uncertainty coming from the calibration to a different stellar mass in the two-population approach by Birnstiel et al. (2012).

The growth of the dust grains in the disc is limited by drift, fragmentation, and drift-induced fragmentation. The relative ratio between large and small grains changes depending on the growth-limiting mechanism. It is f = 0.75 for the fragmentation limit and f = 0.97 for the drift limit. The grains are then advected with the average velocity of the two components weighted by their relative fractions. Enrichment of the disc occurs through evaporating volatiles from inward drifting pebbles at the particular ice lines of each species (Schneider & Bitsch 2021a).

2.1 Viscous evolution

The disc’s viscous evolution is described by the viscous disc equation, which is given as the time evolution of the gas surface density Σ_gas and can be derived from the conservation of mass and angular momentum (Pringle 1981; Armitage 2013), $$\frac{\mathrm{\partial }{\mathrm{\Sigma }}_{\text{gas,}Y}}{\mathrm{\partial }t}-\frac{3}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[\sqrt{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\sqrt{r}\nu {\mathrm{\Sigma }}_{\text{gas,}Y}\right)\right]={\dot{\mathrm{\Sigma }}}_Y.$$(1)

ν is the disc's viscosity, and ${\dot{\mathrm{\Sigma }}}_Y$ denotes the source term of the molecular species Y. For a full list of molecules included in the code, we refer to Schneider & Bitsch (2021a) and Paper I. The source term results from the evaporation and condensation of pebbles and is given by $${\dot{\mathrm{\Sigma }}}_Y=\{\begin{array}{c}{\dot{\mathrm{\Sigma }}}_Y^{\text{evap}}r<r_{\text{ice,}Y}\\ {\dot{\mathrm{\Sigma }}}_Y^{\text{cond}}r\ge r_{\text{ice,}Y}.\end{array}$$(2)

Here, ${\dot{\mathrm{\Sigma }}}_Y^{\text{evap}}$ and ${\dot{\mathrm{\Sigma }}}_Y^{\text{cond}}$ are the source terms for evaporation and condensation of species Y. They originate from the evaporation and condensation of volatiles of molecule Y at the respective ice line, r_ice,Y.

The viscosity, ν, in Equation (1) is given by $$\nu =\alpha \frac{c_{\text{s}}^2}{{\mathrm{\Omega }}_{\text{K}}},$$(3)

where α describes the strength of the turbulence, c_s is the isothermal sound speed and ${\mathrm{\Omega }}_{\text{K}}=\sqrt{\frac{G{M}_{\star }}{r^3}}$ is the Keplerian angular frequency. Here, M_★ denotes the mass of the host star and r the radial distance to it. The dimensionless turbulence parameter α is kept fixed in time and has a radially constant value of α = 10⁻⁴. The speed of sound can be linked to the mid-plane temperature of the disc, which does not evolve in time for simplicity and is thus calculated only at the initialisation of the simulation via an equilibrium between viscous and stellar heating with radiative cooling.

One has to keep in mind that in reality, low-mass stars (M_★ = 0.1-0.5 M_⊙) experience a decrease in their luminosity over time, leading to lower irradiation temperatures (Baraffe et al. 2015). As a result, the evaporation lines for different molecules move inwards. However, the innermost evaporation lines might not be affected as their position is determined by viscous heating, which evolves slowly if the viscosity is low (α = 10⁻⁴). In addition, the inward flux of pebbles occurs on much faster timescales than the luminosity evolution of the host star (Brauer et al. 2008). This means that the pebbles have already reached the inner disc and evaporated their volatile content before the shift in the evaporation lines would affect their evaporation. Therefore, keeping the mid-plane temperature constant in time is a reasonable assumption for our simulations, where we used a small turbulence parameter of α = 10⁻⁴, which results in a slow viscous evolution.

2.2 Internal photoevaporation

For the simulations in this paper, the viscous disc is subject to internal photoevaporation. In this case, an additional term ${\dot{\mathrm{\Sigma }}}_{\text{w}}$, describing the photoevaporative loss rate of the gas surface density, is subtracted from the right-hand side of the viscous disc Equation (1), $$\frac{\mathrm{\partial }{\mathrm{\Sigma }}_{\text{gas,}Y}}{\mathrm{\partial }t}-\frac{3}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[\sqrt{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\sqrt{r}\nu {\mathrm{\Sigma }}_{\text{gas,}Y}\right)\right]={\dot{\mathrm{\Sigma }}}_Y-{\dot{\mathrm{\Sigma }}}_{\text{w}}.$$(4)

This additional term acts on the total gas surface density, not separately on each included molecule. The description we employed for photoevaporation due to X-rays follows the work of Picogna et al. (2019, 2021) and Ercolano et al. (2021). Only the soft X-ray regime was used, as it has shown to be the most efficient part of the X-ray spectrum.

The photoevaporative gas surface density loss rate is given by $$\begin{array}{rl}{\dot{\mathrm{\Sigma }}}_{\text{w}}& =(\frac{6a\ln (r{)}^5}{r\ln (10{)}^6}+\frac{5b\ln (r{)}^4}{r\ln (10{)}^5}+\frac{4c\ln (r{)}^3}{r\ln (10{)}^4}+\frac{3d\ln (r{)}^2}{r\ln (10{)}^3}\\ & +\frac{2e\ln (r)}{r\ln (10{)}^2}+\frac{f}{r\ln (10)})\frac{\ln (10)}{2\pi r}{\dot{M}}_{\text{w}}(r)[{\text{M}}_{\odot }{\text{AU}}^{-2}{\text{yr}}^{-1}],\end{array}$$(5)

with a mass loss rate of $$\begin{array}{rl}{\dot{M}}_{\text{w}}(r)=& {10}^{a{\log }_{10}(r{)}^6+b{\log }_{10}(r{)}^5+c{\log }_{10}(r{)}^4+d{\log }_{10}(r{)}^3}\\ \cdot & {10}^{e{\log }_{10}(r{)}^2+f{\log }_{10}(r)+g}\cdot {\dot{M}}_{\text{w}}(M_{\star },L_{\text{X,soft}})[{\text{M}}_{\odot }{\text{yr}}^{-1}],\end{array}$$(6)

where the fit parameters a-g for our stellar sample are presented in Table 1. L_X,soft denotes the soft part of the X-ray luminosity of the star.

Equation (5) is visualised in Figure 1, with the photoevapo-rative gas surface density loss rates plotted as a function of disc radius and stellar mass. The mass loss peaks between 1 AU and 10 AU, depending on the mass of the central star, with the general trend of a shift in the peak location towards larger radii with increasing stellar mass. The disc mass is generally lost between 0.7 AU and 300 AU, enabling the opening of a gap in this disc region. Similar to the trend of the peak location, the gap position moves to larger disc radii as a function of stellar mass.

An exception from both these trends is seen for the 0.3 M_⊙ star. The reason is that the photoevaporation rates result from a combination of two changes, namely the stellar mass and the X-ray luminosity. A change in stellar mass results in a different gravitational potential, which is responsible for keeping the gas close to the host star. The X-ray luminosity, on the other hand, is responsible for the repulsion of the gas. Changing it affects the internal physics of the underlying photoevaporation model. A lower X-ray luminosity, as in the case of the 0.3 M_⊙ star, is associated with a harder spectrum in the soft X-ray regime, which leads to a shift of the peak location towards larger radii (Picogna et al. 2021). Additionally, the X-ray luminosity affects the thermal structure of the disc. Here, lower X-ray luminosities lead to a larger flaring of the disc, which allows mass removal from a larger region (Ercolano et al. 2021). Combining all these effects leads to the non-uniform dependence of the photoevaporative gas surface density loss rates on the stellar mass that we see in Figure 1, because a change in stellar mass is always accompanied by a change in X-ray luminosity. A linear trend for both the peak and gap position of the photoevaporation rates is observed when fixing the X-ray luminosity across the stellar sample; see Figure 9 from Picogna et al. (2021), where the mass loss rate as a function of distance is then just dependent on the star’s gravitational potential.

As we see from Figure 1, the photoevaporation rates vary in radius but are constant in time, leading to the same amount of disc material being taken away in each time step of the simulations. The evolution of other disc areas unaffected by photoevaporation is then solely determined by the viscous transport of gas.

The photoevaporative mass loss rate Ṁ_w(M_★, L_X,soft) in Equation (6) depends on the stellar mass and X-ray luminosity. According to Güdel et al. (2007), the X-ray luminosity L_X of a star can be calculated from its mass by using the following relation, $${\log }_{10}\left(L_{\text{X}}\right)=(1.54\pm 0.12){\log }_{10}(M_{\star })+(30.31\pm 0.06)[{\text{erg s}}^{-1}],$$(7)

where the stellar mass has to be given in multiples of the solar mass M_⊙. Employing Equation (7), we obtain the X-ray luminosities for our stellar sample, as given in the right column of Table 1 (see also Table 1 in Picogna et al. 2021).

Since only the soft part of the X-ray spectrum (0.1-1 keV) is important for photoevaporation, we used the following scaling relation to calculate the soft X-ray luminosity from the full X-ray luminosity, $${\log }_{10}\left(L_{\text{X,soft}}\right)=0.95{\log }_{10}\left(L_{\text{X}}\right)+1.19[{\text{erg s}}^{-1}].$$(8)

This relation is a result of interpolating the values of Table 4 in Ercolano et al. (2021).

The photoevaporative mass loss rate Ṁ_w(M_★, L_X,soft) in Equation (6) consists of two parts. One is its stellar massdependent part and the other is its X-ray-dependent part. To calculate the first, stellar mass-dependent part, we employed Equation (5) from Picogna et al. (2021), $$\dot{M}(M_{\star })=3.93\cdot {10}^{-8}{M}_{\star }[{\text{M}}_{\odot }/\text{yr}].$$(9)

The second, X-ray-dependent part was taken from Ercolano et al. (2021), $${\log }_{10}(\dot{M}(L_{\text{X,soft}}))=A\exp \left(\frac{{(\ln ({\log }_{10}\left(L_{\text{X,soft}}\right))-B)}^2}{C}\right)+D[{\text{M}}_{\odot }/\text{yr}].$$(10)

The parameters are given as A = −1.947 · 10¹⁷, B = −1.572 · 10⁻⁴, C = −2.866 · 10⁻¹, D = −6.694.

The combined mass loss rate is then given by a scaling relation because the stellar mass dependence needs to be rescaled by the difference in the soft X-ray component adopted with respect to the mean X-ray luminosity for a given stellar mass, $$\dot{M}(M_{\star },L_{\text{X,soft}})=\dot{M}(M_{\star })\frac{\dot{M}(L_{\text{X,soft}})}{\dot{M}(L_{\text{X,soft,mean}})}[{\text{M}}_{\odot }/\text{yr}].$$(11)

Here, L_X,soft is the soft part of the observed X-ray luminosity of any given input star, with the observational value taken from Figure 2 and the soft part calculated via Equation (8). To obtain L_X,soft,mean, we used the mass of our input star and calculate its mean X-ray luminosity by employing Equation (7). The soft part of that mean X-ray luminosity was then calculated once more using Equation (8).

Since there is a large spread in X-ray luminosities for one specific stellar mass, the actual X-ray luminosity of a star can differ massively from the mean value calculated with Equation (7); see Figure 2. This results in a spread in the photoevapora-tive mass loss rates as well, which might change the lifetimes of the protoplanetary discs. To investigate how and to what extent the disc lifetimes change, we reduced the X-ray luminosities by a factor of [5.0, 3.5, 2.2] compared to the mean values for our stellar sample of [0.5, 0.3, 0.1] M; see Table 2. These reduced X-ray luminosities are well within the spread of the observed X-ray luminosities of Güdel et al. (2007); see simulation data in Figure 2. The decrease in stellar X-ray luminosity leads to a reduction of the photoevaporative mass loss rates by a factor of 3-4 compared to the nominal values. A decrease in the strength of internal photoevaporation prolongs the lifetimes of the protoplanetary discs, which is more in agreement with observations.

The simulations by Sellek et al. (2024) show lower photoevaporation rates compared to Picogna et al. (2021) due to the inclusion of additional cooling processes. However, these studies have so far only been applied to solar-type stars. We thus probed reduced photoevaporation rates in discs around lower-mass stars by changing the X-ray luminosity, as explained above.

An additional factor not considered here is the time variability of X-ray luminosities. A star’s luminosity is not constant over its lifetime, leading to a time-dependent photoevaporative mass loss rate. Furthermore, assuming a different disc structure, i.e. a puffed-up inner part, results in the extinction of X-ray radiation, hindering internal photoevaporation. The puff-up is a result of viscous heating, which is accounted for in our model chemcomp, but not in the photoevaporation models by Picogna et al. (2019), who only incorporate stellar heating. Therefore, we cannot study how much photoevaporation would be hindered by puffed-up discs.

Fig. 1 Photoevaporative gas surface density loss rate as a function of disc radius for different stellar masses, as given in Equation (5), adopted from Picogna et al. (2021). Here, the nominal values for the photoevap-orative mass loss rate, as given in Table 1, are used.

Fig. 2 X-ray luminosity, derived from the measured X-ray flux, as a function of stellar mass for all detected XEST sources, re-plotted from Güdel et al. (2007) using their original data. The distance to Taurus is estimated at 140 pc, although we note that Taurus consists of subgroups with varying distances (Galli et al. 2019). The flux of stars with more than one measurement is averaged to retrieve only one data point for the plot. The circle, star, and diamond symbols indicate the stellar class, as defined for XEST sources; and the straight line gives a linear regression for the logarithmic values, as described by the following equation, log(LX) = 1.54 log(M★) + 30.31. Synthetic values from our simulations from Sections 3.2 and 3.3 are added as purple-coloured squares, with the colour indicating the gap opening time of the disc, tgap.

2.3 Pebble evaporation and chemistry model

Protoplanetary discs in our model consist of gas and dust grains. The grains can grow to pebbles with sizes ranging from millimetres to centimetres. Drift, turbulent fragmentation, and drift-induced fragmentation limit their growth (Brauer et al. 2008). Both gas and dust move inwards through the disc, with the dust pebbles drifting faster due to gas drag in smooth discs without substructures (Mah et al. 2024). In discs containing pressure bumps, pebbles move towards pressure maxima. On their inward journey through the disc, they cross several ice lines of different molecular species and evaporate their respective volatile content. Consequently, the inner disc is enriched with the evaporated material. When gas gets closer to the central star, it will eventually be accreted onto it. On the other hand, re-condensation at evaporation fronts is also possible for gaseous species when they are moving outwards rather than inwards (Ros & Johansen 2013). This can lead to large enhancements of solids compared to the solar composition (e.g. Aguichine et al. 2021; Mousis et al. 2022; Mah & Bitsch 2023). Inward drifting pebbles are self-consistently stopped at pressure bumps.

For the initial composition of the disc, we assumed solar abundances. The exact values can be found in Paper I and are taken from Asplund et al. (2009). The elements were then distributed into volatile and refractory molecules, with CO, CO₂, and CH₄ as the main carbon-bearing species and water as the main oxygen-bearer. We note that C₂H₂ may be an important molecular carrier of carbon in the gas phase as well (Grant et al. 2025), but is currently not included in our model. The distribution of elements X into molecules Y follows a simple chemical partitioning model from Schneider & Bitsch (2021a), which is based on Bitsch & Battistini (2020). The disc composition is dependent on the disc radius since the disc temperature varies with the orbital distance from the central star. Depending on their condensation temperature and position in the disc, molecules from a certain species Y are then either present in gaseous or solid form. The orbital distance where the disc’s mid-plane temperature is equal to the evaporation (condensation) temperature of a certain species is called the evaporation (ice) line of that species.

The number ratio of two elemental species X₁ and X₂ is defined via $$X_1/X_2=\frac{m_{X_1}}{m_{X_2}}\frac{{\mu }_{X_2}}{{\mu }_{X_1}},$$(12)

where m_X1 and m_X2 are the mass fractions of the two elements and μ_X1 and μ_X2 are their atomic masses. In this paper, this definition is used to calculate the C/O, C/H, O/H, and N/H ratios.

2.4 Initial conditions

Our standard simulations were carried out using host star masses of M_★ = [0.5, 0.3, 0.1] M_⊙. The corresponding stellar luminosities were retrieved from the stellar evolution models of Baraffe et al. (2015); for this work, we chose the luminosities at 2.5 Myr to be consistent with Paper I. The resulting values for the luminosities are L_★ = [0.34, 0.17, 0.03] L¹. The initial disc in our standard simulations is characterised by a turbulence parameter of α = 10⁻⁴, an initial disc mass of M_disc = [0.05,0.03,0.01] M_⊙, an initial disc radius of R_disc = [95, 65, 30] AU, and an initial dust-to-gas ratio of 2%, as listed in Table 3. The initial disc mass was fixed at 10% stellar mass under the assumption that our simulated systems began to evolve at an earlier stage when the disc is thought to be more massive than what is seen in observations. The values for the initial disc radius were obtained using the scaling relation of Mah et al. (2023) based on Andrews et al. (2018), where R_disc = 50 AU (M_★/0.1 M_⊙)^0.7. To match the parameters from our first paper, we scaled down the ones calculated with the above-mentioned relation to arrive at the ones listed before. The initial disc radius is generally defined as the point in the disc where the exponential cut-off of the initial gas surface density sets in. The latter is described by a profile with an exponential decay in the outer disc regions, starting at R_disc; see Appendix A for the gas surface density evolution of a purely viscous disc. We ran all our simulations for 10 Myr with a time step of 10 yr. The disc composition was retrieved every 0.1 Myr.

3 Results

In this work, we studied the influence of internal photoevaporation on the (chemical) evolution of discs around low-mass stars, i.e. stars with masses of M_★ = 0.1-0.5 M_⊙. We present our results in this section, focusing mainly on the contrast between the nominal values for photoevaporative mass loss from Picogna et al. (2021) and reduced mass loss rates. The results for the nominal values are found in Section 3.1, whereas the results for the reduced rates are presented in Section 3.2. Section 3 closes with comparing the disc lifetimes for different photoevaporation rates in Section 3.3.

For comparison, the appendices contain additional material on non-photoevaporative discs; see Appendix A for the case of a purely viscously evolving disc or Appendix B for a disc with a giant planet forming in it. Additionally, we provide more information on photoevaporative discs in Appendix C, where we study a reduction factor of 10 in the photoevaporative mass loss, as well as a photoevaporative disc with a viscosity increased by a factor of 10.

Fig. 3 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from Table 1. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius, time evolution is shown in colour - from black, which corresponds to 0 Myr, to dark red, which corresponds to 10 Myr. Bottom : gaseous C/O ratio as a function of disc radius and time (colour-coded). The evaporation lines for the different molecules are given as dashed grey lines. Note that the C/O ratio is calculated from number densities and that, by definition, we have no specified C/O ratio in the gas phase beyond the CO evaporation front. We use our standard parameters for this simulation, as given in Table 3.

3.1 Nominal values for photoevaporative mass loss

3.1.1 Disc evolution

All results shown in this subsection are obtained using the standard simulation parameters as described in Section 2.4 and shown in Table 3. For the photoevaporative mass loss rates, we use the nominal values given in Table 1.

Figure 3 shows the gas surface density in the top row and the gaseous C/O ratio in the bottom row, both as a function of disc radius, for a viscous disc with active internal photoevaporation. Stellar masses are varied across the different columns, with the mass decreasing from M_★ = 0.5 M_⊙ on the left to M_★ = 0.1 M_⊙ on the right. The time evolution spans over 10 Myr, indicated via colour coding. The C/O ratio is calculated from number densities.

Our results for the gas surface density align with those from Paper I, where a similar study was conducted for solar-mass stars. The gas surface density decreases over time in the inner disc due to viscous accretion of disc material onto the host star. In the outer disc regions, viscous spreading is the dominating mechanism. Both of these processes occur on relatively long timescales as the viscous parameter α is small, α = 10⁻⁴.

Additionally, internal photoevaporation opens a gap between 1-2 Myr. Its inner edge moves inwards with time for all host star masses shown here. The exact gap size and position depend on the particular shape of the gas surface density loss rate; see Figure 1. Photoevaporation does not necessarily have the strongest visible effect at the position of the peak of the photoevaporative mass loss rate. Since the disc’s surface density decays exponentially in its outer regions, it is important to check the relative strength of photoevaporation compared to the disc’s gas surface density.

In comparison to paper I, the inner disc’s gas surface density decreases faster for lower-mass stars, especially for the star with M_★ = 0.1 M_⊙. This is a direct result of the lower initial disc masses.

3.1.2 Evolution of the C/O ratio

The gap opened by photoevaporation has a large impact on the inner disc’s chemical evolution. It blocks the pebbles from the outer disc regions, and at the same time, photoevaporative winds carry away inward-moving gas. The direct results of this mechanism are seen in the C/O ratio; see the lower panels of Figure 3.

The initial C/O ratio in Figure 3 shows the classical steplike behaviour (e.g. Öberg et al. 2011; Mollière et al. 2022), with the different volatiles evaporating at different distances. These positions in the disc are determined by the disc temperature. The result is either an increase in the C/O ratio at the evaporation lines of carbon-rich molecules or a decrease in the C/O ratio at those of oxygen-rich molecules.

As the disc starts to evolve, a rapid drop in the C/O ratio in the inner disc, interior to the water-ice line, is seen during the first one million years for all host star masses shown here. The C/O ratio then continues to stay subsolar for the rest of the time evolution, in contrast to the cases of a purely viscously evolving disc or a disc hosting a giant planet, where the C/O ratio reaches supersolar values at the end of the time evolution (see Appendices A, B and the study by Mah et al. 2023). However, apart from the C/O ratio generally being subsolar for all host star masses in the case of the photoevaporative discs that we investigate here, the rest of the time evolution differs between the three different central star masses.

M_★ = 0.5 M_⊙: in this case, the C/O ratio continues to decrease after the initial rapid drop. Its behaviour is a direct consequence of pebble drift and evaporation. The position of the water evaporation line, marked with number 2 in Figure 3, is at a smaller disc radius than those of the carbon-bearing molecules CO₂, CH₄, and CO (marked with 4, 5, and 6). Consequently, pebbles from outer disc regions first reach the evaporation lines of the carbon-bearing species before crossing the water-ice line when moving inwards. Therefore, the carbon-rich vapour is created farther out in the disc than the water vapour. Additionally, pebbles move much faster than gas, resulting in water-ice pebbles reaching the inner disc much earlier than the carbon-rich gas. As a result, the inner disc is first enriched with water vapour, leading to the rapid drop in the C/O ratio in the first one million years due to the large amount of oxygen contained in water. With time, the carbon-rich gas that formed in the outer disc regions enters the inner parts. At the same time, inner disc material is accreted onto the host star. Before a visible effect of these processes can occur, a gap opens in the disc between 1 Myr and 2 Myr. The aforementioned processes stop since neither pebbles nor gas can pass the deep gap carved by photoevaporation. Instead, another mechanism comes into play.

After the gap is fully opened, the flux of carbon vapour from the outer disc is blocked. On the other hand, the gas already present in the inner disc is accreted onto the host star, but with one exception. This is due to the occurrence of a water equilibrium cycle at the water evaporation front (see also Paper I for more details). Since the water-ice line lies at the outer edge of the inner disc (see gas surface density in the upper left panel of Figure 3), water vapour diffuses outwards due to the strong pressure gradient caused by the photoevaporative gap. It can then recondense at the water evaporation front and form water-ice pebbles. These then drift inwards and evaporate again. This process enables water to stay in the inner disc on much longer timescales than carbon-bearing species that have their evaporation lines farther out and not within the inner disc (for comparison see the position of the CO₂ ice line in the gas surface density plot in the upper left panel of Figure 3). Carbon-bearing molecules instead are either accreted onto the central star or removed via outward diffusion into the photoevaporative gap. As a result, the C/O ratio decreases over time as oxygen remains very abundant in the inner disc but carbon continues to be removed.

M_★ = 0.3 M_⊙: in this case, the C/O ratio increases after the initial rapid drop and continues to do so for the rest of the time evolution until it reaches a value of 0.5. Whereas the explanation for the initial drop is the same as for the 0.5 M_⊙ star, the following increase differs in its behaviour and explanation. Instead of the water equilibrium cycle, the 0.3 M_⊙ star has a CO₂ equilibrium cycle since the CO₂ evaporation front instead of the water-ice line now lies at the outer edge of the inner disc (see upper middle panel of Figure 3), where a strong pressure gradient leads to the outward diffusion of material. Due to its position closer to the photoevaporative gap, the CO₂ equilibrium cycle is dominating over the equilibrium cycle for water, which does not exist in this case. Therefore, the CO₂ is kept in the inner disc, leading to a continuous increase in the C/O ratio in this region. With time, all other molecules except CO₂ are accreted onto the host star, leading to a final value of 0.5 for the C/O ratio. This value is a direct consequence of the distribution of atoms in the CO₂ molecule.

It should be noted that the CO₂ equilibrium cycle, although essential in keeping the carbon in the inner disc over longer time periods, is very sensitive to the relative position of the inner edge of the photoevaporative gap and the evaporation lines. Changing, for example, the disc’s viscosity, as done in Appendix C.2, shifts the CO₂ evaporation line for the M_★ = 0.3 M_⊙ star into the gap, leading to the non-existence of the CO₂ equilibrium cycle. The gap position itself is, in turn, dependent on the photoevaporation model, as seen in Figure 1. Another important factor is the disc’s thermal structure, which depends on viscous heating and the evolution of the stellar luminosity, and which also affects the positions of the evaporation lines.

M_★ = 0.1 M_⊙: In this case, the C/O ratio increases to supersolar values after the initial drop, where the initial drop is again a result of the inner disc first being enriched with water vapour. The following increase is a consequence of carbon-rich vapour being carried into the inner disc. This happens faster in comparison to the other stellar masses because the evaporation lines for the 0.1 M_⊙ star are closer to the host star. Additionally, between 2 Myr and 3 Myr, the existence of a CO₂ equilibrium cycle ensures that carbon-rich gas stays in the inner disc. At the same time, carbon-rich gas from the outer disc can no longer reach the inner disc because the gap is already fully opened. After 3 Myr, the CO₂ equilibrium cycle ceases to exist. The reason is the inward movement of the inner edge of the gap, shifting the CO₂ ice line from the inner disc to inside the widened gap. In turn, the water-ice line now lies at the outer edge of the inner disc and a water equilibrium cycle occurs as in the case of the 0.5 M_⊙ star. As a result, the C/O ratio decreases rapidly after 3 Myr.

3.1.3 C/H, O/H, and N/H

In addition to the C/O ratio, we similarly investigate the C/H, O/H, and N/H abundances in the inner disc; see the top, middle, and bottom rows of Figure 4. As for the C/O ratio, we show the time evolution of these element ratios as a function of disc radius for different stellar masses, indicated at the top of each column.

As expected, there is not much carbon and nitrogen in the case of the 0.5 M_⊙ and 0.1 M_⊙ stars, but a lot of oxygen in their inner discs at the end of the time evolution. However, the 0.1 M_⊙ star initially has a CO₂ and an NH₃ equilibrium cycle, resulting in an elevated C/H and N/H ratio after 3 Myr compared to the C/H and N/H ratios at the same time for the 0.5 M_⊙ star. These equilibrium cycles, however, stop once the inner edge of the pho-toevaporative gap reaches the position of their ice lines, resulting in a decrease in the C/H and N/H ratio. For the 0.3 M_⊙ star, the carbon and nitrogen fractions are generally much higher, which is in agreement with our results for the C/O ratio and the position of the CO₂ and the NH₃ ice line. In contrast to the 0.5 M_⊙ and 0.1 M_⊙ stars, these two ice lines in the case of the 0.3 M_⊙ star keep their position in the inner disc throughout the full time evolution. Due to their position at the outer edge of the inner disc and the consequential existence of equilibrium cycles for both molecules, we end up with an elevated carbon and nitrogen content in the inner disc.

3.2 Reduced photoevaporation rate

3.2.1 Gas surface density and C/O ratio

The results shown in this subsection are obtained using the standard simulation parameters as described in Section 2.4 and shown in Table 3. For the photoevaporative mass loss rates, we do not use the nominal values, as in the previous subsection, but reduced rates as given in Table 2, with a reduction factor of 3-4. A description of how these are obtained can be found in Section 2.2.

Figure 5 shows the gas surface density in the top row and the gaseous C/O ratio in the bottom row, both as a function of disc radius, for discs with less effective internal photoevaporation. As in Figure 3, time evolution is indicated via colour coding and the different host star masses are given at the top of each column.

The gas surface density generally behaves similarly to that shown in the previous subsection, where nominal values are used for the photoevaporative mass loss. We see a decrease in the inner disc due to accretion onto the central star and viscous spreading in the outer disc regions over time. However, the gap opened by photoevaporation occurs much later, at around 56 Myr. Again, the 0.1 M_⊙ star is the one with the highest mass loss in the inner disc.

The delayed gap opening directly impacts the C/O ratio; see the lower panels of Figure 5. As before, the initial phase is dominated by water-rich pebbles reaching the inner disc first, resulting in an excess of water vapour and therefore oxygen, which leads to a low C/O ratio. This phase is followed by carbon-rich vapour moving inwards, elevating the C/O ratio again. Since the gap opening is delayed in these simulations, the second phase is much longer than in the nominal case, resulting in much larger C/O ratios. The evolution of the C/O ratio during this phase follows that of a purely viscous disc; see Figure A.1. After the gap has fully opened, the behaviour of the C/O ratio switches from following that of a purely viscous disc to that described in the previous subsection. In the case of the 0.5 M_⊙ and 0.1 M_⊙ stars, the inner disc is dominated by the water equilibrium cycle. As a result, their C/O ratios drop again quickly after gap opening. However, the drop of the C/O ratio is not as extreme as for the case of higher photoevaporation rates, compare the bottom row of Figure 5 for the reduced photoevaporation rates with that of Figure 3 for the nominal values. The reason for this is that, in the case of reduced photoevaporative mass loss rates, the disc has more time to evolve without gap opening, resulting in a decay of the water content of the inner disc due to accretion onto the host star. Consequently, the effect of the water equilibrium cycle is not as strong as for our nominal simulations; see Figure 3. For the 0.3 M_⊙ star, on the other hand, the inner disc is dominated by a CO₂ equilibrium cycle after the gap is fully opened. The C/O ratio in this case decreases very slowly, approaching a value of 0.5.

Simulations of discs around low-mass stars that follow a pure viscous evolution (see Mah et al. 2023) and observations of such systems (see e.g. Kanwar et al. 2024a) indicate high C/O ratios of C/O >1 for the inner disc, which is not in agreement with the low C/O ratios we obtain from our simulations using the nominal photoevaporation rates from Section 3.1. However, already the reduction of the photoevaporative mass loss rates by a factor of 3-4, as studied in this subsection, results in C/O ratios that align much better with the observational values. Kanwar et al. (2024a) find a C/O ratio >1 for the system Sz 28, which has a host star mass of M_★ = 0.12 M_⊙. This result is in very good agreement with the C/O ratio we obtain for the inner disc of the 0.1 M_⊙ star between 2-6 Myr of its evolution; see the lower right panel of Figure 5. From this, we conclude that a photoevaporation rate reduced by a factor of 3-4 can already solve the discrepancy between models and observations of discs around low-mass stars.

Fig. 4 Different element ratios in the gas phase as a function of disc radius and time for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from Table 1. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: carbon over hydrogen. Middle: oxygen over hydrogen. Bottom: nitrogen over hydrogen. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

Fig. 5 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, but with mass loss rates reduced by a factor of 3-4. For the exact values, see Table 2. The reduced rates result in later gap opening and therefore longer disc lifetimes. Host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

Fig. 6 Gaseous C/O ratios as a function of disc radius and time for discs with internal photoevaporation due to X-rays and a host star mass of 0.5 M⊙. This plot shows a comparison between different α values and different photoevaporation rates, with the left column depicting α = 10−4, the right column depicting α = 103, the top row showing our nominal photoevaporation rate of Ṁw = 1.90460 · 108 M⊙yr−1 (see also Table 1), and the bottom row showing a reduced photoevaporation rate of Ṁw = 0.51324 · 10−8 M⊙ yr−1 (see also Table 2). Colour coding, plotting, and remaining simulation parameters are indicated as in Figure 3.

3.2.2 Influence of photoevaporation rate versus disc viscosity on the C/O ratio

In Figure 6, we plot exemplarily the C/O ratio for the 0.5 M_⊙ star for different scenarios to study the effects of the different photoevaporation strengths in comparison with the influence of the turbulence parameter α or the disc’s viscosity, respectively. The left column shows the C/O ratio for a disc with a viscous parameter of α = 10⁻⁴, which is the nominal value, and the right column depicts the same disc but with an α parameter of α = 10⁻³. Nominal values of the photoevaporation rate are used for the top row, whereas the bottom row is simulated with the aforementioned reduced photoevaporation rates.

The comparison shows that higher disc viscosities and higher photoevaporation rates lead to lower C/O ratios in the later evolution stages of the inner disc. For higher disc viscosities, the low C/O ratio is a result of the evaporation lines being shifted outwards in the disc due to enhanced viscous heating. Consequently, no water equilibrium cycle occurs as the waterice line is now lying in the gap instead of in the inner disc; see Figure C.2 for more details on the photoevaporative disc with a higher viscosity. Additionally, the inner disc is accreted much faster onto the host star as compared to the low-viscosity case. Therefore, the C/O ratio remains low after gap opening. For higher photoevaporation rates, the C/O ratio is low due to a different reason. In this case, photoevaporative gaps open fast, leaving no time for carbon-rich vapour to reach the inner disc to enhance the C/O ratio. A reduced photoevaporative mass loss rate, in combination with a viscous parameter of α = 10⁻⁴, and the subsequent delayed opening of gaps is therefore the only case where an elevated C/O ratio in the later disc evolution stages is possible; see the lower left panel of Figure 6.

Fig. 7 Gap opening time versus photoevaporation rate for the low-mass stars studied in this paper as well as the solar-mass star from Paper I (colour coded). Plotted are data points from our simulations, both from nominal (rightmost point of each line) and reduced photoevaporative mass loss rates with reduction factors of 2, 3-4, and 5, and a reciprocal fit through those points. The fit parameters for each stellar mass are listed in Table 4.

3.3 Disc lifetimes

For a more sophisticated analysis and comparison of the gap opening times depending on the photoevaporative mass loss rates, we added simulations to our dataset with general reduction factors of 2, 5, and 10 compared to the nominal values. In Figure 7, we show these together with the nominal values from Section 3.1 and the individual reduced mass loss rates from Section 3.2 for our stellar sample of low-mass stars studied in this paper as well as the solar-mass star from Paper I. However, only four data points instead of the expected five are visible in Figure 7 for each stellar mass because a reduction of the photoevapora-tive mass loss rates by a factor of 10 results in discs that do not show any signs of a photoevaporative gap. Therefore, they are not included in this figure.

We define the gap opening time as the time when the gas surface density drops below a specific threshold of Σ_thresh = 10⁻⁸ g/cm² in the area between 0.1 AU and the position of the CO ice line of the respective disc. Figure 7 depicts a comparison between the different host star masses, which are indicated via colour coding. The gap opening times are plotted as a function of the photoevaporation rate, and disc lifetimes can be inferred from the gap opening times. Our simulation results are shown as data points in the figure. Additionally, we plot a reciprocal fit function through the data points for each stellar mass. The corresponding fit parameters can be found in Table 4.

In Figure 2, we added the nominal photoevaporation rates and the ones with reduction factors of 2, 3-4, 5, and 10 to the observational data, after converting these rates to X-ray luminosities. The gap opening times for each case are highlighted with colour.

It becomes evident that the photoevaporation rate is a critical parameter as the gap opening time varies between around 1-8 Myr when changing said photoevaporation rate, depending on the host star mass. The gap opening time then directly influences the chemical evolution of the inner disc, resulting in different C/O ratios, with a later gap opening resulting in higher C/O ratios of the inner disc.

4 Discussion

4.1 Dependence on model assumptions

As seen in Sections 3.1 and 3.2, the results of our simulations strongly depend on the photoevaporative mass loss rate. However, this is not the only parameter that influences the outcome of our simulations. In this section, we qualitatively discuss the dependence of some of our model assumptions. For more details, we refer to the discussion in Paper I.

4.1.1 Photoevaporation rate

The comparison of our results in Section 3.1, where we use the nominal values from Picogna et al. (2021) for the photoevaporative mass loss, to our results in Section 3.2, where the photoevaporative mass loss is reduced by a factor of 3-4, already demonstrates how sensitive the inner disc’s chemistry is to this parameter. The reduction of the photoevaporation rate is motivated by a large spread in the X-ray luminosity of stars, as indicated by observations; see Figure 2. A result of the lower photoevaporation rates is a delay in the gap opening time, leading to higher C/O ratios in the inner disc. This is essential because the nominal photoevaporation rates overestimate the mass loss, resulting in lower C/O ratios that are not in agreement with observations of discs around low-mass stars (see e.g. Tabone et al. 2023; Kanwar et al. 2024a).

In addition to the indication by observations, simulations by Sellek et al. (2024) show that lower photoevaporation rates should be favoured over those by, for example, Picogna et al. (2021). Lower photoevaporative mass loss rates in the simulations of Sellek et al. (2024) result from additional cooling processes, which are not considered in the model of Picogna et al. (2021).

The sensitivity of the inner disc’s chemistry to the photo-evaporative mass loss rates becomes even more evident when examining the simulation results of a disc with a photoevaporation rate reduced by a factor of 10; see Appendix C. We observe that a reduction factor of 10 already reduces the photoevaporation rates so much that we do not observe a significant variation in the disc structure during our studied evolution of 10 Myr. The discs in these simulations do not show any signs of gap opening, with the results therefore looking almost identical to those of purely viscously evolving discs, compare Figures C.1 to A.1.

4.1.2 Photoevaporation of different molecules

Internal photoevaporation is implemented in our code chemcomp in a way that it acts uniformly on all molecules present in the code. As a result, all species are evaporated equally and leave the disc via photoevaporative winds on the same timescales. However, in reality, it is physically more accurate that lighter molecules are removed more efficiently than molecules with a larger mass, for example, it is easier for a photoevaporative wind to blow away H₂ rather than H₂O. This results in different disc leaving timescales for different molecules. However, it is very difficult to assess how this change in timescales would affect the evolution of the inner disc chemistry. To evaluate this, a detailed calculation of the disc leaving timescales for the different molecules present in chemcomp would be needed.

4.1.3 Refractory carbon grains

Refractory carbon is not included in the molecule list used for the simulations done in this paper. Its inclusion might, however, impact the C/O ratio of the inner disc. One possibility to take refractory carbon into account is the existence of a so-called soot line (~ 300 K), where carbon grains then sublimate at much lower temperatures than the actual sublimation temperature of carbon (2000 K). This has been proposed by Van ’T Hoff et al. (2020) in their work. However, the soot line requires heavily modified carbonaceous material normally occurring in meteoritic parent bodies and does not reflect the state of carbonaceous grains in discs. A second possibility is the destruction of carbon grains via pyrolysis and oxidation, releasing carbon into the gas phase at lower temperatures as well (Lee et al. 2010; Gail & Trieloff 2017; Wei et al. 2019). The third possibility is the reaction of carbon grains with hydrogen to form small molecules such as CH₄ or C₂H₂, which would also occur at lower temperatures (Nakano et al. 2003; Li et al. 2021; Kanwar et al. 2024b; Raul et al. 2025). Lenzuni et al. (1995) and Borderies et al. (2025) have shown that the latter is the most efficient of these processes. In the first two cases of a soot line and carbon destruction via pyrolysis and oxidation, the freed carbon will also react with hydrogen. The formation of molecules such as CH₄ or C₂H₂ might hinder the appearance of an equilibrium cycle for water or CO₂ . On the other hand, including refractory carbon via one of the abovementioned processes might lead to an elevation of the C/O ratio of the inner disc as more carbon is then available via CH₄ and C₂H₂ (Houge et al. 2025).

4.1.4 Initial disc radius

Banzatti et al. (2023) indicate that the water abundance in the inner disc depends on the outer radius of the disc, which is defined by the location of the outermost gap. For this study, they analyse the spectra of four discs observed with the JWST Mid-Infrared Instrument (MIRI): two are compact discs (1020 AU), and two are large discs (100-150 AU) with multiple gaps. The two compact discs show a water excess in comparison to the extended discs.

On the other hand, different JWST observations reveal a large water reservoir in an extended disc with substructures (Gasman et al. 2023), which may be an indication that gap structures play a more dominant role in creating water reservoirs than the disc radius (for a discussion of an extended set of sources, see Gasman et al. 2025). As discussed in Paper I for solar-mass stars, our simulations support the latter claim. This is also true for low-mass stars. On the one hand, we find in the simulations done in this paper that the disc radius does not affect the evolution of the gas surface density and the C/O ratio of lower-mass stars; see Appendix A and compare Figure A.1 with Figure A.2. On the other hand, pebbles in the outer areas of extended discs generally grow more slowly and reach the inner disc later than those in discs with smaller radii. The inner disc of an extended disc is therefore enriched with less water vapour over a longer period in comparison to the inner region of a compact disc. Smaller discs have in turn a higher water content over a shorter period because the pebbles grow faster in the outer disc regions due to a higher density and then have an enhanced drifting speed towards the inner disc. However, in both cases, the inner disc is enriched with water vapour, contradicting the argument that the disc radius is the sole indicator of a water-rich disc. Far more important for the chemical evolution of the inner disc is the presence of gaps and their opening timescales as shown in Sections 3.1 and 3.2.

4.2 Implications of our results

JWST observations of the MIRI Mid-Infrared Disk Survey (MINDS) collaboration of the system Sz 28 (M_★ = 0.12 M_⊙) suggest a high gaseous C/O ratio of C/O > 1 for the inner disc; see Kanwar et al. (2024a). Constraining the C/O ratio in the inner disc of PDS 70 (M_★ = 0.76 M_⊙) also predicts a super-stellar C/O ratio; see Portilla Revelo (2023). Other papers on low-mass stars, see Tabone et al. (2023) or Arabhavi et al. (2024), show results going in the same direction although these do not include detailed thermochemical calculations. Atacama Large Millime-ter/Submillimeter Array (ALMA) observations on the other hand mostly give estimates on the C/O ratio of the outer disc (>10 AU) (see e.g. Bosman et al. 2021; Le Gal et al. 2021). For most discs, there is no ALMA data available probing the gas in the inner disc (<10 AU).

These observational results imply that our nominal values for the photoevaporative mass loss are indeed overestimating the actual mass loss, leading to C/O ratios in the inner disc that are too low to match observations; see Section 3.1. The reduced rates discussed in Section 3.2 on the other hand fit the observational values much better. Especially, the supersolar C/O ratio of the 0.1 M_⊙ star between 2-6 Myr matches very well the observed values for the system Sz 28 (Kanwar et al. 2024a).

5 Summary and conclusions

In the second paper of our series, we performed 1D semi-analytical simulations of protoplanetary discs around low-mass stars, with masses ranging from 0.1 M_⊙ to 0.5 M_⊙. Our model includes pebble drift and evaporation. Additionally, internal photoevaporation due to X-rays is active. Combining both these mechanisms opens up new perspectives for understanding the composition and evolution of inner discs. We compared our results, where we used nominal values for the photoevaporative mass loss from Picogna et al. (2021), to discs with less effective photoevaporation.

Our results clearly indicate that internal photoevaporation strongly affects the evolution and chemical composition of pro-toplanetary discs. Due to the consequential gap opening, gas diffusing into such gaps is carried away by photoevaporative winds and pebbles from the outer regions are blocked in their inward motion.

For the nominal photoevaporation rates, this results in early gap opening and a low C/O ratio in the inner disc, where the latter contrasts with the observations of discs around low-mass stars. Our model can naturally explain the observed high C/O ratios of C/O > 1 when using photoevaporative mass loss rates reduced by a factor of 3-4, which lead to delayed gap opening times. Using the nominal mass loss rates, on the other hand, results in the removal of the discs after 1-2 Myr, which is in agreement with the observations of stars that have already lost their discs at such a young age (Pfalzner et al. 2022).

This dichotomy leads to the following implications for low-mass stars: Our model predicts that the inner discs of young low-mass stars (<2 Myr) should be oxygen-rich and carbon-poor, while the older discs (>2 Myr) should be carbon-rich. The latter require lower photoevaporation rates, which can originate from either the large spread in observed X-ray luminosities or from the photoevaporation model used in this study (Picogna et al. 2021). This model likely overestimates the photoevaporative mass loss at a given X-ray luminosity, leading to a mismatch between the calculated and the observed C/O ratios in discs around low-mass stars. A reduction of the photoevaporation rate by a factor of 3-4 brings the elemental abundances from our simulations into better agreement with observations.

Acknowledgements

We thank Giulia Perotti for helpful discussions and Manuel Güdel for providing the data to create Figure 2. Th.H. acknowledges the support of the ERC Origins grant number 832428. J.L.L. is a fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD).

Appendix A Additional material: Pure viscous disc

A.1 Nominal case

The results for the pure viscous disc are obtained using the standard simulation parameters as described in section 2.4 and shown in table 3. In this case, internal photoevaporation is switched off to have a comparison disc, which is dominated solely by viscous evolution.

In figure A.1, we show the gas surface density in the top row and the C/O ratio in the bottom row, both as a function of disc radius and time. The data are generated and plotted as in the case with photoevaporation; see section 3.1 for a detailed description. Variation in stellar mass is indicated at the top of each column, with the masses decreasing from 0.5 M_⊙ on the left to 0.1 M_⊙ on the right.

Our results for the gas surface density in the top row of Figure A.1 do not show a significant difference between the different host star masses. We see a decrease in the gas surface density in the inner disc for all three cases, resulting from the viscous accretion of disc material onto the central star. Additionally, the outer disc is dominated by viscous spreading. These two effects increase in their intensity with decreasing host star mass and happen on relatively long timescales as the viscous parameter α is small, α = 10⁻⁴. The results obtained for the gas surface density align with those from Paper I, where we studied solarmass stars.

The C/O ratio, depicted in the bottom row of Figure A.1, also shows similar behaviour to that of the inner discs around solarmass stars, except for the 0.1 M_⊙ star. For all three host star masses studied here, the C/O ratio drops during the first one million years, followed by an increase. In the case of the 0.5 M_⊙ and 0.3 M_⊙ stars, this rise in the C/O ratio continues until the end of the time evolution. For the 0.1 M_⊙ star, however, this rise is followed by a very slow decrease, leading to a more or less constant C/O ratio over time.

As in the case of a photoevaporative disc, see Section 3.1, the behaviour of the C/O ratio is a direct consequence of pebble drift and evaporation. Due to the water-ice line being closer to the star and the fact that pebbles drift much faster through the disc than gas, the inner disc is first enriched with water vapour. This results in high amounts of oxygen being present in the beginning, leading to a low C/O ratio. With carbon-rich vapour, which forms in the outer disc, arriving later in the inner regions, this trend is slowly reversed. As a result, the C/O ratio increases, reaching super-solar values after a few million years. However, in the case of the 0.1 M_⊙ star, the C/O ratio decreases again after around 4 Myr while remaining supersolar. This occurs because the viscous evolution is faster for smaller discs. The ice lines are much closer to the host star for the 0.1 M_⊙ star than for the 0.5 M_⊙ star. The former is therefore much faster dominated by CH₄ and CO, leading to a C/O ratio of about 2, as CH₄ and CO exist in the same fraction in our simulations. On the other hand, the 0.5 M_⊙ star is at the same time still dominated by solely CH₄, leading to a higher C/O ratio.

A.2 Smaller initial disc radius

This section compares our nominal viscous discs to discs with an initial disc radius reduced by a factor of 2. This decreases R_disc from [95, 65, 30] AU to [47.5, 32.5, 15] AU for our host star masses of [0.5, 0.3, 0.1] M_⊙. The rest of the simulation parameters are kept as before; see Section 2.4 and table 3 for more details.

figure A.2 shows the gas surface in the top row and the gaseous C/O ratio in the bottom row, both as a function of disc radius and time. The columns correspond to different host star masses, varying from 0.5 M_⊙ on the left to 0.1 M_⊙ on the right.

The general evolution of the gas surface density in discs with smaller initial disc radii is very similar to that of discs with larger initial radii, compare the top rows of Figure A.1 and Figure A.2. However, smaller initial radii lead to faster evolution, which can be seen in the accretion of inner disc material onto the central star and in the viscous spreading in the outer disc.

The behaviour of the C/O ratio in a viscous disc with a smaller initial radius, as shown in the bottom row of Figure A.2, is also very similar to that of a disc with a larger initial radius; see the bottom row of Figure A.1.

Fig. A.1 Disc evolution for a viscous disc without internal photoevaporation. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

Fig. A.2 Disc evolution for a viscous disc without internal photoevaporation. The initial disc radii are reduced by a factor of 2 compared to the nominal values given in Table 3. Their values are now [47.5, 32.5, 15] AU for our host star masses of 0.5 M⊙ (left), 0.3 M⊙ (middle) and 0.1 M⊙ (right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and remaining simulation parameters are indicated as in Figure 3.

Appendix B Additional material: Viscous disc with a planet

In addition to comparing with a viscous disc, it is important to also compare the photoevaporative disc to a viscous disc hosting a giant planet, especially since both of the latter cases are characterised by a gap dominating the disc evolution. In Figure B.1, we therefore plot the gas surface density in the top row and the C/O ratio in the bottom row for a viscous disc with a growing planet, again as a function of both disc radius and time. We vary the central star’s mass, as indicated at the top of each column.

Fig. B.1 Disc evolution for a viscous disc without internal photoevaporation, with a planet seed placed at [2.0, 1.8, 1.0] AU at 0.05 Myr, using an envelope opacity of κenv = 0.5 cm2/g. The planet has a final mass of about [400,165,10] MEarth and reaches pebble isolation mass at [0.12, 0.14, 0.27] Myr for our host star masses of [0.5, 0.3, 0.1] M⊙. The stellar mass decreases from the left to the right. Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

The results in this section are obtained using the standard simulation parameters as described in Section 2.4 and shown in table 3. Additionally, we use an envelope opacity of κ_env = 0.5 cm²/g. The value of the envelope opacity determines the duration of the envelope contraction phase of the planet, with high envelope opacities allowing for slower gas accretion compared to low values because high envelope opacities decrease the cooling rate of the envelope. Planets in our model then grow via pebble and gas accretion, migration is turned off.

For each case considered here, a planetary seed is placed in the disc between the H₂O and the CO₂ evaporation line at 0.05 Myr. This corresponds to [2.0, 1.8, 1.0] AU for our host star masses of [0.5, 0.3, 0.1] M_⊙. The planetary seeds do not migrate but grow over time due to pebble and gas accretion, where pebble accretion stops at [0.12, 0.14, 0.27] Myr, respectively, when the planets reach pebble isolation mass. At this time, the pebbles are blocked from reaching the inner disc and their inward flux is stopped. However, gas from the outer disc is still able to diffuse through the gap, (see e.g. Paardekooper & Mellema 2006; Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018; Weber et al. 2018), moving towards the inner disc regions. It is only slightly hindered by the growing planet. This distinguishes gaps opened by planets from those caused by internal photoevaporation. In both cases, the pebble flux to the inner disc is stopped once the gap is fully opened or the planet reaches pebble isolation mass. However, in the case of a planetary gap, gas can still move inwards, whereas photoevaporative winds carry away the gas completely, resulting in a cut-off of the inner disc from the outer disc’s gas supply.

The planets in our simulations reach final masses of about [400, 165, 10] M_Earth, for a list of their exact parameters see table B.1. Their position is seen as a dip in the gas surface density; see the upper row of Figure B.1. With time, the dip deepens as the planetary mass grows.

The C/O ratio evolves analogously to that of a pure viscous disc without a planet, compare the bottom row of Figure B.1 to the bottom row of Figure A.1. The only difference that occurs when a planet is present is a generally higher C/O ratio in the inner disc throughout the full time evolution. This is because the planets in all scenarios studied here reach pebble isolation mass relatively early, after [0.12, 0.14, 0.27] Myr for our host star masses of [0.5, 0.3, 0.1] M_⊙. After reaching pebble isolation mass, all pebbles from the outer disc are blocked. This is especially true for the water-rich pebbles since the planets are located between the water evaporation front and the ice lines of carbon-bearing molecules. As a result, the water-rich pebbles are trapped in the outer disc with no chance of evaporating their water-ice. Carbon-rich pebbles on the other hand are still able to evaporate their volatile content. Therefore, carbon-rich gas is created outside of the planet’s position, can pass the planetary gap and move to the inner disc. Subsequently, water-rich pebbles have much less time than in the purely viscous case to enrich the inner disc with water vapour before gap-opening. Consequently, the carbon-rich gas, which arrives later, needs to balance less oxygen, ultimately leading to a higher C/O ratio. In the planetary disc, it stays continuously higher than in the pure viscous disc because only water-rich pebbles are blocked but carbon-rich gas can still pass to the inner disc. The result of this process is a supersolar C/O ratio of the inner disc, which is even higher than in the pure viscous disc for the 0.5 M_⊙ and 0.3 M_⊙ stars. However, the 0.1 M_⊙ star shows the same C/O ratio regardless of whether a planet is forming in the disc or not. The reason is that the C/O ratio in this case towards the end of the time evolution at 10 Myr is entirely determined by the CH₄ and CO gas moving inwards from the outer disc regions; see the discussion in section A.1 for more details.

Finally, note that a variation in the envelope opacity κ_env only changes the final mass of the planets, but not the C/O ratio in the disc. Once a planet reaches pebble isolation mass, the pebbles are blocked by the planetary gap. However, gas can still pass through the gap and move into the inner disc. Depending on the mass of the planet, more or less gas is hindered in its flow, but this does not change the overall C/O ratio as both carbon and oxygen are hindered equally from moving inwards.

Appendix C Additional material: Photoevaporative disc

C.1 Factor 10 reduced mass loss rates

Here, we show the results for a viscous disc with active internal photoevaporation. The photoevaporative mass loss rate is reduced globally by a factor of 10 for all host star masses considered here. Our results are obtained using the standard simulation parameters as described in Section 2.4 and shown in Table 3. figure C.1 shows the gas surface density in the top row and the C/O ratio in the bottom row, both as a function of disc radius and time. The mass of the central star is varied across the different columns, indicated at the top of each of them, with the masses decreasing from 0.5 M_⊙ on the left to 0.1 M_⊙ on the right.

The results for both gas surface density and C/O ratio are consistent with those from a purely viscously evolving disc up to 10 Myr; see Appendix A. This consequently follows from the fact that such a low photoevaporation rate does not open any gaps in the disc, leaving the disc no other choice but to follow a pure viscous evolution. However, we expect a gap to open at later times, as we already see a small reduction in the gas surface density out to 100 AU in the case of the photoevaporation rates reduced by a factor of 10, when comparing the top rows of figures A.1 and C.1.

We therefore do not discuss the results for the disc with a pho-toevaporative mass loss rate reduced by a factor of 10 here and refer to Appendix A for a detailed description. The missing gaps are also the reason why the discs studied here do not appear in the plot in Section 3.3, where we discuss disc lifetimes, defined via gap opening times, as a function of photoevaporative mass loss rates; see figure 7.

C.2 Larger viscous parameter α

This section compares our nominal photoevaporative discs to discs with a turbulence parameter, α, that is larger by a factor of 10; we increase α from the nominal value of 10⁻⁴ to 10⁻³. The rest of the simulation parameters are the same as before; see Section 2.4 and Table 3 for more details.

figure C.2 shows the gas surface density in the top row and the gaseous C/O ratio in the bottom row, both as a function of disc radius and time. The columns correspond to different host star masses, varying from 0.5 M_⊙ on the left to 0.1 M_⊙ on the right. The general evolution of the gas surface density in a disc with a larger turbulence parameter, α, and a therefore higher viscosity is similar to that of a disc with a lower α, compare the top row of figure C.2 to that of figure 3. However, larger viscosities lead to smaller particles and a faster evolution, both for the accretion of inner disc material onto the central star and the viscous spreading in the outer disc.

The C/O ratio in a photoevaporative disc with α = 10⁻³, as shown in the bottom row of figure C.2, is solar (M_★ = 0.5 M_⊙) or even supersolar (M_★ = 0.3/0.1 M_⊙) before gap-opening but quickly and drastically decreases after the photoevaporative gap has opened. The amount of oxygen compared to carbon is so large that after 1 Myr, for the 0.5/0.1 M_⊙ stars, or 3 Myr, for the 0.3 M_⊙ star, no C/O ratio can be calculated within the parameter range shown here. The high C/O ratios in the first one million years are a result of the much faster evolution in the case of higher viscosities. At this point, the water-rich phase in the inner disc is already over and carbon-rich vapour has already enriched it enough for the C/O ratio to reach such high values. Additionally, high viscosities result in smaller pebbles that then drift inwards more slowly than larger pebbles. Therefore, the water enrichment of the inner disc in the case of α = 10⁻³ is not as strong as in the case of α = 10⁻⁴, resulting in less water vapour that needs to be balanced by carbon-rich gas.

Fig. C.1 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, but with a mass loss rate reduced by a factor of 10. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in figure 3.

Fig. C.2 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from table 1 and an α parameter of α = 10−3, corresponding to an increase by a factor of 10 compared to our nominal viscosity. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and remaining simulation parameters are indicated as in figure 3.

References

All Tables

All Figures

	Fig. 1 Photoevaporative gas surface density loss rate as a function of disc radius for different stellar masses, as given in Equation (5), adopted from Picogna et al. (2021). Here, the nominal values for the photoevap-orative mass loss rate, as given in Table 1, are used.
In the text

Fig. 2 X-ray luminosity, derived from the measured X-ray flux, as a function of stellar mass for all detected XEST sources, re-plotted from Güdel et al. (2007) using their original data. The distance to Taurus is estimated at 140 pc, although we note that Taurus consists of subgroups with varying distances (Galli et al. 2019). The flux of stars with more than one measurement is averaged to retrieve only one data point for the plot. The circle, star, and diamond symbols indicate the stellar class, as defined for XEST sources; and the straight line gives a linear regression for the logarithmic values, as described by the following equation, log(LX) = 1.54 log(M★) + 30.31. Synthetic values from our simulations from Sections 3.2 and 3.3 are added as purple-coloured squares, with the colour indicating the gap opening time of the disc, tgap.

In the text

Fig. 3 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from Table 1. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius, time evolution is shown in colour - from black, which corresponds to 0 Myr, to dark red, which corresponds to 10 Myr. Bottom : gaseous C/O ratio as a function of disc radius and time (colour-coded). The evaporation lines for the different molecules are given as dashed grey lines. Note that the C/O ratio is calculated from number densities and that, by definition, we have no specified C/O ratio in the gas phase beyond the CO evaporation front. We use our standard parameters for this simulation, as given in Table 3.

In the text

Fig. 4 Different element ratios in the gas phase as a function of disc radius and time for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from Table 1. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: carbon over hydrogen. Middle: oxygen over hydrogen. Bottom: nitrogen over hydrogen. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

In the text

Fig. 5 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, but with mass loss rates reduced by a factor of 3-4. For the exact values, see Table 2. The reduced rates result in later gap opening and therefore longer disc lifetimes. Host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

In the text

Fig. 6 Gaseous C/O ratios as a function of disc radius and time for discs with internal photoevaporation due to X-rays and a host star mass of 0.5 M⊙. This plot shows a comparison between different α values and different photoevaporation rates, with the left column depicting α = 10−4, the right column depicting α = 103, the top row showing our nominal photoevaporation rate of Ṁw = 1.90460 · 108 M⊙yr−1 (see also Table 1), and the bottom row showing a reduced photoevaporation rate of Ṁw = 0.51324 · 10−8 M⊙ yr−1 (see also Table 2). Colour coding, plotting, and remaining simulation parameters are indicated as in Figure 3.

In the text

Fig. 7 Gap opening time versus photoevaporation rate for the low-mass stars studied in this paper as well as the solar-mass star from Paper I (colour coded). Plotted are data points from our simulations, both from nominal (rightmost point of each line) and reduced photoevaporative mass loss rates with reduction factors of 2, 3-4, and 5, and a reciprocal fit through those points. The fit parameters for each stellar mass are listed in Table 4.

In the text

	Fig. A.1 Disc evolution for a viscous disc without internal photoevaporation. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.
In the text

Fig. A.2 Disc evolution for a viscous disc without internal photoevaporation. The initial disc radii are reduced by a factor of 2 compared to the nominal values given in Table 3. Their values are now [47.5, 32.5, 15] AU for our host star masses of 0.5 M⊙ (left), 0.3 M⊙ (middle) and 0.1 M⊙ (right). Top: gas surface density as a function of disc radius and time. Bottom: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and remaining simulation parameters are indicated as in Figure 3.

In the text

Fig. B.1 Disc evolution for a viscous disc without internal photoevaporation, with a planet seed placed at [2.0, 1.8, 1.0] AU at 0.05 Myr, using an envelope opacity of κenv = 0.5 cm2/g. The planet has a final mass of about [400,165,10] MEarth and reaches pebble isolation mass at [0.12, 0.14, 0.27] Myr for our host star masses of [0.5, 0.3, 0.1] M⊙. The stellar mass decreases from the left to the right. Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in Figure 3.

In the text

Fig. C.1 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, but with a mass loss rate reduced by a factor of 10. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and simulation parameters are indicated as in figure 3.

In the text

Fig. C.2 Disc evolution for a viscous disc with internal photoevaporation due to X-rays, using the nominal photoevaporative mass loss rates from table 1 and an α parameter of α = 10−3, corresponding to an increase by a factor of 10 compared to our nominal viscosity. The host star masses vary from 0.5 M⊙ (on the left) to 0.1 M⊙ (on the right). Top: Gas surface density as a function of disc radius and time. Bottom: Gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting, and remaining simulation parameters are indicated as in figure 3.

In the text