© ESO 2019

1. Introduction

The asymptotic giant branch (AGB) is a late evolutionary stage of low-to-intermediate mass stars (1 M_⊙ ≲M_⋆ ≲ 8 M_⊙) that is largely dominated by mass-loss processes. AGB stars can shed significant portions of their outer atmospheric layers in a dust-driven wind, with mass-loss rates of up to Ṁ ∼ 10⁻⁴ M_⊙ yr⁻¹ (e.g. Habing 1996; Höfner & Olofsson 2018). The material expelled by the central star (with very low effective temperatures of T_eff ∼2000 − 3000 K) forms a cool, dense circumstellar envelope (CSE) that is rich in dust grains and a large variety of molecules. After a significant decrease of the mass-loss rate, the AGB phase ends and the “star+CSE” system begins to evolve to the planetary nebula (PN) phase, at which the CSE is fully or almost fully ionized as a consequence of the much higher central star temperatures (T_eff ≈10⁴–10⁵ K) and more diluted envelopes. During the AGB-to-PN transition, shocks – resulting from the interaction between slow and fast winds at the end of the AGB phase (or early post-AGB) – also play an important role in changing the morphology, dynamics, and chemistry of the CSEs (e.g. Kwok 2000; van Winckel 2003; Bujarrabal 2006).

The CO¹ molecule is an excellent tracer of the CSEs of AGB stars, post-AGB objects, and PNe (e.g. Groenewegen et al. 1996; Schöier & Olofsson 2001; Schöier et al. 2002; Teyssier et al. 2006). The rotational transitions of the ground vibrational level over a wide range of excitation energies sample from cold (∼10 K) to hot gas (∼1000 K). Literature on observations of the cold, extended CO J = 1 − 0 to J = 6 − 5 emission around many evolved stars is abundant at millimetre/submillimetre wavelengths (e.g. Knapp et al. 1982; Knapp & Morris 1985; Knapp & Chang 1985; Olofsson et al. 1993; Bujarrabal et al. 1989; Jorissen & Knapp 1998; Castro-Carrizo et al. 2010; Sánchez Contreras & Sahai 2012; Ramstedt & Olofsson 2014). Studies looking at far-infrared (FIR) observations of even higher J CO transitions probing the warmest gas (∼200–1000 K) in deeper layers of CSEs are much more scarce. Pioneering works based on observations with the Infrared Space Observatory (e.g. Justtanont et al. 2000; Schöier et al. 2002) have continued more recently with Herschel (ISO; e.g. Groenewegen et al. 2011; Bujarrabal et al. 2012; Khouri et al. 2014; Danilovich et al. 2015; Nicolaes et al. 2018).

This is the second of a series of papers (Paper I Ramos-Medina et al. 2018a) in which we analyse Herschel/PACS spectra of a large sample of evolved stars from the THROES catalogue in a uniform and systematic way (Ramos-Medina et al. 2018b) to study their warm inner envelope regions using high-J CO transitions at FIR wavelengths. As in other previous studies, we divide our sample into O-rich and C-rich targets (Papers I and II, respectively) since these two major chemistry classes correspond to progenitor stars with different masses that follow somewhat different evolutionary paths and also have a dissimilar dust composition; both of these characteristics potentially affect the mass-loss process. We include targets with the following different evolutionary stages: AGB, post-AGB (or pre-PNe), and young planetary nebulae (yPNe).

The goal of our study (Papers I and II) is to obtain a first estimate of the average excitation temperature (T_rot) and mass (M_H2) of the warm envelope layers traced by the PACS CO lines in a uniform way using a simple analysis technique, the well-known rotational diagram (RD) method. The RD technique is useful to analyse large data sets rapidly (large number of lines and/or large samples) and to provide some constraints on these fundamental parameters. With the aim of benchmarking the results of the RD approximation, we obtain rough estimates of mass-loss rates (Ṁ) and compare these to values in the literature, paying particular attention to studies including at least a few high-J CO (FIR) transitions.

The impact of possible non-LTE effects on the results from the simple RD analysis was investigated in Paper I. It was concluded that, although they are expected to be minor in general and probably only affect the highest-J CO transitions studied in this work (J ≳ 27) at most, their existence cannot be ruled out in the lowest mass-loss rate stars and/or the outermost layers of the PACS CO-emitting volume. We also showed that even under non-LTE conditions, the masses derived from the RDs are approximately correct (or, at the very least, not affected by unusually large uncertainties) since the average excitation temperature describes rather precisely the molecular excitation (i.e., the real level population). This is also corroborated by the good agreement found between our estimates of Ṁ and those from detailed non-LTE excitation and radiative transfer (nLTEexRT) studies that exist for a number of targets. We stress that the RD method enables us to characterize the warm CSEs of evolved stars in a first approximation. For a more robust and detailed study of the radial structure of the density, temperature and velocity in the CSEs, as well as for establishing potential mass-loss rate variations with time, more sophisticated analysis is needed (e.g. Ryde et al. 1999; Schöier et al. 2002; De Beck et al. 2012).

Upon submission of this manuscript (and after our Paper I was accepted for publication) we became aware of a recent work by Nicolaes et al. (2018) who independently presented Herchel/PACS (and SPIRE) range spectroscopy of a sample of 37 AGB stars. These authors performed a similar RD analysis of the CO spectra and focused on deriving excitation temperatures; in contrast to our study, estimates of the envelope mass or mass-loss rates are not reported. Other differences with respect to the work by Nicolaes et al. (2018) is that we introduce a canonical opacity correction in the RDs and that we include post-AGBs and PNe.

2. Observations

2.1. Observations and data reduction

The PACS is a photometer and medium-resolution grating spectrometer (Poglitsch et al. 2010) on board the Herschel Space Telescope (Pilbratt et al. 2010) probing the FIR wavelength range. The PACS spectrometer covers the wavelength range from 51 μm to 210 μm in two different channels, which 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, structured in an array of 5 × 5 spatial pixels (“spaxels”) with 9.4″ × 9.4″. The PACS provides a resolving power between 5500 and 940, i.e., a spectral resolution of approximately 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). The technical details of the instrument can be found in the PACS Observer’s Manual².

The PACS (1D) spectral data were taken from the THROES (caTalogue of HeRschel Observations of Evolved Stars) website³, which contains fully reduced PACS spectra of a collection of 114 stars, mostly low-to-intermediate mass AGB stars, post-AGB, and PNe. The data reduction is explained in Ramos-Medina et al. (2018b) in full detail. The catalogue also consists of a compilation of previous photometry measurements at 12, 25, 60, 100 μm with the Infrared Astronomical Satellite (IRAS).

Table B.2 offers a description of the observations used in this work in which we provide the target name as listed in the header of the PACS FITS files and an alternative name for which some of the stars are more well known in the literature. We only analysed the spectra between [55–95] and [101–190] μm because the flux densities are unreliable above 190 μm, below 55 μm and between 95–101 μm due to spectral leakage. Two observation identifiers (OBSIDs) per target, corresponding to the bands B2A, B2B and R1, are necessary to cover the full PACS wavelength range. In the case of IRC+10216, band B2A data exists (see Decin et al. 2010), but these were acquired in a nonstandard spectroscopy mode and have restricted access in the Herschel Science Archive and, thus, are not included in the THROES catalogue. Instead, we used 4+3 OBSIDs corresponding to seven different operational days: OD = [745, 1087, 1257, and 1296] covering the spectral range [69–95, 140–190] μm and OD = [894, 1113, and 1288] covering a narrower interval [77–95, 155–190] μm.

2.2. Sample overview

We searched for CO emission amongst the entire THROES catalogue, but in this paper we focus on C-rich CSEs (29% of the entries) observed in PACS range mode. We found 15 evolved stars with at least three CO emission lines with signal-to-noise ratio above 3.

This sample contains bright infrared targets spanning a range of evolutionary stages from the AGB to the PN phase, but sharing similar carbon chemistry with strong CO emission at high excitation temperatures (up to E_u/k ≈ 5688 K). The AGBs, which are mostly Mira variables, are known for their high mass-loss rates compared to the mean value of Ṁ ∼ 1.5 × 10⁻⁷ M_⊙ yr⁻¹ derived from studies with large samples of carbon stars (Olofsson et al. 1993). We also include two mixed-chemistry post-AGBs (Red Rectangle and IRAS 16594−4656; Waters et al. 1998; Woods et al. 2005) and two mixed-chemistry yPNes (Hen 2-113 and CPD-56°8032; De Marco & Crowther 1998; Danehkar & Parker 2015, and references therein), that is, they show both C-rich and O-rich dust grains. In Table B.3 we provide a summary of some relevant properties such as distance, effective temperature, and gas mass-loss rate, along with additional references.

Figure 1 shows the classic IRAS colour-colour diagram (van der Veen & Habing 1988) featuring the colours [25]−[60]= and of all the stars in THROES. The stars studied are highlighted with coloured filled symbols, which we use consistently throughout the paper. This diagram is known to be a good indicator of the evolutionary stage of low-to-intermediate evolved stars; AGBs populate the lower left corner and more advanced stages are located on the diametrically opposed, so-called cold side of the diagram. It also shows an evolution in terms of the mass-loss rate and/or progressively increasing optical depths (Bedijn 1987). C-rich AGBs clearly constellate in a different box compared to the O-rich stars in Paper I, which has been interpreted as a consequence of emissivities of the different grains (Zuckerman & Dyck 1986). The AGB star that falls outside the expected box with a clear 25 μm excess ([12]−[25]=0.1) is AFGL 3068 (LL Peg), which is an “extreme carbon star”: very dust-obscured by optically thick shells due to high mass-loss rate (Volk et al. 1992; Winters et al. 1997). The post-AGB object HD 44179, best known as the Red Rectangle, is also outlying in this same box with respect to objects in a similar evolutionary status beyond the AGB, which typically show much larger 25 μm excess indicative of detached cold dust (and gas) envelopes.

Fig. 1. IRAS colour diagrams for the stars in the THROES catalogue. The colours are defined from the infrared fluxes at 12, 25, and 60 μm. The boxes in the left panel are those defined in van der Veen & Habing (1988), and the highlighted points correspond to the stars studied in this paper.

3. Observational results

3.1. Features in PACS spectra

The PACS continuum-subtracted spectra are plotted in Fig. 2, sorted in terms of stellar effective temperature (T_eff increases from top to bottom), and arbitrarily scaled for better comparison of the CO spectra. The continuum was fitted using a non-parametric method after identifying the line-free regions of the spectrum for each target. A brief line-detection statistics summary is given in Table B.4. Broad emission or absorption features due to instrumental artefacts are occasionally visible in some of the targets (e.g. near 62 μm). The spectra of IRC+10216 taken at seven different epochs are shown in the top panel of Fig. 3, and in more detail in Fig. C.2.

Fig. 2. Continuum-subtracted PACS spectra of our sample. The vertical segments (at the top) and bars indicate the rotational transitions of 12CO (light grey dotted lines), 13CO (dark grey), HCN (red), CS (yellow), OH (cyan), CH+ (green), and forbidden lines (blue) of [C II] 157.7 μm, [O I] 63.2, 145.5 μm, [N II] 121.9 μm.

Fig. 2. continued.

Fig. 3. Line and continuum variability of IRC+10216 in the FIR. Top panel: continuum-subtracted spectra at different epochs (see also Fig. C.2); the dotted lines represent the CO transitions. Bottom panel: sine-wave fit to the continuum variability at two selected wavelengths for a fixed pulsation period of 630 days (solid lines) and a free-period fit corresponding to 680 and 660 days for the blue and red dashed curves, respectively.

Straightaway we see significant differences between the spectra of these targets such as the larger line density, resulting from a richer molecular content, in AGB stars compared to post-AGBs and yPNe as expected. In AGB stars the strongest spectral lines are due to CO emission, while post-AGBs and yPNe show much more prominent forbidden lines, except for AFGL 2688 and AFGL 618. Additional molecular features with intense emission are attributed to rotational transitions of common C-bearing species such as ¹³CO, HCN, and CS.

Not so common in C-rich stars is the presence of OH, yet we identify several OH doublet lines (at 79, 84, 119, and 163 μm) in all yPNe and in the post-AGB star IRAS 16594−4656, which has the highest T_eff∼10 000 K in its class and where the presence of shocks is plausible. Two additional OH doublets at 65 and 71 μm, which have high-excitation energy of about 500 and 600 K, respectively, are also observed in IRAS 16594−4656 and in the yPNe Hen 2-113. In those five (T_eff > 10 000 K) targets, we also identify pure rotational transitions of CH⁺; the J = 3–2 line at 119.86 μm is one of the most clearly detected and best isolated lines⁴. The 179.61 μm CH⁺ line could be blended with the water 2₁₂-1₀₁ transition. The highest J lines of CH⁺ are also present in IRAS 16594−4656 and Hen 2-113. Interestingly, we do not see CH⁺ in the spectrum of the Red Rectangle, where sharp emission features near 4225 Å were discovered and assigned to this ion by Balm & Jura (1992). The detection of CH⁺ in AFGL 618 and the more evolved PNe NGC 7027 had been previously reported by Wesson et al. (2010, Herschel/SPIRE) and Cernicharo et al. (1997, ISO data). As we see below in this section, the targets with CH⁺ and OH emission also show the strongest atomic/ionic fine-structure lines.

We identify lines due to H₂O (orto- and para-) transitions in all AGBs, in the post-AGB AFGL 2688 and in the yPNe AFGL 618. The study of water lines in the C-rich AGB stars of our sample (except for AFGL 2513) was already conducted by Lombaert et al. (2016) who suggested that both shocks and UV photodissociation may play a role in warm H₂O formation. Water lines had also been previously identified in the SPIRE spectrum of AFGL 618 and AFGL 2688 Wesson et al. (2010).

The fine-structure lines [O I] at 63.18 μm and 145.53 μm and [C II] at 157.74 μm are very prominent in the spectra of all post-AGBs and yPNe in our sample with the exception of AFGL 2688. The post-AGB object IRAS 16594−4656 demonstrates clear signs of ionization since it shows a [C II] line even stronger than that of AFGL 618, which has a hotter central star (Table B.1). Weak [N II] 121.89 μm emission is also detected in Hen 2-113 and IRAS 21282+5050, and, tentatively in IRAS 16594−4656 and HD 44179.

The 69 μm band of forsterite (Mg₂SiO₄) is visible in CPD-56 ° 8032, Hen 2-113, and HD 44179, as already reported by Cohen et al. (2002) and Blommaert et al. (2014). The presence of O-rich dust in such C-rich objects corroborates a mixed chemistry nature.

3.2. Line flux variability in IRC+10216

Variability of the continuum at optical and IR wavelengths due to stellar pulsations is a common property of AGB stars. We present PACS spectra at multiple epochs of the Mira-type C-rich star IRC+10216 displaying periodic variations both in the continuum and in the lines (see Fig. 3 and Fig. C.2).

Cernicharo et al. (2014) first reported on the discovery of strong intensity variations in high-excitation lines of abundant molecular species towards IRC+10216 using Herschel/HIFI and IRAM 30 m data. Line variability was attributed to periodic changes in the IR pumping rates and also possibly in the dust and gas temperatures in the innermost layers of the CSE. From the analysis of a three-year-long monitoring of the molecular emission of IRC+10216 with Herschel (including HIFI, SPIRE, and PACS data), Teyssier et al. (2015) concluded that intensity changes of CO lines with rotational numbers up to J = 18 are within the typical instrument calibration uncertainties. But in the higher PACS frequency range (J ≥ 28), line strength variations of a factor ≳1.6 and scaling with J were found. More recently, He et al. (2017) reported 5%–30% intensity variability of additional millimetre lines with periodicities in the range 450–1180 days. In sect. 5.3 we study in greater detail the impact of CO line variability on the estimate of T_rot and M_H2 from the RD analysis using PACS data for different epochs.

3.3. CO line fluxes

We now focus on the purely rotational spectrum of CO in the ground vibrational state (v = 0), which we use to study the physical properties of the warm regions of the molecular CSEs of our targets. In the PACS range we can potentially find very high CO rotational transitions from J = 14 − 13 (E_u ∼581 K) to J = 45 − 44 (E_u ∼5688 K). The lack of data between 95 and 101 μm means we cannot detect the transitions J = 26–25 and J = 27–26. In the case of IRC +10216, the much larger gap in the PACS coverage (∼95–140 μm, Fig. 3) prevents detection of CO transitions with upper-level rotational number between J_u = 19 and J_u = 27.

The line fluxes are given in Table B.5. The quoted uncertainties correspond to the propagated statistical errors and do not contain absolute flux calibration errors (typically of ∼15%–20%) or underlying continuum subtraction uncertainties.

As expected, the resolving power of PACS (∼80–300 km s⁻¹) does not allow us to resolve the CO profiles spectrally. This is true not only for AGB CSEs with full width at half maximum (FWHM) similar to or smaller than the terminal expansion velocity of the envelopes (FWHM  ∼  10 − 25 km s⁻¹; Table B.3), but also for post-AGB objects and yPNe, even in targets that are known to have fast (≈100 km s⁻¹) molecular outflows like AFGL 618 (e.g. Bujarrabal et al. 2010). For this reason we measured CO fluxes by simply fitting a Gaussian function to the PACS lines.

Because of insufficient spectral resolution some reported line fluxes are affected by line blend. These are identified with asterisks in the tables and figures. Some of the well-known line blends are CO J = 30–29 with HCN J = 39–38 and CO J = 20–19 with HCN J = 26–25 at 87.2 and 130.4 μm, respectively. Also, the CO transitions J = 21–20 and J = 22–21 are blended with ¹³CO J = 22–21 and J = 23–22 at 124.2 μm and 118.6 μm, respectively.

Figure 4 compares the integrated flux of the CO (J = 15–14) line (F_CO (15−14)) 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. The J = 15 − 14 transition is a strong, non-blended line, detected towards all of our targets, and it was also used by Lombaert et al. (2016) who analysed PACS data of most of the AGBs studied in this work, which facilitates comparison.

Fig. 4. Correlation between CO line flux and continuum intensities. Top panel: CO J = 15 − 14 vs. IRAS 100 μm flux. Middle panel: CO J = 15 − 14 vs. PACS 170 μm continuum flux. Bottom panel: continuum normalised CO J = 15 − 14 flux vs. IRAS [12]–[25] colour. The symbols and colours are the same as in Fig. 1.

As for the O-rich sample (Paper I), there is a positive correlation between the CO line strength and the IRAS₁₀₀ and PACS₁₇₀ continuum fluxes. The correlation between the CO J = 1–0 line and the IRAS fluxes of evolved stars of various chemical types had been previously reported (Olofsson et al. 1987, 1988; Bujarrabal et al. 1992). We also confirm the anticorrelation between the line-to-continuum (F_CO 15 − 14/PACS₁₇₀) ratio and the IRAS [12]–[25] colour (both distance-independent) for the AGB stars, which was noted in Paper I. For the more evolved targets the trend is not so obvious, but the sample size is small. We also see that, in general, the ratio of the molecular emission to dust emission is higher in less evolved objects than in the most evolved objects, which could be partially attributed to more prominent CO photodissociation as the objects evolves along the AGB-to-PNe track. AFGL 618 and IRAS 16594−4656 are two clear outliers in this relation since they show a line-to-continuum emission ratio as large as that of the AGB class. The Red Rectangle (HD 44179) is well isolated in all the panels owing to its comparatively weak CO emission and low CO-to-dust ratio. This surely reflects the different nature of this object with respect to the rest of post-AGB and yPNe in our sample, which is well known from previous works. The Red Rectangle belongs to a special class of post-AGB objects with relatively weak CO emission coming from large (∼1000–2000 AU) circumbinary rotating disks; these disks have very prominent IR emission by warm/hot dust in the disk, but lack massive molecular outflows found in many other evolved stars (e.g. Bujarrabal et al. 2016, and references therein).

4. Rotational diagram analysis

Following the same approach of Paper I, we used the well-known and widely employed RD technique (e.g. Goldsmith & Langer 1999) to obtain a first estimate of the (average) excitation temperature (T_rot) and total mass (M_H2) of the warm inner layers of the molecular envelopes of our sample. A canonical opacity correction factor, C_τ, as defined by Goldsmith & Langer (1999), has been included to take into account moderate optical depth effects (Appendix A.1). We refer to Paper I for a more detailed description of the method.

4.1. Characteristic size of the CO-emitting layers

To compute the optical depth of the line and to make the corresponding C_τ correction, we need an estimate of the column density of CO, which is computed in a simplified manner dividing the total number of CO molecules (N_CO) by the projected area of the CO-emitting volume on the sky. The characteristic size of the envelope regions where the CO PACS emission is produced (r_CO) is one of the main sources of uncertainty since, except for a few targets (see below), these high-J CO-emitting layers are unresolved by PACS. Therefore, the value of r_CO needs to be adopted based on several criteria. These criteria are described in much detail in Paper I, Appendix B. In the following, we provide a brief summary of the general method and provide additional arguments particularised to the C-rich targets under study.

For AGB CSEs, a first estimate of the size of the CO-emitting volume can be derived from the envelope temperature structure, T(r), estimated from detailed non-LTE molecular excitation and radiative transfer (nLTEexRT) calculations in the literature. These have been done for many AGB CSEs using low-J CO transitions (typically J_u ≲ 6), and, for a few cases, also using certain CO transitions from higher J levels (see, e.g. Decin et al. 2010; Khouri et al. 2014; Danilovich et al. 2014; Maercker et al. 2016; Van de Sande et al. 2018). As deduced from these studies, the gas temperature is approximately 1000–2000 K close to the dust condensation radius (∼5–15 R_⋆), and decreases gradually towards the outermost layers approximately following a power law⁵ of the type ∼1/r^α, with α ∼ 0.5 − 1.0. Such models also exist in the literature for three targets in our C-rich sample: IRAS 15194−5115, CIT 6, and IRC+10216 (Ryde et al. 1999; Schöier et al. 2002; De Beck et al. 2012). According to this, the high-excitation transitions observed with PACS require relative proximity to the central star and, in particular, we anticipate that regions with a few hundred kelvin, as deduced from our RDs (Figs. 5 and 6), should not extend much farther than a few 10¹⁵ cm (≲10² R_⋆) in AGB stars.

Fig. 5. Rotational diagrams of the CO molecule. The grey line corresponds to a single least-squares fit to the full range of transitions from where a rotational temperature, Trot, and total gas mass, MH2, is computed. The characteristic radius of the CO-emitting volume (rCO) adopted is indicated. The red and blue lines correspond to a two-component model consisting of a “warm” and “hot” region, respectively. Asterisks indicate line blends that were excluded in the fit.

Fig. 5. continued.

An additional constraint on the radius can be imposed from the fact that the deepest layer traced by the observed CO emission must be such that τ <  1 because of the almost null escape probability from deeper, very optically thick regions (Appendix A.1). For all our AGB CSEs, we explored a range of radii around 1 × 10¹⁵ cm (Fig. A.1). We found that values around r_CO ∼ [1−4] × 10¹⁵ cm result in line optical depths close to, but smaller than unity (typically τ_{J = 14 → 13} ∼ 0.5–0.9), which yield moderate C_τ opacity correction factors for lines with J_u  <   19 and negligible for higher J transitions. For r_CO < 1 × 10¹⁵ cm, the opacity of the CO J = 14 − 13 line, which is the optically thickest transition in our sample, becomes larger than 1 in all our targets (also including post-AGBs and yPNes).

We checked that the range of plausible radii found for the AGB CSEs in our sample, r_CO ∼ [1−4] × 10¹⁵ cm, is consistent with the upper limits to the size of their envelopes deduced from the PACS spectral cubes and/or photometric maps by da Silva Santos (2016) and with any other information on the molecular envelope extent from the literature. IRC+10216 is the only AGB CSE that is partially resolved in the PACS cubes, and it is also the closest in our sample. The PACS cubes show a slightly extended source (with a half-intensity size 4σ above the instrument PSF in both bands) that implies a deconvolved Gaussian radius of about 2″, that is ∼4 × 10¹⁵ cm at d = 150 pc. This value is in good agreement with the lower limit to r_CO needed to satisfy the τ_{J = 14 → 13}< 1 criteria in this target, r_CO ∼ 2 × 10¹⁵ cm. In this case, we then rather confidently use an intermediate value of r_CO = 3 × 10¹⁵ cm.

Contrary to AGBs, for post-AGBs and yPNe there are no model temperature profiles in the literature of the molecular gas in the CSEs (except for the rotating, circumbinary disk of the Red Rectangle, Sect. 6.1.2). The range of representative radius adopted for post-AGBs and yPNe is r_CO = [0.4−4] × 10¹⁶ cm (Table B.1) based on a moderate opacity criteria, the extent of the emission in the PACS cubes and photometric maps, and on additional information on the extent of the intermediate-to-outer molecular envelope from the literature. In particular, in AFGL 2688, a deconvolved diameter of ∼4″ in the PACS blue band suggests that r_CO ≲ 1 × 10¹⁶ cm (at d = 340 pc). Previous CO J = 2–1 mapping observations identified a compact shell of radius ∼2″ (∼5 × 10¹⁵ cm) around the centre of the nebula (Cox et al. 2000). Since the CO J = 14–13 emission is optically thick at r_CO ≲ 6 × 10¹⁵ cm (as shown in Fig. A.1), we adopt as representative radius an intermediate value of r_CO = 8 × 10¹⁵ cm. HD 44179 (the Red Rectangle) is a point source in the PACS photometric maps, therefore, for a distance of 710 pc, r_CO should be on the same order as in AFGL 2688, which is roughly consistent with interferometric observations (e.g. Bujarrabal et al. 2016). For IRAS 16594−4656, only a very loose upper limit to the radius of r_CO < 3 × 10¹⁶ cm is inferred from optical images and H₂ emission maps in this object (e.g. Hrivnak et al. 2008). We explored the range, r_CO ∼ 0.8 − 2 × 10¹⁶, similar to AFGL 2688, and adopted as a reference value the midpoint value where τ_{J = 14 → 13} ∼ 0.7.

All yPNe in our sample are point sources in the PACS spectral cubes, which means that the upper limit to the radius is of about r_CO ≲ 1″ ∼[1−2] × 10¹⁶ cm for AFGL 618, and ∼[4 − 6] × 10¹⁶ cm for the rest. They are known to have central H II regions that have recently formed as the star has become progressively hotter along the PNe evolution. Because the CO envelope surrounds the ionized nebula, a lower limit to r_CO can be established from the extent of the latter. Taking this into account, we set a representative radius to r_CO = 1 × 10¹⁶ cm for AFGL 618 (Sánchez Contreras et al. 2017; Lee et al. 2013) and r_CO = 4 × 10¹⁶ cm for the rest (see, e.g. Danehkar & Parker 2015; Castro-Carrizo et al. 2010).

Given the uncertainty in r_CO, we systematically explored a range of radii around optimal/plausible values of r_CO to assess the impact of this parameter in our results (Fig. A.1). The opacity correction increases the smaller r_CO is, therefore the slope and y-intercept of the RD increase as a direct result of the frequency-dependence of C_τ (Appendix A.1), which results in lower values of T_rot and higher values of N_CO, thus M_H2. This also means that, in practice, only the lowest frequency points are affected by C_τ, while the highest frequency transitions (i.e., highest excitation energies) are unaltered regardless of the radius that we chose within the reasonable constraints that we placed.

4.2. Non-LTE effects

In Paper I we examine and discuss extensively the impact of non-LTE excitation effects (if present) on the values of T_rot and M_H2 derived from the RDs in a sample of 26 non-C-rich evolved stars with mass-loss rates in the range Ṁ ∼ 2 × 10⁻⁷–1× 10⁻⁴ M_⊙ yr⁻¹. The C-rich targets studied in this work have on average higher mass-loss rates than those in Paper I. Therefore, using similar reasoning, the CO population levels are also most likely close to thermalization in the inner dense regions of the CSEs under study.

This is further supported by nLTEexRT computations of a selection of high-J CO transitions observed with PACS (from J_u = 14 to 38) by Lombaert et al. (2016). Their sample included all our targets except for AFGL 2513 and IRC+10216. These authors conclude that over a broad range of mass-loss rates (Ṁ ∼ 10⁻⁷–2 × 10⁻⁵ M_⊙ yr⁻¹) the CO molecule is predominantly excited through collisions with H₂; a minor effect of FIR radiative pumping is from the dust radiation field. The role of dust excitation on FIR CO lines was also investigated by Schöier et al. (2002) and found to be of minor importance for AGBs with typical mass-loss rates of ∼ 10₋₅ M_⊙ yr⁻¹. We assess the FIR pumping effect further in Sect. 5.3, where we study multi-epoch RDs of IRC+10216, which is a source with well-known CO line variability.

We stress that even under non-LTE conditions, for a simple diatomic molecule such as CO, the RD method provides a reliable measure of the total mass within the emitting volume. Although T_rot may deviate from the kinetic temperature in regions where the local density is lower than the critical densities of the transitions considered (n_crit ∼ 5 × 10⁵–3 × 10⁶ cm₋₃, for J_u = 14 and 27, respectively, and n_crit ≈10⁶ − 10⁷ cm₋₃ for J_u > 27), it describes very precisely the molecular excitation, i.e., the real level population. Therefore, the total number of emitting molecules (and, thus, the mass) is robustly computed by adding up the populations of all levels. This is also supported by the good agreement (within uncertainties) between the mass-loss rates derived from this (and other works) using similar LTE approximations, and those obtained from nLTEexRT models (including PACS lines for a few targets); see Paper I and Sect. 6.2.

As discussed in Paper I, conceivable LTE deviations would have its largest impact on the excitation temperature of the hot component, since high-J levels have the lowest critical densities (see above). If this is the case, in low mass-loss rate objects, the value of T_rot deduced for the hot component could deviate from the temperature of the gas and approach that of the dust within the CO-emitting volume. We note that, in any case, the gas and dust temperatures, although not equal, are not excessively divergent in the warm envelope regions around ∼10¹⁵ cm under study (e.g. Danilovich et al. 2014; Schöier et al. 2002).

5. Results

5.1. Gas temperatures and masses

The opacity-corrected RDs of the CO molecule are plotted in Figs. 5 and 6, where multi-epoch RDs for IRC+10216 are separately shown, together with the best-fit parameters using a single-temperature component (grey line) and double-temperature component (red and blue lines). For 8 out of 14 sources, the RDs include CO transitions with upper-level energies that range from E_u ∼580 K to E_u ≲2000 − 3000 K. For 6 targets, CO transitions with upper-level energies of up to E_u ∼5000 K are also detected. Our temperatures for the group of AGBs are consistent with those reported in Nicolaes et al. (2018) within the uncertainties, although for many AGBs we find lower values because of the opacity correction term that we included.

Fig. 6. Rotational diagram of the CO molecule in IRC+10216 at different epochs. Analogous to Fig. 5.

There are cases for which it is clear by the analysis of residuals that a single straight line does not fit the entire range of excitation energies. That is the case of V Hya, IRAS 16594−4656, AFGL 3116, CIT 6, and AFGL 618. Because it is not always clear by eye whether the slope changes and at which point it occurs, we used the Bayesian information criterion (BIC) to help us deciding where to split the diagram and to quantify significance. This is explained in the Appendix A.2 and illustrated in the supplementary Fig. A.2.

We provide the fitting parameters in Table B.1 corresponding to a single-fit and a double-component fit in targets where a single line does not equally fit all data points. As in Paper I, we call these two components “warm” and “hot” with T_w <  T_single <  T_h. Their mean temperatures are K and K, respectively. The corresponding masses are M_w and M_h and the former is 4–10 times larger than the latter. We find single-fit rotational temperatures in the range T_rot ∼ 200–700 K, where some of the post-AGBs and yPNe are the targets with the coolest gas, except for the yPNe AFGL 618 which has the largest rotational temperature in our sample similar to that of the warmest AGB CSEs.

The total number of CO molecules is in the range N_CO ∼ 10⁴⁹ − 10⁵¹, resulting in column densities of ∼ 10¹⁶− 10¹⁹ cm⁻² for the adopted radii (Table B.1). To estimate the total gas mass from CO we assumed the same fractional abundance X_CO = 8 × 10⁻⁴ (e.g. Teyssier et al. 2006) with respect to H₂ for all targets. The single-fit values of the total mass of the CO-emitting volume range between M_H2 ∼ 3 × 10⁻⁵ M_⊙ (V Cyg and HD 44179) and ∼8 × 10⁻³ M_⊙ (IRAS 21282+5050), with a median value of M_H2 ∼ 4 × 10⁻⁴ M_⊙.

Figure 7 shows the single-fit temperature versus mass for the opacity-corrected diagram (coloured symbols) and uncorrected (grey symbols). The opacity correction results in changes in T_rot of 10–15% in AGBs, and lower than 5% in post-AGBs and yPNe. In mass this typically corresponds to 60% in AGBs, lower in the post-AGBs and AFGL 618 (< 10% in M_H2), and negligible in the other three yPNes. Figure A.1 shows that τ_{J = 14 − 13} is close to unity but it quickly falls off with increasing J meaning that C_τ → 0 for E_u/k >  2000 K, and that M_h, in particular, is not underestimated by opacity effects.

Fig. 7. Summary of the RD results. Top panel: coloured symbols correspond to the opacity corrected values (single fit), which are connected by dotted lines to the corresponding uncorrected values in grey; on the sides we show the histograms of temperature and mass of the single fits; in the case of IRC+10216 we show weighted averages of four observations. Bottom panel: temperature and mass of the warm and hot components are connected by dashed lines for the same target (if applicable).

As in Paper I, we find an anticorrelation between M_H2 and T_rot, especially if we only consider the group of AGBs (Fig. 7). In general, post-AGBs and yPNe have the highest masses. One clear exception to trend is the Red Rectangle, which is one of the least massive targets (with only a few ∼10⁻⁵ M_⊙ yr⁻¹) in contrast to the rest of post-AGBs and yPNe. This is not surprising given the nature of this object, which is the prototype of a special class of post-AGB objects with hot rotating disks and tenuous winds that are very different from the massive and fast (high-momentum) outflows of standard pre-PNe (Bujarrabal et al. 2016, and references therein).

We also investigated the correlation between the CO J = 15 − 14 flux and the gas mass and find that the strongest CO emitters have tendentially more massive envelopes, although there is a significant scatter (Fig. C.1). Also, the targets with the highest temperatures have the highest line-to-continuum (F_CO 15 − 14/PACS₁₇₀) ratios. IRAS 21282+5050, which has the most massive warm envelope in our sample, defies this trend since it shows relatively strong CO J = 15 − 14 emission, although its T_rot ∼ 170 K is the lowest in the sample. The Red Rectangle appears isolated in a region of rather weak CO emission in spite of relatively high temperatures (∼400–500 K). We compare these results to Paper I in Sect. 6.3.

5.2. Mass-loss rates

The mass-loss rates are estimated by simply dividing the total mass by the crossing time of the CO-emitting layers, that is,

(1)

where v_exp is the expansion velocity of the gas, which has been taken from literature (Table B.3). For the characteristic radius of the CO-emitting region we use the same value (r_CO) adopted for the opacity correction. We note that this estimate represents a mean or “equivalent” mass-loss rate assuming constant-velocity, spherically symmetric mass loss during the time when the warm inner envelope layers where the CO PACS lines arise were ejected, that is, during the last ∼20–50 yr and ≲300 yr for AGBs and post-AGBs/yPNe, respectively, given the values of v_exp and r_CO.

The mass-loss rates are listed in Table B.1. We find a range of values of Ṁ ∼ 10⁻⁷−10⁻⁵ M_⊙ yr⁻¹, and a median value in our AGB stars of Ṁ ∼ 6 × 10⁻⁶ M_⊙ yr⁻¹. These values are not to be taken as representative of the whole class of C-rich evolved stars, since our sample is small and not unbiased. This is because the objects in the THROES catalogue were originally selected for Herschel observations for various reasons, probably including their strong CO emission.

As in Paper I, we investigated a possible correlation between Ṁ, and T_rot, v_exp. We see little evidence of an anticorrelation between Ṁ and T_rot, although the relation is strongly influenced by AFGL 618, which is a strong outlier in this parameter space (Fig. C.1). In this, and maybe other objects (mainly post-AGB/yPNe), we expect departures from the simple (constant mass-loss rate, spherically symmetric) model adopted to estimate the “equivalent” Ṁ. We compare the results obtained in this work and in Paper I in Sect. 6.3.

Figure 8 shows the logarithm of the integrated flux of the CO J = 15 − 14 line versus the logarithm of the mass-loss rate of the single component fit. The upper and lower limit of the error bars in Ṁ correspond to a range of radii (see Fig. A.1) around the representative one (see Table B.1). We find a positive trend that is consistent with a power-law relation similar to that found by Lombaert et al. (2016) in their sample of C-rich AGB CSEs with H₂O FIR emission lines.

Fig. 8. Logarithm of the mass-loss rate vs. total CO J = 15 − 14 flux.

Separate values of Ṁ for the hot and warm components are computed for completeness, but the difference found (M_h <  M_w) should not be overinterpreted as a recent decrease of the mass-loss rate. The hot and warm components most likely trace adjacent layers of the inner winds of our targets; the hot component presumably best samples regions closer to the center. However, for simplicity and since we ignore the true CO excitation structure, we use the same radius to formally compute Ṁ for both components. We note that owing to the 1/r_CO dependence, the values of M_h and M_w can be brought closer to the single-fit value if the warm and hot correspond to different r_CO. We refrain from discussing M_h and M_w separately since a more sophisticated analysis, including nLTEexRT modeling, is needed to assess mass-loss time variability. For this reason we only compare our single-component mass-loss rates to the literature (Sect. 6.2).

5.3. Influence of line variability on T_rot and M_H2

We have shown in Fig. 3 the temporal variability of the continuum and line fluxes in the case of the Mira-type variable AGB star IRC+10216. The higher J CO lines (E_u > 2000 K) are those that show the strongest variations with time. In this case, we are interested in studying how CO line variability affects the values of M_H2 and T_rot derived using the RD method but using data acquired at different epochs. We use the same seven available OBSIDS⁶ of IRC+10216 as Teyssier et al. (2015), which span a time period of 551 days (Table B.2).

The RDs of IRC+10216 for these seven different observing epochs are shown in Fig. 6, and the results of the RD analysis are tabulated in Table B.1, together with the remaining targets. Since three of the OBSIDs corresponding to the ODs 894, 1133, and 1288 have a more restricted wavelength coverage, the fits to the RDs have an inherently larger uncertainty because the fit is more sensitive to the low number statistics. The error-weighted mean (single-fit) rotational temperature and mass are T_rot∼520 K and M_H2 ∼ 4 × 10⁻⁴ M_⊙, respectively.

In Fig. 9, we plot the temperature and mass for single- and double-temperature components versus the operational day of the observations. The bottom panel shows that the total gas mass deduced from the single-fit of the RD or for the warm and hot components does not reflect the line flux variability since it stays essentially constant with time about the average value (dotted lines), well within the estimated uncertainties.

Fig. 9. Rotational temperature and mass vs. operational day in IRC+10216. Top panel: sinusoidal fit to the Trotvariation with fixed period of 630 days (solid line). Bottom panel: total gas mass over time. In each panel, the dotted lines indicate the average of each component and the asterisks denote unused data points in the fit (see text). The colour code is the same as in Fig. 6.

In a similar manner, the temperature of the warm component does not clearly reflect the CO line flux variations, since all multi-epoch values are in good agreement within uncertainties. The hot component shows the largest variations (perhaps periodic) of the temperature, where T_h is ∼100 K (∼16%) above the average (∼ 640 K) at OD 1087 and then relaxes back to normal values in the remaining epochs. This variation is echoed in the single-fit value of T_rot (∼12%).

Temperature variations are not necessarily expected to be periodic, but if we assume that to be true by fitting a sinusoidal function with the known pulsation period P = 630 days to the RD parameters, we find that such model accommodates reasonably well the data points corresponding to T_h (and also the single-fit T_rot). It is also possible that IRC+10216 underwent an abrupt change of the physical conditions in its inner wind layers at epoch OD 1087, since the remaining data points by themselves do not justify/indicate periodic variability.

In summary, probably because of a compensation between the changes in the y-intercept and the slope of the RD from CO line flux variations, which are largest for transitions with the highest J, the total gas mass M_H2 appears to be robust to FIR pumping (non-LTE) effects. These, however, could have a measurable, yet moderate impact on the rotational temperatures.

6. Discussion

6.1. Gas temperatures

The detection of high-J CO rotational lines is an indication of a significant amount of molecular gas under relatively high temperature conditions. From our simple RD analysis we inferred that the average gas temperatures of the layers sampled by FIR CO lines are much larger (T_rot ∼ 200 − 900 K) than those typically derived from millimetre/submilimetre observations, which are sensitive to ≲100 K gas from the intermediate-to-outer layers of the envelopes of evolved stars (at ≈10¹⁶–10¹⁷ cm; see, e.g. De Beck et al. 2010; Schöier et al. 2011).

In a number of targets we identified a double-temperature (“warm” and “hot”) component. Deviations from a single straight line fit to the RD have also been found in some of the O-rich objects in Paper I and in other previous works using, for example, ISO and/or Herschel/SPIRE CO spectra in a number of AGB and post-AGB CSEs (e.g. Justtanont et al. 2000; Wesson et al. 2010; Matsuura et al. 2014; Cernicharo et al. 2015b; Cordiner et al. 2016). AFGL 618 and IRAS 16594−4656 are the targets whose RDs show the most obvious departure from linearity. IRAS 16594−4656 is particularly interesting since we found a much cooler warm component of just T_w ∼ 220 K, and a breakpoint at lower energies (E_u/k ≈ 1049.9 K) compared to other targets typically with T_w ∼ 450 K up to E_u/k ≈ 1794.3 K.

6.1.1. Temperature stratification

As in Paper I, in order to investigate if the double-temperature component is consistent with resulting from the temperature stratification within the inner layers of the CSEs, we compared the hot-to-warm M_H2 and T_rot ratios (Fig. 10), and find that they are correlated. If the temperature profiles in the envelope follow a power law of the type T(r)∝r^−α, where α is a constant, then the trend in Fig. 10 should also follow approximately a power-law function. We find α ∼ 0.4, which is similar to the value found for the O-rich targets studied in Paper I. The value of α is in agreement with past works that suggested that the kinetic temperature distribution is shallower, with values of α down to ∼0.4–0.5, for the inner (∼5 × 10¹⁴–3 × 10¹⁵ cm) CSE layers (De Beck et al. 2012; Lombaert et al. 2016; Matsuura et al. 2014) than for the outer regions, where the steepest temperature variations (α ∼ 1–1.2) are found (≳10¹⁶ cm; Teyssier et al. 2006).

Fig. 10. Hot/warm temperature ratio vs. hot/warm mass ratio for AGBs and post-AGBs and the power-law index. Top panel: lines correspond to power-law fits with index α = 0.4 for the opacity corrected RDs (solid) and α = 0.44 for the uncorrected RDs (dashed); the open symbols correspond to the O-rich AGB stars in Paper I. Bottom panel: power-law index for each individual target vs. mass-loss rate (single fit).

In addition to the average power-law exponent obtained by fitting all the targets simultaneously, we can also derive a value for each individual case applying

(2)

as shown in the bottom panel of Fig. 10. In the case of IRC+10216, we find α = 0.45 ± 0.06, which is in good agreement with the power-law exponent in the inner CSE between 9 and 65 stellar radii (i.e., up to ∼3 × 10¹⁵ cm) deduced from detailed non-LTE excitation and radiative transfer models (De Beck et al. 2012).

Therefore, the empirical relation found between the hot-to-warm ratio of M_H2 and T_rot is consistent with the double-T_rot component in some of our targets stemming (at least partially) from the temperature stratification across the inner envelope layers. The two components in the RD do not necessarily imply two distinct/detached shells of gas at different temperatures, but they most likely reflect the temperature decay laws. As explained in Paper I, in case of LTE deviations (not impossible in the lowest mass-loss rate objects), the value of α obtained from this simple approach would more closely represent the dust (rather than the gas) temperature distribution. This requires confirmation by detailed nLTEexRT models to the individual targets, which will be performed in a future publication.

6.1.2. Individual targets: Comparison with previous works

It is not possible to directly compare most of our results with literature because past studies have been focusing on the cold, outer components of CSEs. Prior to Herschel there was a study based on ISO LWS data in roughly the same wavelength range by Justtanont et al. (2000) who also performed RD analysis (without opacity correction). They found T_rot ∼ 700(±90) K and T_rot ∼ 380(±30) K for AFGL 618 and AFGL 2688, respectively. We obtained similar results without opacity correction, but introducing this effect lowered these values to T_rot ∼ 690(±30) K and T_rot ∼ 300(±20) K. In AFGL 618 the central star is hot enough (T_eff ∼ 33 000 K) to produce FUV photons that heat the gas. In AGFL 2688 this is probably not the case since the central star is much cooler (T_eff ∼ 7250 K), so low-velocity shocks are the most likely heating mechanism.

Also partially based on ISO data, Ryde et al. (1999) used a number of spectral lines of CO between J = 1 − 0 and J = 21 − 20 to infer the kinetic temperature profile across the CSE of IRAS 15194−5115. According to their model, T_kin ∼ 400 K at r_CO ∼ (1 − 2)×10¹⁵ cm, which is in excellent agreement with our opacity-corrected single component T_rot ∼ 410 K for r_CO = 2 × 10¹⁵ cm. This is because, as already pointed out by Ryde et al. (1999), these high-J levels are mainly populated by collisions, therefore they are proxies for the kinetic temperature, at least out to a few ∼10¹⁵ cm. Also using ISO data, Schöier et al. (2002) presented a kinetic temperature model for CIT 6 that shows that T_kin ∼ 400 − 500 K at approximately r_CO ∼ (1 − 2)×10¹⁵ cm, which is consistent with the warm component (T_w ∼ 460 K) that we infer and the representative radius r_CO = 2 × 10¹⁵ cm adopted. The hot component (T_w ∼ 830 K) that we found would imply that regions closer to the star have an important contribution to the emission of the highest J lines.

The only star whose inner/warmer (gaseous) CSE had been studied in detail before using Herschel/PACS data is IRC+10216. This was done by Decin et al. (2010) using high-J CO spectral lines to infer the kinetic temperature profile as a function of radial distance. They found T_kin ∼ 500 − 600 K at r_CO ∼ (1 − 2)×10¹⁵ cm. In this case, we applied the opacity correction for r_CO = 3 × 10¹⁵ cm, which seems to be the layer at which T_kin ∼ 300 K in their model. This is also the average value of T_w, which we found by fitting only the lowest J transitions; this value does not change appreciably with time despite strong line flux variability (Fig. 9). The hot component (∼600–700 K) may correspond to r_CO ∼ 1 × 10¹⁵ cm according to their model.

Using Herschel/SPIRE data, Wesson et al. (2010) also performed RD analysis of CO spectra and derived T_rot ∼ 70–230 K and T_rot ∼ 100–200 K for AFGL 618 and AFGL 2688, respectively. The rotational temperatures we find are thus higher than those obtained in Wesson et al. (2010) as expected. Further large-velocity gradient (LVG) calculations by those authors suggested that hot material at approximately 1000 K might exist at the high-velocity wind region of AFGL 618, which could be that traced by PACS since we found T_h ∼ 900 K.

For HD 44179 (the Red Rectangle) we obtained T_rot ∼ 440 K, which is about two to three times larger than the range of values inferred by Bujarrabal & Alcolea (2013) from Herschel/HIFI lower J CO observations. It is not surprising that we find a higher value since the higher J lines probed by PACS are probably formed deeper inside the rotating circum-binary disk at the core of this object. Meanwhile follow-up analysis has more clearly shown an outflow with ∼500 K (Bujarrabal et al. 2016), but also seems that such temperature conditions could exist in a region of the inner disk with a radius of about r_CO ∼ 2 × 10¹⁵ cm, unresolved by PACS. The opacity correction would still be moderate for this value (τ_{J = 14 − 13} ∼ 0.67), and it lowers the rotational temperature to T_rot ≈ 410 K which is still within the uncertainties. For the remaining yPNes (IRAS 21282+5050, CPD-568032, Hen 2-113) there are no kinetic temperature models in the literature to which we could compare our results.

6.2. Mass-loss rate

As in Paper I, we compared the values of the mass-loss rates derived from our simple RD analysis with other values found in the literature mostly from low-J observations. We paid special attention to a few targets with detailed non-LTE excitation and radiative transfer analysis of CO data including at least some high-J transitions observed with Herschel. This is also a way of ascertaining the robustness of the RD method.

Figure 11 shows our estimate of the mass-loss rate versus that found in the literature (see Table B.3). For each target the markers correspond to the single temperature component for the radius mentioned in Table B.1, and the error bars represent the uncertainty in the radius for the adopted v_exp. We show the computed opacities for the considered radii in each target in the supplementary Fig. A.1. The range of values found in the literature are shown by the grey shaded area whose bounds are set by the maximum and minimum Ṁ plus uncertainties when reported (factor of ∼3 in AGBs, Ramstedt et al. 2008; De Beck et al. 2010). Those values have been rescaled to the same distance, v_exp, and X_CO adopted in this work.

Fig. 11. Comparison between the mass-loss rate in this work and the literature. The markers correspond to the representative radii listed in Table B.1 for each target and the error bars are mass-loss rates for a given range of radii. The shaded area delimits the range of values found in the literature scaled to the same parameters assumed in this work (see text).

Similar to the non-C-rich THROES targets in Paper I, our mass-loss rates are in good agreement with values in the literature within the large uncertainties. In our case, these are dominated by the uncertainty in r_CO. We see that the range of radii we explored yields a Ṁ that fits within the shaded area. In many cases, the error bars are truncated at the upper limit above which the line opacities would be too large to allow reliable estimates of the masses and mass-loss rates (see Fig. A.1). For example, in the case of AFGL 3068 this seems to imply smaller radius than we adopted to better match the values in the literature.

In the case of IRC+10216, the opacity correction for r_CO = 2 × 10¹⁵ cm (instead of r_CO = 3 × 10¹⁵ cm) would result in a mass-loss of Ṁ ≈ 1.4 × 10⁻⁵ M_⊙ yr₋₁ (not displayed), which would best match the estimate from millimetre observations. However in these circumstances the expected opacities in the lowest J lines would be too large (τ ≫ 1). Nonetheless, the derived T_w and T_h are consistent with kinetic temperature profile models as explained in Sect. 6.1. We find an average value of four OBSIDs of Ṁ ∼ 5 × 10⁻⁶ M_⊙ yr⁻¹, which is lower than the range Ṁ ∼ (1−3) × 10⁻⁵ M_⊙ yr⁻¹ scaled from Teyssier et al. (2006), De Beck et al. (2010, 2012), Guélin et al. (2018). However the value obtained for the warm component agrees with the lower limit of this range (because of the higher M_H2). We obtain T_w ∼ 300 K and Ṁ_w ∼ 1 × 10⁻⁵ M_⊙ yr⁻¹ for a radius of r_CO = 3 × 10¹⁵ cm, which is in good agreement with the results of radiative transfer modeling of the same high J CO lines by Decin et al. (2010). Scaling their Ṁ to the same d and X_CO gives Ṁ ∼ 1.2 × 10⁻⁵ M_⊙ yr⁻¹ with an uncertainty of a factor 2.

For V Hya we obtain Ṁ ∼ 8 × 10⁻⁶ M_⊙ yr⁻¹, which is lower than Ṁ ∼ 3 × 10⁻⁵ M_⊙ yr⁻¹ from Camps (2011) who performed radiative transfer calculations using the same PACS spectrum. In this case, however, the spatio-kinematic structure of the molecular outflow is more complex than assumed in this case. In particular, multiple kinematic (fast and slow) components seem to be present (Hirano et al. 2004; Sahai et al. 2009), which not only translates into a larger uncertainty in the characteristic value of v_exp in the PACS CO-emitting layers, but also implies that the “equivalent” mass-loss rate is particularly questionable.

As for V Hya, for pPNe and yPNe, the assumption of constant-velocity, spherically symmetric mass loss may also not hold. For completeness, the mass-loss rates estimates for these targets are shown Table B.3, but they are subject to larger uncertainties and have to be interpreted with caution. We checked that, even in these cases, our values are in good agreement with previous estimates (making similar simplifying assumptions) in the literature. For example, for the Red Rectangle, we obtain Ṁ ∼ 2 × 10⁻⁷ M_⊙ yr⁻¹, which is within the (scaled) range of Ṁ ∼ (0.2−1.4) × 10⁻⁷ M_⊙ yr⁻¹ reported by De Beck et al. (2010).

In summary, our results are consistent with the literature within the typical uncertainties, but it is hard to tell what is the exact cause of the slight discrepancies from case to case. One obvious reason is that the simple RD method and our assumption of a characteristic value of r_CO (unknown, but crudely constrained from first principles and observations) only provides a rough estimate of the mass-loss rate. We also note, that the bulk of the CO emission under study is produced in the warm inner layers of the CSEs of our targets down to a location where τ ∼ 1. For very optically thick CSEs, there may be an additional amount of gas that is not fully recovered after the moderate opacity correction applied. Another reason for Ṁ discrepancies is the different number of transitions and range of E_u covered by different studies. Non-LTE excitation and radiative transfer models of the CO emission including a wide range of J-transitions is needed to obtain accurate estimates of the mass-loss rates and, in particular, to address Ṁ time modulations.

6.3. Comparison between the sample of O-rich and C-rich stars

In Fig. C.1 we plot together the results presented to those obtained for the sample of O-rich and S-type stars in Paper I.

We find that the range of T_rot is approximately the same among C-rich and O-rich stars, but the O-rich AGBs have less massive CSEs by typically one order of magnitude; these AGBs also have lower expansion velocities. This in turn reflects lower mass-loss rates on average. The scenario is reversed in the groups of PNes with the few O-rich PNes having more warm gas than the carbon counterparts. Most likely, the O-rich AGB stars studied in Paper I are typically low massive stars with very large evolution times, while the O-rich post-AGBs and PNe are very massive objects that have undergone through the hot bottom burning stage.

We stress that even though these trends are indicative of clear differences in the properties of the CSEs of the targets in the THROES catalog, they may not be a general property of C-rich versus non-C-rich targets, since our samples are not necessarily unbiased and they are definitely not statistically significant.

7. Summary

In this paper (Paper II), we use Herschel/PACS FIR spectra of a sample of 15 C-rich evolved stars, including AGBs, post-AGBs, and yPNe, from the THROES catalogue (Ramos-Medina et al. 2018b). These data contain valuable information about the physical-chemical properties of evolved stars as shown, for instance, by the striking differences of spectral features (molecular, atomic/ionized and solid state) as functions of evolutionary stage. In this work, we focus on the rotational spectrum of CO (up to J = 45 − 44), which was used as a proxy for the molecular component of the gas in the warm regions of the CSEs. Our findings can be summarized as follows:

• Because of the higher sensitivity of Herschel compared to ISO, the range of detected ¹²CO transitions has been extended to high rotational levels of up to J_u=45 in low-to-intermediate mass evolved stars. Rotational diagrams using high-excitation CO (v = 0) rotational emission lines, with upper-level energies E_u ∼ 580 to 5000 K, have been plotted to estimate 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 range of temperatures found in our sample, T_rot ∼200− 700 K, is larger than what had been deduced from millimetre/submillimetre observations, and even Herschel/HIFI and SPIRE observations, confirming that PACS CO lines probe deeper layers yet poorly studied to date (typically, ≈10¹⁵ cm for AGBs and ≈10¹⁶ cm for post-AGBs and yPNes).)

• The total gas mass of the warm envelope layers sampled by PACS data are between M ∼ 10⁻⁵ − 10⁻³ M_⊙; post-AGBs and yPNe are overall more massive.

• We find clearly different temperature distributions for the different classes with AGBs having typically hotter gas (up to T_rot  ∼  1000 K) than post-AGBs (T_rot ≲ 500 K) and yPNes (T_rot ≲ 400 K). The yPN AFGL 618 is a clear outlier with a very high amount (M_H2  ∼  2 × 10⁻³ M_⊙) of rather hot (up to T_h  ∼  900 K) gas, similar to the most massive AGBs in the sample.

• For AFGL 3116, CIT 6, AFGL 2513, V Hya, IRAS 16594−4656, and AFGL 618 a double temperature (hot and warm) component is inferred from the RDs. The mean temperatures of the warm and hot components are ∼400 K and ∼820 K, respectively. The mass of the warm component (∼10⁻⁵ − 8 × 10⁻³ M_⊙) is always larger than that of the hot component by a factor ∼4–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 T ∝ r^−0.4, that is, slightly shallower than in the outer envelope layers, in agreement with recent studies.

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

• We investigated the impact of CO line flux variability on the values of M_H2 and T_rot derived from the simple RD analysis. We studied in detail the case of the Mira-variable AGB star IRC+10216, for which multi-epoch PACS data exist. In spite of strong line flux variability we find that the total gas mass and average temperature derived from the RDs at different epochs are minimally affected. Only the hot component shows a sign of line variability (δT/T ∼ 16%), roughly in-phase with the continuum periodicity.

• Similar to Paper I, we find an anticorrelation between T_rot and M_H2, which may result from a combination of CO line cooling and opacity effects, and we find a correlation between Ṁ and v_exp. This correlation is consistent with the wind acceleration mechanism being more efficient the more luminous/massive the star. These trends had been reported in previous studies using low-J CO transitions.

We show that high-J CO emission lines probed by Herschel/PACS are good tracers of the warm gas (T  ∼  200 − 900 K) surrounding evolved carbon stars. Using the simple RD technique, we provided systematic and homogeneous insight into the deepest layers of these CSEs, although it relies on several approximations. Detailed non-LTE excitation and radiative transfer calculations are needed to determine the temperature stratification of the CSEs, to infer mass-loss rates, and to address their time variability.

Acknowledgments

We thank the referee for useful comments and remarks. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). This publication makes use of data products from the THROES catalog, which is a project of the Centro de Astrobiología (CAB-CSIC) with the collaboration of the Spanish Virtual Observatory (SVO), funded by the European Space Agency (ESA). J.M.S.S. acknowledges financial support from the ESAC Faculty and the ESA Education Office under the ESAC trainee program. The Institute for Solar Physics is supported by a grant for research infrastructures of national importance from the Swedish Research Council (registration number 2017-00625). 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.

References

Appendix A: Methods

A.1. Opacity correction

In case of optically thick emission, an optical depth correction factor C_τ should be added to the RD to compensate for an underestimation of the mass within the CO-emitting volume. The optical depth at the line centre (τ) of a given CO (J_u–J_l) transition is

(A.1)

where with v_exp being the expansion velocity of the gas, is the column density of the upper level, A_ij is the Einstein coefficient for spontaneous emission, and λ₀ is the peak wavelength. It follows that the correction factor is written as

(A.2)

To compute C_τ, we perform a first fit of the RD data points, , starting with null opacity correction (lnC_τ = 0), which yields initial (or opacity-uncorrected) values of T_rot and N_CO. We use these values to calculate τ and C_τ and to apply the opacity correction, that is, = . A second fit to is then performed, which renders the so-called opacity-corrected values of T_rot and N_CO (Table B.1).

We want to highlight that, as explained in Goldsmith & Langer (1999), the opacity correction is only reliable for moderate values of the optical depth. For this reason, for all objects in our sample, the minimum acceptable value of r_CO used to compute the CO column density is always chosen so that it results in values of τ close to, but smaller than, unity (Fig. A.1). Indeed, the envelope layers where τ ∼ 1 are the deepest regions observationally accessible because of the almost null escape probability (1/C_τ →0) from deeper, optically thicker regions. As seen in Fig. A.1, for r_CO ≲ 1 × 10¹⁵ cm the opacity of the CO J = 14–13 line (E_u = 580 K), which is the optically thickest transition in our sample, becomes larger than 1 in all our targets. The expected size of the CO-emitting volume in our sample is discussed in detail in Sect. 4.

The optical depth of the line is also sensitive to the expansion velocity of the gas (Eq. (A.1)), which is unconstrained from PACS data and has been assumed to be the terminal expansion velocity of the AGB CSE from the literature (values and references are given in Table B.2). In the case of pAGBs and yPNe, we adopt the average expansion velocity of the bulk of the envelope from different previous works. The uncertainty of v_exp is normally less than 10%.

Fig. A.1. Optical depth at line centre (τ) as a function of the energy of the upper level (Eu). τ was computed for a range of radii; the characteristic radius adopted (rCO, Table B.1) is indicated with a thick line. In each panel, the inset shows the values of Trot and MH2 after the Cτ correction has been applied for the corresponding radii, along with the opacity-uncorrected results (open triangle).

Fig. A.1. continued.

Fig. A.2. Bayesian information criterion and residual sum of squares as functions of number of breakpoints. The BIC is given by the solid line and the RSS is given by the dashed line. The statistic was not computed for IRAS 21282+5050 due to insufficient number of points.

A.2. Model selection with BIC

In this appendix, we describe the method used to automatically find a break point in the linear relation in the RD using the BIC⁷.

Our approach concerns the minimization of the RSS (residual sum of squares) by computing the BIC for a range of possible locations of break points. The BIC can be regarded as a rough estimate of the Bayes factor (Schwarz 1978), and it is given by

(A.3)

where a penalty term, kln(n), penalizes model complexity depending on the number of parameters k and data points n (see also Maindonald & Braun 2010). We also tested the Akaike information criterion (AIC; Akaike 1974), but we concluded that for this particular application it tends to overfit essentially because the penalty term is smaller. In fact, because of the small number of data points, we are not interested in splitting the RD too much in order to obtain robust results. It can happen that BIC suggests more than one breakpoint due to the smaller accuracy of the higher J line fluxes, and in that case only the first break point is taken into account. Obvious line blends were excluded upon the computation of both statistics.

Figure A.2 shows the BIC test in a graphical way where we see that the BIC curve is usually minimal at 0 or 1 breakpoints. For example for AFGL 2513 it is shown that the residuals of the fit keep decreasing with two breakpoints, although this is strongly penalized by the BIC, so a two-component fit is favored. On the contrary, in AFGL 3068, IRAS 15194−5115, and IRC+10216 a single temperature component suffices to reproduce the data.

We note that these methods do not prove that a double temperature component is physically true and they do not provide an alternative explanation for the trends in the data. They simply highlight hidden patterns in the residuals, which can be caused by many effects, of which line blend or heteroscedasticity are the most obvious among them.

Appendix B: Tables

Appendix C: Supplementary figures

Fig. C.1. Supplementary plots. Temperature, mass, mass-loss-rate, line to continuum ratio, and IRAS colour of our sample of C-rich (filled symbols) and O-rich stars in Paper I (open symbols and magenta pentagrams). The colour code is the same as in Fig. 1.

Fig. C.2. Spectra of IRC+10216 at different epochs. Each row corresponds to a different observing day (from top to bottom: OD = 745, 894, 1087, 1113, 1257, 1288, 1296). Same colour code for the top tick marks as in Fig. 2.

All Tables

All Figures

	Fig. 1. IRAS colour diagrams for the stars in the THROES catalogue. The colours are defined from the infrared fluxes at 12, 25, and 60 μm. The boxes in the left panel are those defined in van der Veen & Habing (1988), and the highlighted points correspond to the stars studied in this paper.
In the text

	Fig. 2. Continuum-subtracted PACS spectra of our sample. The vertical segments (at the top) and bars indicate the rotational transitions of 12CO (light grey dotted lines), 13CO (dark grey), HCN (red), CS (yellow), OH (cyan), CH+ (green), and forbidden lines (blue) of [C II] 157.7 μm, [O I] 63.2, 145.5 μm, [N II] 121.9 μm.
In the text

	Fig. 2. continued.
In the text

Fig. 3. Line and continuum variability of IRC+10216 in the FIR. Top panel: continuum-subtracted spectra at different epochs (see also Fig. C.2); the dotted lines represent the CO transitions. Bottom panel: sine-wave fit to the continuum variability at two selected wavelengths for a fixed pulsation period of 630 days (solid lines) and a free-period fit corresponding to 680 and 660 days for the blue and red dashed curves, respectively.

In the text

	Fig. 4. Correlation between CO line flux and continuum intensities. Top panel: CO J = 15 − 14 vs. IRAS 100 μm flux. Middle panel: CO J = 15 − 14 vs. PACS 170 μm continuum flux. Bottom panel: continuum normalised CO J = 15 − 14 flux vs. IRAS [12]–[25] colour. The symbols and colours are the same as in Fig. 1.
In the text

Fig. 5. Rotational diagrams of the CO molecule. The grey line corresponds to a single least-squares fit to the full range of transitions from where a rotational temperature, Trot, and total gas mass, MH2, is computed. The characteristic radius of the CO-emitting volume (rCO) adopted is indicated. The red and blue lines correspond to a two-component model consisting of a “warm” and “hot” region, respectively. Asterisks indicate line blends that were excluded in the fit.

In the text

	Fig. 5. continued.
In the text

	Fig. 6. Rotational diagram of the CO molecule in IRC+10216 at different epochs. Analogous to Fig. 5.
In the text

Fig. 7. Summary of the RD results. Top panel: coloured symbols correspond to the opacity corrected values (single fit), which are connected by dotted lines to the corresponding uncorrected values in grey; on the sides we show the histograms of temperature and mass of the single fits; in the case of IRC+10216 we show weighted averages of four observations. Bottom panel: temperature and mass of the warm and hot components are connected by dashed lines for the same target (if applicable).

In the text

	Fig. 8. Logarithm of the mass-loss rate vs. total CO J = 15 − 14 flux.
In the text

	Fig. 9. Rotational temperature and mass vs. operational day in IRC+10216. Top panel: sinusoidal fit to the Trotvariation with fixed period of 630 days (solid line). Bottom panel: total gas mass over time. In each panel, the dotted lines indicate the average of each component and the asterisks denote unused data points in the fit (see text). The colour code is the same as in Fig. 6.
In the text

Fig. 10. Hot/warm temperature ratio vs. hot/warm mass ratio for AGBs and post-AGBs and the power-law index. Top panel: lines correspond to power-law fits with index α = 0.4 for the opacity corrected RDs (solid) and α = 0.44 for the uncorrected RDs (dashed); the open symbols correspond to the O-rich AGB stars in Paper I. Bottom panel: power-law index for each individual target vs. mass-loss rate (single fit).

In the text

	Fig. 11. Comparison between the mass-loss rate in this work and the literature. The markers correspond to the representative radii listed in Table B.1 for each target and the error bars are mass-loss rates for a given range of radii. The shaded area delimits the range of values found in the literature scaled to the same parameters assumed in this work (see text).
In the text

Fig. A.1. Optical depth at line centre (τ) as a function of the energy of the upper level (Eu). τ was computed for a range of radii; the characteristic radius adopted (rCO, Table B.1) is indicated with a thick line. In each panel, the inset shows the values of Trot and MH2 after the Cτ correction has been applied for the corresponding radii, along with the opacity-uncorrected results (open triangle).

In the text

	Fig. A.1. continued.
In the text

	Fig. A.2. Bayesian information criterion and residual sum of squares as functions of number of breakpoints. The BIC is given by the solid line and the RSS is given by the dashed line. The statistic was not computed for IRAS 21282+5050 due to insufficient number of points.
In the text

	Fig. C.1. Supplementary plots. Temperature, mass, mass-loss-rate, line to continuum ratio, and IRAS colour of our sample of C-rich (filled symbols) and O-rich stars in Paper I (open symbols and magenta pentagrams). The colour code is the same as in Fig. 1.
In the text

	Fig. C.2. Spectra of IRC+10216 at different epochs. Each row corresponds to a different observing day (from top to bottom: OD = 745, 894, 1087, 1113, 1257, 1288, 1296). Same colour code for the top tick marks as in Fig. 2.
In the text