© ESO 2019

1 Introduction

Asymptotic giant branch (AGB) stars are among the most important contributors of dust and heavy elements in the universe. These stars enrich the interstellar medium (ISM) and galaxies by ejecting a large fraction of their material through strong stellar winds. An extended circumstellar envelope (CSE) will be created around the star as a consequence of the intense mass loss (e.g. Habing & Olofsson 2003). Understanding the complex chemical networks in the CSE of AGB stars is required for a better understanding of the enrichment and chemical evolution of the ISM and galaxies.

Carbon monoxide (CO), the most abundant molecule after molecular hydrogen (H₂), has been used to constrain the physical properties and chemical composition of the ISM and CSEs (e.g. Goldreich & Scoville 1976; Scalo & Slavsky 1980; Millar et al. 1987; Garrod & Herbst 2006; Morata & Herbst 2008). Photodissociation by ultraviolet (UV) radiation is the dominant process destroying CO and determining its abundance distribution. Therefore, a precise estimation of the CO photodissociation rate is important in both chemical and physical modelling of CSEs.

Generally speaking, molecular photodissociation by UV radiation can dominate by direct or indirect photodissociation, depending on the molecular structure. In direct photo-dissociation, the photodissociation cross section is continuous as a function of photon energy (continuum photodissociation). Thus, all absorptions lead to molecular dissociation. For indirect photodissociation, the photodissociation cross section contains a series of discrete peaks (line photodissociation). Therefore, absorption at only certain wavelengths leads to molecular dissociation. At wavelengths shorter than that of the H Lyman limit (911.7 Å), CO photodissociation occurs entirely in a set of discrete UV wavelength lines. This process makes CO strongly subject to self shielding, in cases of high abundance, and to mutual shielding by other species which are dissociated at the same wavelengths such as atomic and molecular hydrogen, atomic carbon, and dust. The dissociating photons will be absorbed by species closer to the UV source and thus molecules at the deeper regions will be shielded. The amount of shielding depends on the UV intensity, the geometry of the cloud which determines the photon penetration probability, and the column density of species with the same dissociating wavelengths. Bally & Langer (1982) and Glassgold et al. (1985) have shown the importance of the CO self-shielding in molecular clouds based on anomalous intensity ratios of various CO isotopologues which are selectively photodissociated in the edge of molecular clouds due to various column densities.

The most updated CO unshielded photodissociation rate in the interstellar medium of the solar neighborhood is estimated to be 2.6 × 10⁻¹⁰ s⁻¹ by Visser et al. (2009). We derive the same unshielded photodissociation rate in the outflows of evolved stars assuming that the Draine (1978) radiation field penetrates the envelope. This rate depends only on the following: the radiation field, the accuracy of the CO spectroscopic data and the surface temperature of the astrophysical region.

However, the geometrical distances over which the photodissociation is significant depends strongly on the physicalconditions of the penetrated environment such as its geometry, temperature distribution, and H₂ density. Thus, the environmental properties should be taken into account in calculations of the self- and mutual-shielding functions in various environment such as interstellar clouds and CSEs.

There is a good understanding of the depth dependency of the CO photodissociation mechanism in interstellar clouds and the photodissociation rate has been regularly updated with new laboratory data (e.g. Solomon & Klemperer 1972; Bally & Langer 1982; Glassgold et al. 1985; van Dishoeck & Black 1986, 1988; Visser et al. 2009). The shielding functions for several sets of input parameters for interstellar clouds can be downloaded¹ from the work by Visser et al. (2009).

In case of CSEs, Morris & Jura (1983, hereafter MJ83) present the theory of calculating the depth dependency of the CO photodissociation rate by considering CO self-shielding and H₂ mutual-shielding in a spherical expanding envelope through a “one-band approximation”. In this approximation, they assume that CO dissociates only at 1000 Å. After higher resolution laboratory measurement of far-UV absorption and fluorescence cross sections by Letzelter et al. (1987), the rate was updated by Mamon et al. (1988, hereafter MGH88) considering 34 dissociating-bands. Afterwards, Visser et al. (2009) collected the latest CO laboratory measurements of 855 UV dissociating lines arising in levels J = 0–9 in the lowest vibrational state v = 0. The updated CO photodissociation rate in a CSE is presented in two recent works by Li et al. (2014) and Groenewegen (2017). However, both works use the shielding functions calculated for an interstellar cloud (Visser et al. 2009) not a CSE. In the current work, we aim to calculate the CO photodissociation rate in a CSE using the most updated laboratory data and following the shielding functions developed for CSEs by MJ83 and MGH88. We discuss the differences of the depth dependency of the photodissociation rate in interstellar clouds and CSEs in Sect. 4.3. In addition, we have investigated the effect of several more free parameters on the CO dissociation rate. The new treatment of the CO photodissociation has been incorporated into a chemical network describing the CSEs of AGB stars.

2 Circumstellar chemistry

We use an extended version of the publicly available circumstellar envelope chemical model rate13-cse code² (McElroy et al. 2013). We modified calculations of the CO photodissociation rate from a one-band approximation to a treatment where all known lines are taken into consideration (see Sect. 2.2).

The code assumes a spherically symmetric envelope which is formed due to a constant mass loss Ṁ. The envelope expands with a constant radial velocity V_exp. The H₂ density falls as 1∕r² where r measures the distance from the central star. We used the gas temperature profile given by MGH88 which is derived for CW Leo as follows: $$T_{\text{kin}}(r)=14.6{\left(\frac{r_0}{r}\right)}^{\beta }(\text{K}),$$(1)

where r₀ = 9 × 10¹⁶ cm, β = 0.72 for r < r₀ and β = 0.54 for r > r₀. We assume a minimum temperature of 10 K in the outer CSE.

We incorporated an extended version of the chemical network from the UMIST Database for Astrochemistry (McElroy et al. 2013) which includes the ¹³C and ¹⁸O isotopes, all corresponding isotopologues, their chemical reactions and the properly scaled reaction rate coefficients (Röllig & Ossenkopf 2013). The chemical network includes 933 species and 15108 gas-phase reactions. The isotopologue chemistry will be discussed in detail in a forthcoming paper.

2.1 CO spectroscopic data

The dissociation energy of the CO ground state is 11.09 eV, thus the photodissociation occurs in the wavelength range 911.75–1117.8 Å. The CO photodissociation occurs entirely through line absorptions into pre-dissociating states (e.g. Visser et al. 2009). We use the latest laboratory measurements of CO data (Visser et al. 2009) which includes 855 UV transitions containing rotational excitations J = 0−9 to calculate the photodissociation rate. The higher excitation levels up to J = 30 were examined in the calculations and since the effect was marginal, we excluded the higher transitions to increase the time efficiency of the code.

2.2 CO photodissociation rate

The total CO photodissociation rate at radius r is the summation of the photodissociation rates of all discrete contributing lines i as follows: $$k(r)={\displaystyle \sum _{i=1}^{855}}{k}_i^0{\beta }_i(r){\gamma }_i(r)({\text{s}}^{-1}),$$(2)

where $k_i^0$ is the CO unshielded photodissociation rate at the edge of the cloud, β_i counts the CO self-shielding efficiency and γ_i counts the mutual shielding from other species. $k_i^0$ is estimated to be 2.6 × 10⁻¹⁰ s⁻¹ and the detailed calculations are presented in the Appendix A.1. We note that in calculations of the photodissociation rates we assumed the CO excitation temperature T_ex to be the same as the gas kinetic temperature T_kin given in Eq. (1).

Dissociating radiation can be absorbed by CO (self-shielding), H, C, H₂ and dust (mutual-shielding) (e.g. MGH88; Visser et al. 2009). The mutual shielding by different species depends on their column density and the amount of line overlaps in the relevant region of the spectrum. The amount of mutual shielding by dust is assumed to depend on the total number of protons [n(H)+2n(H₂)] which determine the dust extinction. In our models, the column densities of C and H are insufficient to produce very much blocking. Thus, the shielding by dust and H₂ dominates the CO mutual shielding. We note that in the presence of extra UV radiation from a hot binary companion and/or stellar chromospheric activity there could be an enhancement of the atomic abundances. The investigation of how and whether these would impact the CO abundance distribution is, however, beyond the scope of the current study.

Calculations of the depth dependency of the CO photodissociation rate require accurate information on the line wavelength, oscillator strengths, the pre-dissociation probabilities and the line widths of CO and H₂. We used the new compiled H₂ line list (J. Black, priv. comm.) based on energy levels and transition probabilities computed by Abgrall et al. (1994, 1997). In the Appendix A, we review the underlying physics of calculations of shielding functions β and γ in an expanding CSE.

Fig. 1 CO fractional abundance distributions for simulations with different shielding functions that regulate CO photodissociation: shielding by H2, shielding by dust, CO self-shielding, and the total shielding for the reference model.

3 Results

To study the significance of the different shielding processes we use a reference model with the physical parameters listed in Table 1. This model serves as a direct comparison to the work by MGH88. Figure 1 shows the fractional CO abundance profile calculated by considering different shielding contributions from CO, H₂, dust, and the total shielding for the reference model. As we can see, the CO self-shielding plays a major role in the total shielding.

In Fig. 2 we illustrate how the CO abundance distribution varies for models with different initial CO abundances f₀ and mass-loss rates Ṁ while preserving the total amount of CO ejected by the star. It is clear that the H₂ density (set by Ṁ and V_exp) dictates the size of the envelope, while the initial CO abundance affects the steepness of the slope. This differentiation provides a way to break the degeneracy between f_CO and Ṁ encountered when both low-J and high-J CO transitions are used in CO radiative transfer (RT) modelling.

Figure 3 presents the variation of the CO distribution profiles with Ṁ, keeping all other parameters the same as the reference model. An increase in Ṁ translates into a stronger shielding and hence a larger CO envelope with a sharper drop-off at the outer edge.

The CO abundance profiles derived from our simulations can be fitted by the analytical formula derived by MGH88 as follows: $$f_{\text{CO}}=f_0\text{exp}\left(-\text{ln}(2){\left(\frac{r}{r_{1/2}}\right)}^{\alpha }\right),$$(3)

where f₀ is the initial CO abundance, α determines the steepness of the profile and r_1∕2 marks the radius where the CO abundance drops to half of its initial value. Figure 4 shows the accuracy of the fitting formula for a range of mass-loss rates; all other parameters are the same as the reference model given in Table 1. In general, f₀ depends on the chemical type of AGB star and is commonly assumed to have average values (2 − 6 − 10) × 10⁻⁴ for M-, S-, and C-type AGB stars, respectively (Ramstedt & Olofsson 2014).

We derive α and r_1∕2 for a grid containing 390 models with varying f₀, Ṁ, and V_exp. We considered ten values for f₀ (1, 2, ⋯, 10 × 10⁻⁴), 13 values for Ṁ ([1, 2, 5] ×10⁻⁸, [1, 2, 5] × 10⁻⁷, …, 1 × 10⁻⁴ M_⊙ yr⁻¹), and three values for V_exp (7.5, 15, 30 km s⁻¹). Table B.1 lists the resulting α and r_1∕2 values for all models. Figure 5 shows how variations of f_CO, Ṁ, and V_exp affect α and r_1∕2. Increasing Ṁ and decreasing V_exp leads to an enhancement of the $n_{{\text{H}}_2}$ and thus more shielding and a larger CO envelope. Similarly an increase in f_CO enhances the CO self-shielding. The most drastic changes in α and r_1∕2 happen for oxygen-type AGB stars with lower CO abundance, where the CO self-shielding is not very efficient.

Fig. 2 CO fractional abundance distributions for models with a constant expansion velocity Vexp = 15 (km s−1) and various Ṁ (M⊙ yr−1) and fCO in a way to keep the same amount of effective CO for all models.

Fig. 3 CO fractional abundance distributions for models with fCO = 8 × 10−4, Vexp = 15 (km s−1), and a range of mass-loss rates which are marked in the figure.

Fig. 4 Comparison of the CO fractional abundance distributions from modelling (solid black lines) and the analytic fitting formula (dashed blue lines) for models with fCO = 8 × 10−4, Vexp = 15 (km s−1) and a range of 10−8 < Ṁ < 10−4 (M⊙ yr−1).

4 Discussion

4.1 Influence of the temperature on the CO abundance distribution

We examined the influence of the gas kinetic temperature profile on the CO envelope size for the reference model. The CO excitation temperature is assumed to be the same as the gas kinetic temperature profile which is given in Eq. (1). We assumed the temperature and radius at the inner envelope to be T₀ = 2000 K and r₀ = 10¹⁴ cm. We considered 0.4 < β < 1.0 which reasonably covers the gas temperature profile of AGB CSEs (De Beck et al. 2012; Danilovich et al. 2014; Khouri et al. 2014; Maercker et al. 2016; Ramos-Medina et al. 2018; Van de Sande et al. 2018). β determines the slope of the profile. Figure 6 shows the considered temperature profiles and their corresponding CO abundance distributions. The fitting parameters vary in the ranges 1.65 × 10¹⁷ < r_1∕2 < 2.03 × 10¹⁷ cm and 3.00 < α < 3.24. The temperature profiles with β ≥ 0.6 give rise to similar CO abundance distributions whereas those with β < 0.6 lead to smaller envelopes. This is due to a lower shielding effectiveness in the hotter envelope (see Appendix A).

We also examined the influence of constant CO excitation temperature T_ex = 5, 20, 50, 100 K as assumed by Groenewegen (2017) versus T_ex = T_kin for the reference model. The two fitting parameters vary in ranges 1.53 × 10¹⁷ < r_1∕2 < 2.25 × 10¹⁷ cm and 2.88 < α < 3.41. As shown in Fig. 7 a lower CO excitation temperature results in a bigger CO envelope. The reduced T_ex leads to higher lower-level populations x_l (see Appendix A) and thus more efficient shielding and therefore a bigger envelope.

4.2 Influence of the ISRF on the CO abundance distribution

The Draine interstellar radiation field (ISRF) that we considered in this work has been measured in the solar neighborhood. McDonald & Zijlstra (2015) show that the ISRF is considerably higher in globular clusters. This will have a significant impact on the size of the CO envelope of evolved stars which are located in clusters (McDonald et al. 2015). On the other hand, objects which lie above the Galactic plane, for example IRC+10216, are possibly exposed to a weaker ISRF. We investigated the influence of the strength of the ISRF on the CO envelope size. We scaled the Draine (1978) radiation field that we used for the reference model by factors χ = 0.25, 0.5, 1, 2, and 4. Table 2 presents the results for models with the mass-loss rates Ṁ = 10⁻⁸, 10⁻⁷, 10⁻⁶, 10⁻⁵, 10⁻⁴ M_⊙ yr⁻¹. We considered the expansion velocity V_exp = 15 km s⁻¹ and the initial CO abundance of f_CO = 8 × 10⁻⁸ for all models. As expected, the effect of the ISRF strength on the CO envelope size is more prominent in stars with low mass-loss rates with weaker shielding efficiency. Increasing the ISRF by a factor of two reduces the CO envelope size by 40, 34, 28, 20, 14% for stars with Ṁ = 10⁻⁸, 10⁻⁷, 10⁻⁶, 10⁻⁵, 10⁻⁴ M_⊙ yr⁻¹, respectively.

We notethat in some cases, there is extra UV radiation which internally penetrates into the CSE. The inner UV radiation can arises from stellar chromospheric activity and/or a hot binary companion (e.g. Montez et al. 2017). In hot post-AGB stars, the high temperature of the star itself also generates UV photons. In such cases, the inner UV radiation will be quickly absorbed by inner dust and thus is not expected to affect the extent of the CO envelope. However this likely affects the CSE chemistry (e.g. Saberi et al. 2017, 2018; Van de Sande & Millar 2019) in the dust formation region, which is beyond the scope of this paper.

4.3 Environmental dependency of the CO photodissociation rate

The differences in the physical and chemical properties between interstellar clouds and CSEs can affect the CO shielding functions and thus the depth dependency of CO photodissociation rates. Here we list the environmental differences between interstellar clouds (Visser et al. 2009) and CSEs (this work) that enter in calculations of shielding functions:

• (a)

  Model geometry: in interstellar cloud models, plane-parallel geometry has been considered while in the CSE model spherically symmetric geometry is assumed. The difference in the geometry affects the photon penetration probability.

• (b)

  CO excitation temperature: a constant temperature that can be chosen among the values of 5, 20, 50 and 100 K is available in interstellar clouds model. In the modelling of the CSE, we consider T_ex (CO) to be the same as T_kin which varies in a range ~10−2000 K. This affects the fractional population of the lower-level especially in the inner CSE. Assuming a too low T_ex leads to an overestimation of the CO self-shielding as shown in Fig. 7.

• (c)

  Atomic and molecular line broadening: the Doppler and natural broadenings are the dominant processes which control the CO, H₂ and H line widths in the interstellar clouds. Since the CSE around AGB stars is expanding at velocities of typically a few up tosome 30 km s⁻¹, the expansion velocity should also be considered in the calculations of the line widths.

Fig. 5 Variation of two fitting parameters α and r1∕2 for all models presented in Table B.1.

Fig. 6 Left panel: tested gas kinetic temperature profiles with different β values. Theprofile with label “reference” corresponds to the temperature of the reference model. Right panel: CO abundance distribution profile derived using the temperature profiles from the left panel.

Fig. 7 CO fractional abundances distributions for the reference model when different CO excitation temperature profiles are considered in calculations of the CO photodissociation rate.

4.4 Comparison with literature

The most commonly used method to derive the size of the CO envelope is the one presented by MGH88. They used the Jura (1974) radiation field to calculate the CO photodissociation rate. The Jura radiation field in addition to more UV observations at longer wavelengths up to 2000 Å has been later used by Draine (1978) to derive an analytical formula for the standard ISRF (Lee 1984). Thus, in principle in the wavelength range 930–1125 Å, both Draine and Jura radiation fields represent the same UV observational data. MGH88 presented the fitting parameters of the CO envelope size for a grid of models with a constant CO abundance of $f_{{\text{CO/H}}_2}=8\times {10}^{-8}$ and 13 varying mass-loss rates and three expansion velocities. However, in addition to now outdated low-resolution CO laboratory data which causes underestimation of the unshielded photodissociation rate by 30%, a major drawback of this work is that it is based on calculations that assume one fixed value of f₀ = 8 × 10⁻⁴. We clearly demonstrate in the previous section that the role of f₀ is non-negligible in the overall dissociation efficiency. Table 3 and Fig. 8 compare the MGH88 results for the r_1∕2 and α fitting parameters with those derived in this work. The difference in fitting parameters ranges from 0.6 to 15% for α and 1 to 39% for r_1∕2 between two works. In almost all tested cases, our models predict the steepness parameter α to be larger than that derived by MGH88, indicating a less efficient shielding of CO in our study. In line with this, we calculate consistently smaller values for r_1∕2. This is consistent with observational data for W Hya (Khouri et al. 2014), TX Cam (Ramstedt et al. 2008), and R Dor (Maercker et al. 2016) for example. Models of these objects reproduce the observed line emission better when they assume a smaller CO envelope than is predicted by MGH88, in line with the results of this paper.

Our derived CO radius is also smaller by 11–60% compared to Groenewegen (2017) depending on the H₂ density of the envelope. This comes from different shielding functions used to calculate the depth-dependent CO photodissociation rate.

4.5 Predicting the CO line fluxes

We have tested the effect of CO abundance profiles derived in this work and by MGH88 on the CO line intensities for two models in Table 3. We use a non-local thermodynamic equilibrium (non-LTE) RT code based on the Monte Carlo programme (see e.g. Schöier & Olofsson 2001) for the excitation analysis.

Model 1. The first model is the reference model with its CSE properties are given in Table 1. We assumed that the reference star is at a distance of 1 kpc in order to cover the entire envelope with a 12-m telescope beam. The fitting parameters of the CO envelope size are r_1∕2 = 2.02 × 10¹⁷ cm and α = 3.10 from this work and r_1∕2 = 2.35 × 10¹⁷ cm and α = 2.79 from MGH88. Figure 9 presents the results of RT modelling for CO(1–0, 2–1, 3–2) transitions for this model. The 14% difference in r_1∕2 size only becomes visible in J = 1−0 spectrum in this case. We find a 5% difference in the total intensity of J = 1−0 transition which is less than the commonly assumed 20% uncertainty of single-dish observations. However, comparing the integrated intensity at radial offset points from the central star shows bigger discrepancy between two models. The differences in the integrated intensity at distance 22″ = 3.2 × 10¹⁷ cm reaches to 34%. This indicates that CO high-resolution observations which provide the integrated intensity at radial offset points are more powerful to constrain photodissociation models.

Model 2. We selected the model with the biggest discrepancy with MGH88 work. This model has Ṁ = 10⁻⁴ M_⊙ yr⁻¹ and V_exp = 7.5 km s⁻¹. The star is assumed to be at distance of 3 kpc from the earth to cover the entire envelope with a 12-m telescope beam. The fitting parameters of the CO envelope size are r_1∕2 = 9.91 × 10¹⁷ cm and α = 3.80 from this work and r_1∕2 = 1.64 × 10¹⁸ cm and α = 3.71 from MGH88. Figure 10 shows the results of the CO RT modelling for these two models. The clearest difference is again apparent for CO(1–0), where the emission is substantially different in the velocity range [−10,−5] km s⁻¹. In the blue model, the emission is entirely self-absorbed. This can be explained by the absence of a significant amount of cold CO gas in the red model, which is characterised by a smaller r_1∕2 than the blue model. For this model, the integrated intensities at radial offset points of the red model (this work) are higher due to the absorption in the blue model (MGH88).

5 Summary

We have presented detailed calculations of the CO photodissociation rate in a spherically symmetric CSE which is expanding with a constant velocity. The standard Draine (1978) radiation field is assumed to penetrate into the CSE from all directions. We used the latest CO spectroscopic data to calculate the shielding functions. We examined the impact of variation of five primary important factors Ṁ, V_exp, f_CO, T_ex (CO), and the strength of the ISRF χ on the CO abundance distributions. The effect of varying parameters on the CO envelope size is more prominent for either low-mass loss stars or the ones with low initial CO abundance. This can be explained by lower shielding efficiency.

Assuming the same ISRF and CSE properties, our derived CO envelope size is smaller than the commonly used radius presented by MGH88. We show that having the same amount of effective CO with a different set of physical parameters does not necessarily give the same CO abundance distribution. High-resolution ALMA observations, for example the DEATHSTAR³ project (Ramstedt et al. In prep), can, together with our new formalism of determination of the CO envelope size, be used to decrease theuncertainty in mass-loss determinations.

Although we listed two fitting parameters of the CO abundance distribution for a large grid of models, it is recommended to run models for individual stars separately considering their individual physical parameters. Moreover, we strongly recommend running optimised models in case there are clear indications for a locally weak or strong ISRF based on other observations.

Fig. 8 Comparison of CO abundance distribution parameters to MGH88 results. The left panel shows $r_{1/2}/r_{1/2}^{*}$ ratio and theright panel shows α∕α* with $r_{1/2}^{*}$ and α* are the MGH88 values.

Fig. 9 Results of CO RT modelling for model 1 (the reference model) with the CO abundance distribution from this work (red) and MGH88 (blue). The transitions and the beam size are marked in each panel. Bottom right panel: compares the radial distribution of the CO(1–0) integrated intensities at radial offset points.

Fig. 10 Results of CO RT modelling for model 2 with extreme discrepancy between the CO abundance distribution from this work (red) and MGH88 (blue). The transitions and the beam size are marked in each panel. The bottom right panel compares the radial distribution of the CO(1–0) integrated intensities at radial offset points.

Acknowledgements

This work was supported by ERC consolidator grant 614264. E.D.B. acknowledges financial support from the Swedish National Space Agency. We thank John Black for fruitful discussions on the detailed calculations of photodissociation rates, providing the H₂ data list, and commenting on the manuscript, Tom Millar for several discussions and explanations of the original circumstellar chemistry code and commenting on the manuscript, Theo Khouri for useful discussions, Evelyne Roueff and Franck Le Petit for useful discussions on PDR models. We are grateful to the anonymous referee for their constructive feedback on the manuscript.

Appendix A CO shielding functions

A.1 Unshielded CO photodissociation rate

The contribution of each individual transition to the unshielded photodissociation rate at the outer radius can be calculated using (van Dishoeck & Black 1988): $$k_i^0=\frac{\pi e^2}{m_e{c}^2}{f}_i{\eta }_i{x}_{\text{l}}{\lambda }_i^2{I}_{\text{ISRF}}({\lambda }_i)(s^{-1}),$$(A.1)

where f is the absorption oscillator strength that expresses the probability of absorption of electromagnetic radiation in transitions between energy levels of a molecule, η is the probability for dissociation of the upper level, x_l is the fractional populations of the lower level and I_ISRF is the mean intensity of the ISRF in unit (photons cm⁻² s⁻¹ Å⁻¹). The constant factor πe²∕m_ec² takes the value 8.85 × 10⁻²¹ if λ is in Å. We have assumed x_l has a Boltzmann distribution profile as follow: $$x_{\text{l}}=\frac{g_{\text{l}}\exp (-E_{\text{l}}/k{T}_{\text{ex}})}{{\sum }_{l=0}^9{g}_{\text{l}}\exp (-E_{\text{l}}/k{T}_{\text{ex}})},$$(A.2)

where g_l and E_l are the degeneracy and energy of the lower level. Since in our model T_ex (CO) = T_kin, to derive the fractional population of lower levels at the edge of the CSE, we considered T_ex = 10 K which is the minimum acceptable temperature in our chemical network.

We considered the standard Draine radiation field Draine (1978) as the ISRF which has the form: $$I_{\text{ISRF}}(\lambda )=\left(\frac{6.36\times {10}^7}{{\lambda }^4}-\frac{1.0237\times {10}^{11}}{{\lambda }^5}+\frac{4.0812\times {10}^{13}}{{\lambda }^6}\right),$$(A.3)

here I is in unit (ergs cm${}^{-2}{\text{s}}^{-1}{Å}^{-1}$). We multiplied I by a factor (5.03 × 10⁷ × λ) to convert the intensity unit to (photons cm⁻² s⁻¹ Å⁻¹) which is needed in Eq. (A.1).

A.2 CO self-shielding

The CO self-shielding in an expanding spherical envelope was approximated by MJ83 to be: $${\beta }_i(r)=\frac{1-\exp (-1.5{\tau }_{\text{CO}}({\nu }_i,r))}{1.5{\tau }_{\text{CO}}({\nu }_i,r)},$$(A.4)

where β_i(r) indicates the escape probability for a photon which is generated at radius r. When the absorption occurs primarily in the Doppler core of the line and by assuming that the CSE expansion velocity V_exp dominates the line width, the optical depth at the centre of each dissociating line ν_i can be written: $${\tau }_{\text{CO}}({\nu }_i,r)=\frac{\pi e^2}{m_{e}c}{N}_{\text{CO}}(r)x_l(r)f_{ij}{\lambda }_i\frac{1}{V_{\text{exp}}}+{\displaystyle \sum _j}{\tau }_{\text{CO}}({\nu }_j,r),$$(A.5)

where N_CO is a column density integrated from outside to r. The first term counts the CO optical depth at the centre of each dissociating line i and the second term counts the effect of all blended lines j if λ_i− λ_j < Δλ_j. Here, Δλ_j is the CO line broadening. In the first term πe²∕m_ec = 0.0265 and λ is in cm and V_exp is in cm s⁻¹. The fractional population of the lower level x_l varies by the radius here.

To derive τ_CO(ν_j, r) we have to estimate the line width of each dissociating line. Since the effect of the expansion velocity is equivalent for all lines, we need to onlyconsider the thermal and natural broadenings.

Thermal broadening

The thermal or Doppler broadening due to the random thermal motions of molecules has a gaussian profile as follows: $${\varphi }_{\nu }=\frac{1}{\sqrt{\pi }}\frac{1}{\Delta {\nu }_{\text{th}}}\text{exp}\left\{-{(\frac{\nu -{\nu }_{\text{0}}}{\Delta {\nu }_{\text{th}}})}^{\text{2}}\right\},$$(A.6)

where the line of sight thermal width is defined as $$v_{\text{th}}\equiv \sqrt{\frac{2k{T}_{\text{k}}}{m_{\text{CO}}}},$$(A.7)

and in frequency units: $$\Delta {\nu }_{\text{th}}=\frac{{\nu }_0}{c}\sqrt{\frac{2k{T}_{\text{k}}}{m_{\text{CO}}}}.$$(A.8)

The peak value of ϕ(ν = ν₀) is given by $${\varphi }_0=\frac{1}{\sqrt{\pi }}\frac{1}{\Delta {\nu }_{\text{th}}},$$(A.9)

and the full-width at half maximum (FWHM): $$\Delta v_{1/2}=1.665v_{\text{th}},$$(A.10)

so the CO thermal broadening across the temperature range 10–2000 K is equivalent to 0.13–1.8 km s⁻¹. These give the marginal line broadenings of Δλ_1∕2 = (Δv_1∕2∕c) λ ~ 0.4 × 10⁻³−0.6 × 10⁻² at the wavelength of 1000. Thus, we can ignore the thermal line broadening of CO.

Natural broadening

The natural line broadening resulting from the Heisenberg uncertainty principle gives a Lorentzian profile: $${\varphi }_{\nu }=\frac{\Gamma }{4{\pi }^2{(\nu -{\nu }_0)}^2+{(\Gamma /2)}^2},$$(A.11)

where Γ is the quantum-mechanical damping constant and represents the total radiative decay probability or the inverse radiative lifetime of the upper level in s⁻¹. The peak value of ϕ(ν = ν₀) at the line centre is given by $${\varphi }_0=\frac{4}{\Gamma },$$(A.12)

and the FWHM in frequency units will be: $$\Delta {\nu }_{1/2}=\frac{\Gamma }{2\pi }.$$(A.13)

CO has Γ ~ 10⁸−10¹³ s⁻¹ which leads to the line broadening Δλ_1∕2 = (λ²∕C)Δν_1∕2 ~ (0.5 × 10⁻⁵−0.5). Thus, for the lines with large Γ the natural broadenings are significant.

As we now have the precise line positions and line widths of each individual line, we can estimate the number and contribution of blended lines in the opacity of each dissociating lines. We assumed that if λ_i − λ_j < Δλ_j, the j line is be considered in the second term of Eq. (A.5).

A.3 Mutual-shielding by H₂ and dust

The shielding of CO by molecular hydrogen, H₂, and dust in an expanding CSE is numerically approximated by MJ83 to be: $${\gamma }_i(r)=\exp (-\alpha {({\tau }_{\text{dust}}(r)+{\tau }_{{\text{H}}_2}({\nu }_i,r))}^b),$$(A.14)

where α = 1.644, b = 0.86. ${\tau }_{{\text{H}}_2}$ is the H₂ opacity at each CO dissociating line i and can be calculated as: $${\tau }_{{\text{H}}_2}({\nu }_i,r)=\frac{\pi e^2}{m_{e}c}{f}_m{\varphi }_{{\nu }_i}{N}_{{\text{H}}_2},$$(A.15)

here ${\varphi }_{{\nu }_i}$ is the H₂ line shape at the CO line i. Since the H₂ lines become optically thick, the absorption occurs in the radiative wings of the Lyman and Werner transitions. Thus, CO mutual-shielding by H₂ mostly occurs at the Lorentzian damping wings of the H₂ line profile (MJ83, MGH88) which is presented in Eq. (A.11). This gives the

H₂ opacity at each CO dissociating line i: $${\tau }_{{\text{H}}_2}({\nu }_i,r)={\displaystyle \sum _m}\frac{\pi e^2}{m_{e}c}{f}_m{N}_{{\text{H}}_2}\frac{{\Gamma }_m}{4{\pi }^2{({\nu }_m-{\nu }_i)}^2+{({\Gamma }_m/2)}^2},$$(A.16)

where m sums over all H₂ dissociating lines.

In Eq. (A.14), we assumed that the dust absorption is independent of the wavelength in the spectral region of interest. We also ignored dust scattering and assume that the dust absorption dominates the dust extinction. Therefore, we considered a constant dust extinction at 1000, as MJ83, to be: $${\tau }_{\text{dust}}(r,1000\text{Å})=\frac{4.65\times 2\times N_{{\text{H}}_2}(r)}{1.87\times {10}^{21}},$$(A.17)

where $N_{{\text{H}}_2}$ is the H₂ column density to infinity.

Appendix B Grid models

Table B.1 presents the two fitting parameters of the CO envelope size for the full grid of models.

References

All Tables

All Figures

	Fig. 1 CO fractional abundance distributions for simulations with different shielding functions that regulate CO photodissociation: shielding by H2, shielding by dust, CO self-shielding, and the total shielding for the reference model.
In the text

	Fig. 2 CO fractional abundance distributions for models with a constant expansion velocity Vexp = 15 (km s−1) and various Ṁ (M⊙ yr−1) and fCO in a way to keep the same amount of effective CO for all models.
In the text

	Fig. 3 CO fractional abundance distributions for models with fCO = 8 × 10−4, Vexp = 15 (km s−1), and a range of mass-loss rates which are marked in the figure.
In the text

	Fig. 4 Comparison of the CO fractional abundance distributions from modelling (solid black lines) and the analytic fitting formula (dashed blue lines) for models with fCO = 8 × 10−4, Vexp = 15 (km s−1) and a range of 10−8 < Ṁ < 10−4 (M⊙ yr−1).
In the text

	Fig. 5 Variation of two fitting parameters α and r1∕2 for all models presented in Table B.1.
In the text

	Fig. 6 Left panel: tested gas kinetic temperature profiles with different β values. Theprofile with label “reference” corresponds to the temperature of the reference model. Right panel: CO abundance distribution profile derived using the temperature profiles from the left panel.
In the text

	Fig. 7 CO fractional abundances distributions for the reference model when different CO excitation temperature profiles are considered in calculations of the CO photodissociation rate.
In the text

	Fig. 8 Comparison of CO abundance distribution parameters to MGH88 results. The left panel shows $r_{1/2}/r_{1/2}^{*}$ ratio and theright panel shows α∕α* with $r_{1/2}^{*}$ and α* are the MGH88 values.
In the text

	Fig. 9 Results of CO RT modelling for model 1 (the reference model) with the CO abundance distribution from this work (red) and MGH88 (blue). The transitions and the beam size are marked in each panel. Bottom right panel: compares the radial distribution of the CO(1–0) integrated intensities at radial offset points.
In the text

	Fig. 10 Results of CO RT modelling for model 2 with extreme discrepancy between the CO abundance distribution from this work (red) and MGH88 (blue). The transitions and the beam size are marked in each panel. The bottom right panel compares the radial distribution of the CO(1–0) integrated intensities at radial offset points.
In the text