© The Authors 2024

1 Introduction

A large number of exoplanets have now been characterised, showing very diverse planetary system architectures and exoplanet properties. However, their exact formation pathway and composition are not fully understood (for a review, see Drążkowska et al. 2023). It is therefore important to study their formation environment via simulations in order to gain greater insight into the structure and composition of protoplanetary discs and their time evolution. This in turn helps us to constrain the properties and formation scenarios of planets.

In the core accretion scenario, the planetary core forms via solid accretion. This solid accretion can come in the form of planetesimals (Ida & Lin 2004; Guilera et al. 2014; Mordasini 2018; Miguel et al. 2019; Emsenhuber et al. 2021) or pebbles (Bitsch et al. 2015, 2019; Ndugu et al. 2018; Liu et al. 2019; Izidoro et al. 2021), where pebble accretion is especially more efficient than planetesimal accretion in the outer disc regions (e.g. Johansen & Bitsch 2019). Planetesimal accretion models have traditionally been favoured, but models via pebble accretion have recently become popular. Most of these population synthesis studies tried to reproduce only observed planetary masses and orbital distances. However, it is now possible to constrain planet formation models even further with newly measured atmospheric compositions of observed exoplanets (e.g. Mollière et al. 2020; Line et al. 2021; Pelletier et al. 2021; Mollière et al. 2022; August et al. 2023). Recent theoretical studies from Mordasini et al. (2016), Schneider & Bitsch (2021a,b) or Bitsch et al. (2022) tried to additionally incorporate these atmospheric abundances.

Planetary compositions originate from the composition of the disc they form in. Therefore, it is very important to simulate the full evolution of the structure and chemical composition of protoplanetary discs. The latter is determined by the disc temperature and density, radiation fields, and the position of evaporation lines of different species (see e.g. Öberg et al. 2011; Henning & Semenov 2013; Schneider & Bitsch 2021a; Eistrup & Henning 2022; Mollière et al. 2022).

Even if the temperature of the disc were to remain constant over time, the chemical composition of the disc would still evolve. Pebbles drifting through the disc evaporate their volatile content at the evaporation lines, therefore altering its composition (e.g. Piso et al. 2015; Booth et al. 2017; Schneider & Bitsch 2021a; Kalyaan et al. 2023). On the other hand, chemical reactions play only a minor role in the chemical evolution of pebbles in discs because their drift timescales are much shorter than those of the chemical reactions (Booth & Ilee 2019; Eistrup & Henning 2022).

The composition and evolution of the disc are influenced by pressure bumps that can block inward drifting pebbles (e.g. Pinilla et al. 2012). Such traps can be caused by planets reaching their pebble isolation mass and creating gaps in the discs (e.g. Paardekooper & Mellema 2006; Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018). These gaps influence the chemical composition of inner discs (see e.g. Bitsch et al. 2021; Schneider & Bitsch 2021a,b; Kalyaan et al. 2023). Banzatti et al. (2023), Grant et al. (2023), Perotti et al. (2023), and Tabone et al. (2023) studied the chemical composition of discs via James Webb Space Telescope (JWST) observations, showing that pressure bumps can indeed influence the chemical composition of inner discs.

However, disc evolution is not only determined by viscous effects. Evaporation via magnetohydrodynamic (MHD) winds (for a review, see Lesur et al. 2023) or photoevaporation (for a review, see Pascucci et al. 2023) plays a significant role as well because these processes carry away disc material, which, in the case of photoevaporation, carves deep gaps into the disc structure. Thus, these effects not only influence disc evolution and their late stages of dispersal but also change the disc’s chemical composition. Photoevaporation becomes especially important as it determines the final stages of disc evolution (e.g. Ercolano et al. 2009; Pascucci & Sterzik 2009; Ercolano & Clarke 2010; Owen et al. 2012, 2013). Such a disc dissolving process is necessary since it is known from observations that discs only live fora few million years (Mamajek et al. 2009; Fedele et al. 2010), even though the disc lifetime can vary significantly depending on the stellar type of the system (e.g. Michel et al. 2021; Pfalzner et al. 2022).

Photoevaporation is a process in which high-energy radiation transfers so much of its energy to gas particles in the disc that their velocities become high enough to escape. Energetic radiation can come from nearby stars outside the system, in which case the evaporation process is referred to as external photoevaporation (for a review, see Winter & Haworth 2022). Internal photoevaporation on the other hand is where high-energy radiation from the host star is used to energise particles in the disc (for a review, see Pascucci et al. 2023). Gaps created by photoevaporation block inward-moving pebbles. At the same time, gas is carried away by photoevaporative winds, thus cutting off the inner regions of the disc from their supply from outer areas.

In this work, we combine a viscous disc evolution model – that includes pebble drift and evaporation – with a model of internal photoevaporation. We study protoplanetary discs around solar-like stars (M_* = 1 M_⊙) and compare our results first to a purely viscous disc and second to a disc with a forming giant planet to distinguish between the two kinds of gaps that form for each case.

This paper is structured as follows: Section 2 includes an outline of our model and an explanation of the numerical methods used in chemcomp. The obtained results are shown in Section 3, and their implications are discussed in Section 4 before we conclude with Section 5.

2 Methods

For our numerical simulations, we use the chemcomp code as described in Schneider & Bitsch (2021a). The code is a 1D, semi-analytical model of protoplanetary discs coupled to a planetary growth model, where the latter is only used for one case where we compare the effects of a planet on the disc structure. The details of the planet formation model can be found in Schneider & Bitsch (2021a). The disc evolution follows a classic viscous evolution model (see e.g. Lynden-Bell & Pringle 1974), making use of the alpha-viscosity description from Shakura & Sunyaev (1973).

In the disc, dust grains can grow to pebbles, with their growth being modelled by the two-population approach developed by Birnstiel et al. (2012), where the grains exist in two representative sizes, the monomer size and the size of a large grain population. Their growth is being limited by drift, fragmentation and drift-induced fragmentation. The disc is enriched by volatiles evaporating from the inward drifting pebbles at the particular ice lines of each species (Schneider & Bitsch 2021a).

2.1 Viscous evolution

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

where α is a dimensionless factor describing the strength of the turbulence, c_s is the isothermal sound speed and ${\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{\frac{G{M}_{\star }}{r^3}}$ is the Keplerian angular frequency. The speed of sound can be linked to the mid-plane temperature of the disc. The temperature is calculated from an equilibrium between viscous and stellar heating with radiative cooling. In our model, it does not evolve in time for simplicity and is thus calculated only at the initialisation of the simulation.

The time evolution of the gas surface density in our disc is given by the viscous disc equation, which can be derived from the conservation of mass and angular momentum (Pringle 1981; Armitage 2013), $$\frac{\mathrm{\partial }{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{a}\mathrm{s},\mathrm{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 }}_{\mathrm{g}\mathrm{a}\mathrm{s},\mathrm{Y}}\right)\right]={\dot{\mathrm{\Sigma }}}_{\mathrm{Y}},$$(2)

where ${\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}$ is the source term of the molecular species Y, as given in Table 1. It originates from the evaporation and condensation of pebbles and is given by $${\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}=\{\begin{array}{cc}{\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}^{\text{evap\hspace{0.17em}}}& r<r_{\text{ice\hspace{0.17em}},\mathrm{Y}}\\ {\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}^{\text{cond\hspace{0.17em}}}& r\ge r_{\text{ice\hspace{0.17em}},\mathrm{Y}}.\end{array}$$(3)

Here, ${\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}^{\text{evap\hspace{0.17em}}}$ and ${\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}^{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}$ are the evaporation and condensation source terms of species Y, originating from the evaporation and condensation of volatiles of species Y at the respective ice line r_ice,Y.

2.2 Disc lifetime and internal photoevaporation

When a viscous disc is subject to internal photoevaporation, an additional term ${\dot{\mathrm{\Sigma }}}_{\mathrm{w}}$, describing the photoevaporative gas surface density loss rate, is subtracted from the right-hand-side of the disc’s viscous Equation (2), $$\frac{\mathrm{\partial }{\mathrm{\Sigma }}_{\text{gas\hspace{0.17em}},\mathrm{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\hspace{0.17em}},\mathrm{Y}}\right)\right]={\dot{\mathrm{\Sigma }}}_{\mathrm{Y}}-{\dot{\mathrm{\Sigma }}}_{\mathrm{w}}.$$(4)

We use a description for internal photoevaporation due to X-rays following the work of Picogna et al. (2019, 2021) and Ercolano et al. (2021). They only use the soft X-ray regime as it has shown to be the most efficient for internal photoevaporation. The surface density loss is given by $$\begin{array}{rl}{\dot{\mathrm{\Sigma }}}_{\mathrm{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}}_{\mathrm{w}}\left[M_{\odot }{\mathrm{A}\mathrm{U}}^{-2}{\mathrm{y}\mathrm{r}}^{-1}\right],\end{array}$$(5)

with $$\begin{array}{rl}\frac{{\dot{M}}_{\mathrm{W}}(R)}{{\dot{M}}_{\mathrm{W}}\left(L_{\mathrm{X}}\right)}=& {10}^{a{\log }_{10}(R{)}^6+b{\log }_{10}(R{)}^5+c{\log }_{10}(R{)}^4+d{\log }_{10}(R{)}^3}\\ & \times {10}^{e{\log }_{10}(R{)}^2+f{\log }_{10}(R)+g},\end{array}$$(6)

where the fit parameters for a solar-mass star are presented in Table 2.

The loss rate of the gas surface density is visualised in Figure 1 as a function of disc radius. It is constant in time but varies in radius, meaning that the same amount of disc material is taken away in each time step of our simulations. The peak of photoevaporative mass loss is around 5 AU, with the mass generally being lost between 1.7 AU and 119 AU, enabling the creation of a gap in this area. Disc evolution in regions unaffected by photoevaporation is then purely determined by the viscous transport of gas and dust.

Fig. 1 Gas surface density loss rate as a function of disc radius, as given in Equation (5), taken from Picogna et al. (2021).

2.3 Pebble drift and evaporation

In our model, protoplanetary discs consist of gas and dust grains. The latter can grow to millimetre- to centimetre-sized pebbles, with their growth being limited by drift, turbulent fragmentation and drift-induced fragmentation. Gas and pebbles will move through the disc, with the pebbles drifting faster due to gas drag and therefore reaching the inner disc first. On their way through the disc, the pebbles will cross several evaporation fronts of different molecular species and evaporate their respective volatile content. Enrichment of the inner disc is therefore a natural outcome of our model. Gas close to the host star will eventually be accreted onto it. It is also possible for gaseous species to recondense at such evaporation lines when crossing them in the outwards direction, which can lead to large enhancements of these particular solids compared to the solar composition (e.g. Aguichine et al. 2021; Mousis et al. 2022; Mah & Bitsch 2023).

Dust grains move through the disc with a radial velocity u_Z (Weidenschilling 1977; Brauer et al. 2008; Birnstiel et al. 2012) $$u_{\mathrm{Z}}=\frac{1}{1+{\mathrm{S}\mathrm{t}}^2}{u}_{\mathrm{g}\mathrm{a}\mathrm{s}}-\frac{2}{{\mathrm{S}\mathrm{t}}^{-1}+\mathrm{S}\mathrm{t}}\mathrm{\Delta }v,$$(7)

where St is the Stokes number of a particle, and u_gas and Δv are the velocity of the gas and its azimuthal speed, respectively. The latter is defined by $$\mathrm{\Delta }v=v_{\mathrm{K}}-v_{\phi }=-\frac{1}{2}\frac{\mathrm{d}\ln P}{\mathrm{d}\ln r}{\left(\frac{H_{\mathrm{g}\mathrm{a}\mathrm{s}}}{r}\right)}^2{v}_{\mathrm{K}},$$(8)

where v_K = Ω_Kr is the Keplerian velocity, P is the gas pressure and $H_{\text{gas\hspace{0.17em}}}=\frac{c_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{K}}}$ is the scale height of the disc.

The two-population approach by Birnstiel et al. (2012) used here divides the grain size distribution into two groups: small and large grains. This model then takes the average of their velocities. The code self-consistently stops the inward drifting pebbles at pressure bumps.

Equation (7) implies that small dust grains growing larger increase their radial velocity. On the other hand, their growth can be limited by their velocity. When exceeding a certain velocity limit, the dust grains will fragment when colliding (Birnstiel et al. 2009). This fragmentation velocity has been measured in laboratory experiments to values of 1 m/s to 10 m/s (Poppe et al. 2000; Gundlach & Blum 2014). We use here a variable fragmentation velocity that changes its value from 1 m/s to 10 m/s at the water-ice line in the radial direction. This is motivated by the idea that water-ice grains can reach larger fragmentation velocities before they fragment compared to silicate grains (Gundlach & Blum 2014, but see also Musiolik & Wurm 2019).

2.4 Chemistry model and initial conditions

We assume that the disc’s original chemical composition is similar to the one of its host star. Since we investigate solar-like stars, it is reasonable to use solar abundances [X/H] for our initial disc composition, see Table 3 for the exact values which are taken from Asplund et al. (2009). Since the disc temperature varies according to the orbital distance from the star, the disc composition is likewise dependent on the disc radius.

Elements are distributed into molecules Y following a simple chemical partitioning model from Schneider & Bitsch (2021a), based on Bitsch & Battistini (2020). Molecules from a certain species Y are then either present in gaseous or solid form, dependent on their condensation temperature and position in the disc. If they are at a location where the disc temperature is higher than their condensation temperature, the species are available as vapour, and if they are at a location where the disc temperature is lower than their condensation temperature, the species are available in condensed form. The point where the mid-plane temperature is equal to the condensation temperature of a certain molecule species is the evaporation line of that species.

The number ratio of two elemental species X₁ and X₂ is defined via $${\mathrm{X}}_1/{\mathrm{X}}_2=\frac{m_{{\mathrm{X}}_1}}{m_{{\mathrm{X}}_2}}\frac{{\mu }_{{\mathrm{X}}_2}}{{\mu }_{{\mathrm{X}}_1}},$$(9)

where $m_{{\mathrm{X}}_1}$ and $m_{{\mathrm{X}}_2}$ are the mass fractions of the two elements in a medium of mass m (e.g. m = Σ_gas) and ${\mu }_{{\mathrm{X}}_1}$ and ${\mu }_{{\mathrm{X}}_2}$ are the atomic masses of the specific elements. In this paper, this definition is used to calculate the C/O ratio.

For our standard simulations, we use a star with a mass of M_* = 1 M_⊙ and a luminosity of L_* = 1 L_⊙. Our disc is determined by a viscous α parameter of α = 10⁻⁴, an initial disc mass of M_disc = 0.1 M_⊙, an initial disc radius of R_disc = 75 AU and an initial dust-to-gas ratio of 2%, as listed in Table 4. Our initial gas surface density profile is determined by an exponential decay in the outer disc regions. The initial radius is the point where the exponential cut-off sets in.

3 Results

In this work, we discuss three different disc evolution scenarios: a purely viscously evolving disc (Scenario 1), a viscous disc with additional internal photoevaporation due to X-rays (Scenario 2), and a viscous disc with a giant planet forming in it (Scenario 3). As this work primarily focuses on internal photoevaporation and its effects on disc evolution and composition, the main part of Section 3 is devoted to Scenario 2.

In Section 3.1, we discuss the results for our standard parameter case, where we compare a purely viscously evolving disc (Scenario 1) to one with additional internal photoevaporation (Scenario 2) and a viscous disc with a forming planet (Scenario 3). Section 3.2 extends the standard parameter case to a parameter study of a disc with photoevaporation. We conclude Section 3 with a special case, where we compare our standard case of a disc with internal photoevaporation to a disc with a lower photoevaporative mass loss rate in Section 3.3.

3.1 Standard case

All results shown in this subsection are obtained using the standard simulation parameters as described in Section 2.4 and shown in Table 4.

3.1.1 Scenario 1: Pure viscous disc

In the top row of Figure 2, we show the gas surface density on the left and the gaseous C/O ratio on the right, both as a function of disc radius, for a viscous disc without internal photoevaporation. The time evolution over a span of 10 Myr is given by colour coding. The C/O ratio is calculated from number densities, normalised to solar values.

Our simulations confirm the expected behaviour for a viscously evolving disc with regard to the gas surface density. It decreases over time in the inner disc due to the viscous accretion of disc material onto the host star. Additionally, the outer parts of the disc are dominated by viscous spreading. Both of these processes happen relatively slowly as the alpha parameter is small, α = 10⁻⁴.

The initial C/O ratio shows the classical step-like behaviour (e.g. Öberg et al. 2011; Mollière et al. 2022), where the different volatiles evaporate at different distances, resulting in a change in the C/O ratio. The positions of these evaporation fronts are determined by the disc temperature. The C/O ratio either increases at the evaporation lines of carbon-rich molecules or decreases at those of oxygen-rich molecules.

As the disc starts to evolve, the C/O ratio decreases rapidly in the first one million years in the inner disc, interior to the water-ice line. With time, it increases again, reaching values even higher than the initial ones at the end of the time evolution (see also Mah et al. 2023). This is a direct consequence of pebble drift and evaporation. Since the position of the water evaporation line is at a smaller disc radius than those of carbon-bearing molecules, pebbles from the outer regions moving inwards through the disc therefore first reach the evaporation lines of CO, CO₂ and CH₄ and evaporate their carbon content before evaporating water vapour. This means that the water vapour is formed closer to the star. In addition, pebbles move much faster than gas, water-ice pebbles thus reach the inner parts of the disc much earlier than the carbon-rich gas and evaporate their water content there. Therefore, the inner disc is enriched with water vapour first. The large amount of oxygen in water leads to a low C/O ratio. With time, more and more carbon-rich gas enters the inner disc region, elevating the C/O ratio again. At the same time, the material of the inner disc is accreted onto the host star. In the end, the C/O ratio reaches supersolar values.

The high C/O ratio we see in the inner disc at the end of the time evolution is a consequence of the large methane abundance in our code, see Table 1. If methane were less abundant, the C/O ratio would not be able to increase as much as it does here and reach values larger than 1 (see also Mah et al. 2023).

Our results for a viscously evolving disc essentially reproduce studies by e.g. Piso et al. (2015), Booth et al. (2017), Bitsch et al. (2021), Kalyaan et al. (2021), Eistrup & Henning (2022) or Bitsch & Mah (2023), who have also shown that it is important to study the chemical evolution of protoplanetary discs caused by pebble flow.

3.1.2 Scenario 2: Viscous disc with additional internal photoevaporation

In a disc, where in addition to a pure viscous evolution, as in the top row of Figure 2, internal photoevaporation is active as well, the behaviour of the gas surface density and the C/O ratio changes significantly, see middle row of Figure 2.

Between 2–3 Myr, a large gap opens in the gas surface density, spanning from 5.5 AU to 117 AU. With time, its inner edge moves even further inward in the disc, which is a result of the particular shape of the gas surface density loss rate by photoevaporation (see Figure 1). The opened gap leads to a blocking of the pebbles from the outer disc regions. At the same time, inward-moving gas is carried away by photoevaporative winds. The gap itself is a direct result of internal photoevaporation, where gas is taken from the disc at the location where photoevaporation has the strongest impact. This is not necessarily at the position of the peak of the photoevaporative mass loss rate, as we see in our data. The left panel of the middle row of Figure 2 shows that the gap starts to open around 110 AU, but the peak of the photoevaporative surface density loss rate is at 5 AU, see Figure 1. The reason is that the disc’s surface density decays exponentially in its outer regions resulting in the fact that the relative impact of photoevaporation is strongest here and not at its peak.

The gap has a large impact on the chemical evolution of protoplanetary discs, as can be seen from the C/O ratio (right panel of the middle row of Figure 2). It remains subsolar in the inner disc throughout the whole evolution and never reaches supersolar values as in the case of a purely viscously evolving disc. More specifically, the C/O ratio in the case with internal photoevaporation shows similar behaviour to that of the pure viscous case in the first two million years since the gap has not opened yet. During that time, it decreases rapidly in the inner disc compared to its initial state, resulting from the evaporation of water-rich pebbles. After some time, carbon-rich gas slowly starts to enrich the inner disc as well, leading to a rise in the C/O ratio. However, before this process can lead to a supersolar C/O ratio, as in the viscous case, photoevaporation takes away too much disc material and a gap opens, as can be seen in the left panel of the middle row of Figure 2. The pebbles from outer disc regions are blocked and gas is carried away by photoevaporative winds. Therefore, no more carbon-rich gas reaches the inner disc, the C/O ratio remains low.

To understand why the C/O ratio drops again after the gap-opening, we take a closer look at the water content in the inner disc. In Figure 3, the relative water abundance with respect to the total gas surface density is shown as a function of disc radius. We see that after an initially low abundance, it increases in the first one million years when the water vapour, evaporated from pebbles, arrives. With the arrival of carbon vapour, the water abundance decreases for the next two million years. After that, the gap is fully opened and the flux of carbon vapour is blocked, leading to an increase in the relative water abundance. The reason for this is that the gas already present in the inner disc is accreted onto the host star, with one exception. Since the water evaporation front lies within the inner disc, a water equilibrium cycle is present. Water vapour can recondense at the water-ice line and form water-ice pebbles, which in turn drift inwards and evaporate again. This enables water to stay in the inner disc for much longer than other species like CO, CO₂, and CH₄ that have their evaporation lines further out and not within the inner disc. These species are then either accreted onto the central star or removed via outward diffusion into the photoevaporation zone. The water equilibrium cycle is so effective that the relative water abundance in the inner disc increases until it makes up more than 90% of the inner disc mass at the end of its lifetime.

This equilibrium cycle is also the reason for the decrease in C/O in discs with internal photoevaporation after three million years, as seen in the right panel of the middle row of Figure 2. The carbon-bearing molecules are accreted onto the host star, while water, and therefore oxygen, remains present in the inner disc for much longer due to the above-mentioned equilibrium cycle. Consequently, the C/O ratio decreases over time again. This is the exact opposite of what happens in the pure viscous case.

Additionally, we investigated the time evolution of the C/H, O/H and N/H abundances in the inner disc for the case with internal photoevaporation, see Figure 4. As expected, there is not much carbon and nitrogen but a lot of oxygen in the inner disc at the end of the disc’s lifetime. This is in agreement with our study of the water content in the inner disc. For the calculation of the elemental abundance ratios, no distinction is made as to which molecule the elements belong to. Therefore, C/H and N/H decrease rapidly at the water-ice line because of the hydrogen excess due to water evaporation. Additionally, the hydrogen that is bound in water molecules remains in the inner disc due to the water equilibrium cycle, in contrast to carbon and nitrogen which are much more volatile and are being accreted onto the host star.

Fig. 2 Disc evolution for different disc configurations. Top: viscous disc without internal photoevaporation. Middle: viscous disc with internal photoevaporation due to X-rays. Bottom: viscous disc without internal photoevaporation where a planet seed is put at 3 AU at 0.05 Myr. The planet has a final mass of about 2860 MEarth, and reaches pebble isolation mass at 0.1 Myr. Left: 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. Right: gaseous C/O ratio as a function of disc radius and time (colour-coded). The evaporation lines for different molecules are given as grey dashed lines. Note that the C/O ratio is calculated from number densities scaled to the solar value and that we, by definition, 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 4.

Fig. 3 Relative water content as a function of disc radius and time. Internal photoevaporation is active in this simulation. Note that the relative water content is calculated as a ratio of gas surface densities, and is therefore shown as a mass fraction here. Colour coding, plotting and simulation parameters as in Figure 2.

3.1.3 Scenario 3: Viscous disc with a planet

Photoevaporation is just one possible mechanism causing gaps in protoplanetary discs, while another one is growing giant planets. The main difference is that photoevaporation takes away gas from the disc leading to a complete blocking of the pebble flow as soon as a gap is opened. Additionally, gas diffusing into the gap is carried away by photoevaporative winds. A planet on the other hand only blocks the pebbles from the outer parts, but not the gas (see e.g. Paardekooper & Mellema 2006; Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018; Weber et al. 2018). In this case, the gas can still move inwards through the disc, being only slightly hindered, but not blocked by the planet.

We study how the gas surface density and C/O ratio change for our standard parameter case when instead of internal photoevaporation, a planet is causing the gap (see bottom row of Figure 2). A planet seed is therefore placed in the disc at 3 AU, which is between the water and the CO₂ evaporation front in our model, at 0.05 Myr. It does not migrate but grows over time due to pebble and gas accretion, where pebble accretion stops at 0.1 Myr when the planet reaches the pebble isolation mass. At this time, also the inward flux of pebbles is stopped and the inner disc cannot be reached by water-ice pebbles anymore. In the end, the planet has a mass of about 2860 M_Earth.

The position of the planet is seen in the gas surface density in the left panel of the bottom row in Figure 2. At its position, the gas surface density shows a dip, deepening with time as the planetary mass increases.

The C/O ratio evolves analogously to the case of a viscous disc without a planet, with the only difference that here, it is generally higher in the inner disc over the complete time evolution. The reason is that the planet reaches pebble isolation mass already after 0.1 Myr, blocking all pebbles from the outer disc at this time. Since the planet is located between the water evaporation front and the ice lines of carbon-bearing molecules, waterrich pebbles have much less time than in the purely viscous disc to enrich the inner disc with water vapour before gap-opening. Therefore, the later arriving carbon-rich gas needs to balance less oxygen, leading to a higher C/O ratio. It stays continuously higher than in the pure viscous case because only water-rich pebbles are blocked. Carbon-rich pebbles have their evaporation lines further out than the planet’s position. The carbon-rich gas passes the planet and enriches the inner disc. The result is a supersolar C/O ratio in the inner disc at the end of the time evolution, which is even higher than in the pure viscous case.

For a side-by-side comparison between the three cases discussed in this section, a viscous disc, a disc with internal photoevaporation and a disc with a forming planet, see Appendix A. We show there the gas surface density and the C/O ratio at 5 Myr.

3.2 Parameter study

The results of our standard case in the previous section are embedded in a much larger parameter study. In this section, we investigate a range of alpha parameters, initial disc masses and radii. We use α = {10⁻⁴, 10⁻³, 10⁻²}, M_disc = {0.01 M_⊙, 0.1 M_⊙} and R_disc = {75 AU, 150 AU), the stellar mass in all simulations remains M_* = 1 M_⊙. The results of the parameter study for the pure viscous case are shown in Appendix B, the ones for the viscous case with internal photoevaporation are shown here.

In Figure 5, we plot the gas surface density as a function of disc radius for all the parameters mentioned above. α is increased along the x-axis of the grid. Each row corresponds to a specific combination of initial disc mass and radius, as depicted in the first panel of each row, respectively. Our standard case from Section 3.1.2 is shown in the third panel of the first column.

We see that the influence of α is strongest in the inner disc and the very outer parts. The larger α is, the faster the accretion of inner disc material onto the host star. For the largest value of α = 10⁻², the inner disc has in all cases of initial disc mass and radius completely vanished after 10 Myr. Additionally, a stronger viscous spreading in the outer disc is noticed for larger α values.

An impact of the initial disc mass on the simulations is seen in the gap-opening time. For larger initial disc masses of M_disc = 0.1 M_⊙, the gap caused by internal photoevaporation opens 1–2 Myr later than in the M_disc = 0.01 M_⊙ case. The reason is that for larger initial disc masses, the gas surface densities are higher and can thus withstand the photoevaporative mass loss for a longer time.

The effect of the initial disc radius on our simulations seems rather small because it changes the initial gas surface density and the corresponding accretion rate just by a factor of 2. However, comparing the third and the fourth panel in the first column of Figure 5, we see that the position of gap-opening due to photoevaporation changes when using a different initial disc radius. For an initial radius of 75 AU, the hole starts to open around 200 AU, see the third panel of Figure 5, whereas for an initial radius of 150 AU, the position of gap-opening lies much further inward at 15 AU, see the fourth panel of Figure 5. Since the initial disc radius is defined as the location where the exponential cut-off in the initial gas surface density sets in, the two discs compared here differ in their rate of gas surface density versus photoevaporative mass loss, especially in the outer regions. This leads to photoevaporation being dominant in very different regions of those discs and thus in different gap-opening locations.

Figure 6 shows the C/O ratio as a function of disc radius for a viscous disc with active internal photoevaporation. The parameter grid used here is the same as in Figure 5.

A similar trend as in the standard case, which is depicted in the first panel of row three, is seen for a similar parameter configuration, α = 10⁻⁴, M_disc = 0.1 M_⊙ and R_disc = 150 AU, see first panel of row four.

In general, all parameter configurations show a similar behaviour, with subsolar C/O ratios in the inner disc regions, some of them being even completely depleted of carbon and oxygen, and therefore showing no C/O ratio at all. The latter is seen only in discs with larger values of α, because in those discs, the viscous transport is so fast and dominant that the inner disc is fully accreted onto the host star, with no carbon or oxygen left to calculate a C/O ratio from.

There is only one configuration where we have a supersolar C/O ratio in the inner disc after 2 Myr (α = 10⁻³, M_disc = 0.1 M_⊙ and R_disc = 75 AU), which then vanishes during the rest of the disc lifetime as well. This trend is only seen in this specific parameter configuration because the α value is high enough to transport a lot of carbon-rich material into the inner disc before the gap opens, compared to the parameter configuration with a lower value of α, see panel three of column one in Figure 6 (α = 10⁻⁴, M_disc = 0.1 M_⊙ and R_disc = 75 AU). This leads to an elevated amount of carbon-rich material and a therefore high C/O ratio in the inner disc during the first 2 Myr. With a smaller α (α = 10⁻⁴), the material transport is just not fast enough to result in a supersolar C/O ratio. Furthermore, the initial disc radius plays an important role as well since it defines the location of exponential cut-off in the initial gas surface density. The disc with the supersolar C/O ratio (panel three of column two in Figure 6) therefore has another rate of gas surface density versus photoevaporative mass loss than the same disc but with a larger initial disc radius (α = 10⁻³, M_disc = 0.1 M_⊙ and R_disc = 150 AU), see panel four of column two in Figure 6. As stated before, these two discs then also differ in their location of gap-opening. Since the gap for the disc with R_disc = 75 AU opens further out than that of the disc with R_disc = 150 AU, carbon-rich material can still evaporate and enrich the inner disc inside of the hole. This is not possible in the case of a larger initial disc radius because the gap starts to open in the inner disc, that is, further in than the ice lines of carbon-bearing molecules, and therefore directly hindering evaporated carbon-rich vapour from reaching the inner disc. The special behaviour of having a supersolar C/O ratio during the first 2 Myr (see panel three in column two of Figure 6) of course matches the behaviour of the other panels once the gap is fully opened, reaching inner disc regions and thus hindering carbon-rich vapour from entering the inner disc.

We also note that in the low-viscosity case (α = 10⁻⁴), a larger initial disc mass leads to an elevated C/O ratio in the inner disc, compared to smaller initial disc masses. This is because higher initial disc masses come along with higher gas surface densities that can withstand photoevaporative mass loss for a longer time, which causes the gap to open later than in the case of a disc of lower mass. Consequently, carbon-rich material from the outer disc has more time to move towards the inner disc and enrich it.

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. Left: carbon over hydrogen. Middle: oxygen over hydrogen. Right: nitrogen over hydrogen. Colour coding, plotting and simulation parameters as in Figure 2.

Fig. 5 Gas surface density for a viscous disc with internal photoevaporation due to X-rays. The gas surface density is shown as a function of disc radius and time. A variety of parameter configurations is used for this study: along the x-axis of this figure, α is varied from α = 10−4 to α = 10−2, whereas on the y-axis, different sets of disc masses and radii are used, varying in the range of Mdisc = 0.01–0.1 M⊙ and Rdisc = 75–150 AU. Colour coding, plotting and stellar mass as in Figure 2.

Fig. 6 Gaseous C/O ratio for a viscous disc with internal photoevaporation due to X-rays. The C/O ratio is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.

Fig. 7 Disc evolution for a viscous disc with internal photoevaporation due to X-rays but with a smaller mass loss rate of 3.6 · 10−8 M⊙/yr instead of 3.864 · 10−8 M⊙/yr. Left: gas surface density as a function of disc radius and time. Right: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting and simulation parameters as in Figure 2.

3.3 Lower photoevaporation rate

Since the X-ray luminosity for a given stellar mass is not unambiguously defined, the photoevaporative mass loss rate due to X-rays can vary. A lower limit for solar-mass stars can be read off the lower end of the error bar in Figure 6 of Picogna et al. (2021). For the study in this section, we therefore use a mass loss rate of 3.6 · 10⁻⁸ M_⊙/yr instead of 3.864 · 10⁻⁸ M_⊙/yr for our standard simulation to study the effects of a lower mass loss rate on our previous results. The disc parameters used here are given in Table 4.

In the gas surface density in the left panel of Figure 7, we see that the gap has fully opened after 4 Myr, one million years later than in the case with a slightly higher mass loss rate. For the C/O ratio in the right panel, only minor differences are noticeable. It is slightly higher than in the case of a disc with stronger photoevaporation. This is a direct consequence of the later gap-opening, ensuring that more carbon-rich gas can move into the inner disc before the gap is fully opened.

3.4 Graphical representation of our results

In Section 3.1, we studied three different disc evolution scenarios. To sum up our results, we picture them in a graphical representation, see Figure 8. The left column depicts a purely viscously evolving disc (Scenario 1), the middle one a viscous disc with additional internal photoevaporation due to X-rays (Scenario 2), and the right one a viscous disc with a giant planet forming in it (Scenario 3).

The protoplanetary disc is always shown as a trapezoidal shape with three areas coloured differently. This should highlight the enrichment of those disc regions with different molecules, being separated by the evaporation lines of these molecules. Those ice lines are depicted as black, curved lines, their curvature being the result of a temperature gradient between the mid-plane and surface layers of the disc (Dullemond et al. 2001; Dullemond & Dominik 2004). As an example, we focus on the two evaporation lines that are most important for our study, namely the water-ice line and the combined evaporation front of all carbon-bearing volatiles. Pebbles, displayed as black circles, are distributed evenly across the full disc, showing icy layers according to their position.

All three different evolution scenarios start with the same initial condition, shown in the top panel of Figure 8. The inner disc is initially slightly enriched with water vapour, the middle part is slightly enriched with carbon vapour and the outer part is not enriched at all. Pebbles in the outer part still contain all their volatile content in ice form.

With time evolving continuously, pebbles and gas drift through the disc. When crossing an evaporation line, pebbles evaporate their respective volatile content, enriching the disc there. Gas reaching the very inner parts of the disc is accreted onto the host star. Depending on the scenario chosen, disc enrichment and pebble drift differ, see the lower panels of Figure 8.

In a purely viscously evolving disc (Scenario 1), which is depicted in the left column of Figure 8, inward drifting pebbles evaporate their volatile content whenever crossing an evaporation front. Due to the position of the ice lines and the fact that pebbles drift faster than gas, the inner disc is first enriched with water vapour (top panel). With time, carbon vapour arrives (middle panel) and slowly enriches the inner disc (lower panel), while the gas is simultaneously accreted onto the host star. This can be seen in our results in Section 3.1.1 as well.

In the case of a viscous disc with additional internal photoevaporation (Scenario 2), see middle column of Figure 8, the disc evolves similarly to the pure viscous case in the beginning (first and second panel). When photoevaporation opens a gap, gas diffusing into it is carried away by photoevaporative winds while pebbles are blocked by the gap itself (third panel). The inner disc is then enriched with water vapour again due to the water equilibrium cycle (fourth panel). All volatiles, except water, are accreted onto the host star. Water contrarily recondenses at the water-ice line and forms water-ice pebbles that evaporate again, contributing to a full cycle of evaporation and condensation that allows the water vapour to survive for a long time, see also our results in Section 3.1.2.

A planet forming between the water evaporation front and the ice lines of carbon-bearing molecules (Scenario 3), as shown in the right column of Figure 8, leads to the same chemical evolution of the disc as in the pure viscous case. When the planet is massive enough, it opens a gap but contrary to the gap caused by photoevaporation, a gas flow to the inner disc is still possible (middle panel). Only the pebbles are blocked once the planet reaches pebble isolation mass. The only differences to a purely viscously evolving disc (see left column of Figure 8) are less water enrichment in the beginning due to the planetary embryo already blocking some inflow (top panel) and a higher carbon enrichment in late stages (lower panel), see also our results in Section 3.1.3.

Fig. 8 Phases of disc evolution, with the left column depicting a pure viscous disc, the middle column a viscous disc with internal photoevaporation, and the right column a viscous disc with a growing planet. Top panel: initial state of the protoplanetary disc, with its inner part being slightly enriched with water vapour, its middle part being slightly enriched with carbon vapour and its outer regions showing no enrichment. The evaporation line of water and the combined one of all carbon-bearing molecules are marked with black, curved lines. Their curvature is the result of a temperature gradient between the mid-plane and surface layers of the disc. Pebbles are distributed evenly across the full disc, showing icy layers according to their position in it. Left column: during the pure viscous evolution, pebbles drifting inwards evaporate their volatile content whenever crossing an evaporation front. Due to the position of the ice lines and the fact that pebbles drift faster than gas, the inner disc is first enriched with water vapour (top panel). With time, carbon vapour arrives (middle panel) and slowly enriches the inner disc (lower panel), while the gas is accreted onto the host star. Middle column: with internal photoevaporation, the disc evolves similarly to the pure viscous case in the beginning (first and second panel). When photoevaporation opens a gap, gas diffusing into it is carried away by photoevaporative winds while pebbles are blocked by the gap itself (third panel). The inner disc is then enriched with water vapour again due to the water equilibrium cycle (fourth panel). All volatiles, except water, are accreted onto the host star. Water contrarily recondenses at the water-ice line and forms water pebbles that evaporate again, contributing to a full cycle of evaporation and condensation that allows the water vapour to survive for a long time. Right column: a planet forming between the water evaporation front and the ice lines of carbon-bearing molecules leads to the same chemical evolution of the disc as in the pure viscous case. When the planet is massive enough, it opens a gap, but contrary to the gap caused by photoevaporation, a gas flow to the inner disc is still possible. The only differences are less water enrichment in the beginning due to the planetary embryo already blocking some inflow (top panel) and a higher carbon enrichment in the late stages (lower panel).

4 Discussion

4.1 Dependence on model assumptions

As shown in Section 3.2, the results of our simulations strongly depend on the chosen disc parameters. In this section, we qualitatively discuss some of the parameters that we have not investigated in detail in this work and hint at future investigations.

4.1.1 Chemistry model

Our results clearly depend on the underlying chemistry model, as presented in Table 1. Varying the initial composition will result in a changed significance of individual evaporation lines (Schneider & Bitsch 2021a,b). Models that include different fractions of water and carbon- and oxygen-bearing molecules will therefore show a different water content in the inner disc as well as different C/O ratios.

For our initial elemental abundances, we used solar values. But to be able to compare our results to observations of protoplanetary discs, it is beneficial to use observed stellar abundances since stars can differ quite a lot in their C/O ratio which is used as an initial condition (see e.g. Buder et al. (2018) for elemental abundances and Bitsch & Battistini (2020) or Cabral et al. (2023) for molecular abundances). Additionally, the stellar abundance is not determined by the stellar mass but rather by the formation region of the system. This implies that to determine correctly the initial C/O ratio for individual systems, one should keep in mind that the stellar compositions are not necessarily solar, even for solar-mass stars.

Additionally, a different disc composition leads to a change in the mid-plane temperature. However, changing only the composition does not affect the latter very much. Stronger effects are only seen when choosing a larger disc metallicity (Bitsch & Savvidou 2021). Since the mid-plane temperature influences the position of the evaporation lines, a change will have a significant impact on the chemical composition of the inner disc, but only if the water-ice line is shifted so far out that it does not lie in the inner disc anymore. This implies that the temperature evolution of the disc is important for the continuous water evaporation/recondensation cycle in the inner disc and would depend on the specific system.

4.1.2 Fragmentation velocity

How fast a pebble can become before a collision leads to fragmentation instead of coagulation is determined by the fragmentation velocity. Since collision properties of dust aggregates are not easy to quantify theoretically, fragmentation velocities are measured in laboratory experiments (e.g. Poppe et al. 2000; Blum & Wurm 2008; Gundlach & Blum 2014; Musiolik & Wurm 2019). Typical results range from 1 m/s to 10 m/s.

High fragmentation velocities result in larger grains (higher Stokes numbers) and high pebble velocities. Consequently, lower fragmentation velocities lead to a decrease in the grain size and reduced pebble velocities (Brauer et al. 2008; Birnstiel et al. 2012).

In our simulations, we use a variable fragmentation velocity that changes its value from 1 m/s to 10 m/s at the water-ice line in the radial direction. This is motivated by the idea that waterice grains can reach larger fragmentation velocities before they fragment compared to silicate grains (Gundlach & Blum 2014).

Using instead a constant fragmentation velocity of 1 m/s throughout the whole disc would slow down all initial processes where material is transported radially inwards through the disc. A slower drift of the pebbles leads to less carbon-rich gas being evaporated and thus even less carbon vapour in the inner disc before the gap opens, compared to the case of a variable fragmentation velocity. This results in a smaller C/O ratio. As the same holds for the transport of water, its equilibrium cycle will most likely just start slower. The results for the relative water abundance and the C/O ratio in the inner disc will therefore not change considerably.

Increasing the overall fragmentation velocity to 10 m/s will not influence the transport from the outer to the inner disc, because the change of the fragmentation velocity just happens at the water-ice line. Therefore, the outer disc is unchanged compared to our initial assumption, and the inner disc does not change, because we only look at the gas phase (Bitsch & Mah 2023).

4.1.3 Stellar mass

In this paper, all simulations have been carried out with a stellar mass of 1 M_⊙. Changing the star’s mass will also change the structure of its surrounding disc, resulting in a different initial disc mass and radius. Additionally, the disc’s temperature profile changes, leading to a shift of the evaporation lines, as well as the photoevaporative mass loss rate due to X-ray radiation from the star. Since these factors cannot be considered separately due to their complex interplay, a detailed prediction of the water content and the C/O ratio in the inner disc for a non-solar mass star with active internal photoevaporation is nearly impossible. We will therefore conduct a similar, extensive study for M dwarfs in a follow-up paper. For studies of intermediate-mass stars (1–3 M_⊙) with active internal photoevaporation, see Ronco et al. (2024), where they investigate the effect of the latter on the disc structure without focusing on the chemical evolution.

4.1.4 Radial position of the gap

For the comparison between a gap caused by internal photoevaporation and one that is caused by a planet, we place a planetary seed at 3 AU. The planet then forms between the water and the CO₂ ice line. The resulting C/O ratio in the inner disc in this case is supersolar, see Section 3.1.3 for more details.

To study the effect of the radial position of such a planetary gap, we simulate two additional scenarios. In the first one, the planetary seed is placed at 10 AU, between the CO₂ and the CH₄ ice line. This results in a visible gap in the protoplanetary disc at the position of the planet. The C/O ratio in the inner disc is supersolar, as in our initial case. In the second scenario, the planetary seed is placed at an even larger disc radius of 50 AU, which is between the CH₄ and the CO ice line. In this case, no visible gap opens up in the disc since the planet never reaches the pebble isolation mass. The resulting gas surface density and the C/O ratio show no significant deviation from the gas surface density and the C/O ratio of a purely viscously behaving disc without a planet, where the latter is shown in the top row of Figure 2. In the end, the C/O ratio of the second scenario reaches supersolar values as well.

This simple comparison shows that the C/O ratio in the inner disc is independent of the planetary location. It is supersolar for all cases studied here. The results of this comparison of course strongly depend on the underlying chemistry model. For a further discussion of the latter, see Section 4.1.1.

4.2 Implications of our results

4.2.1 Disc observations

Our results clearly indicate that gaps caused by planets and gaps caused by internal photoevaporation differ strongly in how they affect the protoplanetary disc structure and composition. Gaps caused by planets only block the pebble flow from the outer parts once the planet reaches pebble isolation mass but not the gas. Photoevaporative gaps on the other hand block not only the pebble flow but also lead to gas diffusing into them being carried away by photoevaporative winds. These mechanisms leave no room for further enrichment of the inner disc once a gap is fully opened. This results in very different chemical compositions of the inner disc regions, where discs hosting planets can have supersolar C/O ratios, whereas discs with active internal photoevaporation have difficulties reaching supersolar C/O ratios, see Sections 3.1.2 and 3.1.3. As a consequence, we should be able to distinguish observationally between these two scenarios when measuring the C/O ratio. This implies that we can infer the root cause of a gap structure when observing protoplanetary discs. Since these gap structures are equally seen in the gas and dust phases in our simulations, observations are not limited to one or the other.

Banzatti et al. (2023) find in their study of the spectra of four discs observed with the JWST Mid-Infrared Instrument (MIRI) that the water abundance is mainly dependent on the outer disc radius. They show that the analysed spectra of the two compact discs (10–20 AU) they studied hint towards a water excess compared to the two large discs with multiple gaps (100–150 AU). Another JWST observation shows a big water reservoir in a large disc with substructures (Gasman et al. 2023). This may hint towards the fact that gap structures play a more important role in creating water reservoirs than the disc radius. In agreement with the observation by Gasman et al. (2023), we find that the disc radius does not play a major role in the time evolution of the gas surface density and the C/O ratio and therefore the water content of the inner disc, see Section 3.2. Far more important than the disc radius is the presence of gaps in protoplanetary discs and whether they are caused by planets or photoevaporation, see Sections 3.1.2 and 3.1.3. Our results point to a planetary origin of the gaps in the two extended discs studied by Banzatti et al. (2023). This stems from the fact that we see a water excess in our results only in discs where the substructures are caused by internal photoevaporation. However, the two extended discs studied by Banzatti et al. (2023) show no water excess.

Grant et al. (2023) have shown that the disc around the T Tauri star GW Lup is very abundant in CO₂, using observations from the JWST-MIRI programme: MIRI Mid-Infrared Disk Survey (MINDS). As a possible explanation, they proposed among other scenarios an inner cavity between the H₂O and CO₂ ice lines. With our model, we are now able to identify the cause of this possible gap. Using our arguments above, a high abundance in CO₂ points strongly towards a planet in this inner cavity, rather than the gap being formed by internal photoevaporation. Although GW Lup has a stellar mass of M_* = 0.46 M_⊙ (Alcalá et al. 2017; Andrews et al. 2018), the main trend of our study for solar-mass stars, namely a clear compositional difference of the inner disc depending on which scenario is the main cause for the gap-opening, will still be valid in this mass regime. The reason is that here the exact positions of the ice lines do not matter, but instead their relative locations.

Measured CO₂ abundances of Grant et al. (2023) are higher than what we find is possible to achieve in discs with active internal photoevaporation. However, GW Lup has a lower stellar mass than the stars in our simulations. Studies of Mah et al. (2023) for discs without active internal photoevaporation have shown that less massive stars achieve supersolar C/O ratios in the inner disc faster since all of the ice lines are closer to the host star. Furthermore, the work of Picogna et al. (2021) for discs with active internal photoevaporation points out that inner disc lifetimes increase for less massive stars, indicating that more carbon-rich gas can move to the inner disc. The combination of those two factors might then allow for higher C/O ratios in photoevaporative discs around low-mass stars compared to those around solar-like stars.

The lifetimes of the discs studied by Banzatti et al. (2023) and Grant et al. (2023) are shorter than the timescale in our model that is necessary to achieve a supersolar C/O ratio in the inner disc. However, our standard disc, with which we compare here, is modelled using a small value for the turbulence parameter α and has therefore a low viscosity. Changing the disc’s viscosity directly impacts the timescale for the evaporation of molecules and the subsequent enrichment of the gas with vapour. Larger viscosities result in smaller formation timescales of supersolar C/O ratios (see Figure B.2 and also Mah et al. 2023), making them again comparable to the disc lifetimes in Banzatti et al. (2023) and Grant et al. (2023).

4.2.2 Planetary compositions

From the work done in this paper, we can infer implications for the atmospheric composition of planets forming in discs with active internal photoevaporation.

If a planet forms in the inner disc, it will accrete its atmosphere from the volatile content present there. One would naively assume that such a planet in our model will have a low atmospheric C/O ratio and a high water abundance. However, we have to distinguish between different formation timescales.

On the one hand, when a planet forms early on, it might stay in the inner disc for a longer time. Due to the existence of the water equilibrium cycle and the water thus being present there for a long time as well, the planet could accrete a high amount of water, possibly forming a very water-rich planet. This is supported by observational data from PDS 70 showing that there are indeed systems with water being abundant in the inner disc (Perotti et al. 2023). However, the presence of a planet in the inner disc will most likely perturb the water equilibrium cycle. With the planet accreting part of the water vapour, less water remains to go through the equilibrium cycle, thus reducing it. In this case, the planet will most likely end up with a less water-abundant atmosphere than in the previously discussed case. On the other hand, a late-forming planet might therefore have better chances to accrete a water-rich atmosphere. The reason is that the water equilibrium cycle has not been perturbed yet, giving it more time to keep the water in the inner disc. To conclude, the composition of planetary atmospheres is strongly dependent on the point in time when the planet starts to accrete gas efficiently.

In addition to the formation timescale, the formation position and migration of a planet are additional factors that play a major role in the planet’s ability to accrete certain molecular species. This makes it overall very difficult to precisely predict the planet’s final atmospheric composition in a disc with active internal photoevaporation.

4.2.3 Debris discs

Our simulations show that a disc with active internal photoevaporation still has material left in the region outside of the gap carved by photoevaporative winds at the end of its lifetime. However, the total remaining dust mass of 1.3 M_Earth in the outer disc, when converted to planetesimals, is not enough to account for debris discs, for which several tens of Earth masses are needed.

5 Summary and conclusions

We performed 1D semi-analytical simulations of protoplanetary discs with active internal photoevaporation due to X-rays for solar-like stars. The model we used includes pebble drift and evaporation. Combining a pebble drift and evaporation model with internal photoevaporation opens up new perspectives for understanding the composition of inner discs. We compare our results to discs that evolved purely viscously without photoevaporation.

Our results clearly indicate that internal photoevaporation has significant effects on the evolution of protoplanetary discs and their chemical composition. The consequential opening of gaps strongly influences the distribution of gas and dust in such discs. Due to the diffusion of gas into these gaps – which is then carried away by photoevaporative winds – and the consecutive blocking of pebbles by the gaps themselves, it is impossible to have solar or even supersolar C/O ratios in the inner disc. Additionally, the present equilibrium water cycle allows the formation of very water-rich inner disc regions. Potential planets forming in these areas could therefore easily have large water abundances in their atmospheres.

The above-mentioned gaps due to internal photoevaporation are very different from those caused by planets forming in protoplanetary discs. The latter block only the pebble flow and not the gas, leading to a disc evolution similar to that of a pure viscous disc and resulting in a supersolar C/O ratio for the inner disc region. Future observations could therefore easily infer the root cause of a gap structure when observing the C/O ratio in the inner regions of protoplanetary discs with gaps. As these structures are equally seen in the gas and dust phases in our simulations, observations are not limited to one or the other.

We conclude that future simulations designed to study the evolution of protoplanetary discs should never neglect internal photoevaporation, as it appears to have a significant impact on the evolution and chemical composition of the discs, and therefore on the composition of the planets forming in them.

Acknowledgements

We thank Giovanni Picogna, Remo Burn, Manuel Güdel, James Owen and Sho Shibata for helpful discussions. J.L.L. and B.B. acknowledge the support of the European Research Council (ERC Starting Grant 757448-PAMDORA). Th.H. acknowledges the support of the ERC Origins grant number 832428.

Appendix A Side-by-side comparison of the three discs presented in 3.1

Figure A.1 shows a side-by-side comparison of the three disc versions that are discussed in detail in section 3.1. The solid line represents the purely viscous disc, the dashed line the viscous disc with internal photoevaporation via X-rays and the dotted line the viscous disc with a forming planet. In the left and right panels, the gas surface density and the C/O ratio are plotted for the time of 5 Myr, respectively.

Fig. A.1 State of the disc at 5 Myr for the three different scenarios studied in this paper: a viscous disc (solid line), a viscous disc with internal photoevaporation due to X-rays (dashed line), and a viscous disc with a planet sitting at 3 AU (dotted line). Top: gas surface density as a function of disc radius. Bottom: gaseous C/O ratio as a function of disc radius. Plotting and simulation parameters as in Figure 2.

Appendix B Parameter study of the viscous case

We present here the results of a parameter study of a purely viscously evolving disc. An analogous analysis for a viscous disc with additional internal photoevaporation is done in section 3.2. Investigated here is a range of alpha parameters, initial disc masses and radii. We use α = {10⁻⁴, 10⁻³, 10⁻²}, M_disc = {0.01 M_⊙, 0.1 M_⊙} and R_disc = {75 AU, 150 AU), the stellar mass in all simulations remains M_* = 1 M_⊙.

In figure B.1, we plot the gas surface density as a function of disc radius for all the parameters mentioned above. α is increased along the x-axis of the grid. Each row corresponds to a specific combination of initial disc mass and radius, as depicted in the first panel of each row, respectively. Our standard parameter case, as presented in section 3.1.1, is shown in the third panel of the first column.

We note that the largest visible difference between the plots is seen with varying α. Its influence is strongest in the inner disc and the very outer parts. The larger α is, the faster the accretion of inner disc material onto the host star. Additionally, a stronger viscous spreading in the outer disc is seen for bigger values of α.

An impact of the initial disc mass on the simulations is seen in the accretion of inner disc material onto the host star. For larger initial disc masses of M_disc = 0.1 M_⊙, disc material seems to be accreted slightly faster than in the M_disc = 0.01 M_⊙ case. The reason is that the accretion rate correlates with the initial disc mass.

The effect of the initial disc radius on our simulations seems rather small because it changes the initial gas surface density and the corresponding accretion rate just by a factor of 2. However, comparing the first with the second and the third with the fourth row of figure B.1, we see a slight difference in the viscous spreading in the outer disc with varying initial disc radius. For larger initial radii, the viscous spreading is less dominant. This is because the initial disc radius is defined as the location where the exponential cut-off in the initial gas surface density sets in. Discs with larger initial radii therefore have their cut-off further out in the disc, and can therefore better balance the effect of viscous spreading in the outer parts due to more gas being present, leading to less viscous spreading.

Figure B.2 shows the C/O ratio as a function of disc radius for a disc that evolves purely viscously. The parameter grid used here is the same as in figure B.1.

We see again that varying the α parameter has the largest visible effect on our results. For α = 10⁻⁴, the C/O ratio behaves similarly to our standard case, which is depicted in the third panel of the first column (see also section 3.1.1). However, for larger values of α, the initial drop of the C/O ratio is not as high as in the low-viscosity case. Additionally, it increases very rapidly to supersolar values from one time step to the next, only to decrease again, but very slowly this time, until the end of the time evolution. As in the low-viscosity case, pebbles drift faster than gas and therefore, water-rich pebbles reach the inner disc first before carbon-rich gas arrives. They evaporate their water content, enriching the inner disc with oxygen, which leads to an initial decrease in C/O. However, gas is transported faster through the disc due to the higher viscosity. Carbon-rich gas therefore arrives only shortly after, the decrease is not as big as in the low-viscosity case, meaning less water vapour needs to be balanced. As more and more carbon-rich gas arrives quickly, the C/O ratio increases. With still evaporating water-rich pebbles, it balances out over time but stays supersolar.

Varying the initial disc mass and radius seems to have only a minor influence on the C/O ratio in the inner disc. In general, all parameter configurations show similar behaviour, with super-solar C/O ratios in the inner disc regions at the end of the time evolution.

Fig. B.1 Gas surface density for a viscous disc without internal photoevaporation. The gas surface density is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.

Fig. B.2 Gaseous C/O ratio for a viscous disc without internal photoevaporation. The C/O ratio is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.

References

All Tables

All Figures

	Fig. 1 Gas surface density loss rate as a function of disc radius, as given in Equation (5), taken from Picogna et al. (2021).
In the text

Fig. 2 Disc evolution for different disc configurations. Top: viscous disc without internal photoevaporation. Middle: viscous disc with internal photoevaporation due to X-rays. Bottom: viscous disc without internal photoevaporation where a planet seed is put at 3 AU at 0.05 Myr. The planet has a final mass of about 2860 MEarth, and reaches pebble isolation mass at 0.1 Myr. Left: 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. Right: gaseous C/O ratio as a function of disc radius and time (colour-coded). The evaporation lines for different molecules are given as grey dashed lines. Note that the C/O ratio is calculated from number densities scaled to the solar value and that we, by definition, 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 4.

In the text

	Fig. 3 Relative water content as a function of disc radius and time. Internal photoevaporation is active in this simulation. Note that the relative water content is calculated as a ratio of gas surface densities, and is therefore shown as a mass fraction here. Colour coding, plotting and simulation parameters as in Figure 2.
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. Left: carbon over hydrogen. Middle: oxygen over hydrogen. Right: nitrogen over hydrogen. Colour coding, plotting and simulation parameters as in Figure 2.
In the text

Fig. 5 Gas surface density for a viscous disc with internal photoevaporation due to X-rays. The gas surface density is shown as a function of disc radius and time. A variety of parameter configurations is used for this study: along the x-axis of this figure, α is varied from α = 10−4 to α = 10−2, whereas on the y-axis, different sets of disc masses and radii are used, varying in the range of Mdisc = 0.01–0.1 M⊙ and Rdisc = 75–150 AU. Colour coding, plotting and stellar mass as in Figure 2.

In the text

	Fig. 6 Gaseous C/O ratio for a viscous disc with internal photoevaporation due to X-rays. The C/O ratio is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.
In the text

	Fig. 7 Disc evolution for a viscous disc with internal photoevaporation due to X-rays but with a smaller mass loss rate of 3.6 · 10−8 M⊙/yr instead of 3.864 · 10−8 M⊙/yr. Left: gas surface density as a function of disc radius and time. Right: gaseous C/O ratio as a function of disc radius and time. Colour coding, plotting and simulation parameters as in Figure 2.
In the text

Fig. 8 Phases of disc evolution, with the left column depicting a pure viscous disc, the middle column a viscous disc with internal photoevaporation, and the right column a viscous disc with a growing planet. Top panel: initial state of the protoplanetary disc, with its inner part being slightly enriched with water vapour, its middle part being slightly enriched with carbon vapour and its outer regions showing no enrichment. The evaporation line of water and the combined one of all carbon-bearing molecules are marked with black, curved lines. Their curvature is the result of a temperature gradient between the mid-plane and surface layers of the disc. Pebbles are distributed evenly across the full disc, showing icy layers according to their position in it. Left column: during the pure viscous evolution, pebbles drifting inwards evaporate their volatile content whenever crossing an evaporation front. Due to the position of the ice lines and the fact that pebbles drift faster than gas, the inner disc is first enriched with water vapour (top panel). With time, carbon vapour arrives (middle panel) and slowly enriches the inner disc (lower panel), while the gas is accreted onto the host star. Middle column: with internal photoevaporation, the disc evolves similarly to the pure viscous case in the beginning (first and second panel). When photoevaporation opens a gap, gas diffusing into it is carried away by photoevaporative winds while pebbles are blocked by the gap itself (third panel). The inner disc is then enriched with water vapour again due to the water equilibrium cycle (fourth panel). All volatiles, except water, are accreted onto the host star. Water contrarily recondenses at the water-ice line and forms water pebbles that evaporate again, contributing to a full cycle of evaporation and condensation that allows the water vapour to survive for a long time. Right column: a planet forming between the water evaporation front and the ice lines of carbon-bearing molecules leads to the same chemical evolution of the disc as in the pure viscous case. When the planet is massive enough, it opens a gap, but contrary to the gap caused by photoevaporation, a gas flow to the inner disc is still possible. The only differences are less water enrichment in the beginning due to the planetary embryo already blocking some inflow (top panel) and a higher carbon enrichment in the late stages (lower panel).

In the text

Fig. A.1 State of the disc at 5 Myr for the three different scenarios studied in this paper: a viscous disc (solid line), a viscous disc with internal photoevaporation due to X-rays (dashed line), and a viscous disc with a planet sitting at 3 AU (dotted line). Top: gas surface density as a function of disc radius. Bottom: gaseous C/O ratio as a function of disc radius. Plotting and simulation parameters as in Figure 2.

In the text

	Fig. B.1 Gas surface density for a viscous disc without internal photoevaporation. The gas surface density is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.
In the text

	Fig. B.2 Gaseous C/O ratio for a viscous disc without internal photoevaporation. The C/O ratio is shown as a function of disc radius and time. We vary the disc parameters as in Figure 5. Colour coding, plotting and stellar mass as in Figure 2.
In the text