© ESO 2018

1. Introduction

The asymptotic giant branch (AGB) is one of the late stages in the evolution of low and intermediate mass stars (1 M _⊙ ≲ M ≲ 8 M_⊙). At this phase, stars are very luminous (∼10⁴ L_⊙) and cool, with effective temperatures ranging from T_eff ∼ 2000 to 4000 K. The stellar radius of a typical AGB star is of about 1 AU, i.e. R_⋆ ≈ 10¹³ cm. AGB stars are dominated by an intense mass loss process with rates that range, typically, between 10⁻⁷ and 10⁻⁴ M_⊙ yr⁻¹. For complete reviews see, for example, Habing (1996), Herwig (2005), Höfner & Olofsson (2018).

The vigorous mass loss generates a circumstellar envelope (CSE) of gas and dust around these evolved objects. The mass loss is thought to be produced by the combination of two mechanisms: pulsations and radiation pressure on dust (Bowen 1988; Simis et al. 2001). Both mechanisms work together driving a slow, dense stellar wind. The wind is strongly accelerated from the dust condensation region (at a few stellar radii, ≈10¹⁴ cm), reaching terminal velocities of about 5 to 25 km s⁻¹ at distances from the star of ≳10–20 R_⋆. The temperature of the gas decreases along the CSE from 2000–2500 K, close to the stellar atmosphere, to ∼1000 K in the dust condensation region, and progressively falling down to just a few tens of K in the most external envelope layers, at ≈10¹⁶ − 10¹⁷ cm (e.g. Decin 2012).

The material ejected due to the mass loss process, is expected to enrich the interstellar medium (ISM) with stellar nucleosynthesis products. It is believed that about half of the elements heavier than Fe are formed in the AGB phase (Herwig 2005). The analysis of the molecular products found in the AGB CSEs is very useful to derive information on the nucleosynthesis history of the central star and to improve our current understanding of the complex circumstellar chemistry (e.g. Millar 2016).

The relative abundance of carbon and oxygen in the CSE completely determines the dominant chemistry of these envelopes. According to the ratio of carbon over oxygen atoms, we identify envelopes with a carbon-rich (C/O > 1) or oxygen-rich (C/O < 1) chemistry, with a few objects (known as S-stars) showing C/O ratios very close to unity. This chemical branching is very important and depends on the nucleosynthesis processes that have taken place in the interior of these stars during their AGB evolution, which are basically determined by the mass of the progenitor star and its initial metallicity. All stars are born with a C/O ∼ 0.5 and therefore O-rich. However, this chemical composition may change during their AGB evolution. According to theoretical models (Marigo et al. 2003; Karakas et al. 2009), low-mass progenitors (M < 2 M_⊙) remain as O-rich objects throughout the whole AGB evolution, while those stars with higher initial masses, between 2 and 4 M_⊙, may become C-rich objects during the AGB phase due to the thermal pulses and convective processes (dredge-ups) that enrich the surface of the AGB star with carbon, produced as a result of nuclear reactions. Finally, the most massive objects (AGB stars with initial masses in the range 4 M_⊙ < M < 8 M_⊙) remain also O-rich throughout their evolution in the AGB as they do not produce so much carbon, instead they convert carbon into nitrogen via the CN cycle (in a process known as hot bottom burning).

After the AGB phase the convective envelope is expelled and the hot core starts to ionize the surrounding material. When the bulk of the gas in the CSE is ionized the object is recognized as a planetary nebula (PN), see Kwok (2005).

The transition phase between AGBs and PNe is called the post-AGB phase. It is a very fast (≈10³ − 10⁴ yr) phase when the shell, formed in the AGB phase, detaches from the central star and the spherical symmetry is usually broken due to the effect of the onset of fast bipolar or multi-polar winds (Balick & Frank 2002).

The CO molecule is one of the most abundant species found in these objects as it is quickly and easily formed from the O and C atoms present in the envelope and it is also the most stable at the typical conditions found in these envelopes, CO molecules develop strong rotational emission lines which are very useful to trace the physical properties of the molecular envelopes around evolved stars in different evolutionary stages, from the AGB to the PN phase.

The molecular emission of CO has been studied over the last decades mainly at mm and sub-mm wavelength ranges using the low-J transitions that trace outer (> 10¹⁶ cm) relatively cold (10–50 K) regions of the CSE (Bujarrabal 1999; Schöier & Olofsson 2001; Teyssier et al. 2006; Ramstedt et al. 2009). Studying the hotter layers of the envelope, presumably closer to the central star, traced by high-J CO transitions has only been possible with space telescopes like the Infrared Space Observatory (ISO; see, for example Justtanont et al. 2000) and in more recent times with Herschel, particularly, making use of PACS spectroscopy data (see, Groenewegen et al. 2011; Ueta et al. 2014).

Here, we analyse, in a uniform and systematic way, PACS spectroscopy data of a sample of O-rich (and a few S-type) evolved stars from the THROES catalogue (Ramos-Medina et al. 2018). We focus our study on the high-excitation rotational transitions of the CO molecule detected in these stars. In Sect. 2, a general description of our sample is provided. After that, in Sect. 3, we present the main characteristics of the PACS instrument and provide a brief description of the major data reduction steps. In Sect. 4 we show the principal observational results, including our measurements of the CO line fluxes. The rotational diagram (RD) analysis method used in this work and the results derived can be found in Sects. 5 and 6, respectively. The analysis of our results is presented in Sect. 7. Finally, in Sect. 8 we summarize our main conclusions.

2. Sample information

In this work we have studied 26 evolved stars at different evolutionary stages, from the AGB to PN phase, taken from the THROES catalogue (Ramos-Medina et al. 2018). In its current version, the THROES catalogue contains PACS spectroscopic cubes for a total of 114 low-to-intermediate mass evolved stars. In this study, we focus on non C-rich targets, with most of our sources being O-rich objects but also including a few S-type AGB stars (Table D.1). Sources with a dominant C-rich chemistry are studied separately by da Silva Santos et al. (in prep., Paper II). Our sources are grouped in five main categories as originally defined in the THROES catalogue: classical O-rich and S-type AGB stars, OH/IR stars (strong far-IR emitters with intense OH maser lines, typically late-AGB or early post-AGB stars heavily obscured by dust at optical wavelengths), post-AGB objects, and PNe. For this study, we have selected objects with a minimum of three CO transitions detected with a ratio of S/N > 3, a total of 26 out of 66 non C-rich targets in the THROES catalogue.

In Fig. 1 we present IRAS fluxes and colours of the sources of our sample in comparison with the rest of objects (C-rich or CO non-detections of any chemical type) contained in the THROES catalogue. IRAS colours in this diagram are defined as: $[12]-[25]=2.5{log}_{10}\frac{{\mathrm{IRAS}}_{25}}{{\mathrm{IRAS}}_{12}}$ and $[25]-[60]=2.5{log}_{10}\frac{{\mathrm{IRAS}}_{60}}{{\mathrm{IRAS}}_{25}}$. In the colour–colour diagram, the sources included in the THROES catalogue follow a well-known evolutionary sequence, with the AGBs located at the bottom left region (boxes: I, II, IIIa and VII) and the post-AGBs and PNe occupying the upper right area (boxes: V and VIII); more details can be found in van der Veen & Habing (1988). We note that the only OH/IR star located close to other PNe in region V is OH 231.8+4.2, which is a peculiar AGB star with a massive outflow with physical and chemical particularities uncommon for the AGB stage but approaching those of post-AGB objects (e.g. Alcolea et al. 2001; Sánchez Contreras et al. 2015).

All the objects of our sample except two, IRAS 17347-3139 and NGC 6537, have previous estimations of the mass loss rate but using low-J CO lines detected in the sub-millimetre and millimetre ranges (see Table D.2). Only for six objects namely, NML Tau, TX Cam, R Dor, R Cas, W Hya and W Aql, the mass loss rates provided come from recent works (Danilovich et al. 2014; Khouri et al. 2014; Maercker et al. 2016; Van de Sande et al. 2018) which are based on detailed non-local thermal equilibrium (non-LTE) radiative transfer models which make use of a combination of different sets of CO lines taken from several telescopes and including HIFI and PACS observations that, eventually, cover CO rotational lines up to J = 24 − 23. We have paid particular attention to the mass loss rates estimations derived in these works in order to test the results obtained with the rotational diagram method, used in this paper.

Fig. 1. IRAS colour–colour diagram (left) and 100 μm flux vs. [12]−[25] IRAS colour (right) of the THROES targets with good quality IRAS data (Quality Flag = 3) in the 12, 25, 60 and 100 μm bands. The colour–colour diagram is divided in different boxes where sources with common characteristics and evolutionary stages are located (see van der Veen & Habing 1988, for details). Objects with CO line detections, studied in this work, are highlighted using the symbol and colour code shown in the figure. In the rest of the plots presented throughout this paper, AGBs, OH/IRs, post-AGBs and PNe are coloured and symbol-coded as indicated in this figure.

3. Observations and data reduction

PACS (Poglitsch et al. 2010) is an instrument onboard Herschel (Pilbratt et al. 2010) with photometric and spectroscopic capabilities in the FIR wavelength region. The PACS spectrometer covers the wavelength range from 51 to 210 μm in two different channels that operate simultaneously in the blue (51–105 μm) and red (102–220 μm) bands. The Field of View (FoV) covers a 47″×47″ region in the sky, it is structured in an array of 5 × 5 spatial pixels (“spaxels”), each spaxel is a 9.4″×9.4″ square. PACS provides a resolving power between 5500 and 940 (i.e. a spectral resolution of 55–320 km s⁻¹) at short and long wavelengths, respectively, and the point spread function (PSF) of the PACS spectrometer ranges from ∼9″ in the blue band to ∼13″ in the red band. The PSF is described in da Silva Santos (2016) and Bocchio et al. (2016). Complete details about PACS instrument can be found in the PACS Observer’s Manual¹.

In this work we have used one-dimensional (1D) PACS spectra, integrated over the emitting source, taken from the THROES catalogue (Ramos-Medina et al. 2018). THROES is a web-based catalogue that contains fully reduced PACS spectra, in 3D cubes structure (FITS files) and 1D spectra format (FITS and CSV files), for 114 low-to-intermediate mass evolved stars, covering a range of evolutionary stages from the AGB to the PN phase. The THROES catalogue is public and can be accessed online².

The main goal of the THROES project was the production of enhanced PACS spectroscopic data through the systematic careful interactive reprocessing, for full science exploitation of the spectra. Apart from the standard data reduction steps (most of them) common to the pipeline, some extra corrections were added such as: flat-field, background correction and other post-processing corrections applied in particular cases to extract the science-ready 1D PACS spectra from the 3D data cubes. In this process, the extent of the sources and their precise position in the FoV of PACS were taken into account. Particular attention was given to the mispointed cases (those objects not perfectly centred on the central spaxel of the PACS field of view). Full details on the data reduction process of the THROES targets can be found in Ramos-Medina et al. (2018).

The complete 1D spectra reported here are formed by four individual sub-spectra that cover four different wavelength subranges resulting in a final spectral coverage that ranges from 55 to 95 μm (blue band) and from 105 to 190 μm (red band). All the objects included in this work have a full spectral coverage except RR Aql and T Cep, whose 1D PACS spectra cover only two of the spectral subranges: from 55 to 73 μm and from 105 to 147 μm.

4. Observational results

In Figs. 2 and 3 we display the continuum subtracted PACS spectra of the sources studied in this paper from 55 to 95 μm and from 110 to 190 μm, respectively. The continuum was fitted using a non-parametric method after identifying the line-free regions of the spectrum for each target.

Fig. 2. PACS continuum subtracted spectra of our sample in the 55–75 μm range. The 12CO rotational lines are indicated by vertical dotted lines (see Table D.4). The upper-level rotational quantum number Ju of a few selected CO transitions are shown as a reference. The sources are sorted from top to bottom in these five groups: AGBs, S-type stars, extreme OH/IR stars, post-AGBs, and PNe. Within each category, targets are ordered based on the stellar effective temperature from the coolest source to the hottest one.

Fig. 3. The same as in Fig. 2 but covering the spectral range from 110 to 190 μm.

The sources are sorted in Figs. 2 and 3 attending to their classification with the O-rich AGBs at the top, followed by S-type stars, OH/IR stars, post-AGB stars and PNe. Inside each group, the objects are sorted attending to the effective temperature of the central star, with the coolest object at the top and the hottest one at the bottom. The way the objects have been sorted evidences how different the spectra are along the different evolutionary stages from AGB phase to the PN stage. While the spectra of the envelopes around AGB stars show hundreds of molecular emission lines the density of emission lines in post-AGB objects is clearly lower. Finally, the most evolved objects, with highest temperatures of central stars, the PNe, show intense atomic emission lines in their envelopes and a lack of molecular features compared to the AGBs.

Rotational CO transitions in the ground vibrational state (v = 0) are the most prominent molecular lines found in all the spectra. We have detected and measured the integrated flux of CO v = 0 rotational lines with high-J levels, up to J = 42 (E_up ∼ 4700 K). We have checked that, as expected, the resolving power of PACS (of a few 100 km s⁻¹) does not allow us to spectrally resolve the line profiles, which are expected to have full widths at half intensity similar to the terminal expansion velocity of the envelopes, FWHM ∼ 20–40 km s⁻¹ in the majority of our sample. The PACS CO line profiles are also unresolved for a few post-AGB objects and PNs in our sample that are known to have fast outflows sampled by low-J transitions, that is OH 231.8+4.2 (Alcolea et al. 2001).

Apart from CO, we have also detected molecular emission lines from H₂O, including orto- and para- configurations, OH, HCN or SiO amongst others. Some of these lines have been reported for some objects of our subsample in (e.g. Danilovich et al. 2014; Maercker et al. 2016). Atomic lines such as OI (63 μm), OIII (88 μm), NII (122 μm) and CII (158 μm) are also identified in post-AGBs and PNe, but not in AGBs. Finally, broad dust emission features associated to forsterite at 69 μm are also found in 4 objects: NGC 6537, NGC 6302, OH 231.8+4.2 and OH 26.5+0.6. Previous detections of forsterite in evolved stars with Herschel are shown, for example, in de Vries et al. (2015) and Blommaert et al. (2014). The integrated fluxes of the CO rotational emission lines detected in our sample are shown in Table D.4.

Certain CO rotational lines are known to be blended with transitions by other abundant species that lie within the same PACS spectral resolution element. We have identified, for example, the blend of CO (J = 23 − 22) with o–H₂O (414–303), and that of CO(J = 22 − 21) with the ¹³CO(J = 23 − 22) and HNC (J = 28 − 27) transitions. Line blends were excluded in our RD analysis, but these lines are listed and integrated fluxes are reported in Table D.4 for completeness. In this table, we also list the lines that present a poor estimation of the underlying continuum level, although they have not been used for the rotational diagram analysis given their very large flux error bars. Finally, we note that the CO (J = 27 − 26) and (J = 26 − 25) transitions are missing for all targets since there is a gap in the 95–102 μm region of the PACS detector where these lines lie. The errors in the CO line fluxes are derived using the error associated to the fitting, estimated using the rms value of the continuum adjacent to the line, and assuming a flux calibration error of 15%.

In Fig. 4 we compare the integrated flux of the CO (J = 15−14) transition at 173.6 μm (F_{CO 15−14}) with other observational parameters, in particular, with the IRAS 100 μm flux (IRAS₁₀₀), the PACS continuum at 170 μm (PACS₁₇₀), i.e. near the CO (J = 15−14) line, and the [12]–[25] IRAS colour. As expected, there is a correlation between the CO line strength and the IRAS₁₀₀ and PACS₁₇₀ continuum fluxes, which confirms that objects with stronger molecular emission show, typically, stronger continuum emission from dust as well. Of course, in addition to an intrinsic relation between these two observational parameters, there is also a dependence with the distance to the sources, which is partially responsible of the observed relation. The scatter of the points is significantly larger in the F_{CO 15−14} vs. IRAS₁₀₀ distribution, which partially reflects the variability of the majority of the targets (pulsating AGB stars) given that PACS and IRAS observations were obtained at different epochs, as well as different instrument calibration uncertainties, background contamination sources, beam size and pointing, etc. These effects are discussed for the whole sample of the THROES catalogue in Ramos-Medina et al. (2018). We find an anticorrelation between the line-to-continuum (F_{CO 15−14}/PACS₁₇₀) ratio and the IRAS [12]–[25] colour, which are both distant-independent parameters. In general, we find that the ratio between the molecular emission and the dust emission is higher in less evolved objects than in the most evolved ones.

Fig. 4. Top: integrated flux of the CO J = 15−14 transition at 173.6 μm (FCO 15−14) against the IRAS 100 μm flux. Targets with no PACS data beyond ∼145 μm (i.e. not covering this CO transition; RR Aql and T Cep – see Fig. 3) or with bad quality (Quality Flag≠3) IRAS100 data (OH 26.5+0.6) are excluded. Middle: fluxes of the CO J = 15−14 transition and the adjacent continuum (PACS170). In this diagram, OH 26.5+0.6 is included. Bottom: line-to-continuum (F CO 15−14-to-PACS170) flux ratio against the IRAS colour [12]−[25].

5. CO emission analysis: rotational diagrams (RD)

The CO rotational lines provide a fundamental diagnostic of the physical conditions of the molecular gas. As a first approximation, we calculated the rotational population diagrams of the CO (J_u−J₁) v = 0 transitions observed with PACS assuming optically thin emission under LTE conditions. We have used the approach presented by Justtanont et al. (2000), where the total flux observed is used to derive the total mass of the CO-emitting volume adopting the distance to the source indicated in Table D.2.

The concept of the RD analysis is very simple, see Goldsmith & Langer (1999). Assuming LTE, the population of each level N_u can be described by a single rotation temperature T_rot, which is expected to be the same as the kinetic temperature, and given by the Boltzmann distribution:

$$\begin{array}{c}\hfill N_{\mathrm{u}}=\frac{N_{\mathrm{CO}}}{Q(T_{\mathrm{rot}})}{g}_{\mathrm{u}}{\mathrm{e}}^{-E_{\mathrm{u}}/(\kappa T_{\mathrm{rot}})},\end{array}$$(1)

where N_CO is the total number of CO molecules in the emitting volume, Q(T_rot) is the rotational partition function and g_u are the statistical weights for each upper level. Equation (1) can be rewritten in the following form:

$$\begin{array}{c}\hfill ln\left(\frac{N_{\mathrm{u}}}{g_{\mathrm{u}}}\right)=ln\left(\frac{N_{\mathrm{CO}}}{Q(T_{\mathrm{rot}})}\right)-\frac{E_{\mathrm{u}}}{\kappa T_{\mathrm{rot}}}-ln{C}_{\tau },\end{array}$$(2)

where C_τ is the opacity correction term defined as $C_{\tau }=\frac{{\tau }_{\lambda }}{1-{\mathrm{e}}^{-{\tau }_{\lambda }}}$ by, e.g, Goldsmith & Langer (1999). The left hand side of Eq. (2) can be estimated using the relation between the observed flux of a single line and the number of molecules at the upper level (N_u):

$$\begin{array}{c}\hfill N_{\mathrm{u}}=\frac{4\pi F{\lambda }_{\mathrm{ul}}{d}^2}{A_{\mathrm{ul}}hc},\end{array}$$(3)

where F is the integrated line flux, λ is the rest wavelength of the transition, d is the distance to the source, A_ul is the Einstein’s spontaneous emission coefficient associated to the transition, h is the Planck’s constant and c represents the speed of light.

In our RDs we represent $ln\left(\frac{N_{\mathrm{u}}}{{\mathit{g}}_{\mathrm{u}}}\right)$ against $\frac{E_{\mathrm{u}}}{\kappa }$ and we apply a linear regression that allow us to obtain T_rot (derived from the slope) and N_CO (estimated from the intercept in the y-axis) following Eq. (2).

With the values of N_CO and T_rot derived from this first fitting we estimate those of τ and C_τ from Eq. (4) in Sect. 5.1. This initial estimation is used to obtain new values of $ln\left(\frac{N_{\mathrm{u}}}{{\mathit{g}}_{\mathrm{u}}}\right)$ applying the correction shown in Eq. (2). Then we can apply a new linear fit to this new set of values to derive the so-called opacity corrected values for T_rot and N_CO that are the ones adopted for the rest of the analysis. We want to highlight that the opacity correction derived is moderate being the τ values estimated for each object close to 1 (see Sect. 5.1). Very high values of τ would lead to unreliable results from the RD analysis, see Goldsmith & Langer (1999) for details.

5.1. Opacity correction

The optical depth (τ_λ) of a CO line is given by:

$$\begin{array}{c}\hfill {\tau }_{\lambda }=\frac{A_{\mathrm{ul}}{\lambda }_{\mathrm{ul}}^3{N}_{\mathrm{u}}^{\mathrm{col}}}{8\pi v_{\exp }\sqrt{\pi }/(2\sqrt{log2})}\times ({\mathrm{e}}^{(hc/({\lambda }_{\mathrm{ul}}\kappa T_{\mathrm{rot}}))}-1),\end{array}$$(4)

where $N_{\text{u}}^{\text{col}}$ is the column density of CO in the upper level of the transition.

We need an estimation of T_rot and $N_{\text{CO}}^{\text{col}}$ to calculate τ. T_rot is directly extracted from the RDs and $N_{\text{CO}}^{\text{col}}$ is derived from the total number of CO molecules, N_CO, through the following, simplified, equation: $N_{\mathrm{CO}}^{\mathrm{col}}=\frac{N_{\mathrm{CO}}}{\pi {r_{\mathrm{CO}}}^2}$. Where r_CO represents the characteristic radius of the region where the bulk of the CO emission is produced.

The value of r_CO is one of the main sources of uncertainty to compute the optical depth and the subsequent opacity correction term. Unfortunately, as mentioned in Sect. 3, the inner envelope layers where the CO PACS lines are produced are unresolved or barely resolved in the PACS spectroscopy cubes for most of the sources given the large PSF (∼9″–13″) and pixel size (9″.4) of the instrument (see Sect. 3). Therefore, the value of r_CO needs to be adopted based on several criteria. These criteria are described in detail in Appendix B, and a summary is offered in the following paragraphs.

For AGB stars, an estimate of r_CO is obtained by comparing the values of T_rot deduced from our analysis, ∼200–700 K, with the radial temperature distribution T(r) in the envelopes of a number of AGB stars with detailed non-LTE excitation and radiative transfer studies of their CO emission. The gas temperature can be up to 2000 K near the dust condensation radius (∼10¹⁴ cm), and decreases gradually towards the outermost layers following, approximately, a power-law of the type ∼1/r^α, with α ∼ 0.5−1.0, see Van de Sande et al. (2018); Maercker et al. (2016); Danilovich et al. (2014); Khouri et al. (2014). For our sample, we deduce a range of values of r_CO from about 2×10¹⁴ cm to [2–4]×10¹⁵ cm, with the highest mass-loss rate objects requiring the largest values of r_CO. The values of r_CO adopted for our targets are given in Table D.3.

The CO column densities derived from our RD analysis adopting a range of r_CO from ∼10¹⁴ to ∼10¹⁵ cm are moderate, $N_{\text{CO}}^{\text{col}}$ < ∼ 10¹⁹cm2, resulting in line optical depths of the order of ∼0.9, and even lower in some cases, for the CO J = 14 → 13 transition, which is the line with the largest opacity in our study. For AGBs, we have tuned r_CO to a value slightly larger than the radius of the wind layer where τ_{J = 14 → 13} ∼ 1; the latter is indeed the deepest region of the envelope observationally accessible, because of the almost null escape probability ($1/C_{\tau }=\frac{(1-{\text{e}}^{-\tau })}{\tau }\to 0$) from much deeper, optically thicker (τ ≫ 1) layers of the wind.

Not only for the AGB stars but also for the OH/IR stars, we find that values of the radius around r_CO ∼ [2–4]×10¹⁵ cm result in line optical depths close to, but smaller than, unity (typically τ _{J = 14 → 13}∼0.5–0.7) that result in moderate C_τ opacity correction factors. For all the objects studied in this work, the opacity correction is typically negligible for lines higher than J = 19−18.

Unlike for AGBs and OH/IR stars, for post-AGBs/PNe there are no theoretical or observationally constrained temperature profiles of the molecular gas component available in the literature. For these targets, we find that τ_{J = 14 → 13} approaches a value of ∼1 at distances of r_CO ∼ 1×10¹⁶ cm, which therefore have been taken as a lower limit to the radius of the CO emitting volume detected with PACS. The PACS spectroscopy cubes indeed show extended, partially resolved nebulosities consistent with a linear diameter of ≲[1–2]×10¹⁷ cm (da Silva Santos 2016), in some of these cases. Using additional information on the molecular envelope extent from the literature and aiming to treat the sources as homogeneously as possible, we have adopted a common value of r_CO = 6×10¹⁶ cm for the 5 post-AGBs/PNe in our sample, resulting in column densities of N_CO ≲ 1–5×10¹⁷ cm⁻² and small line optical depths of τ_{J = 14 → 13}∼0.02–0.07.

It must be noted that, in addition, we have completed our analysis by systematically exploring a range of radii around the value of r_CO adopted in each case to investigate the impact of this uncertain parameter in our results (see Fig. C.1). As expected, the smaller the r_CO, the larger the opacity correction, which results in lower values of T_rot and larger values of N_CO. We note that after applying the C_τ correction, the slope of the RD increases, as a direct result of the frequency-dependence of C_τ and the typical values of T_rot in our sample, see Appendix B. For example, in AGBs and OH/IR stars, values of r_CO ≳ 5×10¹⁵ cm always result in negligible opacity corrections, whereas for values of r_CO ≲ 8×10¹⁴ cm the opacity becomes larger than one for most targets, except for the sources with the lowest mass-loss rates (Ṁ ∼ 10⁻¹ M_⊙ yr⁻¹, see Table D.2) where the CO lines remain optically thin at these radii.

6. Results

6.1. Rotational temperatures and total emitting mass of H₂

The CO RDs of our targets, from which we derive the rotational temperature (T_rot) and the total number of CO molecules (N_CO) in the emitting volume, are shown in Figs. 5a–c. For most objects, the RDs include CO transitions with upper-level energies that range from E_up ∼ 580 K to E_up ≲ 2500 K and are well fitted using a single temperature component. The single-fit rotational temperatures obtained range from ∼200 K (NGC 6537) to ∼700 K (π Gru), with an average value of T_rot ∼ 450 K.

Fig. 5. RDs of CO for our targets (continues in 5b and 5c). For each source, the values of Trot and MH2 derived from the fits (solid lines) as well as the characteristic radius of the emitting envelope layers adopted (RCO) are shown (Sect. 5). The errors of Trot and MH2 represent the uncertainty of the fit considering the error bars of the individual points, which include formal errors due to the noise and absolute flux uncertainty of the PACS spectra. For a few sources, a double component fit was necessary. For these cases, the values of Trot and MH2 for the so-called warm (blue) and hot (red) components are also indicated.

Fig. 5. continued.

Fig. 5. continued.

For a number of strong CO emitters (NML Tau, χ Cyg, W Aql, NML Cyg and IRAS 17347-3139) transitions with upper-level energies of up to E_up ∼ 5000 K are also detected. In these cases, the data points in the RD over the full E_up covered are not well reproduced by a single straight-line fit. A double-slope in the RD is expected if a range of temperatures occur along the line of sight.

We have formally tested that a double-slope fit is indeed most appropriate in these cases using statistical tools for model selection: F statistics and Bayesian Information Criterion (BIC), implemented in R programming language in the strucchange ³ package (see Paper II for more details). We refer to the CO emitting regions with the lowest and highest T_rot values derived from the double-fit as the “warm” and “hot” CO components, respectively (Figs. 5a–c). The rotational temperatures of the warm and hot components are in the range ∼370–480 K and ∼570–900 K, respectively. As expected, the values of T_rot obtained from a single-fit to the RDs are intermediate to those of the warm and hot components.

The total number of CO molecules in the emitting volume, N_CO, derived from the RDs has been converted to total mass of molecular hydrogen using:

$$\begin{array}{c}\hfill M_{{\mathrm{H}}_2}=\frac{N_{\mathrm{CO}}}{{\chi }_{\mathrm{CO}}}{m}_{{\mathrm{H}}_2},\end{array}$$(5)

where m_H2 is the mass of a molecule of hydrogen (m_H2 = 2.0159 a.u. = 3.3475 × 10⁻²⁴ g) and χ_CO is the CO-to-H₂ fractional abundance, assumed to be 2 × 10⁻⁴ for O-rich objects (Ramstedt & Olofsson 2014) and 6 × 10⁻⁴ for S-type stars (Ramstedt & Olofsson 2014; Danilovich et al. 2014).

The resulting values of M_H2 and T_rot are listed for each target in Table D.3, together with a few other relevant parameters (such as the adopted r_CO for the moderate opacity correction applied, Sect. 5.1). In Fig. 6 we show the M_H2 and T_rot values for each target, with and without the opacity correction, and the final histogram distributions.

Fig. 6. Values of Trot and MH2 from our CO rotational diagram analysis (single fit; Sect. 6) before and after applying the opacity correction (grey and coloured symbols, respectively).

The single-fit values of the total mass of the CO-emitting volume range between M_H2 ∼ 3 × 10⁻⁶ M_⊙ (X Her) and ∼5 × 10⁻² M_⊙ (IRAS 17347-3139), with an average value of M_H2 ∼ 4 × 10⁻³ M_⊙. The opacity correction applied in all cases is small-to-moderate, which results in a mass-correction factor less than 20% in most sources. The source with the largest opacity correction (of a factor ∼1.5) is NML Cyg, the AGB star with the largest M_H2, and highest mass-loss rate, in our sample (Table D.2).

The mass of the warm CO component is always larger than the mass of the hot component. For the warm component, these values vary from M_H2 ∼ 1.3 × 10⁻⁵ to 7.8 × 10⁻² M_⊙, and for the hot between 3.7 × 10⁻⁶ and 1.6 × 10⁻² M_⊙. The warm-to-hot M_H2 and T_rot ratio, and their relation, is further examined in Sect. 7.

We have investigated if there are any trends between T_rot and M_H2, derived for our RD analysis, and other observational parameters, such as the terminal expansion velocity of the envelope (v_exp) and the strength of the CO (J = 15−14) transition (Fig. 7). First, we notice that objects in the most advances stages of the evolution, namely, post-AGBs and PNe, systematically have the lowest temperatures, T_rot ∼ 200–300 K, and the highest masses, M_H2 > 0.01 M_⊙ (see also Fig. 6). IRAS 17347-3139 is the only object with an exceptionally large T_rot for its class. The distribution of T_rot and M_H2 for post-AGB and PNe is also much narrower than for AGBs and OH/IR. The latter show a broad range of T_rot ∼ 200–700 K and M_H2 ∼ 10⁻⁶ − 10⁻³ M_⊙. This segregation of the sources exists before and after the opacity correction has been applied; the opacity correction is indeed marginal for post-AGBs and PNe (grey symbols in Fig. 6). The observed bias in favour of massive post-AGB/PNe could reflect that, in our sample, post-AGB/PNe are on average more distant (d > 1 kpc) than AGBs and OH/IRs stars, which hampers CO line detections in the less massive post-AGB/PNe.

Fig. 7. Relation between MH2 and Trotderived from our CO RD analysis and other magnitudes. Top-left: MH2 vs. Trot. For objects (a total of five) with a double-Trot component, both the single-fit and double-fit values are connected with dashed lines for an easier identification. Top-right: values of MH2 (single-fit) vs. vexp taken from the literature. Bottom panels: integrated fluxes of CO (J = 15−14) line vs. MH2 (left) and the same line normalized to its adjacent continuum vs. Trot (right).

We observe an anticorrelation between M_H2 and T_rot (Fig. 7, top-left panel). This is most clearly noticed amongst AGB and OH/IRs, which are more numerous and cover a broader range of T_rot and M_H2 values, but it is also followed by post-AGB/PNe. We note that for objects for which a double-T_rot component has been identified (NML Tau, χ Cyg, W Aql, NML Cyg and IRAS 17347-3139), both the double- and the single-fit values (connected with dashed lines in Fig. 7, top-left panel) follow the trend. The object IRAS 17347-3139 is an outlier, since it has a much larger T_rot than anticipated based on the overall properties of the rest of our sample. This object is known to be at the earliest stages of the post-AGB evolution and making a fast transition to the PN phase at present (Table D.2, see Gómez et al. 2005; Tafoya et al. 2009).

For AGBs and OH/IRs stars, there is also an evident linear correlation between M_H2 and the terminal expansion velocity of the envelope (from the literature – Table D.3). As we will see in Sect. 7.4, this is linked to an analogue Ṁ-to-v_exp relation, which is known to exist based on previous low-J CO studies of AGB CSEs (e.g. Danilovich et al. 2015) and that we corroborate in this work using higher-J CO lines. The post-AGB objects and PNe in our sample, however, show very similar values of the mass over a the full range of velocities.

As expected, the integrated flux of the CO lines correlate with the total mass of the CO-emitting volume and, as we will see later, with the average mass-loss rate of the envelopes. As an example we show this relation for the CO J = 15−14 transition (E_up ∼ 660 K) in Fig. 7 (bottom-left). There is a significant scatter in the line flux vs. M_H2 relation since the latter also depends on the excitation temperature (T_rot), the velocity field in the envelope, the optical depth, etc. (Sect. 5).

Finally, we investigated if a correlation exist between the CO J = 15−14 line integrated flux normalized to the underlying ∼170 μm-continuum observed with PACS (F _{CO 15−14}/PACS₁₇₀) and the parameters derived from the RD analysis. The strongest correlation is found with T_rot (Fig. 7, bottom-right), which is somewhat expected since the mass and the distance dependence is removed in the line-to-continuum flux ratio. The relation between F _{CO 15−14}/PACS₁₇₀ with M_H2 (not shown) has a much larger scatter and suggests a certain anticorrelation (given that T_rot anticorrelates with M_H2 in our sample).

In Sect. 7 we further discuss the observed trends.

6.2. Mass-loss rate

Using the values of the total emitting mass, we calculate the mass-loss rates (Ṁ) for each source in a simplified manner dividing the total mass by the crossing time of the CO-emitting layers, that is:

$$\begin{array}{c}\hfill \dot{M}=\frac{M_{{\mathrm{H}}_2}{v}_{\exp }}{r_{\mathrm{CO}}},\end{array}$$(6)

Apart from M_H2, derived from the RDs, we need to adopt a value for the expansion velocity (v_exp), which has been taken from literature, as well as for the characteristic radius of the CO emitting region, for which we use the same value used for the estimation of the opacity correction, r_CO. The values of Ṁ are given in Table D.2.

In Sect. 5, we already mentioned the uncertainties associated to the adoption of a particular linear radius for the emitting region of each source and the importance of exploring a wide range of values to determine a reasonable estimation of the opacity correction and, thus, of M_H2. The same kind of analysis is even more important now, as the impact of r_CO in the mass-loss rate estimation is larger than in the calculation of the opacity correction (it affects M_H2 but also directly Ṁ as ∼ 1/r_CO). Another uncertainty arises in v_exp, however, this parameter is rather well known (from previous works), to an accuracy of ∼10–20% typically.

In Fig. C.2, we show, for every object, the different values of Ṁ and T_rot that correspond to three different values of the radius around the adopted value. For each source, we have also plotted a grey region that represents the broad interval of mass-loss rates found in the literature. The latter are usually derived from the analysis of low-J CO lines and, therefore, they represent the average mass-loss rate in the outer envelope layers. These values of Ṁ from the literature have been re-scaled after matching the expansion velocity (v_exp), the distance (d) and the relative abundance of CO ( χ _CO) to the values adopted here, which is needed for a proper comparison with our estimates of Ṁ. The values of the mass-loss rates reported in the bibliography for each object are displayed in Table D.2 together with the associated interval range of mass-loss rates after scaling and relevant references. To our knowledge, for IRAS 17347-3139 and NGC 6537 there are no previous estimates of the gas mass-loss rate in the literature.

Our estimates of the mass-loss rates are shown in Fig. 8 plotted against the rotational temperature, the expansion velocity, and the CO (J = 15−14) line strength to explore any trends that may be present. In the Ṁ vs. T_rot diagram we also include the mass-loss rates associated to the warm and hot component for objects with a double-slope in their RDs. There is an anti-correlation between Ṁ and T_rot, which was expected based on a similar trend observed directly between M_H2 and T_rot (Fig. 7). This trend together with the clear correlation between Ṁ and v_exp and, more scattered, trend with F_{CO 15−14} are briefly discussed in Sect. 7.

Fig. 8. Top: relation between Ṁ and Trot. Single- and double-fit (warm and hot) values, whenever these exist, are connected using dashed lines. Middle: relation between Ṁ and vexp. The distribution approaches a power-law of the type Ṁ ∝ vexp2.5 (dashed line). Bottom: integrated flux of the CO J = 15−14 line vs. Ṁ in a logarithmic scale.

7. Analysis and interpretation of the results

In this section, we analyse and discuss the results obtained from the RDs as well as the correlations found between different magnitudes.

7.1. Rotational temperature and total mass

For most of our sources, a single-T_rot component was sufficient to obtain a good fit of the CO RD, leading to single-fit rotational temperatures and total (H₂) masses of the CO-emitting region of T_rot ∼ 200–700 K and M_H2 ∼3 × 10⁻⁶ − 5 × 10⁻² M_⊙, respectively (Sect. 6).

The rotational temperatures of a few 100 K derived in this work using CO-PACS data (J_u ∼ 14–27, and up to J_u = 40 lines in some cases) are clearly larger than the typical values of T_rot ∼ 10–100 K found from CO low-J line studies (most commonly using J_u < 4 lines at mm and sub-mm wavelengths) in the external layers of the CSE around evolved stars, at ≈10¹⁶ − 10¹⁷ cm from the center (e.g. Teyssier et al. 2006, and references therein). This corroborates that, as anticipated, the CO-PACS lines sample the warm, presumably inner, regions of the CSEs, probably at a typical/average distance from the central star of about ≈10¹⁵ cm for AGB stars (e.g. Danilovich et al. 2014).

We find that the total mass of the warm inner layers where the PACS CO (J_u > 14) lines mainly arise is typically a factor ∼100–1000 lower than the masses deduced from low-J CO line studies in the literature, reinforcing the idea that CO PACS lines sample mainly the inner part of the envelopes, which contains a small fraction of the total envelope mass.

For a total of five targets in our sample (NML Tau, π Gru, W Aql, χ Cyg and IRAS 17347-3139), the RD diagram shows a double slope and, thus, a double-T_rot fit is most appropriate. Deviations from a single straight line fit of the RDs, most likely denoting a temperature stratification along the line of sight (see below), are not uncommon and have been previously found in a number of AGB and post-AGB CSEs using ISO and/or SPIRE CO spectra (e.g. Justtanont et al. 2000; Wesson et al. 2010; Matsuura et al. 2014; Cernicharo et al. 2015; Cordiner et al. 2016).

The rotational temperatures and masses of the hot and warm components of these five objects are displayed in Fig. 7 (together with the single-fit T_rot and M_H2 values derived for the whole sample). The mass of the warm component is always larger (by a factor 2–10) than that of the hot component. Adopting a constant mass-loss rate and a gas temperature that decreases with the distance, the difference between the mass of the warm and hot component suggests a CO-emitting volume smaller for the hot component than for the warm one.

We compare in detail the hot-to-warm M_H2 and T_rot ratios in Fig. 9. We find that the mass and temperature ratios are well correlated and follow approximately a power-law function. For AGB CSE with an inverse-square density profile, the latter is indeed expected if the temperature stratification within the emitting layers is described by a power-law function, T(r)∼1/r^α. This is because, in this case, the total mass within a given CO-emitting volume is roughly proportional to its radius (M_H2 ∝ r) and, therefore, it is easy to derive that:

$$\begin{array}{c}\hfill \frac{T_{\mathrm{rot}}^{\mathrm{hot}}}{T_{\mathrm{rot}}^{\mathrm{warm}}}\sim {\left(\frac{M_{{\mathrm{H}}_2}^{\mathrm{hot}}}{M_{{\mathrm{H}}_2}^{\mathrm{warm}}}\right)}^{-\frac{1}{\alpha }},\end{array}$$(7)

Fig. 9. Ratio of the mass and rotational temperatures for the hot and warm components, $M_{{\text{H}}_{\text{2}}}^{\text{hot}}/M_{{\text{H}}_{\text{2}}}^{\text{warm}}\hspace{0.17em}\text{vs}.\hspace{0.17em}{T}_{\text{rot}}^{\text{hot}}/T_{\text{rot}}^{\text{warm}}$, for sources with a double-Trot component, namely, the AGB stars NML Tau, π Gru, W Aql, and χ Cyg and the yPN IRAS 17347-3139. The solid line is a fit by a power-law function (y ∼ x−1/α) to the AGB stars’ data – see text in Sect. 7.1. A similar trend is found in our analogue study of the THROES C-rich sample reported in Paper II.

Fitting a power-law function to the set of points corresponding to our AGB stars⁴, we obtain α = 0.5 ± 0.1. This value is at the low end of the range of exponents found in the literature (from a combination of theoretical modelling and observational works) for AGB CSEs, typically α ∼ 0.5 − 1 (e.g. De Beck et al. 2010; Maercker et al. 2016; Danilovich et al. 2015, and references therein). There is some indication from recent works that the kinetic temperature distribution is shallower, with values of α down to ∼0.4, for the inner CSE layers (∼5×10¹⁴–3×10¹⁵ cm; De Beck et al. 2012; Lombaert et al. 2016) than for the outer regions, where the steepest temperature variations (α∼1–1.2) are found (≳10¹⁶ cm; Teyssier et al. 2006). This is indeed in very good agreement with the small value of α found in the inner envelope layers from our analysis. (The hot-to-warm mass vs. T_rot trend with a similar value of α is also deduced in our sample of C-rich targets studied in Paper II). We believe that the empirical relation between the hot-to-warm ratio of M_H2 and T_rot therefore supports that the double-T_rot component identified in some of our objects (and possibly also present in more targets for which CO J_u > 27 transitions remain undetected), is a natural consequence of the temperature stratification across the inner envelope layers. Detailed non-LTE excitation and radiative transfer of the CO emission including high-J transitions, J_u ∼ 14–40, are needed for a proper characteriztion of the temperature structure in the warm inner-wind regions of AGB stars. See Appendix A for a discussion of non-LTE effects on the RD analysis.

7.2. Comparison with other works: Ṁ_THROES vs. Ṁ_Bibliography

In Fig. 10 we show, for each source, the mass-loss rate estimated in this work for the adopted radius (r_CO) against a mean value of gas mass-loss rates found in the bibliography (Ṁ_Bibliography). We assign error bars of a factor 3 to our estimates of the mass-loss rate; this is a common uncertainty value adopted in most molecular line-emission studies (e.g. Maercker et al. 2016). The error bars of Ṁ_Bibliography represent the low and high ends of the broad range of values found in bibliography after scaling to our d, v_exp and χ_CO values, and is always ∼3 or more. For objects with only one value of Ṁ from the literature, an uncertainty factor of 3 has also been assumed after scaling.

Fig. 10. Comparison of the mass-loss rates from this work and from the literature (Tables D.2 and D.3, respectively). The solid line represents a 1:1 ratio. Target colour and symbol code as in Fig. 1. The big circles are used to locate a few targets with previous ṀBibliography estimates obtained from non-LTE excitation and radiative transfer analysis of CO lines, including some of the PACS transitions.

The comparison between our mass-loss rates and those found in the bibliography is in quite good agreement for all the objects, see also Fig. C.2. As already mentioned, we are particularly interested in the comparison of our Ṁ with those obtained from comprehensive non-LTE excitation and radiative transfer models that make use of some CO rotational lines (normally, J < 20) measured with PACS, which exist for a small sample of objects: for R Cas, NML Tau, TX Cam, R Dor, W Hya and W Aql (Van de Sande et al. 2018; Maercker et al. 2016; Khouri et al. 2014; Danilovich et al. 2014). We find that, also in these cases, discrepancies with our values are always within 50% for all sources except W Hya, demonstrating that the simple and approximate RD approach and the assumptions (e.g. the adopted size of the emitting layers) used in this work yield reasonable first estimate results.

Considering again our sample as a whole, our mass-loss rate estimations of post-AGBs and PNe seem to be systematically lower than those in the literature by a ∼20%. In principle, this small difference could suggest that the radius of the CO-emitting region is slightly smaller than the value of r_CO adopted by us (at least in some objects). As it is shown in Fig. C.2, a fine tuning of r_CO can be performed to obtain a perfect match between Ṁ_Bibliography and our estimates for each individual object. However, we prefer to treat the sources as homogeneously as possible and to be conservative avoiding the use of very small values of r_CO, which would result in large opacities ( τ ≫ 1) and would consequently lead to unreliable opacity corrections.

A similar discrepancy is reported by Danilovich et al. (2015) from a CO-line study (including non-LTE radiative transfer modelling) of AGB CSEs in the SUCCESS programme including CO transitions from J = 1 − 0 to J = 9 − 8, as well as J = 14 − 13 in some targets. Their models give mass-loss rates that are on average 40% lower than those derived from past studies (see their Fig. 7). Danilovich et al. (2015) note that, as already mentioned here, most past studies are based on observations of a small number of low-J CO transitions (with similar excitation conditions, i.e., E_up), and conclude that taking a large sample of CO transitions, including high-J lines, properly into account when modelling mass-loss rates, results in lower derived values.

It is important to stress that the mass-loss rates derived in this work represent average values during the period of time when the warm, inner envelope layers were ejected, that is, during the last < 50–100 years, for AGBs, and the last ∼600–2000 years, for PNe (considering the values of r_CO and v_exp given in Table D.3). Since the low-J lines commonly used in the previous studies to derive Ṁ_Bibliography mainly trace the outer envelope layers, the discrepancies observed with respect to past studies could reflect a real time variation of the mass-loss rate in some cases. This is indeed the case of OH 231.8+4.2, where mass-loss rate when the bulk of the nebula was ejected during the last ≈1000 yr (∼2×10⁻⁴ M_⊙, Alcolea et al. 2001) is much larger than the present-day mass-loss rate (< 2×10⁻⁵ M_⊙), Sánchez Contreras et al. 2002). The upper limit to the present-day mass-loss rate of OH 231.8+4.2 is indeed in very good agreement with our estimate using high-J CO transitions.

Finally, as a word of caution we note that the mass-loss rates estimated for post-AGB objects and PNs are particularly uncertain as it is not easy to define essential parameters like r_CO or v_exp for these spatio-kinematically and morphologically complex objects. This makes the comparison with values from past studies subject to larger uncertainties. In spite of this, the agreement is reasonably good for the two post-AGBs and one PN in our sample with values of Ṁ reported in the literature.

7.3. Mass and mass-loss rate vs. T_rot

As shown in Fig. 7, we find a statistically significant anticorrelation between the total mass of the PACS-CO emitting region, M_H2, and the rotational temperature, T_rot, derived from our RD analysis. As expected, a similar trend is present between the mass-loss rate, Ṁ, and T_rot (Fig. 8).

Prior to this work, Teyssier et al. (2006) and Danilovich et al. (2015) noticed a tentative Ṁ vs. T_rot anticorrelation based on their studies of the intermediate-to-outer regions of the envelopes of AGB and post-AGB stars using submm/mm-wavelength CO transitions. Particularly, Danilovich et al. (2015) reported Ṁ and T_kin at a representative distance of ≈10¹⁶ cm (i.e. ∼100 times the dust condensation radius) for a sample of 53 AGB stars. From their analysis, using self-consistent non-LTE radiative transfer models of the dust continuum emission and the CO emission transitions over a range of excitations (from J = 1 − 0 to J = 9 − 8, and in a few cases, up to J = 14 − 13), Danilovich et al. (2015) found that in general high mass-loss rate objects show relatively low temperatures (see their Fig. 6).

The temperatures of the regions studied by Danilovich et al. (2015), at distances ∼[1–2]×10¹⁶ cm, range between ∼10 and 80 K to be compared with the larger temperatures, ∼200–600 K, found by us at distances typically a factor ∼10 smaller in our AGB+OH/IR targets. This indicates that this trend is preserved from the inner envelope layers, traced by J_u ≥ 14 CO lines, through the outer parts of the CSE, sampled by lower-J transitions.

As discussed by Danilovich et al. (2015), a possible reason for the anticorrelation observed is the more efficient line cooling expected for high mass-loss rate objects (see also Groenewegen 1994, and references therein). In addition to this, we believe that optical depth effects can partially contribute to, and reinforce, the observed anticorrelation. This is because for high mass-loss rate objects, the τ _CO∼1 surface layer, which is the deepest region that can be observationally sampled by a given CO line, is located at larger distances from the star where the gas kinetic temperature is lower.

7.4. Mass and mass-loss vs. expansion velocities

In Sect. 6 we already mentioned the presence of evident correlations between M_H2 and v_exp and between Ṁ and v_exp (Figs. 7 and 8, respectively).

The correlation found between M_H2 and v_exp is stronger if we consider only the AGB and OH/IR stars in our sample. In AGB and OH/IRs stars, the observed correlation could be telling us about the processes which are responsible for the acceleration of the stellar wind. Since AGB winds are primarily driven by radiation pressure on dust, the most luminous objects (which are also the most massive) will exhibit the largest expansion velocities.

The fact that post-AGB and PNs deviate from the observed trend is probably related with the different mass-loss mechanism (not driven by radiation pressure on dust) in these latest stages (Bujarrabal et al. 2001).

Moreover, the expansion velocities of post-AGBs and PNs do not represent the terminal velocity of the present-day stellar wind but the result from the complex hydrodynamical interaction between the stellar high-velocity winds, that appear beyond the AGB phase, and the slow ambient circumstellar material (expelled in the previous AGB stage).

Based on the above conclusion that, in AGB stars, v_exp correlates with M_H2 via the stellar luminosity, the trend found between v_exp and Ṁ is not surprising (Fig. 8). In this case, however, the post-AGBs and PNs in our sample do not deviate significantly from the general trend, although, as we mentioned above, the values and interpretation of v_exp as well as the Ṁ derived (under simplifying assumptions) have to be regarded with caution.

The Ṁ-to-v_exp relation has been observed before in other CO line studies in AGB CSEs. For example, Danilovich et al. (2015) report a correlation between Ṁ and v_exp, including stars of different chemical types (O-rich, C-rich and S-stars), which is in very good agreement with our data-points (in spite of the lower-J transitions, J < 14, used by these authors). We have fitted a power-law function to the observed distribution (dashed line in Fig. 8-middle panel), the best fit gives an index of about 3 (i.e., Ṁ ∝v_exp^2.5 ) although the uncertainties associated to the fitting are important (2.5 ± 1).

7.5. Other trends

As shown in Sect. 6, we find a clear and expected correlation between the strength of the CO line emission, in particular, of the J = 15−14 line (F _{CO 15−14}), with M_H2 and Ṁ (Figs. 7 and 8, respectively). The correlation found by us is in very good agreement with that reported recently by Lombaert et al. (2016) using a number of CO transitions observed with PACS to model the envelopes of a sample of C-rich AGB CSEs with H₂O emission. The dependence with the distance to the sources of F _{CO 15−14}, M_H2 and Ṁ partially explain the observed trend; additionally, the line flux is proportional to the amount of mass in the emitting volume (Eq. (3)), although other factors, such as the opacity and the excitation temperature, have an important influence.

The influence of the excitation temperature in the resulting strength of the CO lines relative to the continuum becomes evident in Fig. 8-bottom panel, where we see a clear correlation between F_{CO 15−14}/PACS₁₇₀ and T_rot, which are both distance independent parameters. We find that objects with the lowest values of T_rot show small line-to-continuum (F_{CO 15−14}/PACS₁₇₀) ratios. We have found from this work that post-AGBs and PNs have the most massive envelopes in our sample (> 0.01–0.1 M_⊙), but in spite of this, they are relatively weak high-J CO emitters as a consequence of the relatively low excitation conditions prevailing in their envelopes (T_rot ∼ 200–300 K). Other factors, such as a different gas-to-dust mass fraction at different evolutionary stages, may also lie behind the observed trend.

8. Summary

We have analysed Herschel-PACS 1D spectra of 26 evolved low-to-intermediate mass stars, in different evolutionary stages from the AGB to the PN phase, contained in the THROES catalogue⁵. Our sample includes 23 O-rich and 3 S-type objects (C-rich targets are separately studied in Paper II).

RDs using high-excitation CO v = 0 rotational emission lines, with upper-level energies E_up ∼ 580–5000 K, have been created to study the warm inner layers of the CSEs of these objects. We derive fundamental physical parameters such as: rotational temperatures (T_rot), total molecular mass in the CO-emitting layers (M_H2) and average mass-loss rates during the ejection of these layers (Ṁ).

The rotational temperatures are found to vary from ∼200 K to ∼700 K. For AGB and OH/IR stars, which represent ∼80% of our sample, the rotational temperatures are around T_rot ∼ 500 K and systematically larger than for the most evolved objects, post-AGBs and PNs (T_rot ∼ 200 K). These temperatures are one order of magnitude higher than estimated from past low-J CO line studies of the outer envelope. This is expected as, in the PACS far-infrared range, we are tracing gas that is closer to the central star (≈10¹⁵ cm for AGBs and ≈10¹⁶ cm for post-AGBs and PNs).

The total mass of the inner CSE regions where the PACS CO lines arise ranges from M_H2 ∼ 10⁻⁶ to 10⁻²M_⊙ in our sample, with the most evolved objects being the most massive. The values of M_H2 found in this work are typically 100–1000 times lower than the total envelope mass obtained from previous low-J CO line studies.

The mass-loss rates estimated are in the range Ṁ ≈ 10⁻⁷ − 10⁻⁴ M_⊙ yr⁻¹, in agreement (within uncertainties) with values found in the literature.

For some strong CO emitters in our sample (NML Tau, χ Cyg, W Aql, NML Cyg and IRAS 17347−3139) a double temperature (hot and warm) component is inferred. The temperatures of the warm and hot components are ∼400–500 K and ∼600–900 K, respectively. The mass of the warm component (∼10⁻⁵ − 8 × 10⁻² M_⊙) is always larger than that of the hot component, by a factor ∼2–10.

The warm-to-hot M_H2 and T_rot ratios in our sample are correlated and are consistent with an average temperature radial profile of ∝ r^{−0.5 ± 0.1}, that is, slightly shallower than in the outer envelope layers, in agreement with recent studies.

We have explored the distributions of the different magnitudes and physical parameters of these objects and if any trends exist. We find that M_H2 and Ṁ are both anticorrelated with T_rot. This is interpreted as a combination of two different processes: CO line cooling and opacity effects. Previous works have also reported similar trends based on low-J CO line observations.

A strong correlation has been found, also, between M_H2 and v_exp, particularly for AGB stars. This correlation is also present in the Ṁ vs. v_exp distribution, as expected. This trend, which is observed in past studies using low-J CO transitions, is consistent with the wind acceleration mechanism being more efficient in the more luminous and massive stars.

Acknowledgments

We acknowledge support from the Faculty of the European Space Astronomy Centre (ESAC). C.S.C. acknowledges financial support by the Spanish MINECO through grants AYA2016-75066-C2-1-P and by the European Research Council through ERC grant 610256: NANOCOSMOS. We also acknowledge our referee, Kay Justtanont, for her useful and precise comments which have helped to improve our work.

References

Appendix A: Non-LTE effects

The rotational diagram (RD) method is a classical and widely used line-emission analysis technique useful to obtain a first estimate of the excitation temperature and mass of the gas in the line emitting region (Goldsmith & Langer 1999; Justtanont et al. 2000). As explained in Sect. 5, it is based on two major assumptions: optically thin emission and excitation under LTE conditions. A canonical opacity correction factor, C_τ, has been included to take into account moderate optical depth effects (Sect. 5). Here, we examine whether non-LTE excitation effects are expected to be significant in our sample and, if this is the case, what is the impact on the values of M_H2 and T_rot inferred from our RD analysis. We present several arguments that suggest that although small LTE deviations cannot be ruled out (particularly, in the lowest mass-loss rate objects), we expect a negligible impact on the derived values of M_H2. Additionally, the double slope in the RD observed for some targets predominantly reflects the temperature stratification in the emitting layers.

A first simple evaluation of LTE deviations is obtained from the comparison of the critical densities of the CO transtions studied here and the densities expected in the CO-emitting regions in our sample. For most (80%) of our targets, the far-IR CO transitions detected and used in our RD analysis arise from states ∼500–2000 K above the ground, i.e. with upper-level rotational quantum numbers J_u ∼ 14–27. According to most recent collisional rate coefficients computed by Yang et al. (2010)⁶, these transitions have critical densities of n_crit ∼ 5×10⁵ and 3×10⁶ cm⁻³, for J_u = 14 and J_u = 27, respectively, for a range of temperatures consistent with our estimates, ∼200 and 700 K, see Fig. 10 in Yang et al. (2010).

From the lowest to the highest mass-loss rate objects in our sample (Ṁ ∼ 2×10⁻⁷–1×10⁻⁴ M_⊙ yr⁻¹, Fig. 8), densities of n_H2 ≳10⁶ − 10⁸ cm⁻³ are expected at the typical radius adopted (from r_CO ∼ 4×10¹⁴ to r_CO ∼ [1–4]×10¹⁵ cm, respectively – see Table D.2, Appendix B) consistent with CO population levels predominantly thermalized in most cases.

As already pointed out in pioneering studies of ISO-LWS observations of high-J CO lines using the RD approach, see Justtanont et al. 2000, the mere fact that we detect CO J_u ∼ 14–27 transitions together with a roughly linear distribution in the RD is a hint of high-density gas and level populations in this range of J_u close to thermal equilibrium.

For a few targets, we also detect CO lines with upper-levels above J_u > 27, which have critical densities of n_crit ∼ 10⁶ − 10⁷ cm⁻³ and, therefore, some LTE deviations cannot be rule out (but we note that these transitions presumably arise from regions closer to the center, within the dust formation layers).

The simple calculations above are indeed corroborated by non-LTE excitation and radiative transfer analysis computations of a selection of high-J CO transitions observed with PACS (from J_u = 14 to 38) on a sample of C-rich AGB CSEs by Lombaert et al. (2016). These authors conclude that over the broad range of mass-loss rates studied in their work, Ṁ ∼ 10⁻⁷–2×10⁻⁵ M_⊙ yr⁻¹, which is comparable to the range in our sample, the CO molecule is predominantly excited through collisions with H₂, with a minor effect of far-IR radiative pumping due to the dust radiation field.

We stress that even under non-LTE conditions, for a simple diatomic molecule like CO, the RD method provides a reliable measure of the total mass within the emitting volume. This is because, although T_rot may deviate from the kinetic temperature, T_kin, in regions where the local density is lower than the critical densities of the transitions considered, it does describe quite precisely the molecular excitation (i.e. the real level population). To test this idea and to check the reliability of our analysis, we have compared our values of M_H2 and Ṁ with those obtained from state-of-the-art line excitation and radiative transfer models of high-J CO transitions, sometimes including PACS data, performed by other authors in objects of our sample. To our knowledge, these have been published for a total of 6 sources: R Cas, R Dor, TX Cam, NML Tau, W Hya and W Aql (see Van de Sande et al. 2018; Maercker et al. 2016; Danilovich et al. 2014; Khouri et al. 2014 and references therein). For these objects, the CO emission over a broad range of E_up has been modelled to derive the temperature, density and velocity profile across their ∼10¹⁴ − 10¹⁷ cm CSE layers, as well as the mass-loss rate of the stellar wind. As shown in Sect. 6, discrepancies between our values of Ṁ and those obtained from comprehensive non-LTE radiative transfer analysis are very small, always within 50%, demonstrating that the simple and approximate RD approach and the assumptions (e.g. the adopted size of the emitting layers) used in this work yield reasonable results.

An additional test of the reliability and uncertainties of the RD method is presented in Paper II, where we perform our RD analysis using multi-epoch CO-PACS spectra for the C-rich AGB star IRC+10216. This object shows prominent line variability (mainly of high-excitation transitions of CO and other molecules) denoting radiative farIR (ro-vibrational) pumping via the dust radiation field modulated by the stellar pulsation (Cernicharo et al. 2014; Teyssier et al. 2015; He et al. 2017). This object, thus, serves as a benchmark for testing the effect of radiative vs. collisional excitation and its impact of the values of M_H2 and T_rot derived from our RD analysis. We corroborate that, in spite of strong CO line flux variations and confirmed changes in the relative level population, the total number of emitting molecules (adding up all levels) from the RD analysis remains constant (within uncertainties, < 15%) at multiple epochs. This is true for both the warm and hot components (both identified in IRC+10216) and, thus, for the single-fit value of M_H2.

We also find that line variability has also a small impact on the value of T_rot derived at multiple epochs in IRC+10216. The values of T_rot from the single-fit and warm component (using J_u < 27 transitions) remain essentially constant within uncertainties, consistent with negligible impact of non-LTE effects (if any). The highest-excitation CO lines (J_u > 27) show stronger variations and, as a consequence, the variations of T_rot for the hot component are larger, of ≲±15% about the average, maybe denoting some repercussion of LTE deviations on T_rot in this case.

Although it is improbable given our previous discussion, let us now consider whether the double slope in the RD diagram found for a few objects in our sample could predominantly result from non-LTE effects, with T_rot deviating most significantly from T_kin for the highest-J CO transitions, and what would be the implications in this case.

Among our targets (O-rich and S-type), we find a double-T_rot component in a total of five sources, covering a broad range of mass-loss rates, from Ṁ ∼ 2×10⁻⁷ to 8×10⁻⁵ M_⊙ yr⁻¹. The mass-loss rate distribution of these sources is the same as the sources for which only one T_rot component is found. Moreover, we do not find any correlation between Ṁ and the properties of the hot and warm components (also including C-rich targets, see e.g. Fig. A.1), which makes improbable the double T_rot component being predominantly a non-LTE effect (expected to be most remarkable in low mass-loss rate objects).

Moreover, to explain the observed double slope as a purely non-LTE effect, two conditions need to be simultaneously met in all double-T_rot targets. First, a very specific value of the number density, in a narrow range around ∼10⁷ cm⁻³, is needed to have J_u < 27 levels close to thermal equilibrium but J_u > 27 levels deviating perceivably from LTE (this would generate the double slope in the RD diagram). Second, a strong background radiation source T_bg > T_kin is required to reproduce a T_rot that is higher for J_u > 27 transitions than for J_u < 27 lines. Dust mixed with the gas in the CO-emitting volume and heated by the stellar radiation to 800 K could produce the required T_bg. If this is the case, then it would also imply similarly high gas temperatures in the emitting region since the difference between the gas and dust temperature across these inner wind layers is expected to be moderate (of < 100 K) and to follow a similar radial dependence based on detailed dust-continuum and non-LTE radiative transfer CO emission models (e.g. for W Aql, T_kin ∼ T_dust ∼ 300 K at r_CO = 2×10¹⁵ cm Danilovich et al. 2014). In the extreme case (e.g. maybe for very low mass-loss rate objects with densities < 10⁶ cm⁻³) that T_rot is totally decoupled from T_kin but very close to T_bg for all transitions from J_u = 14 to 40, then the double-slope of the RDs would most likely reflect the stratification of the dust temperature (rather than the gas) across the envelopes.

Fig. A.1. Ratio of the hot-to-warm rotational temperature vs. the mass-loss rates for objects with a double-Trot component: targets from this work (Paper I) are represented by open symbols, C-rich objects (from Paper II) are plotted using filled symbols.

Summarizing, a minor impact of small LTE deviations on the values of M_H2 and Ṁ is expected, even for the targets with the lowest mass-loss rates, Ṁ ≈ 10⁻⁷ M_⊙ yr⁻¹. The double slope of the RD is unlikely to be primarily due to large LTE deviations and points to a temperature stratification in the warm inner layers of the envelopes probed by PACS CO line emission. The temperature profile deduced is further discussed in Sect. 7.1.

Appendix B: The characteristic radius of the CO-emitting volume

In this appendix, we provide a more detailed description of the procedure that we have followed to decide on a “characteristic” radius of the CO-emitting volume (r_CO) to be used to compute, and correct for, the CO line opacities (Sect. 5.1) and to obtain an estimate of the mean mass-loss rate of our targets (Sect. 6.2).

Fig. B.1. Power-law radial distribution of the gas kinetic temperature estimated from non-LTE excitation and radiative transfer models of the CO emission in a number of AGB envelopes (indicated at the top of the box using the same colour for the temperature laws) from the literature. References: Van de Sande et al. (2018); Maercker et al. (2016); Khouri et al. (2014); Decin et al. (2010). The solid and dotted lines of R Dor and IK Tau are different temperature profiles derived by different authors (solid from Maercker et al. 2016, and dotted from Van de Sande et al. 2018 – R Dor and Decin et al. 2010 – IK Tau). The thick black line is an average temperature law.

For AGB CSEs (also including OH/IRs), a first estimate of the size of the CO-emitting volume is derived using the envelope temperature structure estimated from detailed non-LTE excitation and radiative tranfer calculations in the literature. These exist for many AGB CSEs using low-J CO transitions (typically J_u ≲ 6), and, for more reduced samples, also using certain CO transitions from higher J levels.

In Fig. B.1 we show temperature power-law radial profiles in the form

$$T(r)=T_i\times {\left(\frac{r}{R_i}\right)}^{\alpha }$$ that have been deduced in a number of previous works (references are given in the figure caption) for the few targets in our sample that include at least one high-J CO line observed with any of the three Herschel instruments (HIFI, SPIRE, or PACS). As we see in this figure, the full range of rotation temperatures deduced from our RD analysis, T_rot ∼ 200–700 K (Sect. 6), is consistent with the bulk of the PACS CO emission studied here arising from layers between ∼2×10¹⁴ and ∼5×10¹⁵ cm, if we take together all the T(r) distributions in Fig. B.1 (coloured lines), or between ∼3×10¹⁴ and ∼3×10¹⁵ cm, if we use an average T(r) power-law (thick black line, in Fig. B.1). The quadratic mean of these two ranges is r_CO ∼ 10¹⁵ cm, which we uniformly use as the reference value in most AGB stars to compute M_H2, T_rot, and Ṁ (but see more details below). Since the radius is, in any case, uncertain, we offer additional estimates of these parameters for a larger range of radius around the default r_CO (Appendix C). As we will see next, for the lowest (larger) mass-loss rate AGB stars in our sample, the characteristic CO-radius could be a factor of ∼2–4 smaller (larger) than the r_CO ∼ 10¹⁵ cm reference value.

In the following, we perform a case-by-case RD analyis of three targets (R Dor, W Aql, and TX Cam) that are representative of three main mass-loss regimes in our sample: low, medium and high (Ṁ ≈ M_⊙ yr⁻¹, ≈10⁻⁶ M_⊙ yr⁻¹, and ≳10⁻⁵ M_⊙ yr⁻¹, respectively) and for which detailed non-LTE excitation and radiative transfer studies of the CO emission including some mid-to-high J_u-CO transitions exist in the literature. The purpose is twofold: (1) to better illustrate our method to constrain r_CO in our sample (needed to decide on a plausible reference value for a uniform analysis), and (2) to compare the mass-loss rates from our approximate RD analysis with those computed using more sophisticated modelling tools.

R Dor. Precise CO-modelling studies of the molecular envelope of this object, including CO transitions up to J_u = 16 (E_up = 750 K), have been carried out by Maercker et al. (2016) and Van de Sande et al. (2018), references to other molecular line studies using lower-J_u transitions are given in Table D.2;); . The gas temperature structure deduced from these works, which use the same set of CO lines, are not fully coincident (dotted and solid red lines in Fig. B.1). Van de Sande et al. (2018) claim that their model achieves a better fit to more lines, and therefore, we adopt the T(r) power-law determined by these authors, with an exponent α = 0.65.

Fig. B.2. Rotational diagrams of the low mass-loss rate O-rich AGB star R Dor to illustrate the different steps that we have routinely followed to constrain rCO according to our temperature-opacity combined criteria (full details are given in the text, Appendix B).

Our general procedure to determine an optimal/acceptable value of r_CO is as follows. From a linear fit to the CO datapoints in the RD diagram (dotted line and black dots, respectively, in the top panel of Fig. B.2) we find a first or iteration-zero estimate of T_rot(i = 0) before any opacity correction is applied. Assuming that T_rot(i = 0) ∼ 570 K represents (or is close to) the average gas temperature within the CO-emitting volume and using the T(r) function depicted in Fig. B.1, we derive an initial value of r_CO(i = 0) ∼ 2.5×10¹⁴ cm. We use this radius to compute the opacity of the different transitions. In this case, we obtain τ_{14 → 13}∼1.8 for the CO J = 14 − 13 line, and smaller for higher-J lines (see Eq. (4)). After applying the opacity correction, through the frequency-dependent factor C_τ = τ/(1 − −e^−τ), a second linear fit to the opacity-corrected datapoints is performed (pink line and dots, respectively, in top panel of Fig. B.2), resulting in a new estimate of T_rot(i = 1). The latter always is smaller than the iteration-zero value because after applying the C_τ correction, the slope of the RD increases, as a direct result of the frequency-dependence of C_τ and the typical values of T_rot that we measure in our sample. Comparing again the value of T_rot(i = 1) ∼ 480 K with the T(r) law, we deduce a new value of r_CO(i = 1) ∼ 3.4×10¹⁴ cm, where the gas would have such a temperature.

Fig. B.3. Rotational diagram of W Aql and NML Tau (bottom). Left: results of the RD for rCO = 1×1015 cm, which results in a perfect match between the opacity- and temperature-radius (labelled as rCO and rT, respectively, in the plot) but results in an undesired large opacity correction, τ14 → 13> 1. Right: results of the RD for rCO = 1.4×1015 cm, which satisfies our τ14 → 13 < 1 criteria.

This exercise indicates that an intermediate value between r_CO(i = 0) and r_CO(i = 1) is needed to reach full consistency between the r_CO value used for the opacity correction and the one that reproduces the opacity-corrected T_rot derived (according to the T(r) adopted). In the particular case of R Dor, we find a perfect match between r_CO(i = 0) and r_CO(i = 1) for r_CO ∼ 3.15×10¹⁴ cm (leading to T_rot(i = 1) ∼ 500 K at this radius – Fig. B.1, middle panel).

An additional constrain to the radius is obtained from the fact that, as already explained in Sect. 5, the deepest layer traced by the observed CO emission must be such that τ< 1. In order to satisfy this moderate-opacity condition, a value of r_CO slightly larger than the “perfect-match” r_CO above is preferred. For R Dor, a value of r_CO ∼ 3.5×10¹⁴ cm satisfies the opacity criteria (τ_{14 → 13}< 1) and is, within uncertainties, still fully consistent with the radius where the final value of T_rot is attained, adopting the T(r) law from previous works – Fig. B.2, bottom panel). The mass-loss rate that results adopting this value of the radius r_CO ∼ 3.5×10¹⁴ cm, which guarantees line opacities lower than 1, is Ṁ ∼ 4×10⁻⁷ M_⊙ yr⁻¹.

In this and in most targets, in principle we cannot rule out a larger radius (leading to even lower line opacities), even if T_rot deviates from the gas temperature predicted at that radius based on the T(r) law adopted. This is because there is no guarantee that the adopted T(r) is a 100% accurate representation of the real temperature structure in these inner-wind regions and because T_rot may not be either a 100% accurate representation of the temperature of the gas (especially for low mass-loss rate objects). For R Dor, if we adopt the default radius of r_CO ∼ 1×10¹⁵ cm previously established as a reasonable characteristic radius for the overall AGB class, we deduce a temperature of T_rot ∼ 570 K and a mass-loss rate of Ṁ ∼ 1×10⁻⁷ M_⊙ yr⁻¹.

For R Dor, the mass-loss rates derived from Van de Sande et al. (2018) and Maercker et al. (2016) from their non-LTE CO excitation and radiative transfer analysis are very similar in both cases (∼2×10⁻⁷ M_⊙ yr⁻¹) and lie within the range of values estimated by us from our simple RD approach and the plausible range of r_CO explored, Ṁ ∼ [1–4]×10⁻⁷ M_⊙ yr⁻¹.

W Aql. Non-LTE excitation and radiative transfer of the CO emission, including CO transitions up to J_u = 25 (E_up = 1796 K), have been carried out for this S-type AGB star by Danilovich et al. (2014). The mass-loss rate obtained in this work is Ṁ ∼ 4×10⁻⁶ M_⊙ yr⁻¹ (as quoted in Table D.2, where we provide additional references to other molecular line studies using lower-J_u transitions). Danilovich et al. (2014) derive the gas and dust temperature profile in the 2×10¹⁴–10¹⁷ cm envelope layers (see their Fig. 3). Both the gas and dust temperatures are quite similar throughout the whole envelope, with the dust being slighly cooler than the gas below 2×10¹⁵ cm (i.e. the region sampled by the PACS CO data). In this region, the gas and dust temperatures do not differ by more than 50% (in particular, at ∼10¹⁵ cm, values of T_dust∼ 400 K and T_kin ∼ 500 K are expected).

In Fig. B.1 we reproduce the gas temperature structure deduced by Danilovich et al. (2014); the exponent of this power-law function is α = 0.65. We have used this profile to constrain r_CO following the same procedure explained in the previous paragraphs for R Dor. As we can see in the left panel of Fig. B.3, a radius of r_CO ∼ 1×10¹⁵ cm produces a perfect match between the final (opacity-corrected) temperature derived (T_rot ∼ 500 K) and the radius where this temperature is expected (r_T) according to Danilovich et al. (2014). This value of the radius, however, results in line opacities larger than desirable and a mass-loss rate of Ṁ ∼ 9×10⁻⁶ M_⊙ yr⁻¹. In this case, the radius that produces τ_{14 → 13}< 1 and simultaneously satisfies (within uncertainties) the temperature criteria is r_CO ∼ 1.4×10¹⁵ cm (right panel of Fig. B.3. This results in a somewhat smaller mass-loss rate of Ṁ ∼ 5×10⁻⁶ M_⊙ yr⁻¹, which is in excellent agreement with previous estimates.

IK Tau. The physical properties of the circumstellar gas, including the temperature, velocity and density structure, have been determined from the radiative transfer modelling of the multi-transitional (sub)millimetre CO line observations by Decin et al. (2010) and Maercker et al. (2016). The temperature distributions from these two works differ substantially (as for R Dor) and are shown in Fig. B.1 (solid and dotted cyan lines).

Adopting the α = 0.6 distribution derived by Decin et al. (2010), we find that a radius of r_CO ∼ 5×10¹⁴ cm produces a perfect match between the final (opacity-corrected) value of the temperature (T_rot = 420 K) and the radius r_T where this temperature is expected, however, it results in opaque lines (Fig. B.3, bottom-left). For IK Tau, a radius of r_CO ∼ 9×10¹⁴ cm or larger is needed to meet the τ_{14 → 13}< 1 criteria (Fig. B.3, bottom-rigth). This results in a somewhat smaller mass-loss rate of Ṁ ∼ 1×10⁻⁵ M_⊙ yr⁻¹, also in agreement with the range of values from previous estimates (Ṁ ∼ [0.2–3.5]×10⁻⁵ M_⊙ yr⁻¹, Table D.2).

Rest of the targets. We have followed this procedure to constrain the range of r_CO values that are plausible (based on our opacity- and temperature- combined criteria) for each target individually. Whenever a T(r) distribution was available in the literature (this happened for a few targets, see Fig. B.1) we used that specific law, and for the rest we used a mean T(r) distribution represented by the thick black line in Fig. B.1 as a guide.

Remarks on the temperature profiles. The T(r) profiles used and represented in Fig. B.1, although obtained from detailed non-LTE excitation and radiatiave transfer models of the CO emission, are not free from large uncertainties. This is denoted by the discrepancies between the temperature structure derived for the same targets when computed by different teams (notable cases are R Dor and IK Tau, in spite of these being among the most studied and best characterized AGB envelopes). Moreover, the power-law distribution is only an approximation to the real temperature structure. Indeed, the latter can deviate significantly from a simple power-law in some regions, in particular, T(r) can have multiple-slopes as well as local minima and maxima as shown in De Beck et al. (De Beck et al. (2010), see their Fig. 9, )) and (Maercker et al. (2016), see their Fig. 5,)).

Additionally, given the moderate range of J_u in the CO transitions predominantly modelled and published to date, the curves represented in Fig. B.1 are expected to be most reliable in the intermediate-to-outer envelope layers, beyond ∼10¹⁵ cm. Recently, Lombaert et al. (2016) performed non-LTE excitation and radiative transfer models of a sample of C-rich AGB CSEs using a selection of seven high-J CO transitions observed with PACS and SPIRE covering levels from J_u = 15 to 38, that is, higher than in the studies represented in Fig. B.1. One interesting result from this study is that the temperature of the gas falls off less steeply, with a power-law exponent close to α ∼ 0.4 − 0.5, in the inner-warm envelope regions sampled by these high-J CO transitions. This is indeed consistent with the crude estimate of α deduced in this work by comparing the hot and warm envelope components for a few targets in our sample (see Sect. 7). More dedicated non-LTE excitation and radiative transfer models of the high-J CO emission in larger samples is needed to confirm this result and to obtain a more precise description of the temperature distribution in these inner and warm regions of the envelopes of AGB stars.

Final remark on the reference values of r_CO. It is evident that the value of r_CO derived using a temperature criteria is very model-dependent since it fully relies on T(r) and on T_rot being a perfect representation of T_kin. For this reason, we prefer to be cautious and, after having explored carefully the range of plausible r_CO values based on a temperature-opacity combined criteria for all targets, to use a set of (order-of-magnitude) estimates of r_CO for the whole sample, rather than using the best-match radius individually found for each target. As shown in Table D.2, we chose a reference radius: r_CO = 1×10¹⁵ cm for most AGB CSEs, r_CO = [2–4]×10¹⁵ cm for high-Ṁ AGB CSEs (> 1×10⁻⁵ M_⊙ yr⁻¹) and r_CO = [4–6]×10¹⁴ cm for low-Ṁ AGB CSEs (< 1×10⁻⁶ M_⊙ yr⁻¹).

As discussed in Sect. 5, the temperature profile in the molecular envelopes of post-AGB objects and PN is extremely poorly characterized due to a much more limited number of CO emission cartography and, most importantly, modelling studies (especially using high-J CO transitions). Based on the systematically low values of T_rot in post-AGB objects and PNs found in this work (∼200 K) as well as their generally larger dimensions, we chose a reference value of the radius of r_CO = 6×10¹⁶ cm in these cases (Table D.3), although we explored a larger range of r_CO (examples are given in Figs. C.1 and C.2). This range of r_CO resulted in negligible RD opacity corrections in all these sources.

Appendix C: Additional opacity and mass-loss rate plots

In this appendix we show, for each source, the values of the optical depth (τ, defined in Sect. 5.1) of the different CO transitions used to build the RDs for a range of values of the characteristic radius (r_CO) of the CO-emitting volume (Fig. C.1). As we can see in this figure, the radius adopted implies a moderate value of τ and, therefore, a moderate opacity correction. The radius of the PACS CO-emitting region is the major source of uncertainty in the computation of τ and, therefore, of the emitting mass and mass-loss rates. It also has an effect on the rotational temperature. The values of M_H2 and T_rot for a range of r_CO are plotted inside the boxes in Fig. C.1. The resulting mass-loss rates, Ṁ, are shown in Fig. C.2, together with the broad range of mass-loss rates found in the literature for each object, typically derived from low-J CO line studies (references are provided in Table D.2).

Fig. C.1. Plots of τ vs. Eup corresponding to different linear radii. In each chart there is a subgraphic showing MH2 vs. Trot. The thickest curve and the biggest point are those associated to the rCO shown in Table D.3.

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.2. Plots of Ṁ vs. Trot. In each panel three different points, corresponding to three different rCO are shown. The point of the rCO shown in Table D.3 is bigger than the rest. The grey region represents the interval of mass-loss rates covered by previous works found in bibliography. For more details see Sect. 6. The linear radii explored are the same as the ones shown in Fig. C.1.

Fig. C.2. continued.

Appendix D: Tables

All Tables

All Figures

Fig. 1. IRAS colour–colour diagram (left) and 100 μm flux vs. [12]−[25] IRAS colour (right) of the THROES targets with good quality IRAS data (Quality Flag = 3) in the 12, 25, 60 and 100 μm bands. The colour–colour diagram is divided in different boxes where sources with common characteristics and evolutionary stages are located (see van der Veen & Habing 1988, for details). Objects with CO line detections, studied in this work, are highlighted using the symbol and colour code shown in the figure. In the rest of the plots presented throughout this paper, AGBs, OH/IRs, post-AGBs and PNe are coloured and symbol-coded as indicated in this figure.

In the text

Fig. 2. PACS continuum subtracted spectra of our sample in the 55–75 μm range. The 12CO rotational lines are indicated by vertical dotted lines (see Table D.4). The upper-level rotational quantum number Ju of a few selected CO transitions are shown as a reference. The sources are sorted from top to bottom in these five groups: AGBs, S-type stars, extreme OH/IR stars, post-AGBs, and PNe. Within each category, targets are ordered based on the stellar effective temperature from the coolest source to the hottest one.

In the text

	Fig. 3. The same as in Fig. 2 but covering the spectral range from 110 to 190 μm.
In the text

Fig. 4. Top: integrated flux of the CO J = 15−14 transition at 173.6 μm (FCO 15−14) against the IRAS 100 μm flux. Targets with no PACS data beyond ∼145 μm (i.e. not covering this CO transition; RR Aql and T Cep – see Fig. 3) or with bad quality (Quality Flag≠3) IRAS100 data (OH 26.5+0.6) are excluded. Middle: fluxes of the CO J = 15−14 transition and the adjacent continuum (PACS170). In this diagram, OH 26.5+0.6 is included. Bottom: line-to-continuum (F CO 15−14-to-PACS170) flux ratio against the IRAS colour [12]−[25].

In the text

Fig. 5. RDs of CO for our targets (continues in 5b and 5c). For each source, the values of Trot and MH2 derived from the fits (solid lines) as well as the characteristic radius of the emitting envelope layers adopted (RCO) are shown (Sect. 5). The errors of Trot and MH2 represent the uncertainty of the fit considering the error bars of the individual points, which include formal errors due to the noise and absolute flux uncertainty of the PACS spectra. For a few sources, a double component fit was necessary. For these cases, the values of Trot and MH2 for the so-called warm (blue) and hot (red) components are also indicated.

In the text

	Fig. 5. continued.
In the text

	Fig. 5. continued.
In the text

	Fig. 6. Values of Trot and MH2 from our CO rotational diagram analysis (single fit; Sect. 6) before and after applying the opacity correction (grey and coloured symbols, respectively).
In the text

Fig. 7. Relation between MH2 and Trotderived from our CO RD analysis and other magnitudes. Top-left: MH2 vs. Trot. For objects (a total of five) with a double-Trot component, both the single-fit and double-fit values are connected with dashed lines for an easier identification. Top-right: values of MH2 (single-fit) vs. vexp taken from the literature. Bottom panels: integrated fluxes of CO (J = 15−14) line vs. MH2 (left) and the same line normalized to its adjacent continuum vs. Trot (right).

In the text

	Fig. 8. Top: relation between Ṁ and Trot. Single- and double-fit (warm and hot) values, whenever these exist, are connected using dashed lines. Middle: relation between Ṁ and vexp. The distribution approaches a power-law of the type Ṁ ∝ vexp2.5 (dashed line). Bottom: integrated flux of the CO J = 15−14 line vs. Ṁ in a logarithmic scale.
In the text

Fig. 9. Ratio of the mass and rotational temperatures for the hot and warm components, $M_{{\text{H}}_{\text{2}}}^{\text{hot}}/M_{{\text{H}}_{\text{2}}}^{\text{warm}}\hspace{0.17em}\text{vs}.\hspace{0.17em}{T}_{\text{rot}}^{\text{hot}}/T_{\text{rot}}^{\text{warm}}$, for sources with a double-Trot component, namely, the AGB stars NML Tau, π Gru, W Aql, and χ Cyg and the yPN IRAS 17347-3139. The solid line is a fit by a power-law function (y ∼ x−1/α) to the AGB stars’ data – see text in Sect. 7.1. A similar trend is found in our analogue study of the THROES C-rich sample reported in Paper II.

In the text

Fig. 10. Comparison of the mass-loss rates from this work and from the literature (Tables D.2 and D.3, respectively). The solid line represents a 1:1 ratio. Target colour and symbol code as in Fig. 1. The big circles are used to locate a few targets with previous ṀBibliography estimates obtained from non-LTE excitation and radiative transfer analysis of CO lines, including some of the PACS transitions.

In the text

	Fig. A.1. Ratio of the hot-to-warm rotational temperature vs. the mass-loss rates for objects with a double-Trot component: targets from this work (Paper I) are represented by open symbols, C-rich objects (from Paper II) are plotted using filled symbols.
In the text

Fig. B.1. Power-law radial distribution of the gas kinetic temperature estimated from non-LTE excitation and radiative transfer models of the CO emission in a number of AGB envelopes (indicated at the top of the box using the same colour for the temperature laws) from the literature. References: Van de Sande et al. (2018); Maercker et al. (2016); Khouri et al. (2014); Decin et al. (2010). The solid and dotted lines of R Dor and IK Tau are different temperature profiles derived by different authors (solid from Maercker et al. 2016, and dotted from Van de Sande et al. 2018 – R Dor and Decin et al. 2010 – IK Tau). The thick black line is an average temperature law.

In the text

	Fig. B.2. Rotational diagrams of the low mass-loss rate O-rich AGB star R Dor to illustrate the different steps that we have routinely followed to constrain rCO according to our temperature-opacity combined criteria (full details are given in the text, Appendix B).
In the text

	Fig. B.3. Rotational diagram of W Aql and NML Tau (bottom). Left: results of the RD for rCO = 1×1015 cm, which results in a perfect match between the opacity- and temperature-radius (labelled as rCO and rT, respectively, in the plot) but results in an undesired large opacity correction, τ14 → 13> 1. Right: results of the RD for rCO = 1.4×1015 cm, which satisfies our τ14 → 13 < 1 criteria.
In the text

	Fig. C.1. Plots of τ vs. Eup corresponding to different linear radii. In each chart there is a subgraphic showing MH2 vs. Trot. The thickest curve and the biggest point are those associated to the rCO shown in Table D.3.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.2. Plots of Ṁ vs. Trot. In each panel three different points, corresponding to three different rCO are shown. The point of the rCO shown in Table D.3 is bigger than the rest. The grey region represents the interval of mass-loss rates covered by previous works found in bibliography. For more details see Sect. 6. The linear radii explored are the same as the ones shown in Fig. C.1.
In the text

	Fig. C.2. continued.
In the text