Open Access
Issue
A&A
Volume 701, September 2025
Article Number A239
Number of page(s) 22
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202555195
Published online 22 September 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

The composition of dust and gas in protoplanetary discs determines the composition of any planets forming therein. Consequently, elemental abundances such as C/H and O/H and particularly ratios such as C/O have been suggested as tools to trace planetary formation locations and mechanisms. At overall solar composition, the C/O ratio of both the solids and gas in a static disc would rise with distance from the central star as successive snowlines of H2O, CO2, and CO are crossed (Öberg et al. 2011b). However, there is a developing understanding that the picture is more complex than implied by such static models. Disc evolution models show significant redistribution of molecules resulting from the rapid radial drift of dust grains, the sublimation of molecular ices as these grains cross snowlines, and the subsequent slow accretion of this enriched gas inwards (Öberg & Bergin 2016; Booth et al. 2017; Krijt et al. 2018). Further complications include evolution of the disc temperature causing the snowlines to migrate (Miley et al. 2021), migration of planets themselves (Turrini et al. 2021; Pacetti et al. 2022), trapping of volatile species in less-volatile ices (Collings et al. 2004; Ligterink et al. 2024), vertical temperature gradients that lead to ‘snow surfaces’, and chemical evolution of the gas and ice (Eistrup 2023).

Confirming how abundances vary across a disc thus requires comparing such models to measurements from observations. While the Atacama Large Millimeter/submillimeter Array (ALMA) is sensitive to material in the outer disc (∼10s au) (e.g. Kama et al. 2016a; Bergner et al. 2019; Miotello et al. 2019; Öberg et al. 2021; Sturm et al. 2022; Bergin et al. 2024), observations in the infrared, such as with JWST’s Mid-Infrared Instrument (MIRI) in the Medium Resolution Spectroscopy (MRS) mode, are better suited to characterising warm gas in the inner disc (<few au).

So far, the inner disc chemistry depends strongly on stellar mass. In discs around stars ≳0.3 M – for example, T Tauri stars (TTS) – the most frequently detected molecules are H2O, CO2, C2H2, and HCN (Grant et al. 2023, 2024; Banzatti et al. 2023, 2025; Gasman et al. 2023, 2025; Henning et al. 2024; Pontoppidan et al. 2024; Romero-Mirza et al. 2024b; Schwarz et al. 2024; Temmink et al. 2024; Salyk et al. 2025; Vlasblom et al. 2025; Arulanantham et al. 2025), with retrieved emitting areas that imply emission mostly inside ≲1 au. This pattern, with strong O-bearing species, is thought to occur when O is in excess with respect to C, i.e. C/O<1. Conversely, discs around very low mass stars (VLMSs) ≲ 0.3 M show stronger C2H2 relative to the other molecules (Grant et al. 2025), strengthening earlier Spitzer analyses (Pascucci et al. 2009, 2013). In addition, MIRI-MRS has revealed and resolved a rich hydrocarbon chemistry with column densities orders of magnitude larger than previously thought (Tabone et al. 2023; Arabhavi et al. 2024; Kanwar et al. 2024; Arabhavi et al. 2025a). Thermochemical models imply that C/O>1 is required to achieve such conditions (Najita et al. 2011; Kanwar et al. 2024), which are found at a wide range of ages (Long et al. 2025). Notable exceptions to this pattern are DoAr 33 (Colmenares et al. 2024), a 1.1 M star with a hydrocarbon-rich chemistry consistent with C/O=2–4, and the H2O-rich inner disc of the VLMS Sz 114 (Xie et al. 2023).

Two broad explanations have been proposed for the high C/O cases (Tabone et al. 2023; Arabhavi et al. 2024). On the one hand, C/H could be increased by the destruction or erosion of a refractory carbon reservoir. Conversely, the O/H could be depleted in the inner disc gas as a result of the chemical evolution expected due to transport (Mah et al. 2023) or retention of O-rich ices in dust traps or planetesimals (Najita et al. 2013; McClure et al. 2020). Distinguishing between these cases is important, as they predict distinct C/H and O/H for a given C/O.

In the erosion scenario, refractory C could be liberated at high temperatures by pyrolysis of organic solids (their thermal decomposition at ≳500 K) at the ‘soot line’ (Kress et al. 2010; Li et al. 2021; Houge et al. 2025), or chemical erosion, such as oxidation reactions (Lee et al. 2010) or hydrogen chemisputtering (Lenzuni et al. 1995; Borderies et al. 2025). The sublimation of amorphous C is less likely, as it occurs at considerably higher temperatures (≳1100 K, Gail & Trieloff 2017). Alternatively, the radiation field can lead to the destruction of grains or Polyaromatic Hydrocarbons (PAHs) by ultraviolet photons (photolysis, Anderson et al. 2017; Bosman et al. 2021; Vaikundaraman et al. 2025) or X-rays (Siebenmorgen & Krügel 2010).

In the O-depletion scenario, ices desorbing from inwardly drifting dust deliver less-volatile, more oxygen-rich molecules closer to the star. Consequently, the inner disc is expected to first go through a H2O-enriched phase with high O/H and low C/O. This is followed by an increase in C/O as the H2O accretes onto the star, and successively carbon-rich molecules - CO2 and later CO and CH4 - arrive (Mah et al. 2023, 2024; Sellek et al. 2025). The strength and duration of these phases can depend on the formation, location, and leakiness of any dust traps that block drift. Recent observations of VLMSs suggest that the inner disc may indeed become increasingly carbon rich when comparing older discs to younger ones (Long et al. 2025; Arabhavi et al. 2025a).

Crucially, Mah et al. (2023) argue that the colder temperatures of discs around VLMSs result in the snowlines lying closer to the star and thus a faster progression through this evolutionary pattern. Therefore, at a given age, discs around VLMSs are likely to show a higher C/O. The disc transport models of Mah et al. (2023) do cross C/O=1 (required for rapid hydrocarbon formation) but only due to the inclusion of abundant CH4 (as a species containing C but not O). However, those transport models never suggest a C/O>1.5, despite such values being inferred for some discs (Colmenares et al. 2024; Kanwar et al. 2025; Long et al. 2025). This is despite making generous assumptions for the initial abundance of CH4 in the disc (see also discussion in Houge et al. 2025) that vastly exceed the amount inherited from ices in protostellar envelopes (Öberg et al. 2011a; Boogert et al. 2015) or found in cometary ices (Mumma & Charnley 2011; Altwegg et al. 2019; Schuhmann et al. 2019). The question thus remains open whether transport alone explains the difference in C/H, O/H, and C/O between VLMSs and TTSs or if in situ C/H enhancement is required.

The fundamental barrier to achieving a high C/O ratio is that it is buffered by CO being the most abundant C-bearing molecule in these disc models. With C atoms thus always being associated with around one O atom, transport alone is unable to avoid C/O∼1. Raising C/O further thus requires destroying CO in a way such that the C atoms are liberated in the process.

It is indeed often found that CO and its isotopologues – as observed with ALMA – are less abundant in the outer disc around TTSs than the ISM, typically by factors 10–100 (Favre et al. 2013; Kama et al. 2016b; McClure et al. 2016; Anderson et al. 2022; Trapman et al. 2022; Paneque-Carreño et al. 2025). This so-called CO depletion sets in on timescales of ∼1 Myr (Harsono et al. 2014; Zhang et al. 2020; Bergner et al. 2020). While some effect is expected due to freeze-out and photodissociation, the term ‘volatile CO depletion’ refers to a lack of CO on top of these processes (Miotello et al. 2017). Models have explored physical routes – involving the sequestering of CO ice onto pebbles that settle to the midplane and drift (Kama et al. 2016b; Krijt et al. 2020; Powell et al. 2022) – and chemical routes – involving the conversion of CO molecules to less volatile species (Yu et al. 2016; Molyarova et al. 2017; Bosman et al. 2018; Schwarz et al. 2019; Krijt et al. 2020; Van Clepper et al. 2022; Zwicky et al. 2025) – to achieving this depletion1. For our purposes, however, physical sequestration does not liberate C, so we must also consider chemical routes.

Conversely, liberating C from CO cannot by itself alter C/O (or lead to C/O>1) if the liberated O is not removed. While other works also argue for the importance of liberating the C from CO in order to produce hydrocarbons, for example as a result of ionising X-rays (Walsh et al. 2015), they employ physically static models. However, any O-bearing species produced are generally less volatile than the hydrocarbons. Thus, the differential transport of vapour and ice due to dust radial drift, which is not considered in the thermochemical models, could separate the liberated C and O reservoirs and produce conditions where C/O>1.

In this work, we implement simple prescriptions for chemical transformation of CO in a disc evolution code (Booth et al. 2017), as described in Section 2, in order to study whether this can prevent CO transport to the inner disc in a way that results in C/O>1. We thus explore when depletion of O as a result of global-scale disc evolution could explain high C/O measurements in the inner disc via grids of such models. Section 3 presents the evolution of C/H, O/H, and C/O for these grids, while Section 4 considers which conditions would need to be met for discs around VLMSs to have a higher C/O than those around TTSs and how to go about determining if CO-depleted outer disc gas is indeed responsible for the elevated C/O. We summarise our conclusions in Section 5.

2 Methods

2.1 Disc evolution model

We used the 1D disc evolution code by Booth et al. (2017), including the updates to the dust transport in Sellek et al. (2020) and the chemistry in Sellek et al. (2025), to model the evolution of discs around 1 M and 0.1 M stars (Table 1). We modelled the discs for 10 Myr, which spans the ages of the nearby starforming regions studied in Cycles 1 & 2 with JWST, including the older Upper Sco region; the endpoint must be defined manually since we do not include any process to disperse the discs (see Lienert et al. 2024, 2025, for models including photoevaporation). In brief, the code computes the viscous evolution of the disc’s gas according to the Shakura & Sunyaev (1973) α prescription where ν = αcS H for sound speed cS and scale height H and we assume α = 10−3. On top of this, the dust evolution is computed using the two-population model (Birnstiel et al. 2012) implemented following Laibe & Price (2014); Dipierro et al. (2018). In this approach, the majority of mass is contained in large grains, which have sizes limited by either drift or fragmentation and which undergo significant radial drift, while there is always a remnant population of small grains that remain coupled to the gas. The grains are assumed to fragment at a velocity of 1 m s−1 for bare dust grains and 10 m s−1 for icy dust grains.

Several gas and ice phase molecular tracers are advected following the motion of the gas or dust components respectively. At each time step, the equilibrium between freeze-out and thermal desorption is calculated to determine the fraction of each molecule in the gas and ice phases. The species included in this work are those used by Sellek et al. (2025) (H2O, CO2, CO, CH4, CH3OH, O2, C-grains, and Si-grains), plus we additionally investigate C2H2, C2H4, and C2H6. All prefactors and binding energies used (that for the hydrocarbons are for pure multilayer ices from Behmard et al. 2019), as well as resulting snowline locations for each M*, are given in Table A.1.

For the initial conditions, we assumed a surface density profile given by (Lynden-Bell & Pringle 1974) Σt=0(R)=MD,02πRC,01Rexp(RRC,0),\Sigma_{t=0}(R) = \frac{M_{\rm D,0}}{2\pi R_{\rm C,0}} \frac{1}{R} \exp\left(-\frac{R}{R_{\rm C,0}}\right) \label{eq:LBP} ,(1)

where we explore several values for the initial mass, MD,0, and radius, RC,0 (Table 1). All of our MD,0 are ≲ 0.1 M*, so they represent discs in the gravitationally stable regime. This regime may follow a period of gravitational instability maintained by infall, where non-axisymmetric substructures and temperature profiles that cannot be captured in 1D form; although these effects can already leave an imprint in the distributions of different molecules (Molyarova et al. 2025), we use uniform initial abundances (Appendix A) in this work. The RC,0 follow the range of values inferred by Trapman et al. (2023) from CO gas disc radii and are used for both M* investigated. The Σ(R) profiles for t = 01 and 5 Myr, are shown for fiducial parameters in Fig. A.1.

The temperature is prescribed as a power-law proportional to R−1/2 and scaled based on the aspect ratio at 1 au (h0), where we assumed h0M*−0.425 (Sinclair et al. 2020): T0=GMμRRh02.T_0 = \frac{GM_*\mu}{R\mathcal{R}}h_0^2 .(2)

This profile is assumed to be constant. We do not model the effects of stellar evolution (unlike e.g. Miley et al. 2021) or heating from hydrodynamic processes in the disc (unlike e.g. Molyarova et al. 2025). The profile also applies to both dust and gas since they couple thermally at midplane densities.

While Equation (1) defines a ‘smooth’ density profile that is free of substructures or dust traps, we also considered perturbations to the gas surface density. These were introduced by varying the viscosity locally, centred around a radius, Rgap, and with an amplitude, Agap, where the perturbation is fully developed by time, tf (see Sellek et al. 2025, for further details). We chose either Rgap = 5 au, which creates a ‘warm’ trap, or Rgap = 50 au, which leads to a ‘cold’ trap. These terms reference the fact that they lie at respectively warmer and colder temperatures than the freeze-out of CO (and CH4). Discs have been observed to show differences in the strength of their C2H emission (related to their outer disc C/O ratio) between these cases (van der Marel et al. 2021). While the most easily resolved substructures lie at 10s au (Bae et al. 2023), high-resolution surveys have shown they can form at distances as close as 5 au in a small fraction (~10%) of compact discs (Guerra-Alvarado et al. 2025). Moreover, we modelled cases where traps a) are present from the start (tf = 0), which maximises their chances of impacting pebble drift and volatile delivery (Mah et al. 2024; Sellek et al. 2025), or b) only form by tf = 1 Myr, as they are more rarely observed in the earlier embedded stages (Ohashi et al. 2023; Hsieh et al. 2024). This lack of substructures at earlier times may be due to observational limitations or the fact that planets (if they are indeed the origin of traps) may take megayear timescales to form and grow to the pebble isolation mass, especially since it is harder for them to open a gap in younger, warmer, discs (Nazari et al. 2025).

Table 1

Parameters explored in our model grids. Fiducial values are highlighted in bold.

2.2 Reactions

The chemistry in disc midplanes is driven by ionisation (Sections 2.2.3, 4.1.1), which is often assumed to result from cosmic rays (CRs). Notably, CRs also produce a secondary UV field (Prasad & Tarafdar 1983). Molecules can be photodissociated by this UV – or react with the ions or radicals created by the CRs or secondary photons – with rates that imply changes on megayear timescales for typical ionisation rates. We do not include a detailed chemical network, but instead a toy model based on results from thermochemical disc models (e.g. Walsh et al. 2015; Yu et al. 2016; Eistrup et al. 2016, 2018; Bosman et al. 2018). Our approach is summarised in Figure 1 and involves simple prescriptions for the net effect of major chemical pathways, based on the initial rate-limiting steps, similar to that used by Krijt et al. (2020). For each time step of the disc evolution code (typically a few years), the rates were calculated and the abundances of molecules updated in an Eulerian manner; the time step for the disc evolution was taken as the minimum of the shortest time step required to model the transport processes and the shortest time required to deplete any molecular tracer by a factor e in order to ensure that no abundance ever becomes ≤0.

thumbnail Fig. 1

Summary of the reaction scheme(s) used in this work. Here, ‘s-’ indicates that the species is in the ice phase. The process in bold is the most critical for raising the inner disc C/O.

2.2.1 CO transformation

Bosman et al. (2018) characterised CO destruction in disc midplanes by three main routes. In terms of the net consumption of our molecular reservoirs, they are as follows:

  1. Dissociative ionisation of gas phase CO by He+, where we assume that the liberated C+ reacts to form CH4 CO+3H2He+CH4+H2O,\ce{CO + 3 H2 ->[He^+] CH4 + H2O}(3)

  2. Hydrogenation of CO ice on grain surfaces, which forms CH3OH via HCO CO(s)+2H2H​CH3OH(s),\ce{CO(s) + 2 H2 ->[H] CH3OH(s)}(4)

  3. Reaction of CO with OH on grain surfaces to form CO2 CO(s)+H2O(s)OHCO2+H2.\ce{CO(s) + H2O(s) ->[OH] CO2 + H2}(5)

Above each arrow are intermediaries, created by ionisation or photodissociation, required for the reactions (Appendix A).

We treated grain-surface reactions as tunnelling processes following Hasegawa et al. (1992), where the top two surface layers are active (Bosman et al. 2018). The diffusion was assumed to be dominated by H (for reactions it participates in) at the maximum of the thermal hopping and tunnelling rates. The barrier to diffusion was set as 0.3 of the binding energy (for H we took the amorphous solid water value in Minissale et al. 2022). The barrier width for the reactions was taken to be the standard 1 Å, and the barrier heights are equal to those in gas-grain codes used for the thermochemical modelling of protoplanetary discs (Walsh et al. 2015, see Table A.2). The 300 K barrier for O2 + H→ HO2 used in that network (Lamberts et al. 2013; Taquet et al. 2016) results in a very short reaction timescale. For computational efficiency, we thus assumed that any frozen-out O2 is entirely consumed (otherwise this process seriously throttles the code time step). We checked that this assumption had negligible impact on our results: the reaction rate of O2 is fast enough that in either case, very little O2 survives for the megayear timescales that we are interested in. Gas-phase reaction rates were taken from the UMIST RATE22 database (Millar et al. 2024).

2.2.2 Further reactions

CH3OH ice, for example formed by hydrogenation of CO, may undergo UV photodissociation. For simplicity, we followed Krijt et al. (2020) and assumed that this always splits the C–O bond, with the CH3 and O being rapidly hydrogenated to CH4 and H2O: CH3OH(s)+H2uVCH4(s)+H2O(s).\ce{CH3OH(s) + H2 ->[UV] CH4(s) + H2O(s)}(6)

We note, however, that the branching ratios for photodissociation of CH3 OH are uncertain: other works (Laas et al. 2011; Drozdovskaya et al. 2014) adopt a ‘standard’ assumption of 60% CH3 production, while laboratory experiments suggest it could be even lower (≲14% Öberg et al. 2009); in this sense, our rate is an upper limit for the conversion of CH3OH ice to CH4.

Moreover, the most commonly observed hydrocarbon in the inner regions of protoplanetary discs is C2H2. We assumed this forms after photodissociation of gas-phase CH4 CH4uv12C2H2+32H2.\ce{CH4 ->[UV] 1/2 C2H2 + 3/2 H2}(7)

On grain surfaces, C2H2 may be destroyed by hydrogenation (Aikawa et al. 1999; Yu et al. 2016; Bosman et al. 2018). The radicals C2H3 and C2H5 are highly reactive, with no barrier for hydrogenation (Hasegawa et al. 1992), and hence they are very short lived. Thus we considered only the more stable C2H4 and C2H6: C2H2+H2H​C2H4\begin{align} \ce{C2H2 + H2 &->[H] C2H4} \\ \end{align}(8) C2H4+H2H​C2H6.\ce{C2H4 + H2 &->[H] C2H6}.(9)

The initial abundances of these hydrocarbons were assumed to be 0 everywhere so we could focus on their production routes.

Another way to liberate additional C is by destroying or eroding carbonaceous grains. In Appendix B.3 we estimate UV photolysis of grains in the outer disc to have little effect on our models, and so henceforth neglect this process in this work. We refer the reader to Houge et al. (2025) for an exploration of the role of thermal decomposition of grains in the inner disc.

2.2.3 Ionisation rates

Every process discussed above depends indirectly on the ionisation rate. At the midplane, CRs are often assumed to be the principal source of ionisation – with a fiducial rate assumed in databases such as umist (Millar et al. 1991, 2024) of ζH2 = 1.3 × 10−17 s−1 (a dense cloud value; Prasad & Huntress 1980) – and thus drive the chemistry (Eistrup et al. 2016, 2018). CRs are subject to attenuation, with an average penetration depth of ΣCR = 96 g cm−2 (Umebayashi & Nakano 1981). For simplicity, we assumed they are attenuated by the vertical column of gas2.

When CRs are highly attenuated, ionisation does not cease entirely. Short-lived radionuclides (SLRs) can provide ionisation rates of up to 7.6 × 10−19 (Umebayashi & Nakano 2009), though this value will decay over time as the SLRs have half-lives of around 1 Myr (e.g. Cleeves et al. 2014; Eistrup et al. 2018). This rate is too low to produce much chemical evolution within typical protoplanetary discs lifetime (Eistrup et al. 2018) but we included it as a floor value nonetheless.

The ionising particle flux incident on a disc can also vary from the fiducial value. For example, the flux of galactic CRs reaching Earth is modulated by solar magnetic activity; the increased activity of young stars may reduce the ionisation rate to as little as ∼10−20 s−1 even before attenuation in the disc (Cleeves et al. 2013). Conversely, energetic particles may be accelerated locally by flares from the central star (Feigelson et al. 2002; Rab et al. 2017; Rodgers-Lee et al. 2017, 2020a), reconnection events in the interaction region between the stellar and disc magnetic fields (Orlando et al. 2011; Colombo et al. 2019; Brunn et al. 2023), or shocks in protostellar jets (Padovani et al. 2015, 2016). Finally, hard X-rays also produce ionisation, with midplane ionisation rates typically of the order of 10−20−10−19 s−1 (Cleeves et al. 2013, 2014; Rab et al. 2017). To encompass these possibilities in a general way, agnostic to the origin of the variations, we defined a single ionisation rate incident on the disc ζH2 and explored values raised or lowered uniformly by a factor ten with respect to the fiducial CR value (Table 1).

thumbnail Fig. 2

Evolution of different disc models for 10 Myr. a) Models that do not include reactions of C-bearing species for various initial disc-to-star mass ratios and initial disc sizes. b) Models that include reactions of CO. All models have a fixed initial disc-to-star mass ratio of 0.05 but different combinations of ionisation rate and disc size. In each panel, the tracks start at ‘S’ and evolve in the direction of the arrows, ending at ‘F’. The solid sections represent ages of 1–5 Myr appropriate for most discs studied so far with MIRI-MRS. In each block, the left-hand columns show models around a 1.0 M star, representative of TTSs, while the right-hand columns show models around a 0.1 M star, representative of VLMS. The top rows plot the discs in the plane of the inner-disc gas-phase C/H vs gas-phase O/H. Dashed lines indicate where C/O=0.32,1,3,10. The bottom rows plot the gas-phase C/O in the inner disc over time. The horizontal dashed lines show ratios of 0.32 (initial ratio of volatile abundances), 0.43 (ratio of total abundances) and 1 and the lightly shaded regions indicate the typically observed range of C/O at 1–5 Myr for each stellar mass regime.

3 Results

3.1 No C-chemistry

As a baseline, we consider a model set without any reactions of C-bearing species; the models only include the oxidation of O2 ice to H2O (otherwise O2 could provide a large, unconstrain-able volatile-O reservoir). Figure 2a shows the evolution of these models in the plane of the gas-phase C/H and O/H in the inner disc and the evolution of the resulting C/O over time. We measured this at the inner boundary of our disc model (0.1 au), but it should be representative of the conditions inside the H2O snowline. We shade regions where C/O<1 for TTSs and C/O>1 for VLMSs to illustrate the typical constraints inferred without placing any weight on exact individual ages or C/O values, which vary depending on the method used to derive them.

Figure 2a shows that initially O/H is enhanced (with a corresponding drop in C/O) as radial drift delivers O-rich H2O ice to the inner disc. As the H2O vapour is accreted onto the star, lowering O/H again, CO2 becomes the dominant molecule in the inner disc (Sellek et al. 2025), leading to an increase in C/H and an overall rise in C/O. The initial disc size RC,0 has a small effect on the evolution. For more extended discs, a greater fraction of the molecular budget is initially frozen out, resulting in slightly larger enhancements in these phases. The evolution happens faster for lower-mass stars, reflecting the snowlines lying closer to the central star (Mah et al. 2023). The initial disc-to-star mass ratio has negligible effect.

Over time, CO becomes the dominant molecular species, and the C/O reaches ∼1. This is the end state after 10 Myr for 1.0 M stars, while for VLMSs, the closer-in snowlines result in enough time to deplete the CO vapour significantly through accretion onto the star, resulting in lower final C/H and O/H. At late times, a small number of dust grains, that originate at large radii and never grew large enough to undergo significant radial drift, remain in the disc. These, by definition, experience the least decoupling from the gas and so produce little segregation of C- and O-bearing molecules. Consequently, after the remnants of the earlier phases of strong radial drift completely accrete onto the star (only for VLMS discs does this happen within 10 Myr), these small grains reset the elemental composition to the initial volatile C/O=0.32 assumed in the models. A gradual fall in C/O towards this value is thus seen at the latest times ≳7 Myr. Other than in this very last phase, at a given time the VLMS discs show higher C/O as a result of their faster evolution (Mah et al. 2023).

The C/O of the disc can be calculated as CO=CO+CH4+CO2+CH3OHCO+2×O2+2×CO2+CH3OH+H2O.\frac{\rm C}{\rm O} = \frac{\rm CO+CH_4+CO_2+CH_3OH}{\rm CO+2\times O_2+2\times CO_2+CH_3OH+H_2O} .(10)

In the absence of all O-bearing species other than CO, which all have shorter lifetimes in the disc than CO (including O2 due to the efficient oxidation to H2O), this reduces to CO=CO+CH4CO.\frac{\rm C}{\rm O} = \frac{\rm CO+CH_4}{\rm CO} .(11)

As CH4 is only marginally less volatile than CO, their ratio is only weakly altered by the relative dust and gas evolution. Thus, for our assumed inherited initial CH4/CO=0.02, C/O cannot rise significantly above 1. Somewhat greater C/O is seen by Mah et al. (2023), as their models use a much higher initial CH4 /CO=0.5. However (as also shown in Fig. 8 of Long et al. 2025), C/O still rarely exceeds ~1.2. This is comfortably below the inner disc C/O≳ 2 inferred for several VLMSs (Long et al. 2025; Kanwar et al. 2025) or DoAr 33 (Colmenares et al. 2024). We note that CH4/CO would be further lowered if most of the CH4 were trapped in H2O ices such that it only desorbed at higher temperatures (e.g. Collings et al. 2004).

We conclude that the destruction of CO, along with a considerable additional source of hydrocarbons, would be required to reach C/O>1. Next, we explore how to achieve such conditions.

3.2 CO transformation

As the dominant molecule at late times in the models above is CO, which imposes C/O=1, it is imperative that it is depleted in favour of other species if C/O>1 is to be reached. Reactions alone only reportion a fixed C and O budget between different molecules, but cannot alter the (total) C/O ratio. Instead, C and O must be separated into different molecules with the C remaining in a volatile species and the O becoming sequestered in a less volatile species (e.g. H2O) and transported more rapidly out of the disc by radial drift. In the gas phase, this can be achieved by the formation of CH4. Where CO is comfortably below its freeze-out temperature, it is predominantly destroyed by hydrogenation to CH3OH, which if followed by UV photodissociation, may potentially also lead to CH4. In Appendix B.1 we clarify the role played by the different reactions in our CO destruction scheme by adding them one by one for an ionisation rate ζH2 = 1.3 × 10−17 s−1. We demonstrate that the most crucial reaction is the gas phase route as it directly liberates C. CH3 OH formation alone does not liberate C and cannot produce C/O>1, and in smooth discs, the timescale to photodissociate the CH3 OH ice is too short compared to the drift timescale to rectify this.

Here, we focus on the results of a grid of models including all CO transformation routes. As the models without reactions involving C showed a much stronger dependence on the initial disc size than the initial disc-to-star mass ratio, in Fig. 2b we only show models with a fiducial MD,0/M* = 0.05 but note that models with other mass ratios show very similar behaviours, with a tendency towards slightly faster chemical evolution for lower mass discs as they attenuate CRs slightly less.

For a given set of properties (MD,0/M*, RC,0, ζH2), the evolutionary pathways of discs around 1.0 and 0.1 M stars are similar (Fig. 2b). As C/O>1 is reached by destroying CO, such that the O ends up in less-volatile molecules with a shorter lifetime in the disc, it happens primarily through mild O depletion at roughly solar or slightly enhanced C/H. The maximum C/O occurs when the gas is slightly depleted in C, but very depleted in O. As seen when reactions are neglected, the evolution is more rapid around the VLMSs as predicted by Mah et al. (2023). Thus, as before, by the arbitrary end time of 10 Myr, VLMSs show significant depletion of both C/H and O/H, while the TTSs end with somewhat depleted O/H, but still relatively high C/H, and, consequently often a higher C/O.

3.2.1 Dependence on ionisation rate

The clearest trend in Fig. 2b is that higher C/O is achieved for higher ionisation rates. Essentially all models for both stellar masses remain under C/O=1 when the ionisation rate is decreased to ~10−18 s−1. The fiducial ionisation rate is able to create a maximum 1 ≲ C/O ≲ 2. However, a highly elevated ionisation rate of ~ 10−16 s−1 can reach 2 ≲ C/O ≲ 10, consistent with some of the more extreme cases inferred from observations.

This trend can be understood by considering the various timescales at play, relative to the 10 Myr for which we run our models (only a few rare exceptions, for example J0446B as studied by Long et al. (2025), live longer than this). The average timescale that a molecule, m, spends in the gas phase before accreting onto the star after being released at its snowline can be estimated following Eq. (18) of Sellek et al. (2025) as tacc,m=0.28Myr(MM)3/2(α103)1(h00.021)2(Tsnow,m100K)2.t_{\rm acc,m} = 0.28\,\mathrm{Myr} \left(\frac{M_*}{M_{\odot}}\right)^{3/2} \left(\frac{\alpha}{10^{-3}}\right)^{-1} \left(\frac{h_0}{0.021}\right)^{2} \left(\frac{T_{\rm snow,m}}{100\,\mathrm{K}}\right)^{-2} .(12)

Given that we scale h0M*−0.425 (Sinclair et al. 2020), then for our parameters and assuming Tsnow,cO = 20 K, we find tacc,cO = 7Myr(M*/M)0.65, i.e. 1.6Myr and 7Myr for 0.1 and 1.0M stars respectively. For comparison, the timescale for the most critical reaction, COHe+CH4$\ce{CO ->[He+] CH4}$, is 4Myr(ζH2/1017s1)1${\sim}4\,\mathrm{Myr}\left(\zeta_{\rm \ce{H2}}/10^{-17}\,\mathrm{s^{-1}}\right)^{-1}$ assuming CO dominates the consumption of He+.

Thus, if the ionisation rate is lowered, the reaction rate becomes longer than both the 10 Myr for which we run the models, and the advection timescales, so most CO molecules reach the star without a chance of being converted to CH4 and maintain a C/O~ 1. Conversely if the ionisation rate is raised, CO molecules can react in the gas phase in <1 Myr, much faster than it would take to accrete onto the star around stars of any mass. Consequently, the destruction of CO and formation of CH4 is efficient in these high ζH2 models, resulting in a high C/O within 2-3 Myr for any M*. Only for the fiducial ionisation rate of ~ 10−17 s−1 do we find that, within the timescale of most observations so far, C/O<1 for the TTSs and >1 for the VLMSs.

3.2.2 Dependence on disc size

The second trend in Fig. 2b is a dependence on initial disc size, which becomes especially clear at high ionisation rate. Intrigu-ingly, initially smaller discs around VLMSs reach higher C/O, while for the TTSs, initially larger discs are more C-rich for much of their early evolution. The sensitivity to this parameter is such that at the fiducial ionisation rate, whether discs around VLMSs cross the C/O=1 boundary depends on their initial size. Larger (≳50 au) discs around VLMSs do not reach C/O>1 without an elevated ionisation rate.

The short accretion timescales at the CO snowline for the VLMSs mean that they are observed over several such timescales. Thus, their behaviour is controlled by how quickly the reservoir of C-rich gas is able to be depleted by accretion to a level where sublimation from the remnant small dust grains comes to dominate and resets C/O to ~0.32. The larger the gasphase reservoir of CO relative to the abundance of small grains, the later this happens and the more time there is to convert CO to CH4 and raise C/O. Smaller discs have a larger fraction of CO (94% for 5 au discs compared to 24% for 50 au discs) initially in the gas phase (which prevents it from instead being converted to CH3OH). Moreover, they have a shorter radial drift timescale such that less dust is left at late times. Combined, these mean that the discs reach higher later peaks in C/O.

Conversely, the TTSs are assumed to have denser discs, leading to more significant attenuation of ionising particles. In larger discs, the mass is more spread out, making the disc optically thinner. Consequently, a larger fraction of the disc CO gas (46% for 50 au discs compared to 8% for 5 au discs) experiences ionisation rates sufficiently above the ∼10−18 s−1 floor (that alone is insufficient to destroy CO within a disc lifetime).

thumbnail Fig. 3

C/O ratio as a function of the inner disc H2O vapour mass (calculated as the total integrated H2O vapour mass inside of the midplane snowline) for the same models as Fig. 2b. The line styles are the same as in that figure, and the arrows illustrate the evolution with time.

3.2.3 Relationship with the presence of water

As outlined above, small amounts of O-rich ices sublimating from remnant dust can contribute to limiting the growth of C/O. Figure 3 shows the C/O ratio as a function of the total integrated H2O vapour mass in the inner disc. This reduces the spread between the tracks for different initial radii, demonstrating the evolution of the disc’s solids is important even though the critical destruction of CO happens in the gas phase. The higher the ionisation rate, the sooner CO is broken up, and the sooner the O (in the form of H2O) can be removed by the action of drift. Thus we see that while fiducial ionisation rates make H2O slightly more abundant at 10 Myr as they produce some H2O from CO, shifting the tracks to the right compared to the lowest ionisation rates, the highest ionisation rates result in much lower H2O after 10 Myr.

We note, however, that this may not imply an anti-correlation between C/O and the observable H2O emission. The observable H2O emission is limited by photodissociation in the UV-exposed upper layers, and the dust continuum optical depth, where the amount of dust may be correlated with the amount of H2O delivered to the inner disc resulting in little sensitivity of the column densities to the H2O abundance (Sellek et al. 2025; Houge et al. 2025). Indeed, H2O is detected in a few VLMS discs, and at line-to-continuum ratios similar to TTS discs, but is harder to identify due to the strong C2H2 emission (Arabhavi et al. 2025b).

3.2.4 Formation of larger hydrocarbons

We note that CH4 has rarely been detected so far in MIR spectra of discs around stars >0.3 M, with tentative exceptions of DoAr33 (Colmenares et al. 2024) and CY Tau (Temmink et al. 2025). Rather, one of the most commonly detected species in MIR spectra of protoplanetary discs is C2H2. VLMS disc spectra do show CH4 (Arabhavi et al. 2024; Kanwar et al. 2024), but also show larger hydrocarbons including C4H2, C6H6, C3H4, C2H6 and C2H4 (in order of detection frequency, Arabhavi et al. 2025b). Thus, for completeness, we henceforth use models in which the hydrocarbon formation progresses further to C2H2, and beyond, via successive hydrogenation, to C2H4 and C2H6.

These species are all less volatile than CH4 but more so than CO2. These less-volatile hydrocarbons are transported inwards by drifting ices to a closer snowline, shortening the lifetime of C-rich gas per Eq. (12) and limiting the maximum C/O somewhat. Nevertheless, we find that the picture of C/O evolution remains broadly similar to models that only include CO transformation, so we defer a more detailed discussion to Appendix B.2.

3.3 The role of dust traps

Dust traps have been connected to the depletion of C (e.g. McClure et al. 2020; Sturm et al. 2022) – as well as refractory elements such as Ca (Micolta et al. 2023, 2024) – in accreting material. Moreover, some works (Banzatti et al. 2020, 2023; Romero-Mirza et al. 2024a) relate the strength of MIR H2O emission to dust disc size – as a proxy for substructures at large R that hold back the dust – although other works (Gasman et al. 2025; Temmink et al. 2025) find less clear trends or distinctions. Thus, we build on other modelling works (Kalyaan et al. 2021, 2023; Easterwood et al. 2024; Mah et al. 2024; Sellek et al. 2025) and explore the impact of dust traps in the outer disc, which reduce the inward pebble flux and subsequent delivery of volatiles to the inner disc. We consider a model grid (using the hydrocarbon formation scheme) with four options for pressure bumps that trap dust: forming either early or late and at warm or cold temperatures (Table 1). Appendix D shows the results for all models. The results are similar for 1.0 M stars and all RC,0, so here we focus on VLMS models with RC,0 = 50 au, which is a) the only size where significant mass lies outside both of the trap locations considered and b) the least effective size at raising C/O in smooth discs, to see if traps can remedy this.

Figure 4 shows how, at a given time, any of the models with traps have a lower O/H than the corresponding model with no trap. The trap must therefore reduce the pebble flux to the inner disc and thus the supply of H2O (a major O-bearing molecule, especially in the presence of CO destruction, e.g. Fig. 3), which would be frozen out at both traps. Even the cold gap at 50 au, despite only being able to trap a small fraction of the dust mass, increases C/O to above 1. This trapping results in a similar flux of dust to the inner disc as in the more compact models without traps (which had simply lost their dust more quickly) and consequently the evolution of the C/O is very similar.

Conversely, warm dust traps at 5 au show a distinct behaviour, where they have significantly higher C/H than a model with no traps, despite a lower O/H; consequently they also sustain a high C/O. At the location of these traps, H2O, CH3 OH, CO2 and C2Hx hydrocarbons will be frozen out and trapped, while CH4 and CO are in the gas phase and can pass through normally. As the O/H is lowered, the ices cannot be leaking through significantly, yet the trapped material must be the source of additional C. Instead, the trap results in CH3 OH ice residing long enough in the disc (unlike in smooth discs, Appendix B.1) to be photodissociated to CH4 (by the secondary UV created by CRs)3. This then sublimates and passes through. We confirmed this by running extra models in which this reaction was switched off, and in which the sustained high C/H was no longer seen. Warm traps that form early enough to trap a significant mass of CH3 OH before most of the pebbles drift past 5 au, can reach C/O>10, while later trap formation leads to lower C/O∼ 3, though these values are upper limits due to the uncertain branching ratios for CH3 OH photodissociation (Öberg et al. 2009). This scenario involves a large change in the total (gas+ice) C/O across the dust trap, from >10 inside, to <1 outside (see Fig. D.1). This loosely resembles the model proposed for J1605 by Kanwar et al. (2025), where strong O depletion is needed in the hottest gas, but C/O falls again further out (in order to reproduce the CO2 fluxes), with the two reservoirs separated by a gas gap. Cold traps do not produce the same effect: any CH4 produced remains frozen out.

thumbnail Fig. 4

Same as for the right-hand panels of Fig. 2b but comparing models with and without dust traps for a RC,0 = 50 au disc around a 0.1 M star with an ionisation rate of ∼10−17 s−1.

thumbnail Fig. 5

Maximum inner-disc gas-phase C/O achieved for different parameters in smooth disc models with hydrocarbon formation around 1.0 M TTSs (left) or 0.1 M VLMSs (right). Darker green colours indicate more C-rich discs, while darker purple colours indicate more O-rich discs. We note that the colour bar is neither linear nor symmetric. The central circles show the maxima within the first 5 Myr, while the outer rings show the maxima across the full 10 Myr model run.

4 Discussion

4.1 Exploration of what is needed for higher C/O around VLMSs

Figure 5 summarises the maximum C/O achieved in the first 5 Myr and the full 10 Myr in the smooth disc models with hydrocarbon formation. Our models suggest that for VLMSs, CO transformation combined with dust and gas transport reliably produces C/O>1, as observed for most of their discs, for ζH2 ≳ 10−17 s−1 and RC,⊙ ≲ 16 au. Conversely, for more massive TTSs where C/O is usually <1, our models match this for the first 5 Myr if ζH2 ≲ 10−17 s−1 - though ζH2 ≲ 10−18 s−1 is required to keep discs below C/O=1 for the full 10 Myr – with weaker dependence on initial disc size. Figure 6 thus summarises the differences hypothesised for discs around VLMSs compared to the typical values measured for TTSs that would promote a higher C/O: a colder disc with close in snowlines, a more compact disc, and a higher ionisation rate (e.g. due to weaker CR exclusion). In the following, we discuss the viability of such conditions.

thumbnail Fig. 6

Summary of differences between VLMS discs and TTS discs in this work. For the purposes of this comparison, we assumed a TTS model with the typically measured ionisation rate of ~10−18 s−1 and an average RC,0 = 16 au. We assumed that the VLMS has a more compact disc of RC,0 = 5 au subject to an unexcluded galactic CR ionisation rate of ~10−17 s−1.

4.1.1 High ionisation rates

A key ingredient in our model is ionisation, which must be at a rate ≳10−17 s−1 for the VLMS discs to become C-rich and ≲10−17 s−1 for TTSs to remain O-rich for the first few megayears. ALMA observations have constrained ionisation in protoplane-tary discs (around M* ~ 0.5-2.1 M* stars) using the molecular ions HCO+, N2H+, and their isotopologues (Cleeves et al. 2015; Seifert et al. 2021; Aikawa et al. 2021; Long et al. 2024; Kashyap et al. 2024). These efforts have produced values ranging from 10−20−10−18 s−1, which are compared with the rates assumed in our 1 M models in Fig. A.1. While some works quote midplane rates and others incident rates, the observations are mostly sensitive to scales of 10s au where, with the ‘standard’ attenuation of CRs, the ionisation would not be reduced at the midplane compared to the incident rate.

As discussed in Section 2.2.3, there are several contributions to the ionisation in discs, including galactic CRs (which have the highest energy so most easily penetrate to the midplane), SLRs, and energetic particles (protons) and radiation (X-rays) produced by stellar magnetic activity. The low ionisation rates ≲10−18 s−1 inferred by the aforementioned works suggest that the galactic CR flux is reduced. CRs may be more strongly attenuated than assumed e.g. due to non-vertical propagation (Fujii & Kimura 2022; Long et al. 2024). Alternatively, the incident rate may be severely reduced (Cleeves et al. 2015) due to the exclusion of galactic CRs by the strong, highly magnetic stellar winds of young, rapidly rotating stars (Cleeves et al. 2013). However, while the observed rates meet our ≲10−17 s−1 criterion to avoid C-rich TTSs, they would be insufficient to drive the chemistry to C/O>1 around VLMSs. This suggests that for their discs to undergo sufficient chemical evolution, VLMSs must have higher ionisation rates, either due to negligible exclusion of CRs or a significant additional source of ionisation that compensates.

X-rays could be one such additional ionisation source. At the ages of interest, most stars ≲1.0 M are still on a Haysashi track and thus are fully convective (e.g. Sokal et al. 2020; Getman & Feigelson 2021). They are also relatively rapid rotators. Thus, both TTSs and VLMSs are highly magnetically active and produce strong X-rays. Both lie in the ‘saturated regime’ where LX/Lbol~10−3, (Preibisch et al. 2005), such that LX increases steeply with M (see also Güdel et al. 2007; Flaischlen et al. 2021). Assuming a passively heated irradiated disc, material at a given temperature – e.g. at a given snowline – is heated by a given stellar flux (Lbol/r2). Given similar LX/Lbol, the incident X-ray fluxes (LX/r2) should thus also be similar at a given snowline4. While flares are common and strong, their time-averaged luminosity is of the order of ∼10% of the background ‘characteristic’ level for all spectral types (Audard et al. 2000; Getman & Feigelson 2021; do Amaral et al. 2022) and so are unlikely to alter this trend. Moreover, X-rays seem unlikely to dominate ionisation or drive chemistry at the midplane (Cleeves et al. 2013, 2014; Rab et al. 2017). This is distinct from the disc surface, where VLMSs have a higher LX/LUV that changes the dominant chemical pathways (Walsh et al. 2015, we note that their M dwarf model used a roughly maximal LX for a 0.1 M star). Nevertheless, the X-ray-dominated layer may penetrate closer to the midplane in lower mass discs.

There are next to no direct constraints on stellar energetic particles, though it has been suggested that they could explain potentially higher ionisation rates close to DM Tau (Long et al. 2024). Young stars, along with other active stars and individual flares, follow a Güdel-Benz relation (Guedel & Benz 1993; Benz & Guedel 1994): a correlation between their X-ray and gyrosynchotron luminosities (Dzib et al. 2013, 2015; Pech et al. 2016). Gyrosynchotron radiation traces the acceleration of relativistic electrons by magnetic activity, with protons, which dominate the stellar energetic particle flux, likely also accelerated. This correlation could tentatively suggest that VLMSs, which have lower LX, would produce fewer energetic particles than earlier type stars, thus going against the trend we would need.

Finally, both solar-type and late M-dwarf stars are expected to have stellar winds that can exclude CRs (e.g. Rodgers-Lee et al. 2020b; Mesquita et al. 2022). The CR energies that are excluded at a distance r are those for which the advective timescale in the wind tadv = r/vwind is shorter than the diffusive timescale tdiff = r2, where the diffusivity κ ∝ 1/B, the magnetic field strength. Thus, if VLMSs launch slower or more weakly magnetic winds, they would exclude CRs less effectively. B declines at least as fast as r−1, making tadv /tdiff an increasing function of radius. Thus, CR exclusion should be weaker further from the star (Langner & Potgieter 2005; Mesquita et al. 2022), as suggested to explain the increase in IM Lup’s ionisation rate beyond 100 au (Seifert et al. 2021). Therefore, on ∼10 au scales, for example, VLMSs may have to make up less than the three orders of magnitude difference implied by comparing directly to the ζH2,1 au ~ 10−20 s−1 that Cleeves et al. (2013) estimated by extrapolating from the CR fluxes reaching Earth. Overall, it seems possible that VLMSs could exclude CRs less efficiently than TTSs but modelling of CR propagation as a function of distance around pre-MS stars of different masses is required to confirm if the exclusion is sufficiently negligible to reach ζH2 ~ 10−17 s−1 within the CO snowline.

4.1.2 Small initial disc sizes

Compact initial discs are needed around VLMSs to allow them to deplete enough dust to reach C/O significantly greater than 1 (unless ζH2 is high enough that TTSs would also host carbon-rich discs with 5 Myr). Thus, a dependence of disc size on stellar mass would help ensure that the VLMSs reach a consistently higher inner disc C/O than the more massive TTSs.

Theory predicts that discs may form more compact around lower mass stars (e.g. Hennebelle et al. 2016). This could place them towards the lower end of our adopted range, which was used for both M but based on a sample skewed towards higher masses (Trapman et al. 2023). Observations of embedded proto-stellar discs may hint in this direction (e.g. Yen et al. 2017) but at the Class II stage, only the few brightest (and hence largest) VLMS discs have been studied so far (Kurtovic et al. 2021). More measurements of gas disc sizes around VLMSs are needed to confirm if they are indeed systematically smaller.

4.1.3 Warm dust traps

Section 3.3 shows that CH3 OH photodissociation at a dust trap inside the CH4 snowline could provide a supply of additional C to the inner disc, potentially producing quite extreme C/O in discs that do not otherwise reach C/O>1. The known gaps on ≲5 au scales (Guerra-Alvarado et al. 2025) span spectral types K7 to M5 so include both TTS discs and VLMS discs. While this scenario could explain unusual or extreme cases such as DoAr 33 or J1605 respectively, we consider it an unlikely explanation for the C-rich – but less extreme – chemistry of most VLMS inner discs given that gap opening is not expected to be easier for lower M (Sinclair et al. 2020) and substructured discs indeed appear rarer around lower mass stars (in line with the occurrence rates of sufficiently massive exoplanets van der Marel & Mulders 2021). Thus it is more probable that such a scenario would occur around a TTS, rather than a VLMS.

4.2 The need for more inner disc C/O measurements

At present, there are few robust MIR measurements of inner disc C/H, O/H, and C/O; mostly what has been inferred is whether C/O<1 or C/O>1. More accurate values would help discriminate whether, for example, a ‘standard’ cosmic ray ionisation rate around VLMSs of ζH2,CR ~ 10−17 s−1 suffices or whether the C/O is extreme enough that their discs would need to receive enhanced fluxes of energetic particles to explain the observations. Likewise, it would be useful to know how close the discs that do not show C/O>1 come to doing so, and to know elemental abundances in order to understand if chemical transformation is nevertheless at work in some way. A reliable method must account for effects such as radial gradients in molecular abundances and variations in their respective emitting areas/locations; developing a reliable proxy in terms of retrieved column densities that averages over these effects would be useful.

Moreover, measurements across a wider range of ages would break some degeneracies. We have argued that a minimum ionisation rate of 10−17 s−1 is needed to explain the C-rich inner discs of most VLMSs while a maximum of 10−17 s−1 is needed to explain the O-rich inner discs of most TTSs. It is thus not clear whether we require higher ionisation rates around the VLMSs, as 10−17 s−1 could apply to all discs around stars of all masses. However, this value would lead to discs around TTSs becoming C-rich (i.e. reaching C/O>1) on timescales of ∼10 Myr, while lower ionisation rates would keep these disc below C/O=1. Studying discs around TTSs in older star-forming regions such as Upper Scorpius and determining whether they remain below C/O=1 would help distinguish between these two outcomes. JWST Cycle 2 GO programs 2970 (PI: Pascucci) and 3034 (PI: Zhang) should shed light on this matter; finding that the discs are C-rich would suggest that the lower C/O seen so far around TTSs is purely a question of the slower physical evolution of the gas and dust, whereas if the discs remain below C/O=1, that would require a slower chemical evolution driven by weaker ionisation.

4.3 The need for complementary observations with ALMA

As the models presented here involve the transport of chemically evolved gas from the outer disc to the inner disc, another way to test the origin of the high C/O could be to directly probe the outer disc chemical evolution around VLMSs and compare to that around TTSs using ALMA. Figure E.1 shows how our models predict that inward transport would cause the inner disc chemistry to lag behind that of the outer disc, and thus we expect C/Oin<C/Oout. So far, however, only a small survey of discs around five M4-M5 stars was conducted by Pegues et al. (2021), finding that the disc chemistry was fairly similar to TTSs, though with stronger C2H described as reminiscent of the hydrocarbon-rich inner discs.

More measurements of the chemistry of the outer discs of VLMSs are required to build on those initial results. For example, one could attempt to constrain C/O using CO and C2H (e.g. Miotello et al. 2019), C/H using CO and N2H+ (Anderson et al. 2022; Trapman et al. 2022) or [C 1] (Sturm et al. 2022), or the ionisation rate using N2H+ or HCO+ (Cleeves et al. 2015; Seifert et al. 2021; Aikawa et al. 2021; Long et al. 2024; Kashyap et al. 2024). These would measure, respectively, whether there is high C/O gas in the outer disc that can be transported to the inner disc, whether (chemical) depletion of CO is operating more or less efficiently around VLMSs, and whether there is sufficient ionisation to drive the chemistry needed to destroy CO. At present, there are no measurements of whether the ionisation rate is higher for VLMSs, as we suggest is required (Section 4.1.1).

5 Conclusions

Distinct differences have been identified between the chemistry of the inner regions of protoplanetary discs around VLMSs (≲0.3 M) and more massive TTSs (≳0.3 M). The former tend to have strong hydrocarbon emission, implying a more C-rich chemistry with gas-phase elemental C/O likely >1, while the latter are more O rich. Various suggestions have been put forward to explain these findings. For example, models for inner disc chemical evolution have previously predicted that the C/O should rise after the first ∼1 Myr (within which large amounts of H2O are delivered to the inner disc) and do so more rapidly for the colder discs around VLMSs (Mah et al. 2023). However, such models are typically limited to C/O∼1 (Fig. 2a) by having CO as the dominant C-carrying molecule. In this work, we explored such models further by demonstrating that in order to reach a C/O substantially greater than 1, the CO must be chemically depleted in such a way that the C is eventually liberated. The chemical transformation follows pathways commonly discussed in the context of underabundant CO (with respect to ISM levels) in the outer disc (Bosman et al. 2018; Krijt et al. 2020, Fig. 1). In the cold, dark disc midplane, such pathways are ultimately driven by the presence of ionising energetic particles and/or radiation.

  • The most critical reaction to allow C/O>1 is the dissociative ionisation of CO in the gas phase by He+, which directly liberates C, making it available for forming simple hydrocarbons (with O locked up in H2O ice). The less volatile the formed hydrocarbons are, the shorter the lifetime of C-rich gas in the disc. This initially raises C/O slightly faster, but limits the maximum C/O reached (Figs. 2b, B.2).

  • Whether a disc becomes C-rich and when depends on the relative physical and chemical timescales (Fig 5). In a disc subject to a low ionisation rate ≲10−18 s−1, the chemical timescales are longer than typical disc lifetimes, and C/O remains ≲1. For a very high ionisation rate, ≳10−16 s−1, the conversion of CO to hydrocarbons and the corresponding O depletion of the gas happens much more quickly than the timescale on which gas moves inwards through the disc such that any disc would reach C/O>1 within 1-3 Myr.

  • While, due to shorter physical timescales for their gas transport, discs around VLMSs increase in C/O more quickly and reach their maximum C/O sooner than those around TTSs, they nevertheless need ζH2 ≳ 10−17 s−1 to transform enough CO into hydrocarbons to reach C/O>1. Given that observations suggest ζH2 ≲ 10−18 s−1 around stars 0.5-2.0 M (Fig. A.1), inline with needing to be ≲10−17 s−1 to remain at C/O< for TTSs, VLMSs likely need to experience higher ionisation rates and thus faster chemical conversion, potentially requiring less-effective exclusion of CRs.

  • To reach a high C/O, it is crucial to ensure a very low resupply of H2O vapour at late times (Fig. 3). Discs around VLMSs therefore need to be at least as compact as the average TTS disc (≲16 au), since discs that are initially more compact lose their dust more quickly, or have effective trapping of dust and ices in pressure bumps in order to reach C/O>1.

  • The presence of a warm dust trap between the CH4 and CH3 OH snowlines allows the CH3 OH ice that formed outside the CO snowline to remain in the disc for long enough to be photodissociated by the secondary UV field produced by CRs (Fig. 4). Assuming that CH4 is the primary product of this process, it can sublimate and escape through the gap, leading to a sharp C/O gradient with a sustained high C/O inside the gap compared to outside (Fig. D.1). This could potentially explain discs with more extreme C/O, such as J1605.

We conclude that CO destruction and the subsequent transport of O-depleted material from the outer disc is a viable way to produce the higher C/O seen in the inner disc around VLMSs, provided that the ionisation rates are ≳10−17 s−1. Constraints on inner disc chemistry for older discs should demonstrate whether discs around TTSs also reach these conditions – and the differences are just a matter of their longer physical timescales – or whether the chemical transformation is also faster in the VLMS discs. More observations of VLMSs with ALMA could help confirm whether they have suitable disc sizes, ionisation rates, and chemistries for our proposed route to C/O>1.

Acknowledgements

We thank the anonymous reviewer for helpful suggestions that improved the clarity of the work. We also thank L. Trapman and A. Bosman for sharing details of the CO chemical network, C. Walsh for discussions about the grain surface processes and the role of X-rays in chemistry, J. Kanwar & I. Kamp for insights about determination of C/O from spectra, N. Kurtovic for discussions about ALMA observations, and D. Rodgers-Lee & L. I. Cleeves for discussions about ionisation processes. A.D.S. and E.v.D. acknowledge support from the ERC grant 101019751 MOLDISK. E.v.D. also acknowledges support from the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150) and grant TOP-1 614.001.751 from the Dutch Research Council (NWO).

Appendix A Additional details of chemical model

A.1 Molecular data and tunnelling barriers

Table A.1

Binding energies, prefactors, and initial abundances (with respect to H) assumed for the species considered in this work and the resulting snowlines around 0.1 and 1.0 M stars (defined as the locations where the gas and ice abundances are equal at t = 0).

Table A.2

Tunnelling barriers assumed for grain-surface reactions.

A.2 Equilibrium processes for reactive intermediaries

For each of the reactive intermediaries that react with CO - He+, H, and OH - we assume equilibrium abundances between a formation rate driven directly by cosmic ray ionisation (CR) or their secondary UV field (CRPHOT) and all major destruction routes, following those in (Krijt et al. 2020): He + CR  He+He+ + {CO, N2, H2, dust}  He + {CO+, N2++, H2+, dust+},H2O(s) + CRPHOT  OH(s) + H(s)OH(s) + {CO(s), H(s)}  {CO2(s) + H(s), H2O(s)},H2+CR2.4HHH(s)H(s) + H2S(s) 3H2S2(s).\ce{He + CR &-> He+} \\ \ce{He+ + \{CO,N2,H2,dust\} &-> He + \{CO+,N2+,H2+,dust+\}} ,

Here, H2S provides a baseline sink of H atoms on grain surfaces (we do not follow S-bearing species explicitly so assume a constant abundance of 2 × 10−8) and the numbers above the reactions scale the rates of the reactions shown to account for other pathways following Krijt et al. (2020).

A.3 Ionisation

A key input of our model is the ionisation, which creates the intermediaries that reach with CO (and other species). Figure A.1 shows the ionisation rate profile at the midplane for our models: for a fixed incident ionisation rate at the disc surface (10−17 s−1 is shown here), the midplane rate varies as a function of time as the disc becomes less optically thick and the attenuation reduces - for comparison, the bottom panel shows the evolving surface densities. Observational estimates (Cleeves et al. 2015; Seifert et al. 2021; Aikawa et al. 2021; Long et al. 2024) are included for the TTSs. In addition, we include the median and range of the ionisation rate posteriors derived from CO and N2H+ fluxes when kinematically derived masses are used as priors, as derived in the exoALMA program by (Trapman et al. 2025).

thumbnail Fig. A.1

Top: Midplane ionisation rate profile at t=0, 1, and 5 Myr (increasingly dark blues, evolution indicated by arrow) for fiducial models around TTS (left) and VLMSs (right). The TTS models are compared to observationally estimated values for individual sources (red and yellow colours) and the range and median for the exoALMA program (Trapman et al. 2025). The black horizontal line on each plot represents the floor due to SLRs and the vertical dotted line shows the location of the CO snowline. Bottom: Surface densities of the gas at the times indicated in the top panels. The horizontal line shows the assumed 96 g cm−2 penetration depth of CRs.

Appendix B Role of different processes

B.1 Impact of each molecular reaction

To test which reactions matter most, we added them into our model one by one, for a fiducial case with RC = 16 au, MD/M* = 0.05, and ζH2 = 1.3 × 10−17 s−1. The results are shown in Fig. B.1 and compared to a reference model where no reactions were included and C/O reaches a maximum of 0.98 (black). The maximum C/O reached in each case is tabulated in Table B.1. For the TTS models, the C/O always rises monotonically over time after the first Myr so a greater maximum is reached after 10 Myr than within 5 Myr; the VLMS models typically peak just after 5 Myr so behave likewise.

The first reaction considered is the conversion of CO to CH4 (red). We find that this is sufficient for discs around either stellar mass to cross C/O=1 within 10 Myr, although only the discs around VLMSs cross within the first 5 Myr. This reaction is key to achieving high C/O levels.

Next we add hydrogenation of CO to CH3OH (dark green). While this does not itself separate C and O from one another, it does help flush the most volatile reservoir of O out of the disc more quickly, as shown by the model tracks moving up and to the right in the region where C/O<1. This prevents much of the delivery of CO to its snowline and thus reduces the amount available to buffer the C/O ratio to 1, enhancing the effectiveness with which the conversion of gas phase CO to CH4 raises C/O. We also tested a model where this is the only reaction (orange), confirming that while it helps to accelerate the reduction in C/H and O/H, it never leads to C/O>1 on its own.

The third reaction we add is the photodissociation of the CH3 OH ice forming CH4 (light green). While this process can play the crucial role of breaking the CO bond and separating out C and O into species with different volatilities, in practice it makes barely any discernible difference to the evolution of our models, raising C/O by only a percent or two. This is because the timescale for these reactions is on the order of Myr, while the majority of the ice has a much shorter lifetime in smooth discs due to efficient radial drift.

The final route for CO transformation is the reaction with OH radicals forming CO2 (light blue). We find this mostly has no impact on our results as a very limit range of radii at the midplane are within the right temperature range for this to dominate over the other gas or ice phase routes.

In addition, we add three reactions forming the larger hydrocarbons C2H2, C2H4, and C2H6 (dark blue). At intermediate times, they slightly accelerate the delivery of C-bearing molecules to the inner disc, raising C/O. Conversely, at late times they have very small impact in the direction of lowering the C/O ratio, as they convert CH4 into less volatile species, thus reducing the lifetime of the carbonaceous component in the disc and preventing such high C/O from building up.

B.2 Impact of hydrocarbon formation

For ionisation rates < 10−17 s−1, the evolution of C/O (likewise C/H and O/H) is very similar when we include larger hydrocarbons (compare Figs. B.2 and 2b). This is because, due to the long chemical timescales, large reservoirs of C & O remain in CO and are unaffected by the addition of a route out of CH4. The fiducial ionisation rates are more visibly affected, but the conclusions remain the same. The timescale for C2H2 production from CH4 in our model is ~ 3 Myr - comparable to the accretion timescale at the CH4 snowline for the TTS model and longer than that for the VLMS model - so there is hardly any time for a CH4 molecule once formed to be converted to larger hydrocarbons.

thumbnail Fig. B.1

Same as in Fig. 2a but for disc models with different numbers or types of reactions included.

Table B.1

Maximum C/O achieved within the time range specified for models including different numbers or types of reactions for a fiducial case with RC = 16au, MD/M* = 0.05, and ζH2 = 1.3 × 10−17s−1.

A much larger difference emerges for high ionisation rates. Previously, CH4 was the dominant molecule but now it is efficiently converted to C2H2 , which is much less volatile. This results in C moving from the gas phase into ices and being transported inwards much faster by drift such that it is ultimately removed from the disc much sooner. This prevents such a large gaseous C/O ratio from building up, though values ~ 3 are still achieved for both M*. The formation of less volatile hydrocarbons still allows high ionisation rate models to explain DoAr 33, but they would not be able to reach the extreme levels of O depletion inferred for J1605 by Kanwar et al. (2025) without also having considerable depletion of C and thus too low C/O. We note that trapping of CH4 in mixed ices would have a similar effect as it would desorb at high temperatures (Collings et al. 2004), but as much of the CH4 produced in our models is inside of its snowline, not all CH4 would be affected by this.

For VLMSs, the endpoint after 10 Myr is similar, as the gas accretion of CH4 onto the star around cooler stars (Eq. 12) is fast enough to have removed it anyway, regardless of conversion into less volatile species. Conversely, for TTSs, loss of CH4 via accretion is insufficiently fast, so C2H2 formation lowers the C/O - and underlying C/H - by over an order of magnitude. Such C/H values are comparable to the 1-2 orders of magnitude depletions measured from NIR atomic lines (tracing accreting gas inside the sublimation radius McClure 2019; McClure et al. 2020; Sturm et al. 2022), but only occur in our models on timescales longer than the ages of most of those sources.

thumbnail Fig. B.2

Same as in the lower panels of Fig. 2b but for disc models that include formation of C2Hx hydrocarbons.

B.2.1 Molecular ratios involving C2H2

Figure B.2 focuses on the inner-disc gas-phase C/O ratio, rather than individual molecules, as significant chemical reprocessing may happen in the observable upper layers, with the complexity of the products depending on the relative roles of UV or X-ray radiation (Walsh et al. 2015; Raul et al. 2025). For reference, however, B.3 shows the evolution of averaged CO2/C2H2 and C2H6/C2H2 column density ratios for these models.

The column density ratios were approximated as the ratio of area-averaged observable column densities: NxNy=Nx,obsNy,obs,\frac{N_x}{N_y} = \frac{\langle N_{x,\rm obs}\rangle}{\langle N_{y,\rm obs}\rangle},(B.1)

where Nx,obs=0RsnowNxfobs2πdRπRsnow2.\langle N_{x,\rm obs}\rangle = \frac{\int_0^{R_{\rm snow}} N_x f_{\rm obs} 2\pi{\rm d}R}{\pi R_{\rm snow}^2} .(B.2)

Over time, as more C2H2 is created and as the less-volatile CO2 advects onto the star on a shorter timescale than C2H2 , the CO2/C2H2 falls. This is in line with the tentative observational trends, where Arabhavi et al. (2025a) find that VLMSs in the older Upper Sco star forming region have lower CO2/C2H2 flux ratios. Moreover, Long et al. (2025) note an average increase in the NC2H2 /NCO2 column density ratio over time - with the extremes being the young H2O-rich VLMS Sz 114 and the 34 Myr old J0446B - that they attributed to an increase in C/O. The increase we see is not driven by the fact that equilibrium chemistry forms more C2H2 and less CO2 at higher C/O as in thermochemical models, but simply as a result of the secular evolution of the disc as the molecules are formed and transported. At late times, the hydrocarbon-rich gas starts to accrete onto the star to be replaced with C-depleted gas, and the CO2/C2H2 ratio rises again.

The C2H6/C2H2 ratio is strongly dependent on the ionisation rate, and rather weakly on the disc size. The values are similar for both stellar masses at a given ionisation rate and in all cases are very high, as hydrogenation is very efficient and as in our simplified approach, C2H6 acts as a sink. Photodissociation reactions of C2H6 (and C2H4) could help reverse this balance as could the formation of more complex hydrocarbons or organic molecules via other reactions that were not considered.

B.3 Impact of carbonaceous grain erosion

Carbonaceous grains can be eroded by various processes, as outlined in the introduction. While most of these processes need high temperatures to become efficient and are thus restricted to the inner disc, previous modelling has suggested that the destruction of carbonaceous grains by UV at large distances from the star could be responsible for the high C/O observed there in some discs (Bergin et al. 2016; Anderson et al. 2017; Bosman et al. 2021). As the focus of this work is the connection between the outer disc chemical processing and the inner disc, we focus on whether this ‘UV photolysis’ would affect our results, or if grain erosion can only raise the inner-disc C/O if it happens in situ.

We follow the approach used by Anderson et al. (2017); Bosman et al. (2021) who based their rates of destruction by UV on measurements from photolysis experiments of hydrogenated amorphous carbon grains, in which CH4 was the primary product, expressed in terms of the UV photon flux FUV = 108G0 and the yield of C atoms per photon YUV = 8 × 10−4 (Alata et al. 2014, 2015): kphot=σYUVFUVNC,k_{\rm phot}=\frac{\sigma Y_{\rm UV} F_{\rm UV}}{N_{\rm C}} ,(B.3)

where σ = πa2 is the geometric cross-section of the grains, and NC = (4πρsa3)/(3µCmH) is the number of carbon atoms per grain where µC = 12 is the relative atomic mass of C and ρs = 2.24 g cm−3 is the bulk density of carbonaceous material. These factors combine to give (3µCmH)/(4aρs) showing that even if the large grains are exposed to the UV (e.g. in the outer disc), their low surface area-to-volume ratio means they are inefficiently eroded. We therefore consider only small grains (assumed to be 0.1 µm) for our calculation.

thumbnail Fig. B.3

Evolution of molecular gas-phase column density ratios averaged across the disc inside their respective snowlines. Top: C2H6 and C2H2. Bottom: CO2 and C2H2.

The disc evolution code is a 1D, vertically integrated, code, that tracks the mass fractions ϵi=Σi/Σ$\epsilon_i = \Sigma_i/\Sigma$. We express the rate of change of the carbonaceous grain mass fraction as dϵidt=1ΣdΣidt\begin{align} \frac{{\rm d}\epsilon_i}{{\rm d}t} &= \frac{1}{\Sigma} \frac{{\rm d} \Sigma_i}{{\rm d}t} \\ \end{align}(B.4) =2Σ0dρi(z)dtdz&=\frac{2}{\Sigma} \int_0^\infty \frac{{\rm d} \rho_i(z)}{{\rm d}t} {\rm d}z \\(B.5) =2ΣσYUVNC0ρi(z)FUV,zdz,&=\frac{2}{\Sigma} \frac{\sigma Y_{\rm UV}}{N_{\rm C}} \int_0^\infty \rho_i(z) F_{\rm UV, z} {\rm d}z ,(B.6)

We assume the UV is attenuated by the vertical column of gas as FUV,z=FUV,0eτUV,z$F_{\rm UV, z}=F_{\rm UV,0}e^{-\tau_{\rm UV, z}}$ where the optical depth to UV τUV,z=σH(ϵd/0.01)znH(z)dz$\tau_{\rm UV,z} = \sigma_H (\epsilon_d/0.01) \int_z^\infty{n_H(z') {\rm d}z'}$ for UV cross section per hydrogen atom σH = 2.6 × 10−21 cm−2 (Facchini et al. 2016). By writing μimHni=ρi=ϵiρ=ϵiμmHn$\mu_i m_H n_i = \rho_i = \epsilon_i \rho = \epsilon_i \mu m_H n$ assuming the composition is vertically uniform, and taking nH = 2nH2, we have τUV,z=μσH(ϵd/0.01)ϵH2zn(z)dz=σeffNz$\tau_{\rm UV,z} = \mu \sigma_H (\epsilon_d/0.01) \epsilon_{\rm H_2} \int_z^\infty{n(z') {\rm d}z'} = \sigma_{\rm eff} N_z$ where σeff=μσH(ϵd/0.01)ϵH2$\sigma_{\rm eff} = \mu \sigma_H (\epsilon_d/0.01) \epsilon_{\rm H_2}$ and Nz=zn(z)dz$N_z = \int_z^\infty{n(z') {\rm d}z'}$. Noting that n(z)dz = −dNz, we can then write dϵidt=2ϵiμmHΣσYUVNCz=z=0FUV,0eσeffNzdNz\begin{align} \frac{{\rm d}\epsilon_i}{{\rm d}t} &= \frac{2\epsilon_i \mu m_H}{\Sigma} \frac{\sigma Y_{\rm UV}}{N_{\rm C}} \int_{z=\infty}^{z=0} F_{\rm UV, 0} e^{-\sigma_{\rm eff} N_z} {\rm d}N_z \\ \end{align}(B.7) =ϵiσYUVNC1eτUV,0τUV,0,&=\epsilon_i \frac{\sigma Y_{\rm UV}}{N_{\rm C}} \frac{1-e^{-\tau_{\rm UV, 0}}}{\tau_{\rm UV, 0}} ,(B.8)

i.e. the vertically averaged rate is modulated by a factor (1 -e−τUV,0)/τUV,0. By using the vertically averaged abundance and rate, we assume that the material is well-mixed vertically, which is reasonable for our purposes as the timescale for vertical mixing is shorter than the timescale for accretion from the outer disc to the inner disc by a factor h2. This approach is in contrast to the works above, which are static thermochemical models with a vertical C/O gradient that only rises above 1 in the surface layers that are optically thin to the UV. As any grain effectively spends only (1 - e−τUV,0)/τUV,0 of its time in the optically thin layers due to vertical mixing, the impact on raising C/O is essentially diluted by this factor under our assumptions. We note that this ‘unrestricted photolysis regime’ (Klarmann et al. 2018; Vaikundaraman et al. 2025) is a best case scenario: the average erosion may be less efficient if the residence time of individual grains in the upper layers is longer than their destruction timescale 1/kphot, such that erosion becomes limited by the rate of vertical mixing of fresh grains to the surface (Anderson et al. 2017; Klarmann et al. 2018). In order to fully deplete the carbon from the grains, efficient collisions are also required to recycle material from larger grains into smaller grains that are able to reach the exposed surface (Vaikundaraman et al. 2025).

thumbnail Fig. B.4

Comparison of the timescales at t = 0 yr around a 1.0 M star (left) and 0.1 M star (right) for the erosion of carbonaceous grains by UV photolysis (blue, solid) with the timescales of gas accretion (orange, dotted; see Equation 12), dust grain growth (green, dashed) and radial drift (red, dashed). The carbon depletion timescale at 1 au in the ‘fiducialfc025’ model from Vaikundaraman et al. (2025) is shown as the blue V.

As we argue that photolysis is only relevant for the small dust grains, and as these are the only grains that can reside in the upper layers where the UV can penetrate, we scale the rate with the fraction of dust in small grains which, in the two population model, is fsmall = 0.03 - 0.25 of the total grain mass (in the drift limit and fragmentation limit respectively Birnstiel et al. 2012). This factor is assumed to modify both the optical depth to UV τUV,0 and the final rate dϵi/dt${\rm d}\epsilon_i/{\rm d}t$, though in the limit of large τUV,0, the two contributions will cancel, as, essentially, the UV penetrates (and can thus erode) a fixed number of grains

The incident UV field FUV,0 is given by a 12000, K blackbody with luminosity Lfuv,0 equal to the accretion luminosity: LFUV,0=fUV(12000K)0.8GMM˙R,L_{\rm FUV, 0} = f_{\rm UV}(12000\,\mathrm{K}) \frac{0.8 G M_* \dot{M}}{R_*} ,(B.9)

where fFUV(12000 K) = 0.245 is the fraction of a 12000 K blackbody’s luminosity radiated at FUV energies. The UV radiation (in Habing units) intercepted per unit area of the disc is then G0=Lacc1.63×103ergs1cm2ϕ4πr2,G_0 = \frac{L_{\rm acc}}{1.63\times10^{-3}\,\mathrm{erg\,s^{-1}\,cm^{-2}}} \frac{\phi}{4\pi r^2} ,(B.10)

where ϕ~h is the flaring angle of the disc surface.

The model is also included in Fig. B.1 as the model involving C-grains (purple) and in the final row of Table B.1. A very small (sub-percentage) increase is seen for the VLMS model. This can be understood by comparing the grain destruction timescale to the timescales of the dust and gas evolution (Fig. B.4). Although the timescales for grain destruction are shorter than the typical disc lifetime, in these models, over much of the disc, the dust can grow to large grains, drift, and accrete onto the star faster than this timescale. Thus a very small fraction of the dust is eroded before it is accreted onto the star, and so although there was initially a large reservoir of refractory grains, it is not very effectively liberated for either M*. This is in agreement with the results of Klarmann et al. (2018). We note that our UV photolysis timescale estimate is somewhat longer, but of the same order of magnitude as, that of Vaikundaraman et al. (2025): the differences could be due to differences in UV luminosity and scattering or grain size distributions, compositions, and structure, but nevertheless our conclusions would still hold. Models with less efficient growth and drift might allow a larger fraction of the carbonaceous grains to be destroyed and to contribute to enriching the gas, for example if they are retained at a dust trap (similar to the scenario for CH4 production from CH3OH in Section 3.3).

Appendix C Observational outliers

So far, two discs observed by MIRI-MRS stand out as outliers to the trend of C-rich VLMS discs and O-rich TTS discs. DoAr 33 has a C-rich inner disc with C/O = 2 - 4 (Colmenares et al. 2024) around a 1.1 M star in the 0.3 - 6Myr ρ Ophiuchi region. Explaining this source, on these timescales, via transport of outer disc material would require an ionisation rate elevated by a factor ~ 10 above the galactic rate. This depletes O sufficiently, while C remains somewhat enhanced inside the CO snowline following sublimation of drifting ices, to reproduce the inferred abundance for this disc (without additional mechanisms e.g. erosion of carbonaceous grains). Conversely, Sz 114, a H2O-rich disc around a VLMS Xie et al. (2023), has an estimated RC = 18.2 au (Trapman et al. 2023), making it larger than average, consistent with its status as one of the brightest discs at millimetre wavelengths around a mid-to-late M-type star in the Lupus star-forming region (with millimetre flux being correlated with millimetre disc size, Tripathi et al. 2017; Andrews et al. 2018; Kurtovic et al. 2021). The large size of Sz 114 could play a role in its relatively C-poor chemistry, as explored in Section 3.2.2. Alternatively, Sz 114 may simply me too young to have yet crossed to C/O> 1 (Xie et al. 2023; Long et al. 2025). In our models, this could happen if it is ≲ 2 Myr, consistent with Lupus’ age range. Finally, the H2 O-rich phase could be prolonged past that time if the disc hosts leakier traps (e.g. Mah et al. 2024) than those we explored in Section 3.3 & Appendix D and so long as gas-phase CH4 is not produced there (e.g. if the traps are cold enough for it to remain frozen out, or the photolysis of CH3 OH is less effective than assumed).

Appendix D Further examples of models with traps

thumbnail Fig. D.1

C/O as a function of R at 2.5 Myr for models with no trap (black) or a warm trap (red) in a 50 au disc around a 0.1 M star with an ionisation rate of 1.3 × 10−17 s−1. Vertical lines mark the location of the minimum gas density (’Gap location’) and the pressure maximum where dust accumulates (’Dust ring’).

Figure D.1 shows the C/O ratio as a function of radius for the warm trap model in a 50 au disc around a 0.1 Mo star as explored in Section 3.3. Since the C/O gradient is similar to the case inferred for J1605, we plot this model at 2.5 Myr, the dynamical traceback age of the β Sco subregion of Upper Sco where J1605 is located (Miret-Roig et al. 2022).

Figure D.2 shows equivalent warm and cold-trap models to Fig. 4 for 1.0 M stars and Fig. D.3 shows warm-trap models for all initial radii. A similar effect is seen for all RC,0 investigated, for example with warm traps photodissociating CH3 OH and leaking CH4 into the inner disc. Most notably, more material starts outside the trap for larger discs, and therefore in these cases, there is a larger trapped CH3 OH reservoir that can sustain a higher C/H and C/O.

thumbnail Fig. D.2

Same as in Fig. 4 but for 1 M.

thumbnail Fig. D.3

Same as in Fig. 4 and D.2 but comparing models with and without warm dust traps for a range of initial disc radii.

Appendix E Comparing the inner and outer disc

Figure E.1 compares the inner and outer disc chemistry for the model grid with hydrocarbon formation. As before, the inner disc is measured at the inner boundary of the model in order to represent the hottest gas in the disc, since most lines are observed with excitation conditions of several 100s K. For the outer disc, the choice of is more arbitrary (as we are mostly looking to demonstrate the lag between the two), but since we are interested in the volatile C and O abundances, we needed to measure at a location where both of the main volatile C carriers - CO and CH4 - are in the gas phase and thus chose just inside the CH4 snowline. At this location, the C/H initially rises a little as CO and CH4 ices desorb at their snowlines but eventually starts to deplete by up to two orders of magnitude (consistent with typical TTS disc measurements Favre et al. 2013; Kama et al. 2016b; McClure et al. 2016; Miotello et al. 2017; Anderson et al. 2022; Trapman et al. 2022; Paneque-Carreño et al. 2025) in line with the gas outside the snowlines.

thumbnail Fig. E.1

Relationship between the inner disc C/O (measured at the 0.1 au inner boundary of the models) and the outer disc (measured at 90% of the CH4 snowline) C abundance (left) and C/O (right) for discs around VLMSs.

More interesting is the relationship between inner and outer disc C/O. Initially, we measure C/O< 1 in both the inner and outer disc, due to the presence of O2 gas in the outer disc. The inner disc C/O initially drops quite rapidly as O-rich H2O ice is deposited in the inner disc by drift. Conversely, the enhancement at the CH4 and CO snowlines takes slightly longer to occur, as the dust grains take longer to grow and start to drift further out in the disc, and as a larger fraction of the molecules are already in the gas-phase any potential enhancement is weaker. Thus, the model tracks initially drop almost vertically on the plot of inner disc vs outer disc C/O (Fig. E.1). As the O-rich gas starts to accrete, they reverse this course, the timescale for this accretion is very short due to the proximity of the H2O snowline to the star, so the rise in inner disc C/O happens faster than that in the outer disc. Nevertheless, the discs do cross to C/O in the outer disc before they can do so in the inner disc, as once the remnants of the initial drift phase are accreted, the chemistry is driven by the inward spreading of the outer-disc gas. Moreover, the inner disc C/O is kept somewhat lower than around the CH4 snowline due to the sublimation of O-rich ices between those locations, especially H2O (Fig. 3). Thus, we see all the model tracks pass through the lower right quadrant of Fig. E.1. Eventually, the inner disc may also reach C/O> 1 and then the discs move into the upper right of the plot, but remain below the line where the inner and outer disc C/O are equal, though discs with sufficient resupply of H2O may just remain in the bottom right (e.g. the lightest green track).

We note that if C/O is instead reached by the release of additional C in the inner disc, without O-depletion in the outer disc, we might expect discs to populate the upper left hand quadrant, where our models are never seen. This may be a way to distinguish between the two proposed scenarios of high C/O due to inner disc grain erosion vs O-depleted gas. For example, in the low-turbulence (α = 10−4) models of Houge et al. (2025), a strong, early pebble flux causes the irreversible thermal decomposition of refractory organics into C2H2 on drift timescales (< 0.1 Myr), producing an early high C/O in the inner disc, while the composition of the outer disc has barely evolved (chemical reactions of CO are neglected in their work). On longer timescales, however, the diffusion of this C-rich gas as far as the C2H2 snowline raises C/O above 1 outside the H2O snowline and creates a longer-lived volatile C reservoir that can resupply the inner disc on longer timescales once the initial C-rich reservoir accretes onto the star. Thus, we might expect these models to approach similar conditions to ours at late times, but future work should compare the two directly. Finally, we also caution that Fig. D.1 shows an alternative way to achieve a higher C/O in the inner disc than the outer disc, but we do find that in the plane shown in Fig. E.1, our models stay in the upper right as the CH4 gas produced can also diffuse outwards towards its snowline.

References

  1. Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. M. 1999, ApJ, 519, 705 [Google Scholar]
  2. Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, ApJS, 257, 13 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alata, I., Cruz-Diaz, G. A., Muñoz Caro, G. M., & Dartois, E. 2014, A&A, 569, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Alata, I., Jallat, A., Gavilan, L., et al. 2015, A&A, 584, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Altwegg, K., Balsiger, H., & Fuselier, S. A. 2019, ARA&A, 57, 113 [NASA ADS] [CrossRef] [Google Scholar]
  6. Anderson, D. E., Bergin, E. A., Blake, G. A., et al. 2017, ApJ, 845, 13 [NASA ADS] [CrossRef] [Google Scholar]
  7. Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927, 229 [NASA ADS] [CrossRef] [Google Scholar]
  8. Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157 [Google Scholar]
  9. Arabhavi, A. M., Kamp, I., Henning, T., et al. 2024, Science, 384, 1086 [NASA ADS] [CrossRef] [Google Scholar]
  10. Arabhavi, A. M., Kamp, I., Henning, T., et al. 2025a, A&A, 699, A194 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Arabhavi, A. M., Kamp, I., van Dishoeck, E. F., et al. 2025b, ApJ, 984, L62 [Google Scholar]
  12. Arulanantham, N., Salyk, C., Pontoppidan, K., et al. 2025, AJ, 170, 67 [Google Scholar]
  13. Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396 [Google Scholar]
  14. Bae, J., Isella, A., Zhu, Z., et al. 2023, ASP Conf. Ser., 534, 423 [NASA ADS] [Google Scholar]
  15. Banzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, ApJ, 903, 124 [NASA ADS] [CrossRef] [Google Scholar]
  16. Banzatti, A., Pontoppidan, K. M., Carr, J. S., et al. 2023, ApJ, 957, L22 [NASA ADS] [CrossRef] [Google Scholar]
  17. Banzatti, A., Salyk, C., Pontoppidan, K. M., et al. 2025, AJ, 169, 165 [Google Scholar]
  18. Behmard, A., Fayolle, E. C., Graninger, D. M., et al. 2019, ApJ, 875, 73 [Google Scholar]
  19. Benz, A. O., & Guedel, M. 1994, A&A, 285, 621 [NASA ADS] [Google Scholar]
  20. Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101 [Google Scholar]
  21. Bergin, E. A., Booth, R. A., Colmenares, M. J., & Ilee, J. D. 2024, ApJ, 969, L21 [NASA ADS] [CrossRef] [Google Scholar]
  22. Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25 [Google Scholar]
  23. Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020, ApJ, 898, 97 [NASA ADS] [CrossRef] [Google Scholar]
  24. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541 [Google Scholar]
  26. Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994 [Google Scholar]
  27. Borderies, A., Commerçon, B., & Bourdon, B. 2025, A&A, 694, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Bosman, A. D., Alarcón, F., Zhang, K., & Bergin, E. A. 2021, ApJ, 910, 3 [NASA ADS] [CrossRef] [Google Scholar]
  30. Brunn, V., Marcowith, A., Sauty, C., et al. 2023, MNRAS, 519, 5673 [NASA ADS] [CrossRef] [Google Scholar]
  31. Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5 [Google Scholar]
  32. Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123 [NASA ADS] [CrossRef] [Google Scholar]
  33. Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204 [CrossRef] [Google Scholar]
  34. Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133 [NASA ADS] [CrossRef] [Google Scholar]
  35. Colmenares, M. J., Bergin, E. A., Salyk, C., et al. 2024, ApJ, 977, 173 [NASA ADS] [CrossRef] [Google Scholar]
  36. Colombo, S., Orlando, S., Peres, G., et al. 2019, A&A, 624, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Dipierro, G., Laibe, G., Alexander, R., & Hutchison, M. 2018, MNRAS, 479, 4187 [NASA ADS] [CrossRef] [Google Scholar]
  38. do Amaral, L. N. R., Barnes, R., Segura, A., & Luger, R. 2022, ApJ, 928, 12 [NASA ADS] [CrossRef] [Google Scholar]
  39. Drozdovskaya, M. N., Walsh, C., Visser, R., Harsono, D., & van Dishoeck, E. F. 2014, MNRAS, 445, 913 [NASA ADS] [CrossRef] [Google Scholar]
  40. Dzib, S. A., Loinard, L., Mioduszewski, A. J., et al. 2013, ApJ, 775, 63 [Google Scholar]
  41. Dzib, S. A., Loinard, L., Rodríguez, L. F., et al. 2015, ApJ, 801, 91 [Google Scholar]
  42. Easterwood, W., Kalyaan, A., & Banzatti, A. 2024, ApJ, 977, 21 [Google Scholar]
  43. Eistrup, C. 2023, ACS Earth Space Chem., 7, 260 [NASA ADS] [CrossRef] [Google Scholar]
  44. Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Facchini, S., Clarke, C. J., & Bisbas, T. G. 2016, MNRAS, 457, 3593 [NASA ADS] [CrossRef] [Google Scholar]
  47. Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38 [Google Scholar]
  48. Feigelson, E. D., Garmire, G. P., & Pravdo, S. H. 2002, ApJ, 572, 335 [Google Scholar]
  49. Flaischlen, S., Preibisch, T., Manara, C. F., & Ercolano, B. 2021, A&A, 648, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Fujii, Y. I., & Kimura, S. S. 2022, ApJ, 937, L37 [Google Scholar]
  51. Gail, H.-P., & Trieloff, M. 2017, A&A, 606, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Gasman, D., van Dishoeck, E. F., Grant, S. L., et al. 2023, A&A, 679, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Gasman, D., Temmink, M., van Dishoeck, E. F., et al. 2025, A&A, 694, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Getman, K. V., & Feigelson, E. D. 2021, ApJ, 916, 32 [NASA ADS] [CrossRef] [Google Scholar]
  55. Grant, S. L., Kurtovic, N. T., van Dishoeck, E. F., et al. 2024, A&A, 689, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Grant, S. L., Temmink, M., van Dishoeck, E. F., et al. 2025, A&A, in press https://doi.org/10.1051/0004-6361/202555862 [Google Scholar]
  57. Grant, S. L., van Dishoeck, E. F., Tabone, B., et al. 2023, ApJ, 947, L6 [NASA ADS] [CrossRef] [Google Scholar]
  58. Güdel, M., Briggs, K. R., Arzner, K., et al. 2007, A&A, 468, 353 [Google Scholar]
  59. Guedel, M., & Benz, A. O. 1993, ApJ, 405, L63 [NASA ADS] [CrossRef] [Google Scholar]
  60. Guerra-Alvarado, O. M., van der Marel, N., Williams, J. P., et al. 2025, A&A, 696, A232 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Harsono, D., Jørgensen, J. K., van Dishoeck, E. F., et al. 2014, A&A, 562, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [Google Scholar]
  63. Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8 [Google Scholar]
  64. Henning, T., Kamp, I., Samland, M., et al. 2024, PASP, 136, 054302 [NASA ADS] [CrossRef] [Google Scholar]
  65. Houge, A., Johansen, A., Bergin, E., et al. 2025, A&A, 699, A227 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Hsieh, C.-H., Arce, H. G., Maureira, M. J., et al. 2024, ApJ, 973, 138 [NASA ADS] [CrossRef] [Google Scholar]
  67. Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96 [EDP Sciences] [Google Scholar]
  68. Kalyaan, A., Pinilla, P., Krijt, S., Mulders, G. D., & Banzatti, A. 2021, ApJ, 921, 84 [CrossRef] [Google Scholar]
  69. Kalyaan, A., Pinilla, P., Krijt, S., et al. 2023, ApJ, 954, 66 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kama, M., Bruderer, S., Carney, M., et al. 2016a, A&A, 588, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016b, A&A, 592, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Kanwar, J., Kamp, I., Jang, H., et al. 2024, A&A, 689, A231 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Kanwar, J., Kamp, I., Woitke, P., et al. 2025, A&A, submitted [Google Scholar]
  74. Kashyap, P., Majumdar, L., Dutrey, A., et al. 2024, ApJ, 976, 258 [Google Scholar]
  75. Klarmann, L., Ormel, C. W., & Dominik, C. 2018, A&A, 618, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Kress, M. E., Tielens, A. G. G. M., & Frenklach, M. 2010, Adv. Space Res., 46, 44 [NASA ADS] [CrossRef] [Google Scholar]
  77. Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78 [Google Scholar]
  78. Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134 [NASA ADS] [CrossRef] [Google Scholar]
  79. Kurtovic, N. T., Pinilla, P., Long, F., et al. 2021, A&A, 645, A139 [EDP Sciences] [Google Scholar]
  80. Laas, J. C., Garrod, R. T., Herbst, E., & Widicus Weaver, S. L. 2011, ApJ, 728, 71 [Google Scholar]
  81. Laibe, G., & Price, D. J. 2014, MNRAS, 444, 1940 [Google Scholar]
  82. Lamberts, T., Cuppen, H. M., Ioppolo, S., & Linnartz, H. 2013, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 15, 8287 [Google Scholar]
  83. Langner, U. W., & Potgieter, M. S. 2005, ApJ, 630, 1114 [Google Scholar]
  84. Lee, J.-E., Bergin, E. A., & Nomura, H. 2010, ApJ, 710, L21 [Google Scholar]
  85. Lenzuni, P., Gail, H.-P., & Henning, T. 1995, ApJ, 447, 848 [NASA ADS] [CrossRef] [Google Scholar]
  86. Li, J., Bergin, E. A., Blake, G. A., Ciesla, F. J., & Hirschmann, M. M. 2021, Sci. Adv., 7, eabd3632 [NASA ADS] [CrossRef] [Google Scholar]
  87. Lienert, J. L., Bitsch, B., & Henning, T. 2024, A&A, 691, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Lienert, J. L., Bitsch, B., & Henning, T. 2025, A&A, 700, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Ligterink, N. F. W., Kipfer, K. A., & Gavino, S. 2024, A&A, 687, A224 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Ligterink, N. F. W., Pinilla, P., van der Marel, N., et al. 2025, Nat. Astron., 8, 1257 [Google Scholar]
  91. Long, D. E., Cleeves, L. I., Adams, F. C., et al. 2024, ApJ, 972, 88 [Google Scholar]
  92. Long, F., Pascucci, I., Houge, A., et al. 2025, ApJ, 978, L30 [Google Scholar]
  93. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
  94. Mah, J., Bitsch, B., Pascucci, I., & Henning, T. 2023, A&A, 677, L7 [CrossRef] [EDP Sciences] [Google Scholar]
  95. Mah, J., Savvidou, S., & Bitsch, B. 2024, A&A, 686, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. McClure, M. K. 2019, A&A, 632, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167 [Google Scholar]
  98. McClure, M. K., Dominik, C., & Kama, M. 2020, A&A, 642, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Mesquita, A. L., Rodgers-Lee, D., Vidotto, A. A., Atri, D., & Wood, B. E. 2022, MNRAS, 509, 2091 [Google Scholar]
  100. Micolta, M., Calvet, N., Thanathibodee, T., et al. 2023, ApJ, 953, 177 [Google Scholar]
  101. Micolta, M., Calvet, N., Thanathibodee, T., et al. 2024, ApJ, 976, 251 [Google Scholar]
  102. Miley, J. M., Panic´, O., Booth, R. A., et al. 2021, MNRAS, 500, 4658 [Google Scholar]
  103. Millar, T. J., Bennett, A., Rawlings, J. M. C., Brown, P. D., & Charnley, S. B. 1991, A&AS, 87, 585 [NASA ADS] [Google Scholar]
  104. Millar, T. J., Walsh, C., Van de Sande, M., & Markwick, A. J. 2024, A&A, 682, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ACS Earth Space Chem., 6, 597 [NASA ADS] [CrossRef] [Google Scholar]
  106. Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Miret-Roig, N., Galli, P. A. B., Olivares, J., et al. 2022, A&A, 667, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130 [Google Scholar]
  110. Molyarova, T., Vorobyov, E., & Akimkin, V. 2025, PASA, 42, e004 [Google Scholar]
  111. Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471 [NASA ADS] [CrossRef] [Google Scholar]
  112. Najita, J. R., Ádámkovics, M., & Glassgold, A. E. 2011, ApJ, 743, 147 [NASA ADS] [CrossRef] [Google Scholar]
  113. Najita, J. R., Carr, J. S., Pontoppidan, K. M., et al. 2013, ApJ, 766, 134 [Google Scholar]
  114. Nazari, P., Sellek, A. D., & Rosotti, G. P. 2025, A&A, 697, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19 [Google Scholar]
  116. Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 504, 891 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011a, ApJ, 740, 109 [Google Scholar]
  118. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011b, ApJ, 743, L16 [Google Scholar]
  119. Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1 [CrossRef] [Google Scholar]
  120. Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8 [NASA ADS] [CrossRef] [Google Scholar]
  121. Orlando, S., Reale, F., Peres, G., & Mignone, A. 2011, MNRAS, 415, 3380 [Google Scholar]
  122. Pacetti, E., Turrini, D., Schisano, E., et al. 2022, ApJ, 937, 36 [NASA ADS] [CrossRef] [Google Scholar]
  123. Padovani, M., Hennebelle, P., Marcowith, A., & Ferrière, K. 2015, A&A, 582, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., Rosotti, G., & Tabone, B. 2025, A&A, 698, A231 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  126. Pascucci, I., Apai, D., Luhman, K., et al. 2009, ApJ, 696, 143 [Google Scholar]
  127. Pascucci, I., Herczeg, G., Carr, J. S., & Bruderer, S. 2013, ApJ, 779, 178 [NASA ADS] [CrossRef] [Google Scholar]
  128. Pascucci, I., Skinner, B. N., Deng, D., et al. 2023, ApJ, 953, 183 [CrossRef] [Google Scholar]
  129. Pech, G., Loinard, L., Dzib, S. A., et al. 2016, ApJ, 818, 116 [Google Scholar]
  130. Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2021, ApJ, 911, 150 [NASA ADS] [CrossRef] [Google Scholar]
  131. Pontoppidan, K. M., Salyk, C., Banzatti, A., et al. 2024, ApJ, 963, 158 [NASA ADS] [CrossRef] [Google Scholar]
  132. Powell, D., Gao, P., Murray-Clay, R., & Zhang, X. 2022, Nat. Astron., 6, 1147 [CrossRef] [Google Scholar]
  133. Prasad, S. S., & Huntress, Jr., W. T. 1980, ApJS, 43, 1 [NASA ADS] [CrossRef] [Google Scholar]
  134. Prasad, S. S., & Tarafdar, S. P. 1983, ApJ, 267, 603 [Google Scholar]
  135. Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401 [NASA ADS] [CrossRef] [Google Scholar]
  136. Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Raul, E., Alarcón, F., & Bergin, E. A. 2025, ApJ, 982, 155 [Google Scholar]
  138. Rodgers-Lee, D., Taylor, A. M., Ray, T. P., & Downes, T. P. 2017, MNRAS, 472, 26 [Google Scholar]
  139. Rodgers-Lee, D., Taylor, A. M., Downes, T. P., & Ray, T. P. 2020a, MNRAS, 491, 4742 [Google Scholar]
  140. Rodgers-Lee, D., Vidotto, A. A., Taylor, A. M., Rimmer, P. B., & Downes, T. P. 2020b, MNRAS, 499, 2124 [Google Scholar]
  141. Romero-Mirza, C. E., Banzatti, A., Öberg, K. I., et al. 2024a, ApJ, 975, 78 [NASA ADS] [CrossRef] [Google Scholar]
  142. Romero-Mirza, C. E., Öberg, K. I., Banzatti, A., et al. 2024b, ApJ, 964, 36 [CrossRef] [Google Scholar]
  143. Ruaud, M., Gorti, U., & Hollenbach, D. J. 2022, ApJ, 925, 49 [NASA ADS] [CrossRef] [Google Scholar]
  144. Salyk, C., Pontoppidan, K. M., Banzatti, A., et al. 2025, AJ, 169, 184 [Google Scholar]
  145. Schuhmann, M., Altwegg, K., Balsiger, H., et al. 2019, A&A, 630, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2019, ApJ, 877, 131 [Google Scholar]
  147. Schwarz, K. R., Henning, T., Christiaens, V., et al. 2024, ApJ, 962, 8 [NASA ADS] [CrossRef] [Google Scholar]
  148. Seifert, R. A., Cleeves, L. I., Adams, F. C., & Li, Z.-Y. 2021, ApJ, 912, 136 [NASA ADS] [CrossRef] [Google Scholar]
  149. Sellek, A. D., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 498, 2845 [NASA ADS] [CrossRef] [Google Scholar]
  150. Sellek, A. D., Vlasblom, M., & van Dishoeck, E. F. 2025, A&A, 694, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  151. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
  152. Siebenmorgen, R., & Krügel, E. 2010, A&A, 511, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Sinclair, C. A., Rosotti, G. P., Juhasz, A., & Clarke, C. J. 2020, MNRAS, 493, 3535 [Google Scholar]
  154. Sokal, K. R., Johns-Krull, C. M., Mace, G. N., et al. 2020, ApJ, 888, 116 [Google Scholar]
  155. Sturm, J. A., McClure, M. K., Harsono, D., et al. 2022, A&A, 660, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  156. Tabone, B., Bettoni, G., van Dishoeck, E. F., et al. 2023, Nat. Astron., 7, 805 [NASA ADS] [CrossRef] [Google Scholar]
  157. Taquet, V., Furuya, K., Walsh, C., & van Dishoeck, E. F. 2016, MNRAS, 462, S99 [NASA ADS] [CrossRef] [Google Scholar]
  158. Temmink, M., van Dishoeck, E. F., Grant, S. L., et al. 2024, A&A, 686, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  159. Temmink, M., Sellek, A. D., Gasman, D., et al. 2025, A&A, 699, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  160. Trapman, L., Zhang, K., van’t Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022, ApJ, 926, L2 [NASA ADS] [CrossRef] [Google Scholar]
  161. Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41 [NASA ADS] [CrossRef] [Google Scholar]
  162. Trapman, L., Longarini, C., Rosotti, G. P., et al. 2025, ApJ, 984, L18 [Google Scholar]
  163. Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44 [Google Scholar]
  164. Turrini, D., Schisano, E., Fonte, S., et al. 2021, ApJ, 909, 40 [Google Scholar]
  165. Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617 [NASA ADS] [Google Scholar]
  166. Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69 [Google Scholar]
  167. Vaikundaraman, V., Dra˛z˙kowska, J., Binkert, F., Birnstiel, T., & Miotello, A. 2025, A&A, 696, A215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  168. Van Clepper, E., Bergner, J. B., Bosman, A. D., Bergin, E., & Ciesla, F. J. 2022, ApJ, 927, 206 [NASA ADS] [CrossRef] [Google Scholar]
  169. van der Marel, N., & Mulders, G. D. 2021, AJ, 162, 28 [NASA ADS] [CrossRef] [Google Scholar]
  170. van der Marel, N., Bosman, A. D., Krijt, S., Mulders, G. D., & Bergner, J. B. 2021, A&A, 653, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  171. Vlasblom, M., Temmink, M., Grant, S. L., et al. 2025, A&A, 693, A278 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  172. Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  173. Xie, C., Pascucci, I., Long, F., et al. 2023, ApJ, 959, L25 [NASA ADS] [CrossRef] [Google Scholar]
  174. Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017, ApJ, 834, 178 [Google Scholar]
  175. Yu, M., Willacy, K., Dodson-Robinson, S. E., Turner, N. J., & Evans, Neal J. I. 2016, ApJ, 822, 53 [NASA ADS] [CrossRef] [Google Scholar]
  176. Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJ, 891, L17 [Google Scholar]
  177. Zwicky, L., Molyarova, T., Kóspál, Á., & Ábrahám, P. 2025, A&A, 699, A377 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

1

Some other works use ‘depletion’ only if CO is less abundant than predicted given chemical transformation processes as well as freeze-out and photodissociation, which has lead to competing claims (e.g. Ruaud et al. 2022; Pascucci et al. 2023) of whether CO is depleted (see discussion in Section 4.1 of Paneque-Carreño et al. 2025). In the works that find CO ‘depletion’, grain surface processes are not modelled in detail, so CO remains abundant at the 10−4 level. The volatile C abundance is instead artificially reduced in order to approximate the transformation of CO to less volatile species and thus reduce the CO abundance. Ultimately, all works agree on the need to transform CO to a less volatile species in order to reproduce the observed line fluxes.

2

If CRs detour along sheared magnetic field lines, the surface density they traverse to reach the midplane at a distance R can be greater than Σ(R) by around two orders of magnitude (Fujii & Kimura 2022).

3

Ligterink et al. (2025) also suggested that warm, UV-irradiated dust traps could trap ices long enough for the formation of organic macromolecules.

4

This is different to the main sequence, where M dwarfs have a higher level of magnetic activity than solar-type stars due to their faster rotation periods and deeper convective zone, resulting in a higher LX/Lbol. Then, orbiting planets at a given temperature, for example in the habitable zone, would experience a higher X-ray flux (Johnstone et al. 2021).

All Tables

Table 1

Parameters explored in our model grids. Fiducial values are highlighted in bold.

Table A.1

Binding energies, prefactors, and initial abundances (with respect to H) assumed for the species considered in this work and the resulting snowlines around 0.1 and 1.0 M stars (defined as the locations where the gas and ice abundances are equal at t = 0).

Table A.2

Tunnelling barriers assumed for grain-surface reactions.

Table B.1

Maximum C/O achieved within the time range specified for models including different numbers or types of reactions for a fiducial case with RC = 16au, MD/M* = 0.05, and ζH2 = 1.3 × 10−17s−1.

All Figures

thumbnail Fig. 1

Summary of the reaction scheme(s) used in this work. Here, ‘s-’ indicates that the species is in the ice phase. The process in bold is the most critical for raising the inner disc C/O.

In the text
thumbnail Fig. 2

Evolution of different disc models for 10 Myr. a) Models that do not include reactions of C-bearing species for various initial disc-to-star mass ratios and initial disc sizes. b) Models that include reactions of CO. All models have a fixed initial disc-to-star mass ratio of 0.05 but different combinations of ionisation rate and disc size. In each panel, the tracks start at ‘S’ and evolve in the direction of the arrows, ending at ‘F’. The solid sections represent ages of 1–5 Myr appropriate for most discs studied so far with MIRI-MRS. In each block, the left-hand columns show models around a 1.0 M star, representative of TTSs, while the right-hand columns show models around a 0.1 M star, representative of VLMS. The top rows plot the discs in the plane of the inner-disc gas-phase C/H vs gas-phase O/H. Dashed lines indicate where C/O=0.32,1,3,10. The bottom rows plot the gas-phase C/O in the inner disc over time. The horizontal dashed lines show ratios of 0.32 (initial ratio of volatile abundances), 0.43 (ratio of total abundances) and 1 and the lightly shaded regions indicate the typically observed range of C/O at 1–5 Myr for each stellar mass regime.

In the text
thumbnail Fig. 3

C/O ratio as a function of the inner disc H2O vapour mass (calculated as the total integrated H2O vapour mass inside of the midplane snowline) for the same models as Fig. 2b. The line styles are the same as in that figure, and the arrows illustrate the evolution with time.

In the text
thumbnail Fig. 4

Same as for the right-hand panels of Fig. 2b but comparing models with and without dust traps for a RC,0 = 50 au disc around a 0.1 M star with an ionisation rate of ∼10−17 s−1.

In the text
thumbnail Fig. 5

Maximum inner-disc gas-phase C/O achieved for different parameters in smooth disc models with hydrocarbon formation around 1.0 M TTSs (left) or 0.1 M VLMSs (right). Darker green colours indicate more C-rich discs, while darker purple colours indicate more O-rich discs. We note that the colour bar is neither linear nor symmetric. The central circles show the maxima within the first 5 Myr, while the outer rings show the maxima across the full 10 Myr model run.

In the text
thumbnail Fig. 6

Summary of differences between VLMS discs and TTS discs in this work. For the purposes of this comparison, we assumed a TTS model with the typically measured ionisation rate of ~10−18 s−1 and an average RC,0 = 16 au. We assumed that the VLMS has a more compact disc of RC,0 = 5 au subject to an unexcluded galactic CR ionisation rate of ~10−17 s−1.

In the text
thumbnail Fig. A.1

Top: Midplane ionisation rate profile at t=0, 1, and 5 Myr (increasingly dark blues, evolution indicated by arrow) for fiducial models around TTS (left) and VLMSs (right). The TTS models are compared to observationally estimated values for individual sources (red and yellow colours) and the range and median for the exoALMA program (Trapman et al. 2025). The black horizontal line on each plot represents the floor due to SLRs and the vertical dotted line shows the location of the CO snowline. Bottom: Surface densities of the gas at the times indicated in the top panels. The horizontal line shows the assumed 96 g cm−2 penetration depth of CRs.

In the text
thumbnail Fig. B.1

Same as in Fig. 2a but for disc models with different numbers or types of reactions included.

In the text
thumbnail Fig. B.2

Same as in the lower panels of Fig. 2b but for disc models that include formation of C2Hx hydrocarbons.

In the text
thumbnail Fig. B.3

Evolution of molecular gas-phase column density ratios averaged across the disc inside their respective snowlines. Top: C2H6 and C2H2. Bottom: CO2 and C2H2.

In the text
thumbnail Fig. B.4

Comparison of the timescales at t = 0 yr around a 1.0 M star (left) and 0.1 M star (right) for the erosion of carbonaceous grains by UV photolysis (blue, solid) with the timescales of gas accretion (orange, dotted; see Equation 12), dust grain growth (green, dashed) and radial drift (red, dashed). The carbon depletion timescale at 1 au in the ‘fiducialfc025’ model from Vaikundaraman et al. (2025) is shown as the blue V.

In the text
thumbnail Fig. D.1

C/O as a function of R at 2.5 Myr for models with no trap (black) or a warm trap (red) in a 50 au disc around a 0.1 M star with an ionisation rate of 1.3 × 10−17 s−1. Vertical lines mark the location of the minimum gas density (’Gap location’) and the pressure maximum where dust accumulates (’Dust ring’).

In the text
thumbnail Fig. D.2

Same as in Fig. 4 but for 1 M.

In the text
thumbnail Fig. D.3

Same as in Fig. 4 and D.2 but comparing models with and without warm dust traps for a range of initial disc radii.

In the text
thumbnail Fig. E.1

Relationship between the inner disc C/O (measured at the 0.1 au inner boundary of the models) and the outer disc (measured at 90% of the CH4 snowline) C abundance (left) and C/O (right) for discs around VLMSs.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.