© ESO, 2017

1. Introduction

Most observed exo-planets orbit close to their parent star (for a review see: Udry & Santos 2007; Winn & Fabrycky 2015). The atmospheres of these close-in planets show a large diversity in molecular composition (Madhusudhan et al. 2014). This diversity in molecular composition must be set during planet formation and thus be representative of the natal protoplanetary disk. Understanding the chemistry of the inner, planet-forming regions of circumstellar disks around young stars will thus give us another important piece of the puzzle of planet formation. Prime molecules for such studies are H₂O, CO, CO₂ and CH₄ which are the major oxygen- and carbon-bearing species that set the overall C/O ratio (Öberg et al. 2011).

The chemistry in the inner disk, that is, its inner few AU, differs from that in the outer disk. It lies within the H₂O and CO₂ icelines so all icy planetesimals are sublimated. The large range of temperatures (100–1500 K) and densities (10¹⁰−10¹⁶ cm^-3) then makes for a diverse chemistry across the inner disk region (see e.g. Willacy et al. 1998; Markwick et al. 2002; Agúndez et al. 2008; Henning & Semenov 2013; Walsh et al. 2015). The driving cause for this diversity is high temperature chemistry: some molecules such as H₂O and HCN have reaction barriers in their formation pathways that make it difficult to produce the molecule in high abundances at temperatures below a few hundred Kelvin. As soon as the temperature is high enough to overcome these barriers, formation is fast and these molecules become major reservoirs of oxygen and nitrogen. An interesting example is formed by the main oxygen bearing molecules, H₂O and CO₂: the gas phase formation of both these molecules includes the OH radical. At temperatures below ~ 200 K the formation of CO₂ is faster, leading to high gas phase abundances, up to ~ 10^-6 with respect to (w.r.t.) total gas density, in regions where CO₂ is not frozen out. When the temperature is high enough, H₂O formation will push most of the gas phase oxygen into H₂O and the CO₂ abundance drops to ~ 10^-8 (Agúndez et al. 2008; Walsh et al. 2014, 2015). Such chemical transitions can have strong implications for the atmospheric content of gas giants formed in these regions if most of their atmosphere is accreted from the surrounding gas.

A major question is to what extent the inner disk abundances indeed reflect high temperature chemistry or whether continuously migrating and sublimating icy planetesimals and pebbles at the icelines replenish the disk atmospheres (Stevenson & Lunine 1988; Ciesla & Cuzzi 2006). Interstellar ices are known to be rich in CO₂, with typical abudances of 25% w.r.t. H₂O ice, or about 10^-5 w.r.t. total gas density (de Graauw et al. 1996; Gibb et al. 2004; Bergin et al. 2005; Pontoppidan et al. 2008; Boogert et al. 2015). Cometary ices show similarly high CO₂/H₂O abundance ratios (Mumma & Charnley 2011; Le Roy et al. 2015). Of all molecules with high ice abundances, CO₂ shows the largest contrast between interstellar ice and high temperature chemistry abundances, and could therefore be a good diagnostic of its chemistry. Pontoppidan & Blevins (2014) argue based on Spitzer Space Telescope data that CO₂ is not inherited from the interstellar medium but is reset by chemistry in the inner disk. However, that analysis used a local thermodynamic equilibrium (LTE) CO₂ excitation model coupled with a disk model and did not investigate the potential of future instruments, which could be more sensitive to a contribution from sublimating planetesimals. Here we re-consider the retrieval of CO₂ abundances in the inner regions of protoplanetary disks using a full non-LTE excitation and radiative transfer disk model, with a forward look to the new opportunities offered by the James Webb Space Telescope (JWST).

The detection of infrared vibrational bands seen from CO₂, C₂H₂ and HCN, together with high energy rotational lines of OH and H₂O, was one of the major discoveries of the Spitzer Space Telescope (e.g. Lahuis et al. 2006; Carr & Najita 2008, 2011; Salyk et al. 2008, 2011b; Pascucci et al. 2009, 2013; Pontoppidan et al. 2010; Najita et al. 2011). These data cover wavelengths in the 10–35 μm range at low spectral resolving power of λ/ Δλ = 600. Complementary ground-based infrared spectroscopy of molecules such as CO, OH, H₂O, CH₄, C₂H₂ and HCN also exists at shorter wavelengths in the 3–5 μm range (e.g. Najita et al. 2003; Gibb et al. 2007; Salyk et al. 2008, 2011a; Fedele et al. 2011; Mandell et al. 2012; Gibb & Horne 2013; Brown et al. 2013). The high spectral resolving power of R = 25 000−10⁵ for instruments like Keck/NIRSPEC and VLT/CRIRES have resolved the line profiles and have revealed interesting kinematical phenomena, such as disk winds in the inner disk regions (Pontoppidan et al. 2008, 2011; Bast et al. 2011; Brown et al. 2013). Further advances are expected with VLT/CRIRES+ as well as through modelling of current data with more detailed physical models.

Protoplanetary disks have a complex physical structure (see Armitage 2011, for a review) and putting all physics, from magnetically induced turbulence to full radiative transfer, into a single model is not feasible. This means that simplifications must be made. During the Spitzer era, the models used to explain the observations were usually LTE excitation slab models at a single temperature. With 2D physical models such as RADLITE (Pontoppidan et al. 2009) and with full 2D physical-chemical models such as Dust and Lines (DALI, Bruderer et al. 2012; Bruderer 2013) or Protoplanetary Disk Model (ProDiMo, Woitke et al. 2009) it is now possible to fully take into account the large range of temperatures and densities as well as the non-local excitation effects. For example, it has been shown that it is important to include radiative pumping introduced by hot (500–1500 K) thermal dust emission of regions just behind the inner rim. This has been done for H₂O by Meijerink et al. (2009) who concluded that to explain the mid-infrared water lines observed with Spitzer, water is located in the inner ~1 AU in a region where the local gas-to-dust ratio is 1–2 orders of magnitude higher than the interstellar medium (ISM) value. Antonellini et al. (2015, 2016) performed a protoplanetary disk parameter study to see how disk parameters affect the H₂O emission. Mandell et al. (2012) compared an LTE disk model analysis using RADLITE with slab models and concluded that, while inferred abundance ratios were similar with factors of a few, there could be orders of magnitude differences in absolute abundances depending on the assumed emitting area in slab models (see also discussion in Salyk et al. 2011b). Thi et al. (2013) concluded that the CO infrared emission from disks around Herbig stars was rotationally cool and vibrationally hot due to a combination of infrared and ultraviolet (UV) pumping fields (see also Brown et al. 2013). Bruderer et al. (2015) modelled the non-LTE excitation and emission of HCN concluding that the emitting area for mid-infrared lines can be ten times larger in disks than the assumed emitting area in slab models due to infrared pumping. Our study of CO₂ is along similar lines as that for HCN.

As CO₂ cannot be observed through rotational transitions in the far-infrared and submillimeter, because of the lack of a permanent dipole moment, it must be observed through its vibrational transitions at near- and mid-infrared wavelengths. The CO₂ in our own atmosphere makes it impossible to detect these CO₂ lines from astronomical sources from the ground, and even at altitudes of 13 km with SOFIA. This means that CO₂ has to be observed from space. CO₂ has been observed by Spitzer in protoplanetary disks through its v₂Q-branch at 15 μm where many individual Q-band lines combine into a single broad Q-branch feature at low spectral resolution (Lahuis et al. 2006; Carr & Najita 2008). These gaseous CO₂ lines have first been detected in high mass protostars and shocks with the Infrared Space Observatory (ISO, e.g. van Dishoeck et al. 1996; Boonman et al. 2003a,b). CO₂ also has a strong band around 4.3 μm due to the v₃ asymmetric stretch mode. This mode has high Einstein A coefficients and thus should thus be easily observable, but has not been seen from CO₂ gas towards protoplanetary disks or protostars, in contrast with the corresponding feature in CO₂ ice (van Dishoeck et al. 1996).

The CO₂v₂Q-branch profile is slightly narrower than that of C₂H₂ and HCN observed at similar wavelengths. These results suggest that CO₂ is absent (or strongly under-represented) in the inner, hottest regions of the disk. Full disk LTE modelling of RNO 90 by Pontoppidan & Blevins (2014) using RADLITE showed that the observations of this disk favour a low CO₂ abundance (10^-4 w.r.t. H₂O, ≈10^-8 w.r.t. total gas density). The slab models by Salyk et al. (2011b) indicate smaller differences between the CO₂ and H₂O abundances, although CO₂ is still found to be 2 to 3 orders of magnitude lower in abundance.

To properly analyse CO₂ emission from disks, a full non-LTE excitation model of the CO₂ ro-vibrational levels must be made, using molecular data from experiments and detailed quantum calculations. This model can then be used to perform a simple slab model study to see under which conditions non-LTE effects may be important. These same slab model tests are also used to check the influences of the assumptions made in setting up the ro-vibrational excitation model. Such CO₂ models have been developed in the past for evolved asymptotic giant branch stars (e.g. Cami et al. 2000; González-Alfonso & Cernicharo 1999) and shocks (e.g. Boonman et al. 2003b), but not applied to disks.

Our CO₂ excitation model is coupled with a full protoplanetary disk model computed with DALI to investigate the importance of non-LTE excitation, infrared pumping and dust opacity on the emission spectra. In addition, the effects of varying some key disk parameters such as source luminosity and gas/dust ratios on line fluxes and line-to-continuum ratios are investigated. Finally, Spitzer data for a set of T Tauri disks are analysed to derive the CO₂ abundance structure using parametrized abundances.

JWST will allow a big leap forward in our observing capabilities at near- and mid-infrared wavelengths, where the inner planet-forming regions of disks emit most of their lines. The spectrometers on board JWST, NIRSPEC and MIRI (Rieke et al. 2015) with their higher spectral resolving power (R ≈ 3000) compared to Spitzer (R = 600) will not only separate many blended lines (Pontoppidan et al. 2010) but also boost line-to-continuum ratios allowing detection of individual P, Q and R-branch lines thus giving new information on the physics and chemistry of the inner disk. Here we simulate the emission spectra of CO₂ and its ¹³CO₂ isotopologue from a protoplanetary disk at JWST resolution. We investigate which subset of these lines is the most useful for abundance determinations at different disk heights and point out the importance of detecting the ¹³CO₂ feature. We also investigate which features could signify high CO₂ abundances around the CO₂ iceline due to sublimating planetesimals.

2. Modelling CO₂ emission

2.1. Vibrational states

The structure of a molecular emission spectrum depends on the vibrational level energies and transitions between these levels that can be mediated by photons. Figure 1 shows the vibrational energy level diagram for CO₂ from the HITRAN database (Rothman et al. 2013). Lines denote the transitions that are dipole allowed. Colours denote the part of the spectrum where features will show up. This colour coding is repeated in Fig. 2 where a model CO₂ spectrum is presented.

CO₂ is a linear molecule with a ground state. It has a symmetric, v₁, and an asymmetric, v₃, stretching mode (both of the Σ type) and a doubly degenerate bending mode, v₂ (Π type) with an angular momentum, l. A vibrational state is denoted by these quantum numbers as: . The vibrational constant of the symmetric stretch mode is very close to twice that of the bending mode. Due to this resonance, states with the same value for 2v₁ + v₂ and the same angular momentum mix. This mixing leads to multiple vibrational levels that have different energies in a process known as Fermi splitting. The Fermi split levels have the same notation as the unmixed state with the highest symmetric stretch quantum number, v₁ and numbered in order of decreasing energy¹. This leads to the vibrational state notation of: where n is the numbering of the levels. This full designation is used in Fig. 1. For the rest of the paper we will drop the (n) for the levels where there is only one variant.

The number of vibrational states in the HITRAN database is much larger than the set of states used here. Not all of the vibrational states are needed to model CO₂ in a protoplanetary disk because some the higher energy levels can hardly be excited, either collisionally or with radiation, so they should not have an impact on the emitted line radiation. We adopt the same levels as used for AGB stars in González-Alfonso & Cernicharo (1999) and add to this set the 03³0 vibrational level.

2.2. Rotational ladders

The rotational ladder of the ground state is given in Fig. 1. All states up to J = 80 in each vibrational state are included; this rotational level corresponds to an energy of approximately 3700 K (2550 cm^-1) above the vibrational state energy. The rotational structure of CO₂ is more complex than that of a linear diatomic like CO. This is due to the fully symmetric wavefunction of CO₂ in the ground electronic state. This means that all states of CO₂ need to be fully symmetric to satisfy Bose-Einstein statistics. As a result, not all rotational quantum numbers J exist in all of the vibrational states: some vibrational states miss all odd or all even J levels. There are also additional selections on the Wang parity of the states (e,f). For the ground vibrational state this means that only the rotational states with even J numbers are present and that the parity of these states is e.

The rotational structure is summarized in Table 1. The states with v₂ = v₃ = 0 all have the same rotational structure as the ground vibrational state. The 01¹0(1) state has both even and odd J levels starting at J = 1. The even J levels have f parity, while the odd J levels have e parity. In general for levels with v₂ ≠ 0 and v₃ = 0, the rotational ladder starts at J = v₂ with an even parity, with the parity alternating in the rotational ladder with increasing J. For v₃ ≠ 0 and v₂ = 0, only odd J levels exist if v₃ is odd, whereas only even J levels exist if v₃ is even. All levels have an e parity. For v₂ ≠ 0 and v₃ ≠ 0, the rotational ladder is the same as for the v₂ ≠ 0 and v₃ = 0 case if v₃ is even, whereas the parities relative to this case are switched if v₃ is odd.

2.3. Transitions between states

To properly model the emission of infrared lines from protoplanetary disks non-LTE effects must be taken into account. The population of each level was determined by the balance of the transition rates, both radiative and collisional. The radiative transition rates were set by the Einstein coefficients and the ambient radiation field. Einstein coefficients for CO₂ have been well studied, both in the laboratory and in detailed quantum chemical calculations (see e.g. Rothman et al. 2009; Jacquinet-Husson et al. 2011; Rothman et al. 2013; Tashkun et al. 2015, and references therein). These are collected in several databases for CO₂ energy levels and Einstein coefficients such as the Carbon Dioxide Spectroscopic Database (CDSD) (Tashkun et al. 2015) and as part of large molecular databases such as HITRAN (Rothman et al. 2013) and GEISA (Jacquinet-Husson et al. 2011). Here the ¹²CO₂ and ¹³CO₂ data from the HITRAN database were used. It should be noted that the differences between the three databases are small for the lines considered here, within a few % in line intensity and less than 1% for the line positions.

The HITRAN database gives the energies of the ro-vibrational levels above the ground state and the Einstein A coefficients of transitions between them. Only transitions above a certain intensity at 296 K are included in the databases. The weakest lines included in the line list are 13 orders of magnitude weaker than the strongest lines. With expected temperatures in the inner regions of disks ranging from 100–1000 K, no important lines should be missed due to this intensity cut. In the final, narrowed down set of states all transitions that are dipole allowed have been accounted for.

Collisional rate coefficients between vibrational states are collected from literature sources. The measured rate of the relaxation of the 01¹0 to the 00⁰0 state by collisions with H₂ from Allen et al. (1980) is used. Vibrational relaxation of the 00⁰1 state due to collisions with H₂ is taken from Nevdakh et al. (2003). For the transitions between the Fermi split levels the rate by Jacobs et al. (1975) for collisions between CO₂ with CO₂ is used with a scaling for the decreased mean molecular mass. Although data used here supersede those in Taylor & Bitterman (1969), that paper does give a sense for the uncertainties of the experiments. The different experiments in Taylor & Bitterman (1969) usually agree within a factor of two, and the numbers used here from Allen et al. (1980) and Nevdakh et al. (2003) fall within the spread for their respective transitions. It is thus expected that the accuracy of the individual collisional rate coefficients is better than a factor of two.

No information is available from the literature for pure rotational transitions induced by collisions of CO₂ with other molecules. We therefore adopt the CO rotational collisional rate coefficients from the LAMDA database (Schöier et al. 2005; Yang et al. 2010; Neufeld 2012). Due to the lack of dipole moment, the critical density for rotational transitions of CO₂ is expected to be very low (n_crit < 10⁴) cm^-3 and thus the exact collisional rate coefficients are not important for the higher density environments considered here. A method similar to Faure & Josselin (2008), Thi et al. (2013), Bruderer et al. (2015) is used to create the full state-to-state collisional rate coefficient matrix. The method is described in Appendix A.

2.4. CO₂ spectra

Figure 2 presents a slab model spectrum of CO₂ computed using the RADEX programme (van der Tak et al. 2007). A density of 10¹⁶ cm^-3 was used to ensure close to LTE populations of all levels. A column density of 10¹⁶ cm^-2 was adopted, close to the observed value derived by Salyk et al. (2011b), with a temperature of 750 K and linewidth of 1 km s^-1. The transitions are labelled at the approximate location of their Q-branch. The spectrum shows that, due to the Fermi splitting of the bending and stretching modes, the 15 μm feature is very broad stretching from slightly shorter than 12 μm to slightly longer than 20 μm for the absorption in the Earth atmosphere. For astronomical sources, the lines between 14 and 16 μm are more realistic targets.

Two main emission features are seen in the spectrum. The strong feature around 4.3 μm is caused by the radiative decay of the 00⁰1 vibration level to the vibrational ground state. As a Σ−Σ transition this feature does not have a Q-branch, but the R and P branches are the brightest features in the spectrum in LTE at 750 K. The second strong feature is at 15 μm. This emission is caused by the radiative decay of the 01¹0 vibrational state into the ground state. It also contains small contribution by the 02²0 → 01¹0 and 03³0 → 02²0 transitions. This feature does have a Q-branch that has been observed both in absorption (Lahuis et al. 2006) and emission (Carr & Najita 2008; Pontoppidan et al. 2010). The CO₂Q-branch is found to be narrow compared to the other Q-branches of HCN and C₂H₂ measured in the same sources.

The narrowness is partly due to the fact that the CO₂Q-branch is intrinsically narrower than the same feature for HCN. This is connected with the change in the rotational constant between the ground and excited vibrational states. A comparison between Q-branch profiles for CO₂ and HCN for two optically thin LTE models is presented in Fig. 3. The lighter HCN has a full width half maximum (FWHM) that is about 50% larger than that of CO₂. The difference in the observed width of the feature is generally larger (Salyk et al. 2011b): the HCN feature is typically twice as wide as the CO₂ feature. Thus the inferred temperature from the CO₂Q-branch from the observations is low compared to the temperature inferred from the HCN feature. The difference is amplified by the intrinsically narrower CO₂Q-branch, making it more striking.

2.5. Dependence on kinetic temperature, density and radiation field

The excitation of, and the line emission from, a molecule depend strongly on the environment of the molecule, especially the kinetic temperature, radiation field and collisional partner density. In Fig. 4 slab model spectra of CO₂ for different physical parameters are compared. The dependence on the radiation field is modelled by including a blackbody field of 750 K diluted with a factor W: ⟨ J_ν ⟩ = WB_ν(T_rad) with T_rad = 750 K. When testing the effects of the kinetic temperature and density, no incident radiation field is included (W = 0).

Figure 4 shows that at a constant density of 10¹² cm^-3 the 4.3 μm band is orders of magnitude weaker than the 15 μm band. The 15 μm band increases in strength and also in width, with increasing temperature as higher J levels of the CO₂v₂ vibrational mode can be collisionally excited. Especially the spectrum at 1000 K shows additional Q branches from transitions originating from the higher energy 10⁰0(1) and 10⁰0(2) vibrational levels at 14 and 16 μm.

In the absence of a pumping radiation field, collisions are needed to populate the higher energy levels. With enough collisions, the excitation temperature becomes equal to the kinetic temperature. The density at which the excitation temperature of a level reaches the kinetic temperature depends on the critical density: n_c = A_ul/K_ul for a two-level system, where A_ul is the Einstein A coefficient from level u to level l and K_ul is the collisional rate coefficient between these levels. For densities below the critical density the radiative decay is much faster than the collisional excitation and de-excitation. This means that the line intensity scales as n/n_c. Above the critical density collisional excitation and de-excitation are fast: the intensity is then no longer dependent on the density. The critical density of the 15 μm band is close to 10¹² cm^-3, so there is little change in this band when increasing the density above this value. However, when decreasing the density below the critical value this results in the a strong reduction of the band strength. The critical density of the 4.3 μm feature is close to 10¹⁵ cm^-3 so below this the lines are orders of magnitude weaker than would be expected from LTE.

Adding a radiation field has a significant impact on both the 4.3 and 15 μm features. The radiation of a black body of 750 K peaks around 3.8 μm so the 4.3 μm/15 μm flux ratio in these cases is larger than the flux ratio without radiation field for densities below the critical density of the 4.3 μm lines. Another difference between the collisionally excited and radiatively excited states is that in the latter case vibrational levels that cannot be directly excited from the ground state by photons, such as the 10⁰0(1) and 10⁰0(2) levels, are barely populated at all.

3. CO₂ emission from a protoplanetary disk

To properly probe the chemistry in the inner disk from infrared line emission one needs to go beyond slab models with their inherent degeneracies. A protoplanetary disk model such as that used here includes more realistic geometries and contains a broad range of physical conditions constrained by observational data. Information can be gained on the location and extent of the emitting CO₂ region as well as the nature of the excitation process. By comparing with observational data, molecular abundances can be inferred as function of location. A critical aspect of the models is the infrared continuum radiation field, which has to be calculated accurately throughout the disk. This means that detailed wavelength dependent dust opacities need to be included and dust temperatures have to be calculated on a very fine grid, since the pumping radiation can originate in a different part of the disk than the lines, for example, the near-infrared for close to the inner rim. The dust is also important in absorbing some of the line flux, effectively hiding parts of the disk from our view.

In this section, the CO₂ spectra are modelled using the DALI (Dust and Lines) code (Bruderer et al. 2012; Bruderer 2013). The focus is on emission from the 15 μm lines that have been observed with Spitzer and will be observable with JWST-MIRI. Trends in the shape of the v₂Q-branch and the ratios of lines in the P- and R-branches are investigated and predictions are presented. First the model and its parameters are introduced and the results of one particular model are used as illustration. Finally the effects of various parameters on the resulting line fluxes are shown, in particular source luminosity and gas/dust ratio. As in Bruderer et al. (2015), the model is based on the source AS 205 (N) but should be representative of a typical T Tauri disk.

3.1. Model setup

Details of the full DALI model and benchmark tests are reported in Bruderer et al. (2012) and Bruderer (2013). Here we use the same parts of DALI as in Bruderer et al. (2015). The model starts with the input of a dust and gas surface density structure. The gas and dust structures are parametrized with a surface density profile (1)and vertical distribution (2)with the scale height angle h(R) = h_c(R/R_c)^ψ. The values of the parameters for the AS 205 (N) disk are taken from Andrews et al. (2009) who fitted both the SED and submillimeter images simultaneously. As the inferred structure of the disk is strongly dependent on the dust opacities and size distribution, the same values from Andrews et al. (2009) are used. They are summarized in Table 2 and the gas density structure is shown in Fig. 5, panel a. The central star is a T Tauri star with excess UV due to accretion. All the accretion luminosity is assumed to be released at the stellar surface as a 10⁴ K blackbody. The density and temperature profile are typical for a strongly flared disk as used here. The temperature, radiation field and CO₂ excitation structure can be found in the Appendix, Fig. C.1.

In setting up the model special care was taken at the inner rim, where optical and UV photons are absorbed by the dust over a very short physical path. To properly get the temperature structure of the disk directly after the inner rim, high resolution in the radial direction is needed. Varying the radial width of the first cells showed that the temperature structure only converges when the cell width of the first handful of cells is smaller than the mean free path of the UV photons.

The model dust structure is irradiated by the star and the interstellar radiation field. A Monte-Carlo radiative transfer module calculates the dust temperature and the local radiation field at all positions throughout the disk. The gas temperature is then assumed to be equal to the dust temperature. This is not true for the upper and outer parts of the disk. For the regions were CO₂ is abundant in our models the difference between dust temperature and gas temperature computed by self-consistently calculating the chemistry and cooling is less than 5%. The excitation module calculates the CO₂ level populations, using a 1+1D escape probablity that includes the continuum radiation due to the dust (Appendix A.2 in Bruderer 2013). Finally the synthetic spectra are derived using the ray tracing module, which solves the radiative transfer equation along rays through the disk. The ray tracing module as presented in Bruderer et al. (2012) is used as well as a newly developed ray-tracing module that is presented in Appendix B which is orders of magnitude faster, but a few percent less accurate. In the ray-tracing module a thermal broadening and turbulent broadening with FWHM ~ 0.2 km s^-1 is used, which means that thermal broadening dominates above ~ 40 K. The gas is in Keplerian rotation around the star. This approach is similar to Meijerink et al. (2009) and Thi et al. (2013) for H₂O and CO respectively. However Thi et al. (2013) used a chemical network to determine the abundances, whereas here only parametric abundance structures are used to avoid the added complexity and uncertainties of the chemical network.

The adopted CO₂ abundance is either a constant abundance or a jump abundance profile. The abundance throughout the paper is defined as the fractional abundance w.r.t n_H = n(H) + 2n(H₂). The inner region is defined by T > 200 K and A_V > 2 mag, which is approximately the region where the transformation of OH into H₂O is faster than the reaction of OH with CO to form CO₂. The outer region is the region of the disk with T < 200 K or A_V < 2 mag, where the CO₂ abundance is expected to peak. No CO₂ is assumed to be present in regions with A_V < 0.5 mag as photodissociation is expected to be very efficient in this region. The physical extent of these regions is shown in panel b of Fig. 5.

As shown by Meijerink et al. (2009) and Bruderer et al. (2015), the gas-to-dust (“G/D”) ratio is very important for the resulting line fluxes as the dust photosphere can hide a large portion of the potentially emitting CO₂. Here the gas-to-dust ratio is changed in two ways, by increasing the amount of gas, or by decreasing the amount of dust. When the gas mass is increased and thus the dust mass kept at the standard value of 2.9 × 10^-4 M_⊙, this is denoted by g/d_gas. If the dust mass is decreased and the gas mass kept at 0.029 M_⊙ this is denoted by g/d_dust.

3.2. Model results

Panel c of Fig. 5 presents the contribution function for one of the 15 μm lines, the v₂1 → 0Q(6) line. The contribution function shows the relative, azimuthally integrated contribution to the total integrated line flux. Contours show the areas in which 25% and 75% of the emission is located. Panel c also includes the τ = 1 surface for the continuum due to the dust, the τ = 1 surface for the v₂1 → 0Q(6) line and surface where the density is equal to the critical density. The area of the disk contributing significantly to the emission is large, an annulus from approximately 0.7 to 30 AU. The dust temperature in the CO₂ emitting region is between 100 and 500 K and the CO₂ excitation temperature ranges from 100–300 K (see Fig. C.1). The density is lower than the critical density at any point in the emitting area.

Panel d of Fig. 5 shows the contribution for the v₃ 1 → 0 R(7) line with the same lines and contours as panel c. The critical density for this line is very high, ~ 10¹⁵ cm^-3. This means that except for the inner 1 AU near the mid-plane, the level population of the v₃ level is dominated by the interaction of the molecule with the surrounding radiation field. The emitting area of the v₃1 → 0 R(7) line is smaller compared to that of the line at 15 μm. The emitting area stretches from close the the sublimation radius up to ~ 10 AU. The excitation temperatures for this line are also higher, ranging from 300–1000 K in the emitting region (see Fig. C.1).

In Fig. 6 the total flux for the 00⁰1−00⁰0 R(7) line at 4.25 μm and the 15 μm feature integrated from 14.8 to 15.0 μm are presented as functions of x_out, for different gas-to-dust ratios and different x_in. The 15 μm flux shows an increase in flux for increasing total CO₂ abundance and gas-to-dust ratio and so does the line flux of the 4.25 μm line for most of the parameter space. The total flux never scales linearly with abundance, due to different opacity effects. The dust is optically thick at infrared wavelengths up to 100 AU, so there will always be a reservoir of gas that will be hidden by the dust. The lines themselves are strong (have large Einstein A coefficients) and the natural line width is relatively small (0.2 km s^-1 FWHM). As a result the line centres of transitions with low J values quickly become optically thick. Therefore, if the abundance, and thus the column, in the upper layers of the disk is high, the line no longer probes the inner regions. This can be seen in Fig. 6 as the fluxes for models with different x_in converge with increasing x_out. Convergence happens at lower x_out for higher gas-to-dust ratios. The inner region is quickly invisible through the 4.25 μm line with increasing gas-to-dust ratios: for a gas-to-dust ratio of 10 000, there is a less than 50% difference in fluxes between the models with different inner abundances, even for the lowest outer abundances. This is not seen so strongly in the 15 μm feature as it also includes high J lines which are stronger in the hotter inner regions and are not as optically thick as the low J lines. There is no significant dependence of the flux on the inner abundance of CO₂ if the outer abundance is >3 × 10^-7 and the gas to dust ratio is higher than 1000. In these models the 15 μm feature traces part of the inner 1 AU but only the upper layers.

Different ways of modelling the gas-to-dust ratio has little effect on the resulting fluxes. Figure 6 shows the fluxes for a constant dust mass and increasing gas mass for increasing the gas-to-dust ratio, whereas Fig. D.1 in Appendix D shows the fluxes for decreasing dust mass for a constant gas mass. The differences in fluxes are very small for models with the same gas/dust ratio times CO₂ abundance, irrespective of the total gas mass: fluxes agree within 10% for most of the models. This reflects the fact that the underlying emitting columns of CO₂ are similar above the dust τ = 1 surface. Only the temperature of the emitting gas changes: higher temperatures for gas that is emitting higher up in a high gas mass disk and lower temperatures for gas that is emitting deeper into the disk in a low dust mass disk.

The grey band in Figs. 6 and D.1 shows the range of fluxes observed for protoplanetary disks scaled to a common distance of 125 pc (Salyk et al. 2011b). This figure immediately shows that low CO₂ abundances, x_out< 3 × 10^-7, are needed to be consistent with the observations. Some disks have lower fluxes than given by the lowest abundance model, which can be due to other parameters. A more complete comparison between model and observations is made in Sect. 4.1.

In Appendix E a comparison is made between the fluxes of models with CO₂ in LTE and models for which the excitation of CO₂ is calculated from the rate coefficients and the Einstein A coefficients. The line fluxes differ by a factor of about three between the models, similar to the differences found by Bruderer et al. (2015, their Fig. 6) for the case of HCN.

3.2.1. The v₂ band emission profile

Figure 7 shows the v₂Q-branch profile at 15 μm for a variety of models. All lines have been convolved to the resolving power of JWST-MIRI at that wavelength (R = 2200, Rieke et al. 2015; Wells et al. 2015) with three bins per resolution element. Panel a shows the results from a simple LTE slab model at different temperatures whereas panels b and c presents the same feature from the DALI models. Panel b contains models with different gas-to-dust ratios and abundances (assuming x_in = x_out) scaled so g/d × x_CO2 is constant. It shows that gas-to-dust ratio and abundance are degenerate. It is expected that these models show similar spectra, as the total amount of CO₂ above the dust photosphere is equal for all models. The lack of any significant difference shows that collisional excitation of the vibrationally excited state is insignificant compared to radiative pumping. Panel c of Fig. 7 shows the effect of different inner abundances on the profile. For the highest inner abundance shown, 1 × 10^-6, an increase in the shorter wavelength flux can be seen, but the differences are far smaller than the differences between the LTE models. Panel d shows models with similar abundances, but with increasing g/d_dust. The flux in the 15 μm feature increases with g/d_dust for these models as can be seen in Fig. D.1. This is partly due to the widening of the feature as can be seen in Panel d which is caused by the removal of dust. Due to the lower dust photosphere it is now possible for a larger part of the inner region to contribute to this emission. The inner region is hotter and thus emits more towards high J lines causing the Q-branch to widen.

Fitting of LTE models to DALI model spectra in Figs. 7b–d results in inferred temperatures of 300–600 K. Only the models with a strong tail (blue lines in 7b and 7d) need temperatures of 600 K for a good fit, the other models are well represented with ~300 K. For comparison, the actual temperature in the emitting layers is 150–350 K (Fig. C.1), illustrating that the optically thin model overestimates the inferred temperatures. The proper inclusion of optical depth effects for the lower-J lines lowers the inferred temperatures. This means that care has to be taken when interpreting a temperature from the CO₂ profile. A wide feature can be due to high optical depths or high rotational temperature of the gas.

A broader look at the CO₂ spectrum is thus needed. The left panel of Fig. 8 shows the P, Q and R-branches of the vibrational bending mode transition at R = 2200, for models with different inner CO₂ abundances and the same outer abundance of 10^-7. The shape for the R- and P-branches is flatter for low to mid-J and slightly more extended at high J in the spectrum from the model with an inner CO₂ abundance of 10^-6 than the other spectra. The peaks at 14.4 μm and 15.6 μm are due to the Q-branches from the transitions between 11¹0(1) → 10⁰1(1) and 11¹0(2) → 10⁰1(2) respectively. These are overlapping with lines from the bending fundamental P and R branches. For the constant and low inner CO₂ abundances, 10^-7 and 10^-8 respectively R- and P-branch shapes are similar, with models differing only in absolute flux. Decreasing the inner CO₂ abundance from 10^-8 to lower values has no effect of the line strengths.

The right panel of Fig. 8 shows Boltzmann plots of the spectra on the left. The number of molecules in the upper state inferred from the flux is given as a function of the upper state energy. The number of molecules in the upper state is given by: , with d the distance to the object, F the integrated line flux, g_u the statistical weight of the upper level and A_ul and ν_ul the Einstein A coefficient and the frequency of the transition. From slope of vs. E_up a rotational temperature can be determined. The expected slopes for 400, 600 and 800 K are given in the figure. It can be seen that the models do not show strong differences below J = 20, where emission is dominated by optically thick lines. Towards higher J, the model with x_in = 10^-6 starts to differ more and more from the other two models. The models with x_in = 10^-7 and x_in = 10^-8 stay within a factor of two of each other up to J = 80 where the molecule model ends.

Models with similar absolute abundances of CO₂ (constant g/d × x_CO2) but different g/d_gas ratios are nearly identical: the width of the Q-branch and the shapes of the P- and R-branches are set by the gas temperature structure. This temperature structure is the same for models with different g/d_gas ratios as it is set by the dust structure. The temperature is, however, a function of g/d_dust, but those temperature differences are not large enough for measurable effects. From this it also follows that the exact collisional rate coefficients are not important: the density is low enough that the radiation field can set the excitation of the vibrational levels. At the same time the density is still high enough to be higher than the critical density for the rotational transitions, setting the rotational excitation temperature equal to the gas kinetic temperature.

The branch shapes are a function of g/d_dust at constant absolute abundance. Apart from the total flux which is slightly higher at higher g/d_dust (Fig. D.1), the spectra are also broader (Panel d, Fig. 7). This is because the hotter inner regions are less occulted by dust for higher g/d_dust ratios. This hotter gas has more emission coming from high J lines, boosting the tail of the Q-branch.

To quantify the effects of different abundance profiles, line ratios can also be informative. The lines are chosen so they are free from water emission (see Appendix F). The top two panels of Fig. 9 shows the line ratios for lines in the 01¹0(1) → 00⁰0(1) 15 μm band: R(37):R(7) and P(15):P(51). The R(7) and P(15) lines come from levels with energies close to the lowest energy level in the vibrational state (energy difference is less than 140 K). These levels are thus easily populated and the lines coming from these levels are quickly optically thick. The R(37) and P(51) lines come from levels with rotational energies at least 750 K above the ground vibrational energy. These lines need high kinetic temperatures to show up strongly and need higher columns of CO₂ at prevailing temperatures to become optically thick. From observation of Fig. 9 a few things become clear. First, for very high outer abundances, it is very difficult to distinguish between different inner abundances based on the presented line ratio. Second, models with high outer abundances are nearly degenerate with models that have a low outer abundance and a high inner abundance. A measure of the optical depth will solve this. In the more intermediate regimes the line ratios presented here or a Boltzmann plot will supplement the information needed to distinguish between a cold, optically thick CO₂ reservoir and a hot, more optically thin CO₂ reservoir that would be degenerate in just Q-branch fitting.

3.2.2. ¹³CO₂v₂ band

An easier method to break these degeneracies is to use the ¹³CO₂ isotopologue. ¹³CO₂ is approximately 68 times less abundant compared to ¹²CO₂, using a standard local interstellar medium value (Wilson & Rood 1994; Milam et al. 2005). This means that the isotopologue is much less likely to be optically thick and thus ¹³CO₂:¹²CO₂ line ratios can be used as a measure of the optical depth, adding the needed information to lift the degeneracies. The bottom panel of Fig. 9 shows the ratio between the flux in the ¹³CO₂v₂Q-branch and the ¹²CO₂v₂P(25) line.

As the Q-branch for ¹³CO₂ is less optically thick, it is also more sensitive to the abundance structure. The Q-branch, situated at 15.42 μm, partially overlaps with the P(23) line of the more abundant isotopologue so both isotopologues need to be modelled to properly account the the contribution of these lines. Figure 10 shows the same models as in Fig. 8 but now with the ¹³CO₂ emission in thick lines. The ¹³CO₂Q-branch is predicted to be approximately as strong as the nearby ¹²CO₂ lines for the highest inner abundances. The total flux in the ¹³CO₂Q-branch shows a stronger dependence on the inner CO₂ abundance than the ¹²CO₂Q-branch. A hot reservoir of CO₂ strongly shows up as an extended tail of the ¹³CO₂Q-branch between 15.38 and 15.40 μm.

3.2.3. Emission from the v₃ band

The v₃ band around 4.25 μm is a strong emission band in the disk models, containing a larger total flux than the v₂ band. Even so, the 4.3 μm band of gaseous CO₂ has not been seen in observations of ISO with the Short Wave Spectrometer (SWS) towards high mass protostars in contrast with 15 μm band that has been seen towards these sources in absorption (van Dishoeck et al. 1996; Boonman et al. 2003a). This may be largely due to the strong solid CO₂ 4.2 μm ice feature obscuring the gas-phase lines for the case of protostars, but for disks this should not be a limitation. Figure 11 shows the spectrum of gaseous CO₂ in the v₃ band around 4.3 μm at JWST-NIRSpec resolving power. The resolving power of NIRSpec is taken to be R = 3000, which is not enough to fully separate the lines from each other. The CO₂ emission thus shows up as an extended band.

The band shapes in Fig. 11 are very similar. The largest difference is the strength of the 4.2 μm discontinuity, which is probably an artefact of the model as only a finite number of J levels are taken into account. The total flux over the whole feature does depend on the inner abundance, but the difference is of the order of ~ 10% for 2 orders of magnitude change of the inner abundance.

Figure 11 also shows the ¹³CO₂ spectrum. The lines from ¹³CO₂ are mostly blended with much stronger lines from ¹²CO₂. At the longer wavelength limit, ¹³CO₂ lines are stronger than those of ¹²CO₂ but there the 6 μm water band and 4.7 μm CO band start to complicate the detection of ¹³CO₂ in the 4–5 μm region.

Average abundances of CO₂ can be derived from observations of the 4.3 μm band. While inferring the abundance structure will be easier from the 15 μm band there are some observational advantages of using the 4.3 μm band. NIRSpec has multi-object capabilities and will thus be able to get large samples of disks in a single exposure, especially for more distant clusters where there are many sources in a single FOV. NIRSpec has the additional advantages that it does not suffer from detector fringing and that it is more sensitive (NIRSpec pocket guide², Rieke et al. 2015; Wells et al. 2015). As both the 4.3 μm and 15 μm bands are pumped by infrared radiation, the flux ratios between these two will mostly contain information about the ratio of the continuum radiation field between the wavelengths of these bands.

3.3. Line-to-continuum ratio

The line-to-continuum ratio is potentially an even better diagnostic of the gas/dust ratio than line ratios (Meijerink et al. 2009). In Fig. 12 we present spectra with the continuum added to it. The spectra have been shifted, as the continua for these models overlap. The models for which the spectra are derived all have a CO₂ abundance of 10^-8 but differ in the gas-to-dust ratio. The gas-to-dust ratio, determines the column of CO₂ that can contribute to the line. A large part of the CO₂ reservoir near the mid-plane cannot contribute due to the large continuum optical depth of the dust. It is thus not surprising that the line-to-continuum ratio is strongly dependant on the gas-to-dust ratio. The precise method for setting the gas-to-dust ratio (by increasing the amount of gas, or decreasing the amount of dust) does not really matter for the line-to-continuum ratio. It does matter for the absolute scaling of the continuum, which decreases if the amount of dust is decreased.

Meijerink et al. (2009) could constrain the gas-to-dust ratio from the data since there is an upper limit to the H₂O abundance from the atomic O abundance. Here it is not possible to make a similar statement as CO₂ is not expected to be a major reservoir of either the oxygen or the carbon in the disk. On the contrary, from Fig. 6 it can be seen that with a gas-to-dust ratio of 100 an abundance of 10^-7 is high enough to explain the brightest of the observed line fluxes. The line-to-continuum ratios from the Spitzer-IRS spectra of 5–10% are also matched by the same models (see Fig. 13). External information such as can be obtained from H₂O is needed to lift the degeneracy between high gas-to-dust ratio and high abundance: if one of the two is fixed, the other can be determined from the flux or line-to-continuum ratio.

The line-to-continuum ratio is very important for planning observations, however, as it sets the limit on how precisely the continuum needs to be measured to be able to make a robust line detection. Figure 13 shows the line-to-continuum ratios for models with a constant abundance. These figures show that high signal-to-noise (S/N) on the continuum is needed to be able to get robust line detections. The ¹²CO₂Q-branch should be easily accessible for most protoplanetary disks. To be able to probe individual P- and R-branch lines of the ¹²CO₂ 15 μm feature as well as the ¹³CO₂Q-branch, deeper observations (reaching S/N of at least 300, up to 1000) will be needed to probe down to disks with CO₂ abundances of 10^-8 and gas-to-dust ratios of 1000.

3.4. CO₂ from the ground

As noted earlier, there is a large part of the CO₂ spectrum that cannot be seen from the ground because of atmospheric CO₂. There are a few lines, however, that could be targeted from the ground using high spectral resolution. The high J lines (J> 70) of the v₁ = 1−0 transition in the R branch around 4.18 μm are visible with a resolving power of R = 30 000 or higher. At this resolution the CO₂ atmospheric lines are resolved and at J> 70 they are narrow enough to leave 20−50% transmission windows between them (ESO skycalc³). The lines are expected to have a peak line-to-continuum ratio of 1:100. So a S/N of 10 on the line peak translates to a S/N of 1000 on the continuum. The FWHM of the atmospheric lines is about 30 km s^-1. So half of the emission line profile should be observable when the relative velocity shift between observer and source is more than 15 km s^-1. Since the Earth’s orbit allows for velocity shifts up to 30 km s^-1 in both directions, observing the full line profile is possible in two observations at different times of year for sources close to the orbital plane of the earth. The exposure time needed to get a S/N of 1000 on the continuum at 4.18 μm, which is 6.7 Jy in our AS 205 (N) model, with a 40% sky transmission on the lines, is about ten hours for VLT-CRIRES⁴. The high JP-branch lines of CO₂ are close to atmospheric lines from N₂O, O₃ and H₂O resulting in a very opaque atmosphere at these wavelengths (Noll et al. 2012; Jones et al. 2013).

The other lines that can be seen from the ground are between 9 and 12 μm (in the N band). These originate from the 01¹1 and 00⁰1 levels. The line to continuum ratios vary from 1:40 to 1:3000 for the brighter lines in the 00⁰1 → 10⁰0(1) band with the most likely models having line to continuum ratios between 1:200 and 1:2000. For a continuum of ~11 Jy a S/N of 2000 for R = 10⁵ could potentially be achieved in about three to ten minutes of integration with the European Extremely Large Telescope (E-ELT). At the location of the CO₂ atmospheric absorption lines in this part of the spectrum the sky transmittance is 50% and the atmospheric lines have a FWHM of ~50 km s^-1

3.5. CO₂ model uncertainties

The fluxes derived from the DALI models depend on the details of the CO₂ excitation processes included in the model. The collisional rate coefficients are particularly uncertain, since the measured set is incomplete. There are multiple ways to extrapolate what is measured to what is needed to complete the model. Modelling slabs of CO₂ using different extrapolations such as: absolute scaling of the rotational collision rate coefficients, including temperature dependence of the vibrational collisional rate coefficients and different implementations of the collision rate coefficients between the vibrational levels with 2v₁ + v₂ = constant, show that fluxes can change by up to 50% for specific combinations of radiation fields and densities. The highest differences are seen in the 4–5 μm band, usually at low densities. The flux in the 15 μm band usually stays within 10% of the flux of the model used here. These uncertainties are small compared to other uncertainties in disk modelling such as the chemistry or parameters of the disk hosting protostar. The main reason that the fluxes are relatively insensitive to the details of the collisional rate coefficients is due to the importance of radiative pumping in parallel with collisions.

The assumption T_dust = T_gas is not entirely correct, since the gas temperature can be up to 5% higher than the dust temperature in the CO₂ emitting regions. This affects our line fluxes. For the 4.3 μm fluxes the induced difference in flux is always smaller than 10%. For the 15 μm fluxes difference are generally smaller than 10%, whereas some of the higher J lines are up to 25% brighter.

A very simplified abundance structure was taken. It is likely that protoplanetary disks will not have the abundance structure adopted here. Full chemical models indeed show much more complex chemical structures (see e.g. Walsh et al. 2015). The analysis carried out in this work here should still be appropriate for more complex abundance structures, and future work will couple such chemistry models directly with the excitation and radiative transfer.

The stellar parameters for the central star and the exact parameters of the protoplanetary disk also influence the resulting CO₂ spectrum. The central star influences the line emission through its UV radiation that can both dissociate molecules and heat the gas. Since for our models no chemistry is included, only the heating of the dust by stellar radiation is important for our models. The CO₂ flux in the emission band around 15 μm scales almost linearly with the bolometric luminosity of the central object (see Fig. 14).

4. Discussion

4.1. Observed 15 μm profiles and inferred abundances

The v₂ 15 μm feature of CO₂ has been observed in many sources with Spitzer-IRS (Pontoppidan et al. 2010; Salyk et al. 2011b). The SH (Short-High) mode barely resolves the 15 μm Q-branch, but that is enough to compare with the models. We used the spectra that have been reduced with the Caltech High-res IRS pipeline (CHIP; Pontoppidan et al. 2010; Pontoppidan 2016). The sources selected out of the repository have a strong emission feature of CO₂ but no distinguishable H₂O emission in the 10–20 μm range. The sources and some stellar parameters are listed in Table 3. The observed spectra are continuum subtracted (Appendix G) and the observed profiles are compared with model profiles by eye (Fig. 15).

Two sets of comparisons are made. For the first set, the model fluxes are only corrected for the distance to the objects. For the other set, the model fluxes are scaled for the distance but also scaled for the luminosity of the central source, using L_CO2 ∝ L_⋆. This relation is found by running a set of models with a range of luminosities, presented in Fig. 14. Aside from the luminosity of the star, all other parameters have been kept the same including the shape of the stellar spectrum. The effective temperature of the star mostly affects the fraction of short wavelength UV photons which can photodissociate molecules, but since no detailed chemistry is included, the use of a different stellar temperature would not change our results. Other tests (not shown here) have indeed shown that the shape of the spectrum does not really matter for the CO₂ line fluxes in this parametric model. All models have a gas-to-dust ratio of 1000 and a constant CO₂ abundance of 10^-8.

Both the total flux in the range between 14.7 and 15.0 μm and that of a single line in this region (the 01¹0Q(6) line) have an almost linear relation with luminosity of the central star. For the 00⁰1R(7) line around 4.3 μm the dependence on the central luminosity is slightly more complex. Below a stellar luminosity of 1 L_⊙ the dependence is stronger than linear, but above that the dependence becomes weaker than linear. Overall, it is reasonable to correct the 15 μm fluxes from our model for source stellar luminosities using the linear relationship. This is because the amount of infrared continuum radiation that the disk produces scales linearly with the amount of energy that is put into the disk by the stellar radiation. It is the infrared continuum radiation that sets the molecular emission the due to radiative pumping, the dominant vibrational excitation mechanism for CO₂.

The model spectra are overplotted on the continuum subtracted observations in Fig. 15. The flux in these models has been scaled with the distance of the source and the luminosity of the central star. A gas-to-dust ratio of 1000 is adopted as inferred from H₂O observations (Meijerink et al. 2009).

An overview of the inferred abundances is given in Table 4. For the DALI models the emitting CO₂ column and the number of emitting CO₂ molecules have been tabulated in Table 4. The column is defined as the column of CO₂ above the τ_dust = 1 line at the radial location of the peak of the contribution function (Fig. 5d). The number of molecules is taken over the region that is responsible for half the total emission as given by the contribution function. The number of molecules shown is thus the minimum amount of CO₂ needed to explain the majority of the flux and a sets lower limit for the amount of CO₂ needed to explain all of the emission. The inferred abundances range from 10^-9−10^-7. They agree with that inferred by Pontoppidan & Blevins (2014) using an LTE disk model appropriate for the RNO 90 disk, demonstrating that non-LTE excitation effects are minor (see also Appendix E).

For GW Lup and SZ 50 the emitting CO₂ columns found by Salyk et al. (2011b) are within a factor of two from those inferred here, whereas for DN Tau and IM Lup our inferred columns are consistent with the upper limits from the slab models (tabulated in Table 4). However, the inferred column for HD 101412 differs greatly. For all disks the number of molecules in our models is at least an order of magnitude higher than the number of molecules inferred from the LTE models. The emitting area used by Salyk et al. (2011b) in fitting the CO₂ feature was fixed and generally taken to be slightly larger than the inner 1 AU. This is very small compared to the emitting area found in this work which extends up to 30 AU. It is thus unsurprising that the total number of CO₂ molecules inferred is lower for the LTE slab models from Salyk et al. (2011b). The high number of molecules needed for the emission in our models is also related to the difference in excitation: the vibrational excitation temperature of the gas in the non-LTE models is lower (100–300 K) than the temperature fitted for the LTE models (~ 650 K). Thus in the non-LTE models a larger number of molecules is needed to get the same total flux. The narrow CO₂ profile is due to low rotational temperatures as emission from large radii >2–10 AU dominates the strongest lines. The visual contrast is enhanced by the fact that the CO₂ feature is also intrinsically narrower at similar temperatures than that of HCN (Fig. 3). For HD 101412 the model feature is notably narrower than the observed feature signifying either a higher CO₂ rotational temperature, or a more optically thick emitting region.

There are of course caveats in the comparisons done here. The standard model uses a T Tauri star that is luminous (total luminosity of 7.3 L_⊙) and that disk is known to have very strong H₂O emission. The sample of comparison protostars consists of 7 T Tauri stars that have luminosities a factor of 2–35 lower than assumed in our model and a Herbig Ae star that is more than three times as luminous. A simple correction for source luminosity is only an approximation. All of these sources have little to no emission lines of H₂O in the mid-infrared. This may be an indication of different disk structures, and the disk model used in this work may not be representative of these water-rich disks. Indeed, Banzatti & Pontoppidan (2015) found that the emitting radius of the CO ro-vibrational lines scales inversely with the vibrational temperature inferred from the CO emission. This relation is consistent with inside-out gap opening. Comparing the CO ro-vibrational data with H₂O infrared emission data from VLT-CRIRES and Spitzer-IRS Banzatti et al. (2017) found that there is a correlation between the radius of the CO ro-vibrational emission and the strength of the water emission lines: the larger the radius of the CO emission, the weaker the H₂O emission. This suggests that the H₂O-poor sources may also have inner gaps, where both CO and H₂O are depleted. There are only two H₂O-poor sources in our sample that overlap with Banzatti et al. (2017).

However, if our analysis is applied to sources that do have water emission, the range of best-fit CO₂ abundances is found to be similar. Figure G.2 shows CO₂ model spectra compared to observations for a set of the strongest water emitting sources. The conclusion that the abundance of CO₂ in protoplanetary disks is around 10^-8 is, therefore, robust.

The inferred CO₂ abundances are low, much lower than the expected ISM value of 10^-5 if all of the CO₂ would result from sublimated ices. This demonstrates that the abundances have been reset by high temperature chemistry, as also concluded by Pontoppidan & Blevins (2014).

The inferred low abundances agree well with chemical models by Walsh et al. (2015). However, the column found for chemical models by Agúndez et al. (2008), ~ 6 × 10¹⁶ cm^-2, is more than an order of magnitude higher. Agúndez et al. (2008) used a different lower vertical bound for their column integration and only considered the inner 3 AU. Either of these assumptions may explain the difference in the CO₂ column.

4.2. Tracing the CO₂ iceline

One of the new big paradigms in (giant) planet formation is pebble accretion. Pebbles, in models defined as dust particles with a Stokes number around one, are badly coupled to the gas, but generally not massive enough to ignore the interaction with the gas. This means that these particles settle to the mid-plane and radially drift inward on short time scales. This pebble flux allows in theory a planetesimal to accrete all the pebbles that form at radii larger than its current location (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Levison et al. 2015).

This flux of pebbles also has consequences for the chemical composition of the disk. These pebbles should at some point encounter an iceline, if they are not stopped before. At the iceline they should release the corresponding volatiles. The same holds for any drifting planetesimals (Ciesla & Cuzzi 2006). As the ice composition is very different from the gas composition, this can in principle strongly change the gas content in a narrow region around the ice line. For this effect to become observable in mid-infrared lines, the sublimated ices should also be mixed vertically to higher regions in the disk.

From chemical models the total gas-phase abundance of CO₂ around the CO₂ iceline is thought to be relatively low, (10^-8, Walsh et al. 2015) similar to the value found in this work. The CO₂ ice content in the outer disk can be orders of magnitude higher. Both chemical models and measurements of comets show that the CO₂ content in ices can be more than 20% of the total ice content (Le Roy et al. 2015; Eistrup et al. 2016), with CO₂ ice even becoming more abundant than H₂O ice in some models of outer disk chemistry (Drozdovskaya et al. 2016). This translates into an abundance up to a few × 10^-5. Here, we investigate whether the evaporation of these CO₂ ices around the iceline would be observable.

To model the effect of pebbles moving over the iceline, a model with a constant CO₂ abundance of 1 × 10^-8 and a gas-to-dust ratio of 1000 is taken. In addition, the abundance of CO₂ is enhanced in an annulus where the midplane temperature is between 70 K and 100 K (grey region in Fig. 5b), corresponding to the sublimation temperature of pure CO₂ ice (Harsono et al. 2015). This results in a radial area between 8 and 15 AU in our case. The abundance is taken to be enhanced over the total vertical extent of the CO₂ in the model, as in the case of strong vertical mixing. The spectra from three models with enhanced abundances of x_CO2,ring = 10^-6,10^-5 and 10^-4 in this ring can be seen in Fig. 16 together with the spectrum for the constant x_CO2 = 10^-8 model. We note that CO₂ ice is unlikely to be pure, and that some of it will likely also come off at the H₂O iceline, but such a multi-step sublimation model is not considered here.

The spectra in Fig. 16 show both ¹²CO₂ and ¹³CO₂ for the model with a constant abundance, an enhanced abundance of 1 × 10^-6, 1 × 10^-5, and 1 × 10^-4 around the CO₂ iceline. The enhanced CO₂ increases the ¹²CO₂ flux up to a factor of 3. The optically thin ¹³CO₂ feature is however increased by a lot more, with the peak flux in the ¹³CO₂Q-branch reaching fluxes that are more than two times higher than peak fluxes of the nearby ¹²CO₂P-branch lines.

Enhanced abundances in the outer regions can be distinguished from an enhanced abundance in the inner regions by looking at the tail of the ¹³CO₂Q-branch. An over-abundant inner region will show a significant flux (10–50% of the peak) from the ¹³CO₂Q-branch in the entire region between the locations of the ¹²CO₂P(21) and P(23) lines and will show a smoothly declining profile with decreasing wavelength. If the abundance enhancement is in the outer regions, where the gas temperature and thus the rotational temperature is lower, the flux between the ¹²CO₂P(21) and P(23) lines will be lower for the low enhancements (0–20% of the peak flux), for the higher enhancements the R(1) line shows up in the short wavelength side of the profile. Other R and P branch lines from ¹³CO₂ can also show up in the spectrum if the abundance can reach up to 10^-4 in the ring around the CO₂ iceline.

4.3. Comparison of CO₂ with other inner disk molecules

With the models presented in this paper, four molecules with rovibrational transitions coming from the inner disk have been studied by non-LTE disk models: H₂O (Antonellini et al. 2015, 2016; Meijerink et al. 2009), CO (Thi et al. 2013), HCN (Bruderer et al. 2015) and CO₂ (this work). Of these CO is special, as it can be excited by UV radiation and fluoresces to excited vibrational states that in turn emit infrared radiation. For the other molecules absorption of a UV photon mostly leads to dissociation of the molecule (Heays et al. 2017). H₂O, HCN and CO all have an permanent dipole moment and can thus also emit strongly in the sub-millimeter. This means that these molecules will have lower rotational temperatures than CO₂ in low density gas but they are actually observed to have broader profiles in the mid-infrared. Thus, our models reinforce the conclusion from the observed profiles that CO₂ comes from relatively cold gas (200–300 K; see panel c of Fig. 5 and panel d of Fig. C.1).

For the disk around AS 205 (N), the emission of both HCN and CO₂ has now been analysed under non-LTE conditions in DALI. As such it is possible to infer the HCN to CO₂ abundance ratio in the disk. The representative, constant abundance model for the CO₂ emission from AS 205 N has an abundance of 3 × 10^-8 with a gas-to-dust ratio of 1000 but the inner abundance can easily vary by an order of magnitude or more while still being in agreement with observations. The models from Bruderer et al. (2015) that best reproduce the data have outer HCN abundances between 10^-10 and 10^-9 for gas-to-dust ratios of 1000. This translates into CO₂/HCN abundance ratios of 30–300 in the region from 2 to 30 AU. The higher abundance of CO₂ in the outer regions of the disk explains the colder inferred rotational temperature of CO₂ compared to HCN.

5. Conclusion

Results of DALI models are presented, modelling the full continuum radiative transfer, non-LTE excitation of CO₂ in a typical protoplanetary disk model, with as the main goal: to find a way to measure the CO₂ abundance in the emitting region of disks with future instruments like JWST and test different assumptions on its origin. Spectra of CO₂ and ¹³CO₂ in the 4–4.5 μm and 14–16 μm regions were modelled for disks with different parametrized abundance structures, gas masses and dust masses. The main conclusions of this study are:

• The critical density of the CO₂00⁰1 state, responsible for emission around 4.3 μm, is very high, >10¹⁵ cm^-3. As a result, in the absence of a pumping radiation field, there is no emission from the 00⁰1 state at low densities. If there is a pumping infrared radiation field, or if the density is high enough, the emission around 4.3 μm will be brighter than that around 15 μm.

• The infrared continuum radiation excites CO₂ up to large radii (10s of AU). The region probed by the CO₂ emission can therefore be an order of magnitude larger (in radius) than typically assumed in LTE slab models. Temperatures inferred from optically thin LTE models can also be larger than the actual temperature of the emitting gas. Differences between LTE and non-LTE full disk model flux are typically a within factor of three.

• Current observations of the 15 μm Q-branch fluxes are consistent with models with constant abundances between 10^-9 and 10^-7 for a gas-to-dust ratio of 1000. Observations of lines corresponding to levels with high rotational quantum numbers or the ¹³CO₂Q-branch will have to be used to properly infer abundances. In particular, the ¹³CO₂Q-branch can be a good indicator of abundance structure from inner to outer disk.

• The gas-to-dust ratio and fractional abundance are largely degenerate. The column of CO₂ above the dust infrared photosphere sets the emission. Models with similar columns have very similar spectra irrespective of total dust and gas mass, due to the excitation mechanism of CO₂. If the gas-to-dust ratio is constrained from other observations such as H₂O the fractional abundance can be determined from the spectra.

• The abundance of CO₂ in protoplanetary disks inferred from modelling, 10^-9–10^-7, is at least two orders of magnitude lower than the CO₂ abundance in ISM ices. This implies that disk chemical abundances are not directly inherited from the ISM and that significant chemical processing happens between the giant molecular cloud stage and the protoplanetary disk stage.

• The ¹³CO₂v₂ Q-branch at 15.42 μm will be able to identify an overabundance of CO₂ in the upper layers of the inner disk, such as could be produced by sublimating pebbles and planetesimals around the iceline(s).

Our work shows that the new instruments on JWST will be able to give a wealth of information on the CO₂ abundance structure, provided that high S/N (>300 on the continuum) spectra are obtained.

---

Acknowledgments

We thank the anonymous referee for his/her suggestions that have improved the paper. 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. This work is based in part on archival data obtained made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

References

Appendix A: Collisional rate coefficients

The collisional rate coefficients are calculated in a way very similar to Bruderer et al. (2015), that is by combining vibrational coefficients with rotational rate coefficients to get the state-to-state ro-vibrational rate coefficients. Only collisions with H₂ are considered, which is the dominant gas species in the regions were CO₂ is expected to be abundant. The vibrational coefficients were taken from the laser physics and atmosperic physics papers. An overview of the final vibrational rate coefficients used are shown in Table A.1. The temperature dependence of the rates is suppressed for the de-excitation collissional rate coefficients and the rate for 300 K are used throughout. The vibrational rate coefficients are not expected to vary much over the range of temperatures considered here. The de-excitation rate coefficient of the bending mode by H₂ (01¹0 → 00⁰0 transition) from Allen et al. (1980) is 5 × 10^-12 cm^-3 s^-1, an order of magnitude faster than the He (Taylor & Bitterman 1969; Allen et al. 1980). This is probably due to vibrational-rotational energy exchange in collisions with rotationally excited H₂ (Allen et al. 1980). For levels higher up in the vibrational ladder we extrapolate the rates as in Procaccia & Levine (1975) and Chandra & Sharma (2001). Combining Eqs. (6) and (8) from the later paper we get, for v>w, the relation (A.1)The rate coefficient measured by Nevdakh et al. (2003) is actually the total quenching rate of the 00⁰1 level. Here we assume that the relaxiation of the 00⁰1 level goes to the three closest lower energy levels (03³0,11¹0(1),11¹0(2)) in equal measure. For the all the rates between the levels of the Fermi degenerate states and the corresponding bending mode with higher angular momentum the CO₂–CO₂ rate measured by Jacobs et al. (1975) was used scaled to the reduced mass of the H₂-CO₂ system. The states with constant 2v₁ + v₂ are considered equal to the pure bending mode with respect to the collisional rate coefficients to other levels.

No information exists on the rotational rate coefficients of CO₂ with H₂. We have decided to use the CO rate coefficients from Yang et al. (2010) instead. Since CO₂ does not have a dipole moment, the exact rate coefficients are not expected to be important since the critical densities for the levels in the rotational ladders are very low, <10⁴ cm^-3.

The method suggested by Faure & Josselin (2008) was employed to calculate the state-to-state de-excitation rate coefficients for initial levels v,J to all levels v′,J′ with a smaller ro-vibrational energy. This is assuming a decoupling of rotational and vibrational levels, so we can write (A.2)where (A.3)with the statistical weights g_i of the levels. All the excitation rates are calculated using the detailed balance.

Appendix B: Fast line ray tracer

For the calculation of the CO₂ lines a new ray tracer was used. The conventional ray tracer used in DALI (Bruderer et al. 2012) can take up to ten minutes to calculate the flux from one line. The CO₂ molecule model used here includes more than 3600 lines. Not all of these lines are directly important but to get the complete spectrum of both the 4.3 μm and the 15 μm bands, a few hundred lines need to be ray traced for each model.

To enable the calculation of a large number of lines a module has been implemented into DALI that can calculate a line flux in a few seconds versus a few minutes for the conventional ray tracer. The module uses the fact that, along a line of sight, the velocity shear due to the finite height of the disk is approximately linear (Horne & Marsh 1986, Eqs. (9) and (10)). Using this, the spectrum for an annulus in the disk can be approximated. At the radius of the annulus in question, the spectra are calculated for different velocities shears. These spectra are calculated by vertically integrating the equation of radiative transfer through the disk and correcting for the projected area for non face-on viewing angles. Then the total spectrum of the annulus is calculated by iterating over the azimuthal direction. For each angle the velocity shear is calculated and the spectrum is interpolated from the pre-calculated spectra. A simple sum over the spectra in all annuli is now sufficient to calculate the total spectrum.

This approximation is a powerful tool for calculating the total flux in a line especially for low inclinations. For the models presented in this paper the fluxes differ by about 4% for the 15 μm lines and 1.5% for the 4.3 μm lines. This is small compared to the other uncertainties in the models.

The approximation breaks down at high inclinations and should be used with care for any inclination larger than 45°. The total line shape is also close to the line shape from the traditional ray tracer, but with the high S/N from ALMA, using the traditional ray tracer is still advised for doing direct comparisons. This is also the case for images for which the errors will be larger than for the integrated flux or line shape as some of the errors made in making the image will cancel out (in first order) when integrating over the annulus.

Appendix C: Model temperature and radiation structure

Figure C.1 shows the model temperature, radiation field and excitation temperature structure corresponding to the model shown in Fig. 5. Panel a shows the dust temperature structure, panels b and c show the excitation temperature of the v₂ 1 → 0 Q(6) line and v₃ 1 → 0R(7) lines. For the excitation temperature only the upper and lower state of the line are used. This is thus a vibrational excitation temperature and can be different from the ground state rotational excitation temperature (that follows the dust temperature) and the rotational excitation temperature within a vibrationally excited state. Where the density is higher than the critical density the excitation temperature is equal to the dust temperature. In the disk atmosphere the Q(6) is mostly subthermally excited, while the R(7) line is superthermally excited. For both lines there is a maximum in the vertical excitation temperature distribution at the point where the gas becomes optically thick to its own radiation. Panel d shows the dust temperature of the region from which most of the CO₂ 15 μm emission originates as function of radius. Most of the emitting gas is at temperatures between 150 and 350 K. Panels e and f show the strength of the radiation field at 15 μm and 4.3 μm is shown as a factor of the radiation field of a 750 K blackbody. This shows where there is a sufficient photon density to radiatively pump CO₂.

Appendix D: Model fluxes g/d_dust

As mentioned in the main text, two different way of changing the gas-to-dust ratio were considered, increasing the gas mass and decreasing the dust mass w.r.t. the gas-to-dust ratio 100 case. Figure D.1 is the counterpart to Fig. 6 showing the modelled fluxes for different inner CO₂ abundances, outer CO₂ abundances and different gas-to-dust ratios. In this case the gas-to-dust ratio is varied by keeping the gas mass of the disk constant and decreasing the amount of dust in the disk.

There are only very slight differences between Figs. D.1 and 6, and all observations made for the figure in the main text are true for this figure as well.

Appendix E: LTE vs. non-LTE

The effects of the LTE assumption on the line fluxes in a full disk model on the v₃,1 → 0,R(7) line and the 15 μm feature are shown in Fig. E.1. Only the models with constant abundance (x_out = x_in) are shown for clarity but the differences between LTE and non-LTE for these models are representative for the complete set of models. For the 15 μm flux the differences between the LTE and non-LTE models is small, of the order of 30%. The radial extent of the emission is, however, different: the region emitting 75% of the 15 μm flux extends twice as far in the non-LTE models (extent of the 15 μm non-LTE emission is seen in Fig. 5, panel f). This a clear sign of the importance of infrared pumping that is included in the non-LTE models. The difference between the fluxes in the 4.25 μm line are greater, up to an order of magnitude. The difference are strongest in the models that have a low total CO₂ content (so low abundance and low gas-to-dust ratio). This is mostly due to the larger radial extent of the emitting region extending up to 20 times further out in the non-LTE models compared to the corresponding LTE model (Extent of the 4.3 μm non-LTE emission is seen in Fig. 5, panel i). This is in line with the higher Einstein A coefficient and upper level energy (and thus higher critical density) of the v₃,1 → 0,R(7) line, giving rise to a large importance of infrared pumping relative to collisional excitation. Figure E.1 uses g/d_dust = 1000, but the plot for g/d_gas = 1000 is very similar.

Appendix F: Line blending by H₂O and OH

One of the major challenges in interpreting IR-spectra of molecules in T Tauri disks are the ubiquitous water lines. H₂O has a large dipole moment and thus has strong transitions. As H₂O chemically favours hot regions (Agúndez et al. 2008; Walsh et al. 2015) there are a lot of rotational lines in the mid infra-red. Figure F.1 shows the H₂O rotational lines near the CO₂ 15 μm feature. The spectra are simulated with a LTE slab model using the same parameters as fitted by Salyk et al. (2011b) for AS 205 (N) as reproduced in Table F.1. It should be noted that AS 205 (N) is a very water rich disk (in its spectra) explaining the large number of strong lines. Fortunately, there are still some regions in the CO₂ spectrum that are not blended with H₂O or OH lines and thus can be used for tracing the CO₂ abundance structure independent of a H₂O emission model. The situation improves at higher resolving power as can be seen in Fig. F.2. The resolving power of 28000 has been chosen to match with the resolving power of the SPICA HRS mode. At this point the line widths are dominated by the assumed Keplerian linewidth of 20 km s^-1. At this resolution the individual Q-branch lines are separable and a lot of the line blends that happen at a resolving power of 2200 are no longer an issue.

Appendix G: Spitzer-IRS spectra

Figure G.1 shows the spectra as observed by Spitzer-IRS reduced with the CHIP software (Pontoppidan et al. 2010; Pontoppidan 2016). Continua that have been fitted to these spectra are also indicated on the figure. The objects have been chosen because their spectra are relatively free of H₂O emission. Even without H₂O lines, it is still tricky to determine a good baseline for the continuum as there are a lot of spectral slope changes, even in the narrow wavelength range considered here. This is especially true for HD 101412 where the full spectrum shows a hint of what looks like R- and P-branches. If these are features due to line emission, it becomes arbitrary where one puts the actual continuum, thus these features are counted here as part of the continuum. Whether these feature are real or not will not matter a lot for the abundance determination as the CO₂Q-branch is separated from the strong R- and P-branch lines.

Figure G.2 shows a comparison between observed spectra of disk with strong H₂O emission and CO₂ model spectra. The spectra are corrected for source luminosity and distance as explained in Sect. 4.1. Assumed distances and luminosities are given in Table G.1.

All Tables

All Figures