© ESO 2018

1 Introduction

To date, a few thousand planetary systems have been found¹. Most of them have system architectures that are very different from our own solar system (Madhusudhan et al. 2014) and explaining the large variety of systems is a challenge for current planet formation theories (see, e.g. Morbidelli & Raymond 2016, and references therein). Thus, the birth environment of planets – protoplanetary disks – are an active area of study. A young stellar system inherits small dust grains from the interstellar medium. In regions with high densities and low turbulence, the grains start to coagulate. In the midplane of protoplanetary disks, where densities are higher than 10⁸ cm⁻³, grain growth can really take off. Grain growth and the interactions of these grown particles with the gaseous disk are of special interest to planet formation (see, e.g. Weidenschilling 1977; Lambrechts & Johansen 2012; Testi et al. 2014). The growth of dust grains to comets and planets is far from straightforward, however.

Pebbles,that is, particles that are large enough to slightly decouple from the gas, have been invoked to assist the formation of planets in different ways (Johansen & Lambrechts 2017). They are not supported by the pressure gradient from the gas, but they are subject to gas drag. As a result pebbles drift on a time-scale that is an order of magnitude faster than the gas depletion time-scale. This flow of pebbles, if intercepted or stopped, can help with planet formation. The accretion of pebbles onto forming giant planetary cores should help these cores grow beyond their classical isolation mass (Ormel & Klahr 2010; Lambrechts & Johansen 2012) while the interactions of gas and pebbles near the inner edge of the disk can help with the formation of ultra compact planetary systems as found by the Kepler satellite (Tan et al. 2016; Ormel et al. 2017). Efficient creation and redistribution of pebbles would lead to quick depletion of the solid content of disks increasing their gas-to-dust ratios by one to two orders of magnitude in 1 Myr (e.g. Ciesla & Cuzzi 2006; Brauer et al. 2008; Birnstiel et al. 2010).

Observations of disks do not show evidence of strong dust depletion. The three disks that have far-infrared measurements of the HD molecule to constrain the gas content show gas-to-dust ratios smaller than 200 (Bergin et al. 2013; McClure et al. 2016; Trapman et al. 2017), including the 10 Myr old disk TW Hya (Debes et al. 2013). Gas-to-dust ratios foundfrom gas mass estimates using CO line fluxes are inconsistent with strong dust depletion as well (e.g. Ansdell et al. 2016; Miotello et al. 2017). Manara et al. (2016) also argue from the observed relation between accretion rates and dust masses, that the gas-to-dust ratio in the 2–3 Myr old Lupus region should be close to 100.

It is thus of paramount importance that a way is found to quantify the inwards mass flux of solid material due to radial drift from observations. This is especially important for the inner ( < 10 AU) regions of protoplanetary disks. Here we propose to look for a signature of efficient radial drift in protoplanetary disks through molecules that are a major component of icy planetesimals.

Radial drift is expected to transport large amounts of ice over the iceline, depositing a certain species in the gas phase just inside and ice just outside the iceline (Stevenson & Lunine 1988; Piso et al. 2015; Öberg & Bergin 2016). This has been modelled in detail for the water iceline by Ciesla & Cuzzi (2006) and Schoonenberg & Ormel (2017), for the CO iceline by Stammler et al. (2017) and in general for the H₂ O, CO and CO₂ icelines by Booth et al. (2017). In all cases the gaseous abundance of a molecule in the ice is enhanced inside of the iceline as long as there is an influx of drifting particles. The absolute value and width of the enhancement depend on the mass influx of ice and the strength of viscous mixing. Such an enhancement may be seen directly with observations of molecular lines. Out of the three most abundant ice species (CO, H₂ O and CO₂ ), CO₂ is the most promising candidate for a study of this nature. Both CO and H₂ O are expected to be highly abundant in the inner disk based on chemical models (see, e.g. Aikawa & Herbst 1999; Markwick et al. 2002; Walsh et al. 2015; Eistrup et al. 2016). As such any effect of radial transport of icy material will be difficult to observe. CO₂ is expected to be abundant in outer disk ices (with abundances around 10⁻⁵ if inherited from the cloud or produced in situ in the disk, Pontoppidan et al. 2008; Boogert et al. 2015; Le Roy et al. 2015; Drozdovskaya et al. 2016), but it is far less abundant in the gas in the inner regions of the disk (with abundances around 10⁻⁸ , Pontoppidan et al. 2010; Bosman et al. 2017). This gives a leverage of three orders of magnitude to see effects from drifting icy pebbles.

Models by Cyr et al. (1998) and Ali-Dib et al. (2014), for example, predict depletion of volatiles in the inner disk. In these models, all volatiles are locked up outside of the iceline in stationary solids. In Ali-Dib et al. (2014) this effect is strengthened by the assumption that the gas and the small dust in the disk midplane are moving radially outwards such as proposed by Takeuchi & Lin (2002). Together this leads to very low inner disk H₂ O abundances in their models. CO₂ is not included in these models, but the process for H₂O would also work for CO₂, but slightly slower, as the CO₂ iceline is slightly further out. However, these models do not include the diffusion of small grains through the disk, which could resupply the inner disk with volatiles at a higher rate than that caused by the radial drift of large grains.

Bosman et al. (2017) showed that an enhancement of CO₂ near the iceline would be observable in the ¹³CO₂ mid infrared spectrum, if that material were mixed up to the upper layer. Here we continue on this line of research by constrainingthe maximal mass transport rate across the iceline and by investigating the shape and amplitude of a possible CO₂ abundance enhancement near the iceline due to this mass transport.

To constrain the mass transport we have build a model along the same lines as Ciesla & Cuzzi (2006) and Booth et al. (2017) except that we do not include planet formation processes. As in Booth et al. (2017) we use the dust evolution prescription from Birnstiel et al. (2012). The focus is on the chemistry within the CO₂ iceline to make predictions on the CO₂ content of the inner disk. In contrast with Schoonenberg & Ormel (2017) a global disk model is used to maintain consistency between the transported mass and observed outer disk masses. Chemical studies of the gas in the inner disk have been presented by, for example Agúndez et al. (2008), Eistrup et al. (2016), and Walsh et al. (2015) but transport processes are not included in these studies. Cridland et al. (2017) use an evolving disk, including grain growth and transport, coupled with a full chemical model in their planet formation models. However, they do not include transport of ice and gas species due to the various physical evolution processes.

Section 2 presents the details of our physical model, whereas Sect. 3 discusses various midplane chemical processes involving CO₂ and our method for simulating infrared spectra. Section 4 presents the model results for a range of model assumptions and parameters. One common outcome of all of these models is that the CO₂ abundance in the inner disk is very high, orders of magnitude more than observed, making it difficult to quantify mass transport. Section 5 discusses possible ways to mitigate this discrepancy and the implications for the physical and chemical structure of disks, suggesting JWST-MIRI observations of ¹³CO₂ that can test them.

2 Physical model

2.1 Gas dynamics

To investigate the effect of drifting pebbles on inner disk gas-phase abundances we build a 1D dynamic model. This model starts with an α-disk model (Shakura & Sunyaev 1973). The evolution of the surface density of gas Σ_gas (t, r) can be described as: $$\frac{\partial {\text{UpSigma}}_{\text{gas}}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\left(\frac{\alpha c_s^2}{\text{UpOmega}}{\text{UpSigma}}_{\text{gas}}{r}^{1/2}\right)\right],$$(1)

where r is the distance to the star, t is time, α is the dimensionless Sakura–Sunyaev parameter, Ω is the local Keplerian frequency and $c_s=\sqrt{k_{b}T/\mu }$ is the local sound speed, with k_b the Boltzmann constant, T the local temperatureand μ the mean molecular mass which is taken to be 2.6 amu. $\alpha c_s^2/\text{UpOmega}$ is also equal to ν_turb, the (turbulent) gas viscosity. The viscosity, or the resistance of the gas to shear, is responsible for the exchange of angular momentum. The origin of the viscosity is generally assumed to be turbulence stirred up by the magneto-rotational instability (MRI; Turner et al. 2014) although Eq. (1) is in principle agnostic to the origin of the viscosity as long as the correct α value is included. For the rest of the paper, we assume that the viscosity originates from turbulence and that the process responsible for the viscous evolution is also the dominant process in mixing constituents of the disk radially.

The evolution of the surface density of a trace quantity has two main components. First, all gaseous constituents are moving with the flow of the gas. This advection is governed by: $$\frac{\partial {\text{UpSigma}}_{x,\text{gas}}}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(F_{\text{gas}}r\frac{{\text{UpSigma}}_{x,\text{gas}}}{{\text{UpSigma}}_{\text{gas}}}\right),$$(2)

where F_gas is the total radial flux, which is related to the radial velocity of the gas due to viscous accretion, u_gas , and is given by: $$F_{\text{gas}}={\text{UpSigma}}_{\text{gas}}{u}_{\text{gas}}=\frac{3}{\sqrt{r}}\frac{\partial }{\partial r}\nu {\text{UpSigma}}_{\text{gas}}\sqrt{r}.$$(3)

Second, the gas is also mixed by the turbulence, smoothing out variations in abundance. This diffusion can be written as (Clarke & Pringle 1988; Desch et al. 2017) $$\frac{\partial {\text{UpSigma}}_{x,\text{gas}}}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_x{\text{UpSigma}}_{\text{gas}}\frac{\partial }{\partial r}\left(\frac{{\text{UpSigma}}_{x,\text{gas}}}{{\text{UpSigma}}_{\text{gas}}}\right)\right),$$(4)

where D_x is the gas mass diffusion coefficient. The diffusivity is related to the viscosity by: $$\text{Sc}=\frac{{\nu }_{\text{gas}}}{D_x},$$(5)

with Sc the Schmidtnumber. For turbulent diffusion in a viscous disk it is expected to be of order unity. As such, Sc = 1 is assumed for all gaseous components.

Advection and diffusion are both effective in changing the abundance of a trace species if an abundance gradient exists. Diffusion is most effective in the presence of strong abundance gradients and strong changes in the abundance gradient, such as near the iceline. At the icelines the diffusion of a trace species will generally dominate over the advection due to gas flow.

2.2 Dust growth and dynamics

The dynamics of dust grains are strongly dependent on the grain size. A grain with a large surface area relative to its mass is well coupled to the gas and will not act significantly different from a molecule in the gas. Solid bodies with a very small surface area relative to its mass, for example, planetesimals, will not be influenced by the gas pressure or turbulence, their motions are then completely governed by gravitational interactions. Dust particles with sizes between these extremes will be influenced by the presence of the gas in various ways. To quantify these regimes it is useful to consider a quantity known as the Stokes number: The Stokes number is, assuming Epstein drag and spherical particles in a vertically hydrostatic disk, given by (Weidenschilling 1977; Birnstiel et al. 2010): $$\text{St}=\frac{a_{\text{grain}}{\rho }_s}{{\text{UpSigma}}_{\text{gas}}}\frac{\pi }{2}.$$(6)

Particles with a very small Stokes number ( ≪1) are very well coupled to the gas and the gas pressure gradients and particles with very large Stokes number ( ≫1) are decoupled from the gas.

The coupling between gas and dust determines both the diffusivity of the dust, that is, how well the dust mixed due to turbulent motions of the disk as well as advection of the dust due to bulk flows of the gas. Youdin & Lithwick (2007) derived that the diffusivity D_dust of a particle with a certain Stokes number can be related to the gas diffusivity by: $$D_{\text{dust}}=\frac{D_{\text{gas}}}{1+{\text{St}}^2}.$$(7)

Similarly the advection speed of dust due to gas advection can be given by: $$u_{\text{dust},\text{adv}}=\frac{u_{\text{gas}}}{1+{\text{St}}^2}.$$(8)

When particles are not completely coupled to the gas they are also no longer completely supported by the gas pressure gradient. The gas pressure gradient, from high temperature, high density material close to the star, to low temperature, low density material far away from the star provides an outwards force, such that the velocity with which the gas has a stable orbit around the star is lower than the Keplerian velocity. Particles that start to decouple from the gas thus also need to speed up relative to the gas to stay in a circular orbit. This induces a velocity difference between the gas and the dust particles. The velocity difference between the gas velocity and a Keplerian orbital velocity is given by: $$\text{UpDelta}u=\text{UpOmega}r-\sqrt{{\text{UpOmega}}^2{r}^2-\frac{r}{{\rho }_{\text{gas}}}\frac{\partial P}{\partial r}},$$(9)

with P the pressure of the gas. We note that in the case of a positive pressure gradient, the gas will be moving faster than the Keplerian velocity

As a result of this velocity difference the particles are subjected to a drag force, which, in the case of a negative pressure gradient, removes angular momentum from the particles. The loss of angular momentum results in an inwards spiral of the dust particles. This process is known as radial drift. The maximal drift velocity can be given as (Weidenschilling 1977): $$u_{\text{drift}}=\frac{2\text{UpDelta}u}{\text{St}+{\text{St}}^{-1}}.$$(10)

The drift velocity will thus be maximal for particles with a Stokes number of unity. Drift velocities of ~1% of the orbital velocity are typical for particles with a Stokes number of unity.

The dynamics of dust are thus intrinsically linked to the size, or rather size distribution of the dust particles. The dust size distribution is set by the competition of coagulation and fragmentation². When two particles collide in the gas they can either coagulate, that is, the particles stick together and continue on as a single larger particle, or they can fragment, the particles break apart into many smaller particles. The outcome of the collision depends on the relative velocity of the particles and their composition. At low velocities particles are expected to stick, while at high velocities collisions lead to fragmentation. The velocity that sets the boundary between the two outcomes is called the fragmentation velocity. For pure silicate aggregates the fragmentation velocities are low (1 m s⁻¹ ) while aggregates with a coating of water ice can remain intact in collisions with relative velocities up to a order of magnitude higher (Blum & Wurm 2008; Gundlach & Blum 2015).

Dust size distributions resulting from the coagulation and fragmentation processes can generally not be computed analytically. Calculations have been done for both static and dynamic disks (Brauer et al. 2008; Birnstiel et al. 2010; Krijt & Ciesla 2016). These models are very numerically intensive and do not lend themselves to the inclusion of additional physics and chemistry nor to large parameter studies. As such we will use the simplified dust evolution prescription from Birnstiel et al. (2012). This prescription has been benchmarked against models with a more complete grain growth and dust dynamics prescription from Brauer et al. (2008). The prescription only tracks the ends of the dust size distribution, the smallest grains of set size and the largest grains at a location in the disk of a variable size. As this prescription is a key part of the model we will reiterate some of the key arguments, equations and assumption of this prescription here, for a more complete explanation, see Birnstiel et al. (2012).

The prescription assumes that the dust size distribution can be in one of three stages at any point in the disk. Either the largest particles are still growing, the largest particle size is limited by radial drift, or the largest particle size is limited by fragmentation. In the first and second case, the size distribution is strongly biased towards larger sizes. In the final case the size distribution is a bit flatter (see, Brauer et al. 2008). The size distribution in all cases is parametrised by a minimal grain size, the monomer size, a maximal grain size, which depends on the local conditions, and the fraction of mass in the large grains.

From physical considerations one can write a maximal expected size of the particles due to different processes. Growth by coagulation is limited by the amount of collisions and thus by the amount of time that has passed. The largest grain expected in a size distribution at a given time is given by: $$a_{\text{grow}}(t)=a_{\text{mono}}\exp \left[\frac{t-t_0}{{\tau }_{\text{grow}}}\right],$$(11)

where a_mono is the monomer size, t₀ is the starting time, t is the current time and τ_grow is the local growth time scale, given by: $${\tau }_{\text{grow}}=\frac{{\text{UpSigma}}_{\text{gas}}}{{\text{UpSigma}}_{\text{dust}}\text{UpOmega}}.$$(12)

Radial drift moves the particles inwards: the larger the particles, the faster the drift. There is thus a size at which particles are removed faster due to radial drift than they are replenished by grain growth, limiting the maximal size of particles. This maximal size is given by: $$a_{\text{drift}}=f_d\frac{8{\text{UpSigma}}_{\text{dust}}}{{\rho }_s}\frac{\pi r^2{\text{UpOmega}}^2}{c_s^2}{\gamma }^{-1},$$(13)

with f_d a numerical factor, ρ_s is the density of the grainsand γ is the absolute power law slope of the gas pressure: $$\gamma =\left|\frac{\text{d}\ln P_{\text{gas}}}{\text{d}\ln r}\right|.$$(14)

As mentioned before, particles that collide with high velocities are expected to fragment. For this model we consider two sources of relative velocities. One source of fragmentation is differential velocities due to turbulence (see, Ormel & Cuzzi 2007, for a more complete explanation). This limiting size is given by: $$a_{\text{frag}}=f_f\frac{2}{3\pi }\frac{{\text{UpSigma}}_{\text{gas}}}{{\rho }_s\alpha }\frac{u_f^2}{c_s^2},$$(15)

where f_f is a numerical fine tuning parameter and u_f is the fragmentation velocity. The other source of fragmentation considered is different velocities due to different radial drift speeds. The limiting size for this process is given by: $$a_{\text{df}}=\frac{u_{f}r\text{UpOmega}}{c_s^2(1-N)}\frac{4{\text{UpSigma}}_{\text{gas}}}{{\rho }_s}{\gamma }^{-1},$$(16)

N is the average Stokes number fraction between two colliding grains, here we use N = 0.5.

The size distribution at a location in the disk spans from the monomer size (a_mono) to the smallest of the four limiting sizes above. In the model the grains size distribution is split into “small” and “large” grains. The mass fraction of the large grains depends on the process limiting the size: if the grain size is limited by fragmentation ( min (a_frag, a_df) < a_drift), the fraction of mass in large grains (f_{m, frag}) is assumed to be 75%. If the grain size is limited by radial drift (a_drift < min(a_frag, a_df)), the fraction of mass in large grains (f_{m, drift}) is assumed to be 97%. These fractions were determined by Birnstiel et al. (2012) by matching the simplified model to more complete grain-growth models.

Using these considerations the advection-diffusion equation for the dynamics of the dust can be rewritten: $$\begin{array}{lll}\frac{\partial {\text{UpSigma}}_{\text{dust}}}{\partial t}\hfill & =\hfill & \frac{1}{r}\frac{\partial }{\partial r}[\text{vphantom}\frac{1}{2}\left(u_{\text{dust},\text{s}}r(1-f_m){\text{UpSigma}}_{\text{dust}}+u_{\text{dust},\text{l}}r{f}_m{\text{UpSigma}}_{\text{dust}}\right)\hfill \\ \hfill & \hfill & +r{\text{UpSigma}}_{\text{gas}}(D_{\text{dust},\text{s}}\frac{\partial }{\partial r}\left(\frac{(1-f_m){\text{UpSigma}}_{\text{dust}}}{{\text{UpSigma}}_{\text{gas}}}\right)\hfill \\ \hfill & \hfill & +D_{\text{dust},\text{l}}\frac{\partial }{\partial r}\left(\frac{f_m{\text{UpSigma}}_{\text{dust}}}{{\text{UpSigma}}_{\text{gas}}}\right))],\hfill \end{array}$$(17)

where $D_{\text{dust},\text{x}}=D_{\text{gas}}/(1+{\text{St}}_{\text{x}}^2)$, $u_{\text{dust},\text{x}}=u_{\text{gas}}/(1+{\text{St}}_{\text{x}}^2)+u_{\text{drift},\text{x}}$ where u_drift is the velocity due to radial drift (Eq. (10)). Here it is assumed that the large grains never get a Stokes number larger than unity, which holds for the models presented here.

The final part of the dynamics concerns the ice on the dust. The ice moves with the dust, so both the large scale movements as well as the mixing diffusion of dust must be taken into account. For the model we assume that the ice is distributed over the grains according to the mass of the grains. This means that if large grains have a mass fraction f_m , then the large grains will have the same mass fraction of f_m of the available ice. This is a reasonable assumption, as long as the grain size distribution is in fragmentation equilibrium and the largest grains are transported on a timescale longer than the local growth timescale. The algorithm used here forces the latter condition to be true everywhere in the disk, while the former condition is generally true for the CO₂ and H₂ O icelines (Brauer et al. 2008), but not for the CO iceline (Stammler et al. 2017). With this assumption the advection-diffusion equation forthe ice can be written as: $$\begin{array}{lll}\frac{\partial {\text{UpSigma}}_{\text{ice},x}}{\partial t}\hfill & =\hfill & \frac{1}{r}\frac{\partial }{\partial r}r(F_{\text{dust}}\frac{{\text{UpSigma}}_{\text{ice},x}}{{\text{UpSigma}}_{\text{dust}}}\hfill \\ \hfill & \hfill & +{\text{UpSigma}}_{\text{dust}}\left(D_{\text{dust},\text{s}}(1-f_m)+D_{\text{dust},\text{l}}{f}_m\right)\frac{\partial }{\partial r}\frac{{\text{UpSigma}}_{\text{ice},\text{x}}}{{\text{UpSigma}}_{\text{dust}}}).\hfill \end{array}$$(18)

Here F_dust is the radial flux of dust, this is given by: $$\begin{array}{lll}{F}_{\text{dust}}\hfill & =\hfill & {\text{UpSigma}}_{\text{dust}}\left(u_{\text{dust},s}(1-f_m)+u_{\text{dust},\text{l}}{f}_m\right)\hfill \\ \hfill & \hfill & +{\text{UpSigma}}_{\text{dust}}(D_{\text{dust},s}\frac{\partial }{\partial r}\left(\frac{(1-f_m){\text{UpSigma}}_{\text{dust}}}{{\text{UpSigma}}_{\text{gas}}}\right)\hfill \\ \hfill & \hfill & +D_{\text{dust},\text{l}}\frac{\partial }{\partial r}\left(\frac{f_m{\text{UpSigma}}_{\text{dust}}}{{\text{UpSigma}}_{\text{gas}}}\right)).\hfill \end{array}$$(19)

2.3 Model parameters

For our study, we pick a disk structure that should be representative of a young disk around a one solar mass star. The initial surface density structure is given by: ${\zeta }_{{\text{H}}_2}$(20)

where p controls the steepness of the surface density profile, r_taper the extent ofthe disk and Σ_{1 AU} sets the normalisation of the surface density profile. For our models we use p = 1 and r_taper = 40 AU. The temperature profile is given by a simple power law: ${\zeta }_{{\text{H}}_2}$(21)

The disk is assumed to be viscous with a constant α, as such there is a constraint on the power law slopes p and q, if we want the gas surface density to be a self similar solution to Eq. (1), it is required to have p + q = 3∕2 (Hartmann et al. 1998).

To calculate the volume densities, a vertical Gaussian distribution of gas with a scale height H_p (r) = h_pr, with h_p = 0.05, is used.

The disk is populated with particles of 0.1 μm at the start of the model, this is also the size of the small dust in the model. The gas-to-dust ratio is taken to be 100 over the entire disk. The density of the grains is assumed to be 2.5 g cm⁻³ and grains are assumed to be solid spheres.

The molecules are initially distributed through simple step functions. These step functions have a characteristic radius “iceline” which differentiates between the inner disk and the outer disk. Within this radius the gas phase abundance of the molecule is initialised, outside of this radius the solid phase abundance is initialised. Water is distributed equally through the disk with an abundance of 1.2 × 10⁻⁴, whereas CO₂ is initially depleted in the inner disk with an abundance of 10⁻⁸ (Pontoppidan et al. 2011; Bosman et al. 2017) and has a high abundance in the outer disk of 4 × 10⁻⁵. This ice abundance of CO₂ is motivatedby the ISM ice abundance (Pontoppidan et al. 2008; Boogert et al. 2015), cometary abundances (Le Roy et al. 2015)and chemical models of disk formation (Drozdovskaya et al. 2016).

A summary of the initial conditions, fixed and variable parameters is given in Table 1.

2.4 Boundary conditions

Due to finite computational capabilities, the calculation domain of the model needs to be limited. The assumptionsmade for the edges can have significant influences on the model results. For the inner edge, it is assumed that gas leaves the disk with an accretion rate as assumed from a self similar solution (p + q = 3∕2) according tothe viscosity and surface density at the inner edge. The accretion rate is given by: $$\dot{M}=3\pi {\text{UpSigma}}_{\text{gas}}{\nu }_{\text{gas}}.$$(22)

All gas constituents advect with the gas over the inner boundary or, in the case of large grains, drift over the boundary. Diffusion over the inner boundary is not possible. For the outer boundary, the assumption is made that the surface density of gas outside the domain is zero. Again viscous evolution or an advection process can remove grains and other gas constituents from the computational domain.

3 Chemical processes

3.1 Freeze-out and sublimation

Freeze-out and sublimation determine the fraction of a molecule that is in the gas phase and the fraction that is locked up on the grains. Freeze-out of a molecule, or a molecule’s accretion onto a grain is given by the collision rate of a molecule with the grain surface times the sticking fraction, f_s, which is assumed to be unity: $$R_{\text{acc},\text{x}}=f_s{\sigma }_{\text{dust}}{n}_{\text{grain}}{v}_{\text{therm},\text{x}},$$(23)

where σ_dust is the average dust surface area, n_grain is the number density of grains and the thermal velocity $v_{\text{therm},\text{x}}=\sqrt{8kT/\pi m_x}$, with k the Boltzmann constant, T the temperature and m_x the mass of the molecule. Molecules that are frozen-out on the grain can sublimate or desorb. For a grain covered with many monolayers of ice, the rate per unit volume for this process is given by: $$f_{\text{des},\text{x}}=R_{\text{des},\text{x}}{n}_{\text{x},\text{ice}}=p_x{\sigma }_{\text{dust}}{n}_{\text{grain}}{N}_{\text{act}}\exp \left[-\frac{E_{\text{bind}}}{kT}\right],$$(24)

where p_x is the zeroth-order “prefactor” encoding the frequency of desorption attempts per unit surface area, σ_dust is the surface area per grain, n_grain is the number density of grains. N_act is the number of active layers, that is, the number of layers that can participate in the sublimation, N_act = 2 is used, T is the dust temperature, which we take equal to the gas temperature. For mixed ices this rate can be modified by a covering fraction χ_x = n_{ice, x}∕∑ _xn_{ice, x}, however, this is neglected here. We note that f_des,x in its current form is independent of the amount of molecules frozen out on the dust grains. Using these rates we get the following differential equation: $$\frac{\partial n_{\text{x},\text{gas}}}{\partial t}=-R_{\text{acc},\text{x}}{n}_{\text{x},\text{gas}}+f_{\text{des},\text{x}},$$(25)

which has an analytical solution: $$\begin{array}{lll}{n}_{\text{x},\text{gas}}(t)\hfill & =\hfill & \min [\text{vphantom}\frac{1}{2}{n}_{\text{x},\text{tot}},\left(n_{\text{x},\text{gas}}(t_0)-\frac{f_{\text{des},\text{x}}}{R_{\text{acc},\text{x}}}\right)\hfill \\ \hfill & \hfill & \times \exp \left(-R_{\text{acc},\text{x}}(t-t_0)\right)+\frac{f_{\text{des},\text{x}}}{R_{\text{acc},\text{x}}}].\hfill \end{array}$$(26)

where n_x,tot is the total number density of a molecule (gas and ice). The number density of ice is given by: $$n_{\text{x},\text{ice}}(t)=n_{\text{x},\text{tot}}-n_{\text{x},\text{gas}}(t).$$(27)

The ice line temperature, defined as the temperature for which n_x,gas = n_x,dust, that is, when freeze-out and sublimation balance, depends on the total number density of the molecule considered. At lower total molecule number densities, the iceline will be at lower temperatures.

For CO₂ we use a binding energy of 2900 K, representative for CO₂ mixed with water (Sandford & Allamandola 1990; Collings et al. 2004). More recent measurements have suggested that the binding energy is lower, around 2300 K (Noble et al. 2012). Using the lower binding energy moves the CO₂ iceline further out from 90 to 80 K or from 6 to 8 AU in our standard model. The prefactor, $p_{x,{\text{CO}}_2}$, of 9.3 × 10²⁶ cm⁻² s⁻¹ from (Noble et al. 2012) is used. The change in iceline location has a minimal effect on the evolution of the abundance profiles. The mixing time-scale becomes longer at larger radii which would increase mixing times. For H₂ O a binding energy of 5600 K and prefactor, $p_{x,{\text{H}}_2\text{O}}$, of 10³⁰ cm⁻² s⁻¹ (Fraser et al. 2001).

3.2 Midplane formation and destruction processes

Radial drift and radial diffusion and advection will quickly move part of the outer disk CO₂ ice reservoir into the inner disk, enhancing the inner disk abundances. To get a good measure of the amount of CO₂ in the innerdisk it is necessary to also take into account the processes that form and destroy CO₂ in gas and ice. The density is highest near the mid-plane, as such this is where formation and destruction process are expected to be most relevant for the bulk of the CO₂. However, in the less dense upper layers, there are UV-photons that can dissociate and ionise molecules, possibly influencing the overall abundance of CO₂.

3.2.1 Gas-phase formation of CO₂

The formation of CO₂ in the inner disk mainly goes through the warm gas-phase route. Here CO₂ forms through the reaction CO + OH →CO₂ + H. The reaction has a slight activation barrier of 176 K (Smith et al. 2004). The parent molecule CO is very stable and is expected to be present at high abundances in the inner disk (10⁻⁴) (Walsh et al. 2015). The OH radical is expected to be less abundant, and it is the fate of this radical that determines the total production rate of CO₂ . OH is formed either directly from H₂O photodissociation (Heays et al. 2017), or by hydrogenation of atomic oxygen, O +H₂ →OH + H, in a reaction that has an activation barrier of 3150 K (Baulch et al. 1992). The atomic oxygen itself also has to be liberated from, in this case, either CO, CO₂ or H₂ O by X-rays or UV-photons. The production rate of CO₂ is thus severely limited if there is no strong radiation field present to release oxygen from one of the major carriers.

The CO₂ formation reaction, CO + OH →CO₂ + H, has competition from the H₂O formation reaction, OH + H₂ →H₂O + H. The hydrogenation of OH has a higher activation energy, 1740 K, than the formation of CO₂ , but since H₂ is orders of magnitude more abundant than CO, the formation of water will dominate over the formation of CO₂ at high temperatures. The rate for CO₂ formation is given by (Smith et al. 2004): ${\zeta }_{{\text{H}}_2}$(28)

The formation rate for H₂O formation is given by (Baulch et al. 1992): ${\zeta }_{{\text{H}}_2}$(29)

This means that the expected $x_{{\text{CO}}_2}/x_{{\text{H}}_2\text{O}}$ fraction from formation is: ${\zeta }_{{\text{H}}_2}$(30)

This function is plotted in Fig. 1 which shows that the formation of CO₂ is faster below temperatures of 150 K, whereas above this temperature formation of water is faster. Above a temperature of 300 K water formation is a thousand times faster than CO₂ formation. The implication is that gaseous CO₂ formation is only effective in a narrow temperature range, 50–150 K, and then only if OH is present as well, requiring UV photons or X-rays to liberate O and OH from CO or H₂O.

Fig. 1 Ratio of the CO2 to H2 O formation rate from a reaction of OH with CO and H2 respectively. Hydrogenation of OH to H2O dominates above 150 K. Formation times of CO2 and H2 O from OH are faster than the inner disk mixing time for all temperatures considered here.

3.2.2 Destruction of CO₂: cosmic rays

Cosmic rays, 10–100 MeV protons and ions, have enough energy to penetrate deeply into the disk. Cosmic-rays can ionise H₂ in regions where UV photons and X-rays cannot penetrate. The primary ionisation as well as the collisions of the resulting energetic electron with further H₂ creates electronically excited H₂ as well as excited H atoms. These excited atoms and molecules radiatively decay, resulting in the emission of UV-photons (Prasad & Tarafdar 1983). These locally generated UV-photons can dissociate and ionise molecules. Here only the dissociation rate of CO₂ is taken into account as ionisationof CO₂ by this process is negligible. Following Heays et al. (2017) the destruction rate of species X is written as: $$k_X=\frac{{\zeta }_{{\text{H}}_2}{x}_{{\text{H}}_2}}{x_X}{\displaystyle \int P}(\lambda )p_X(\lambda )\text{d}\lambda ,$$(31)

where ${\zeta }_{{\text{H}}_2}$ is the direct cosmic ray ionisation rate of H₂, x_X is the abundance of species X w.r.t. H₂ , P(λ) is the photon emission probability per unit spectral density for which we use the spectrum from Gredel et al. (1987). p_X (λ) is the absorption probability of species X for a photon of wavelength λ. This probability is given by: $$p_X(\lambda )=\frac{x_X{\sigma }_X^{\text{destr}}(\lambda )}{x_{\text{dust}}{\sigma }_{\text{dust}}^{\text{abs}}(\lambda )+{\displaystyle \sum _j{x}_j}{\sigma }_j^{\text{abs}}(\lambda )},$$(32)

where ${\sigma }_X^{\text{i}}(\lambda )$ is the wavelength dependent destruction or absorption cross section of species X and $x_{\text{dust}}{\sigma }_{\text{dust}}^{\text{abs}}$ is the dust cross section per hydrogen molecule. For the calculation of the cosmic-ray dissociation rate of CO₂ we assume that the destruction cross section in the UV is equal to the absorption cross section for CO₂ , that is, every absorption of a photon with a wavelength shorter than 227 nm leads to the destruction of a CO₂ molecule. For the calculations the cross sections from Heays et al. (2017) are used³ . These cross sections can also be used to compute destruction rates for CO₂ by stellar UV radiation in the surface layers of the disk.

The dust absorption is an important factor in these calculations and can be the dominant contribution to the total absorption in parts of the spectrum. The dust absorption greatly depends on the dust opacities assumed. A standard ISM dust composition was taken following Weingartner & Draine (2001), the mass extinction coefficients are calculated using Mie theory with the MIEX code (Wolf & Voshchinnikov 2004) and optical constants by Draine (2003) for graphite and Weingartner & Draine (2001) for silicates. Grain sizes are distributed assuming an MRN distribution starting at 5 nm, with varying maximum size are used. The resulting mass opacities and cross sections are shown in Fig. A.1.

Cosmic ray induced destruction rates for CO₂ are calculated for each dust size distribution for a range of CO₂ abundances between 10⁻⁸ and 10⁻⁴. The o-p ratio of H₂, important for the H₂ emission spectrum, is assumed to be 3:1, representative for high temperature gas. The abundances of the other shielding species used in the calculation are shown in Table 2. The destruction rate for CO₂ is plotted in Fig. 2.

The CO₂ destruction is faster for grown grains and low abundances of CO₂. At abundances above 10⁻⁷, the strongest transitions start to saturate, lowering the destruction rate with increasing abundance. Even though the dust opacity changes by more than an order of magnitude, the dissociation rates stay within a factor of 3 for all CO₂ abundances. This is due to the inclusion of H₂O into our calculations. H₂ O has large absorption cross sections in the same wavelength regimes as CO₂. When H₂ O is depleted, such as would be expected between the H₂O and CO₂ icelines, CO₂ destruction rates increase, especially for the largest grains. Even in this optimal case, the destruction rate for CO₂ is limited to 2 × 10⁻¹³ s⁻¹ for typical ${\zeta }_{{\text{H}}_2}$ of 10⁻¹⁷ s⁻¹.

Aside from generating a UV field, cosmic-rays also create ions, the most important of these being He⁺ . Due to the large ionisation potential of He, electron transfer reactions with He⁺ generally lead to dissociation of the newly created ion. For example, CO₂ + He⁺ preferably leads to the creation of O + CO⁺. Destruction of CO₂ due to He⁺ is limited by the creation of He⁺ and the competition between CO₂ and other reaction partners of He⁺. The total reaction rate coefficient of CO₂ with He⁺ is $k_{\text{ion},{\text{CO}}_2}=1.14\times {10}^{-9}$ cm³ s⁻¹ (Adams & Smith 1976). The three main competitors for reactions with He⁺ are H₂ O, CO and N₂ , with reaction rate coefficients of $k_{\text{ion},{\text{H}}_2\text{O}}=5.0\times {10}^{-10}$, k_{ion , CO} = 1.6 × 10⁻⁹ and $k_{\text{ion},{\text{N}}_2}=1.6\times {10}^{-9}$ cm³ s⁻¹ respectively (Mauclaire et al. 1978; Anicich et al. 1977; Adams & Smith 1976). Assuming a cosmic He⁺ ionisation rate of $0.65{\zeta }_{{\text{H}}_2}$ (Umebayashi & Nakano 2009), the destruction rate for CO₂ due to He⁺ reaction can be written as, $$R_{\text{destr},{\text{He}}^{+}}=0.65{\zeta }_{{\text{H}}_2}\frac{x_{\text{ce}He}}{x_{{\text{CO}}_2}}\frac{k_{\text{ion},{\text{CO}}_2}{x}_{{\text{CO}}_2}}{{\displaystyle \sum _X{k}_{\text{ion},\text{X}}}{x}_X},$$(33)

where the sum is over all reactive collision partners of He⁺ of which CO, N₂ and H₂ O are most important. For typical abundances of He of 0.1 and N₂, CO and H₂O of 10⁻⁴, $R_{\text{destr},{\text{He}}^{+}}$ is below $650{\zeta }_{{\text{H}}_2}$ for all CO₂ abundances.Thus cosmic-ray induced photodissociation will always be more effective than destruction due to He⁺ .

Altogether the destruction time-scale for midplane CO₂ by cosmic ray induced processes is long, ~3 Myr for a ${\zeta }_{{\text{H}}_2}={10}^{-17}$ s⁻¹. The latter value is likely an upper limit for the inner disk given the possibility of attenuation and exclusion of cosmic rays (Umebayashi & Nakano 1981; Cleeves et al. 2015).

Fig. 2 Dissociation rate of CO2 due to cosmic ray induced photons, for different dust size distributions. Solid lines are for condition inside of the H2 O iceline, dashed lines for conditions outside of the H2O iceline. The efficiency multiplied by ${\zeta }_{{\text{H}}_2}$ gives the CO2 dissociation rate.

Fig. 3 Time evolution series for a model without grain growth. This model assumes an α of 10−3 . Left: surface density of the gas (solid lines) and solids (dashed lines), the solid surface density is the sum of the dust and ice surface densities, it has been multiplied by a factor of 10 for visualization. Middle: abundance of CO2 . Right: abundance of H2O.

3.2.3 Destruction of CO₂: gas-phase reactions

In warm gas it is possible to destroy CO₂ by endothermic reactions with H or H₂ (Talbi & Herbst 2002). The reaction $${\text{CO}}_2+\text{H}\to \text{CO}+\text{OH}$$(34)

has an activation barrier of 13 300 K and a pre-exponential factor of 2.5 × 10⁻¹⁰ cm³ s⁻¹ for temperatures between 300 and 2500 K (Tsang & Hampson 1986) while the reaction $${\text{CO}}_2+{\text{H}}_2\to \text{CO}+{\text{H}}_2\text{O}$$(35)

has an activation barrier of 56 900 K and has a pre-exponential factor of 3.3 × 10⁻¹⁰ cm³ s⁻¹ at 1000 K. This means that for gas at 300 K, an H₂ density of 10¹² cm⁻³, and a corresponding H density of 1 cm⁻³, the rate for destruction by atomic hydrogen is 1.4 × 10⁻²⁹ s⁻¹, while the destruction rate by molecular hydrogen is 1.4 × 10⁻⁸⁰ s⁻¹. Both are far too low to be significant in the inner disk. However, if the atomic H abundance is higher, destruction of CO₂ by H can become efficient at high temperatures. As such the destruction of CO₂ becomes very dependent on the formation speed of H₂ at high densities and temperatures. The CO₂ abundance as a function of temperature for a gas-phase model is shown in Fig. B.1. Temperatures of >700 K are needed to lower the CO₂ abundance below 10⁻⁷, even with a high atomic H abundance.

3.3 Simulating spectra

To compare the model abundances with observations, infrared spectra are simulated with the thermochemical code DALI (Bruderer et al. 2012; Bruderer 2013). DALI is used to calculate a radial temperature profile, from dust and gas surface density profiles and stellar parameters. The viscous evolution model is initialised using the same surface density distribution as the DALI model. The temperature slope q is taken so that the surface density slope p is consistent with a self-similar solution q = 1.5 − p = 0.9. The temperature at 1 AU is taken from the midplane temperatures as calculated by the DALI continuum ray-tracing module. No explicit chemistry is included in this version of DALI, but the abundances are parametrised using the output of the dynamical model. The viscous model run with α = 10⁻³, u_f = 3 m s⁻¹ and a CO₂ destruction rate varying between 10⁻¹³ and 10⁻⁹ s⁻¹. The results from Sect. 4.4 show that these parameters represent the bulk of the model results.

After simulating 1 Myr of evolution, the gas-phase CO₂ abundance profile is extracted as function of temperature and interpolated onto the DALI midplane temperatures. The abundance is taken to be vertically constant up to the point where A_V = 1. In some variations an abundance floor of 10⁻⁹ or 10⁻⁸ is used for cells with A_V > 1 to simulate local CO₂ production. Using the resulting abundance structure, the non-LTE excitation of CO₂ is calculated using the rate coefficients from Bosman et al. (2017). Finally the DALI line ray-tracing module calculates the line fluxes for CO₂ and its ¹³C isotope.

For calculating the spectra, the same disk model is used as in Bruderer et al. (2015) and Bosman et al. (2017). The model is based on the disk AS 205 N. The parameters for the DALI models can be found in Table 3. For more specifics on the modelling of the spectra and adopted parameters, see Bruderer et al. (2015) and Bosman et al. (2017). To simulate the high gas-to-dust ratios that are inferred from water observations (Meijerink et al. 2009), the overall gas-to-dust ratio is set at 1000 throughout the disk. The artificially high gas-to-dust ratio does not affect any of the modelling, except for the line formation, which is only sensitive to the upper optically thin layers of the disk. High gas-to-dust ratios effectively mimic settled dust near the midplane containing ~90% of the dust mass.

4 Results

4.1 Pure viscous evolution

To start, a dynamical model without grain growth and without any chemistry (except for freeze-out and desorption) is investigated. Figure 3 shows the time evolution for a model with α =10⁻³. As expected, the total gas and dust surface densities barely change over 10⁶ yr. The gas anddust evolve viscously, some mass is accreted onto the star while the outer disk spreads a little bit. There are no changes in gas-to-dust ratio in the disk.

The water abundance in the right panel of Fig. 3 shows no evolution at all. This is to be expected because the only radial evolution is due to the gas viscosity, which affects gas and dust equally. Diffusion of the icy dust grains is equally quick as the diffusion of the water vapour. This means that all the water vapour that diffuses outwards is compensated by icy dust grains diffusing inwards.

The abundance of CO₂ in the middle panel of Fig. 3 does show evolution. Initially there is only a little bit of CO₂ gas in the inner disk, but a large amount of ice in the outer disk. The abundance gradient together with the viscous accretion makes the icy CO₂ move inwards, filling up the inner disk with CO₂ gas. This continues for about 1 × 10⁶ yr until the gaseous abundance of CO₂ is equal to the abundance of icy CO₂. This time-scale directly scales with the assumed α parameter. For α = 10⁻⁴ it takes more than 3 Myr to get a flat abundance profile in the inner disk.

4.2 Viscous evolution and grain growth

For the models including grain growth, we tested a range of fragmentation velocities from 1 to 30 m s⁻¹ . For some ofthese models, especially those with low fragmentation velocities and high α, the results are indistinguishable from the case without grain growth. Figures 4 and 5 show two models where the effect of grain growth and resulting drift can be seen.

Figure 4 shows the surface density and abundance evolution for a model with a fragmentation velocity of 3 m s⁻¹ as appropriate for pure silicate grains. The surface density of the dust shows a small evolution due to radial drift. Instead of decreasing, the surface density of dust in the inner 4 AU slightly increases in the first 300 000 yr due to the supply of dust particles from outside this radius.

The abundance profiles in Fig. 4 show distinct effects of radial drift. Both the gaseous H₂ O and CO₂ abundances are high at all times due to the influx of drifting icy grains. There is also a decrease in the abundance of ices at large radii. This is because the grains carrying the ices have moved inwards, increasing the gas-to-ice ratios.

The models are not in steady state after 1 Myr. If these models are evolved further, at first the abundances of H₂ O and CO₂ will increasefurther. At some point in time the influx of dust will slow down, because a significant part of the dust in the outer disk will have drifted across the snowline. At this point the inner disk abundances of H₂ O and CO₂ start to decrease as these molecules are lost due to accretion onto the star but no longer replenished by dust from the outer disk. For the model shown in Fig. 4 the average gas-to-dust ratio in the disk would be around 500 after 3 Myr.

For higher α the effects of radial drift are limited as the high rate of mixing smooths out concentration gradients created by radial drift. At the same time the maximum grain size is limited due to larger velocity collisions at higher turbulence. As such, for α = 0.01 a fragmentation velocity higher than 10 m s⁻¹ is needed to see significant effects of radial drift and grain growth.

For lower α the effects of grain growth get more pronounced, 90% of the dust mass is accreted onto the star in 2.5 Myr for a fragmentation velocity of 1 m s⁻¹. In this time a lot of molecular material is released into the inner disk and peak abundances of 10⁻³ and 10⁻² are reached for CO₂ and H₂ O, respectively, after 1 Myr of evolution. Due to the low turbulence in the gas, the volatiles released are not well mixed, neither inwards nor outwards, so there is no strong enhancement of the ice surface density just outside of the iceline. An overview of the different model evolutions can be found in Appendix C.

Increasing the fragmentation velocity allows grains in the model to grow to larger sizes, leading again to larger amounts of radial drift. Figure 5 shows the evolution of the same model as shown in Fig. 4, but with a fragmentation velocity of 30 m s⁻¹ , representative for grains coated in water ice. The solid surface density distribution clearly shows the effects of grain growth. At all radii, mass is moved inwards at a high rate. When the silicate surface density has dropped by an order of magnitude, a local enhancement in the surface density of the total solids is seen at the icelines. This enhancement is about a factor of 2 for both icelines.

The gas-phase abundances of CO₂ and H₂ O both show a large increase in the first 3 × 10⁵ yr due to the large influx of icy pebbles. At later times, the abundances are decreased again as the volatiles are accreted onto the star. After 3 × 10⁶ yr the H₂ O and CO₂ abundances are lower than 10⁻⁵.

Fig. 4 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 3 m s−1. Panels as in Fig. 3.

Fig. 5 Time evolution series for a model with grain growth. The fragmentation velocity for this model is 30 m s−1 and α = 10−3. The verticalscales are different from Figs. 3 and 4. Panels as in Fig. 3.

4.3 Viscous evolution and CO₂ destruction

The models shown in Figs. 3–5 without any chemical processing predict CO₂ abundances that are very high in the inner regions, 10⁻⁵–10⁻⁴, orders of magnitude higher that the value of ~10⁻⁸ inferred from observations (Pontoppidan & Blevins 2014; Bosman et al. 2017). There are multiple explanations for this disparity, both physical and chemical, which are discussed in Sect. 5. Here we investigate how large any missing chemical destruction route for gaseous CO₂ would need to be. As discussed in Sect. 3 and Appendix sectionlinking B.1 , midplane CO₂ can only be destroyed by cosmic-ray induced processes in the current networks. However, both cosmic ray induced photodissociation and He⁺ production are an order of magnitude slower than the viscous mixing time. We therefore introduce additional destruction of gaseous CO₂ with a rate that is a constant over the entire disk. This rate is varied between 10⁻¹³ and 10⁻⁹ s⁻¹ to obtain agreement with observations. For comparison, the cosmic ray induced process generally have a rate of the order of 10⁻¹⁴ s⁻¹ (see Sect. 3.2.2). The rate is implemented as an effective destruction rate, so there is no route back to CO₂ after destruction. To get the absolute rate one has to also take into account the reformation efficiency, but that can be strongly dependent on the destruction pathway, of which we are agnostic. These efficiencies are discussed in Sect. 5 and Appendix B.

Figure 6 shows the abundance evolution for a model where only gaseous CO₂ is destroyed, at a rate of 10⁻¹¹ s⁻¹. In this case the innermost parts of the disk are empty of CO₂, whereas the abundance of CO₂ reaches values close to the initial ice abundances around the iceline. Due to the constant destruction of CO₂ near the iceline, the actual ice abundance of CO₂ is also lower than the initial value.

The left panel of Fig. 7 shows the CO₂ abundance distribution for models with α = 10⁻³. The abundance profiles after 1 Myr of evolution are presented, when a semi-steady state has been reached. The peak abundance and peak width of the CO₂ gas abundance profile both depend on the assumed destruction rate: a higher rate leads to a lower peak abundance and a narrower peak, while a lower rate leads to the opposite. A rate of ~ 10⁻¹¹ s⁻¹ or higher is needed to decrease the average gaseous CO₂ abundance below the observational limit (see below). Increasing the rate moves the location of the iceline further out, as the total available CO₂ near the iceline decreases. The viscosity also influences the width of the abundance profile, higher viscosities lead to a broader abundance peak.

Fig. 6 Timeevolution series for a model without grain growth but with a constant destruction rate of gaseous CO2 of 10−11 s−1. This model assumes an α of 10−3 . Panels as in Fig. 3, but with different vertical scales.

Fig. 7 Abundance of CO2 as function of radii for models with different destruction rates of CO2. These modeluse α = 10−3 and uf = 0 m s−1 (left) and uf = 3 m s−1 (right). Abundance profiles after 1 Myr of evolution are shown. At this point the disk has reached a semi steady state. A destruction rate of at least 10−11 s−1 is needed to keep the averaged CO2 abundance below the observational limit. The iceline is further out for models with higher destruction rates as the destruction of gas-phase CO2 lowers the total abundance of CO2 within 10 AU where CO2 can sublimate.

4.4 Viscous evolution, grain growth and CO₂ destruction

The next step is to include CO₂ destruction in viscous evolution models with grain growth. Here only models that retain the disk dust mass, that is, models that have an overall gas-to-dust ratio smaller than 1000 after 1 Myr of evolution are considered since there is no observational evidence for very high gas-to-dust ratios. This means that the grains in our models do not reach pebble sizes such as have been inferred from observations (e.g. Pérez et al. 2012; Tazzari et al. 2016). It is unclear how these large grains are retainedin the disk as the radial drift timescale for these particles is expected to be shorter than the disk lifetime (e.g. Birnstiel et al. 2010; Krijt et al. 2016). The requirement of dust retention restricts our models to α =10⁻⁴ and u_f = 1 m s⁻¹, α =10⁻³ and u_f = 1, 3 m s⁻¹, and α = 10⁻² and u_f = 1, 3, 10 m s⁻¹.

The abundance profiles for different CO₂ destruction rates are shown in the right panel of Fig. 7 for the model with α = 10⁻³ and u_f =  3 m s⁻¹. The abundance profiles are determined after 1 Myr, when the disk is in a semi-steady state. A destruction rate higher than 10⁻¹² s⁻¹ will create an abundance profile with clear abundance maximum at the iceline, but for this rate 10⁻¹² s⁻¹ the disk surface area averaged CO₂ abundance is still higher than the observed value for both cases. Only for a rate of around 10⁻¹⁰ s⁻¹ do the models predict an abundance profile with a disk surface area averaged CO₂ abundance that is consistent with the high end of the observations. Only the models with destruction rate higher than 10⁻⁹ s⁻¹ have a peak in the abundance that is consistent with the maximal inferred abundances from the observations. As shown in Fig. 7 the model with a rate of 10⁻¹¹ s⁻¹ will create a spectrum consistent with the observations, even though the average CO₂ abundance is higher than the observed abundance limit.

For the selected models, the effect of different fragmentation velocities is small, with differences in peak abundances less than a factor of five between similar models, as seen in Fig. 7. This is unsurprising as the selection criteria for the models only consider models that have a transport velocity of solids that is less than an order of magnitude faster than the transport velocity of the model without grain growth and radial drift. For models with faster growth it becomes more arbitrary to give a representative abundance profile as semi steady state is not reached at any point in time before the disk is depleted of most of the dust and volatiles in the inner regions. However, for very young disks that may have very efficient radial drift of grains, a high gas-phase CO₂ abundance is expected within the CO₂ iceline unless the destruction rate is 10⁻¹¹ s⁻¹ or higher.

4.5 Model spectra

The spectralmodelling the focuses on the CO₂ bending mode around centred at 15 μm. This bending mode has a strong Q-branch that has been observed by Spitzer in protoplanetary disks (Carr & Najita 2008; Pontoppidan et al. 2010). This region also has a weaker feature due to the Q-branch of ¹³CO₂ at 15.42 μm which can be used, together with the ¹²CO₂ Q-branch, to infer information on the abundance structure (Bosman et al. 2017).

For the model with α = 10⁻³ and $R_{\text{destr},{\text{CO}}_2}={10}^{-11}$, the 2D abundance distribution assumed for the ray-tracing is shown in Fig. 8. The enhanced CO₂ abundance around the iceline shows up clearly. The results of the continuum ray-tracing can be found in Fig. 5 of Bosman et al. (2017). The temperature at 1 AU is 320 K. This is higher than the value assumed for the models in the previous sections and is in line with the high luminosity of the modelled source of 7.7 L_⊙. The higher temperaturemoves the CO₂ iceline slightly further out. The abundance profiles from the viscous model with dust growth and CO₂ destruction are still very similar to those from the previous section.

Figure 8 also shows the emitting region of the CO₂ Q(6) line for both isotopologues. This is one of the strongest lines in the CO₂ spectrum and part of the Q-branch. As such it is a good representative for the complete Q-branch. Due to its high abundance, the ¹²CO₂ lines have a high optical depth over a larger area in the disk. As such a large part of the disk contributes to the emission of the ¹² CO₂ lines. For ¹³ CO₂ the emission is from a more compact area, close to the peak in CO₂ abundance.

Figure 9 shows CO₂ spectra for different destruction rates and abundance floor of 10⁻⁸. Line-to-continuum ratios for the individual R and P branch linesare all greater than 0.01 and are detectable by JWST (Bosman et al. 2017). Note that only even J levels exist in the vibrational ground state, so P, Q and R branch lines exist only every other J for this band. Both the ¹²CO₂ and the ¹³CO₂ Q-branch fluxes are influenced by the different destruction rates. For a rate of 10⁻¹² s⁻¹ the ¹²CO₂ Q-branch is a factor of two brighter than for the higher rates, which give fluxes very similar in magnitude for the ¹² CO₂ Q-branch. For ¹³ CO₂ the model with a destruction rate of 10⁻¹² s⁻¹ is a least a factor of five brighter than the models with higher rates of 10⁻¹¹ s⁻¹ and 10⁻¹⁰ s⁻¹. Since the latter two ¹²CO₂ spectra are indistinguishable observationally, a conservative lower limit to the destruction rate of 10⁻¹¹ s⁻¹ is chosen. For ¹³CO₂ however, there is a factor of three difference between the two higher rate models; here the ¹³ CO₂ Q-branch flux is determined by subtracting the contribution from the ¹²CO₂ P(23)-line assuming it is the same as that of the neighbouring P(21) line.

The fluxratio of the ¹³CO₂ Q-branch over the ¹²CO₂ Q-branch is shown in Fig. 10 for a range of CO₂ destruction rates. Three regimes can be seen. For ratios larger than ~0.13 the CO₂ released from the ice dominates the abundance profile and spectrum, due to the relatively low destruction rate. The high CO₂ column, dueto the high abundance around the iceline sets the small flux ratio. At ratios smaller than 0.03, corresponding to destruction rates larger than 5 × 10⁻¹¹ s⁻¹, the CO₂ released from the ice at the iceline is destroyed so fast that it cannot leave an imprint on the Q-branch fluxes. In the region in between, both the CO₂ released at the iceline as well as the upper layer CO₂ influence the spectra, as such the flux ratio not only depends on the CO₂ destruction rate, but also on the CO₂ abundance floor used.

Note that Bosman et al. (2017) show that a high ¹³CO₂ Q-branch over ¹²CO₂ Q-branch flux ratio may also be an indication of a high abundance in the inner disk. More generally it can be said that a high flux ratio indicates that the CO₂ responsible for the emission is highly optically thick. In these cases the width of the ¹³ CO₂ Q-branch can be a measure for the temperature of this emitting gas, which can be related to the radial location of the origin of the emission. An analysis like this will probably only be possible with disk specific modelling.

The spectral modelling done here is meant as an illustration to what could be possible with JWST-MIRI observations and to show what kind of spectral features could hold clues to new insights into disk dynamics. As such the model setup is not entirely self consistent. Foremost, the temperature profiles calculated from the radiative transfer are not exactly the same as the temperatureprofile used in the viscous disk model. Furthermore, the starting point for our viscous evolution is rather arbitrary.

Another uncertainty is the vertical distribution of the CO₂ sublimated from the ices. The assumed turbulent viscosity would not only mix material radially, but also vertically, and the CO₂ that comes off the grain near the midplane should make it up to the upper layers of the disk where the infrared emission is generated.

The vertical dispersion of a species injected into the disk midplane grows with time as ${\sigma }_{\text{X}}=\sqrt{D_{\text{gas}}t}$ (Ciesla 2010). A steady state is reached if σ_X = H, the scale height of the disk, at which point a constant abundance is reached. As such we can estimate the time needed for vertical mixing. For a self-similar disk, this time scale becomes: ${\zeta }_{{\text{H}}_2}$(36)

If a similar mixing speed (α) is assumed for the radial and vertical processes, then vertical mixing should happen well within the lifetime of the disk for locations near the CO₂ iceline.

Given these uncertainties, a simple constant abundance up to a certain height has been chosen. However higher up in the disk, photodissociation by stellar radiation and processing by stellar X-rays become important. CO₂ is removed in our model in the harshest environments (A_V < 1 mag), but CO₂ photodissociation is still possible in the upper disk regions where CO₂ is included. The location and thickness of this layer is discussed in Sects. 5.2.5 and 5.2.6.

Fig. 8 2D CO2 abundance distribution from the $R_{\text{destr},{\text{CO}}_2}={10}^{-11}$ model with an abundance floor of 10−9. Black and white lines show the regions where 25% and 75% of the emission of respectively 12 CO 2 and 13 CO2 is emitted from. The 13CO2 emission is more concentrated around the region with the abundance enhancement.

5 Discussion

Viscous disk models including only freeze-out and desorption of CO₂ predict a higher CO₂ abundance than observed. This conclusion holds even for models without grain growth and radial drift. The discussion is divided into two sections, first the gas and grain surface chemistry is discussed to explore alternative chemical destruction routes. Subsequently physical effects are investigated, some of which also imply chemical effects to happen. Figure 11 shows the CO₂ abundance structure in the disk with some of the mechanisms discussed here that could obtain agreement between observations and models.

5.1 Chemical processes

In the scenario that the CO₂ released at the iceline is destroyed, an effective destruction rate of 10⁻¹¹ s⁻¹ or higher is required to match the observed and computed abundances. In Sects. 3.2.2 and 3.2.3 a number of known destruction processes for midplane CO₂ were discussed, but none of these chemical processes is able to account for the destruction rate of 10⁻¹¹ s⁻¹.

5.1.1 Clues from high-mass protostars?

Our finding that a CO₂ ice desorption model is inconsistent with CO₂ gas-phase observations is not the first to do so. For high mass protostars Boonman et al. (2003) noted a similar disagreement between simple in-fall and desorption models, which predict CO₂ gas-phase abundances around 10⁻⁵ whereas observations are more consistent with an abundance of ~10⁻⁷, a two orders of magnitude difference. Assuming that the CO₂ ice sublimates and is destroyed within the dynamical time of 10⁴ yr, a destruction rate of at least 10⁻¹¹ s⁻¹ is needed for protostars as well, of the same order as the rate needed in disks. If the CO₂ destruction mechanism is similar in both types of sources, then this would most likely be a chemical destruction mechanism with no or a weak dependence on total density since densities in protostellar envelopes are orders of magnitude lower than in the inner disk midplanes and the physical environments are very different. However, as argued below the existence of such a chemical pathway is unlikely, except for destruction by UV photons or X-rays.

5.1.2 Alternative gas phase destruction routes

Here we investigate whether a chemical pathway could be missing in the chemical networks. Since ion-molecule reactions initiated by He⁺ have already been discarded, the assumption is that this destruction comes from a neutral-neutral reactive collision. The rate for a two-body reaction can be written as: $$R_{{\text{CO}}_2,\text{destr}}=f_X{k}_{2\text{0body}}{n}_{\text{gas}}{x}_X$$(37)

where n_gas is the gas number density, k_{2- body} is the two-body reaction rate coefficient, x_X is the abundance of the collision partner and f_X is the branching ratio to any product that does not cycle back to CO₂ . Near the iceline densities are around 10¹³ cm⁻³, a typical rate coefficient for a barrier-less neutral-neutral reaction is around 10⁻¹¹ cm³ s⁻¹. As such the partner abundance needs to be higher than 10⁻¹³ to destroy CO₂ at a fast enough rate. The only species that have a nearly barrier-less reaction with CO₂ in either the UMIST (McElroy et al. 2013) or KIDA (Wakelam et al. 2012) databases are C, CH₂ and Si. Reaction rates for CO₂ with C to form two CO molecules and with CH₂ to form CO + H₂CO are 10⁻¹⁵ cm³ s⁻¹ and 3.9 × 10⁻¹⁴ cm³ s⁻¹ (Tsang & Hampson 1986), however, both rates are high temperature estimates and extrapolated to low temperatures (Hébrard et al. 2009). The rate for CO ₂ + Si →CO + SiO has been measured at room temperature to be 1.1 × 10⁻¹¹ cm³ s⁻¹ (Husain & Norris 1978), McElroy et al. (2013) assume a small energy barrier of 282 K, in line with the high temperature measurements. Other species in the databases that have an exothermic reaction with CO₂ are CH and N; both of these reactions are thought to have a barrier of 3000 and 1710 K respectively (Mitchell 1984; Avramenko & Krasnen’kov 1967). Other neutral-neutral reactions in the chemical networks are highly endothermic, for example CO ₂ + S →CO + SO is endothermic by ~20 000 K (Singleton & Cvetanovic 1988), effectively nullifying this reaction in the temperature range we are interested in.

C, CH₂ and Si are, due to their high reactivity, only present in very low quantities in the gas-phase with abundances at the solver accuracy limit of 10⁻¹⁷ for C and CH₂ , and 10⁻¹⁶ for Si and thus can hardly account for the destruction of a significant amount of CO₂ unless it is possible to quickly reform the initial reactant from the reaction products. The reactions products, CO, H ₂ CO and SiO, are again very stable, thus it is unlikely that these molecules quickly react further to form the initial reactants again. This means that the abundance of C, CH₂ and Si will go down with time if more and more CO₂ is added.

5.1.3 Alternative ice destruction routes

The low CO₂ abundance observed could also be explained if CO₂ is converted to other species while in the ice. Bisschop et al. (2007) report that CO₂ bombarded with H does not lead to HCOOH at detectable levels in their experiments. There are no data available on grain surfaces reactions of CO₂ ice with N, Si or C. If any of these reaction pathways were to destroy CO₂ in the ice, they should be very efficient, ~99% conversion, but at the same time be able to explain the high CO₂ ice abundances in comets. This quickly limits the temperature range in which these reactions can be efficient to 30–80 K.

In short, due to the high stability of CO₂, only a handful of highly reactive radicals or atoms can destroy CO₂ , whether in the gas or ice. None of these species is predicted to be abundant enough to quickly destroy the large flux of CO₂ that is brought into the inner disk due to dynamic processes of CO₂-containing icy grains.

5.2 Physical processes

The abundance profiles computed in this work all assume that the CO₂ can move freely through the disk, both in radial and vertical directions, while being shielded from UV and X-ray radiation. It is also assumed that all CO₂ ice sublimates at the CO₂ iceline and that CO₂ is highly abundant in the ice, as found in the ISM. Several physical processes could make these assumptions invalid. An important constraint is that due to the difference between the ice abundances (around 10⁻⁵) and the inferred gas abundance ( < 10⁻⁷) any process needs reduce the CO₂ abundance in the infrared emitting layers by at least two orders of magnitude.

5.2.1 Impact of reduced turbulence

The currentmodels assume a standard viscously spreading and evolving disk with a constant α value. Magneto-hydrodynamical (MHD) simulations currently favour disks with strongly changing viscosity in radial and vertical directions (e.g. Gammie 1996; Turner et al. 2007; Bai & Stone 2013). In these simulations turbulence is only expected to be high at those locations where the ionisation fraction of the gas is also high, that is, the upper and outer regions of the disk. In the disk mid-plane the viscous α may be as low as 10⁻⁶ , especially if non-ideal MHD effects are taken into account. In this regime, radial and vertical mixing due to the turbulence becomes a very slow process, with mixing time-scales close to the disk lifetime. Laminar accretion through the midplane would still enrich the inner disk midplane with CO₂ but the low vertical mixing would greatly reduce the amount of material that is lifted into the IR-emitting layers of the disk. A mixing time of 10⁶ yr would imply an effective vertical α of 10⁻⁶ for a CO₂ iceline at 10 AU (Eq. (36), Ciesla 2010). Note that this still means that CO₂ is mixed upwards within 10⁵ yr in the inner 1 AU, which would be visible as an added warmer component to the CO₂ emission. Observations of another species created primarily on grains at high abundance, > 10⁻⁷, such as CH₄, NH₃ or CH₃OH near their respective icelines can be used to rule out inefficient vertical mixing.

5.2.2 Trapping of CO₂ in large solids

There are ways to trap a volatile in a less volatile substance and prevent or delay its sublimation to higher temperatures. One possibility is to lock up the CO₂ inside the water ice. This can be done in the bulk of the amorphous ice, below the surface layers, or by the formation of clathrates in the water ice at the high midplane pressures. Even in the case that CO₂ is completely trapped in H₂O ice, CO₂ still comes off the grain together with the H₂O at the H₂ O iceline. At this point the CO₂ is in the gasphase and free to diffuse around in the disk. The strong negative abundance gradient outside the H₂ O iceline will transport the CO₂ outwards until the CO₂ starts to freeze-out of its own accord at the CO₂ iceline. This would still result in gas-phase abundances within the CO₂ iceline close to the original icy CO₂ abundances of the outer disk.

Spherical grains between 0.1 μm and 1 mm are coated with 10² to 10¹⁰ layers of water ice, respectively. In the case of small particles, with only 100 layers of ice, it is hard to imagine that the water traps a significant amount of CO₂ . For the larger grains with a smaller relative surface area, however, the water ice could readily trap a lot of CO₂ . In a fully mixed H₂O and CO₂ ice, the sublimation of CO₂ from the upper 100 million ice layers will still keep the fraction of sublimated CO₂ below 1%, however efficient grain fragmentation, such as assumed in our models, will expose CO₂ rich layers allowing for the sublimation of more CO₂. For small grains, clathrates may keep the CO₂ locked-up in the upper layers to prevent sublimation. However, to trap a single CO₂ molecule at least six H₂O molecules are needed (Fleyfel & Devlin 1991).

For amorphous ices similar restrictions exists on the H₂O to CO₂ ratio for trapping CO₂ in the water ice. For ices with a H₂O to CO₂ ratio of 5:1, CO₂ desorption happens at the CO₂ desorption temperature, only at a ratio of 20:1 is the CO₂ fully trapped within the H₂O ice (Sandford & Allamandola 1990). Collings et al. (2004) show that the majority of the trapped CO₂ desorbs at temperatures just below the desorption temperature of H₂O, due to the crystallization of water. Fayolle et al. (2011) note that the fraction of trapped CO₂ also depends on the thickness of the ice. They are able to trap 65% of the CO₂ in a 5:1 ice, which means that a significant fraction of the CO₂ still comes of at the CO₂ iceline.

Both of these mechanisms need to have all the CO₂ mixed with water. For interstellar dust grains there is observational evidence that this is not the case: high quality Spitzer spectra indicate that CO₂ is mixed in both the CO-rich layers of the ice as well as the H₂ O-rich layers of the ice (Pontoppidan et al. 2008). Sublimation of CO from these mixed ices would create a CO₂ rich layer around the mixed H₂ O:CO₂ inner layers.

Another option would be to lock the CO₂ in the refractory material. We note however that the CO₂ ice abundance is of the same order as the Mg, Fe and Si abundances (Asplund et al. 2009), as such it is unlikely that CO₂ can be locked up in the refractory material without chemical altering that material, such as conversion into, for example, CaCO₃ .

Measurements from 67P/Churyumov-Gerasimenko indicate that even a kilometre-sized object outgasses CO₂ at 3.5 AU distance (Hässig et al. 2015). For larger objects, CO₂ retention might be possible, butlocking up 99% of the CO₂ ice in these bodies seems highly unlikely.

5.2.3 Trapping icy grains outside of the iceline

It is possibleto stop the CO₂ from sublimating if the CO₂ ice never reaches the CO₂ iceline, for example by introducing a dust trap (van der Marel et al. 2013). Stopping the CO₂ ice from crossing the iceline would also mean stopping the small dust from crossing the iceline. This would lead to fast depletion of the dust in the inner disk, creating a transition disk, although some small dust grains can still cross the gap (Pinilla et al. 2016). Indeed, most disks for which CO₂ has been measured have strong near- and mid-infrared continuum emission, indicating that the inner dust disks cannot be strongly depleted. Furthermore, measurements by Kama et al. (2015) show that for Herbig Ae stars, the accreting material coming from a full disk has a refractory material content close to the assumed ISM values even though these disks have no detectable infrared CO₂ emission with Spitzer (Pontoppidan et al. 2010).

5.2.4 Shocks

For high-mass protostars, Charnley & Kaufman (2000) proposed that C-shocks with speeds above 30 km s⁻¹ passing through the region within the CO₂ iceline could account for the low observed CO₂ abundance. Dense C-shocks can elevate temperatures above 1000 K and create free atoms such as H, N and C, so all of the CO₂ in the shocked layer can then be efficiently converted to H₂O and CO. At the higher densities in disks, lower velocity shocks can reach the same temperature. For example, Stammler & Dullemond (2014) calculate that a 5 km s⁻¹ shock would reach temperatures upwards of 1000 K in the post shock gas. If shocks with such speeds occur in disks as well this could be an explanation for the low CO₂ abundance. However, assuming the shock only destroys gas-phase CO₂, not the ice, shocks must happen at least once every 3 × 10³ yr (1∕10⁻¹¹ s⁻¹) and the shock speeds have to be of the order of the Keplerian velocity.

These nebular shocks have been invoked to explain the ubiquitous chondrules found in meteorites. For these chondrules to form, dense gas needs to be flash heated, so that grains can evaporate and quickly reform as the gas cools (Hewins et al. 2005). Desch & Connolly (2002) proposed that a slow shock (7 km s⁻¹ ) in the dense midplane can heat up the gas to temperature around 2000 K, enough to start sublimation (see also Stammler & Dullemond 2014). In the case of efficient cooling, it takes the gas a few minutes to cool, but this would be enough for the CO ₂ + H →CO + OH reaction to significantly reduce the CO₂ abundance. Slower shocks, with low cooling rates would also be able to destroy CO₂ efficiently, as long as temperatures of 750 K are reached (see Fig. B.1). To effectively destroy CO₂, these shocksneed to affect the entire inner disk at least once every few thousand years. Each shock could then also induce chondrule formation. There is evidence that chondrules went through multiple creation events, but the frequency of these events is currently not constrained (Ruzicka et al. 2008; Desch et al. 2012).

Further out in the disk, shocks can sputter ices, including CO₂ ice, from the grains (Charnley & Kaufman 2000). However, sputtering does not necessarily destroy the CO₂ . Moreover, full destruction of ices by shocks is unlikely since CO₂ needs to survive in the comet formation regions of the disk.

5.2.5 UV-dominated layers

The inner disk chemical models by Walsh et al. (2015) show that there is a surface layer in the disk where UV destruction becomes important resulting in abundances of 10⁻⁸–10⁻⁷. Depending on the time-scale of the chemical processes and the extent of this layer, it may be optically thick in the ¹² CO₂ infrared lines. This means that low abundances of CO₂ will be measuredin ¹²CO₂ irrespective of the amount of sublimating CO₂ ice in the midplane. If this UV-dominated layer is not too thick, ¹³ CO₂ could still have an imprint of the deeper, more abundant CO₂ that is coming from the iceline. On the other hand, if this layer is thick enough that even ¹³ CO₂ becomes optically thick to its own radiation, then the more abundant CO₂ sublimating from the iceline would be hidden from our view (see Figs. 10 and 12). This would happen at CO₂ columns larger than 10¹⁸ cm⁻².

The (vertical) thickness of the photon dominated layer is very dependent on the UV flux and dust opacity in the upper layers of the disk. For a disk that has only small grains ( < 1 μm) in the upperatmosphere, the UV opacity is high. In this case the UV dominated layer only contains a short column of gas, and thus CO₂ (Fig. 13). Observations of such a disk would lead to the spectra as presented in Sect. 4.5. If it is assumed that grains have grown to mm sizes, even in the upper atmosphere, then the UV can penetrate a lot deeper. This leads to a larger column of UV dominated gas and that could hide the high abundance CO₂ in the midplane, if the CO₂ in this layer becomes optically thick to its own IR radiation or if the dust becomes optically thick at IR wavelengths in the UV dominated layer. The first scenario is sketched in panel C of Fig. 11. To probe the midplane CO₂ abundance one would need to observe disks with a lot of small grains in the upper atmosphere, preferably around UV weak sources such that the CO₂ column in the UV dominated layer is minimised.

Fig. 9 Spectra 12CO2 and 13 CO2 for models with different destruction rates at a spectral resolving power R = 2200 appropriate for JWST-MIRI. The abundance floor is taken to be 10−8. The black line shows the model results for a constant abundance as used in Bosman et al. (2017). The regions were the Q-branch flux is extracted is shown in the black lines. The 13CO2 Q-branch is blended with the 12CO2 P(23) line. The neighbouring P(21) line is indicated. Note that the models with destruction rates of 10−11 and 10−10 s−1 are indistinguishable in their 12CO2 spectra.

Fig. 10 Flux ratios of the 13CO2 and 12 CO2 Q-branches, determined by integrating the region shown in Fig. 9. The 13 CO 2 Q-branch flux is corrected for the emission of the P(23) line by subtraction of the neighbouring P(21) line flux.

Fig. 11 Schematic representation of the CO2 abundance near the iceline in the r–z plane. Panels B, C and D approximately represent the top half of the blue box in panel A. The straight green lines denote the τ = 1 surface for 12CO2 (solid) and 13CO2 (dashed). Panels B, C and D: the abundance structure under different mechanism to explain the difference in CO2 abundance between the predictions of the viscous model and the IR observations. Panel B: the case for a constant gas-phase destruction rate due to some as yet unidentified process (Sect. 4.4). Panel C: the abundance structure for deep UV/X-ray penetration (Sects. 5.2.5 and 5.2.6). In more extreme cases, the 13 CO 2 τ = 1 surface can also be in the UV/X-ray dominated region. Panel D: the abundance structure for negligible vertical mixing (Sect. 5.2.1). Image not to scale.

Fig. 12 Schematic of the vertical CO2 abundance (top panel) and upwards column (bottom panel) as function of disk height. Three different cases areshown in solid black, dashed dark blue and dotted blue lines. The black line shows the case for shallow UV or X-ray penetration, in this case an abundance similar to the midplane CO2 abundance will be measured. The dashed dark blue line shows a slightly deeper penetration of radiation, in which case a lower CO2 abundance will be inferred from the observations because the gas is optically thick and hiding the higher abundance deeper into the disk. The dotted blue line shows the case of deep penetration. In this case a low CO2 abundance will be inferred due to the dust hiding the bulk of the CO2. The vertical lines in the upper panel show where the disk becomes optically thick to CO2 lines, the solid and dashed grey line, show where the gas becomes optically thick for the small and intermediate PDR or XDR respectively, the green dotted line shows where the dust becomes optically thick, which is above where the gas becomes optically thick for the large PDR or XDR.

Fig. 13 Left: depth into the disk to which UV can destroy CO2 at a rate of 10−11 s−1 for a disk with an upper disk gas-to-dust ratio of 1000 and large grains (cyan) and small grains (blue). The black and grey contours show the area from which 25% and 75% of the emission of 12 CO 2 and 13 CO2 respectively originates. The green lines show the dust τ = 1 surface at 15 μm for the large grains (up to 1 mm, dark green) and small grains (up to 1 μm, light green) respectively. The UV photo destruction rate has a stronger dependence on the grain size distribution than this μm dust photosphere. Right: depth into the disk to which X-rays can destroy CO2 at a rate of 10−11 s−1 for different stellar X-ray luminosities (LX in erg s−1 ). The black and grey contours show the area from which 25% and 75% of the emission of 12 CO 2 and 13 CO2 respectively is coming from. The green line shows the dust τ = 1 surface at 15 μm for the large grains (up to 1 mm).

5.2.6 X-ray induced dissociation

T-Tauri stars are known to have strong X-ray emission 10²⁸ –10³⁴ erg s⁻¹ (Feigelson et al. 2002). These X-rays can ionise the gas and create a local UV field in the same way as cosmic-rays, destroying molecules in the gas (see, for example, the models by Stäuber et al. 2005, for high mass protostars). As X-rays are more quickly attenuated than the higher energy cosmic-rays, they are not efficient in destroying CO₂ in the disk midplane, however, X-rays do penetrate further into the disk atmosphere than UV photons. Thus, X-rays can help in decreasing the CO₂ abundance in the atmosphere of the disk, especially if hard X-ray flare, such are thought to occur in TW Hya are common in other T-Tauri stars as well (Kastner et al. 2002; Cleeves et al. 2015, 2017).

For our specific DALI model, a total X-ray flux of 10³⁰ erg s⁻¹, emitted as a 10⁷ K black body, would have a destruction rate of 10⁻¹¹ s⁻¹ in most of the ¹²CO₂ emitting region and could thus have an effect on the ¹²CO₂ fluxes (Fig. 13). A total X-ray flux of 10³² erg s⁻¹ would be needed to have the same destruction rate in most of the ¹³ CO₂ emitting region. Thus, a source that has a time-averaged X-ray luminosity > 10³⁰ erg s⁻¹ could show a low CO₂ abundance in the atmosphere and have a higher abundance underneath. If the X-ray luminosity of the source is lower than 10³² erg s⁻¹, this should be reflected in the ¹³CO₂ measurement. Note however that the X-rays do not reach into the midplane mass reservoir. Since the vertical transport for a fully viscous disk is orders of magnitude faster than the radial transport, the replenishment rate of CO₂ in the upper layers due to mixing from the midplane would also be faster. Thus, these X-ray luminosities should be considered as lower limits. Most disk hosting stars do not have high enough X-ray luminosities to destroy CO₂ efficiently (Bustamante et al. 2016), but X-ray destruction of CO₂ in the surface layers can be an explanation in some sources and lead to similar effects as for the case of UV, such as the CO₂ IR lines becoming optically thick in the X-ray dominated region.

6 Summary and conclusions

We have presented a model for the transport of major gas and ice species. This model predicts that transport of ices from the outer disk towards the inner disk can have a significant effect on the abundance of species in the inner disk. Radial transport is predicted to be so efficient, that even in the case without any radial drift, the modelled CO₂ abundance in the inner disk is orders of magnitude higher than the currently available observations of disks. As such CO₂ cannot be directly used to trace the mass transport rate. The presence of dust in disks older than a Myr does imply that the mass transport rate of dust cannot be larger than $10\times {\dot{M}}_{\text{acc},\text{gas}}/\left(\text{g/d}\right)$. This is inconsistent with the pebble-sized dust observed in disks.

To reconcile model and observation CO₂ abundances,either the midplane CO₂ needs to bedestroyed by some physical or chemical mechanism, or the CO₂ abundance in the upper layers needs to be decoupled from the CO₂ abundance in the midplane. This can be achieved by preventing vertical CO₂ cycling, hiding the CO₂ under the dust IR photosphere or by the creating a low abundance CO₂ column that is optically thick in a UV or X-ray dominated layer. We can summarize this work with the following conclusions.

• Accretion flow and diffusion equilibrate the inner disk CO₂ gas abundance with the high outer disk CO₂ ice abundance of ~10⁻⁵ in less than 1 Myr for α > 10⁻³ (Sect. 4.1).

• For any disk that retains its dust on the same time-scale as it retains the gas, dust diffusion processes dominate over radial drift (Sect. 4.2 and Appendix C).

• If dust can grow efficiently and drift inwards, both CO₂ and H₂O will be enhanced above the outer disk ice abundance as long as the dust is drifting inwards with a mass flux more than 1% of that of the gas, ${\zeta }_{{\text{H}}_2}$ (Sect. 4.4).

• A physical or chemical process that would induce a vertically averaged CO₂ destruction rate of at least 10⁻¹¹ s⁻¹ ( < 10³ yr) is necessary to produce ¹²CO₂ fluxes that are consistent with the current observations of disks (Sect. 4.5).

• Within the current chemical networks, there are no paths to destroy CO₂ fast enough without invoking strong X-ray or UV-fluxes or elevated cosmic ray fluxes. Especially at temperatures between 60 and 150 K, destroying CO₂ is extremely difficult in the current chemical networks since the formation of CO₂ is strongly favoured over the production of H₂O in that temperature range. Above 300 K CO₂ destruction pathways exists, but they only become efficient enough at temperatures above 750 K (Sects. 3.2, 5.1 and Appendix B).

• Shocks could raise temperatures enough to destroy CO₂ efficiently, but effectively lowering the CO₂ abundance would require the entire inner 10 AU of the disk to be processed by a shock every 10³ yr.

• Trapping CO₂ in water ice seems unlikely and is not expected to lower the inner disk gaseous CO₂ abundance significantly. Trapping the grains before they cross the iceline efficiently enough to lower CO₂ abundances in the inner disk would quickly create a transition disk (Sects. 5.2.2 and 5.2.3).

• Deep penetration of UV and strong X-ray fluxes can lower the CO₂ abundance in the ¹²CO₂ emitting region. This can happen if grains in the disk atmosphere have grown to mm sizes and the central source has enough UV luminosity ( > 10²⁷ erg s⁻¹, about 0.1% of the solar UV). Likewise sources that have an X-ray luminosity greater than 10³⁰ erg s⁻¹ efficiently destroy CO₂ in a large part of the IR emitting region. If the UV or X-ray dominated region does not reach down to the dust 15 μm τ = 1 surface, ¹³CO₂ can still contain a information on the amount of sublimated CO₂.

• A high ratio of the ¹³CO₂ Q-branch flux over the ¹²CO₂ Q-branch flux ( >0.05), can hint at an contribution of sublimating CO₂ at the iceline. Sublimating CO₂ is, however, not the only way to get high flux ratios. A more complete, source specific analysis would be needed to exclude other possibilities (Sect. 4.5).

In summary, to explain the discrepancy between our models and observations, we either need to invoke frequent shocks in the inner 10 AU, or need to hide the abundant midplane CO₂ from our view. This can be due to very inefficient mixing or due to deep UV and/or X-ray penetration. In the case of UV penetration, this would imply efficient grain growth and settling, in the case of X-rays, this would be stronger than measured X-ray fluxes for most sources, or frequent X-ray flares. The key to hiding the mid-plane CO₂ is that the low abundance CO₂ in the upper layers has a high enough column density to become optically thick to its own IR radiation. Trapping the CO₂ on the grains, or destroying it in the gas phase near the iceline at ambient temperatures is strongly disfavoured.

Is CO₂ a special case in terms of discrepancy between models and observations? Modelling and observations of other molecules abundant in ices, such as NH₃ and CH₄ , may be used to test whether a chemical or physical reason is the origin of the current discrepancy. Ultimately, such data can constrain the mass flow of solids if the chemical composition of the surface layers is sufficiently representative of that in the midplane.

Acknowledgements

Astrochemistry in Leiden is supported by the European Union A-ERC grant 291141 CHEMPLAN, by the Netherlands Research School for Astronomy (NOVA), by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize. Daniel Harsono, Stefano Facchini, Kees Dullemond and Til Birnstiel are thanked for insightful discussions and Catherine Walsh and Christian Eistrup for sharing and help with the gas-grain chemical network.

Appendix A UV dust cross sections

Figure A.1 shows the UV opacities and cross sections for different grain size distributions used in Sect. 3.2.2.

Fig. A.1 UV opacities and cross section of the dust for different sizes of the largest grains. The MRN power law slope is 3.5 for all cases. Growing grains from 1 μm to 1 mm hasmore than an order of magnitude effect on the opacities.

Appendix B Chemical modelling

B.1 Gas-phase only models

To test the main CO₂ destruction and formation processes, a small grid of chemical models was run. The network from Walsh et al. (2014) was used as basis and modified. The gas-phase chemical network is based on the UMIST 2012 network (McElroy et al. 2013). The network includes freeze-out and desorption. H₂ formation on grains following Cazaux & Tielens (2004), which includes H₂ formation at high temperatures. The grain-surface reactions included in the network were removed for this test. As the main interest is inthe midplane layers, the local UV and X-ray flux are assumed to be negligible. The cosmic ray rate is taken to be 10⁻¹⁷ s⁻¹.

Fig. B.1 CO2 abundance from a gas-phase chemical network as function of temperature after 1 Myr of evolution. The density is 1012 cm−3.

Four series of models were run, two starting from “atomic” initial conditions (all species are in atomic form except for H₂ ), the other two from “molecular” conditions (see Table B.1, the atomic initial conditions use the same elemental abundances as the used for the molecular test). For each set of initial conditions one set of models was run with H₂ formation happening only on grains according to Cazaux & Tielens (2004) leading to an atomic H abundance of 10⁻⁴ with respect to H₂ , while in the other set of models the additional H₂ formation route at high temperature was added to keep the H abundance at 10⁻¹² (1 cm⁻³ ). The density was kept constant at 10¹² cm⁻³ and the temperature was varied between 150 and 800 K.

The CO₂ abundances after 1 Myr of chemical evolution are presented in Fig. B.1. From this plot it can be gathered that there is no efficient destruction mechanism in the current gas-phase network until a temperature of ~450 K is reached.

For temperatures between 150 and 450 K, the different atomic H abundances do not significantly influence the results. The difference between the initial CO₂ abundance and the CO₂ abundance after 1 Myr in the molecular case can be explained by the destruction of CO₂ and H₂ O due to processes driven by cosmic rays, e.g. cosmic ray induced photo dissociation, and the reformation of CO₂ from the fragments.

For models starting from atomic initial conditions the CO₂ content decreases with increasing temperature as expected from the formation speeds of CO₂ and H₂ O in the model with a low H abundance, and this trend continues up to higher temperatures. The local minimum in the CO₂ abundance at 350 K is dueto reactions with atomic N before the nitrogen is locked up in other species. When the atomic nitrogen abundance is high, CO₂ reacts with N to form CO and NO, NO is then further processed to H₂O, releasing a N in the process. This significantly slows down the build-up of CO₂ when both O and OH are available in the gas-phase in high abundances.

In the models with a high H abundance, reactions with atomic hydrogen become important starting at 350 K. This process, which is slightly more efficient for H₂O than for CO₂ , together with the reformation of H₂O and CO₂ from the resulting OH, sets the equilibrium at high temperature. The abundance of H₂ O as function of temperature is fairly constant, at 10⁻⁴, and thus the CO₂ abundance decreases towards the highest temperature, as expected from the formation rate ratio given in Fig. 1. The difference between the CO₂ abundance at high temperature in the case of atomic and molecular abundance can be explained by the lower CO abundance in the case of atomic initial conditions.

For the models with a low H abundance, CO₂ is not efficiently destroyed by H at high temperatures so the CO₂ abundance is nearly constant as function of temperature. From this parameter exploration we can conclude that CO₂ can only be effectively destroyed by radicals like H and N and that this is only efficient if these are available at high abundances > 10⁻⁶. For the model beginning with atomic abundances the low temperature dependence of the abundance is similar to the high H abundance case. Above temperatures of 400 K differences start to be noticeable. The CO₂ abundance does not increase with temperature as it does in the case of high H abundance. This is due to a similar process as in the high H model around 350 K. All the CO₂ that is formed when oxygen is not locked up in water is quickly destroyed by atomic H and N, when atomic H and N have settled to low enough abundances that CO₂ can survive, all the O is already locked up in other species, mostly H₂O and CO leaving no extra O to form CO₂. The small amount of CO₂ that is in the gas-phase is mostly created by cosmic ray induced destruction of H₂ O. The decreasing abundance towards high temperatures is due to the increasingly efficient reformation of H₂ O after destruction due to cosmic rays.

B.2 Grain surface chemistry between the H₂O and CO₂ icelines

Section B.1 investigates the gas-phase chemistry of CO₂ in the region within the H₂ O iceline, that is, above 150K. In this region all species will be primarily in the gas-phase and grain surface chemistry is thus expected to be less important. Between theCO₂ and H₂ O icelines The efficiency of the destruction of gaseous CO₂ is estimatedby investigating the chemical pathways in a full chemical model. The network from Walsh et al. (2014) and Eistrup et al. (2016) was used. A point model representative for the conditions between the H₂ O and CO₂ iceline was run. The number density used is 10¹³ cm⁻³, the temperature 100 K and the cosmic ray ionisation rate 10⁻¹⁷ s⁻¹. There was no external UV in our model. After running the model for 1 Myr to get equilibrium abundances of atomic and ionised species, the chemical pathways of oxygen carrying species are examined.

For each oxygen carrying molecule that could be created from a destruction of CO₂ , we define the destruction efficiency as the fraction the oxygen atoms that is incorporated into a stable molecule, for which we take that a cosmic ray induced process is needed to destroy the molecule. Here, pathways leading to incorporation of oxygen into CO are neglected, as these pathways are dominated by the availability of atomic C, which must originate from an other CO molecule. Pathways towards CO are generally very minor. Destruction efficiencies are tabulated in Table B.2.

From the efficiencies it can be seen that direct ionisation or protonation of CO₂ is not an effective way to destroy itsince CO₂ will be reformed from the products in most cases. Only if an oxygen atom can be removed from the central carbon, is there a chance that the oxygen will be locked up into something other than CO₂ . In these cases the efficiency is 27% (for OH) or higher. Appendix C:

Apppendix C Viscous evolution and grain growth

Figures C.1 through C.9 show model results for the model with grain growth and radial drift, but without any destruction rate. α values of 10⁻², 10⁻³ and 10⁻⁴ are shown, fragmentation velocities of 1, 10 and 100 m s⁻¹ are used. The features seen in ice abundance for the higher fragmentation velocity models is due to a numerical instability in the grain growth algorithm, this does not influence our results.

Fig. C.1 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.

Fig. C.2 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.

Fig. C.3 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.

Fig. C.4 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.

Fig. C.5 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 10 m s−1. Panels as in Fig. 3.

Fig. C.6 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 100 m s−1. Panels as in Fig. 3.

Fig. C.7 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.

Fig. C.8 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 10 m s−1. Panels as in Fig. 3.

Fig. C.9 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 100 m s−1. Panels as in Fig. 3.

References

All Tables

All Figures

	Fig. 1 Ratio of the CO2 to H2 O formation rate from a reaction of OH with CO and H2 respectively. Hydrogenation of OH to H2O dominates above 150 K. Formation times of CO2 and H2 O from OH are faster than the inner disk mixing time for all temperatures considered here.
In the text

	Fig. 2 Dissociation rate of CO2 due to cosmic ray induced photons, for different dust size distributions. Solid lines are for condition inside of the H2 O iceline, dashed lines for conditions outside of the H2O iceline. The efficiency multiplied by ${\zeta }_{{\text{H}}_2}$ gives the CO2 dissociation rate.
In the text

	Fig. 3 Time evolution series for a model without grain growth. This model assumes an α of 10−3 . Left: surface density of the gas (solid lines) and solids (dashed lines), the solid surface density is the sum of the dust and ice surface densities, it has been multiplied by a factor of 10 for visualization. Middle: abundance of CO2 . Right: abundance of H2O.
In the text

	Fig. 4 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 3 m s−1. Panels as in Fig. 3.
In the text

	Fig. 5 Time evolution series for a model with grain growth. The fragmentation velocity for this model is 30 m s−1 and α = 10−3. The verticalscales are different from Figs. 3 and 4. Panels as in Fig. 3.
In the text

	Fig. 6 Timeevolution series for a model without grain growth but with a constant destruction rate of gaseous CO2 of 10−11 s−1. This model assumes an α of 10−3 . Panels as in Fig. 3, but with different vertical scales.
In the text

Fig. 7 Abundance of CO2 as function of radii for models with different destruction rates of CO2. These modeluse α = 10−3 and uf = 0 m s−1 (left) and uf = 3 m s−1 (right). Abundance profiles after 1 Myr of evolution are shown. At this point the disk has reached a semi steady state. A destruction rate of at least 10−11 s−1 is needed to keep the averaged CO2 abundance below the observational limit. The iceline is further out for models with higher destruction rates as the destruction of gas-phase CO2 lowers the total abundance of CO2 within 10 AU where CO2 can sublimate.

In the text

	Fig. 8 2D CO2 abundance distribution from the $R_{\text{destr},{\text{CO}}_2}={10}^{-11}$ model with an abundance floor of 10−9. Black and white lines show the regions where 25% and 75% of the emission of respectively 12 CO 2 and 13 CO2 is emitted from. The 13CO2 emission is more concentrated around the region with the abundance enhancement.
In the text

Fig. 9 Spectra 12CO2 and 13 CO2 for models with different destruction rates at a spectral resolving power R = 2200 appropriate for JWST-MIRI. The abundance floor is taken to be 10−8. The black line shows the model results for a constant abundance as used in Bosman et al. (2017). The regions were the Q-branch flux is extracted is shown in the black lines. The 13CO2 Q-branch is blended with the 12CO2 P(23) line. The neighbouring P(21) line is indicated. Note that the models with destruction rates of 10−11 and 10−10 s−1 are indistinguishable in their 12CO2 spectra.

In the text

	Fig. 10 Flux ratios of the 13CO2 and 12 CO2 Q-branches, determined by integrating the region shown in Fig. 9. The 13 CO 2 Q-branch flux is corrected for the emission of the P(23) line by subtraction of the neighbouring P(21) line flux.
In the text

Fig. 11 Schematic representation of the CO2 abundance near the iceline in the r–z plane. Panels B, C and D approximately represent the top half of the blue box in panel A. The straight green lines denote the τ = 1 surface for 12CO2 (solid) and 13CO2 (dashed). Panels B, C and D: the abundance structure under different mechanism to explain the difference in CO2 abundance between the predictions of the viscous model and the IR observations. Panel B: the case for a constant gas-phase destruction rate due to some as yet unidentified process (Sect. 4.4). Panel C: the abundance structure for deep UV/X-ray penetration (Sects. 5.2.5 and 5.2.6). In more extreme cases, the 13 CO 2 τ = 1 surface can also be in the UV/X-ray dominated region. Panel D: the abundance structure for negligible vertical mixing (Sect. 5.2.1). Image not to scale.

In the text

Fig. 12 Schematic of the vertical CO2 abundance (top panel) and upwards column (bottom panel) as function of disk height. Three different cases areshown in solid black, dashed dark blue and dotted blue lines. The black line shows the case for shallow UV or X-ray penetration, in this case an abundance similar to the midplane CO2 abundance will be measured. The dashed dark blue line shows a slightly deeper penetration of radiation, in which case a lower CO2 abundance will be inferred from the observations because the gas is optically thick and hiding the higher abundance deeper into the disk. The dotted blue line shows the case of deep penetration. In this case a low CO2 abundance will be inferred due to the dust hiding the bulk of the CO2. The vertical lines in the upper panel show where the disk becomes optically thick to CO2 lines, the solid and dashed grey line, show where the gas becomes optically thick for the small and intermediate PDR or XDR respectively, the green dotted line shows where the dust becomes optically thick, which is above where the gas becomes optically thick for the large PDR or XDR.

In the text

Fig. 13 Left: depth into the disk to which UV can destroy CO2 at a rate of 10−11 s−1 for a disk with an upper disk gas-to-dust ratio of 1000 and large grains (cyan) and small grains (blue). The black and grey contours show the area from which 25% and 75% of the emission of 12 CO 2 and 13 CO2 respectively originates. The green lines show the dust τ = 1 surface at 15 μm for the large grains (up to 1 mm, dark green) and small grains (up to 1 μm, light green) respectively. The UV photo destruction rate has a stronger dependence on the grain size distribution than this μm dust photosphere. Right: depth into the disk to which X-rays can destroy CO2 at a rate of 10−11 s−1 for different stellar X-ray luminosities (LX in erg s−1 ). The black and grey contours show the area from which 25% and 75% of the emission of 12 CO 2 and 13 CO2 respectively is coming from. The green line shows the dust τ = 1 surface at 15 μm for the large grains (up to 1 mm).

In the text

	Fig. A.1 UV opacities and cross section of the dust for different sizes of the largest grains. The MRN power law slope is 3.5 for all cases. Growing grains from 1 μm to 1 mm hasmore than an order of magnitude effect on the opacities.
In the text

	Fig. B.1 CO2 abundance from a gas-phase chemical network as function of temperature after 1 Myr of evolution. The density is 1012 cm−3.
In the text

	Fig. C.1 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.2 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.3 Time evolution series for a model with grain growth. This model assumes an α of 10−2 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.4 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.5 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 10 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.6 Time evolution series for a model with grain growth. This model assumes an α of 10−3 . The fragmentation velocity for this model is 100 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.7 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 1 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.8 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 10 m s−1. Panels as in Fig. 3.
In the text

	Fig. C.9 Time evolution series for a model with grain growth. This model assumes an α of 10−4 . The fragmentation velocity for this model is 100 m s−1. Panels as in Fig. 3.
In the text