© ESO 2021

1 Introduction

Essentially all low and intermediate mass stars end their lives on the asymptotic giant branch (AGB), where they lose their envelope (of the order ~0.5 to 6.5 M_⊙ depending on the initial mass) that is returned to the interstellar medium (ISM) in the form of dust and gas. Determining the mass-loss rate (MLR) is of key importance in understanding AGB evolution in detail. One of the main methods to do so is to use carbon monoxide (CO) as it is the most abundant molecule after molecular hydrogen and it is very stable to photodissociation (see Höfner & Olofsson 2018 for an overview). As CO also essentially binds all carbon when C/O < 1 (oxygen-rich stars, or O-stars) and all oxygen when C/O > 1 (carbon-rich stars, or C-stars) the molecule is also important for the chemistry in the circumstellar envelopes (CSEs) around AGB stars, see e.g. Agúndez et al. (2020).

When performing detailed radiative transfer modelling of CO-lines (typically fitting multiple transitions) to determine the MLR (see e.g. Sahai 1990; Groenewegen 1994; Schöier & Olofsson 2001; Decin et al. 2006), a CO abundance profile is required as input. The CO density is usually parameterised as (1)

where f_CO is the photospheric CO abundance, r is the distance to the star, is the distance where the CO abundance is half the photospheric value (and typically called the photodissociation radius), and α describes the shape of the profile.

Recently, Groenewegen (2017) and Saberi et al. (2019) presented results of photodissociation calculations of CO (i.e. values for and α) for a large grid of MLRs, expansion velocities, f_CO values, and strengths of the interstellar radiation field (ISRF). Both papers were updates of Mamon et al. (1988) which was the classical reference for the CO photodissociation radii in circumstellar envelopes for almost thirty years. Both papers used different approaches and the CO photodissociation radii found by Saberi et al. (2019) are 11–60% smaller than found by Groenewegen (2017) (Saberi et al. 2019).

Groenewegen (2017) used the approach introduced in Li et al. (2014, 2016). It is based on the shielding functions of Visser et al. (2009), which depend on the column density of CO and molecular hydrogen (H₂)¹ as well as a numerical integration scheme that takes into account the fact that at any point in the wind UV photons can arrive from all directions (see Appendices A in Li et al. 2014 and Groenewegen 2017). The shielding functions from Visser et al. (2009) were derived for a static environment however.

Saberi et al. (2019) used the escape-probability formalism developed for an expanding envelope by Morris & Jura (1983; herafter MJ83) and later updated by Mamon et al. (1988). They used the latest updated molecular data for CO from Visser et al. (2009) and for H₂ to calculate the CO shielding functions. Both models assume spherical symmetry, a smooth wind (i.e. no clumps), a constant velocity in the outflow, and simplified scattering and absorption properties of the dust in the 912–1076 Å region.

In the escape-probability formalism the relevant integrations over angle are replaced by numerical approximations developed in MJ83. These approximations have been derived for a certain parameter space of line and continuum optical depths, but neither Mamon et al. (1988) nor Saberi et al. (2019) checked whether these conditions are in fact met.

In this paper, the validity of these approximations is investigated and an improved formalism is proposed. The mathematical problem and numerical solution to this problem are presented in Sect. 2. The calculations are presented in Sect. 3 and discussed in Sect. 4, while Sect. 5 concludes this paper.

2 Equations and solutions

Two equations are relevant in the description of the photodissociation process of CO. The first is a term describing the continuum absorption of radiation, (2)

where ϕ is the angle with respect to the radius vector from the central star (MJ83) and its numerical approximation, (3)

with an accuracy of 1% for τ_c < 5 (Eq. (5) in MJ83), where τ_c is the continuum optical depth measured radially outwards from a point at distance r from the central star. In MJ83, only the continuum absorption by dust was considered (τ_d), but also the contribution in the line wings of H₂ transitions should be included, (see the discussion in Mamon et al. 1988 and Eq. (A14) in Saberi et al. 2019).

The second relevant equation describes the escape probability for a line photon, (4)

which is approximated as (5)

with α = 1.5, where Δ is a small correction without which this expression for β is accurate to within 10% (Eq. (8) in MJ83). When (6)

Equation (5) is accurate to <0.5% if τ_d = 0 and a few percent when τ_d ≠ 0 (Eq. (9) in MJ83). Again, the line optical depth τ_i was measured outwards from a point in the envelope at distance r.

It should be pointed out that Mamon et al. (1988) and Saberi et al. (2019) ignore the Δ term (see the Appendix in Mamon et al. 1988 and Eq. (A4) in Saberi et al. 2019). These two papers also use a different notation in the sense that the term between square brackets is denoted with β, corresponding to the CO self-shielding, and γ is used for CO-mutual shielding by H₂ and dust.

It is remarked that Eqs. (2) and (4) ignore the finite size of the central star. Assuming that all radiation from behindthe central star is blocked, and that the UV emission of the central star itself may be ignored, the integration over ϕ effectively runs from 0 to (π − β), where (7)

indicates the angle subtended by the central star (or the inner radius of the envelope) from a point at distance r (see Eq. (5) and Fig. A1 in Groenewegen 2017). This effect was taken into account in Groenewegen (2017), while the approximations in MJ83 assume that the central star is a point source.

The results in the next section have been obtained using routines written in Fortran77 from Press et al. (1992): a Romberg integration schema (qromb) to perform the numerical integrations² and a non-linear least-squares fitting routine (a Levenberg–Marquardt algorithm, mrqmin) to derive the coefficients in Eqs. (3), (5), and (6). The precision of the numerical integration routine has been verified by comparing selected results to those obtained using an online tool for such calculations³.

3 Results

3.1 Calculation of γ

Figure 1 and Table 1 contain the results of the calculation of γ. The upper two panels illustrate the results for τ_c < 5. Panel a shows the exact calculation of γ against τ_c. The approximation of Eq. (3) is also plotted, but they are indistinguishable. Panel b shows the ratio of the exact calculation to that approximation. The maximum deviation is ≈ 1.5%, which is larger than claimed in MJ83. One observes that in this range in optical depth that the approximation is systematically lower (by ≈ 0.5%) than the exact calculation. For larger optical depths the deviations become increasingly larger, up to a factor of two at τ_c ≈ 16 (Table 1).

As the numerical integration scheme is compared to a fitting formulae, there is no observational error as is to be used in a classical χ² analysis. Instead, an error is assigned such that the final fitting formulae (Eqs. (11) and (12)) have a reduced χ² of unity. This allows one to monitor the reduced χ² in the several fitting steps and gives representative error bars in the fitting coefficients.

Comparing in a chi-square sense Eq. (3) (with fixed coefficients − 1.644 and 0.86) to the exact calculations over the range 0 < τ_c < 5 in steps of 0.025 units results in a reduced χ² of 7.3 if an “error” of 0.35% (see below) is assigned to each data point. Fitting for the coefficients and finding the solution with the smallest maximum deviation over the largest possible range in τ_c results in the best fit of: (8)

made in the range τ_c < 5.3. The maximum deviations are ±1.2% and the reduced χ² becomes 2.0, which is a significant reduction, indicating that the fit is much improved although the coefficients are similar. Panel c is arepeat of panel b for this approximation, showing that the deviations are now more symmetric around unity.

3.2 Calculation of β

Figure 2 compares the exact numerical results for the calculation of (β∕γ) in the range τ_i < 5 and τ_c < 5, to the approximation in Eq. (5) with Δ = 0 and γ approximated by Eq. (3). The differences are much larger than the claimed 10%, and almost 25% near τ_i = 2.5 (and for τ_c = 5). Deviations of less than 10% are only reached in a very limited range of τ_i and τ_c. Technically speaking, the claimed accuracies are on β, while β∕γ is considered here; however, as the approximation from Eq. (5) is accurate to about 1.5%, this has little practical effect. The deviations from unity appear to reach a constant level for large τ_i and also depend on τ_c. The top panels in Fig. 3 highlight this by showing the results for τ_i < 15 and the restricted range τ_c < 2.5. At large optical depths, the deviation is at most 10%, but at τ_i ~ 1.5 it is 15%.

The minimum in the ratio of exact-to-approximation curve near τ_i ~ 1 can be largely removed by includingthe term Δ, as shown inthe bottom panels of Fig. 3, showing the results for τ_i < 5, τ_c < 2.5, and using Eq. (6) for Δ. The deviations from unity are now quite uniform with τ_i.

Taking the calculations with limits τ_i < 5, τ_c < 2.5 and the approximations Eq. (3), and Eq. (5) with Δ = 0 as reference (the default in Mamon et al. 1988 and Saberi et al. 2019), the reduced χ² is 45 when assigning an “error” of 1.0% per datapoint (see below) in β∕γ. Performing a fit in the range of τ_i < 15, τ_c < 2.5 and the better approximation from Eq. (8) for γ, the best fit is Eq. (5) with α = 1.5222 ± 0.0001 and (9)

with a reduced χ² of 25. This fit is shown in Fig. 4.

Fig. 1 Natural logarithm of γ plotted against continuum optical depth (panel a), with the ratio of the exact solution to the approximation in MJ83, Eq. (3), in panel b. Panel c shows this ratio using the new best fit, Eq. (8).

Fig. 2 Ratio of β over γ plotted against the line optical depth for the approximation in MJ83: Eqs. (3) and (5) with α = 1.5 and Δ = 0 and calculated for τc < 5. The approximation is plotted with small squares in the upper panel. The ratio of the exact solution to the assumed approximation is shown in the bottom panel. For each discrete value of τi (calculated at 0.025 intervals), the vertical lines indicate the range in the continuum optical depth, from 0 (largest value of β/γ) to the maximum value (lowest value). The small vertical lines at the bottom of the upper panel and at the top of the bottom panel indicate the optical depth where the difference between the exact calculation and the approximation is largest.

4 Average escape probabilities

In the previous section, we show that the approximations in MJ83 have limitations. However, it is not clear how this might affect a full calculation including all CO lines which cover a large range in optical depths.

To investigate this further, we considered the optical depths of the 855 CO lines for the standard model in Saberi et al. (2019), with parameters MLR Ṁ = 10⁻⁵ M_⊙ yr⁻¹, wind expansion velocity V_exp = 15 km s⁻¹, and f_CO = 8 × 10⁻⁴, at a distance of r = 10¹⁶ cm, which implies a kinetic temperature of ~70 K (from Eq. (1) in Saberi et al. 2019). The CO line optical depth τ_i ranges from 0.25 to ~ 900 with a median of 35. The continuum optical depth in the line wings of H₂ lines at the wavelengths of the CO lines ranges from 0.12 to >10 000 with a median of 1.3. The dust optical depth τ_d is 4.14.

These optical depths scale as ~ Ṁ f_CO/V_exp/r. Groenewegen (2017) and Saberi et al. (2019) considered MLRs from 10⁻⁴–10⁻⁸ M_⊙ yr⁻¹, expansion velocities from 3.75 to 30 km s⁻¹, and CO abundances from 2 to 16 × 10⁻⁴. The range in distances was taken from where the abundance (cf. Eq. (1)) changes from (0.999 × f_CO) to (). Using the approximations in Groenewegen (2017) for and α as a function of MLR, f_CO, and V_exp, it was determined that the optical depths range from about unity to 10⁻⁴ times those of the standard case considering the entire parameter space.

Table 2 summarises the calculations for the average escape probability (10)

over all CO lines for different scaling factors of the optical depth. The escape probability was calculated as follows: (a) exactly, (b) using the approximations in MGH and Saberi et al. (2019) (Eqs. (3) and (5) with α = 1.5 and Δ = 0), and (c) with the improved approximations (Eq. (5) with α = 1.522 and Eq. (9)) in Cols. 2, 3, and 4, respectively.

For the standard model, the improved approximation is closer to the exact calculation, but still off by about 30%. In the standard model scaled downwards, one can observe the trend that the simple approximation systematically gives average escape probabilitiesthat are too large, while the improved approximation systematically gives lower values than the exact calculation. The simple approximation actually gives results closer to the exact calculation, which is surprising in view of the results presented in Sects. 3.1 and 3.2.

A closer inspection of the optical depths shows that in the standard case in 769 of the 855 lines τ_i > 15 or τ_c > 5 (and both conditions are violated in 365 lines), that is, values for the optical depths outside the range of the fitting formula. Considering the optical depth of all 855 lines for the five sets of models, and choosing escape probabilities of > 0.001 as being the most relevant, we find that in 90% of the cases the optical depths are τ_i < 7 and τ_c < 1. This indicates that the line and continuum optical depths in realistic cases covers a more restricted range than considered in MJ83 and earlier in this paper. Given this finding, the fits to γ and β∕γ were repeated using these limits, and the results are (11)

a value of α = 1.4376 ± 0.0006 in Eq. (5), and (12)

By choosing the error per datapoint as 0.35% and 1.0% earlier, the reduced χ² s were tuned to be unity in both fits, respectively. With this new approximation, the average escape probabilities (the last column in Table 2) are very close to the exact values, except in the standard case, where the functional forms of Eqs. (3) and (5) do not provide an adequate description.

Fig. 3 Same as Fig. 2, but with τc limited to <2.5 (the top two panels). The calculations in the bottom two panels show the results with the Δ term included (Eq. (6)), which make the ratio of exact calculation to approximation more symmetric around unity.

Fig. 4 Same as Fig. 2, but for τc < 2.5 using the new best fit: Eqs. (8) and (5) with α = 1.522 and Eq. (9).

5 Summary and discussion

Improved numerical approximations are outlined to the formalism presented in MJ83 that were at the basis of the CO photodissociation calculations in Mamon et al. (1988) and Saberi et al. (2019). The results in Table 2 show that the average escape probability in realistic cases deviate by 2–3% from the exact value for low to moderate MLRs (and 35% for large MLRs) in the approximation used by Mamon et al. (1988) and Saberi et al. (2019) as well as by ≲ 0.4% (19%) in the improved approximation. For the quantity (), the differences are 0.002–19% and 0.001–11%, with the largest difference for the smallest MLRs, respectively.

The improved approximations were implemented in the code used in Saberi et al. (2019) and some calculations were performed using V_exp = 15 km s⁻¹, f_CO = 8 × 10⁻⁴, and MLRs Ṁ = 10⁻⁴, 10⁻⁵, 10⁻⁶, 10⁻⁷, 10⁻⁸ M_⊙ yr⁻¹. Figure 5 shows the CO abundance versus distance to the star for the five MLRs, from which was determined. The improved formulae lead to a larger photodissociation radius by 7.2% for Ṁ = 10⁻⁸, 2.2% for 10⁻⁷ M_⊙ yr⁻¹, and a ≲0.5% difference for the other MLRs. A larger photodissociation radius is expected as the escape probability is lower in the improved approximation compared to the classical approximation.

As basic assumptions of spherical symmetry, constant velocity, and smooth outflow are identical in the Groenewegen (2017) and Saberi et al. (2019) models (see the introduction) the differences in the CO photodissociation radii (11–60% smaller in Saberi et al. 2019 than Groenewegen 2017) cannot be attributed to them. These differences are due to the different underlying physical implementation, the CO shielding function derived for a static environment that is not representative of an AGB wind, and the escape-probability formalism in an expanding envelope, respectively. It is shown here that the numerical approximations used in the latter formalism have a small effect on the outcome.

The model by Saberi et al. (2019) is currently the most accurate available and covers a large parameter space (in MLR, V_exp, f_CO, ISRF strength). The results in the present paper show that the CO photodissociation radii are underestimated by a few percent for the lowest MLRs in Saberi et al. (2019), but uncertainties in observational estimates of MLR and f_CO lead to larger uncertainties. Resolved observations of CO shells, such as those carried out within the DEATHSTAR⁴ project (Ramstedt et al. 2020), will allow for a detailed comparison of predicted and observed CO photodissociation radii for a large sample in the near future.

Fig. 5 CO abundance profile for the five indicated MLRs. Blue lines indicate the profiles using the approximations used in Saberi et al. (2019), and black lines indicate the profiles using the improved approximations.

Acknowledgements

MS acknowledges the SolarALMA project, which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 682462), and by the Research Council of Norway through its Centres of Excellence scheme, project number 262622.

References

All Tables

All Figures

	Fig. 1 Natural logarithm of γ plotted against continuum optical depth (panel a), with the ratio of the exact solution to the approximation in MJ83, Eq. (3), in panel b. Panel c shows this ratio using the new best fit, Eq. (8).
In the text

Fig. 2 Ratio of β over γ plotted against the line optical depth for the approximation in MJ83: Eqs. (3) and (5) with α = 1.5 and Δ = 0 and calculated for τc < 5. The approximation is plotted with small squares in the upper panel. The ratio of the exact solution to the assumed approximation is shown in the bottom panel. For each discrete value of τi (calculated at 0.025 intervals), the vertical lines indicate the range in the continuum optical depth, from 0 (largest value of β/γ) to the maximum value (lowest value). The small vertical lines at the bottom of the upper panel and at the top of the bottom panel indicate the optical depth where the difference between the exact calculation and the approximation is largest.

In the text

	Fig. 3 Same as Fig. 2, but with τc limited to <2.5 (the top two panels). The calculations in the bottom two panels show the results with the Δ term included (Eq. (6)), which make the ratio of exact calculation to approximation more symmetric around unity.
In the text

	Fig. 4 Same as Fig. 2, but for τc < 2.5 using the new best fit: Eqs. (8) and (5) with α = 1.522 and Eq. (9).
In the text

	Fig. 5 CO abundance profile for the five indicated MLRs. Blue lines indicate the profiles using the approximations used in Saberi et al. (2019), and black lines indicate the profiles using the improved approximations.
In the text