1. Introduction
The protoplanetary disc matter can be roughly divided into three component: gaseous chemical species, solid dust particles, and icy mantles covering the surface of dust grains. Gas and solid particles become dynamically decoupled, as evolving dust grows and acquires relative velocities leading to the redistribution of elements in the disc and between the phases, and creating the premises for different chemical environments. When planets start to form, their properties, including chemical composition of the atmosphere, are inevitably affected by the location and the mechanism of their formation. This suggests that the origin of (exo)planets might be investigated using their observed chemical composition, and makes understanding the disc chemical evolution vital for creating a consistent planet formation theory.
One of the key parameters that govern the chemical setup of a planetary atmosphere is the relation between the abundances of carbon and oxygen, often referred to as carbon-to-oxygen ratio (hereafter C/O ratio). The variations of C/O ratio in the ice and gas phases at the snowlines of main disc volatiles (CO, CO
$_2$
, and H
$_2$
O) and the prospects of connecting them to planet formation were discussed in Öberg, Murray-Clay, & Bergin (Reference Öberg, Murray-Clay and Bergin2011) within a qualitative freeze-out model. Since then, C/O ratio received a lot of attention in this context. It was thoroughly investigated in modelling (see, e.g. Booth et al. Reference Booth, Clarke, Madhusudhan and Ilee2017; Eistrup et al. Reference Eistrup, Walsh and van Dishoeck2018; Cridland, Eistrup, & van Dishoeck Reference Cridland, Eistrup and van Dishoeck2019a; Cridland et al. 2019b; Cridland, Bosman, & van Dishoeck Reference Cridland, Bosman and van Dishoeck2020a; Cridland et al. 2020b; Krijt et al. Reference Krijt, Bosman, Zhang, Schwarz, Ciesla and Bergin2020; Turrini et al. Reference Turrini2021; Schneider & Bitsch Reference Schneider and Bitsch2021). The connection of disc chemical composition with C/O in exoplanetary atmospheres was modelled using core accretion model (Thiabaud et al. Reference Thiabaud, Marboeuf, Alibert, Leya and Mezger2015) and ‘chain’ planet population synthesis model (Mordasini et al. Reference Mordasini, van Boekel, Mollière, Henning and Benneke2016). Mollière et al. (Reference Mollière2022) considered a simple formation retrieval pipeline and found that this task requires careful consideration of the model assumptions.
The measurements of molecular abundances in the atmospheres of giant exoplanets obtained by a variety of modern facilities, such as HST, Spitzer, VLTI, JWST, Gemini, indicate a diversity of C/O ratios: from low C/O ratio values (
$\approx0.4$
, below the solar value of
$0.54$
; Benneke et al. Reference Benneke2019;
GRAVITY Collaboration et al. 2020; Worthen et al. Reference Worthen2024; Xue et al. Reference Xue, Bean, Zhang, Welbanks, Lunine and August2024) to stellar (
$\approx0.5$
, close to solar; Mollière et al. Reference Mollière2020; Zhang et al. Reference Zhang2021; Smith et al. Reference Smith2024) and close to or above unity (Swain et al. Reference Swain2009; Madhusudhan et al. Reference Madhusudhan2011). A variety of solar and super-solar C/O ratios is observed in four planets within the HR8799 system (Nasedkin et al. Reference Nasedkin2024). Chemical composition of the atmospheres of many hot Jupiters indicates high C/O>1 of the forming material (Moses et al. Reference Moses, Madhusudhan, Visscher and Freedman2013). There is an observational evidence of young planets in PDS 70 disc accreting the material with C/O>1 (Facchini et al. Reference Facchini, Teague, Bae, Benisty, Keppler and Isella2021). The number of exoplanets with constrained atmospheric C/O ratios grows with large studies of multiple planets such as Changeat et al. (Reference Changeat2022), which allows us to make some statistical conclusions. The population study of C/O ratios in exoplanetary atmospheres reveals that there are two populations with different elemental ratio, which are likely formed in different mechanisms (Hoch et al. Reference Hoch2023). Khorshid, Min, & Désert (Reference Khorshid, Min and Désert2023) were able to restrict the formation scenario for WASP-77b based on the measured C/O ratio of the planet (Line et al. Reference Line2021) and the modelling of planet formation and migration.
Elemental abundances in protoplanetary discs can be constrained from observations (Fedele & Favre Reference Fedele and Favre2020), and C/O ratio in the gas can be estimated. Spatially resolved observations can help distinguish between different C/O ratios spectroscopically (Matter, Pignatale, & Lopez Reference Matter, Pignatale and Lopez2020). Cleeves et al. (Reference Cleeves, Öberg, Wilner, Huang, Loomis, Andrews and Guzman2018) report
$\textrm{C/O}\approx0.8$
in the molecular layer of IM Lup disc. ALMA observations of hydrocarbons and sulphur-bearing species indicate
$\textrm{C/O} \gt 1$
in the upper disc layers and in the outer disc in TW Hya and DM Tau (Dutrey et al. Reference Dutrey2011; Bergin et al. Reference Bergin, Du, Ilsedore Cleeves, Blake, Schwarz, Visser and Zhang2016; Semenov et al. Reference Semenov2018) and for a population of discs in Lupus (Miotello et al. Reference Miotello2019). For the nearby discs the solar elemental composition with
$\textrm{C/O}\approx0.54$
is usually expected, thus the observed higher values confirm redistribution of carbon and oxygen in discs. In addition to high C/O in disc atmospheres, there is evidence of both carbon and oxygen depletion from gas (Kama et al. Reference Kama2016; Miotello et al. Reference Miotello2019). However, some of heavy oxygen carriers might not be observable, leading to overestimated C/O in disc observations (Walsh et al. Reference Walsh, Nomura and van Dishoeck2015).
The volatile composition is also used to constrain the origin of bodies in the Solar System. Fraction of CO and CO
$_2$
ices relative to water in cometary comae indicate their formation between the CO and CO
$_2$
snowlines or exterior to the CO snowline (A’Hearn et al. 2012; Seligman et al. Reference Seligman2022). Abundances of CO and N
$_2$
ices were used to analyse the original location of Pluto and Triton (Mousis et al. 2024a). Observed elemental abundances were used to constrain the Jupiter formation scenario (Lodders Reference Lodders2004), relying also on abundances of nitrogen (Öberg & Wordsworth Reference Öberg and Wordsworth2019; Bosman, Cridland, &Miguel Reference Bosman, Cridland and Miguel2019) and chemically inactive species like Ar. However, the model assumptions can lead to different interpretation of the observations: while Öberg &Wordsworth (Reference Öberg and Wordsworth2019) and Bosman et al. (Reference Bosman, Cridland and Miguel2019) suggest that Jupiter formed outside N
$_2$
snowline (at
$ \gt 30$
au), Ohno & Ueda (Reference Ohno and Ueda2021) consider the concept disc shadow, which allows Jupiter to form near its current location.
Chemical processes other than freeze-out and sublimation at the snowlines can alter the composition of ice and gas as well. Due to gas-phase and surface reactions, snowlines can become important for the redistribution of elements. More detailed chemical modelling shows that the C/O ratio in the gas and in the ice depends also on the initial chemical setup and ionisation by cosmic rays and radioactive nuclei (Eistrup, Walsh, & van Dishoeck Reference Eistrup, Walsh and van Dishoeck2016; Eistrup et al. Reference Eistrup, Walsh and van Dishoeck2018). It directly affects the interpretation of observations. Another essential chemical process is CO depletion from the gas, resulting in its transformation to CO
$_2$
ice (Bosman, Tielens, & van Dishoeck Reference Bosman, Tielens and van Dishoeck2018). For example, stellar C/O ratio in the atmosphere of HR 8799e indicates that the planet accreted its material beyond CO snowline (
$\approx$
45 au), but chemical modelling suggests that due to CO depletion, the C/O in the ice already approaches the stellar ratio beyond CO
$_2$
snowline (
$\approx$
20 au), which is closer to the star (Mollière et al. Reference Mollière2020).
Another key process affecting the elemental ratio is dust drift, which leads to spatial segregation between the chemical constituents of the gas and the grains covered with ice. The distribution of CO in the gas and ice phases was studied within dynamical models of dust evolution (Stammler et al. Reference Stammler, Birnstiel, Panić, Dullemond and Dominik2017; Krijt et al. Reference Krijt, Schwarz, Bergin and Ciesla2018). Even without chemical processes, dust evolution and dynamics can substantially alter C/O ratio in the atmospheres of forming planets (Booth et al. Reference Booth, Clarke, Madhusudhan and Ilee2017). Some models combine chemical reactions treatment with dust evolution and transport, usually within 1D viscous models. Dust transport can have a strong effect on the abundances of volatiles in the inner disc regions (Bosman et al. Reference Bosman, Tielens and van Dishoeck2018). However, for discs with low turbulence and high cosmic ray ionisation rate, C/O ratio is rather defined by chemical evolution (Booth & Ilee Reference Booth and Ilee2019).
In our previous work (Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021) we showed that the volatile species tend to concentrate around their snowlines both in the gas and more notably on the dust surface. This accumulation was found to be caused by effective transport of volatiles through the snowlines by azimuthal variations in the gas and dust radial and angular velocity, an effect that cannot be captured in 1D viscous disc models. Such accumulation should immediately affect the local C/O ratio, which suggests the connection between the snowlines of various volatiles and the formation of planets with altered C/O in their atmospheres. In this work, we follow the distribution of the main volatile species in the disc to investigate the distribution of C/O ratio in gas and ice in a 2D thin-disc hydrodynamic model. We study the effect of dust growth and dynamics on the elemental ratios and consider the role of the initial mass of the collapsing core on the distribution of volatiles.
The paper is organised as follows. The main features of the used FEOSAD model are described in Section 2, with the details of the treatment of the volatiles given in Section 2.3. In Section 3 we describe the results of the simulations, focusing on distribution of the volatiles in Section 3.2, the C/O ratios in Section 3.3, and their evolution in Section 3.4. In Section 4, we discuss the implications of our results in the context of planet formation via different mechanisms. The main conclusions are listed in Section 5.
2. Model
We use the global model of protoplanetary disc formation FEOSAD (Vorobyov et al. Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018), which includes disc self-gravity, dust evolution and interaction with gas (including backreaction of dust on gas), turbulent viscosity, adiabatic and radiative cooling and heating. It describes the formation of a protostar and a protoplanetary disc from a collapsing cloud in a 2D thin-disc approach. The model also includes freeze-out of main volatile species as in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021), with the feedback from ice mantles on dust evolution via fragmentation velocity. Here we summarise the key characteristics of the model, more details can be found in the previous works (Vorobyov et al. Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018; Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021; Kadam, Vorobyov, & Basu Reference Kadam, Vorobyov and Basu2022). The main difference from our previous study in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021) is that here we consider the formation of dead zones via variable
$\alpha$
-parameter of Shakura and Sunyaev and also include turbulent diffusion.
2.1 Gas evolution
For the gas component, the hydrodynamic equations for mass, momentum, and internal energy conservation are the following
\begin{align}\frac{\partial \Sigma_\textrm{g}}{\partial t} + \nabla \cdot\left( \Sigma_\textrm{g} \boldsymbol{v} \right) =0, \\[-22pt] \nonumber \end{align}
\begin{align}\frac{\partial}{\partial t} \left(\Sigma_\textrm{g} \boldsymbol{ v } \right) & + [\nabla \cdot \left( \Sigma_\textrm{g} \boldsymbol{v } \otimes \boldsymbol{v } \right)] = - \nabla \mathcal{P} + \Sigma_\textrm{g} \, \boldsymbol{g } + \nonumber\\ & + \nabla \cdot \boldsymbol{\Pi} - \Sigma_\textrm{d,gr} \boldsymbol{f }, \\[-22pt] \nonumber \end{align}
\begin{align}\frac{\partial e}{\partial t} +\nabla \cdot \left( e \boldsymbol{v} \right) = -\mathcal{P}(\nabla \cdot \boldsymbol{v }) -\Lambda +\Gamma +\nabla \boldsymbol{v}:\Pi,\end{align}
where subscripts p and
$p^\prime$
denote the planar components
$(r,\phi)$
in polar coordinates,
$\Sigma_\textrm{g}$
is the gas mass surface density, e is the internal energy per surface area,
$\mathcal{P}$
is the vertically integrated gas pressure calculated via the ideal equation of state as
$\mathcal{P}=(\gamma-1) e$
with
$\gamma=7/5$
, f is the friction force between gas and dust,
$\boldsymbol{ v}=v_r \hat{\boldsymbol{ r}} v_\phi \hat{\boldsymbol{\phi}}$
is the gas velocity in the disc plane, and
$\nabla=\hat{\boldsymbol{r}} \partial / \partial r + \hat{\boldsymbol{\phi}} r^{-1} \partial / \partial \phi $
is the gradient along the planar coordinates of the disc. The gravitational acceleration in the disc plane,
$\boldsymbol{g}=g_r \hat{\boldsymbol{r}}+g_\phi \hat{\boldsymbol{\phi}}$
, includes the gravity of the central protostar when formed and takes into account disc self-gravity of both gas and dust found by solving the Poisson integral (Binney & Tremaine Reference Binney and Tremaine1987).
The consideration of time-dependent energy balance (Equation 3) allows us to accurately calculate the midplane temperature
$T_\textrm{mp}$
and is particularly important to describe the phase state of the volatiles and the level of turbulent viscosity. The terms
$\Lambda$
and
$\Gamma$
describe the rates of dust cooling and heating, respectively, by stellar and background irradiation. They are calculated based on the analytical solution of the radiation transfer equations in the vertical direction (Dong et al. Reference Dong, Vorobyov, Pavlyuchenkov, Chiang and Liu2016; Vorobyov et al. Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018)
\begin{align}\Lambda = \frac{8 \tau_\textrm{P}\sigma T_\textrm{mp}^4}{1+2\tau_\textrm{P}+\frac{3}{2}\tau_\textrm{P}\tau_\textrm{R}}, \textrm{ } \Gamma = \frac{8 \tau_\textrm{P}\sigma T_\textrm{irr}^4}{1+2\tau_\textrm{P}+\frac{3}{2}\tau_\textrm{P}\tau_\textrm{R}}.\end{align}
Here,
$\sigma$
is the Stefan-Boltzmann constant,
$\tau_\textrm{P}$
and
$\tau_\textrm{R}$
are the Planck and Rosseland mean optical depths to the disc midplane, calculated as
$\tau=\kappa \Sigma_\textrm{dust}$
from Planck and Rosseland mean opacities
$\kappa_\textrm{P}$
and
$\kappa_\textrm{R}$
(Semenov et al. Reference Semenov, Henning, Helling, Ilgner and Sedlmayr2003) and total dust surface density
$\Sigma_\textrm{dust}$
. Gas and dust temperatures are assumed to be equal, and the midplane temperature is linked with gas pressure as
$T_\textrm{mp} = \mathcal{P}\mu/\mathcal{R}\Sigma_\textrm{g}$
, where
$\mu=2.3$
is the mean molecular weight of the gas, and
$\mathcal{R}$
is the universal gas constant. The irradiation temperature at the disc surface
$T_\textrm{irr}$
is determined by both stellar and background irradiation. Stellar irradiation includes the luminosity from the photosphere of the protostar and accretion luminosity. The background radiation is assumed as a black body with the temperature of 15 K. For more details on the irradiation we refer to Vorobyov et al. (Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018).
Turbulent viscosity is described using the common
$\alpha$
-parameter approach of Shakura & Sunyaev (Reference Shakura and Sunyaev1973). It is taken into account via the viscous stress tensor
${\boldsymbol{\Pi}}$
(see Vorobyov & Basu Reference Vorobyov and Basu2010, for explicit expressions for the components of the terms with
${\boldsymbol{\Pi}}$
). The magnitude of kinematic viscosity is
$\nu=\alpha c_\textrm{s} H$
, where
$c_\textrm{s}$
is the sound speed and H is the vertical scale height of the gas disc calculated using an assumption of local hydrostatic equilibrium of a self-gravitating disc (see Vorobyov & Basu Reference Vorobyov and Basu2010, Appendix A). Here, we use the adaptive
$\alpha$
approach implying accretion through a layered disc (Gammie Reference Gammie1996; Armitage, Livio, & Pringle Reference Armitage, Livio and Pringle2001; Kadam et al. Reference Kadam, Vorobyov and Basu2022). Turbulence is assumed to be generated by magneto-rotational instability (MRI) which only develops in layers of the disc where the ionisation level is high enough. The MRI-active layer is characterised by its surface density
$\Sigma_\textrm{MRI}$
and relatively high value of turbulent viscosity
$\alpha_\textrm{MRI}=10^{-3}$
. As thermal and photo-ionisation are not efficient enough for the relatively cold and dense matter in the disc at >0.5 au, the main process determining the thickness of the MRI-active layer is ionisation by cosmic rays. It is assumed to be constant
$\Sigma_\textrm{MRI}=100$
g cm-2, which is the typical depth of Galactic cosmic rays penetration in the ISM (Umebayashi & Nakano Reference Umebayashi and Nakano1981) and in protoplanetary discs (Padovani et al. Reference Padovani, Galli and Glassgold2009). The dead zone is characterised by surface density from the midplane
$\Sigma_\textrm{ dz}=\Sigma_\textrm{g}/2-\Sigma_\textrm{MRI}$
with residual turbulence
$\alpha_\textrm{dz}$
. The turbulence in this layer is only hydrodynamic turbulence driven by the Maxwell stress in the active layer, and small value
$\alpha_\textrm{dz}=10^{-5}$
is adopted (Okuzumi & Hirose Reference Okuzumi and Hirose2011). However, if local temperature exceeds the critical value of 1 300 K, thermal ionisation becomes possible, the MRI develops and the dead zone is no longer dead, in which case
$\alpha_\textrm{ dz}=10^{-1}$
(Zhu et al. Reference Zhu, Hartmann, Gammie, Book, Simon and Engelhard2010; Kadam et al. Reference Kadam, Vorobyov, Regály, Kóspál and Ábrahám2019). This value is higher than
$\alpha_\textrm{MRI}$
in the outer disc due to the different ionisation processes and local conditions. In the outer disc, the MRI can be suppressed by non-ideal magneto-hydrodynamic (MHD) effects such as ambipolar diffusion and Ohmic resistivity (Bai & Stone Reference Bai and Stone2013; Gressel et al. Reference Gressel, Turner, Nelson and McNally2015). In the dead zone in the inner disc,
$\alpha$
can reach higher values when the MRI is triggered by thermal ionisation, as shown by 3D MHD simulations (Zhu, Jiang, & Stone Reference Zhu, Jiang and Stone2020). The adopted parameterisation makes use of an effective parameter
$\alpha_\textrm{eff}$
, which at any given location is calculated as
\begin{align} & \alpha_\textrm{eff}=\frac{\Sigma_\textrm{MRI}\alpha_\textrm{MRI}+\Sigma_\textrm{dz}\alpha_\textrm{dz}}{\Sigma_\textrm{ MRI}+\Sigma_\textrm{dz}}, \nonumber \\ & \alpha_\textrm{dz}=\begin{cases}10^{-5}, & \textrm{ if } T_\textrm{mp} \lt 1\,300 \textrm{K}; \\10^{-1}, & \textrm{ if } T_\textrm{mp} \geq 1\,300 \textrm{K}.\end{cases}\end{align}
2.2 Dust evolution
Dust is described as consisting of two components: small grains that are dynamically coupled to the gas, with the mass surface density
$\Sigma_\textrm{d,sm}$
, and grown grains with the mass surface density
$\Sigma_\textrm{d,gr}$
that can move relative to the gas and change in size. The total dust surface density necessary for the calculation of the optical depths for Equation (4) is
$\Sigma_\textrm{dust}=(\Sigma_\textrm{d,gr}+\Sigma_\textrm{d,sm})/2$
. The factor
$1/2$
appears as the optical depths is calculated to the midplane. Each dust population has a power-law size distribution
$f(a)= dN/da = C a^{-p}$
with a normalisation constant C and a fixed exponent
$p=3.5$
. Small dust has sizes between
$a_\textrm{min}=5\times 10^{-7}$
cm and
$a_*=10^{-4}$
cm, grown dust has sizes between
$a_*$
and
$a_\textrm{ max}$
, which can vary due to dust coagulation and fragmentation. Dynamics of these dust components follows the continuity and momentum equations
\begin{align}\frac{{\partial \Sigma_\textrm{d,gr} }}{{\partial t}} + \nabla \cdot\left( \Sigma_\textrm{d,gr} \boldsymbol{u} \right) = S(a_\textrm{max}) + \nabla \cdot \left( D \Sigma_\textrm{g} \nabla\left(\frac{\Sigma_\textrm{d,gr}}{\Sigma_\textrm{g}} \right)\right),\end{align}
\begin{eqnarray}\frac{\partial}{\partial t} \left( \Sigma_\textrm{d,gr} \boldsymbol{u}\right) + \left[\nabla \cdot \left( \Sigma_\textrm{d,gr} \boldsymbol{u} \otimes \boldsymbol{u} \right)\right] &=& \Sigma_\textrm{d,gr} \, \boldsymbol{g}+ \nonumber \\ + \Sigma_\textrm{d,gr} \boldsymbol{f} + S(a_\textrm{max}) \boldsymbol{v},\end{eqnarray}
where
$\boldsymbol{u}$
is the grown dust velocity. The term
$S(a_\textrm{max})$
is responsible for the exchange of matter between the dust populations, as dust is converted from the small to grown component due to coagulation and back due to fragmentation. The details of the dust evolution model are presented in Vorobyov et al. (Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018). The last term in Equations (6) and (7) is responsible for dust turbulent diffusion, similar to Vorobyov et al. (2020a). The coefficient of turbulent diffusion D is related to the kinematic viscosity
$\nu$
as
$D=\nu/(1+\textrm{St}^2)$
(Birnstiel Reference Birnstiel2023). Diffusion affects dust grains along with their ice mantles, as well as the gas-phase species (see Section 2.3).
The innermost regions of the disc are challenging to simulate explicitly due to the Courant criterion: in the highly dynamic inner regions (fraction of au) the timescales are so short that the code demands very small time step in order to preserve stability. Therefore, the inner regions are represented by a sink cell, with a carefully chosen inflow-outflow boundary condition at the sink cell and a parametric description of the accretion onto the star (see Vorobyov et al. Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018, for details). In the simulations presented below, the radius of the sink cell is 0.62 au.
We consider two disc models with different initial mass of the collapsing cloud, 0.66 and 1 M
$_{\odot}$
. We note that about 10% of the gas mass that crosses the sink cell is assumed to be evacuated by jets and outflows, and the other 90% lands on to the star. A small amount of mass remains in the envelope by the end of simulations. In both models, the initial gas temperature
$T_\textrm{init}=15$
K and the ratio of rotational to gravitational energy
$\beta=0.28$
%. The simulations start from the collapse of a molecular cloud, with only small dust grains. The simulation continues until the age of the system becomes equal to 0.5 Myr. Masses of the central protostar and the disc by the end of simulation are
$M_{\star} = 0.4$
M
$_{\odot}$
and
$M_\textrm{disc} = 0.22$
M
$_{\odot}$
in model M1 and
$M_{\star} = 0.58$
M
$_{\odot}$
and
$M_\textrm{disc} = 0.35$
M
$_{\odot}$
in model M2. The disc masses are around
$0.5$
stellar masses, which makes them essentially self-gravitating.
A number of recent studies develop the idea that the mass infall from the ambient ISM continues during the lifetime of the disc, including the Class II stage (Padoan et al. Reference Padoan, Pan, Pelkonen, Haugboelle and Nordlund2024; Pelkonen et al. Reference Pelkonen, Padoan, Juvela, Haugbølle and Nordlund2024; Winter, Benisty, & Andrews Reference Winter, Benisty and Andrews2024). Such models describe a Bondi–Hoyle accretion regime and are in good agreement with the observed properties of the disc population, such as accretion rates, masses and sizes. This input of matter can play an important role in disc evolution and planet formation process (Vorobyov, Lin, & Guedel Reference Vorobyov, Lin and Guedel2015; Vorobyov et al. Reference Vorobyov, Regaly, Guedel and Lin2016). In the FEOSAD model, the mass infall to the disc can be accounted for Vorobyov et al. (2015), (Reference Vorobyov, Regaly, Guedel and Lin2016), but in the present simulation this effect is not included. The mass infall from the envelope continues until the cloud is depleted of matter. Because of the thin-disc geometry, the gravitational contraction of the cloud proceeds in the plane of the disc. The matter is landing on the disc outer edge and is transported towards the star by the combined action of gravitational and viscous torques. The infall on the disc inner regions is therefore neglected. This is a reasonable approximation, considering that most of the matter and angular momentum in a three-dimensional cloud is located at relatively large polar angles and a flared outer edge of the disc intercepts most of them (see, e.g. Visser et al. Reference Visser, van Dishoeck, Doty and Dullemond2009, Figure 1). In our modelling, we do not consider the continuous Bondi-Hoyle accretion and concentrate on the internal disc processes.
Surface density of gas and grown dust, Toomre Q-parameter, maximum dust radius, temperature and viscous
$\alpha$
-parameter in model M1 at selected time moments: 160 kyr,
$80\times80$
au; 350 kyr,
$35\times35$
au; 490 kyr,
$9\times9$
au. The contours indicate the position of the water snowline. Note that at the panels with multiple water snowlines, water is frozen outside the outer line and inside an inner dust ring at 1–2 au.
2.3 Evolution of volatiles
We follow the evolution of four main volatile species: H
$_2$
O, CO
$_2$
, CO, and CH
$_4$
. These are the most abundant carbon- and oxygen-bearing ices observed in protostellar cores (Öberg et al. 2011a). In the model, each of these species can be present in three states: in the gas, in the ice on the surface of small dust, and in the ice on the surface of grown dust. Each species s is described by its surface density in the gas
$\Sigma_{s}^\textrm{gas}$
, on small dust
$\Sigma_{s}^\textrm{sm}$
, and on grown dust
$\Sigma_{s}^\textrm{gr}$
. Their distributions in the disc can change through three main processes: advection together with the corresponding component (gas or small/grown dust); exchange of mantles between dust populations due to grain collisions; phase transitions, including adsorption from gas to dust, and thermal and photo-desorption. Initially, all ices are on small grains and no volatiles present in the gas. The treatment of volatiles is adopted from Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021), who describe the models in more details. Here we recap main features of the chemical model.
Our chemical model only describes phase transitions, i.e. adsorption and desorption, which includes thermal desorption and photo-desorption by interstellar UV radiation. These reactions were shown to be the most important for gas-phase abundances of most species (Ilee et al. Reference Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2011). Due to high computational costs, no other chemical processes, either gas-phase or surface reactions are included, although they also may have significant effect on the composition of both ice and gas (Semenov & Wiebe Reference Semenov and Wiebe2011). The chemical evolution of the surface densities of volatile species is calculated from the system of equations
\begin{align} \frac{\partial \Sigma_{s}^\textrm{gas}}{\partial t} & + \nabla \cdot \left(\Sigma_{s}^{\textrm{gas}} \boldsymbol{v} \right) - \nabla \cdot \left( D \Sigma_\textrm{g} \nabla\left(\frac{\Sigma_{s}^\textrm{gas} }{\Sigma_\textrm{g}} \right)\right) = \nonumber \\ & -\lambda_s\Sigma_{s}^\textrm{ gas}+\eta_s^\textrm{sm}+\eta_s^\textrm{gr}, \\[-20pt] \nonumber \end{align}
\begin{align} \frac{\partial \Sigma_{s}^\textrm{sm}}{\partial t} & + \nabla \cdot \left(\Sigma_{s}^\textrm{sm} \boldsymbol{v} \right) - \nabla \cdot \left( D \Sigma_\textrm{g} \nabla\left(\frac{\Sigma_{s}^\textrm{sm}}{\Sigma_\textrm{g}} \right)\right) = \nonumber \\ & \lambda_s^\textrm{ sm}\Sigma_{s}^\textrm{gas}-\eta_s^\textrm{sm}, \\[-20pt] \nonumber \end{align}
\begin{eqnarray} \frac{{\partial \Sigma_{s}^\textrm{gr} }}{{\partial t}} + \nabla \cdot \left( \Sigma_{s}^\textrm{gr} \boldsymbol{u} \right) -\nabla \cdot \left( D \Sigma_\textrm{g} \nabla\left(\frac{\Sigma_{s}^\textrm{gr}}{\Sigma_\textrm{g}} \right)\right) &=& \nonumber \\ \lambda_s^\textrm{ gr}\Sigma_{s}^\textrm{gas}-\eta_s^\textrm{gr}, \end{eqnarray}
where the mass rate coefficients per disc unit area of adsorption
$\lambda_{s}$
and desorption
$\eta_{s}$
for the species s are calculated for local conditions at every hydrodynamic step, separately for small and grown dust populations. Equations (9)–(11) are solved in two steps: first, the right-hand side is considered without the advection term. The case of pure adsorption/desorption represented by the right-hand side of the equation can be solved analytically (see Appendix A in Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). This is done at every hydrodynamic step, before the dust growth step, when ices on small and grown grains are redistributed proportionally to mass exchange between the dust populations. This is followed by a transport step, when the surface densities of the volatiles are changed according to the fluxes of their respective gas or dust components between the cells. Restricting chemical processes to only adsorption and desorption allows the chemical step to be calculated fast, which is very important for computationally demanding hydrodynamic simulations.
For each dust population, the desorption rate is a sum of thermal desorption and photo-desorption
$\eta=\eta_\textrm{ td}+\eta_\textrm{pd}$
(the indices ‘sm’ and ‘gr’ are omitted for convenience). We operate under the assumption of zeroth-order desorption, i.e. the desorption rate does not depend on the present amount of ice (
$\Sigma_s^\textrm{sm}$
or
$\Sigma_s^\textrm{gr}$
). It implies that only the upper layers of the ice mantle are able to sublimate, which is a more appropriate approach for thick mantles. This assumption better describes desorption of CO and H
$_2$
O in temperature programmed desorption (TPD) experiments (Fraser et al. Reference Fraser, Collings, McCoustra and Williams2001; Öberg et al. Reference Öberg, van Broekhuizen, Fraser, Bisschop, van Dishoeck and Schlemmer2005; Bisschop et al. Reference Bisschop, Fraser, Öberg, van Dishoeck and Schlemmer2006) than first-order desorption, which is more suitable for thin mantles of several monolayers. The rate of thermal desorption is calculated as
\begin{align}\eta_\textrm{td} = \widetilde{\sigma}_\textrm{tot} N_\textrm{ss} \mu_\textrm{s} m_\textrm{p} \sqrt{\frac{2 N_\textrm{ss} E_\textrm{b}k_\textrm{B}}{\pi^2 \mu_\textrm{s} m_\textrm{p}}} \exp{\left(-\frac{E_\textrm{b}}{T_\textrm{mp}}\right)},\end{align}
where
$N_\textrm{ss}=10^{15}$
cm-2 is the surface density of binding sites (Cuppen et al. Reference Cuppen, Walsh, Lamberts, Semenov, Garrod, Penteado and Ioppolo2017),
$E_\textrm{b}$
(K) is the binding energy of the species to the surface,
$\mu_\textrm{s}$
is the species molecular mass,
$m_\textrm{p}$
is atomic mass unit,
$k_\textrm{B}$
is the Boltzmann constant. We follow Hasegawa & Herbst (Reference Hasegawa and Herbst1993) and use binding energy to calculate the pre-exponential factor in Equation (12), the same way it was done in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). However, this approach was recently demonstrated by Minissale et al. (Reference Minissale2022) to underestimate the pre-exponential factor by a few orders of magnitude, as it does not account for the rotational part of the partition functions of desorbing molecules. Total surface area of dust grains (small or grown) per disc unit area
$\widetilde{\sigma}_\textrm{tot}$
(cm
$^2$
cm-2) is calculated for the adopted power-law size distribution with
$p=3.5$
as
\begin{eqnarray}\widetilde{\sigma}_\textrm{tot}^\textrm{sm}&=&\frac{3 \Sigma_\textrm{d, sm}}{\rho_\textrm{s}\sqrt{a_\textrm{min} a_*}},\end{eqnarray}
\begin{eqnarray}\widetilde{\sigma}_\textrm{tot}^\textrm{gr}&=&\frac{3 \Sigma_\textrm{d, gr}}{\rho_\textrm{s}\sqrt{a_* a_\textrm{max}}}.\end{eqnarray}
Here,
$\rho_\textrm{s}=3$
g cm-3 is the material density of silicate cores of the dust grains.
The model includes photodesorption of volatiles by interstellar irradiation, which is mostly relevant in the outer disc regions with lower optical depth. We do not consider the UV radiation field produced by the star and the accretion region as a source of photodesorption, assuming that they do not reach disc midplane due to high optical depth. The photo-desorption rate is calculated as
\begin{align} \eta_\textrm{pd} =\widetilde{\sigma}_\textrm{tot} \mu_\textrm{sp}m_\textrm{p}Y F_\textrm{UV}.\end{align}
where
$F_\textrm{UV}=F_0^\textrm{UV}G_\textrm{UV}$
(photons cm-2 s-1) is the UV photon flux expressed in the units of standard UV field and
$Y = 3.5 \times 10^{-3} + 0.13 \exp\left({-336\,\textrm{K}/T_\textrm{mp}}\right)$
(mol photon-1) is the photodesorption yield adopted from Westley et al. (Reference Westley, Baragiola, Johnson and Baratta1995) for water ice. The intensity of the interstellar UV radiation field with
$G_0=1$
is
$F^\textrm{ UV}_0=4.63\times10^7$
photon cm-2 s-1 (Draine Reference Draine1978). We assume that the disc situated in a star-forming region is illuminated by a slightly elevated unattenuated UV field with
$G_\textrm{env}=5.5G_0$
. For the disc midplane, which is illuminated from above and below, this field scales with the UV optical depth
$\tau_\textrm{UV}$
towards the disc midplane as
\begin{align} G_\textrm{UV}=0.5 G_\textrm{env} e^{-\tau_\textrm{UV}}. \end{align}
The optical depth can be calculated as
$\tau_\textrm{UV}= 0.5 (\varkappa_\textrm{sm} \Sigma_\textrm{d,sm} + \varkappa_\textrm{gr} \Sigma_\textrm{d,gr})$
, where
$\varkappa_\textrm{sm}=10^4$
cm
$^2$
g-1,
$\varkappa_\textrm{gr}=2\times 10^2$
cm
$^2$
g-1 are typical values of absorption coefficients in the UV for small and grown grains (Pavlyuchenkov et al. Reference Pavlyuchenkov, Akimkin, Wiebe and Vorobyov2019, Fig. 1).
We calculate the adsorption rate
$\lambda$
following Brown & Charnley (Reference Brown and Charnley1990). It is proportional to the total cross-section of dust grains per unit volume, which can be obtained from the total surface area of dust grains per unit disc surface
$\widetilde{\sigma}_\textrm{tot}$
. To change the normalisation to the 2D case,
$\widetilde{\sigma}_\textrm{ tot}$
needs to be multiplied by
$1/2H$
. Another factor of
$1/4$
follows from the difference between cross-section and surface area of a sphere. As a result, the adsorption rate is calculated as
\begin{align}\lambda= \frac{\widetilde{\sigma}_\textrm{tot}}{8H} \sqrt{\frac{8k_\textrm{B}T_\textrm{mp}}{\pi \mu_\textrm{sp} m_\textrm{ p}}}. \end{align}
A more detailed derivation of rate coefficients for adsorption and desorption is presented in Section 2.3 of Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021).
Table 1 summarises the molecular parameters used in both models in this work. Binding energies for H
$_2$
O, CO
$_2$
, and CO are based on the experimental data from Cuppen et al. (Reference Cuppen, Walsh, Lamberts, Semenov, Garrod, Penteado and Ioppolo2017) for desorption from crystalline water ice. The binding energy for methane is taken from Aikawa et al. (Reference Aikawa, Miyama, Nakano and Ume-bayashi1996). The values of the initial abundances relative to water
$f_\textrm{s}$
are based on the observations of ices in protostellar cores around the low-mass protostars (Öberg et al. 2011a). They are transformed to the mass fraction relative to gas assuming the water abundance of
$5\times10^{-5}$
relative to gas number density. Total initial mass of ices comprises
$\approx 8.5\%$
of the total mass of refractory grain cores. This is relatively low compared to the estimates suggesting comparable masses of ices and refractories in the discs (Pontoppidan et al. Reference Pontoppidan, Salyk, Bergin, Brittain, Marty, Mousis, Öberg, Beuther, Klessen, Dullemond and Henning2014). However, the lower fraction of ices is more suitable in our approach, that suggests that ice mantles do not change mass and radius of dust grains (Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). As we are mostly interested in the elemental ratios and relative abundances of the considered ices, lower ice fraction is an appropriate simplification. However, we note that dust dynamics can lead to accumulation of ices and ice mantles exceeding masses of silicate cores in some disc regions, as was shown previously in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021).
Binding energies, molecular weights, and initial abundances for the considered volatiles adopted in the modelling. Initial abundances of the species
$f_\textrm{s}$
are shown relative to number density of water molecules in ice phase, and
$\Sigma_{s}^\textrm{sm}/\Sigma_\textrm{g}^\textrm{init}$
is the corresponding initial mass fraction of the ices relative to gas surface density.
Ice mantles also provide feedback on dust evolution. The model includes the effect of ices on fragmentation velocity
$v_\textrm{frag}$
, which is the the maximum collision velocity leading to sticking instead of fragmentation. According to laboratory experiments, icy grains have higher
$v_\textrm{frag}$
than bare silicate grains by an order of magnitude (Wada et al. Reference Wada, Tanaka, Suyama, Kimura and Yamamoto2009; Gundlach & Blum Reference Gundlach and Blum2015). In Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021) we used the values of fragmentation velocity
$v_\textrm{frag}=1.5$
and 15 m s-1 for bare and icy grains, correspondingly. Here we follow Okuzumi & Tazaki (Reference Okuzumi and Tazaki2019) and adopt lower values of
$v_\textrm{frag}=0.5$
and 5 m s-1, which are more relevant for grains consisting of
$\mu$
m-sized monomers. These lower values of
$v_\textrm{frag}$
will lead to higher importance of fragmentation compared to Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). To determine if a dust grain should be considered icy or bare, we compare the local total surface density of all ices on grown dust divided by
$\Sigma_\textrm{d,gr}$
with the threshold value K, which is calculated as
\begin{align} K = \frac{3 a_\textrm{ml}\rho_\textrm{ice}}{\sqrt{a_*a_\textrm{max}} \rho_\textrm{s}}, \end{align}
i.e. an icy grain must have at least one monolayer of ice. Here,
$a_\textrm{ml}$
is the thickness of the ice monolayer estimated as the size of a water molecule
$3 \times 10^{-8}$
cm. The material density of ice
$\rho_\textrm{ ice}=1$
g cm-3 and the mean radius of a grown grain is calculated as
$\sqrt{a_*a_\textrm{max}}$
for the power-law distribution with
$p=3.5$
.
3. Results
To understand the distribution of the species, we first need to consider the global evolution of the disc and its structure. The distribution of volatiles and the C/O ratio is very sensitive to gas and dust substructures appearing in the disc. Variations in temperature and pressure lead to the complex shape and temporal evolution of the snowlines. The dependence of dust fragmentation velocity on the presence of ice mantles implies the feedback from the volatiles on dust and (through back-reaction) on gas.
Our modelling starts with the gravitational collapse of a flattened, slowly rotating molecular cloud. The protoplanetary disc is formed after the formation of the protostar, when the in-spiraling layers of the contracting cloud hit the centrifugal barrier near the inner edge of the sink cell, at a time instance depending on the initial core mass. The disc and the protostar are formed
$\approx53$
kyr after the beginning of the cloud collapse in model M1 and at
$\approx78$
kyr in model M2. If not stated otherwise, times are specified counting from the beginning of the simulation, e.g. the 100 kyr time instance for model M1 refers to a
$\approx50$
kyr old disc, as the first stage includes cloud collapse as well.
An important characteristic of young stellar objects is their variable accretion rate. Our modelling produces accretion bursts with the magnitude of
$\sim100$
L
$_{\odot}$
occurring every
$\sim10^{4}$
yr during the first hundred thousand years of disc evolution. These burst parameters are in line with the episodic accretion scenario (Hartmann & Kenyon Reference Hartmann and Kenyon1985) and resemble the observed phenomenon of FU Ori type stars (see Audard et al. Reference Audard, Beuther, Klessen, Dullemond and Henning2014, for a review). The luminosity outbursts heat up the disc and typically shift the snowlines further away from the star. Although this effect is temporary, it can be reflected in the distributions of the volatiles, and leave long-term imprints in the observed dust properties (Vorobyov et al. Reference Vorobyov2022). The detailed analysis of the effect of such outbursts on the distribution of the elements is worthy of a separate study and lies beyond the scope of this paper.
3.1 Dust and gas structures
During the first hundred thousand years of evolution, protoplanetary discs change from highly asymmetric and dynamic objects to nearly axisymmetric and slowly evolving structures. Fig. 1 shows the examples of different gas and dust substructures that are present in the disc at different stages. The snapshots are shown for model M1, they include a young disc (160 kyr), an intermediate state (350 kyr), and a more evolved and axisymmetric disc (490 kyr). Each of the selected time instances represent some characteristic morphological features addressed below. In model M2, similar structures appear, although sometimes at different evolution times. In this subsection, we consider model M1 as an example and discuss these features, highlighting the properties of dust and gas substructures that are most relevant for the distribution of the volatiles.
The earliest phases of disc evolution are characterised by a large-scale spiral structure in both dust and gas, as well as episodic appearance of clumps. They are the result of gravitational instability (GI) in a massive disc, as in our modelling, the disc mass comprises more than
$0.1$
of the stellar mass, which is roughly a threshold of the disc stability against GI (Vorobyov Reference Vorobyov2013; Kratter & Lodato Reference Kratter and Lodato2016). This is illustrated by the Toomre Q-parameter (Toomre Reference Toomre1964) in the third column of Fig. 1: inside the spirals and the clump
$Q \lt 1$
, which indicates the dominance of self-gravity over Keplerian sheer and gas density. The spiral structures become less prominent with time as the disc loses mass and the Q-parameter increases. However, spirals persist in the model throughout the disc evolution up to 0.5 Myr. For example, at 350 kyr, a very tight spiral is present in the gas, starting at the gas and dust ring at
$\approx10$
au. At 490 kyr, the spiral pattern is weak and exists at
$ \gt 10$
au distances, which are not displayed in Fig. 1.
One- or two-armed grand design spirals are common
$100-200$
kyr after the disc formation in the models. The analogues of such spirals in the observed young protoplanetary discs around low-mass stars are found, for example, in Elias 2-27 or WaOph 6 (Pérez et al. Reference Pérez2016; Huang et al. 2018b). It is not yet clear if the observed spirals are the result of GI or caused by a perturbation from a companion planet or a (sub)stellar object (Meru et al. Reference Meru, Juhász, Ilee, Clarke, Rosotti and Booth2017; Brown-Sevilla et al. Reference Brown-Sevilla2021). A spiral with multiple clumps produced by GI was recently observed in the disc around a FUor V960 Mon (Weber et al. Reference Weber2023).
Another important feature of gas and dust spatial distributions is ring-like structures at various scales. The system of rings starts to form as early as 100 kyr, and develops to the high-contrast multiple ring structure (1–3 orders of magnitude difference between surface densities in rings and gaps), which is evident in the middle and bottom rows of Fig. 1. While gas rings are also common, the annual structures are more prominent in the dust surface density, as well as in dust size. Some of the dust rings correspond spatially to the gas rings, while others are barely reflected in the gas distribution. Overall, the difference between the gas and dust structures develops with time, as the result of dust growth and drift (Testi et al. Reference Testi, Beuther, Klessen, Dullemond and Henning2014). Only some particular substructures, such as the ring between 1–2 au, are present in the distributions of both gas and dust components.
Another location where prominent rings form in both dust and gas is the water snowline. Here, the water snowline is defined as the location in the disc midplane where the amount of water in the gas equals to its amount in the ice (on both dust populations). There are multiple location in the disc where this happens. Water is frozen in most of the disc beyond 5–10 au, and we will refer to the furthest snowline dividing these outer frozen region from the inner one with the gas-phase water as a primary snowline. Generally, the primary snowline is roughly circular, but it can have asymmetries due to the spiral structure and an additional snowline may appear, e.g. around a gravitationally bound clump (upper rows in Fig. 1). Besides, at
$\approx260$
kyr, another region rich in water ice appears in the inner disc, creating a pattern of double or triple water snowline at later times (middle and lower rows in Fig. 1). We will refer to these inner additional snowlines as the secondary snowlines. They circumcise a cold and dense gas-dust ring that forms at 1–2 au. As the disc cools down with time, the primary snowline position moves from
$\approx10$
au distance at 160 kyr, to
$\approx6$
au at 350 kyr and
$\approx4$
au at 490 kyr.
Snowlines are known to be associated with the enhancement of dust and volatiles (Stevenson & Lunine Reference Stevenson and Lunine1988; Cuzzi & Zahnle Reference Cuzzi and Zahnle2004; Drążkowska & Alibert 2017; Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). Dust enhancement was also shown to affect the distribution of gas and its accretion rate through the disc (Gárate et al. Reference Gárate, Birnstiel and Stammler2020) by means of dust back reaction, which is also accounted for in our modelling. In our models, an about 2 au wide dust ring is formed at the inner edge of the water snowline (at
$\approx$
10 au) as soon as 50 kyr after the disc formation. Grown dust grains drifting towards the star through the snowline lose their mantles, their fragmentation velocity drops, rendering them more vulnerable to fragmentation. Consequently, the grain maximum size
$a_\textrm{max}$
decreases, their drift towards the star slows down, which leads to the accumulation of grown dust, as well as small dust as a product of fragmentation. Increase of total dust density also leads to less efficient cooling and results in higher temperature inside the snowline (see the right panels in Fig. 1). Later, at times
$ \gt 200$
kyr, several additional rings form outside the water snowline at distances up to 100 au, the most notable one being
$1-2$
au exterior to the primary snowline. Dust is trapped inside the gas pressure maxima in these rings, which is a self-supporting phenomenon as the temperature also increases inside the dust ring due to high optical depth.
These rings are worthy of attention in the context of possible planet formation. Dust surface density and size are higher in the rings, with
$a_\textrm{max}$
reaching centimetres, dust-to-gas ratio up to
$0.1-0.2$
, and Stokes number up to
$0.05-0.1$
. This could ease the development of the streaming instability (SI), which typically requires pebble-size grains with
$\textrm{St}\gtrsim 0.01$
and dust-to-gas ratio
$\gtrsim0.02$
(Carrera, Johansen, & Davies Reference Carrera, Johansen and Davies2015). Multiple ring-like structures are commonly observed in protoplanetary discs at a range of ages and display a variety of examples, with different widths, contrasts and numbers of rings (Huang et al. 2018a, and many others). However, to directly compare the ring structures in the simulated dust surface density with the observed dust emission, require radiative transfer modelling is needed. Some of the observed ring structures could be produced by radial variation in dust size even in the absence of gaps in dust surface density (Akimkin et al. Reference Akimkin, Vorobyov, Pavlyuchenkov and Stoy-anovskaya2020).
Radial distribution of azimuthally averaged surface densities of the volatiles in the gas and in the ice at various time instances in M1 (
$M_\textrm{core}=0.66$
M
$_{\odot}$
). Pale lines indicate the total surface density of species. The upper panels show surface densities of gas, small dust and grown dust, and the midplane temperature.
Same as Figure but for model M2 (
$M_\textrm{core}=1$
M
$_{\odot}$
).
The most prominent dust rings are located in the vicinity of the water snowline: the ring outside the primary snowline at 5–8 au (depending on the time) and the ring at 1–2 au, inside the primary snowline, which at later times also contains water ice and additional snowlines. The main effect of the snowlines is the change in fragmentation velocity between mantled and non-mantled grains: dust size sharply decreases by
$\approx$
2 orders of magnitude for the latter. Immediately outside the water snowline, the midplane temperature is lower, due to lower dust opacity and thus more efficient cooling. Turbulent
$\alpha$
is on the contrary, higher, providing more efficient radial transport of matter. It leads to lowering the gas surface density in this gap, which in turn increases
$\alpha$
(see Equation 5), creating the positive feedback and further deepening the gap. One of the possible mechanisms to create the initial decrease in the gas surface density is dust back-reaction, which can affect the inward flow of gas at the snowline. This effect was investigated by Gárate et al. (Reference Gárate, Birnstiel and Stammler2020) for different initial dust-to-gas ratios. The radial variation of
$\alpha$
itself lead to the appearance of gas substructures (Tong, Alexander, & Rosotti Reference Tong, Alexander and Rosotti2024). The combined effect of lower temperature and density creates a pressure minimum, which dust tends to avoid.
The dead zone, where
$\alpha$
-parameter values are lower than
$10^{-3}$
, includes the ice-free inner disc and has an approximately two times larger radial span than the iceless region. The distribution of
$\alpha$
-parameter is shown in the right column of Fig. 1. Inside the primary water snowline, the values of
$\alpha$
are the lowest due to high surface density of gas. The dead zone is not axisymmetric and reflects the spiral structures of the gas distribution, as it is sensitive to
$\Sigma_\textrm{g}$
(see Equation 5). The spiral arms of the dead zone span to 15–25 au at 160 kyr, to 10–18 au at 350 kyr and to 10–14 au at 490 kyr. The comparison with the temperature distribution in Fig. 1 indicates that
$\alpha$
and
$T_\textrm{mp}$
are anticorrelated, as higher viscosity provides faster accretion, hence lower surface densities and more efficient cooling. Similarly, a lot of dust accumulates in the dead zone, increasing opacity, which hampers cooling.
The icy dust ring at 1–2 au is especially interesting in the context of planet formation. In this ring, dust-to-gas ratio exceeds 0.1, and surface densities of both gas and dust are increased by more than an order of magnitude compared to the adjacent regions. Due to the ice mantles that protect dust from fragmentation, the dust size reaches several cm, close to the values behind the water snowline. When the ring is established, it is self-supporting, in a sense that without external perturbation (e.g. sharp change in accretion flow from the outer disc), it can be stable for a long time, over tens of kyr.
The presence of water ice inside the dead zone was investigated by Vallet et al. (Reference Vallet, Childs, Martin, Livio and Lepp2023). They showed that in the discs around lower-mass stars, the turbulent heating in the inner disc can be low enough to allow the existence of ices. In our modelling the cooling of the inner region is assisted by the inner dust ring. The freeze-out has a positive feedback on dust growth due to higher
$v_\textrm{frag}$
, which helps dust grow larger and further accumulate towards the pressure maximum in the ring. So in our modelling ices coexist with a lot of grown dust in the dead zone. Such an icy inner region appears to be a promising place for the formation of the volatile rich planets.
3.2 Distribution of volatiles
Even in the absence of chemical reactions, over the years of disc evolution the distribution of volatiles significantly changes compared to the initial one. This is the result of both phase transitions and dust growth and advection. Dust drift brings ices from the outer disc, enriching the inner disc with volatiles. Snowlines provide the conditions favouring accumulation of ices and gas-phase species. The formation of the established disc structures, such as dust and gas rings and spirals (see Section 3.1) leads to a complex pattern of intermittent snowlines. In this Section, we describe the main features of the molecular distributions and their implications for the composition of dust and gas at early stages of protoplanetary disc evolution.
Figs. 2 and 3 show the azimuthally averaged radial profiles of the volatile species in models M1 and M2, respectively. As mentioned earlier, we define the snowline position as the location in the disc midplane where the amount of species in the gas equals to its amount in the ice (on both dust populations). This definition allows the snowline to have complex shape, characterised by different positions for each azimuthal direction. In the azimuthally averaged distributions shown in Figs. 2 and 3, the position would only be approximate. The slope of the species distribution in the vicinity of the snowline reflects the degree of the axial asymmetry in the disc. Flatter distributions appear as the contributions sum up from the snowlines, the radial position of which vary at different azimuthal angles
$\phi$
.
In our model setup, the volatiles can either have no snowline in the disc, or have two or more snowlines depending on the local conditions. Snowlines are absent for more volatile species (CO and CH
$_4$
) at earlier times or during bright luminosity outbursts, when the disc is too hot for them to freeze out. When there are two snowlines, the inner snowline in the disc is the one determined by thermal desorption. We refer to it as the primary snowline. The outer snowline is determined by photo-desorption, it lies in the embedding envelope outside of the body of the disc where optical depth is low. It must be noted that the gas-phase species outside this photo-snowline are vulnerable to photo-dissociation by the UV radiation. This process is not explicitly included in our chemical model, but can be assessed as in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021).
More than two snowlines appear when the disc physical structure becomes more complex, mainly due to the presence of the ring-like structures. Species can freeze inside a cold dense dust ring, creating additional secondary snowlines, also governed by thermal desorption. The concepts of primary and secondary snowline is necessary for H
$_2$
O and CO
$_2$
, which have multiple snowlines at the later stages of disc evolution. For water, the inner icy region appears at
$\approx$
1 au at 490 kyr (see right column of Fig. 3) inside a dense dust ring. For CO
$_2$
, the ring that appears at 7 au, outside the primary water snowline at 4 au, creates the inverted thermal structure in the region with
$T_\textrm{mp}$
close to the typical CO
$_2$
sublimation temperature of 70–90 K. Increase in
$T_\textrm{mp}$
in these ring is caused by higher optical depth and consequently lower cooling on the viscously heated midplane.
Apart from the snowlines, the distributions of volatiles in Figs. 2 and 3 present local radial variations in all of the species components, including total abundance of the species. The initial total abundance is kept only in the envelope. As the matter is being redistributed, abundances of all volatiles in the inner disc grow. Particularly, the process responsible for this is dust drift. It brings the grown ice-covered grains to the inner disc regions, where their ice mantles are sublimated and no longer move with the drifting silicate dust cores. The effect is stronger for less volatile species H
$_2$
O and CO
$_2$
: their abundance in the gas grows by a factor of a few inside their primary snowlines. For CO and CH
$_4$
the effect is weaker, because there is less grown dust outside their snowlines, and those snowlines are not very much established at earlier times. Thus their abundances inside the snowline only grow by a factor of unity, except for the immediate vicinity of the snowline. Moreover, both CO and CH
$_4$
have a bump in the gas-phase distribution just outside the dust ring at 1 au, that is absent in H
$_2$
O and CO
$_2$
. At the inner disc edge the abundances of CO and CH
$_4$
in the gas are lower than the initial value.
The abundance enhancements appear most notably at the snowlines, with local bumps in all three phases at later times (after
$\approx$
350 kyr). They are produced by the combination of the dust radial drift and the azimuthal oscillations of dust and gas radial velocity, described in more detail in Molyarova et al. (Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021) and Molyarova et al. (in preparation). By 500 kyr, the surface density of gas-phase H
$_2$
O exceeds the initial value by a factor of 5, of CO
$_2$
– by a factor of 7, of CH
$_4$
– by 3.5, and of CO – by 2.5. Similar accumulation powered by turbulent diffusion was previously studied for CO molecule in axially symmetric model setup (Stammler et al. Reference Stammler, Birnstiel, Panić, Dullemond and Dominik2017; Krijt et al. Reference Krijt, Schwarz, Bergin and Ciesla2018). Their results indicated similar enhancement in the gas phase by a factor of a few. In our non-axisymmetric approach, diffusion is not a necessary requirement, and the necessary transport is rather provided by azimuthal variations of radial velocity induced by disc self-gravity (Molyarova et al. in preparation).
Distributions of ices on small and grown dust are different as they are affected by dust growth and drift. In general, there is more grown dust than small by mass, and in most of the disc there is more ices on grown dust, particularly at later times. However, in the outer disc and at the earlier times, ices on small grains dominate. As initially all ices are on small dust, it seems inevitable that they will gradually move to the grown grains as small dust coagulates and turns into grown grains. However, near all the snowlines, amount of ices on small dust increases due to the effect of the spiral pattern. Ice abundances are determined by episodic sublimation and freeze-out in the warm spiral arms of the complex density and temperature pattern (see Section 3.3.2. in Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). As adsorption preferably happens to the smaller grains due to their larger total surface area, there is much more ices on small dust in the wake of the spiral arms. This concerns the photo-desorption snowlines as well.
3.3 C/O ratios
Here we consider the relative amount of carbon and oxygen in the gas, on the surface of small and grown grains, and in total for all for all phases and disc locations. C/O ratio is seen as a perspective tracer of planet formation mechanisms (see Öberg et al. 2011b; Thiabaud et al. Reference Thiabaud, Marboeuf, Alibert, Leya and Mezger2015). Therefore, we are especially interested in the regions of the disc and the volatile phase component where the C/O ratio declines from the total initial value (see below). Particularly, we are interested in C/O ratio noticeably above the initial value, as it was observed in atmospheres of exoplanets and suggests the formation of these planets in similarly carbon-enriched environments (Swain et al. Reference Swain2009; Madhusudhan et al. Reference Madhusudhan2011; Moses et al. Reference Moses, Madhusudhan, Visscher and Freedman2013; Facchini et al. Reference Facchini, Teague, Bae, Benisty, Keppler and Isella2021).
The C/O ratio is calculated as the total number of carbon atoms contained in the molecules, divided by the total number of oxygen atoms in the molecules. This calculation can be done separately for the species contained in the gas phase, species on dust surface, or for the species in all phases. Thus, we calculate the C/O ratio in the gas, in the ice, and total C/O, respectively. For the ice species, we include both ices on grown and small dust grains.In some disc regions, the amount of carbon and oxygen in in a particular phase is very low, e.g. in the ice phase inside the water snowline, where all the molecules are in the gas. For such cases, C/O ratio would not be representative of the chemical composition, so we exclude the corresponding computational cells from the consideration. We only display the C/O ratio if the mass of the volatiles in the phase is higher than
$10^{-3}$
of the total mass of volatiles at a given location.
For the adopted molecular abundances based on Öberg et al. (2011a) and also used by Eistrup et al. (Reference Eistrup, Walsh and van Dishoeck2018), the baseline value is
$\textrm{C/O}=0.34$
, which is in line with the gas-phase abundances in the ISM (
$\approx0.36$
, Przybilla, Nieva, & Butler Reference Przybilla, Nieva and Butler2008). This value is different from the typical solar value of
$\approx0.5$
(Przybilla et al. Reference Przybilla, Nieva and Butler2008), which is commonly used as a reference (e.g., Öberg et al. 2011b; Semenov et al. Reference Semenov2018). First, local galactic abundances changed since the Solar system formation 4.6 billion years ago (see, e.g. Appendix A in Bergin et al. Reference Bergin, Booth, Colmenares and Ilee2024, for the comparison). Second, the difference also stems from inclusion or non-inclusion of the dust component. The stellar atmospheric abundances comprise all the elements present in the medium where the star formed, while the elements available for the volatiles are only a fraction of that. The values of the elemental abundances can be determined separately for the volatile (gas and ices) and the refractory (rocky dust grain cores) components of the ISM (Hensley & Draine Reference Hensley and Draine2021). Here, we do not take into account carbon and oxygen from the rocky cores of dust grains. Due to the model assumptions, the initial ice-to-rock mass ratio in our model is quite low (0.08 compared to 2–4 in Pontoppidan et al. Reference Pontoppidan, Salyk, Bergin, Brittain, Marty, Mousis, Öberg, Beuther, Klessen, Dullemond and Henning2014). Therefore, adding the elements from solid dust grains to those contained in ices would be misleading, as the former would dominate in the resulting C/O ratio. In this work, we concentrate on considering the C/O ratios of the volatile component (ice and gas), and compare them with the initial value of 0.34.
Fig. 4 shows radial profiles of the C/O ratios averaged over the azimuthal angle
$\phi$
by the end of the simulations, at 490 kyr. The key values of C/O ratio relevant in the context of planet formation are the initial value of
$0.34$
, motivated by the comparison with the initial abundances, and
$1.0$
, which is a boarder between carbon- and oxygen-dominated chemical regimes in the ISM and (exo)planetary atmospheres (see, e.g. Tielens Reference Tielens2005; Seager et al. Reference Seager, Richardson, Hansen, Menou, Cho and Deming2005; Madhusudhan Reference Madhusudhan2012).
Radial profiles of the C/O ratio at 490 kyr in models M1 (left) and M2 (right). The plots show C/O in total (black), in the gas (red), and in the ice (blue). The C/O ratio in the ice (gas) is only shown for radial distances where the mass of the volatiles in the ice (gas) is larger than 0.1% of the total mass of the volatiles in the gas (ice). The grey horizontal line indicates the baseline
$\textrm{C/O}=0.34$
. Positions of the snowlines are highlighted by vertical dashed lines. The regions where water (thus all other species) is in the gas are shaded with purple.
By the age of 490 yr, the C/O ratio in both models is significantly different from the initial value. This concerns the total C/O ratio, as well as the C/O ratio in the gas and in the ice. Let us analyse the distributions of C/O ratios that we can see from Fig. 4, moving inwards from the envelope to the centre. The main features of the C/O distributions are the following:
-
• the C/O ratio in the envelope is close to the initial value
$0.34$
in the gas and in total, and it grows up to unity from the envelope to the CO snowline in the disc; -
• around the CO snowline, the C/O approaches unity;
-
• between CO
$_2$
and CH
$_4$
snowlines the C/O in the gas is
$ \gt $
1, and the C/O in the ice is below the initial; -
• at the distances of tens of au, there are variations in the total C/O ratio that are not connected with any snowlines, or variations in gas- or ice-phase C/O;
-
• ice-phase C/O ratios peak around all snowlines except water;
-
• inside the primary water snowline is the region rich in volatiles, with both gas- and ice-phase C/O are the lowest in the disc.
As expected, the key changes in the C/O ratio distribution are associated with the positions of the snowlines. There are two main processes. First, the freeze-out and desorption at the snowlines transfer the elements between phases, altering the C/O in the gas and in the ice. Second, the snowlines favour the accumulation of the respective volatiles both in the gas and in the ice, as was discussed in Section 3.2, pumping up the amount of both components and altering their proportion. Besides the snowlines, radial drift of grown grains (Weidenschilling Reference Weidenschilling1977) transports the volatiles inwards. We discuss below which processes are responsible for the formation of the listed features.
C/O in the outer disc. In the surrounding envelope, outside the disc, all volatiles are in the gas due to photo-desorption, and C/O in the gas is close to the initial value of 0.34. Around the CO snowline and beyond, the C/O in the gas is close to unity, and in the gas C/O is higher that the initial, which requires explanation. The region beyond the CO snowline is usually described as the place where all species are frozen, thus having a stellar C/O ratio in the solid phase and practically no carbon or oxygen in the gas (e.g. Öberg et al. 2011b; Öberg &Wordsworth Reference Öberg and Wordsworth2019; Mollière et al. Reference Mollière2020). In our simulations, only in model M2 there is a region where less than
$10^{-3}$
of C and O is in the gas, and in model M1 such region is absent. Due to the asymmetric spiral structure that persists even at 490 kyr, even though most of CO is frozen beyond the snowline, there is a significant (
$ \gt 10^{-3}$
) fraction of it in the gas. Additionally, the position of CO snowline itself is significantly affected by the dust drift, as it declines from the equilibrium between adsorption and desorption.The indistinctness of the CO snowline also helps CO to persist in the outer regions: all other species are ultimately frozen, so they are efficiently carried away to the inner disc via radial drift, while this works worse for only partially frozen CO. It creates relative overabundance of CO in the outer disc, elevating the total C/O ratio and later the ice-phase C/O ratio, when the preserved CO freezes out. Between the CO ice line and the envelope, there is a gradient of C/O in the gas due to photo-desorption of the ices. The last molecule to be photo-desorbed is water, which returns oxygen to the gas phase at the farthest radial distance.
$C/O\approx1$
around the CO snowline. Beyond the CH
$_4$
snowline, CO dominates the composition of volatiles in the gas phase, leading to the gas-phase C/O close to unity, the value characteristic of the CO molecule. At the CO snowline, the CO dominates the ice-phase composition as well. While dust drift substantially lowered the abundances of CO
$_2$
and H
$_2$
O by 490 kyr in these regions, CO ice accumulated at the snowline. This leads to the ice-phase C/O closer to 1, too, making the vicinity of the CO snowline a region where total amounts of C and O are similar.
High C/O in the gas, low C/O in the ice. In the disc regions beyond the primary CO
$_2$
snowline, only CO and CH
$_4$
are in the gas, meaning the dominance of carbon and C/O
$ \gt 1$
. Consequently, the C/O in the ice is generally lower, around
$0.2$
, as the ices are mainly H
$_2$
O and CO
$_2$
rich in oxygen. This is consistent with the classical step-like picture of Öberg et al. (2011b), with the addition of carbon-rich methane allowing the
$\textrm{C/O} \gt 1$
. The midplane C/O ratios we simulate are difficult to directly compare with observations, which mostly trace the molecular layer above the disc midplane. High C/O ratios in the gas are indeed observed in many protoplanetary discs (e.g. Miotello et al. Reference Miotello2019), but they are considered to be a natural consequence of dust settling (Krijt et al. Reference Krijt, Schwarz, Bergin and Ciesla2018, Reference Krijt, Bosman, Zhang, Schwarz, Ciesla and Bergin2020). Dust inward drift could also enhance this effect in the outer disc regions, creating the radial gradient of total C/O ratio in the disc. Observations of CS and SO emission coming from close to the midplane layers potentially indicate the presence of such radial gradient of the gas-phase C/O in the PDS 70 disc (Rampinelli et al. Reference Rampinelli2024).
Variations of total C/O. Throughout the disc, there are sharp changes in total C/O ratio not connected with any snowline. Most noticeable are the variations between CO
$_2$
and CH
$_4$
snowlines, where gas-phase and ice-phase C/O ratios are stable. These variations are associated with the disc substructures, particularly with the dense dust rings described in Section 3.1. The total C/O changes due to radial variations in dust-to-gas ratio: when it is higher, the total C/O is closer to the ice-phase C/O, and vice versa. Inside the dust-rich rings, H
$_2$
O and CO
$_2$
ices are abundant, due to high surface density of dust relative to gas. At the same time, CO and CH
$_4$
in the gas phase have similar surface densities inside dust rings and between them. Thus, in the dust rings, CO
$_2$
and water are overabundant, leading to lower total C/O ratio. We consider this effect in more detail in Section 3.5. Variations of the total C/O ratio due to dust substructures are also present in the cold ring at
$\approx1$
au.
C/O peaks at the snowlines. There are peaks of the C/O ratio in the ice right outside the snowlines of CO
$_2$
, CH
$_4$
and CO, produced by the accumulation of the respective ices (see Figs. 2 and 3). In M2 model, the amount of CH
$_4$
and CO ices at their respective snowlines becomes comparable with or even larger than those of CO
$_2$
and H
$_2$
O (compare the right panels of Fig. 3), leading to C/O in the ice
$\approx$
0.6–0.9. In model M1, the accumulation is more prominent, so the C/O ratio in the ice approaches unity at CO snowline and
$ \gt 1$
at the methane snowline. These peaks distort the pattern of generally low C/O ratio in the ice and preset additional regions where carbon-rich planetesimals could be formed, and carbon-rich pebbles could be accreted onto forming protoplanets. At the snowlines of H
$_2$
O and CO
$_2$
, the C/O in the ice approaches the C/O ratios of these molecules, 0 and 0.5, respectively.
Lower C/O inside the water snowline. Inside the primary water snowline, the C/O ratio is generally lower than outside of it, close to the initial 0.34. Contrary to the outer regions depleted of ices due to the radial drift, the disc parts inside and around the water snowline are enriched in volatiles, and particularly of oxygen-rich water and CO
$_2$
. Total and gas-phase C/O ratios vary from
$\approx0.2$
to
$0.6$
, as the secondary water snowlines add more substructure to the C/O distribution. The regions where the only ice is water and the ice-phase
$\textrm{C/O}=0$
are the inner cold dust ring at 1 au and the narrow annuli between water and CO
$_2$
snowlines. The enrichment of the inner disc regions with oxygen as a result of dust radial drift is suggested by the resolved observations of molecules in protoplanetary discs (Banzatti et al. Reference Banzatti2020).
These characteristic features of the C/O distributions are similar in the two presented models. The main difference is the radial distances where the borders between the zones are located; they are closer to the star in the less massive and thus colder model M1. This is mainly due to slightly different masses of the central star accumulated throughout 490 kyr of non-identical protostellar accretion history, which lead to different luminosity and thermal structure (see upper panels of Fig. 5). Particularly, the stellar masses and luminosities at this time instance are:
$1.07$
L
$_{\odot}$
and
$0.34$
M
$_{\odot}$
(for M1), 1.89 L
$_{\odot}$
and 0.58 M
$_{\odot}$
(for M2). The less massive model M1 demonstrates overall higher C/O ratio in both phases.
Evolution of central source luminosity and C/O ratio in models M1 (
$M_\textrm{core}=0.66\,{\rm M}_{\odot}$
, left) and M2 (
$M_\textrm{core}=1\,{\rm M}_{\odot}$
, right). The upper panels show stellar and accretion luminosity depending on time. Below are, successively, total C/O ratio, C/O in the gas, and C/O in the ice, depending on time. The C/O values above and below the initial value of 0.34 are coloured in shades of red and blue, respectively. The regions with low abundances of both carbon and oxygen, either in the gas or ice phases, are shown in white. Coloured contours correspond to the positions of the snowlines. Photo-dissociation snowlines are not shown.
3.4 Evolution of the snowlines and the C/O ratios
The positions of the snowlines are crucial for the values of C/O ratios, both because of the direct change through the phase transitions and the associated accumulation of volatiles. They depend on the local gas and dust properties, particularly on temperature. They can also be shifted inward due to dust drift (Piso et al. Reference Piso, Öberg, Birnstiel and Murray-Clay2015). Snowlines evolve as the disc structure changes with time. In this Section, we consider the co-evolution of the snowlines and the C/O ratios and discuss the mechanisms of species redistribution over the disc.
The temporal evolution of the azimuthally averaged C/O ratios and the equilibrium positions of the snowlines is shown in Fig. 5. It is evident from Fig. 5 that the C/O ratio indeed follows the snowlines, particularly inside
$\approx10$
au. The C/O structure changes throughout the disc evolution, and some key features only appear at later times.
One of the key factors affecting the disc thermal structure is the luminosity of the central source. In the upper panels of Fig. 5, we show the evolution of total luminosity, which directly affects the positions of the snowlines. The luminosity is the sum of stellar and accretion components, coming from the protorstar itself and gravitational potential energy of the accreted matter. The stellar luminosity gradually decreases as the protostar becomes more compact. Accretion luminosity depends on the accretion rate, which is highly variable as a result of magnetorotational and gravitational instabilities in the disc (Kadam et al. Reference Kadam, Vorobyov, Regály, Kóspál and Ábrahám2019, Reference Kadam, Vorobyov, Regály, Kóspál and Ábrahám2020; Vorobyov et al. 2020b). The simulated episodic luminosity outbursts are similar to those occurring in the observed YSOs (Connelley & Reipurth Reference Connelley and Reipurth2018), with their amplitudes of tens and hundreds
$L_{\odot}$
. In the more massive model M2, the outbursts are more frequent, brighter and occur until later times. The reasons behind this difference is the massive disc being more prone to both MRI and GI, which needs to be investigated in more detail in a separate study.
Snowlines of the least volatile of the considered species, H
$_2$
O and CO
$_2$
, exist in the model since the earliest phases of the disc formation. Even during bright luminosity outbursts (
$\sim100$
L
$_{\odot}$
) they do not disappear, but move farther away from the star. During the first
$\approx50$
kyr the disc is spreading out, it is highly asymmetric and dynamic, so the snowline positions oscillate. Later on, the disc generally cools down, and the snowlines gradually move towards the star (except during the outbursts). In model M2, between 130 and 500 kyr, the primary water snowline moves from 12 to
$4.5$
au; the primary CO
$_2$
snowline moves from 23 to 7 au. In model M1, the water snowline moves from 9 to 4 au, and the CO
$_2$
snowline moves from 17 to 4.5 au.
Despite a factor of two different binding energies of H
$_2$
O and CO
$_2$
(5 770 and 2 360 K, respectively), the locations of their primary snowlines do not differ much due to the steep radial temperature gradient around these distances (see upper panels of Figs. 3 and 2). Inside 10–20 au,
$T_\textrm{mp}$
is determined by heating mechanisms other than external irradiation: viscous heating, gas work (PdV heating), heating by shocks and energy transport with advection. This also means that the water snowline is less sensitive to the level of irradiative heating, thus only slightly affected by the luminosity outbursts. The temperature change is particularly sharp at the water snowline(s), with the absolute value of the approximated power-law slope of 1.5–6. High surface density of dust in water-ice-free regions makes cooling less efficient and leads to locking up the produced heat and consequently higher temperatures. Besides, both these species have multiple snowlines due to the formation of ring-like substructures with the conditions close to the borderline between their frozen and gaseous state. These additional snowlines also affect the C/O distributions.
Methane and CO are the more volatile species in our model. They either have zero or two snowlines in the disc. The snowlines are absent during the outbursts brighter than
$\approx100$
L
$_{\odot}$
for methane and
$\approx200$
L
$_{\odot}$
for CO. Until approximately 200–250 kyr, there is no established CO snowline inside the disc. For example, at 160 kyr (see lower left panel in Figs. 3 and 2) there is CO ice on both small and grown dust, but their amount is an order of magnitude lower than that of the gas. Additionally, most of this ice is located at the outer disc edge, where gas and dust surface densities sharply drop. Similar distribution appears for CO and CH
$_4$
during bright outbursts: their ices are present in some disc regions, but due to non-axisymmetric disc structure, they do not dominate in the averaged profiles, and there is no common snowline for the whole disc. There are no secondary snowlines of CO or CH
$_4$
, because there is no prominent gas and dust structures in the outer disc where these species are frozen.
Snowlines divide the disc into several zones with different characteristic C/O ratios. However, the chemical composition and C/O ratios in these zones change with time. One of the distinct zones is the region where water is not frozen, shaded with purple in Fig. 4 and circumscribed by the dark purple dotted line in Fig. 5. In this zone, all the species are in the gas phase, so the C/O ratio is initially close to 0.34. However, as the disc evolves, dust brings more volatiles from the outer disc. Particularly, the abundance of water grows most significantly, making C/O in the gas decrease with time. This happens because the main mechanism of the transport of the volatiles is dust radial drift, which works best in the regions where grown grains can sustain their mantles. Banzatti et al. (Reference Banzatti2020) suggest that dust growth and drift can be responsible for the observed anticorrelation between disc radius and H
$_2$
O emission, implying that the inner regions of small discs with are enriched in water brought by efficiently drifting grains. Water is the least volatile species in our model, it is frozen in the largest part of the disc, thus its distribution is most strongly affected by the dust drift. At later times, this zone is divided into two, when a cold dense dust ring forms at
$1-2$
au, which happens at
$\approx460$
kyr in model M2 and at
$\approx270$
kyr in model M1. Interior to the ring, water abundance in the gas is around an order of magnitude lower than outside of it in both models. This happens because the inward flow of gas-phase H
$_2$
O is ‘blocked’ by the cold ring where it freezes and accumulates with the grown grains in the pressure maximum. So the C/O ratio in the gas at
$r\lessapprox2$
au is determined by CO
$_2$
and thus close to 0.5.
The C/O ratio in the ice is not defined in the envelope, where all ices are photo-desorbed, and in the warmest inner disc, where even water is in the gas. In disc regions where only water is frozen, the C/O ratio in the ice is zero. Before 200–250 kyr, when the disc cools down enough for CO to freeze out, the C/O ratio in the ice in the rest of the disc is close to 0.2. It is determined by CO
$_2$
and H
$_2$
O ices, with a small contribution from the low-abundance CH
$_4$
. After that, when CO freezes out and CO snowline appears, the C/O ratio in the outer disc region becomes close to the initial value both in the ice and in total. This region beyond the CO snowline is frequently referred to in relation to giant exoplanets with stellar C/O ratios (e.g. Mollière et al. Reference Mollière2020; Öberg & Wordsworth Reference Öberg and Wordsworth2019; Ohno & Ueda Reference Ohno and Ueda2021). This region would be a perfect location for planets to accrete pebbles covered with icy mantles with the primordial elemental ratio, which would directly become part of the planetary atmosphere. However, pebbles are not initially present in the disc, and their existence in these outer regions is not guaranteed. Pebbles are dust grains large enough to move relative to gas (e.g. Lambrechts & Johansen Reference Lambrechts and Johansen2012; Lenz, Klahr, & Birnstiel Reference Lenz, Klahr and Birnstiel2019), and a fraction of the grown dust in our modelling can be classified as pebbles. The properties of pebbles and composition of their ice mantles in the same setting of the FEOSAD model were studied by Topchieva et al. (Reference Topchieva, Molyarova, Akimkin, Maksimova and Vorobyov2024). They show that pebbles appear in the disc as early as 50 kyr after its formation, and exist in a wide region of the disc. This partially includes the region beyond the CO snowline, but only the area of CO enhancement, where CO dominates in the ice composition and thus the C/O ratio in the ice is close to unity (see lower panels in their Fig. 3).
As shown by Topchieva et al. (Reference Topchieva, Molyarova, Akimkin, Maksimova and Vorobyov2024), ices on pebbles are dominated by H
$_2$
O and CO
$_2$
. In this case, relatively high values of the C/O in the ice (
$\approx0.5$
) correspond to the regions where there is more CO
$_2$
, i.e. around the CO
$_2$
snowlines. This is also the region where a prominent dust ring forms under the influence of the primary water snowline. It is characterised by accumulation of CO
$_2$
and to lesser degree H
$_2$
O ice, as well as vapours, and relatively high amount of grown dust in the ring. The dust ring is situated between the two regions with C/O
$\approx0.34$
in the ice. It presents another favourable location for accreting the ice content with close-to-initial C/O ratio.
The total C/O ratio in the disc also changed significantly from the initial value due to dust drift that redistributes the ices. The matter becomes more carbon-rich as the grains bring CO and CH
$_4$
from the outer disc parts. This effect was previously investigated by Stammler et al. (Reference Stammler, Birnstiel, Panić, Dullemond and Dominik2017) and Krijt et al. (Reference Krijt, Schwarz, Bergin and Ciesla2018) in their modelling of CO dynamics and dust evolution. However, as was shown by Krijt et al. (Reference Krijt, Schwarz, Bergin and Ciesla2018, Reference Krijt, Bosman, Zhang, Schwarz, Ciesla and Bergin2020), vertical settling of grown grains towards the midplane is responsible for depleting the upper layers in the outer disc of gas-phase oxygen, which cannot be captured within our thin-disc modelling. The panels in the second row of Fig. 5 demonstrate strong enhancement of total C/O ratio in the intermediate disc regions. In model M2, the total C/O ratio between 10–100 au becomes
$\approx0.7$
after
$\sim300$
kyr, which is two times higher than the initial value. At the snowlines of CH
$_4$
and CO
$_2$
it approaches unity, mostly due to the accumulation of these species in the gas. In model M1, this process begins
$\approx200$
kyr earlier and consequently leads to even higher total C/O ratio. At 400–500 kyr, most of the disc between 5 and 100 au has total C/O
$\gtrsim1$
in model M2, which demonstrates the powerful impact of dust drift.
A distinctive feature of the produced C/O distributions is the peaks of total and ice-phase C/O ratios around the CO and CH
$_4$
snowlines. In model M2, the increase of the C/O ratios becomes noticeable only after
$\approx$
450 kyr, while in model M1 it starts to form around
$\approx$
300 kyr. Initial abundances of CH
$_4$
and CO are lower than that of CO
$_2$
and H
$_2$
O. As these species accumulate at their snowlines, the abundances become comparable, leading to C/O in the ice
$\approx$
0.6–0.8 in model M2 and up to 1 in model M1. Similar accumulation is seen around CO
$_2$
snowline, but unlike CO and CH
$_4$
snowlines, it is also connected to the interaction with the dust ring structure and is considered in more detail in Section 3.5.
3.5 Two-dimentional structure
The disc is not axisymmetric even at later stages of its evolution (see Section 3.1 and Fig. 1), which is also reflected in the distributions of volatiles and C/O ratios. Examples of 2D distributions of C/O in model M1 at two time instances are shown in upper panels of Fig. 6. The left-hand side group of panels shows the structure that is characteristic of the earlier phases, the same time instance as the upper row in Fig. 1. The prominent spiral structure and a clump both have their reflection in the C/O ratios. The snowlines of CO
$_2$
and CH
$_4$
have clearly non-regular shape affected by the spiral pattern of the gas, with blobs of frozen methane inside the main snowline. Inside the clump, the temperature is higher (see Fig. 1), and both CO
$_2$
and water are in the gas. The separation between gas-phase and ice-phase C/O ratios is clear, but the total C/O ratio in the clump is similar to the surroundings and is only slightly above the primordial value. The gas-phase C/O ratio inside the clump is around
$0.7$
, lower than in the surrounding gas at the same radial distances. Similar decrease of C/O ratio indicated by lowered CS/CO ratio was observed in the disc around DR Tau (Huang et al. Reference Huang, Bergin, Le Gal, Andrews, Bae, Keyte and Sturm2024). Lifetime of the clump (several orbital periods) is too short for significant differentiation between gas- and ice-phase composition to develop, mainly because of insufficient numerical resolution. Focused studies with higher resolution are needed to explore the C/O ratio in the clumps as precursor of giant planets formed via disc fragmentation.
Distributions of C/O ratios and gas/dust surface densities. Model M1 160 kyr (upper left) and 300 kyr (upper right); model M2 250 kyr (lower left) and 490 kyr (lower right). Dotted lines mark the positions of the snowlines for H
$_2$
O (dark purple), CO
$_2$
(magenta), CH
$_4$
(green), and CO (yellow).
The upper right-hand side group of panels in Fig. 6 features a later stage of the disc evolution in model M1. By 300 kyr, the accretion rate and the average positions of the snowlines are stabilised (see Fig. 5). The spiral structure in the gas is much weaker but still present at
$r \gt 10$
au. The CO snowline is established at
$\approx80$
au. The 2D shape of the snowlines is more circular at 300 kyr, particularly for the less volatile CO
$_2$
and H
$_2$
O. The effect of spirals is still evident in the contours of the CO and CH
$_4$
snowlines. At the outer side of these snowlines, at approximately 40 and 80 au, respectively, the C/O in the gas has local maxima, which are also seen in Fig. 4. The one at 40 au also has a clearly spiral shape, repeating the pattern of the gas distribution. The radial span of these peaks is of the same scale as the dispersion of the respective snowline distance from the star. For example, at the CH
$_4$
snowline the C/O peak is
$\approx10$
au wide, and the snowline is at 28–37 au distances. Similar accumulation powered by diffusion of water vapour was shown, e.g. by Drążkowska & Alibert (2017). In addition to diffusion, in our modelling, the species are delivered outward from the snowline due to the dynamic shape of the snowline and two-dimensional movement of gas and dust (see Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021, and Molyarova et al. in preparation).
In more massive model M2, the shape of the snowlines remains complex for longer times. Examples of C/O distributions in M2 are shown in the lower panels of Fig. 6. The radial scale is chosen so that CO
$_2$
and H
$_2$
O snowlines are seen in more detail. At 250 kyr, even the grown dust distribution still has spiral pattern, and CO
$_2$
snowline is following its complex multi-armed shape. The accumulation of CO
$_2$
at the snowline is also very efficient, but it does not strongly affect the C/O ratios, as it drives them to the value 0.5 of the CO
$_2$
molecule, which is close to the intrinsic value. The complex-shaped region between these snowlines has the ice-phase C/O
$\approx 0.7$
. This pattern moves and changes its shape, thus affecting all radial distances between 10 and 20 au. The total C/O ratio is slightly above the ambient value, approaching 0.5, because CO
$_2$
in both phases begins to accumulate in this region.
At later times, the C/O distribution becomes more complex due to the presence of disc substructures, particularly the dust rings. The lower right-hand side group of panels in Fig. 6 show this in more details. First, the temperature and density variations across the dust rings lead to the formation of multiple CO
$_2$
and H
$_2$
O snowlines. An additional annulus of icy CO
$_2$
appears in a relatively cold region between the dense dust rings at 6–7 and 8–13 au. The dust rings are warmer due to higher optical depth and active heating by disk internal sources, and the inner edge of the 6 au ring is warm enough to sustain gas-phase CO
$_2$
. Second, the accumulation of icy dust grains in the rings alters the total C/O ratio. The total C/O ratio anticorrelates with the distribution of grown dust grains: inside the dense rings it is close to the initial value of 0.34, while between the rings, it is higher and reaches
$0.8-0.9$
. At the same time, neither ice-phase, nor gas-phase C/O ratio displays any noticeable variations at 12–25 au, but they do have variations at
$r \lt 12$
au, following the dust ring pattern.
The total C/O ratio is defined by the combination of the ice- and the gas-phase component. Their relative contribution is proportional to the dust-to-gas mass ratio. This is illustrated in Fig. 7. Within the rings, the total C/O is dominated by ices of oxygen-rich species CO
$_2$
and H
$_2$
O, particularly on grown dust accumulated in the pressure maxima. For example, there is a lot of water ice at 6–7 au, and its contribution to the total C/O is weighted with high dust-to-gas ratio of almost
$0.1$
, so the resulting total C/O is lower than the initial value. At the same time, at
$\approx5$
au, dust-to-gas ratio is around
$10^{-4}$
. The dominant species there are in the gas, ice surface densities are
$1-2$
orders of magnitude lower, so the total C/O ratio goes up. In the wide dust ring at 9–13 au, both H
$_2$
O and CO
$_2$
contribute, although at the warmer inner edge of the dust ring CO
$_2$
is sublimated. The value of total C/O is approaches 0.34, also elevated by the presence of gas-phase CO and CH
$_4$
. Beyond
$\approx10$
au, both H
$_2$
O and CO
$_2$
are frozen, and the variations of the total C/O ratio are clearly anticorrelated with the dust-to-gas ratio and the position of the rings. Between the dust rings, where contribution of ices is low, the C/O is determined by CO and CH
$_4$
gases, and reaches values of
$\approx1$
.
Averaged radial profiles of dust-to-gas ratio, the total C/O ratio, and CO
$_2$
and H
$_2$
O in the gas and in the ice in model M2 at 490 kyr. The horizontal lines in the upper panel show the reference values for the C/O ratio (0.34) and dust-to-gas ratio (0.01).
The anticorrelation between the total C/O ratio and dust-to-gas mass ratio is an interesting finding. It is illustrated in Fig. 8 for both models. Only the points within the disc are shown, with
$\Sigma_\textrm{ gas} \gt 0.1$
g cm-2.The pattern obviously changes with time, but the anticorrelation persists. Model M1 demonstrates wider variety of C/O and dust-to-gas ratios. The disc points are grouped in tangled curved ‘branches’, some of them steeper than others. The ‘width’ (or spread) of these branches is determined by the azimuthal substructures. In the axially symmetric parts of the disc the values of dust-to-gas ratio and the total C/O are similar at a given radial distance. Different branches, which can be closer to vertical or horizontal orientation, are the result of the radial variations in the ice-phase C/O ratio and in ice fraction relative to the dust silicate cores. Horizontal branches correspond to weak or absent anticorrelation. For example, near the snowlines of carbon-rich species, C/O in the ice is high and close to that of the gas, which decreases the effect of dust mass fraction on the total C/O. In areas with no ices, the anticorrelation is also irrelevant, because the total C/O is determined entirely by the gas phase. Vertical branches, on the contrary, correspond to the strongest anticorrelation effect, which is expected in the regions between the snowlines, where ice-phase C/O is the lowest. The identified anticorrelation in Fig. 8 is similar to the results of chemical population synthesis modelling (Cridland et al. 2019b, 2020b) showing that the more solids a planet accreted in the disc, the lower the C/O ratio is in its atmosphere.
Dependence between total C/O ratio and dust-to-gas mass ratio in models M1 (upper panel) and M2 (lower panel). Three time instances are shown. The dashed line shows the fitted log-linear dependence for 490 kyr.
We fit the data for
$t=490$
kyr with a linear law (taking the logarithm of dust-to-gas ratio) and obtain the following fits:
$\textrm{C/O} = -0.056 -0.25\log_{10} \left(\xi\right)$
for model M1, and
$\textrm{C/O} = -0.18 -0.27\log_{10} \left(\xi\right)$
for model M2. Here,
$\xi$
is the dust-to-gas mass ratio. The correlation coefficients are
$-0.54$
for M1 and
$-0.57$
for M2. At the shown earlier times (160 and 350 kyr), the correlation coefficient changes between approximately
$-0.9$
and
$-0.5$
, which indicates noticeable anticorrelation throughout the disc evolution.
We note that these C/O ratios only include the volatile component, without the contribution from the refractory material. Although the refractory material is typically considered as silicates, which are rich oxygen, it contains a significant amount of carbon, with the resulting
$\textrm{C/O}\approx0.5$
(see Table 2 in Hensley & Draine Reference Hensley and Draine2021). This solid carbon can be subject to carbon grain destruction (Lee, Bergin, & Nomura Reference Lee, Bergin and Nomura2010; Gail & Trieloff Reference Gail and Trieloff2017; Wei et al. Reference Wei, Nomura, Lee, Ip, Walsh and Millar2019), but this process should be treated separately, as it also affects the gas-phase carbon abundance. Taking refractory cores into account should affect the dependence between total C/O ratio and dust-to-gas ratio, as adding more rock would make the C/O ratio closer to
$0.5$
. Thus the degree of the anticorrelation must be affected by the composition of rocky cores, but the anticorrelation itself should remain even when refractories are included, because the
$\textrm{C/O}\approx0.5$
is still lower than typical C/O of the gas in most of the disc (
$\gtrapprox$
1).
Dust rings are detected in the majority of the observed protoplanetary discs (Long et al. Reference Long2018; Huang et al. 2018a). They are considered as a plausible sites of planet formation (Carrera et al. Reference Carrera, Johansen and Davies2015; Yang, Johansen, & Carrera Reference Yang, Johansen and Carrera2017; Li & Youdin Reference Li and Youdin2021; Lee, Fuentes, & Hopkins Reference Lee, Fuentes and Hopkins2022; Jiang & Ormel Reference Jiang and Ormel2023). The anticorrelation between the total C/O ratio of the volatiles and dust-to-gas mass ratio that we point out is a logical consequence of ices having typically lower C/O ratios and being attached to dust grains. If planets are formed in the dust rings with high dust-to-gas ratios (
$ \gt 10^{-2}$
), either exclusively from solids, or with the inclusion of the dust component, this would imply that their material initially has lower C/O ratio of
$\approx0.5$
and below. To reach higher C/O ratios up to unity and above, which are observed in many exoplanets, these planets would need to migrate and accrete carbon-rich gas from regions other than their immediate formation sites inside the dust rings. In case if planet formation occurs independently of the dust rings, e.g. in the GI, their material is not determined by this anticorrelation.
4. Discussion
Our simulations present a wide range of C/O ratios in the disc in different phases evolving with time. For the atmospheres of giant exoplanets, a variety of C/O ratios were retrieved, too. Here we can compare them to identify the disc regions and times where the chemical and physical conditions for planet formation are consistent. Most of the exoplanets for which the atmospheric composition was retrieved have super-stellar C/O ratios (Hoch et al. Reference Hoch2023; Weiner Mansfield et al. Reference Weiner Mansfield2024), which draws more attention to carbon-enriched areas. They are suggested to form by core accretion, and accreting mostly the gas, which is typically more carbon rich (beyond water snowline). A lot of planets are observed to have stellar or slightly super-stellar C/O (e.g. Mollière et al. Reference Mollière2020; Zhang et al. Reference Zhang2021; Smith et al. Reference Smith2024; Sing et al. Reference Sing2024; Nortmann et al. Reference Nortmann2024, and many others). One way to form such planets is GI, which includes solids and gas together, thus undifferentiated matter is suitable for producing planets with stellar C/O. Disc fragmentation to clumps due to GI requires particular conditions (Meru & Bate Reference Meru and Bate2010; Vorobyov Reference Vorobyov2013), and the direct collapse of gravitationally unstable clumps tends to produce rather massive objects (e.g.
$\approx5$
$M_\textrm{J}$
planets and
$\approx$
60–70
$M_\textrm{J}$
brown dwarfs, see Figure 4 in Vorobyov et al. Reference Vorobyov, Zakhozhay and Dunham2013) at larger radial distances (
$ \gt $
10—100 au, see Vorobyov et al. Reference Vorobyov, Zakhozhay and Dunham2013; Kratter & Lodato Reference Kratter and Lodato2016). GI can also assist the assemblage of planetary cores (Nayakshin Reference Nayakshin2010a,b; Nayakshin, Helled, & Boley Reference Nayakshin, Helled and Boley2014; Vorobyov & Elbakyan Reference Vorobyov and Elbakyan2019). A planet formed through core accretion can also accrete planetesimals, which can be covered with ice, and enrich the atmosphere with oxygen, making the C/O ratio close to the initial stellar value. There are particular exoplanets, where lower than stellar C/O ratio is observed in the atmosphere, such as
$\beta$
Pic b (GRAVITY Collaboration et al. 2020; Worthen et al. Reference Worthen2024), HD 209458 b (Xue et al. Reference Xue, Bean, Zhang, Welbanks, Lunine and August2024), or HD 189733b (Fu et al. Reference Fu2024). Such planets need even more enrichment in ices with low C/O, which makes the regions with low C/O in the ice also more attractive sites for planet formation.
Gravitational instability implies that the planet forms from a mix of gas and dust (Bodenheimer Reference Bodenheimer1974), this is why it is suitable to explain the formation of planets with solar, or unaltered C/O ratios. In our modelling, GI would be associated with the total C/O ratio, which we find to be significantly variable, too. For GI to form a planet with a primordial C/O ratio, it has to occur during the first 100 kyr after the disc formation. At later times, the total C/O ratio changes, and the only region with the primordial C/O ratio is the very outer disc parts, at
$ \gt 100$
au, which is the part of a protoplanetary disc, where conditions for GI are the most consistent with the observed properties of these objects (Rafikov Reference Rafikov2005). Planet formation through GI is indeed more likely at earlier evolutionary stages, when gas surface density is higher (Armitage Reference Armitage2010). We can highlight the areas where planet formation via GI is possible in our modelling as the regions where
$Q_\textrm{Toomre}\leq1$
. They are shown in the upper panel of Fig. 9 for model M1. These regions appear between
$\approx$
10–100 au before
$\approx300$
kyr. However, at later times, clumps could appear in the disc as a result of an external perturbation, such as a stellar flyby (Thies et al. Reference Thies, Kroupa, Goodwin, Stamatellos and Whitworth2010). In this case, the planet would be formed from the material with altered C/O ratio, most probably with elevated amount of carbon, as the regions outside 5–10 au are typically more gravitationally unstable. This means that GI can produce planets with super-solar C/O ratios, if it is induced by external influence at later stages of disc evolution.
The disc regions where the conditions for GI and SI are fulfilled in model M1. The regions and times where there is no instability are shaded in white. In the upper panel, the colour indicates the minimum value of
$Q_\textrm{ Toomre}$
at a given radius, if
$Q_\textrm{Toomre}\leq1$
. In the lower panel, the colour indicates the fraction of mass at a given radius where SI can be triggered according to Li & Youdin (Reference Li and Youdin2021) criterion. Positions of the snowlines are shown for reference in dashed lines.
Core accretion is another most widely discussed scenario of giant planet formation. Accretion of gas should produce atmospheres with the C/O ratio close to the one in the gas phase of protoplanetary disc. However, dust grains are also accreted, so pebble and planetesimal accretion can enrich the atmosphere in volatile components (Mordasini et al. Reference Mordasini, van Boekel, Mollière, Henning and Benneke2016; Danti, Bitsch, & Mah Reference Danti, Bitsch and Mah2023). This makes atmospheric C/O ratio closer to the ice-phase C/O, but in case of gas giants, the amount of the solids needed to compensate the prevalence of carbon in the gas should be quite high, up to hundreds of Earth masses (GRAVITY Collaboration et al. 2020). The C/O ratios in exoplanetary atmospheres are often interpreted in terms of pebble accretion, so planets with stellar C/O ratios are assumed to form in the environments where solid phase C/O is unprocessed and thus close to the initial value. One of such locations is beyond CO snowline, where most of the carbon- and oxygen- bearing material is in the ice, e.g. for Jupiter (e.g. Öberg & Wordsworth Reference Öberg and Wordsworth2019; Ohno & Ueda Reference Ohno and Ueda2021). In our models, this is rather the vicinity of the CO
$_2$
snowlines, and the pebbles beyond the CO snowline are mostly covered with carbon-rich CO ice (Topchieva et al. Reference Topchieva, Molyarova, Akimkin, Maksimova and Vorobyov2024). Additionally, the CO snowline is typically very far from the star (
$ \gt 40$
au), so such scenarios must rely on planet migration to obtain their current location. Interpretations relying on chemical modelling extend this region to include the area beyond CO
$_2$
snowline due to additional chemical processing of CO in this region, e.g. HR 8799e (Mollière et al. Reference Mollière2020). This puts milder constraints on the original distance from the star where pebbles should be accreted and requires less migration. Modelling of planet formation and migration including pebble and gas accretion puts the formation location of planets with super-solar C/O ratios beyond water and CO
$_2$
snowlines (Bitsch, Schneider, & Kreidberg Reference Bitsch, Schneider and Kreidberg2022).
In our modelling, the C/O in the ice is close to initial value in the regions beyond CO snowline, excluding the area of CO accumulation. Between CO and CO
$_2$
snowlines, it is lower, as our model does not include chemical processes apart from adsorption and desorption. However, there is another region with C/O in the ice close to initial. The vicinity of primary water snowline and the ring induced outside of it has values of C/O in the ice only slightly above the initial value of 0.34. It is surrounded by the snowlines of CO
$_2$
. This region could be another favourable location for forming planets with the stellar C/O. As it is situated closer to the star, it would imply less migration.
Rare planets with lower than stellar C/O ratio, such as
$\beta$
Pic b (GRAVITY Collaboration et al. 2020; Reggiani et al. Reference Reggiani, Galarza, Schlaufman, Sing, Healy, McWilliam, Lothringer and Pueyo2024), HD 209458 b (Xue et al. Reference Xue, Bean, Zhang, Welbanks, Lunine and August2024), HD 189733b (Fu et al. Reference Fu2024), or KELT-1 b, Kepler-13A b andWASP-79 b (less precisely determined, see Hoch et al. Reference Hoch2023), need to have accreted a lot of oxygen-rich ice. Therefore, they are more likely to accrete solid material in the regions with the lowest ice-phase C/O ratios. The most suitable region would be at the distances between H
$_2$
O and CO
$_2$
snowlines, where ice mantles are made of pure water. However, in our modelling results, this region is very small, typically only a few au wide, as the snowlines are close to each other. This is because of steep temperature profile in this region, which is a result of the significant contribution of non-irradiation heating mechanisms, particularly viscous heating. Beyond CO
$_2$
snowline, there are also regions with relatively low (0.2–0.3) C/O in the ice, but much more solids need to be accreted in such areas to compensate for the excess of carbon from the gas.
Let us summarise the above constraints on planet formation locations and mechanisms implied by our simulated C/O ratios. Core accretion is suitable for forming planets with high C/O (
$\approx1$
) in the atmosphere around the snowlines of CO, CH
$_4$
and CO
$_2$
, or anywhere beyond CO
$_2$
snowline if they did not accrete much solids. Planets with stellar or slightly super-stellar C/O ratio need to accrete (a lot of) oxygen-rich solids to compensate their initially high C/O inherited from the gas. The locations where this is possible is between CO
$_2$
and CH
$_4$
and between CO and CH
$_4$
snowlines. Planets with low C/O ratio could accrete ices between H
$_2$
O and CO
$_2$
snowlines. Snowlines are favourable planetesimal formation sites, so the planets that accreted planetesimals/pebbles there can have altered C/O ratios. The values will be lower than the initial if they form at the water snowline, and higher if they form at the snowlines of carbon-rich species. At the same time, to obtain planets with stellar C/O formed at the snowlines, these planets would need to migrate and accrete matter in different regions of the disc to make their C/O ratio close to the initial stellar value. Alternatively, planets with stellar C/O ratio can form via disc fragmentation through GI at earlier stages. Dedicated modelling of planet formation accounting for evolution of dust and volatiles is necessary to put more particular constraints on planet formation scenarios.
Planetesimals play an important role in delivering the ice-phase elements to planetary atmospheres. To form planetesimals, additional physical process is needed, such as SE (see Youdin & Goodman Reference Youdin and Goodman2005), which is not explicitly included in our modelling because of insufficient numerical resolution and simplified vertical disc structure. However, we can post-process the simulation results to check if the conditions for SI are fulfilled in some regions of the disc where dust-to-gas ratio and dust size are enhanced, following Vorobyov et al. (Reference Vorobyov, Skliarevskii, Guedel and Molyarova2024). Dense dust rings forming at later stages (see Section 3.1) seem to be an ideal location for triggering the SI, which would ultimately lead to formation of planetesimals and then planets in the disc. Triggering the SI requires specific relations between local dust-to-gas ratio and Stokes number (Yang et al. Reference Yang, Johansen and Carrera2017). The criteria vary depending on the model, we adopt them from Li & Youdin (Reference Li and Youdin2021). Another criteria would be the requirement of volume density of dust to exceed that of gas in the midplane (Youdin & Goodman Reference Youdin and Goodman2005). We do not apply it, as in our modelling, the residual value of
$\alpha$
is
$10^{-3}$
which makes this condition unreachable outside of the dead zone. The regions in model M1 where the conditions of Li & Youdin (Reference Li and Youdin2021) are satisfied are shown in the lower panel of Fig. 9. Most of the suitable regions are in the inner disc (
$r \lt 20$
au) inside the dust rings, and appear after 200 kyr. However, there are some suitable regions between 10–100 au at earlier times, where SI could be triggered in the spirals.
The regions where GI and SI are possible shown in Fig. 9 are separated in space and time, and they have different characteristic C/O ratios. We can sum up all the volatiles in these regions (throughout the disc lifetime) to assess typical C/O ratios of the planet-forming material. For GI, we exclude the pre-disc phase (
$t \lt 53$
kyr) and consider the total C/O ratio, assuming both gas and solids are included in the forming planet. For SI, we separate the gas- and ice-phase C/O ratios. The formed planetesimals would only include the ices, however, if they form the planetary cores, these cores would also accrete gas. We note that the composition of the rocks, which are typically carbon-rich, is not included in our assessment. The resulting distributions of the C/O ratios in planet-forming regions are shown in Fig. 10. For GI regions, the C/O distribution has a distinct and relatively narrow peak around 0.5. It is slightly higher than the initial value of
$0.34$
. For SI regions, the ice-phase C/O is below 0.5, with major peaks at 0 and
$\approx0.2$
, and the gas-phase C/O has a broad distribution with multiple peaks between
$\approx$
0.2–1.4. Distributions of C/O ratios in the regions where GI and SI can be triggered are noticeably different.
Distribution of C/O ratios in the regions where gravitational and streaming instabilities are triggered. For GI, total C/O ratio is shown, for SI, the C/O ratios in the ice and in the gas. Black and grey points show the observed C/O ratios in two populations of exoplanets, the data is adopted from Hoch et al. (Reference Hoch2023).
It was shown by Hoch et al. (Reference Hoch2023) that there are two different populations of C/O ratios observed in giant exoplanets. They find that directly imaged exoplanets have
$\textrm{C/O}\approx$
0.5–0.8, while transiting hot Jupiters have a wider variety of C/O ratios (
$\approx$
0.3–1.7, see Figures 12 and 13 in Hoch et al. Reference Hoch2023), and suggest that these two populations could have different formation pathways. We add the C/O data of exoplanetary atmospheres compiled in Table 3 of Hoch et al. (Reference Hoch2023) to Fig. 10 (with arbitrary position at the y-axis). The narrow distribution of C/O ratios in the regions with GI is in step with the distribution of directly imaged exoplanets, albeit with a slightly shifted value due to our assumed initial conditions, while the wide range of C/O values in the regions of SI matches the variety of C/O ratios in transiting exoplanets. This could suggest that directly imaged exoplanets could form as a result of GI, which is also in line with their typically higher masses and orbital separations. At the same time, the transiting hot Jupiters could have experienced a lot of migration (Lin, Bodenheimer, & Richardson 1996; Dawson & Johnson Reference Dawson and Johnson2018, during which they accrete material with a variety of C/O ratios both from the gas and solid phase. It suggests that they could also form in the core accretion scenario. The origin of wide separation planets was also investigated by Bergin et al. (Reference Bergin, Booth, Colmenares and Ilee2024), based on the comparison with the observed
$\textrm{C/O} \gt 1$
in protoplanetary discs (including the full composition of solids). They conclude that both core accretion and GI can work as the formation mechanism of these planets.
Apart from (exo)planets, the C/O ratios can be measured for the comets, which present the best preserved sample of the primordial composition of the ices in the Solar System. Spectroscopic measurements of molecular composition in the comae suggest that the C/O ratio of cometary ice is quite low, typically below
$0.1$
due to the dominance of water ice (A’Hearn et al. 2012; Seligman et al. Reference Seligman2022; Harrington Pinto et al. Reference Harrington Pinto, Womack, Fernandez and Bauer2022). Although most comets are carbon-depleted, there are individual measurements of C/O in comets above 0.5, for example in C/2006 W3 Christensen and 29P/Schwassmann–Wachmann (Ootsubo et al. Reference Ootsubo2012; Seligman et al. Reference Seligman2022), or even close to 1 in C/2016 R2 (PanSTARRS) (Wierzchos & Womack Reference Wierzchos and Womack2018; McKay et al. Reference McKay2019). Additionally, high value of
$\textrm{C/O}\approx 1$
was observed in the interstellar object 2I/Borisov (Bodewits et al. Reference Bodewits2020). Our modelling results show the ice-phase
$\textrm{C/O}=0$
in the vicinity of the water snowline, as well as low values between the CO
$_2$
and CH
$_4$
snowlines (
$\approx0.2$
) and between the CH
$_4$
and CO snowlines (
$\approx0.3$
). These are the locations where comets could originate from. However, in the vicinity of the CO
$_2$
, CH
$_4$
and CO ice lines themselves, the C/O ratio in the ice phase is much higher. The fact that carbon-rich cometary ices are extremely rare in the Solar System may indicate that planetesimals formed on ice lines from carbon-rich volatiles do not persist throughout the evolution of a planetary system. This means that they are likely to be included in larger bodies, which favours the snowline-aided planet formation scenarios (Drążkowska & Alibert 2017; Hyodo et al. Reference Hyodo, Guillot, Ida, Okuzumi and Youdin2021). This is also consistent with the abundance of exoplanets with high C/O (Weiner Mansfield et al. Reference Weiner Mansfield2024), which could be formed around the snowlines of carbon-rich species. In the giant planets of the Solar System, the C/O ratios are not well constrained (Mousis et al. 2024b). However, the existing data suggest rather super-solar values for all giant planets except for Neptune (Cavalié et al. Reference Cavalié, Lunine, Mousis and Hueso2024); for Jupiter, the C/O ratio is assessed as
$\approx0.9$
(Wong et al. Reference Wong, Mahaffy, Atreya, Niemann and Owen2004; Li et al. Reference Li2024).
Our model only considers four most abundant chemical species. However, there can be other more complex molecules in protoplanetary discs, which could affect the balance of carbon and oxygen. The most obvious candidate is methanol CH
$_3$
OH, which has the abundance similar to methane in the protostellar cores (Öberg et al. 2011a). It was also observed in a protoplanetary disc around an erupting star V883 Ori (Lee et al. Reference Lee2019). We do not consider it in the model as its binding energy is close to that of water, thus the snowlines would have similar positions, but the abundance is an order of magnitude lower. However, it could somewhat increase the local C/O ratio in the inner regions where there are no other carbon-bearing species, such as the ice in the ring at 1 au. Including methanol would alter the distribution of the C/O ratio. Interactions between the ices considered in the model could also affect the results. As was recently shown by Ligterink, Kipfer, & Gavino (Reference Ligterink, Kipfer and Gavino2024), trapping of volatile species inside the mantles of less volatile ices could have a significant impact on the distribution of C/O ratios.
Another important process missing in our modelling is gas-phase and surface chemical reactions. They could significantly affect the distribution of C/O ratio in the gas and in the ice, particularly with high level of cosmic ray ionisation (Eistrup et al. Reference Eistrup, Walsh and van Dishoeck2016) or if carbon grain destruction is considered (Cridland et al. 2019a). One particular mechanism is the transformation of CO to CO
$_2$
on the surface of dust grains, which can lead to the depletion of CO from both gas and ice phases between CO and CO
$_2$
snowlines (Molyarova et al. Reference Molyarova, Akimkin, Semenov, Henning, Vasyunin and Wiebe2017; Bosman et al. Reference Bosman, Tielens and van Dishoeck2018). Considering this mechanism can change the conclusions about planet formation location (Mollière et al. Reference Mollière2020). Nevertheless, radial variations of the C/O ratio are necessary to explain molecular emission of discs with gaps (Leemker et al. Reference Leemker, Booth, van Dishoeck, Wölfer and Dent2024), and they can only be result of dust dynamics. In order to more consistently describe the distribution of molecules and elements in the disc, the models combining dust evolution and dynamics with more complex chemistry treatment are necessary.
Our simulations adopt the thin-disc approximation and focus on the midplane of the protoplanetary discs, in order to capture the essential physics of self-gravity, thermal balance, and dust evolution in a global modelling within reasonable computational times. This means that some relevant processes connected with the vertical structure are inevitably excluded. For example, vertical mixing and dust settling affect the C/O ratio in the upper layers of the disc (Krijt et al. Reference Krijt, Schwarz, Bergin and Ciesla2018, Reference Krijt, Bosman, Zhang, Schwarz, Ciesla and Bergin2020). Dust settling is implicitly included in our modelling through separate scale heights of drown dust and gas (as well as small dust), affecting dust number density in the midplane. However, this approach does not allow to reproduce vertical stratification in dust properties and chemical composition, which is particularly relevant for the interpretation of molecular observations. Vertical structure is also relevant for the accretion of matter on forming giant planets, which should proceed in 3D manner through meridional flows (Morbidelli et al. Reference Morbidelli, Szulágyi, Crida, Lega, Bitsch, Tanigawa and Kanagawa2014). Cridland et al. (2020a) showed that the C/O ratio in the atmospheres of giant planets is rather affected by the composition of the molecular layer than that of the midplane.
5. Conclusions
In this work, we studied the distribution of volatiles in a viscous self-gravitating protoplanetary disc with dust evolution using a thin-disc hydrodynamic code FEOSAD (Vorobyov et al. Reference Vorobyov, Akimkin, Stoyanovskaya, Pavlyuchenkov and Liu2018; Molyarova et al. Reference Molyarova, Vorobyov, Akimkin, Skliarevskii, Wiebe and Güdel2021). We calculated the C/O elemental ratio in the gas, in the ice, and in total, identified the key properties of the distribution of elements over 500 kyr of disc evolution and considered their implications for planet formation theory. Our main findings can be summarised as follows.
-
• The simulated C/O ratios in the regions where GI and SI conditions are fulfilled are consistent with the C/O ratios in two populations of exoplanets possibly formed in different mechanisms pointed out by Hoch et al. (Reference Hoch2023). We show that narrow C/O distribution of directly imaged planets is consistent with their formation via GI, while a variety of C/O in transiting hot Jupiters is in line with their migration through varying C/O conditions after the formation via either core accretion or GI.
-
• The lower C/O ratio in the ice between the CO
$_2$
, CH
$_4$
and CO snowlines is consistent with the typical composition of Solar System comets, while the higher value of
$\textrm{C/O}\approx 0.5-1$
at these snowlines corresponds to the composition of rare carbon-rich comets. This may indicate that matter from the snowlines is hardly preserved during the evolution of the disc and planetary system, possibly due to the inclusion in to planets. -
• The distribution of volatiles is affected by the disc substructures, such as rings and spirals, as well as by dust radial drift. Variations of physical conditions create multiple snowlines of CO
$_2$
and H
$_2$
O inside 10 au. Dust drift of icy grains brings the volatiles from the outer to the inner disc, enriching the inner disc with both C and O. It also creates a radial gradient of the total C/O ratio: its value is around 0.2 where water is not frozen, and 0.6–0.9 where it is icy, compared to the initial value of
$0.34$
. -
• Volatiles accumulate at their snowlines in both ice and gas phases due to the combined effect of dust drift and azimuthal variations of gas and dust radial velocities in a self-gravitating, non-axisymmetric disc. The species with low initial abundances, such as CH
$_4$
(or methanol not considered here), can significantly affect C/O ratio, as their accumulation at the snowline creates a bump in C/O in all phases above
$1.0$
. The total mass of the model affects the timescales and the magnitude of the accumulation by a factor of two. -
• Forming planets can accrete gas with
$\textrm{C/O} \gt 1$
beyond CO
$_2$
snowline, ices with
$\textrm{C/O}\approx0.5-1$
at the CO, CH
$_4$
and CO
$_2$
snowlines, ices with
$\textrm{C/O}\approx 0.2-0.3$
between these snowlines and ices with
$\textrm{C/O}=0$
between H
$_2$
O and CO
$_2$
snowlines. Planets with stellar C/O would need to migrate through these regions to acquire necessary composition or form via GI at earlier stages from the mixture of gas and dust with unaltered C/O ratio. -
• Dust-to-gas mass ratio and the total C/O ratio are systematically anticorrelated, because in dust-rich regions the volatile composition is close to that of the ice (which is lower), and in dust-poor regions, gas determines the C/O ratio.
The connection between protoplanetary disc components and exoplanets based on their composition should be more thoroughly investigated in the models focused on the planet formation process. We emphasise that these models should also take into account the effect of dust evolution and dynamics on the distribution of the elements in the planet-forming material. Inclusion of chemical processes and more accurate consideration of the bulk composition of dust grains could also affect the C/O ratios of the planet-forming environment.
Acknowledgement
We are thankful to the anonymous referee for useful comments that helped to improve the manuscript. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC) and the local computing facility of the Southern Federal University.
Data availability statement
The data underlying this article will be shared on reasonable request to the corresponding author.
Funding statement
The work is supported by Russian Science Foundation grant 22-72-10029, https://rscf.ru/project/22-72-10029/
Competing interests
None.
Open access
