© M. Jeste et al. 2022

1 Introduction

Low- and intermediate-mass stars will reach the asymptotic giant branch (AGB) phase in a late stage of their evolution, during which they will play a crucial role in the cosmic matter cycle (e.g. Kippenhahn et al. 2012; Höfner & Olofsson 2018). Due to the mass-loss processes in this phase, circumstellar envelopes (CSEs) are created around these AGB stars. Depending on the element ratio [C]/[O] in the stellar photosphere, AGB stars can belong to one of three types: oxygen-rich stars (also known as M-type stars, [C]/[O] < 1), S-type stars ([C]/[O] ≈ 1), or carbon-rich stars (also known as C-type stars, [C]/[O] > 1) (Habing 1996; Höfner & Olofsson 2018). In light of different chemical compositions for different AGB types (Tsuji 1964; Cherchneff 2006; Schöier et al. 2013), different tracers should be carefully selected to study the physical conditions of CSEs towards different types of AGB stars.

The HCN molecule is of great interest for studying the CSEs of carbon-rich stars because it is one of the parent molecules that is formed in the stellar atmosphere, and is also the most abundant molecule after H₂ and CO in such stars, as suggested by chemical models (e.g. Tsuji 1964; Cherchneff 2006). Observations have also confirmed that the HCN fractional abundances relative to H₂ are as high as ~3×10⁻⁵ for carbon-rich stars (e.g. Fonfría et al. 2008; Schöier et al. 2013). The HCN molecule is thus widely used to study the CSEs of carbon-rich stars.

The HCN molecule has three vibrational modes, corresponding to the CH stretching mode, ν₁ , the doubly degenerate bending mode, ν₂, and the CN stretching mode, ν₃, together denoted as (ν₁, ν₂, ν₃). The wavelengths of the fundamental transitions of these modes are 3.0, 14.0, and 4.8 µm, which correspond to excitation temperatures of 4764, 1029, and 3017 K, respectively (e.g. Adel & Barker 1934; Barber et al. 2014; Menten et al. 2018). These infrared ro-vibrational lines are difficult to observe with ground-based telescopes due to the Earth’s atmosphere. In contrast, numerous HCN rotational transitions from many vibrational states are easily accessible with ground-based millimetre and sub-millimetre telescopes. Many such HCN transitions have been widely detected in CSEs of carbon-rich stars (Cernicharo et al. 2011; Bieging 2001), the preplanetary nebula CRL 618 (Thorwirth et al. 2003), star-forming regions (Rolffs et al. 2011), and even external galaxies (Sakamoto et al. 2010; Martín et al. 2021). Because of the high extinction and density, the innermost regions of CSEs are difficult to study in lower-J rotational transitions from the vibrational ground state, as they become optically thick; where J is the angular momentum quantum number. The HCN rotational transitions in vibrationally excited states usually have much lower opacities. These properties, together with their high energy levels and critical densities, make them exceptional tracers for investigating the innermost regions of CSEs around carbon-rich stars, which are still poorly constrained. Similar to the observations of vibrationally excited SiO and H₂O lines towards oxygen-rich stars (Wong et al. 2016), observations of vibrationally excited HCN lines may pave the way to understanding the acceleration zone and the atmosphere around the photosphere of carbon-rich stars. Furthermore, detections of multiple transitions of HCN will place constraints on the excitation conditions of the maser lines arising in these regions.

Studies of vibrationally excited HCN transitions have been limited. A large number of HCN J = 3−2 rotational transitions from 28 vibrational states have been investigated towards IRC+10216 (Cernicharo et al. 2011). Using Herschel/PACS and SPIRE data, Nicolaes et al. (2018) performed a population diagram analysis of HCN vibrationally excited lines for eight carbon-rich stars, but they could not resolve the line profiles due to the coarse spectral resolution of their data. Both studies infer high excitation temperatures, supporting the idea that these transitions arise from the innermost region of CSEs. High spectral resolution observations towards more carbon-rich stars can provide more statistical constraints on the role of vibrationally excited HCN transitions for studying the entire class of carbon-rich AGB stars.

Masers are important tools for studying CSEs (e.g. Habing 1996; Reid & Honma 2014). The OH, SiO, and H₂O masers are common in CSEs of oxygen-rich stars, but these masers have not yet been detected towards carbon-rich stars (Gray 2012). Instead, HCN and SiS transitions are commonly found to show maser actions towards carbon-rich stars (e.g. Lucas et al. 1988; Bieging 2001; Henkel et al. 1983; Fonfría Expósito et al. 2006; Schilke et al. 2000; Gong et al. 2017; Fonfra et al. 2018; Menten et al. 2018). Compared with SiS masers, HCN masers are found to be more widespread (Menten et al. 2018). However, previous maser studies have been limited to certain HCN transitions and a low number of targets (see Table A.1). Extending the search for HCN masers to other vibrationally excited HCN transitions and towards more carbon-rich stars may shed light on the pumping mechanism.

With the motivations mentioned above, studies of HCN rotational lines in vibrationally excited states with high spectral resolutions are indispensable in order to investigate the innermost region of carbon-rich CSEs. We therefore undertake new observations of multiple vibrationally excited transitions with different rotational quantum numbers, J, towards a selected sample of carbon-rich stars. In this paper, we start with an overview of our sample of carbon-rich stars in Sect. 1, followed by the observations and data reduction of the stellar sample in Sects. 2 and 3. We present our results in Sect. 4 and discuss all the sources (some individually) in detail in Sect. 5, followed by a brief discussion on the aspect of summary and outlook of the project in Sect. 6.

2 Sample of carbon stars

Table 1 lists the information of the 16 carbon-rich stars in our survey. These sources were selected by applying the following two criteria: (1) the targets should be accessible with the Atacama Pathfinder Experiment¹ telescope (APEX); (2) the carbon-rich nature of these stars and the detection of HCN have already been confirmed by previous studies. We ended up with a 16-object sample, three of which were from Menten et al. (2018), where the authors explored 177 GHz HCN maser emission towards a sample of 13 carbon-rich AGB stars, another three from Rau et al. (2017), and the remaining ten from Massalkhi et al. (2018). These stars have mass-loss rates of the order of 10⁻⁷ to 10⁻⁴ M_⊙ yr⁻¹, spanning about three orders of magnitude. Additionally, we estimated the pulsation phase of the observed stars for different observational epochs (see Table B.1), which is discussed further in Appendix B.

3 Observations and data reduction

We carried out spectral line observations towards the selected 16 carbon-rich stars with the APEX 12-m telescope (Güsten et al. 2006) from 2018 June 8 to 2018 December 20 under the projects M-0101.F-9519A-2018 and M-0102.F-9519B-2018. The total observing time is 185.06 h. We observed a total of 26 rotational transitions of HCN, J = 2−1, J = 3−2, and J = 4−3, in various vibrationally excited states. Table 2 shows the list of the observed HCN lines. We made use of SEPIA (a Swedish-ESO PI receiver; Billade et al. 2012; Immer et al. 2016), PI230, and FLASH (both MPIfR receivers; Klein et al. 2014), as the frontend instruments to observe HCN transitions in J = 2−1, 3−2, and 4−3, respectively. Table 3 gives a brief summary of the instrumental setups of our observations. All HCN lines in the same frequency band were observed and calibrated simultaneously in a single setup. The Fast Fourier Transform Spectrometers (FFTSs) were connected as backend instruments to analyse the spectral data (Klein et al. 2006). These FFTS modules provide us with spectral resolutions of 76 kHz (or 0.13 km s⁻¹), 61 kHz (or 0.07 km s⁻¹), and 38 kHz (or 0.03 km s⁻¹) for lines at 177 GHz, 267 GHz, and 355 GHz, respectively.

Additional observations of IRC +10216, II Lup, and RAFGL 4211 (i.e. the three targets from Menten et al. 2018), were observed again with the SEPIA receiver (η_mb = 0.85) in October 2021 under the project M-0108.F-9515B-2021 to study the variability of the vibrationally excited HCN lines in J = 2−1 (Sect. 4.3).

The spectra were obtained with a secondary mirror, wobbling at a rate of 1.5 Hz and a beam throw of ± 60″ in azimuth. Total integration times for different targets are summarised in Table 1. The chopper-wheel method was used to calibrate the antenna temperature, $T_{\text{A}}^{*}$ (Ulich & Haas 1976), and the calibration was done every 8–10 min. We converted the antenna temperatures to main-beam temperatures, Tmb, with the relationship, $T_{\text{mb}}-T_{\text{A}}^{*}{\eta }_{\text{f}}/{\eta }_{\text{mb}}$, where η_f and η_mb are the telescope’s forward efficiency and main-beam efficiency, respectively. We obtained the values of η_f, η_mb, and the Jansky to Kelvin (Jy K⁻¹) conversion factor for HCN J = 3−2 from the telescope’s report². To get a precise η_mb for the J = 2−1 and J = 4−3 transitions of HCN, which was not included in the telescope’s report, we calculated the antenna temperature from the cross scans of Mars taken during the observations and compared it with the modelled brightness temperature of Mars at that time³. Table 3 lists the main-beam efficiency and Jy K⁻¹ factor for each HCN line. The telescope’s focus was checked after sunrise and sunset, while the pointing, based on nearby pointing sources, was accurate to within ~2″. The full width at half maximum (FWHM) beam sizes, are 36″, 24″, and 18″ at 177 GHz, 267 GHz, and 356 GHz, respectively. Velocities are given with respect to the local standard of rest (LSR). A brief summary of the observed HCN rotational lines in various vibrational states is given in Table 2.

Data reduction was performed with the GILDAS⁴ software package (Pety 2005). The first-order polynomial baseline was subtracted for all observed spectra. The spectra were smoothed at the expense of velocity resolution to improve the signal-to-noise ratios (S/Ns).

4 Results

Figures 1–2, 3, and 4 present the observed spectra towards the three exemplary sources, IRC +10216, RAFGL 4211, and II Lup, respectively, while the spectra of the other 13 sources are presented in Appendix C. In our study, the detection criteria is set to be a S/N of at least three.

4.1 Ground-state HCN transitions

HCN J = 2−1, 3−2, and 4−3 transitions in the ground state are detected towards all 16 sources. These spectral line profiles are parabolic, which is characteristic of a saturated line from an expanding CSE (e.g. Morris 1975; Olofsson et al. 1982). We fit the line profiles in the GILDAS/CLASS software using the shell method to obtain line parameters including velocity-integrated intensity (f T_mbdV), systemic LSR velocity (Vlsr), and expansion velocity (V_exp). The results are given in Table 4. Compared with the integrated HCN J = 2−1 (0, 0, 0) intensities of IRC +10216, II Lup, and RAFGL 4211 in Menten et al. (2018), our values are lower by a factor of ~50% on the T_mb scale. Regular measurements of calibration sources with APEX in the HCN J = 2−1 line find a scatter of intensities of about 15%. Given the large spatial extent of the HCN emission (Dayal & Bieging 1995; Menten et al. 2018), variability of the vibrational ground state HCN lines is unlikely. We therefore attribute this difference to calibration uncertainties⁵. Comparing with the published results (see Table 1), we find that our derived V_LSR are generally higher than the systemic LSR velocities in the literature by 1−3 km s⁻¹, except for V636 Mon. As already explained by previous studies (Olofsson et al. 1982; Menten et al. 2018), this velocity difference is attributed to the presence of blueshifted self-absorption caused by low-excitation gas in the outer parts of CSEs. Our observations demonstrate that such HCN blueshifted self-absorption is ubiquitous in carbon-rich stars. Figure 5 presents the distribution of expansion velocities derived from the HCN J = 3−2 ground-state line towards our sample. The median value of the expansion velocity is 16.1 km s⁻¹, higher than previously reported medians in carbon-rich stars in Ramstedt et al. (2009) because our sample includes objects with higher mass-loss rates than the Ramstedt et al. study.

4.2 Vibrationally excited HCN transitions

HCN rotational transitions in vibrational states with upper energy levels of up to 4789 K (see Table 2) were targeted in our study. We detected HCN rotational lines in the following 11 excited vibrational states: (1, 0, 0), (0, 0, 1), (0, 1^1e, 0), (0, 1^1f, 0), (0, 2^2e, 0), (0, 2^2f, 0), (0, 2⁰, 0), (0, 3^1e, 0), (0, 3^1f, 0), (0, 3^3e, 0), and (0, 3^3f, 0). Table 5 summarises the detection statuses in each source. Vibrationally excited HCN transitions are found to show more complex line profiles, because both thermal and maser emission may contribute to the observed line profiles. The properties of the thermal emission and the maser emission are discussed in Sect. 4.2.1 and Sect. 4.2.2, respectively.

4.2.1 Thermal emission

We detect 39 vibrationally excited HCN thermal lines towards IRC +10216, II Lup, RAFGL 4211, S Aur, W Ori, RV Aqr, IRC +30374, and CRL 3068, while thermal lines are not detected towards the other sources. Most of these sources have high mass-loss rates of ≳1×10⁻⁵ M_⊙ yr⁻¹, except for RV Aqr and W Ori, which have lower mass-loss rates of <3×10⁻⁶ M_⊙ yr⁻¹. These thermal lines display Gaussian-like profiles. In order to derive the observed parameters of the thermal emissions, we assumed a single Gaussian component to fit the observed line profile. The results are shown in Table 6. The peak intensities are within a range of 0.001−0.5 K for these lines. The fitted velocities are generally in a good agreement with the systemic velocities in the literature, within ~3σ (see also Table 1).

All of these lines have FWHM line widths of >8 km s⁻¹. We also note that the velocity spread of thermal vibrationally excited HCN lines can be as broad as that of the ground-state line (see Figs. 1 and C.4 for example). Because we only detected enough thermal vibrationally excited HCN lines in IRC +10216 and CRL 3068, we investigated the fitted line widths as a function of upper energy levels for these two sources. The results are shown in Fig. 6. The line widths do not have an apparent dependence on either their upper energy levels or their Einstein A coefficients (see Table 2). Even though a couple of transitions have similar upper energy levels, their line widths can vary by a factor of 1.5. Particularly, the ν₂ = 3 lines are found to be the narrowest, with expansions of ~5 km s⁻¹ in IRC +10216 due to their high energy levels and critical densities. The line widths are comparable to those of other high energy transitions from parent molecules (e.g. HCN, CS, SiS; Highberger et al. 2000; He et al. 2008; Patel et al. 2009). Comparing with the models of the velocity structure in IRC +10216 (e.g. Fig. 5 in Decin et al. 2015), we find that the ν₂ = 3 lines might arise from the wind acceleration zone, which is about 8 to 10 R_* from the star.

The innermost regions of CSEs have high extinctions, volume densities, and temperatures. Such regions are poorly studied because they are not easily probed with the standard dense-gas tracers in the ground-state lines due to effects such as photon trapping and self-absorption. Our observations extend the spectrally resolved detection of the vibrationally excited thermal HCN transitions towards carbon-rich stars beyond IRC +10216, which paves the way towards a better understanding of the innermost regions of carbon-rich CSEs with vibrationally excited HCN thermal transitions.

4.2.2 Maser emission

We identified HCN maser action following the criteria that the line profiles display either higher peak intensities than expected, or strong asymmetry (i.e. narrow spike features) overlaid on a broader line profile. For instance, different intensities of the J = 2−1 (0, 1^1e, 0) and J = 2−1 (0, 1^1f, 0) lines can be taken as evidence for maser action. Under local thermodynamic equilibrium (LTE) conditions, these two lines are expected to emit at nearly identical intensities, so the intensity difference of these two lines indicates that at least one of them deviates from LTE. In our observations, the J = 2−1 (0, 1^1e, 0) line is significantly brighter than the J = 2−1 (0, 1^1f, 0) line (see Figs. 1−4), which is suggestive of maser actions. Symmetric profiles are generally expected from thermal emission of typical AGB CSEs and narrow asymmetric spikes are not expected in thermal line profiles by the classic models (e.g. Morris 1975; Olofsson et al. 1982). Previous observations suggested that such asymmetric spikes arise from maser action (e.g. Henkel et al. 1983; Schilke et al. 2000; Fonfría Expósito et al. 2006; Gong et al. 2017; Fonfría et al. 2018), because these masers are only formed at specific regions in CSEs, which correspond to specific velocities in their line profiles. Hence, the presence of narrow asymmetric spikes are regarded as an indicator of maser action in this work.

In this study, we detect 29 HCN maser lines from 12 carbon-rich stars, whereas HCN masers are not detected towards Y Pav, IRC +30374, CRL 2513, or CRL 2477 in our observations. Except for six maser lines already reported by previous observations (e.g. Menten et al. 2018), 23 transitions are newly discovered HCN masers. Hence, our observations expand the total number of previous known HCN masers towards carbon-rich stars by 47% (see Table A.1). Our results also confirm that HCN masers are common in carbon-rich stars. In particular, eight CSEs, namely CRL 3068, RV Aqr, R Vol, AQ Sgr, V Aql, V636 Mon, S Aur, and IRC +20370, are newly discovered hosts of HCN masers. The 29 HCN masers arise from eight different transitions, namely J = 2−1 (0, 1^1e, 0), J = 2−1 (0, 1^1f, 0), J = 3−2 (0, 1^1e, 0), J = 3−2 (0, 3^1e, 0), J = 3−2 (0, 2, 0), J = 4−3 (0, 1^1e, 0), and J = 4−3 (0, 1^1f, 0). Among them, the J = 3−2 (0, 3^1e, 0), J = 3−2 (0, 2, 0), J = 4−3 (0, 1^1e, 0), and J = 4−3 (0,1^1f, 0) masers are first detected in the interstellar medium.

IRC +10216. Five HCN masers, including J = 2−1 (0, 1^1e, 0), J = 2−1 (0, 1^1f, 0), J = 3−2 (0, 1^1e, 0), J = 3−2 (0, 1^1f, 0), and J = 4−3 (0, 1^1f, 0), are detected in IRC +10216. The J = 2−1 (0, 1^1e, 0) maser was reported in previous observations (Lucas & Cernicharo 1989; Menten et al. 2018). Our 2018 observations show at least three peaks at −32.2, −27.9, and −26.3 km s⁻¹, which have peak flux densities of about 145, 295, and 351 Jy, respectively. The intensity of this line is nearly 20 times that of the J = 2−1 (0, 1^1f, 0) line (see Fig. 1). Comparison with other measurements at different epochs are presented in Sect. 4.3. Furthermore, we find that the flux density is lower for the components with LSR velocities farther from the systemic velocity of −26.5 km s⁻¹. The J = 2−1 (0, 1^1f, 0) line was believed to show maser actions due to its intensity variation (Menten et al. 2018). This maser peaks close to the systemic velocity and has a flux density of about 17 Jy. The J = 3−2 (0, 1^1e, 0) line shows as a narrow spike overlaid on a broad thermal profile, while such a feature is not seen in the J = 3−2 (0, 1^1f, 0) line. The J = 3−2 (0, 1^1e, 0) maser peaks at a redshifted LSR velocity of −24.0 km s⁻¹. Subtracting the broad component, we obtained a flux density of 17.5 Jy. We note that this line was observed with NRAO-12 m and IRAM-30 m (Ziurys & Turner 1986; Cernicharo et al. 2011), but they did not detect a spike feature. On the other hand, the J = 3−2 (0, 1^1e, 0) maser was detected by previous extended Submillimeter Array (eSMA) observations (Shinnaga et al. 2009). We also see a presence of an asymmetrical profile feature in the J = 3−2 (0, 1^1f, 0) line, which is marked as a maser. He et al. (2017) monitored the line shape of this maser and reported such a variability in its blue-shifted peak. The J = 4−3 (0, 1^1f, 0) line appears to show a prominent asymmetry narrow spike on the blushifted side, peaking at an LSR velocity of −35.9 km s⁻¹, while no asymmetric features are found in the J = 4−3 (0, 1^1e, 0) line. The flux density is estimated to be 68 Jy by subtracting the broad component. The flux densities of our detected HCN masers are at least a factor of two lower than those (>800 Jy) of sub-millimetre HCN lasers (Schilke et al. 2000; Schilke & Menten 2003). Furthermore, we find that the intensity of the J = 3−2 (0, 3^1e, 0) line may also deviate from the values expected under LTE, because its peak intensity is about 1.4 times that of the J = 3−2 (0, 3^1f, 0) line. However, the tentative maser nature of the J = 3−2 (0, 3^1e, 0) line still needs to be confirmed by follow-up observations because the signal-to-noise ratios of the two lines are not good enough for the difference to be significant.

RAFGL 4211. We detected four HCN masers. Among them, the J = 2−1 (0, 1^1e, 0) maser was first reported by Menten et al. (2018). Our J = 2−1 (0, 1^1e, 0) spectrum shows at least two peaks at −8.1 km s⁻¹ and −5 km s⁻¹ of flux density of 24.2 Jy and 32.5 Jy, respectively. These peak flux densities are about 20 times that of the J = 2−1 (0, 1^1f, 0) line (see Fig. 3), strongly supporting its maser nature. The flux density of the J = 2−1 (0, 1^1e, 0) maser is comparable to that of the J = 1−0 (0, 2, 0) maser in this source (see Fig. 1 in Smith et al. 2014). Although the J = 3−2 (0, 2, 0) and J = 3−2 (0, 3^1e, 0) masers only have flux densities of 2.4 Jy and 11.8 Jy, respectively, the detection of the two masers are exceptional, because these two masers have only been discovered towards RAFGL 4211 so far. Both of the masers peak at −2.1 km s⁻¹, close to the stellar systemic velocity. Unlike the other three masers, the J = 4−3 (0, 1^1e, 0) profile exhibits a narrow spike overlaid on a broad component, while such a narrow spike does not appear in the J = 3−2 (0, 1^1e, 0) spectrum. The J = 4−3 (0, 1^1e, 0) maser peaks at a slightly blueshifted side of −3.7 km s⁻¹, and its flux density is about 18.5 Jy, which is fainter than the J = 2−1 (0, 1^1e, 0) maser.

II Lup. We detected the J = 2−1 (0, 1^1e, 0) and J = 4−3 (0, 1^1f, 0) masers in II Lup. The latter one is newly detected in this source. The J = 2−1 (0, 1^1e, 0) maser shows at least three peaks at −22.7, −18.1, and −15.2 km s⁻¹, which have flux densities of 21.3, 39.3, and 39.3 Jy, respectively. This maser is much brighter than the J = 2−1 (0, 1^1f, 0) line (see Fig. 4). Two peaks at −25.4 and −6.2 km s⁻¹ are also identified in the J = 4−3 (0, 1^1f, 0) spectrum. Their corresponding flux densities are 14.3 Jy and 8.2 Jy, which are about four times the peak flux densities at the respective velocity positions of the J = 4−3 (0, 1^1e, 0) spectrum.

R Vol. No HCN masers were reported in this star prior to our observations. In R Vol, our observations result in the first detection of four HCN masers that are the J = 2−1 (0, 1^1e, 0), J = 2−1 (0, 1^1f, 0), J = 3−2 (0, 1^1e, 0), and J = 4−3 (0, 1^1e, 0) lines. The J = 2−1, J = 3−2, and J = 4−3 lines in the (0, 1^1e, 0) state have peak flux densities of 8.4 Jy, 3.8 Jy, and 7.5 Jy at −12.9, −14.5, −11 km s⁻¹, respectively, which are much brighter than their corresponding transitions in the (0, 1^1f, 0) state (see Fig. C.1). This is suggestive of their maser nature. Although the J = 2−1 (0, 1^1f, 0) line only has a flux density of 1.9 Jy, this line is still regarded as a maser in this work, because its velocity coverage is much narrower than the ground state-line and the LSR velocity of its peak is slightly shifted to the blueshifted side at−13.2 km s⁻¹.

AQ Sgr. The J = 3−2 (0, 1^1e, 0) line has peak flux density of 1.2 Jy at 22.6 km s⁻¹. Despite its low flux density, this line is brighter than the J = 3−2 (0, 1^1f, 0) line, which is not detected by the same sensitivity. Furthermore, the peak velocity is slightly redshifted with respect to the systemic velocity (~20 km s⁻¹, see Fig. C.3). Therefore, the J = 3−2 (0, 1^1e, 0) line is regarded as a maser in this study, which makes it the first HCN maser towards this source.

CRL 3068. Two HCN masers are newly discovered in the source: the J = 2−1 (0, 1^1e, 0) and J = 4−3 (0, 1^1f, 0) masers. Asymmetric features indicate that these two transitions show maser actions (see Fig. C.4). The J =2−1 (0, 1^1e, 0) maser displays five distinct peaks at −38.5 km s⁻¹, -35.3 km s⁻¹, −33.8 km s⁻¹, −29.4 km s⁻¹, and −26.8 km s⁻¹, and their corresponding peak flux densities are 12.9 Jy, 23.8 Jy, 18.2 Jy, 9.9 Jy, and 8.6 Jy, respectively. Stronger peaks are observed on the blueshifted side. They are at least four times higher than the peak flux density of the J = 2−1 (0, 1^1f, 0) line. The J = 4−3 (0, 1^1f, 0) maser shows a prominent peak of flux density 70.6 Jy at −41.7 km s⁻¹, which is brighter than any other transitions in the (0, 1^1e, 0) state, including even the J = 2−1 (0, 1^1e, 0) maser. The J = 4−3 (0, 1^1f, 0) maser is about 10 km s⁻¹ away from the systemic velocity of −31.5 km s⁻¹, but still does not reach its terminal velocity of 14.5 km s⁻¹.

V636 Mon. The J = 2−1 (0, 1^1e, 0) line is the only vibrationally excited HCN maser detected in V636 Mon. Its spectrum exhibits a narrow spike at 5.7 km s⁻¹, and its peak flux density is 2.2 Jy (see Fig. C.6). This line is at least three times the peak flux density of the J = 2−1 (0, 1^1f, 0) line. Furthermore, the narrow spike is about 4 km s⁻¹ away from the systemic velocity. These two features support that the J = 2−1 (0, 1^1e, 0) line is a maser, which make it the first maser to be detected in V636 Mon.

IRC +20370. Two HCN masers from the J = 2−1 (0, 1^1e, 0) and J = 3−2 (0, 1^1e, 0) transitions are detected towards this source. Our results are the first maser detection towards this source. Both maser transitions are brighter than their corresponding lines in the (0, 1^1f, 0) state by a factor of roughly ten (see Fig. C.7). The J = 2−1 (0, 1^1e, 0) maser hosts two peaks at −4.6 km s⁻¹ and −2 km s⁻¹, and their corresponding flux densities are 6.7 and 7.3 Jy. The J = 3−2 (0, 1^1e, 0) maser exhibits at least one peak at −4 km s⁻¹ with a flux density of 7.1 Jy, which is slightly weaker than the J = 2−1 (0, 1^1e, 0) maser. These velocity components are slightly blueshifted from the systemic velocity of −0.8 km s⁻¹. The peak flux density at −4.6 km s⁻¹is lower than that at −1.7 km s⁻¹ in the J = 2−1 (0, 1^1e, 0) maser, which is in stark contrast to the J = 3−2 (0, 1^1e, 0) maser. The different velocities indicate that these masers could arise from different zones or different pumping mechanisms.

W Ori. In this study, we report three new HCN masers from the J = 2−1 (0, 1^1e, 0), J = 3−2 (0, 1^1e, 0), and J = 4−3 (0, 1^1e, 0) transitions that display single-peak profiles. Their maser nature is strongly supported by the fact that these three transitions have peak intensities about 20 times those of their corresponding (0, 1^1f, 0) transitions (see Fig. C.8). The respective LSR velocities of the three masers are around 0.13 km s⁻¹, −1.1km s⁻¹, and −2.5 km s⁻¹, while the corresponding flux densities are 14.5 Jy, 16 Jy, and 6.5 Jy. All of these maser peaks are within 3 km s⁻¹ away from the systemic velocity of −1.0 km s⁻¹. W Ori has been found to host the HCN J = 1−0 (0, 0, 0) maser by spike emission at an LSR velocity of ~8 km s−1 and shows temporal variations in this maser (Olofsson et al. 1993; Izumiura et al. 1995). The LSR velocities indicate that our newly detected masers arise from the regions closer to the star when compared with the ground-state J = 1−0 maser. Furthermore, the J = 2−1 (0, 1^1e, 0) maser is more luminous than the HCN J = 1−0 (0, 0, 0) maser.

S Aur. In this source, the J = 2−1 (0, 1^1e, 0) transition exhibits a narrow component with a high flux density of 9.1 Jy at −17.0 km s⁻¹ (see Fig. C.9). Its peak intensity is at least 30 times that of the J = 2−1 (0, 1^1f, 0) transition, and is even twice the peak intensity of the J = 2−1 (0, 0, 0) transitions. Furthermore, its line width is much narrower than that of the J = 3−2 (0, 1^1e, 0) transition. These facts support that the J = 2 − 1 (0, 1^1e, 0) transition is a maser, which is the first maser reported towards this source. We also tentatively detect the J = 3−2 (0, 0, 1) line exhibiting a peak intensity of 8.5 Jy with a very narrow line width of 0.19 km s⁻¹. The emission feature peaks close to the systematic velocity of −17 km s⁻¹. Since this emission feature is only detected in two native spectral channels above the noise level, we conservatively do not identify it as a maser.

V Aql. The J = 3−2 (0, 1^1e, 0) transition is the only vibrationally excited HCN line detected in V Aql. This transition displays a single-peaked spectrum with a peak flux density of ~4.9 Jy at 55.4 km s⁻¹ (see Fig. C.10). The peak flux density is at least three times that of the J = 3−2 (0, 1^1f, 0) transition. Furthermore, the LSR velocity appears to be slightly redshifted with respect to the systemic velocity of 53.5 km s⁻¹. These facts support that the J = 3−2 (0, 1^1e, 0) transition is a maser, which is the first maser detected in V Aql.

RV Aqr. In this study, we report three new HCN masers from the J = 2−1 (0, 1^1e, 0), (0, 1^1f, 0), and J = 3−2 (0, 1^1e, 0) transitions. Figure C.13 shows that the J = 2−1 (0, 1^1e, 0) and J = 3−2 (0, 1^1e, 0) transitions are at least ten times brighter than their corresponding (0, 1^1f, 0) transitions. Compared with the vibrationally ground-state lines, the J = 2−1 (0, 1^1f, 0) line, and also the two detected (0, 1^1e, 0) lines, have much narrower velocity spreads. Furthermore, their peak velocities are about −1 km s⁻¹, which are blueshifted with respect to the systemic velocity of 0.5 km s⁻¹. Therefore, all three transitions are identified as masers in this work, which are the first maser detection in RV Aqr. Their flux densities are 9.5 Jy, 1.9 Jy, and 14.7 Jy for the J = 2−1 (0, 1^1e, 0), (0, 1^1f, 0), and J = 3−2 (0, 1^1e, 0) masers, respectively. On the other hand, the J = 4−3 (0, 1^1e, 0) transition is only slightly brighter than the J = 4−3 (0, 1^1f, 0) transition; the latter has a S/N of approximately two at the peak channel and does not formally meet our detection criteria. We conservatively do not identify the J = 4−3 (0, 1^1e, 0) transition as a maser in this work.

Non-detection. HCN masers are not detected towards Y Pav, IRC +30374, CRL 2513, or CRL 2477 in our observations. Hence, we are only able to give a 3σ upper limit of 0.9 Jy, 1.8 Jy, 0.8 Jy, and 2.6 Jy for the J = 2−1 (0, 1^1e, 0) transitions towards Y Pav, IRC +30374, CRL 2513, and CRL 2477, respectively. Among these sources, we note that the HCN J = 1−0 (0, 2, 0) maser was detected towards CRL 2513 by Lucas et al. (1988), but the higher rotational transitions (J = 2−1 and J = 3−2) in the (0, 2, 0) state are not detected in this work (see Fig. C.11).

In summary, the J = 2−1 (0, 1^1e, 0) maser is detected towards 12 out of 16 sources. The detection rate of 75% is similar to that of the previous J = 2−1 (0, 1^1e, 0) maser study (85%; Menten et al. 2018) and is much higher than that of the J = 1−0 (0, 2, 0) maser in carbon-rich stars (7%; Lucas et al. 1988). This suggests that the J = 2−1 (0, 1^1e, 0) maser is more common than the J = 1−0 (0, 2, 0) maser in carbon-rich stars.

The peak flux density of the J = 2−1 (0, 1^1e, 0) maser is brighter than that of the J = 2−1 (0, 0, 0) line in W Ori and S Aur, while the peak flux densities of these two lines are comparable in RAFGL 4211, R Vol, and RV Aqr. In contrast, the J = 2−1 (0, 1^1f, 0) masers are much weaker and are only detected in three sources. The J = 3−2 (0, 1^1e, 0) maser is detected towards seven sources. The flux density of the J = 3−2 (0, 1^1e, 0) maser is much lower than that of the J = 2−1 (0, 1^1e, 0) maser towards IRC +10216 and R Vol, but is roughly comparable to that of the J = 2−1 (0, 1^1e, 0) maser towards IRC +20370, W Ori, and RV Aqr. The J = 4−3 (0, 1^1f, 0) masers are only detected in IRC +10216, CRL 3068, and II Lup, and all are found to show a prominent blueshifted component with a peak at ~10 km s⁻¹ from their systemic LSR velocities.

Moreover, we find that none of the detected HCN maser peaks reach the terminal velocities, suggesting that these HCN masers arise from shells that are not yet fully accelerated. This supports the notion that these masers can be used to investigate the wind acceleration of carbon-rich stars. Furthermore, most of the bright HCN maser components appear on the blueshifted side with respect to their respective systemic LSR velocities, which may be attributed to their underlying pumping mechanisms (see Sect. 5.2).

Fig. 1 From left to right: Spectra for J = 2−1 (first) and 3−2 (second and third) transitions of HCN towards IRC +10216. The vibrational state of each transition is indicated in the top right corner of the respective spectrum. The vertical dashed blue lines denote the range of the emission in VLSR, the vertical dotted line shows the systemic VLSR of the star, and the horizontal dotted line is the baseline. The number below the vibrational quantum numbers in the panel is a code assigned to the date of observation of the particular transition. The corresponding dates of the codes are listed in the notes of Table 4.

Fig. 2 Spectra for J = 4−3 transitions of HCN towards IRC +10216. The description of the figure is the same as in the Fig. 1 caption.

Fig. 3 Spectra for J = 2−1 (left), 3−2 (centre) and 4−3 (right) transitions of HCN towards RAFGL 4211. The description of the figure is the same as in the Fig. 1 caption.

Fig. 4 Spectra for J = 2−1 (left), 3−2 (centre) and 4−3 (right) transitions of HCN towards II Lup. The description of the figure is the same as in the Fig. 1 caption.

Fig. 5 Histogram distribution of expansion velocities derived from J = 3−2 ground-state HCN lines. The magenta line marks the median for HCN J = 3−2 from this work. The median values calculated by Ramstedt et al. (2009) towards M-type, S-type, and C-type stars are marked in orange, green, and red, respectively.

Fig. 6 FWHM line widths as a function of upper energy levels for IRC +10216 (left) and CRL 3068 (right). The colours represent the rotational transitions and the symbols represent the various thermal vibrationally excited lines, as indicated in the legends.

4.3 Time variability of HCN masers

Three of our targets, IRC +10216, II Lup, and RAFGL 4211, have been observed in the J = 2−1 (0, 1^1e, 0) and J = 2− 1 (0, 1^1f, 0) lines at least twice in previous works (Lucas & Cernicharo 1989; Menten et al. 2018), allowing us to study the time variability of the HCN masers. Figure 7 presents the HCN masers obtained at different epochs on a flux density scale. It can be seen that the line shape of the J = 2−1 (0, 1^1e, 0) line has changed significantly over time. For IRC +10216, the spectra of the (0, 1^1e, 0) line show four peaks at −32.6 km s⁻¹, −29.8 km s⁻¹, −25.1 km s⁻¹, and −22.6 km s⁻¹ in 1989, two peaks at −33 km s⁻¹, and −27.2 km s⁻¹ in 2015, three peaks at about −32.2 km s⁻¹, −26.3 km s⁻¹, and −27.9 km s⁻¹ in 2018, and five peaks at −33.5 km s⁻¹, −32.1 km s⁻¹, −29.7 km s⁻¹, −26.7 km s⁻¹, and −24.1 km s⁻¹ in 2021 (see the top left panel in Fig. 7).

Towards II Lup, this line shows two peaks at −23.6 km s⁻¹ and −15.9 km s⁻¹ in 2015, three peaks at −22.6 km s⁻¹, −18.3 km s⁻¹, and −15.1 km s⁻¹ in 2018, and two peaks again in 2021 at −25.1 km s⁻¹ and −15.3 km s⁻¹ (see the bottom left panel in Fig. 7). Similarly, RAFGL 4211 shows three peaks at −8.2 km s⁻¹, −5.1 km s⁻¹, and −3.7 km s⁻¹ in 2018, two peaks at −5.1 km s⁻¹ and −2.9 km s⁻¹ in 2015, and three peaks in 2021 at −9.7 km s⁻¹, −3.7 km s⁻¹, and −2.6 km s⁻¹ (see the bottom right panel in Fig. 7).

Based on the measurements of the ground-state HCN lines, we note that there might be about a 50% uncertainty in flux calibration. In order to minimise the effect of calibration uncertainty, we used the integrated intensity ratio between the maser line and the ground-state J = 2−1 (0, 0, 0) line to assess the temporal variation between the three epochs of APEX observations (2015−2021). For IRC +10216, the ratio changes from 0.16 → 0.15 → 0.28, suggesting that the flux density significantly increased between 2018 and 2021 (≈86%). For II Lup, the ratio changes from 0.28 → 0.35 → 0.20, indicative of a 25% increase in the flux density of the (0, 1^1e, 0) maser between 2015 and 2018, but a 42% decrease between 2018 and 2021. For RAFGL 4211, the ratio increases from 0.16 → 0.25 → 0.36, indicating that the flux density of the (0, 1^1e, 0) maser increased by 63% in the first two epochs and by 44% in the last two. These suggest that the HCN J = 2−1 (0, 1^1e, 0) masers are variable in the three sources over a timescale of a few years. In contrast, as shown in the top right panel in Fig. 7, the weaker HCN J = 2−1 (0, 1^1f, 0) maser is not quite as variable over the 3 yr in IRC +10216. The spectral profiles are comparable between the three epochs and the intensity ratio between the (0, 1^1f, 0) and (0, 0, 0) lines only changes by ≲10%. However, it is worth noting that the most significant variation in the J = 2−1 (0, 1^1e, 0) transition happened over a longer period in time, from 1989 to 2015 (see Sect. 6.3.4 in Menten et al. 2018).

Although the time variability of HCN masers is confirmed by our study, observations are too sparse to allow for a search for periodic variations or for association with stellar pulsation periods. Previous monitoring observations of the HCN J = 3−2 (0, 1^1f, 0) maser in IRC +10216 indicate a period of ~730 days (He et al. 2017), which is at least 100−200 days longer than the periods in the infrared (~630 ± 3 days; Le Bertre 1992; Menten et al. 2012) and centimetre bands (535 ± 50 days; Menten et al. 2006). An anti-correlation between HCN maser luminosity and infrared luminosity was found in their single-dish observations (He et al. 2017), but was not confirmed by the follow-up monitoring observations with the Atacama Compact Array (He et al. 2019). Previous observations towards RAFGL 4211 indicate that the variability of the HCN J = 1−0 (0, 2, 0) maser may be unrelated to stellar pulsation, although the observations consisted of only three epochs (Smith et al. 2014). Menten et al. (2018) suggest that the apparent absence of correlation between the vibrationally excited HCN lines emission cycle and the stellar cycle may arise from the additional contribution of collisional pumping. Further long-term monitoring observations with higher cadences and consistent flux calibration are needed to identify the exact cause(s) of the time variability of HCN masers, such as inhomogeneity in the outflowing wind or a highly anisotropic infrared field.

Fig. 7 Variability of HCN (2−1) masers seen towards IRC+10216 (top two panels) in (0, lle, 0) and (0, llf, 0), II Lup (bottom left) and RAFGL 4211 (bottom right) in (0, lle, 0) vibrationally excited lines. The black spectra are from observations in 2018, cyan spectra are from 2021, red spectra are from Menten et al. (2018), and blue spectra are from Lucas & Cernicharo (1989) with the date of observations indicated in the legend.

Fig. 8 Rotation diagram of IRC+10216 constructed from vibrationally excited HCN lines. Some masers are also plotted in the diagram, shown by triangles, but are excluded from the linear least-square fit.

5 Discussion

5.1 Population diagram

The vibrationally excited HCN lines detected in our observations occupy energy levels above 1000 Κ (see Table 2). Their spectra, including those of the thermal lines, cover narrow velocity ranges (see Table 6), suggesting that they most likely arise from a hot region where the stellar outflow has not reached its terminal velocity (i.e. close to the stellar surface; Menten et al. 2018). In order to better constrain the properties of the region traced by these transitions, we take advantage of the population-diagram method to estimate the HCN excitation temperature and column density. Assuming LTE conditions and low opacities for the thermal emission of vibrationally excited HCN, we can use the standard formula to derive HCN column density and excitation temperature (e.g. Goldsmith & Langer 1999):

$$N_{\text{u}}/g_{\text{u}}=\frac{N_{\text{tot}}}{Q\left(T_e\right)}{e}^{-E_{\text{u}}/T_{\text{e}}}=\frac{3k\int T_{\text{mb}}\text{d}V}{8{\pi }^{3}v{\mu }^2{S}_{\text{ul}}},$$(1)

where N_tot is the HCN column density, Q(T_e) is the partition function at the excitation temperature T_e, ∫ T_mbdV is the velocity-integrated intensity, k is the Boltzmann constant, ν is the transition’s rest frequency, μ is the permanent dipole moment of HCN, and S_ul is the line strength. We obtained the values of Q, v, and S_ulµ² from CDMS (Müller et al. 2005) via Splatalogue⁶ (see also Table 2). Since the emitting size of vibrationally excited HCN lines should be small, we should correct the observed integrated intensities by dividing by the beam dilution factor, $\frac{{\theta }_{\text{s}}^2}{{\theta }_{\text{s}}^2+{\theta }_{\text{b}}^2}$, where θ_s and θ_b are the source size and FWHM beam size, respectively. However, the source size is still not well constrained. The eSMA observations suggest that the vibrationally excited HCN J = 3−2 (0, l^le, 0) emission can be as extended as 1″ in IRC +10216 (see Fig. 7 in Shinnaga et al. 2009) but with double peaks. Our assumption was to adopt a Gaussian FWHM source size of θ_s = 0.2″ (~28 au) towards IRC +10216. Assuming the physical size to be the same as that of IRC +10216, we can estimate the source sizes for the other sources according to their distances in Table 1. The assumed angular source sizes are given in Table 7.

In order to fit the population diagrams, at least two transitions covering a large enough energy range are needed to meaningfully derive the excitation temperature and column density. Because a large number of vibrationally excited HCN transitions are found to show maser action for which LTE conditions do not hold (see Sect. 4.2.2), linear least-square fits can only be applied to two sources: IRC +10216 and CRL 3068. In the case of IRC +10216, we included lines that are blended with each other: (0, 3^3e, 0) and (0, 3^3f , 0). Assuming LTE conditions, these two lines are expected to show nearly identical integrated intensities. We therefore used half the intensities integrated over the velocity range covered by the two lines for our analysis. Excluding the maser lines and other blended lines, we performed the linear least-square fits to the two sources. The results are shown in Figs. 8–9. In each of these plots, we display, without fitting, the data points of two masers, J = 2–1 (0, 1^1e, 0) and J = 4–3 (0, 1^1f, 0). The J = 2–1 (0, 1^1e, 0) maser lies above the fitted lines in both plots, indicating the highly non-thermal nature of its emission. The J = 4–3 (0, 11f , 0) maser is weaker than the J = 2−1 (0, l^le, 0) maser because its maser component only contributes a small fraction to the total integrated intensity of the transition (see Fig. 1 and Sect. 4.2.2 for more information). The fitted excitation temperatures and HCN column densities from the population diagrams are given in Table 7.

In the range of θ_s=0.1″−1″, we note that the derived HCN column densities are largely dependent on the assumed source size with a relationship of $N_{\text{tot}}\propto 1/{\theta }_{\text{s}}^2$, that is, the derived HCN column densities can vary by at most a factor of 100. On the other hand, the excitation temperature is not sensitive to the adopted value of θ_s. The excitation temperatures of the two sources are consistent, within uncertainties, with the T_ex range of 452−753 Κ as derived from Herschel observations of eight other carbon-rich stars (Nicolaes et al. 2018) and the T_ex values of IRC +10216 (T_ex=1000 K, Schilke et al. 2000; T_ex=410−2465 K, Cernicharo et al. 2011). We also note that Cernicharo et al. (2011) defined three different zones in their population diagram of IRC +10216, based on the upper energy levels of HCN lines. The T_ex value in IRC +10216 derived by our single-component fitting is higher than that of Zone III (410 ± 60 K; E_u≤2000 K), but slightly lower than that of Zone II (1240 ±210 K; 2000<E_u≤5000 K), probably because we fit HCN lines from both zones. Some of our identified thermal lines might in fact show (weak) maser actions because their column densities tend to deviate from other transitions of similar energies in the population diagrams (Figs. 8−9). However, they are not easily identified from their line profiles alone. This may lead to an overestimation of the HCN excitation temperature. Comparing the kinetic temperature model in IRC +10216 (e.g. Fig. 1 in Agúndez et al. 2012), we infer that the vibrationally excited HCN line-emitting regions have a size of <20R_* (O″4), where the stellar radius, R_*, is 22 mas (Monnier et al. 2000). Our assumed source size of O.″2 for IRC +10216 is consistent with this inferred size.

For the other four stars with vibrationally excited HCN lines but without a broad coverage of line excitation energy, we derived their HCN column densities by fixing the excitation temperature to be 700 K, which is a representative value among the results from Herschel data (452−753 K; Nicolaes et al. 2018) and our fitting. The associated uncertainties in column density were derived by Monte Carlo analysis with 10000 simulations. A Gaussian function was used to fit the resulted distribution, and the standard deviation of the Gaussian distribution was regarded as the 1σ error. The results are shown in Table 7. The derived source-averaged HCN column densities range from 0.18−10.8 ×10¹⁹ cm⁻². Among these sources, W Ori has the lowest source-averaged column density, likely due to its lowest mass-loss rate of ~3.l×l0⁻⁷. We note that the uncertainties in the adopted source sizes can lead to even larger uncertainties than those shown in Table 7.

Based on the model of IRC +10216 (e.g. Fig. 1 in Agúndez et al. 2012), the H₂ number density is about 2×10⁸ cm⁻³ at the radius of 0.2″ (i.e. 4×10¹⁴ cm), corresponding to a H₂ column density of 1.6×10²³ cm⁻². The HCN abundance with respect to H₂ is thus determined to be (7.4 ± 5.5) × 10⁻⁵. The HCN abundance is roughly consistent with previous estimates of HCN abundances in carbon-rich stars (Fonfría et al. 2008; Schöier et al. 2013), but it seems to be much higher than the prediction (< 3 × 10⁻⁶ at r >5R_*) in the non-equilibrium chemical model (Cherchneff 2006). However, we also note that the derived HCN abundance largely relies on the assumptions, such as the emission source size, which can result in large uncertainties in the estimated HCN column density and, therefore, its abundance.

Fig. 9 Rotation diagram of CRL 3068. The description of the figure is the same as in the caption of Fig. 8.

Fig. 10 Energy diagram of HCN. The quantum numbers, υ and J, indicate the vibrational state and rotational state, respectively, l-type doubling causes the splitting of every rotational level into two sub-levels, le and 1f, which are also denoted in this figure. Pumping photons from the ground state to highly vibrationally excited states are indicated by cyan arrows. The J = 2−1 (0, lle, 0) and J = 2−1 (0, llf, 0) transitions are indicated by the blue and purple arrows, respectively. The direct l-type transitions are marked by red arrows. We note that the separations of the l-type doublets are not drawn to scale.

5.2 Pumping mechanism of HCN masers

Vibrationally excited HCN masers are generally thought to be explained by pumping through infrared radiations (e.g. Lucas & Cernicharo 1989; Menten et al. 2018). Our detected HCN masers are found in the excited states of the bending mode: v₂ = 1^le or 1^1f, v₂ = 2, and v₂ = 3^1e (Sect. 4.2.2). As illustrated in the energy level diagram in Fig. 10, pumping HCN molecules from the ground state to these excited states requires infrared photons at 14 µm, 7 µm, and 5 µm, respectively (see Fig. 10). Following Menten et al. (2018), we extended the comparison between the isotropic maser luminosities, L_M, and the corresponding infrared luminosities, L_IR, towards all detected HCN masers, as shown in Table 8. Their luminosity ratios, L_IR/L_M, should be more reliable than their luminosities because the ratios are independent of the assumed distances. We find that the luminosity ratios are much greater than 50 for all the detected HCN masers. The infrared flux densities are based on the WISE and IRAS measurements at 4.6 µm and 12 µm (Abrahamyan et al. 2015), but not at the exact wavelengths of HCN vibrational excitation of 14 µm, 7 µm, and 5 µm. While this adds to the uncertainties in the derived luminosity ratios L_IR/L_M, Menten et al. (2018) suggest that the flux densities at 7 µm and 14 µm are within a factor of two of the 12-µm flux density based on the SED towards IRC +10216. Hence, our conclusion is still robust because the luminosity ratios are much higher than two.

The fact that the amount of IR pumping photons is much more than adequate to account for HCN masers indicates that all these masers are unsaturated. The unsaturated nature of detected HCN masers is also supported by their brightness temperature estimates. Assuming a typical maser emission size of 28 au (~0″2 for IRC +10216; Sect. 5.1), the brightness temperatures of these masers are roughly 1000−300 000 K. The map of the HCN J = 3−2 (0, 1^le, 0) emission indicates that the maser can have a smaller emission size of <0.″2 towards IRC+10216 (Shinnaga et al. 2009), which would lead to a higher estimate in brightness temperature. Accounting for the fact that the brightness temperature increases quadratically with decreasing source size and that the emission size may vary in other sources or transitions, the maser brightness temperatures are still too low to reach their saturation levels, even if we assume a source size that is ten times smaller. Unsaturated masers tend to vary on much shorter timescales than saturated masers due to their exponential amplification behaviour. This may explain the observed temporal variability on a relative short timescale of a few years (Sect. 4.3).

Most of the bright HCN maser components appear on the blueshifted side of their spectra with respect to their systemic LSR velocities. This may be explained by having the maser-emitting regions in the front of the stars, assuming an expanding outflow. The observed blueshifted maser may arise from the ‘radial’ amplification of background photospheric emission. This scenario has been invoked to account for the observed blueshifted SiS (1−0) maser in IRC +10216 (Gong et al. 2017). Stellar occultation may cause maser emission arising from the far side of the star to be partially blocked. Alternatively, the blueshift velocities of certain narrow maser profiles are not really significant and the maser spikes appear to peak very close to the systemic velocity. Prominent examples include RAFGL 4211, W Ori, and S Aur. These masers may be better interpreted as ‘tangential’ amplification, where the emission region is on the plane of sky beside the central star. This scenario requires an outward-accelerating velocity field and has been proposed to explain the detection of strong SiO masers near the systemic velocity in oxygen-rich stars (Bujarrabal & Nguyen-Q-Rieu 1981).

The J = 2−1 (0, 1^le, 0) transition is the brightest HCN maser in our observations (see Sect. 4.2.2). This could arise from the fact that the Einstein A coefficient strongly depends on the rotational quantum number J. The HCN molecules rapidly decay to the lower-J levels, and are likely to accumulate at the lowest-J levels. Direct l-type transitions have been detected in emission towards different carbon-rich stars (e.g. Thorwirth et al. 2003; Cernicharo et al. 2011). Emission from direct l-type transitions can cause molecules to cascade down from the (0, 1^1f, 0) level to (0, 1^1e, 0) level. Furthermore, the rotational transitions in the ν₂ = 1 state do not change the parity due to the selection rule. The net effects result in more HCN molecules in the (0, 1^1e, 0) state than in the (0, 1^1f, 0) state. The gain of masers is directly proportional to molecular column densities at their corresponding energy levels, which makes the J = 2−1 (0, 1^1e, 0) transition the brightest HCN masers in the ν₂=1 state. This may explain why HCN masers in the (0, 1^1e, 0) state are typically stronger or more readily detectable than the (0, 1^1f, 0) state, as shown in Figs. 11, 12, and D.1, which compare the (0, 1^1e, 0) and (0, 1^1f, 0) spectra in the same rotation transition at different J levels and in different sources.

The above pumping scenario is only consistent with our results for the J = 2−1 and J = 3−2 masers. For the J = 4−3 transition, it is also consistent with the results in RAFGL 4211, R Vol, W Ori (Fig. 11), and perhaps RV Aqr (Fig. D.1), but fails to explain the fact that (0, 1^1f, 0) maser is brighter than the (0, 1^1e, 0) masers in IRC +10216, II Lup, and CRL 3068 (Fig. 12). There are two common features in these sources showing the brighter (0, 1 ^1f, 0) maser in J = 4−3. First, the peak velocities of these J = 4−3 (0, 1^1f, 0) masers are all blueshifted by about 10 km s⁻¹ from their systemic LSR velocities, regardless of their discrepant ground-state or excited-state line widths (Tables 4 and 6). The peak velocities are different from those of the J = 2−1 (0, 1^1e, 0) masers of the same sources or the J = 4−3 (0, 1^1e, 0) masers of different sources, indicating a different pumping mechanism of this maser in these three sources. Second, all three sources show a broad emission profile in (0, 1^1f, 0), which is roughly a factor of two brighter than the corresponding, apparently thermal, (0, 1^1e, 0) profile. This indicates weak maser activities in the (0, 1^1f, 0) line. The J = 4−3 (0, 1^1f, 0) masers may stem from (a combination of) infrared line overlap, anisotropic infrared radiation fields, and/or non-negligible collisional pumping. Observations have shown a large amount of ro-vibrational HCN and C₂H₂ transitions at ~ 14 µm, and a large fraction of them are overlapped (Fonfría et al. 2008). The ro-vibrational HCN and C₂H₂ lines at 14 µm could partially contribute to the pumping of HCN masers in the v₂=1 state. Such an infrared line overlap has been invoked to explain the observed high-J SiS masers (Fonfría Expósito et al. 2006). Because the HCN and C₂H₂ lines are displaced in frequency, the overlap should give rise to the line asymmetries in the maser profiles. Furthermore, such a mechanism would lead to the presence of masers at nearly the same velocity shift with respect to their systemic velocities for different sources. This is in agreement with our observations, where the most prominent component of the J = 4-3 (0, 1^lf, 0) maser in these three sources peaks at about 10 km s⁻¹, blueshifted from the respective systemic velocities. Alternatively, the asymmetries in the maser profiles can arise from a highly anisotropic infrared radiation field, caused by substructures in the circumstellar envelope such as spirals or broken shells. Anisotropic infrared radiation field has been indicated by the non-spherical infrared morphology towards II Lup (Lykou et al. 2018). Thirdly, collisional pumping might not be negligible, although the vibrationally excited transitions are thought to be dominated by radiative pumping.

As pointed out in previous studies (Smith et al. 2014; Menten et al. 2018), the time variability does not appear to correlate with stellar cycles, indicating that collisional pumping may also contribute to part of the population inversion. This may explain the second feature where the J = 4–3 (0, 1^1f, 0) line exhibits a weak maser over a broad velocity range in these three sources. Inho-mogeneous density structures driven by stellar pulsating shocks may result in the diverse line shapes among the HCN masers.

The HCN J = 3–2 (0, 2, 0) and J = 3–2 (0, 3^1e, 0) masers are only detected in RAFGL 4211 among carbon stars as of today (Table A.1). The uniqueness of these two masers could stem from the peculiar geometry of RAFGL 4211. For instance, the relatively long coherent length of maser amplification along the line of sight may result in a uniquely high gain for this maser. Both J = 2–1 and J = 4–3 (0, 1^1e, 0) masers are detected, with the latter showing a narrow spike on top of a broad component. The higher-J maser lines (3–2 and 4–3) could partly arise from a region closer to the star, where the temperature and H2 number density are higher than in the region producing the J = 2–1 (0, 1 ^1e, 0) maser.

The physical conditions to excite the detected HCN masers can be roughly constrained. Based on the population diagram analysis (see Sect. 5.1), the excitation temperature is found to be around 800–900 K for the vibrationally excited HCN emitting regions. These masers have to be formed in these regions, which should have a kinetic temperature of ≳700 K. Furthermore, the velocity widths of the maser profiles indicate that the masers form in the wind acceleration zones. Comparing with the physical model in IRC +10216 (e.g. Agúndez et al. 2012), the expected gas temperature and velocity correspond to a H2 number density of >10⁸ cm³. This gives a loose constraint on the physical conditions.

Fig. 11 Overlaid spectra for HCN (0, 1le, 0) (in black) and (0, 1lf, 0) (in green) transitions in different rotational transitions towards RAFGL 4211 (left), R Vol (middle), and W Ori (right). These sources show brighter J = 4−3 (0, 1le, 0) emission than J = 4−3 (0, 1lf, 0). The vertical lines are labelled in the same way as for Fig. 1.

Fig. 12 Same as Fig. 11 towards IRC+10216 (left), II Lup (middle), and CRL 3068 (right). These sources show brighter J = 4−3 (0, 1lf, 0) emission than J = 4−3(0, 1le, 0).

5.3 Comparison with masers in oxygen-rich stars

It is known that OH, SiO, and H2O masers are common in oxygen-rich stars. SiO masers are formed in the dust formation zone, OH masers are formed farther out, and H2O masers lie between (e.g. Habing 1996). As discussed above, HCN masers should be formed in the innermost regions (i.e. r < 20 R_*), similar to the formation regions of SiO masers in oxygen-rich stars. Furthermore, both HCN and SiO are linear molecules that show maser actions in their ground vibrational states towards respective classes of stars (see Table A.1 for HCN ground-state masers towards carbon-rich stars, and see Boboltz & Claussen 2004; Dike et al. 2021 for SiO ground-state masers towards oxygen-rich stars). High-excitation masers with lower energy levels above 4000 K are detected in both SiO and HCN, with SiO masers in vibrational states up to v = 4 detected in VY CMa (Cernicharo et al. 1993) and up to v = 7 detected in Betelgeuse (Kervella et al. 2018); and hCn masers in the (0,4⁰,0) state were reported towards IRC +10216, CIT 6, and Y CVn (Schilke et al. 2000; Schilke & Menten 2003). These similar properties suggest that HCN masers in carbon-rich stars may be an analogue of SiO masers in oxygen-rich stars. If confirmed, HCN masers can be regarded as a powerful tool for measuring the accurate proper motions and parallaxes of carbon-rich stars as performed towards oxygen-rich stars with SiO masers (e.g. Zhang et al. 2012), especially towards the high mass-loss carbon-rich stars, which cannot be determined by optical and infrared telescopes due to their high extinctions. Very long baseline interferometry measurements of HCN masers might be the only way to determine their proper motions and parallaxes.

6 Summary and outlook

We performed APEX observations of multiple vibrationally excited HCN transitions towards a sample of 16 carbon-rich stars. The main results are summarised as follows:

1. Our observations resulted in the detection of 68 vibrationally excited HCN transitions from 13 stars, suggesting that vibra-tionally excited HCN lines are ubiquitous in carbon-rich stars. Among the detected transitions, 39 of them tend to be thermal, while 29 of them show maser actions. Our observations extend the spectrally resolved detection of the vibrationally excited thermal HCN transitions in carbon-rich stars beyond IRC +10216. Among the detected masers, 23 are new detections. The results expand the total number of HCN masers towards carbon-rich stars by 47%. The masers arise from 12 different transitions, and the J = 3–2 (0, 3^1e, 0), J = 3–2 (0, 2, 0), J = 4–3 (0, 1^1f, 0) masers are detected in the interstellar medium for the first time.

2. Temporal variations of the J = 2–1 (0, 1^1e, 0) maser on a timescale of a few years have been confirmed by the dramatic change of maser line profiles. However, the cause of the time variability of HCN masers is still unclear.

3. Population diagrams of vibrationally excited thermal HCN transitions suggest an excitation temperature of about 800 K for IRC+10216 and CRL 3068, consistent with the previous studies. Assuming a typical size of 28 au for the regions showing the vibrationally excited HCN transitions, we derived the source-average HCN column densities of (0.18–10.8) × 10¹⁹ cm⁻² for six carbon-rich stars. The corresponding HCN abundance with respect to H2 is roughly ~7.4 ± 5.5 × 10⁻⁵.

4. The ratios between HCN masers and corresponding infrared luminosities suggest that HCN masers are unsaturated. The unsaturated nature is supported by the temporal variation and relatively low brightness temperatures of these HCN masers. The spontaneous decay of rotational transitions and direct l-type transitions in the v₂=1 state may explain the fact that the J = 2–1 (0, 1^1e, 0) transition is the brightest maser in the v₂=1 state. Comparing with the physical model of IRC+10216, we infer that a gas kinetic temperature of ≳700 K and an H₂ number density of > 10⁸ cm⁻³ are required to excite the HCN masers. A comparison between masers in carbon-rich and oxygen-rich stars suggests that HCN masers in carbon-rich stars may be an analogue of SiO masers in oxygen-rich stars.

Our observations suggest that both vibrationally excited HCN thermal and maser line emission are rich in the CSEs of carbon-rich stars and have quite high flux densities. This indicates that these transitions are suitable for follow-up highresolution observations in order to study the entire class of carbon-rich stars. With the very long baselines, the Atacama Large Millimeter/submillimeter Array observations of these vibrationally excited HCN lines can achieve an angular resolution of ≲30 mas (e.g. Wong et al. 2016), which should be able to constrain the physical conditions of the extended atmospheres and inner winds of carbon-rich stars.

Acknowledgements

This publication is based on data acquired with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the MaxPlanck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. This research has made use of NASA's Astrophysics Data System. We thank A. V. Lapinov and the Editor-in-chief of the Astronomicheskii Tsirkulyar for kindly providing the scanned copies of Lapinov et al. (1994) and Zinchenko et al. (1988), respectively. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. This work has made use of Python libraries including NumPy (https://www.numpy.org/) (van der Walt et al. 2011), SciPy (https://www.scipy.org/) (Jones et al. 2001), Matplotlib (https://matplotlib.org/) (Hunter 2007). We thank the anonymous referee for their useful comments.

Appendix A Catalogue of known HCN masers towards carbon-rich stars

Based on our new detections and those in the literature, we have compiled the largest catalogue of HCN masers towards carbon-rich stars, and the results are listed in Table A.1. These masers are assigned to 36 different carbon-rich stars, and these HCN masers can be assigned to ten different vibrational states. Our observations led to the discovery of 23 new HCN masers, increasing the total number to 75.

Appendix B Estimates of the stellar variability phase

The phases of stellar variability on the dates of our 2018 observations were calculated for a number of sources in Table B.1, primarily using the information compiled by the AAVSO. In the used convention, phase zero is at the maximum light. While pulsation periods are available for the majority of sources (cf. Table 1), the epochs of maxima are scarce in the literature, limiting the number of stars with known phases at the epochs of our observations. In addition to AAVSO records, we used data from Menten et al. (2012, for IRC+10216), Feast et al. (2003, for IILup), Whitelock et al. (2006, for IRC+30374), and Menzies et al. (2006, for CRL–2513).

Appendix C Vibrationally excited HCN lines in other sources

In addition to the HCN spectra of IRC +10216, CRL 3068, and II Lup (see Figs. 1–4), the observed APEX HCN spectra in the other 13 sources are presented in Figs. C.1–C.13.

Fig. C.1 Spectra for J =2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards R Vol. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.2 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards Y Pav. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.3 Spectra for J = 3–2 (left) and 4–3 (right) transitions of HCN towards AQ Sgr. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.4 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards CRL 3068. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.5 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards IRC +30374. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.6 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards V636 Mon. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.7 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards IRC +20370. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.8 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards W Ori. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.9 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards S Aur. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.10 Spectra for J = 3–2 (left) and 4–3 (right) transitions of HCN towards V Aql. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.11 Spectra for J = 2–1 (left) and 3–2 (right) transitions of HCN towards CRL 2513. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.12 Spectra for J = 2–1 (left) and 3–2 (right) transitions of HCN towards CRL 2477. The description of the figure is the same as in the Fig. 1 caption.

Fig. C.13 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards RV Aqr. The description of the figure is the same as in the Fig. 1 caption.

Appendix D Comparison of HCN emission in (0,1^1e 0) and (0, 1^1f, 0) vibrational states

Fig. D.1 Same as for Figs. 11 and 12 towards AQ Sgr, V636 Mon, IRC+20370, S Aur, V Aol, and RV Aqr.

References

All Tables

All Figures

Fig. 1 From left to right: Spectra for J = 2−1 (first) and 3−2 (second and third) transitions of HCN towards IRC +10216. The vibrational state of each transition is indicated in the top right corner of the respective spectrum. The vertical dashed blue lines denote the range of the emission in VLSR, the vertical dotted line shows the systemic VLSR of the star, and the horizontal dotted line is the baseline. The number below the vibrational quantum numbers in the panel is a code assigned to the date of observation of the particular transition. The corresponding dates of the codes are listed in the notes of Table 4.

In the text

	Fig. 2 Spectra for J = 4−3 transitions of HCN towards IRC +10216. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. 3 Spectra for J = 2−1 (left), 3−2 (centre) and 4−3 (right) transitions of HCN towards RAFGL 4211. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. 4 Spectra for J = 2−1 (left), 3−2 (centre) and 4−3 (right) transitions of HCN towards II Lup. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. 5 Histogram distribution of expansion velocities derived from J = 3−2 ground-state HCN lines. The magenta line marks the median for HCN J = 3−2 from this work. The median values calculated by Ramstedt et al. (2009) towards M-type, S-type, and C-type stars are marked in orange, green, and red, respectively.
In the text

	Fig. 6 FWHM line widths as a function of upper energy levels for IRC +10216 (left) and CRL 3068 (right). The colours represent the rotational transitions and the symbols represent the various thermal vibrationally excited lines, as indicated in the legends.
In the text

Fig. 7 Variability of HCN (2−1) masers seen towards IRC+10216 (top two panels) in (0, lle, 0) and (0, llf, 0), II Lup (bottom left) and RAFGL 4211 (bottom right) in (0, lle, 0) vibrationally excited lines. The black spectra are from observations in 2018, cyan spectra are from 2021, red spectra are from Menten et al. (2018), and blue spectra are from Lucas & Cernicharo (1989) with the date of observations indicated in the legend.

In the text

	Fig. 8 Rotation diagram of IRC+10216 constructed from vibrationally excited HCN lines. Some masers are also plotted in the diagram, shown by triangles, but are excluded from the linear least-square fit.
In the text

	Fig. 9 Rotation diagram of CRL 3068. The description of the figure is the same as in the caption of Fig. 8.
In the text

Fig. 10 Energy diagram of HCN. The quantum numbers, υ and J, indicate the vibrational state and rotational state, respectively, l-type doubling causes the splitting of every rotational level into two sub-levels, le and 1f, which are also denoted in this figure. Pumping photons from the ground state to highly vibrationally excited states are indicated by cyan arrows. The J = 2−1 (0, lle, 0) and J = 2−1 (0, llf, 0) transitions are indicated by the blue and purple arrows, respectively. The direct l-type transitions are marked by red arrows. We note that the separations of the l-type doublets are not drawn to scale.

In the text

	Fig. 11 Overlaid spectra for HCN (0, 1le, 0) (in black) and (0, 1lf, 0) (in green) transitions in different rotational transitions towards RAFGL 4211 (left), R Vol (middle), and W Ori (right). These sources show brighter J = 4−3 (0, 1le, 0) emission than J = 4−3 (0, 1lf, 0). The vertical lines are labelled in the same way as for Fig. 1.
In the text

	Fig. 12 Same as Fig. 11 towards IRC+10216 (left), II Lup (middle), and CRL 3068 (right). These sources show brighter J = 4−3 (0, 1lf, 0) emission than J = 4−3(0, 1le, 0).
In the text

	Fig. C.1 Spectra for J =2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards R Vol. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.2 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards Y Pav. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.3 Spectra for J = 3–2 (left) and 4–3 (right) transitions of HCN towards AQ Sgr. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.4 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards CRL 3068. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.5 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards IRC +30374. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.6 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards V636 Mon. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.7 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards IRC +20370. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.8 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards W Ori. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.9 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards S Aur. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.10 Spectra for J = 3–2 (left) and 4–3 (right) transitions of HCN towards V Aql. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.11 Spectra for J = 2–1 (left) and 3–2 (right) transitions of HCN towards CRL 2513. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.12 Spectra for J = 2–1 (left) and 3–2 (right) transitions of HCN towards CRL 2477. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. C.13 Spectra for J = 2–1 (left), 3–2 (centre) and 4–3 (right) transitions of HCN towards RV Aqr. The description of the figure is the same as in the Fig. 1 caption.
In the text

	Fig. D.1 Same as for Figs. 11 and 12 towards AQ Sgr, V636 Mon, IRC+20370, S Aur, V Aol, and RV Aqr.
In the text