© ESO, 2013

1. Introduction

Terrestrial planets form in the inner region of protoplanetary discs (typically R < 3−5 AU), where the gas reaches temperatures of a few hundred to a few thousand degrees Kelvin and the dust grains are warm enough such that water is not frozen onto their surfaces (Kamp et al. 2010). The growth of solid cores in the inner disc is slower than in the outer disc, where grains are coated with an icy mantle. As a consequence, the solid bodies never reach the critical mass to attract gravitationally the disc gas before the gas disc has dissipated (Armitage 2010).

At a few hundred to a few thousand degrees, molecules are excited to their ro-vibrational levels and emit in the near- and mid-infrared. Probing the warm gas emissions requires high angular-resolution observations in the infrared. Carbon monoxide is one of the most abundant species in protoplanetary discs, and its ro-vibrational transitions around 4.6 micron (M-band) are the most commonly detected lines with high spectral resolution adaptive-optics (AO) assisted spectrometers around T Tauri (Najita et al. 2003; Rettig et al. 2004; Salyk et al. 2011) and Herbig Ae stars (Blake & Boogert 2004; Carmona et al. 2005; Brittain et al. 2009; van der Plas et al. 2009). The carbon monoxide (CO) peak emission is also seen displaced from the source centre, suggesting the presence of an inner hole (Brittain et al. 2003; Goto et al. 2006; Brown et al. 2012). The spectro-astrometric technique further improves our understanding of the emission by reaching spatial resolutions of the order of 0.1−0.5 AU for nearby (100−140 pc) discs (Pontoppidan et al. 2008). Goto et al. (2012) observed the CO ro-vibrational emission from the HD 100546 disc with a spatial resolution of 10 AU at 100 pc. They found warm CO gas (400−500) out to a distance of 50 AU, which is evidence for the presence of a warm molecular disc atmosphere. Very dense and hot gas located within a few tenths of AU from the stars can be detected by the CO overtone emissions (Tatulli et al. 2008; Berthoud et al. 2007; Thi et al. 2005).

Line profiles are diverse, from broad single peaked to clear double peaked; the latter is a signature of Keplerian rotation of the gas between 1 and 5 AU (Salyk et al. 2009). Bast et al. (2011) focused on 8 out of 50 discs that show broad single-peaked CO profiles. They proposed that those 8 discs have either highly turbulent inner disc gas and/or that a slowly moving disc wind contributes to the emission.

The M-band observations show lines emitted from a large range of excitation levels (from fundamental v = 1 → 0 P(1) and R(0) lines to “hot” lines (v′ > 1, Δv = 1) which make it possible to probe disc regions with different excitation conditions. CO vibrational diagrams of Herbig Ae discs indicate vibrational gas temperatures higher than the rotational ones (a few thousand K instead of a few hundred K). This suggests that the population of high vibrational levels (v >1) may be dominated by infrared (IR) and ultraviolet (UV) fluorescence excitation (Brittain et al. 2007; Brown et al. 2012). The v > 1 ro-vibrational levels can also be efficiently populated by UV fluorescence for T Tauri discs, especially for strongly accreting systems with excess UV emission (Bast et al. 2011). Fluorescence excitation in T Tauri discs is supported by Hubble Space Telescope observations of UV-pumped CO emission (France et al. 2011).

Many models of CO ro-vibrational emission made the assumption of local thermodynamical equilibrium (LTE) population (Hügelmeyer et al. 2009; Regály et al. 2010; Salyk et al. 2011). Nevertheless, a detailed understanding of the CO ro-vibrational lines, especially the spatial location of the lines and the efficiency of the IR/UV fluorescence, requires non local thermodynamic equilibrium (NLTE) modelling and a chemo-physical code that includes detailed continuum and line radiative transfer. The accuracy of NLTE modelling depends on many factors, including the availability of accurate collision rates. The difficulties associated with the measurement and computation of accurate rates result in sparse trustworthy data. Completeness in the rate coefficients can be attained only by using approximate extrapolation rules derived from physical considerations. Krotkov et al. (1980) and Scoville et al. (1980) studied the efficiency of UV and IR pumping but considered only pure-vibrational levels.

In this paper, we describe the implementation of a large CO model molecule, which includes many ro-vibrational levels in the ground electronic state X and the electronic state A, into the code ProDiMo. The radiative chemo-physical code ProDiMo solves the continuum radiative transfer, the gas heating and cooling balance, the gas chemistry, and the disc vertical hydrostatic structure self-consistently (Woitke et al. 2009, Paper I). This paper is the fourth in a series dealing with aspects in protoplanetary disc modelling (Kamp et al. 2010; Thi et al. 2011, Papers II and III).

We will use the code to model two typical discs around an A-type star: a solar-nebula-mass disc (M_disc = 0.01 M_⊙), and a low-mass disc (M_disc = 10^-4M_⊙). The choice of these models is justified by the expectation that the strength of UV/IR fluorescence will depend on the density and optical depth in the discs.

The paper is organized as follows: we present and discuss the CO spectroscopic and collisional data implemented in the ProDiMo code in Sect. 2. The main feature of the ProDiMo code and the disc model parameters are described in Sects. 3.1 and 3.2. In Sect. 3.3, we discuss the model outputs. Finally, we present our conclusions and recommendations on the interpretations of CO ro-vibrational observations in discs in Sect. 4.

2. CO molecular data

2.1. Spectroscopic data

Tashkun et al. (2010) compiled and evaluated existing experimental frequencies for transitions within the ground electronic state. Their work is based on accurate experimental frequency measurements (e.g. George et al. 1994; Gendriesch et al. 2009).

We compared their energy levels and frequencies with those from Chandra et al. (1996) and Hure & Roueff (1993) and found no significant differences. In addition, Chandra et al. (1996) computed Einstein spontaneous probabilities A¹Π for all ro-vibrational transitions between 140 rotational levels and 10 vibrational levels in the ground electronic state. We adopted their transition probabilities, which compared well with other studies (e.g. Okada et al. 2002). The transitions between the electronic excited ro-vibrational levels A and the ground ro-vibrational levels X¹Σ⁺ (4th positive system), which occur around ~1600 Å, were derived from the band-averaged transition probabilities of Beegle et al. (1999) and Borges et al. (2001) using the formula and Höln-London factors in Morton & Noreau (1994), which are given in Table 1. Figure 1 shows a sketch of the fluorescence mechanism involving a ro-vibrational level in the ground electronic state X¹Σ⁺ and a ro-vibrational level in the excited state A¹Π.

2.2. Collisional data

Although CO is one of the most studied interstellar molecules with a large amount of experimental and theoretical collision rates available, the lack of completeness of any individual study obliged us to use scaling laws to fill the gaps. Since protoplanetary discs show a large range of physical and chemical properties, we need to consider collisions of CO with hydrogen atoms (H), helium (He), hydrogen molecules (H₂), and electrons. In the rest of the section, we will first review existing theoretical and experimental collision rates before discussing the strengths and limitations of scaling laws.

2.2.1. CO-H collision rates

The v = 1 → 0 CO-H de-excitation rate coefficients were first measured by Millikan & White (1963) in shock-tube experiments. Subsequent works include the shock-tube experiments by von Rosenberg et al. (1971), Glass & Kironde (1982), and more recently Kozlov et al. (2000) (between 2000 and 3000 K). Shock-tube experimental measurements are de-excitation timescales,which can be fitted in terms of the Landau-Teller rate coefficients in cm³ s^-1 (Ayres & Wiedemann 1989): (A.1)with β = E_{v   =   1}/T and the values for A and B given in Table 2. Neufeld & Hollenbach (1994) gave a simple formula that fits the shock-tube data for T > 300 K: (A.2)This formula supersedes an older one given by Hollenbach & McKee (1979): (A.3)Experimental values for transitions other than the v = 1 → 0 are not available. Balakrishnan et al. (2002) performed quantal scattering computations for rotational and band-averaged vibrational collision rates between CO and hydrogen atoms. They used the close-coupling method for the rotational transitions and the infinite order sudden approximation (IOS, Flower 2007) for the vibrational transitions for a large range of gas temperatures (5 < T < 3000 K). The IOS approximation reduces the computational effort but is strictly valid only for energetic collisions, i.e. ΔE ≪ kT, where ΔE is the energy difference between two levels. The close-coupling computations were done for rotational levels up to J = 7. The vibrational transition rates were computed up to v = 4.

Figure 2 shows the different rate coefficients. The most striking feature is the large dispersion (a few orders of magnitude) between the sources. Ayres & Wiedemann (1989) chose the values of Glass & Kironde (1982), which are much larger than the older values of von Rosenberg et al. (1971). However, recent measurements by Kozlov et al. (2000) reproduced the values of von Rosenberg et al. (1971). Interestingly, the theoretical values calculated by Balakrishnan et al. (2002) are much closer to the values of Glass & Kironde (1982) and are best approximated by the formula of Neufeld & Hollenbach (1994). One possible reason for the large discrepancies is the reactivity of CO with atomic hydrogen (an open shell species), which can falsify the measurements. We chose to use the theoretical rates by Balakrishnan et al. (2002).

2.2.2. CO-He collision rates

Cecchi-Pestellini et al. (2002) computed rotational and vibrational rate de-excitation coefficients of CO in collision with He. The rotation rate coefficients in the ground vibrational level are computed by the quantal close-coupling method. The vibrational rates are computed for 500 ≤ T ≤ 1300 K using the IOS approximation. Krems (2002) calculated quantum-mechanically ro-vibrational rate coefficients for 35 ≤ T ≤ 1500 K. Rates were computed by Flower (2012) and Reid et al. (1995). The experimental fit coefficients of Millikan (1964) are given in Table 2. Wickham-Jones et al. (1987) reported experimental data from the Ph.D. thesis of Chenery (1984), who measured rates down to 35 K. The experimental rate coefficients are compared to the theoretical values in Fig. 3. The theoretical values of Cecchi-Pestellini et al. (2002) are a factor ~2 smaller than the experimental values, whereas the values computed by Krems (2002) match the shock-tube data down to T = 300 K. The formula (A.1) that fits the experimental shock-tube data under-predicts the values for T ≤ 300 K. We adopted the values of Krems (2002) up to T = 1500 K and those of Cecchi-Pestellini et al. (2002) above T = 1500 K. The rate coefficients for CO in collisions with He are around three orders of magnitude smaller than the collision rates with H. This behaviour can be explained by the fact that He is a close-shell species.

2.2.3. CO-H₂ collision rates

In the case of collisional excitation of CO by H₂, we need to consider the collisions with ortho-H₂ and para-H₂ separately. Yang et al. (2010) computed with the close-coupling method rates for temperatures between 1 and 3000 K involving ground vibrational CO transitions and with rotational levels up to J = 40. Rates for vibrational transitions are scarcer and are given averaged over the rotational levels. Between 30 and 70 K, Reid et al. (1997) showed experimental and theoretical rates for the CO v = 1 → 0 transition. Between 70 and 290 K, rates derived from shock-tube experiments are available (Hooker & Millikan 1963; Andrews & Simpson 1975, 1976; Glass & Kironde 1982). At temperatures higher than 300 K, we can assume that the Landau-Teller law is valid and expand the shock-tube rates to a few thousand K (Ayres & Wiedemann 1989). Neufeld & Hollenbach (1994) provided a convenient analytical expression of the shock-tube data (T > 290 K): (A.4)Alternatively, one can use the Landau-Teller formula with the values in Table 2.

In Fig. 4, we compare the two shock-tube experimental rate coefficients, the experimental and theoretical values of Reid et al. (1997), and our adopted fit (A.5)for T ≤ 300 K and the fit to the Hooker & Millikan (1963) data by Neufeld & Hollenbach (1994) above 300 K. We notice a large discrepancy between the experimental and theoretical rates at low temperatures. The shock-tube data are consistent with each other. It is known that the rate extrapolations to lower temperatures with the Landau-Teller law may not be correct. Collisions with H are around two orders of magnitude more effective than collisions with H₂ at all temperatures.

2.2.4. CO-electron collision rates

The collision excitation and de-excitation rates of CO by electrons for rotational transitions are computed using the theory of Dickinson et al. (1977) based on the Born approximation for ΔJ = ± 1 (see also Randell et al. 1996). We adopted a rotational constant B_CO =  1.922529 cm^-1, a centrifugal distortion D_CO = 6.1210779 × 10^-6 cm^-1, and a dipole moment of 0.112 Debye. Ristić et al. (2007) measured the vibrational rates of CO by collisions with electrons for vibrational levels up to ten. Morgan & Tennyson (1993) also computed cross sections for vibrational transitions. Thompson (1973) gave a simple formula for the v = 1 → 0 rate (A.6)where again β = E_{v   =   1}/T. An alternative formula is given by Draine & Roberge (1984)(A.7)Finally, we propose our own analytical formula, which fits the values of Ristić et al. (2007)(A.8)The measurements for the v = 1 → 0 transition are shown in Fig. 5 and compared to the three analytical fitting formulae. The Thompson (1973) formula deviates significantly from the experimental data. We used a combination of formula (A.8), which fits the experimental data relatively well down to 400 K. Below 400 K, we adopted a flat rate coefficient.

2.3. Rate coefficient scaling laws

The current experimental and theoretical works do not provide all the rate coefficients needed to model the CO statistical ro-vibrational population of the ground X¹Σ⁺ and excited A¹Π electronic levels. Therefore, we have to resort to extrapolation and scaling laws to fill the gaps.

2.3.1. Extrapolating rotational transition rates

A few extrapolating laws attempt to estimate rates from existing ones. We used the energy-corrected sudden (ECS) scaling law with a hybrid exponential power-law for the base (k_{J′′ → 0}) rates (Depristo et al. 1979; Depristo 1979; Green 1993). Goldflam et al. (1977) showed that rates between any levels can be derived knowing the rates to the ground level within the IOS approximation. To overcome the limitation of collision energies greater than the level energy differences, Depristo et al. (1979) proposed a correction for the inelasticity and finite collision duration effects using second-order perturbation theory for low-collision energies (i.e. for collisions at low temperatures). A de-excitation (J′ > J′′) rate is given by the ESC-EP (energy-corrected sudden scaling law with a hybrid exponential power-law for the base rates) formula: (A.9)where (A.10)is the Wigner 3-j symbol and E_J′ and E_{J′′   =   J′ − 1} are the upper and lower level energy, respectively. The correction factor introduced by Depristo et al. (1979) is defined by (A.11)with (A.12)where μ is the reduced mass in atomic mass units, ω_J′ is the inelasticity to the next lowest level (E_J′ − E_{J′′   =   J′ − 1}) in cm^-1, and l_c is an adjustable scaling length in Å. We adopted l_c = 3 Å (Schöier et al. 2005). The ESC formula is recovered for l_c = 0 Å.

The scaling law requires the knowledge of the base rates (J → 0). Collision rates with hydrogen atoms to the ground level are known for J up to 7, with H₂ up to 40, and only for J = 1 → 0 for electrons. Further estimates of the rates to the ground level are needed. The ESC-EP model provides an approximation for rates to the ground level (Green 1993; Millot 1990): (A.13)where A(T), γ, and β are free parameters, which are estimated by fitting Eq. (A.13) to the existing data using a least-square solver for over-constrained systems. We performed the extrapolation for collisions with o-H₂, p-H₂, He, and H but not with electrons. An example of extrapolated rate coefficients is shown in Fig. 6.

2.3.2. Extrapolating vibrational transition rates

The extrapolation of an existing rate to low temperature was performed by fitting the available rates to the formula (A.14)where a and b are the two fitting parameters.

The extrapolation of rates between two vibrational levels without any known rate follows physically motivated scaling laws. Based on the initial work of Procaccia & Levine (1975), Chandra & Sharma (2001) proposed using the Landau-Teller relationship for transitions between adjacent states (Δv = 1) (A.15)and for Δv > 1: (A.16)Elitzur (1983) argued for a different formulation (A.17)For Δv = 1 (v′′ = v′ − 1), this last formula reduces to (A.18)Scoville et al. (1980) argued that the probabilities of collision-induced vibrational transitions are proportional to the corresponding radiative transition matrix elements (Born-Coulomb long range interaction). Therefore, the vibrational de-excitation rate coefficients from v′ to v′′ can be derived from the v = 1 → 0 rate coefficients: (A.19)where A_v′v′′ is the vibrational band transition probabilities. The three formulations give quite different rate estimates. In order to choose which one is the most appropriate, we compared the predictions with detailed calculations and/or experimental values.

Figure 7 shows the rate coefficient ratios for transitions between adjacent levels in collisions with atomic hydrogen using the detailed calculation of Balakrishnan et al. (2002). The ratios are grouped around one with a small spread. The rate coefficient ratios suggest that the expansion rules proposed by Elitzur (1983) are more appropriate to expand the detailed rate coefficients.

The rate coefficient ratios for CO-He collisions vary dramatically with temperature (Fig. 8). For temperatures below 1000 K, the ratios are close to 2 × v′, a value that corresponds neither to the Chandra & Sharma (2001) prescription (=v′) nor to that championed by Elitzur (1983) (=1). We chose to use the Chandra & Sharma (2001) for Δv > 1 and the factor 2 × v′ for Δv = 1.

The choice is open concerning the rules to apply to estimate the rate coefficients between CO and H₂ and between CO and electrons other than v′′ = 1 → 0. We chose to use the same combination of scaling rules as for CO-He collisions. We estimate the scaling rules to provide rates correct within an order of magnitude.

2.3.3. Estimating ro-vibrational transition rates

The vibrational rates are averaged over all ro-vibrational transitions. To obtain individual rates between ro-vibrational states, we apply the method described in Faure & Josselin (2008). This method assumes a complete decoupling of rotational and vibrational levels. The state-to-state ro-vibrational rates are related to the corresponding rates for the pure rotational transitions in the ground vibrational state: (A.20)The factor P_v′v′′ is defined as (A.21)This formula ensures detailed balance.

2.3.4. Approximating electronic transition collision rates

In the absence of collision rates connecting the ground X¹Σ⁺ to the A¹Π electronic levels, we adopted the rough values for CO-electron collision rates provided by the approximation. The (van Regemorter 1962) approximation relies on the Born-dipole (also called Born-Coulomb) long-range interaction approximation and relates the collision strength to the transition probability and wavelength for allowed transitions (Itikawa 2007). The Born-dipole approximation is valid for species with a permanent dipole, which is the case for CO (μ = 0.122 Debye). The collision rate for allowed transitions is given by (A.22)For neutrals, the g-bar value can be approximated by (A.23)The collision strength Γ is related to g-bar by (A.24)The approximation provides rate coefficients correct within an order of magnitude.

2.3.5. Combining theoretical, experimental, and estimated rates

The resulting CO rate matrix after combining all existing and estimated rates is complete in the sense that all ro-vibrational levels in the ground electronic state connected by a radiative transition are also connected by collisions. Therefore all levels in the ground electronic levels can be populated in LTE at very high densities and temperatures. We have neglected the radiative and collisional transitions between levels in excited electronic states since radiative transitions to the ground electronic state are very fast. The “completeness” is required to ensure that we don’t artificially create sub-thermally populated levels. The methodology described above was implemented in the radiative photochemical code ProDiMo. The main collisional partner in ro-vibrational transitions is the atomic hydrogen by a few orders of magnitude.

3. CO fundamental and hot-band emission from Herbig Ae discs

3.1. The P roD iM o code

A detailed description of the code can be found in Woitke et al. (2009). Subsequent additions to the code are documented in Kamp et al. (2010), Thi et al. (2011), and Woitke et al. (2011). Although not used in this study X-ray heating and X-ray chemistry are included in ProDiMo (Aresu et al. 2011). We deal in this paper with the excitation mechanisms of the CO ro-vibrational levels and do not attempt to fit actual observations.

3.2. Disc model description

We chose to model two representative discs, which can be found around Herbig Ae stars. Our generic star has an effective temperature of 8600 K, a mass of 2.2 M_⊙, and luminosity of 32 L_⊙. We have no intention to fit a particular object in this study. The input stellar spectrum is taken from the phoenix database of stellar spectra (Brott & Hauschildt 2005). The disc model parameters are summarized in Table 3. The disc geometry is described by a single zone power-law surface-density profile with index ϵ, an inner radius R_in, an outer radius R_out, a disc gas mass M_disc, and a standard gas-to-dust mass ratio Δ of 100 for all the discs. The inner radius was set to two values to model the presence of an inner disc hole and its effects on the CO ro-vibrational emissions.

The dust grains are defined in both discs by their composition (astronomical silicates, Laor & Draine 1993) and their size distribution (a power-law ranging from a_min to a_max). We did not iterate on the vertical hydrostatic profile but used instead a power-law H = H₀(R/R₀)^− p to model the gas scale-height H. We chose to model two flaring discs (flaring index p = 1.2) because observations suggest that the CO levels populated by UV-fluorescence pumping seem to occur in flaring discs.

The two discs differ from each other only in the total (gas+dust) disc mass M_disc. The first disc has a mass consistent with the typical solar nebula mass of 0.01 M_⊙, while the second disc is more akin to an object in transition from a primordial gas-rich disc to a secondary debris disc with a total mass of 10^-4 M_⊙. The choice of the two masses is motivated by the possibility to assess the effects of infrared and UV fluorescence on the CO population in discs with vastly different density structures. The discs are assumed passively heated only (i.e. without viscous heating: α = 0), and we do not consider X-ray heating because the X-ray luminosity is generally low compared to the UV luminosity in Herbig Ae stars. X-ray heating and chemistry can play an important role in discs around T Tauri stars (Aresu et al. 2011, 2012; Meijerink et al. 2012).

The main source of heating for the gas is the stellar UV via photoelectric effects on polycyclic aromatic hydrocarbons (PAHs) assumed here to be circumcoronenes (C₅₄H₁₈). The abundance of PAHs is assumed at 10% of the standard interstellar abundance (f_PAH = 0.1).

The gas and dust were assumed well-mixed with no dust settling. A detailed discussion on the heating and cooling in HerbigAe discs can be found in Kamp et al. (2010).

The continuum radiative transfer and dust thermal balance are calculated in two dimensions using the long-characteristic method (Woitke et al. 2009; Thi et al. 2011). The radiative transfer numerical implementation was benchmarked against other radiative transfer codes (Pinte et al. 2009). The gas temperatures were computed by balancing the heating and cooling (mostly radiative) processes. The CO self-shielding against photodissociation uses factors interpolated from values in look-up tables (Visser et al. 2009).

CO ro-vibrational transitions have moderate Einstein probabilities (A = 1−10 s^-1) so that the line transfer occurs in the Doppler core of the line profile. On the other hand, UV-fluorescence pumping occurs via electronic transitions (A¹Π–X¹Σ⁺) with high Einstein probabilities (A = 10⁴−10⁷ s^-1) and thus with large natural line width. The line core becomes quickly optically thick (line self-shielding) with the UV line transfer happening in the Lorentzian wings of the line. Therefore the use of the Voigt profile in the line transfer is warranted.

The CO levels are computed using a 1+1D escape probability method for Voigt line profiles. We adapted the analytical formulation of Apruzese (1985) to match the formulation in Woitke et al. (2009) at low optical depths and small intrinsic line widths. At high optical depth, the escape probability β varies as τ^− 1/2 (Voigt profile) instead of τ^-1 (Gaussian profile): (A.25)where a is defined as Γ/(4πΔν_D) and τ_c is the critical optical depth and is defined as . The sum of the natural and collisional width is Γ while the effective Doppler width is Δν_D (Rybicki & Lightman 1986). The escape probability function for a Voigt line profile is illustrated in Fig. 9 for the case a = 10^-3. In addition to the UV line shielding, dust grains contribute strongly to the UV flux attenuation. The main limitation in using the escape probability technique is that it does not take overlapping line effects into account.

Excited vibrational levels can also be populated by the absorption of an IR photon emitted by the star and the dust grains, a phenomenon called IR pumping. We also checked the effect of IR pumping on the vibrational ground rotational population by modelling two discs with R_in = 1 AU, one with M_disc = 10^-2M_⊙ (models 1c, 1d, and 1e) and with M_disc = 10^-4M_⊙ (models 3c, 3d, and 3e). For each disc model, we computed the rotational level population by taking into account: (1) ro-vibrational levels up to v = 5 (models 1c and 3c); (2) ro-vibrational levels with v = 0 and v = 1 (models 1d and 3d); and (3) the rotational levels in the ground vibrational only (models 1e and 3e).

In total, we ran four series of models with and without the electronic levels and with different numbers of vibrational levels to assess the importance of UV fluorescence and IR pumping (see Table 4). Complete spectra from four to five microns are generated for discs seen face-on (zero inclination).

3.3. Results and discussion

3.3.1. CO chemistry and location of the CO ro-vibrational emission

We focus on the CO chemistry of model 1a. The disc density, dust temperature, and gas temperature structures are shown in Fig. 10. The panels in Fig. 11 show the enhancement (χ) with respect to the interstellar UV, the abundance relative to the total number of hydrogen-nuclei for He, electron, atomic hydrogen, molecular hydrogen, vibrationally excited H₂ (H2exc), C⁺, C, CO, CH₄, OH, and H₂O. The transition from C⁺ to C occurs at log (χ/n) between −1 and −2, and the transition between C and CO occurs at log (χ/n) ~ −2.5 to − 3.

In Fig. 12, we show the location of the CO v = 1−0 P(19) and v = 4−3 P(20) emissions overlaid with the number density of the three main collision partners (H, H₂, and electrons) for a 10^-2M_⊙ disc model taking UV pumping or not into account (models 1a and 1b). Similar plots are displayed in Fig. 13 for models 3a and 3b. The CO lines are emitted up to a few AU after R_in, where the density is high enough such that the level population is the closest to the LTE population for the fundamental transitions. The emitting area becomes more and more confined to the inner rim as the upper energy of the transition increases. In the 10^-2M_⊙ disc model the fluxes are a few percents higher when UV pumping is switched on. The fundamental P(19) line is emitted in region where the gas is between 500 and 1000 K. In the vertical direction, the CO ro-vibrational lines are mostly emitted in the disc region with z/r ~ 0.2 and R < 5 AU.

In the emitting region, the extinction in the visible A_V is below 1 and some H₂ molecules are excited by UV and/or hot enough to be vibrationally excited (the excited H₂ abundance reaches 10^-12–10^-9). As soon as it is formed, the warm or vibrationally excited H₂ reacts quickly with atomic oxygen to form OH. In turn, OH reacts with H₂ to form H₂O or with atomic carbon to form CO.

The OH and water vapour abundance are high for z/r ~ 0.2 and R < 5 AU. The central role played by OH in hot gas chemistry is discussed in more details in Thi & Bik (2005). The 1a and 1b disc models are massive enough such that in the inner disc midplane (z/r < 0.1) the most abundant carbon-bearing gas-phase species are methane and C₂H₂ (R < 3 AU), whereas most of the oxygen is locked in water vapour up to R ~ 5 AU. The remnant oxygen is in the form of atomic oxygen, CO, and other minor species. The gas in the inner disc is dense enough such that UV fluorescence is quenched by collisional de-excitations between the ro-vibrational levels.

3.3.2. CO spectral line energy distribution

The analysis of the CO ro-vibrational lines can be performed either by comparing directly the observed and modelled line fluxes or by using rotational diagrams. In Figs. 14 and 16 we plotted the continuum-subtracted line fluxes for model 1 to 4 (spectral line energy distribution, SLED) for the fundamental (v = 1−0) and one “hot” transition (v = 4−3). The shape of the SLEDs varies from model to model and between models with and without UV pumping. The large variety of the SLEDs illustrates the sensitivity of the CO ro-vibrational line fluxes to the disc parameters (M_disc, R_in, etc.). The line fluxes range from 10^-19 W m^-2 to a few 10^-16 W m^-2 for a low-mass disc (10^-4M_⊙) with an inner hole to a few 10^-16 W m^-2 to for a massive disc (10^-2M_⊙) at a distance of 140 pc. In comparison, we show in Fig. 15 the v = 1−0 CO line fluxes (Brittain et al. 2003) for the massive disc AB Aur and the low-mass disc HD 141569A.The observed fluxes from HD 141569A should be scaled from its distance of 99 pc (van den Ancker et al. 1997) to 140 pc. Fluxes for the hot lines from the HD 141569A disc are ~10^-18 W m^-2 after scaling the fluxes to 140 pc. When UV pumping is included, the hot line fluxes from our models are of the order of 10^-19–10^-17 W m^-2. The shapes of the SLEDs for models with a small inner radius (R_in = 1 AU, models 1 and 3) reflect the effects of highly optically thick lines. On the other hand, the shapes of SLEDS of models with a large inner radius (R_in = 20 AU, models 2 and 4) are typical of optically thin or moderately optically thick lines. The observed line fluxes are best matched by models with a small inner radius whereas the shapes of the observed SLED are closer to models with a large inner radius. To obtain fluxes that are close to the observed values, the emitting gas should be close to the star, but at the same time lower surface densities are needed to have moderate optical depths for the CO lines. Interestingly, UV pumping will increase the line v = 1−0 fluxes only by a factor of few and the v = 4−3 fluxes by orders of magnitude. In our models, UV pumping does not dramatically change the shape of the SLEDS for v = 1−0 and v = 4−3 transitions.

We ran models with inclination from 1 to 89 degrees with the line ray-tracing module. The line fluxes vary by up to 30% for inclinations below 70 degrees. At high inclination, the cold CO in the outer flaring disc starts to reabsorb the emission from the inner disc.

3.3.3. CO ro-vibrational level excitation

The CO vibrational levels in the ground electronic state can be populated by collisions with the main chemical species (H, H₂, electrons, He), by formation pumping (i.e. formation of CO in an vibrationally excited level), by pumping by dust emission, or by UV-fluorescence pumping. The region of high warm CO abundance is H rich and H₂ poor (Kamp & Dullemond 2004) (see Figs. 12 and 13). We also plotted the excitation temperatures. We chose two typical transitions, one from the v = 1 level and one from the v = 4 level, to illustrate the efficiency of the excitation mechanisms. The relative population between the v = 1,J = 19 and v = 1,J = 18 and between v = 4,J = 20 and v = 4,J = 19 levels are shown as excitation temperatures in Figs. 17 and 18 for the high- and low-mass disc models, respectively. The excitation temperatures between the v = 1,J = 19 and v = 0, J = 18 levels in the line-emitting region are below the gas kinetic temperatures but higher than the radiation temperature at the line frequency 4.844 micron (subthermal population, see Figs. A.2, A.4 and A.5 in the appendix). However, the excitation temperatures between the v = 4,J = 20 and v = 3,J = 19 levels are much higher than the gas kinetic and radiation (at 4.692 micron) temperatures (suprathermal population). We can distinguish between the excitation mechanisms of the v = 1 level and that of the v > 1 levels.

The v = 1 level is mostly populated by collisions and IR pumping because the disc densities in those emitting regions are close to the critical densities (10¹⁰–10¹² cm^-3). The v = 1 level population is weakly affected by UV-fluorescence pumping (see the upper panels of Figs. 17 and 18). Likewise the line fluxes for the fundamental transitions change by at most a factor four when UV-fluorescence pumping is taken into account (Fig. 14).

The collision rates with electrons are ~100 times higher than with H but the electron abundance is ~10^-6 lower than the H-atom abundance. The collision rates with He are of the order of 10^-16 cm³ s^-1 at 500 K. The abundance is 0.075 times less than H+H₂. Therefore He is not a major collision partner. The main collision partner in the CO ro-vibrational emitting region is atomic hydrogen, whose de-excitation collision rates with CO are two orders of magnitude (10^-12 cm³ s^-1 at 500 K) higher than rates with H₂.

The v > 1 levels have higher critical densities than the v = 1 levels and require very high gas temperatures T_gas > 1000 K and high densities to be populated efficiently by collisions. Pumping of the v > 1 levels by UV fluorescence becomes the prominent excitation mechanism. The excitation temperatures of UV-pumped levels (lower left-hand panels in Figs. 17 and 18) tend towards the effective radiation temperature (brightness temperature) in the UV (Fig. A.3, the radiation temperatures at 0.15 micron).

The effect of UV pumping becomes more pronounced in the 10^-4M_⊙ disc models (3a and 3b, Fig. 13), where the gas densities and hence the UV-pumped high-v level quenching efficiencies are lower and the UV photons are less absorbed by the dust grains (Fig. A.3). In the M_disc = 10^-2M_⊙ disc model, the UV pumping is confined to a very thin layer at the inner disc rim, whereas in the low-mass disc models (M_disc = 10^-4M_⊙), the UV pumping occurs till ~0.5 AU inside the disc because the UV photons can penetrate deeper in the disc.

The UV-fluorescence pumping can affect dramatically the line fluxes (Fig. 16). In particular the hot line v = 4 − 3 P(20) is two orders of magnitude stronger with UV pumping. The effect of UV pumping is better captured by computing the vibrational temperature (see Sect. 3.3.5). The relative rotational level populations within a vibrational level are less affected by UV pumping than the relative vibrational level populations (see the shape of the SLED in Fig. 14).

An alternative excitation mechanism is the formation pumping of CO. The CO formation reaction via C + OH has an exothermicity of − 6.5 eV and a rate of k_chem ~ 10^-10 cm³ s^-1 (Zanchet et al. 2007). The excess energy is used to pump CO to high vibrational levels with a maximum population at v = 10 (Bulut et al. 2011). The formation of CO constitutes a pumping mechanism for CO.

We can compare the efficiency of the chemical pumping with the efficiency of the population of the v = 1 level by collisions with atomic hydrogen by noting that OH reaches a maximum abundance of ~10^-6 (Fig. 11). The chemical pumping rate assuming atomic carbon abundance [C] ~ 10^-4 is k_chemn_Cn_OH ~ 10^-20n² cm³ s¹, where n is the gas density. Collision rates with atomic hydrogen are of the order of 10 cm³ s¹, a few orders of magnitude more efficient than the chemical pumping.

3.3.4. Rotational diagram for the fundamental transitions

The analysis of the line fluxes using a rotational diagram can provide more insight than the study of the shape of the SLEDs. We plotted in Fig. 19 the fundamental transition rotational diagram, i.e. the log of the level column density versus B_COJ′(J′ + 1), where B_CO is the rotational constant in Kelvins, for the 10^-2 and 10^-4M_⊙R_in = 1 AU and 20 AU models with (blue symbols) and without UV pumping (red symbols). The shape of a rotational diagram is affected by NLTE (subthermal population and IR/UV pumping) and optical depth effects. To obtain the optically thin emission in LTE (purple solid lines), we populated the CO ro-vibrational levels assuming T_ex = T_gas. The line fluxes were computed in the optically thin approximation, i.e. optically depth effects, which would have prevented the flux from increasing indefinitely, were not taken into account.

The optically thin CO LTE populations probe the actual kinetic temperatures in the discs. The populations were calculated by summing for each level all the CO molecules at that level in the models, assuming that the population is in LTE. For each model, the rotational population of the v = 1 level can be relatively easily matched by a single temperature because the CO fundamental lines probe a small region in discs where the gas is warm and dense. Moreover, the densities are high enough for the low- and medium-J rotational levels within the v = 1 level to be close in rotational LTE. The line profiles will be wider for high-J than for low-J transitions in discs not seen edge-on. The LTE temperatures (fits to the purple lines shown with black crosses) T_LTE are given in Table 4. The gas is warmer for lower mass, hence lower density discs because the molecules, which are efficient gas coolants, are more rapidly formed in dense gases.

Since the stellar radiation decreases with distance squared, it is normal that the gas is cooler for discs with an inner radius starting at 20 AU rather than at 1 AU. Since CO ro-vibrational lines are efficient coolant for the gas, the gas probed by the CO gas is warmer (i.e. with higher ) when the number of ro-vibrational levels decreases (Table 4).

The red and blue lines in Fig. 19 correspond to the CO ro-vibrational NLTE level population computed locally in the disc using the 1+1D escape probability method. Those lines show that the level populations are rotationally subthermal. The effects are stronger for less massive (10^-4M_⊙), hence less dense discs.

The ro-vibrational NLTE population is orders of magnitude smaller than the ro-vibrational LTE population for all models, pointing to a subthermal population of the excited ro-vibrational levels: the rotational level populations within a vibrational level can be in LTE, while the vibrational level populations are subthermal. Assuming a transition probability of 10 s^-1 and a rate coefficient of 10^-16 cm³ s^-1 (H₂) to 10^-12 cm³ s^-1 (H), the critical density is 10¹¹–10¹⁵ cm^-3. Such high densities are reached only close to the star in the inner disc midplane of massive discs.

Radiative transfer of the optically thick lines decreases further the populations derived from the analysis of rotational diagrams. One consequence is that rotational temperatures derived from fitting points from low-J optically thick lines in rotational diagrams would appear lower than the actual gas kinetic temperatures. In Appendix A.1, we provide a simple analytical method to analyse observations of optically thick ¹²CO v = 1 − 0 emission lines. Optical thickness results in flux differences between lines of transitions that have the same upper level but lower levels that differ by ΔJ and different transition probabilities: the P (J′′ = J′ + 1) and the R (J′′ = J′ − 1) branches. As expected, the column density differences between LTE and NLTE population are strong in models 1a and 1b and are weak in model 4a and 4b.

The R-branch transitions generally have higher transition probabilities (i.e. with higher Einstein-A coefficients) than P-branch transitions. Hence, if the lines are optically thick, the column densities derived from R-branch lines are the smallest (N_J′ ∝ F_ν/A_ul and F_ν are almost the same). Lines are optically thick up to very high J, giving rise to a steep slope for the low J and to a shallow slope for the moderate J lines (symbols). Observed rotational diagrams of Herbig Ae discs show similar shapes (Blake & Boogert 2004; Brittain et al. 2007; Brown et al. 2012). Only lines from J′ > 40 (>25 for the 10^-4M_⊙, R_in = 20 AU disc model) are optically thin: their populations join the optically thin NLTE populations and the differences between the P and R branch column densities fade away. The optical thickness of the lines and rotationally subthermal population of the levels can cause the lines to be actually optically thick but effectively thin, especially in discs where the continuum emission at 4−5 micron is optically thin.

UV pumping affects the level population for all levels with v > 1 whatever the inner disc radius. The main effect of UV pumping is to populate the high-J and v > 1 (see discussion below and Fig. 20) levels at the expense of the lower v = 0,J levels. We also ran the models with UV pumping but without collision excitations between the ground and electronic excited states. However, we noticed no difference over a few percents.

The shape of the rotational diagrams (symbols) shows that it would be difficult to fit straight lines to estimate the excitation temperatures. Moreover, they do not reflect the actual CO population due to strong optical depth effects.

3.3.5. Vibrational temperature and UV/IR fluorescence pumping

The vibrational temperature measures the relative population of the vibrational levels. We plotted in Figs. 20 and 21 the log of the vibrational level population (ln(N_v)) as function of the vibration energy E_v = 3122 v, where v is the vibrational quantum number (Brittain et al. 2007). The vibrational level population is the sum of all the rotational level populations at that vibrational level.

We derived the vibrational temperatures by fitting straight lines through the data in models with and without UV pumping (nine and five vibrational levels) apart from the 10^-2M_⊙R_in = 1 AU model where only v = 4, 5, and 6 levels are taken into account because a straight line could not be fitted when those levels were included. The values are reported in the last column of Table 4.

Even without UV pumping, the vibrational temperatures are higher than the LTE rotational temperatures. However, the rotational temperatures derived using points in the shallow slope of the diagrams will yield temperatures higher than the LTE values and closer to the non-pumped vibrational temperatures.

Two clear trends can be seen from the vibrational diagrams and the temperatures in Table 4. First, the vibrational temperature is higher for the lower mass disc. Second, the vibrational temperature increases with the inner hole size if UV pumping is switched on, whereas it decreases if UV pumping is switched off.

When the UV pumping is not present, the high vibrational levels can be populated by collisions and/or by IR pumping. The collision rates for vibrational transitions are lower than for rotational transitions. If the vibrational levels are populated by collision only, then the vibrational temperatures should be lower then the rotational temperatures. Therefore, collisions play a minor role in populating the high-v ro-vibrational levels, leaving IR pumping as the dominating mechanism. In the absence of UV pumping, the vibrational temperatures are still higher than the maximum dust temperatures in the discs (max(T_dust) in Table 4). The IR continuum flux between 1 and 5 μm is dominated by the dust emission for the models with R_in = 1 AU (models 1 and 3) and by the stellar emission for models with an inner hole (R_in = 20 AU; models 2 and 4). The vibrational temperatures are the same whether five or nine vibrational levels are taken into account in the models, suggesting that the effect of IR pumping does not depend on the vibrational level.

As the gas cools with increasing inner radius, the high-v levels are less populated. Absorption of UV photons directly by CO molecules will favour the population of high-v levels unless high gas densities quench the overpopulation. The gas density at the inner edge decreases with the edge’s distance to the star. Higher R_in values and lower mass discs will have less collisional de-excitations of UV-pumped high-v levels.

3.3.6. Effect of IR pumping on the rotational population of the v = 1 level

The left-hand panels in Fig. 22 show the differences between the case with two and five vibrational levels. For the massive disc with M_disc = 10^-2M_⊙, the population decreases as more vibrational levels are included. The CO v > 1 levels can be pumped by absorptions in the hot and overtone bands of stellar or dust IR photons, resulting in lower populations of the v = 1, J levels (upper right-hand panel). The pure rotational levels in LTE cannot be explained by a single excitation temperature contrary to the rotational levels in the first vibrational excited state. The pure rotational transitions probe much larger disc areas with wide-ranging physical conditions than do the fundamental ro-vibrational transitions. This stems from the low critical densities needed to populate the pure-rotational levels.

The right-hand panels in Fig. 22 compare the population in models where only rotational levels with v = 0 are included and where the levels v = 0 and v = 1 are included. In the optically thin cases (solid lines), the differences between the cyan and blue lines reflect the rotationally subthermal populations for the high-J populations. The low-J levels are at LTE. However, the high optical depths for the low- and mid-J lines result in lower derived population and the appearance of a change in the slope at J = 8–12. At high J (E_rot > 2000 K), the overpopulations when only the ground vibrational (v = 0, J) are modelling artifacts because overlapping energies occur between the v = 0, high-J level energies and that of the v = 1, low-J levels. For a gas at temperature T, both the high J and low J of the v = 1 levels can and should be populated. In warm and hot gases, models that do include the pure rotational levels only may overestimate their populations.

4. Conclusion

We have implemented a complete CO ro-vibrational molecular model in the ground X¹Σ⁺ and A¹Π electronic state in the radiative chemo-physical code ProDiMo. We gathered existing collisional rate coefficients and used scaling rules to extrapolate the missing ones. Collision rates of CO with atomic hydrogen are two orders of magnitude higher than rates with molecular hydrogen and three orders of magnitude higher than with He.

We ran models that include continuum, NLTE-line radiative transfer assuming a Voigt line profile, chemistry, and gas thermal balance to understand the effect of UV on the population of high-vibrational and high-rotational levels in protoplanetary discs. CO ro-vibrational lines are emitted within the first few AU and probe disc regions where CO is quickly formed via C + OH, with OH being formed via the fast reaction of vibrationally excited H₂ with atomic oxygen.

The amount of atomic hydrogen is high in the CO line-emitting region. The v = 1 level population is dominated by collisional excitation with atomic hydrogen with contributions from IR-pumping and UV-fluorescence pumping. Excitation of high-J and high-v levels by UV and IR photons can be efficient, especially in low mass discs where collisional quenching is less rapid than in high-mass discs. The molecular gas is rotationally “cool” and subthermal but vibrationally “hot” and suprathermal. However, only the detection of hot lines (v′ ≥ 2, Δv = 1) in discs can provide a test of the efficiency of UV-fluorescence.

The pure rotational transitions probe the entire disc while the ro-vibrational transitions probe the inner warm disc region. Therefore the rotational temperature depends on the vibrational number of the initial level. It is difficult to derive accurate temperatures and column densities from the optically thick ¹²CO fundamental emissions. Spatially and spectrally resolved CO, ¹³CO, C¹⁸O, and C¹⁷O ro-vibrational observations are needed to disentangle the optical depth effects and constrain the location of the line emissions. The low- and medium-J rotational levels on the ground and first excited vibrational levels are populated in LTE in discs due to the low critical densities of the pure rotational transitions.

Future studies will focus on studying the effects of disc parameters (mass, size, flaring index, ...) on the CO ro-vibrational emissions and will include larger numbers of vibration levels and electronic excited states.

Acknowledgments

We thank ANR (contract ANR-07-BLAN-0221) and PNPS of CNRS/INSU, France for support. I.K., W.F.T. and P.W. acknowledge funding from the EU FP7-2011 under Grant Agreement nr. 284405. G.v.d.P. thanks the Millennium Science Initiative (ICM) of the Chilean ministry of Economy. Computations presented in this paper were performed at the Service Commun de Calcul Intensif de l’Observatoire de Grenoble (SCCI) on the super-computer funded by Agence Nationale pour la Recherche under contracts ANR-07-BLAN-0221, ANR-2010-JCJC-0504-01 and ANR-2010-JCJC-0501-01. We acknowledge discussions with Ch. Pinte, F. Ménard, and A. Carmona. We thank the referee and the editor for the useful comments that helped improve the paper.

References

Online material

Appendix A

Appendix A.1: Rotational diagram analysis for levels in LTE and optically thick lines

We adopted the rotational diagram analysis that takes optical depth effects into account (also called population diagram analysis in Goldsmith & Langer 1999). A similar approach was taken by Goto et al. (2011). The line intensity emitted by a slab at a single temperature reads (A.1)B_ν is the Planck function at the excitation T_ex = T_rot = T_vib. The second term in the right-hand side is an optical depth correction factor to the optically thin rotational diagram. The term β = (1 − exp( − τ_ν))/τ_ν is akin to the escape probability. The maximum optical depth occurs at E_rot ~ kT_rot, where E_rot is the rotational upper level energy of an emission line. In the Born-Oppenheimer approximation, E_u = E_rot + E_vib. The rotational energy for CO can be approximated by E_rot = B_rotJ_u(J_u + 1), with the rotational constant B_rot = 2.76 K. The maximum optical depth occurs when . The line optical depth is (A.2)where ν is the line frequency, A the Einstein spontaneous emission probability of the transition, c the speed of light, Δv the turbulent width, N(CO) the CO column density, i the disc inclination, and (A.3)is the CO fractional population in the initial level (v′,J′). For CO, Q_rot ≃ kT_rot/B_CO + 1/3 (B_CO = 2.76 K) and Q_vib = 1/(1 − exp( − 3122/T_vib) are the rotational and vibrational partition function respectively (Brittain et al. 2007).

The rotational diagram of AB Aur (Fig. A.1) shows three parts: a steep slope corresponding to increasing optical depth until E_rot ~ kT_ex, then a shallow slope where the line optical depth decreases because the higher the rotational level, the less they are populated. At very high J, the slope steepens again. This behaviour also appears in our theoretical rotational diagrams, but the second turning point occurs at higher J than in the observations. The second slope change corresponds to lines with τ < 1. Assuming that the population is in LTE, we fitted the AB Aur rotational diagram by a model that takes the optical depth effects into account (Goldsmith & Langer 1999). The model parameters are the CO column density N(CO) = 4.2 × 10¹⁷ cm^-2, the mean excitation temperature T_ex ≃ 600 K (we assume that the vibrational temperature is equal to the rotational temperature), the turbulent width Δv = 0.05, and the inclination i = 30 degree. In the upper panel of Fig. A.1, the analytical solid curve compares well with the observations. The dashed-line curve shows the same model where optical depth effects are not taken into account. The lower panel shows the derived optical depths, which reach ~47.

Appendix A.2: Excitation, radiation, and gas kinetic temperatures

We define the radiation (brightness) temperature T_rad at a given wavelength as the equivalent blackbody temperature that will match the specific intensity computed by the continuum radiative transfer at a given location in the disc (Figs. A.2 and A.3). We also show temperature ratios in Figs. A.4 and A.5: the excitation over the radiation and the excitation over the gas kinetic temperatures.

All Tables

All Figures