© ESO 2018

1. Introduction

The envelope of the asymptotic giant branch (AGB) star CW Leo, IRC+10216, is one of the most popular targets for molecule hunting in space (Cernicharo et al. 2000). The main reason for this is its proximity (~130 pc from the Sun) and high mass-loss rate [$\dot{M}\approx {4\times 10}^{-5}{(d/200\mathrm{p}\mathrm{c})}^2 M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$, (Kwan & Linke 1982; Jura 1986; Cernicharo et al. 2015)].

Dyck et al. (1991) have derived a stellar period of 638 days from observations in the K (2.2 µm) and M (4.64 µm) bands, while (Le Bertre 1992) reported it to be 649 days in Bands J(1.24 µm) and M. Shenavrin et al. (2011) reported on IR observations that (Menten et al. 2012) used to establish that CW Leo is a low-amplitude (~1 mag in the K band) Mira variable star with a period of 630 ± 2.9 days. Cernicharo et al. (2014) reported on the variability of molecular lines in the submillimeter and far-infrared range. The most striking changes were revealed by the CCH radical, which shows intensity fluctuations as large as a factor of 3, while neighboring lines from other species were almost unchanged. Pulsations of the central star and the subsequent IR pumping variations could be the primary explanation for this behavior. He et al. (2017) have reported the variability of SiS J = 14−13 and HCN v₂ = 1^f J = 3−2 (masing lines) with a similar period and phase as the IR light variation in this object, while (Fonfría et al. 2018) have performed an extended study that also included the SiS J = 15−14 line, using data obtained simultaneously with those discussed in this work.

To shed light on the molecular variability of IRC+10216, we started a 3–1 mm time monitoring with the IRAM-30 m telescope with the following aims: low relative calibration uncertainties (typically within ±5% at 3 mm, and scaled accordingly with frequency, see Sect. 3), adequate time sampling (between ~15 and 40 days), large frequency coverage (112 GHz), and high-frequency resolution (~200 kHz).

We present the first results of this monitoring in the 3 mm atmospheric window (where the effect of observational uncertainties is much weaker than at shorter wavelengths), focusing on the most striking effects, such as the clear evidence of periodicity in many thermal lines, similar to the periodicity revealed in the IR emission, and the phase lags between light curves of different species or rotational transitions of the same species. As discovered by Cernicharo et al. (2014), the radical CCH is one of the species showing a clear variability, together with other radicals such as CN and C₄H. All three species are chemically related to the cyanopolyynes HC_2n+1N, which in turn also show some variability. We focus on this chemical family here.

The millimeter-wave observational procedure and complementary IR data are presented in Sect. 2. The strategy for the overall analysis of the millimeterwave data of this study is described in Sect. 3. Results are presented and discussed in Sect. 4. Conclusions and follow-up plans are given in Sect. 5.

2. Observational millimeter-wave procedure and complementary IR data

The observations have been carried out on a regular basis with the IRAM-30 m telescope since JD = 2 457 180 (for convenience, hereafter we use the modified Julian Date MJD = JD– 2 450 000) once every ~14–40 days (depending on schedule and weather constraints) using the wobbler system (to obtain very flat baselines), the EMIR receivers in seven different tunings, providing a quite complete coverage between 80 and 273 GHz, and the BBC, WILMA, and FTS backends. The maximum frequency resolution achieved is 195.3 kHz in all bands, which translates into velocity resolutions from 0.73 to 0.21 km s⁻¹. Image band rejection is −13 DB, so that some strong lines would produce ghost features in the signal band, although they can be well identified and therefore avoided in this study.

So far, at submission/acceptance of this paper, 63 observing runs have been achieved from MJD = 7180 to MJD = 8166 (6 June 2015 to 2 February 2018), although some frequency settings are sometimes missing as a result of technical or meteorological reasons. The average time sampling in round numbers is therefore 16 days, which means ~0.025 stellar cycles. The time sampling is therefore quite satisfactory, although it has not been entirely regular. The largest gap of missing data occurred between MJD 7406 and MJD 7469 (63 days) because of several scheduled runs that could not be carried out due to extremely bad weather conditions. In addition to these observations, we include in our data set eight additional observations achieved at MJD = 6107, 6128, 6257, 6462, 6488, 6538, 6581, and 6604 that are used in this work to check the consistency of some conclusions.

All the observations were reduced in the same way by marking off all detected features above a 5σ limit, fitting a degree−1 baseline to the unmarked channels, and subtracting it from each individual spectrum (~4 GHz wide). The results were stored and then resampled to build data cubes (X = Frequency, Y = MJD, $Z=T_{\mathrm{A}}^{*}$) and time-average spectra that were then used to present and analyze the results.

A global inspection using the data cubes allowed us to identify specific observing runs and frequency settings with potentially bad data that were then carefully checked out. Obvious cases appear to be due to one or a few corrupted scans that were then removed, and the affected frequency setting was reanalyzed. The main and most difficult sources of errors in our monitoring are, however, systemic (telescope) flux gain drifts as well as calibration, pointing, and focus errors. The telescope maintenance and regular flux monitoring measurements, achieved by the telescope staff, are made to ensure users that no significant gain drifts will affect data taken at different epochs. Pointing and focus were directly checked and corrected by us every 60–90 min using the nearby quasar OJ287, so that they were kept within 3″ and 0.2 mm accuracy. Hot-cold and sky calibrations are achieved every 5–7 min. In this paper we show and discuss only data corresponding to the 3 mm receivers (up to 116 GHz) where the pointing, focus, and calibration uncertainties are much smaller than at higher frequencies.

In order to have a reference on the variability of this object, we requested V. I. Shenavrin for the entire set of IR observations of CW Leo in bands H, J, K, L, and M achieved between MJD 1522 and MJD 5494 (Shenavrin et al. 2011). These observations were discontinued after the latter date. Therefore, there is no overlap between them and our own survey. Nevertheless, the IR data of Shenavrin allow a very good determination of the stellar pulsation period and give a valuable reference for studying variabilities in molecular lines. The period, Π^IR, that we established from fitting the IR data (all five bands) to sinusoids (a more sophisticated Fourier analysis is not possible because the number of data points is limited), is Π^IR = 636.28 ± 4.18 days. From the same fit we obtained an extrapolated (to the MJD = 7180–8166 range) value of the Phase 0 (fitted IR curves equaling their averages while the brightness increases), ${\tau }_0^{\mathrm{I}\mathrm{R}}$ of 7600.83 ± 10.38 days. The error bars given in these two parameters correspond to the 3σ values provided by the numerical fitting routine.

In the next section we define systematic routines that are used to search for variabilities in the observed molecular lines. The values (${\mathrm{\Pi }}^{\mathrm{I}\mathrm{R}},{\tau }_0^{\mathrm{I}\mathrm{R}}$) are in this process always used as a reference.

3. Strategy for 3 mm wavelength data analysis

Instead of studying selected millimeter-wave lines in the IRAM- 30 m time variability survey of IRC + 10216 that display the behavior we are searching for even for the naked eye, we conducted a global analysis based on general criteria with the ultimate goal of determining interesting physical effects without favoring any frequency range or molecular line over the others. To do this, we proceeded in several steps.

In order to determine whether we achieved a good absolute calibration through the full length of the monitoring, we defined a tolerance in the measured flux at each epoch τ of a line centered at frequency v_l as follows:$$\frac{\left|({\int }_v{T}_{\mathrm{A}}^{*}\left(\tau \right)-{\int }_v\langle T_{\mathrm{A}}^{*}\rangle )\right|}{{\int }_v\langle T_{\mathrm{A}}^{*}\rangle )}\le {\eta }_{\mathrm{r}\mathrm{e}\mathrm{f}}\left[v_l\left(\mathrm{G}\mathrm{H}\mathrm{z}\right)/80\right]\mathrm{.}$$(1)

By setting a value of η_ref, we assumed that all line flux fluctuations below the resulting tolerance at frequencies around v_l are within the flux variations caused by calibration, pointing, and focus uncertainties. The definition of this tolerance as proportional to the frequency is based on our observing experience at the IRAM 30 m telescope. Pointing, focus, and atmospheric uncertainties seem to combine to behave in this way. The selected value for η_ref is 0.04. In our data set we identified lines with flux variations within this tolerance, for which we considered that there is no real variability. These stable lines are our stronger argument to claim that observed variabilities in other lines, well beyond this tolerance, are real. This chosen η_ref value is rather demanding as it is widely accepted that absolute calibration uncertainties below ±10% are very difficult to achieve in single-antenna millimeter-wave observations.

Figure 1 shows one of these reference lines (¹³CO J = 1–0) that have an integrated line area within the defined tolerance for more than 90% of the observations. Many other lines do not show variability, but this line is the best and most intense example. The ¹³CO J = 1–0 data are very important in this paper as they show that the variabilities displayed by other lines are real and that we can draw important physical conclusions from their study. Despite the stability of the reference lines, we should point out that some changes appear in the emission that correspond to the terminal velocities of this object because of its pulsation. The period of these changes seems to be the same as in the IR and will be studied in detail elsewhere.

Fig. 1. [X = Frequency, Y = MJD, $Z=\mathrm{\Delta }{T}_{\mathrm{A}}^{*}$] 3D light curves of selected 3 mm transitions of C4H, 13CO, and CN in IRC+10216. Each panel also displays the average spectrum through the monitoring (black spectrum in the bottom box [magenta for 13CO J = 1–0]), observed minimum (in blue) and maximum (in red) spectra (except for 13CO J = 1–0), normalized line area (black dots superimposed on the color-coded 3D light curves [magenta dots for 13CO J = 1–0 in the left panel]), and sine function fit to these black dots (thick dashed black curve). The right panel additionally shows the IR Band M light curve (red-blue dashed sinus curve, red upper X-axis) obtained from fitting (Shenavrin et al. 2011) data (shown here as red and blue dots brought from their original observation dates to the MJD range of the figure). Phase delays are also shown.

At this point, it is important to say that the time-averaged spectrum $\langle T_{\mathrm{A}}^{*}\left(v\right)\rangle$ introduced above was built using equal weight for all data sets, although the observation time and system temperatures (typically in the range ~90–250 K in the 3 mm settings) can differ between observing runs. This is the only reasonable way to achieve this if the flux variations caused by calibration, pointing, and focus uncertainties are kept under control (as shown by the reference lines). Otherwise, the time-averaged line would be biased toward the data sets with lower noise levels. The only (limited) effect is for the noise itself in the final spectrum, but what matters is the line intensity, regardless of the noise. The time-averaged spectrum appears as an in-depth integration from which many features arise that are marginal in each individual epoch. However, since the receiver tuning configurations remain the same throughout the monitoring, a few ghost features from the image side band are not canceled out, although they can be identified.

In the time-averaged spectrum, we marked the central frequencies of all features with more than 5σ and then calculated the line area enclosed in the range ±18 km s⁻¹, around each frequency for each individual spectra. This range was chosen to enclose the terminal velocities of IRC+10216 alone [±14.5 km s⁻¹, (Cernicharo et al. 2000)], and avoid picking up signal from other features, and it was adapted in a few cases of blended lines.

We calculated the ratio of the area enclosed within the frequency interval Δv around the reference frequency v_c at each epoch τ_i, with respect to the same area for the average spectrum. If more than 60% of the resulting ratio points,$$\mathrm{\Gamma }(v_c,\mathrm{\Delta }v,{\tau }_i)=\frac{{\int }_{\mathrm{\Delta }v}{T}_{\mathrm{A}}^{*}(v,{\tau }_i)\delta v}{{\int }_{\mathrm{\Delta }v}\langle T_{\mathrm{A}}^{*}\left(v\right)\rangle \delta v)},$$(2)

fall outside the tolerance range, then we considered the feature as potentially variable and proceeded to fit a sinusoidal function Γ(τ) (we replaced by Δv its equivalent velocity interval Δυ),$$\mathrm{\Gamma }(v_c,\mathrm{\Delta }\upsilon ,{\tau }_{})=1+\frac{\mathcal{B}(v_c,\mathrm{\Delta }\upsilon )}{\mathcal{A}(v_c,\mathrm{\Delta }\upsilon )}\mathrm{.}\sin \frac{{\tau -\tau }_0(v_c,\mathrm{\Delta }\upsilon )}{\mathrm{\Pi }(v_c,\mathrm{\Delta }\upsilon )},$$(3)

from which we obtained for each (v_c, Δυ) experimental values for parameters 𝒜 (average area), ℬ (variability amplitude), τ₀ (phase 0), and Π (variability period). The sigma associated with each parameter in this fit was used to check the significance of the parameter. If it was below 15% for parameter A and below 50% for parameter B, we considered the sinus fit acceptable, otherwise we considered that the dispersion of the Γ_i points was due to a low signal-to-noise ratio rather than to variability. On the other hand, if more than 60% of the Γ_i points were within the established tolerance, we considered those spectral features as reference as well.

4. Results

A procedure for a continuous inspection, as the survey continues, of the 3D [X = Frequency, Y = MJD, $Z=T_{\mathrm{A}}^{*}$] data cubes for any desired spectral line has been implemented. The procedure automatically generates figures that display key information in a compact and understandable manner. We here include the very specific Fig. 1 with panels devoted to lines of CCH, C₄H, CN, and ¹³CO that display a well-defined variability. A general view of the results is provided in Table 1, Fig. A.1 and Table A.1.

The 3 mm transitions of the radicals CCH, C₄H, and CN exhibit the clearest periodic variability with the highest contrast of all lines we surveyed. The variability in CCH lines has been pointed out in (Cernicharo et al. 2014). C₄H displays a behavior similar to that of CCH, the periods are almost the same and the variability is in phase for both of them.

The CN lines also display a clear periodic variability, with a similar period, but almost completely in antiphase with the CCH and C₄H lines. In fact, most variable lines in this survey are in phase with CN, i.e. they display a phase lag of about 70 days with respect to the IR curves (Fig. A.1).

The cyanopolyynes HC₃N and HC₅N have several rotational transitions that are observable with the 3 mm receivers. These lines also show variability, but with lower contrast than in the case of the radicals. The periods derived are similar, and the phase of the variability is coincident with that of CN (and opposite to CCH and C₄H).

5. Conclusions and follow-up

For the first time, this work reveals 55 spectral features in the 80–116 GHz range displaying a variability behavior in IRC+10216 with a period that is compatible within the error bars with the average period derived from the IR bands. Therefore, the 3 mm data indicate that variability in the radiation field arising from the central star and dust formation zone affects many molecular species, although the observed contrast in this frequency range is relatively low compared with higher-J lines of molecules such as CCH (Cernicharo et al. 2014). Pointing and focus errors, as well as atmospheric effects, can be ruled out as the main source of the observed flux fluctuations since some reference lines remain quite constant with time, and other lines clearly vary following a smooth pattern and do not vary randomly.

The most striking result of this survey is that molecules with CN (i.e., cyanopolyynes and their radicals) have a common phase 0, while polyyne radicals are in antiphase with them. In addition, neither is in phase with the IR as a delay with respect to the IR of ~60–70 days is the most common. This behavior has also been found by Fonfría et al. (2018) in the masing emission of SiS J = 14–13 and J = 15–14.

In addition to the effect of the central radiation field (which is responsible for the variability period), some dynamical effects result in line profile changes, as well as in internal peculiarities in the ro-vibrational spectrum and pumping mechanisms of different species that introduce the reported phase lags. As an example, we could consider the differences in IR pumping for diatomic, triatomic, and polyatomic species. All diatomics will absorb and emit in the IR in a narrow frequency domain (typically that corresponding to the v = 0 ↔ 1 transition), whereas a triatomic molecule could have absorption around three frequencies, in particular in the bending mode, which involves the lowest energies and frequencies. Polyatomic molecules will have several bending modes. Hence, the IR pumping could have some differences for these three groups of molecules. For HNC, Cernicharo et al. (2014) have shown that although all vibrational modes participate in the pumping of the ground-state rotational levels, the bending mode stronglydominates the IR effects. For other species, with particularly low energy-bending modes such as C₄H and C₆H, collisions could also participate in the population of the rotational levels within vibrationally excited states. That C₂H and C₄H show a behavior different from all other species is therefore not a chemical effect, but probably a different IR pumping scheme. Absorption of IR photons from the ground to vibrationally excited states depopulates low-J levels in favor of high-J levels if the upper vibrational level is unpopulated. However, if its population due to collisions or to a significant permanent IR pumping is large, then the population of the levels of the ground state will follow a much more complex pumping scheme. It is worth noting that C₄H v₇ is easily detected and its phase is similar to that of all other species, while its ground state is in phase opposition. A detailed analysis is beyond the scope of this paper and will be presented in a forthcoming paper.

These results should guide future observations with ALMA as merging the data from different epochs could be a delicate matter. Table A.1 shows that contrasts exceeding 1.5 are hard to find at 3 mm, but they can be much larger at shorter wavelengths. Moreover, for both single-dish and interferometric observations, new radiative transfer models are required to take the variation in population of molecular levels with time into account.

Acknowledgments

This work has been supported by ERC under grant ERC-2013-Syg-610256-NANOCOSMOS, Spanish MINECO under grants AYA2012-32032, AYA2016-75066-C2-1-P and the CONSOLIDER-Ingenio program “ASTROMOL” CSD 2009-00038.

Appendix A:

Fig. A.1. Time-average IRC+10216 spectrum currently achieved in this project at the 3 mm atmospheric window, with periods and phases derived for spectral lines following the procedure described in Sect. 3. The error bars on the periods are also displayed. The sizes of the blue dots for each line are proportional to the percentage of the observed epochs where the flux exceeds the defined tolerance η(vl) (see Sect. 3). Larger dot sizes and smaller error bars mean more and better constrained periodic variability. The period and phase 0 derived from the IR data are plotted as a blue dashed and red dotted lines, respectively, for comparison. The species C2H and C4H are in antiphase with respect to the rest of molecules, which in turn display a delay of ~60–70 days with respect to the average phase 0 of the IR bands.

References

All Tables

All Figures

Fig. 1. [X = Frequency, Y = MJD, $Z=\mathrm{\Delta }{T}_{\mathrm{A}}^{*}$] 3D light curves of selected 3 mm transitions of C4H, 13CO, and CN in IRC+10216. Each panel also displays the average spectrum through the monitoring (black spectrum in the bottom box [magenta for 13CO J = 1–0]), observed minimum (in blue) and maximum (in red) spectra (except for 13CO J = 1–0), normalized line area (black dots superimposed on the color-coded 3D light curves [magenta dots for 13CO J = 1–0 in the left panel]), and sine function fit to these black dots (thick dashed black curve). The right panel additionally shows the IR Band M light curve (red-blue dashed sinus curve, red upper X-axis) obtained from fitting (Shenavrin et al. 2011) data (shown here as red and blue dots brought from their original observation dates to the MJD range of the figure). Phase delays are also shown.

In the text

Fig. A.1. Time-average IRC+10216 spectrum currently achieved in this project at the 3 mm atmospheric window, with periods and phases derived for spectral lines following the procedure described in Sect. 3. The error bars on the periods are also displayed. The sizes of the blue dots for each line are proportional to the percentage of the observed epochs where the flux exceeds the defined tolerance η(vl) (see Sect. 3). Larger dot sizes and smaller error bars mean more and better constrained periodic variability. The period and phase 0 derived from the IR data are plotted as a blue dashed and red dotted lines, respectively, for comparison. The species C2H and C4H are in antiphase with respect to the rest of molecules, which in turn display a delay of ~60–70 days with respect to the average phase 0 of the IR bands.

In the text