The Maser-emitting Structure and Time Variability of the SiS Lines J = 14–13 and 15–14 in IRC+10216^*

The CARMA observations were carried out on 2011 January 5 in its 15-antenna mode under the frame of project c0710 (Fonfría et al. 2014). The antennas were arranged in the B array configuration, with baselines in the range of 70–830 kλ and system temperatures of about 200 K. The weather during the observations was good, with an atmospheric opacity ≃ 0.1 at 230 GHz.

The calibration and data reduction were performed in the standard way using the Miriad package. Bandpass solutions were obtained observing 3C 84, while 0854+201 was periodically monitored every 15 minutes to calibrate phases. The observing track went on during more than 5 hr, and we integrated on-source for more than 3 hr, so that the atmospheric variations are expected to be averaged in the final maps. The phase calibrator was at an angular distance of ≃15° from IRC+10216, which could in principle introduce an error in the phase calibration that is difficult to quantify, although it should be small, as the weather was good. Calibrator 0854+201 was assumed to have a flux density of 3.8 Jy at 257 GHz, as determined by the CARMA and Submillimeter Array (SMA) quasar monitoring programs near the time of our observations. The absolute flux calibration uncertainty is about 15%.

The SiS (14–13) line was observed using a spectral window with a bandwidth of 125 MHz and a spectral resolution of 0.78 MHz (≃0.92 km s⁻¹ at the observed frequency). This observation was taken at the same time as the lower spectral resolution data of this line presented by Fonfría et al. (2014). The obtained maps have a synthesized beam of about 025, resolving partially the structure around the central star. The primary beam of the telescopes was 30''. The dirty point-spread function (PSF) shows a few significant sidelobes with a peak of ≃30% located at ≃1'' from the center of the main beam. These sidelobes are not expected to significantly affect the shape of the lowest-level contours of the brightness distribution considered in this work.

The noise rms determined from the statistical analysis of a velocity channel with no emission is ≃45 mJy beam⁻¹. Following the procedure used by Fonfría et al. (2014), we derive a statistical position uncertainty of ≃2 mas. We estimate the systematic uncertainty to be ≃30 mas after analyzing the rms of the visibility phases of calibrator 0854+201 and the position of the peak continuum emission of IRC+10216 in each observing cycle. Hence, the position uncertainty is clearly dominated by the systematic error. With this error, the derived stellar position is compatible with the very accurate results by Menten et al. (2012).

The ALMA observations were taken on 2012 April 8–23 during Cycle 0 (project 2011.0.00229.S; Cernicharo et al. 2013). The array consisted of 16 antennas with baselines from ≃30 up to 370 kλ. The weather was good, resulting in system temperatures ranging from 100 to 130 K.

The calibration and data reduction were performed in the standard way using the CASA software. Bandpass and absolute flux were calibrated with observations of 3C 279 and 3C 273. J0854+201 and J0909+013 were observed every 10 and 20 minutes to calibrate phase and amplitude. The absolute flux calibration uncertainty is about 8%. See Cernicharo et al. (2013) for a deeper insight into the calibration process.

The observations covered a spectral range of ≃270–274 GHz with a channel width of ≃0.49 MHz (≃0.54 km s⁻¹ at the frequency of the SiS line J = 15–14) and an effective spectral resolution of ≃1.1 km s⁻¹. The data were weighted assuming a robust parameter of 0.5, which produced maps with an angular resolution of ≃055.

The noise rms is ≃5 mJy beam⁻¹. Following the same method to estimate the position uncertainty as for the CARMA observations, we find that the statistical contribution is negligible. Again, the systematic uncertainty dominates and the position error can be considered as ≃80 mas.

In this work, we have used the IRAM 30 m telescope observations of the J = 14–13 and 15–14 SiS lines carried out on 2004 June 19 and published by Fonfría Expósito et al. (2006). We also analyzed a new set of SiS (i.e., ²⁸SiS) and ²⁹SiS lines recently observed as part of a more ambitious project targeting the time variability of the IRC+10216 molecular emission. This monitoring, motivated by the discovery of the time variability of the thermal line emission by Cernicharo et al. (2014), was carried out from 2015 June 7 to the present day with an average sampling time of 16 days. We use the wobbler system to get very flat baselines, the EMIR receivers in seven different tunings, and the BBC, WILMA, and FTS back ends. All this allows us to completely cover 84 GHz within the 80–273 GHz range in each observing run.

Pointing and focus are checked every 90–100 minutes using the nearby quasar OJ 287 (0854+201). The maximum frequency resolution achieved is 195 kHz in all bands, which translates into velocity resolutions from 0.21 to 0.73 km s⁻¹. Typical integration times range from ≃20 to 50 minutes, which ensures ${T}_{A}^{* }$ noise levels below 20 mK for each individual observation (just a few mK at the lowest frequencies). The whole data set will be presented and analyzed in Pardo et al. (2018).

The data imaging of the CARMA and ALMA data sets and the single-dish data reduction were achieved using the GILDAS software.⁸

The comparison between the spectra taken at the central pixel of lines J = 14–13 and 15–14 shows many similarities but also prominent differences (Figure 1). The SiS (14–13) line is mostly composed of four emission components (A, B, C, and D), a deficit of emission (W; see also Figures 2 and 3), and a weak blueshifted absorption mostly noticeable below the continuum in the interferometer observations. The fingerprints of all these components can be found also in line J = 15–14, although two new components, named A* and D* from now on, are noticeably blue- and redshifted with respect to components A and D. It is remarkable that the velocity of the emission peak of line J = 14–13 does not match up with that of the emission peak shown by the same line observed with the IRAM 30 m telescope, contrary to line J = 15–14 (ALMA), which shows the same peak velocity as its single-dish counterpart. It is not clear whether the A and D components comprise the A* and D* components or are just separate components (the ALMA data do not have enough angular resolution to spatially disentangle the emission), strongly affected by the time evolution of the envelope. In addition to the easily visible component A*, hints of components A, B, C, D, D*, and W and the blueshifted absorption can also be found in the single-dish observations acquired with the IRAM 30 m telescope reported by Fonfría Expósito et al. (2006) and the new data set taken at different stellar pulsation phases (Section 3.2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The moment 0 map of the SiS (14–13) line shows a partially resolved compact brightness distribution roughly centered on the star (Figures 2 and 3; Table 1). The lower-level contours are mostly circular. There is an emission deficit to the E of the star that is probably caused by a negative sidelobe of the synthetic PSF close to the main beam, smaller than 10% and arranged along the direction with P.A. ≃ 120°. The higher-level contours are elliptical, with their major axes oriented roughly along the NE–SW direction and a ratio of the minor to the major axes of ≃80%, suggesting that the brightness distribution is elongated along the direction with P.A. ranging from 50° to 70°. The emission peaks at 007 ± 004 to the NW of the star (Fonfría et al. 2014).

Zoom In Zoom Out Reset image size

The position–velocity diagrams show that the SiS maser emission mostly comes from the four relatively isolated components A, B, C, and D presented above (Figures 2 and 3). The strongest component, A, with a peak flux of ≃16.9 Jy beam⁻¹, is located in front of the star at −5.5 km s⁻¹ with respect to the systemic velocity. The following component in brightness, D, displays a maximum flux of ≃8.5 Jy beam⁻¹ and is located behind the star between 8.5 and 10.5 km s⁻¹. The bulk emission of both components A and D comes from the northern hemisphere of the envelope. The resemblance between the position in the plane of the sky and the shape of components A and D (blueshifted and redshifted with respect to the central star) suggests that they are similar emitting structures mirrored with respect to the plane of the sky containing the star.

Component B is a weaker structure around the star elongated along the direction with a P.A. between 70° and 90°, showing an average flux of ≃35% of the cube peak emission (≃16.9 Jy beam⁻¹). Its size is ≃05 × 02 in the plane of the sky and ≃3 km s⁻¹ along the velocity axis. The maximum flux of this structure (≃45% of the peak emission) is found at ≃01 to the E of the star, about 3 times the estimated position uncertainty. The position–velocity diagrams show the structure corresponding to component W, which is located at ≃3.5 km s⁻¹ and where the emission is much weaker than at other velocities. Between component W and the central star is component C, weaker than component B. The whole structure is embedded in an extended, roughly uniform, and weak structure with an emission ≲20% of the cube peak emission.

The integrated emissions coming from every component A, B, C, and D are mapped in Figure 4. The position of their emission peaks and their corresponding elliptical Gaussian fits can be found in Table 1. Regarding the morphology of components A and D, since their lowest contour levels are above the 3σ level, the differences observed in them suggest a slightly different structure of the emitting regions along the E–SE direction. This difference might evidence an anisotropy of the gas density as long as an elongation roughly along the SE direction is also noticeable in the H¹³CO and SiC₂ maps presented by Fonfría et al. (2014). Nevertheless, this elongation is not present in the SiO maps, which could imply that chemistry also plays a role in the SE quadrant of the envelope.

Zoom In Zoom Out Reset image size

Component B is substantially different from components A and D. The lowest-level contours are roughly circular and are centered on the star. However, component B clearly shows a bipolar structure with two peaks at a distance of 013 ± 004 and 017 ± 004 from the star (6.5 ± 2.0 and 8.5 ± 2.0 ${R}_{\star }$) with P.A. ≃ 70° and 230°, respectively. The peak along the direction with P.A. ≃ 70° accounts for 65% of the total emission of this component. It is worth noting that the emission peaks of components A, C, and D lie along the N–S direction, roughly perpendicular to the direction defined by the peaks of component B (NE–SW).

There is a remarkable resemblance between the structure shown by component B and the brightness distribution of the NaCl (21–20) line in the vibrational ground state in the channels around the systemic velocity presented by Quintana-Lacaci et al. (2016; see their Figure 2). The most evident difference between both structures is the position of their emission peaks, which are ≃03 ≃ 15 ${R}_{\star }$ from the star in the NaCl line. This supports the idea that the gas is anisotropically distributed around the star, as was already suggested by Fonfría et al. (2014) after analyzing the high angular resolution emission of H¹³CN, SiO, and SiC₂.

Component C is also different from the other components. It shows an overall circular morphology avoiding the central region and two emission peaks to the N and S. The southern peak is close to the SW maximum displayed by component B, suggesting that the structure of both components could be somehow related. Moreover, the southern emission of this component roughly peaks at the mirrored position of the emission peaks of components A and D with respect to the direction defined by the peaks of component B (see Figure 4). On the other hand, the northern emission peak matches up quite accurately with the emission peaks of components A and D. This could mean that all these components are related to a common structure in the surroundings of the central star.

Figure 5 shows the average Doppler velocity of the observed brightness distribution (moment 1 map) with variations not larger than ≃1.5 km s⁻¹. It is mostly redshifted with respect to the systemic velocity. The average emission that comes from the NW, which coincides with the emission peaks of components A and D, is slightly redshifted with respect to that coming from the south. This means that the emission is roughly symmetric regarding the plane that contains the central star, but the emitting regions behind the star move faster related to Earth than those in front of it. This is particularly evident to the NW of the star, where the emission of component D is weaker than that of component A (Table 1).

Zoom In Zoom Out Reset image size

Figure 6 shows the observations of the SiS v = 0 J = 5–4, 6–5, 14–13, and 15–14 lines and ²⁹SiS v = 0 J = 13–12 and 15–14 lines mentioned above spanning roughly through a single pulsation period. All the data have been arranged in six columns, one per molecular line. In a given column, we have plotted the difference of the corresponding line observed at every pulsation phase with respect to the average line over the whole monitoring.

Zoom In Zoom Out Reset image size

The ²⁹SiS lines show what seems to be a long-term time variation. There are also departures from the long-term behavior that suddenly affect the lines observed at isolated pulsation phases or during a few adjacent phases, e.g., in the ²⁹SiS lines at phase 0.30 or in the SiS (5–4) line at phases 1.10 and 1.25. Since these observations have been corrected from calibration issues (see Appendix A), all of these short-term departures are in principle real.

The J = 5–4 and 6–5 SiS lines show a roughly complementary dependence with time for pulsation phases below 0.8. Above this pulsation phase, however, it is not clear if their time dependences are or are not related. The long-term variation of the SiS (5–4) line seems to show the same period as the NIR light curve and a phase difference of about 0.2 (see Section 5.4), although it displays short-timescale departures from it mostly in the central part of the line. It is worth keeping in mind that these lines were not corrected from pointing or calibration errors, as we explain in Appendix A. The presence of narrow peaks at high expansion velocities in the SiS (5–4) line and their absence in the SiS (6–5) line throughout the whole monitoring are noticeable. In addition, these peaks vary strongly with the pulsation phase, much more than the rest of the line, contrary to what happens with the J = 6–5 line, for which the line profile is approximately constant. These facts suggest that the level J = 5 varies with the pulsation phase in a different way than the levels J = 4 and 6, as was already suggested by Carlström et al. (1990), who claim this difference as due to a possible IR overlap between an SiS line involving the level J = 5 and a line of another molecule.

The profile of the SiS (6–5) line resembles an inverted U, a shape typically associated with optically thick lines coming from spatially resolved regions (e.g., Zuckerman 1987). However, the SiS (5–4) line displays a double-horned profile. Considering that they are formed roughly in the same region of the envelope and the main beam of the 30 m telescope is nearly the same for their frequencies (HPBW ≃ 25''), we suggest that the high-velocity SiS (5–4) peaks are masers in nature, as for the peaks observed in the SiS lines J = 14–13, 15–14, and 1–0 (see below and in Henkel et al. 1983; Gong et al. 2017).

The strongest time variation is shown by the maser peaks of the J = 14–13 and 15–14 SiS lines (labeled as A* and D* in Figure 1). The emission of both lines seems to be phase coherent and also coherent with that shown by the SiS (5–4) line following the NIR light curve during most of the monitoring, although they show a few short-timescale departures from the long-term dependence. As occurs with the SiS (5–4) line, the long-term variation of the SiS lines J = 14–13 and 15–14 shows also a phase difference of ≃0.2 with respect to the NIR light curve.

Figure 7 shows the relative flux of the narrow peak, which is component A*, with respect to the total emission measured in the same velocity range in the SiS lines J = 14–13 and 15–14 against the NIR pulsation curve. The flux below the feature was calculated taking advantage of its small width, which let us determine the shape of the line related to the extended emission. The very good time sampling of a pulsation period achieved during our monitoring allowed us to describe accurately the variation undergone by this component. This variation is periodical, as it occurs with the NIR luminosity. Nevertheless, as we noted above, the SiS and NIR curves are not in phase. We have fitted the data sets with offset sines, finding phase differences of 0.179 ± 0.006 and 0.164 ± 0.007 for lines J = 14–13 and 15–14, respectively, with respect to the NIR light curve. Moreover, the maser light curves cannot be described accurately with a simple trigonometric function: the slope of the dependence after the minimum is steeper than during the decreasing intensity phase. We will discuss the origins of this phase difference in Section 5.4.

Zoom In Zoom Out Reset image size

The observed emission cannot be modeled assuming only positive excitation temperatures. The synthetic emission is far from reproducing the observed one even after assuming excitation temperatures as high as 10⁵ K. The observational uncertainties cannot explain this difference. This fact forced us to consider maser-emitting components in different places of the envelope model that are included only if it is absolutely necessary. Following this strategy, the most important features of the observed brightness distribution were reasonably well reproduced (Figures 8–10). Most of these components are masing clumps elongated along the line of sight that span through different envelope shells. The size of these clumps in the plane of the sky derived from our model is ≃1.5 ${R}_{\star }$ on average. There is also an arc-like masing structure that could be thought to comprise several very close clumps. Although the derived size distribution for the clumps seems to reproduce reasonably well the observations, we note here that higher angular resolution may provide more accurate results. In any case, the actual compact emitting regions are small, which is in good agreement with the morphology of the typical maser-emitting regions found in other AGB stars (e.g., Soria-Ruiz et al. 2007). In addition to the set of clumps, we need an extended maser structure to reproduce the observed extended emission (see below). The high number of different channels to model, each showing complex emission well above the rms noise level, explains the high number of clumps needed to reproduce the observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best way we found to model the extended emission deduced to exist from the observation includes a maser-emitting shell with a radius of 13 ${R}_{\star }$ and a thickness of ≃2 ${R}_{\star }$, i.e., showing an approximate size of $30\,{R}_{\star }\simeq 0\buildrel{\prime\prime}\over{.} 6$. It is not possible to estimate the actual thickness from our data because of the lack of angular resolution (see Section 5 for a discussion about this). We do not expect a bigger structure because, in this case, we would have a contribution to the emission at larger distances from the star that is nonexistent in our observations. We modified the rotational temperature at any position in this shell when needed to reproduce the observation. In fact, the modeling procedure resulted in a maser-emitting shell that is stronger in the velocity channels around 8.0, −0.5, and −5.0 km s⁻¹ and weaker at 11.0, 3.5 (component W), and −9.5 km s⁻¹.

Component A cannot be reproduced with a single clump located at the emission peak because it results in a significantly narrower component than observed. The best fit suggests that component A is mainly described by one arc spanning different velocity channels to the NW of the star (see Figure 10 in the range [−6.9, −3.2] km s⁻¹), which accounts for the emission peak and its close environment. This arc extends spatially from ≃6.5 to 11 ${R}_{\star }$ from the star, and its P.A. varies from ≃300° to 360°. Its emission is more blueshifted to the N of the star than toward the NW. Additional weaker emitting clumps are necessary to reproduce the rest of the A emission in the velocity range [−10.5, −1.5] km s⁻¹. These clumps are placed at different locations of the envelope between 5 and 20 ${R}_{\star }$, i.e., in Region II. In particular, the clumps that reproduce the emission at −1.5 and −10.5 km s⁻¹ are very close to the boundaries of this region. Most of them span up to four velocity channels, although one of them, located at the SE of the star, interestingly spans out 11 channels (i.e., ≃10 km s⁻¹), taking part in component B as well. However, with such a large velocity span, the growth of the emission of this clump along the line of sight cannot be exponential because in this case the emitted flux would be huge, contrary to what is really observed. This means that this structure is divided into several shorter clumps radiatively disconnected. This fact suggests that the coherent velocity is several times lower than 10 km s⁻¹, as expected to occur for the expanding gas close to the star, with a turbulent velocity ≲5 km s⁻¹ (Monnier et al. 2000; Agúndez et al. 2012; Fonfría et al. 2015, 2017).

In Table 3 we have included the fraction of the total flux of each emitting structure accounted for by every component (thermal emission, extended maser emission, and clumps/arcs). Despite the brightness density of the arc and that some clumps in component A are quite strong compared to the thermal and extended maser emissions, each emitting structure shows comparable fractions of the total emission due to the different solid angles subtended. As a result of this, most of the total maser emission comes from the extended maser emission rather than from the clumps and the arc. On the other hand, the maser emission of component A accounts for more than 80% of its total emission (thermal+maser).

Component B is mostly composed of two clumps spanning ≃5 km s⁻¹ and located at 5.7 and 11 ${R}_{\star }$ from the star along the directions with P.A. ≃ 65° and 235°, respectively. The total flux distribution per component is similar to that for component A, although the fraction of the thermal emission is higher.

Component C shows a very simple structure comprising two weak clumps located at 5.7 ${R}_{\star }$ with P.A. ≃ 180° and 295°. Most of its maser emission comes from the extended maser emission structure. This component composes the part of the brightness distribution in which the clumps account for the smallest fraction of the total emission.

Regarding component D, it can be explained with a simpler structure than component A. It also comprises different clumps distributed across Region II, but no arc is needed to explain the redshifted emission. The strongest clumps are located to the W and NW of the star. The clump to the W is particularly interesting because it spans eight velocity channels (≃7.4 km s⁻¹), although the exponential growth along the line of sight is also disabled, as in component A. There are a group of clumps that resemble an arc at 10.1 km s⁻¹ in the same position in the plane of the sky as the arc of component A. The fraction of maser emission accounted for by this maser structure is ≃35%, similar to that for components A and B. However, contrary to the rest of the components, the extended maser emission is less important, and less than one-third of the total emission is thermal.

Putting all these data together, 25% of the total SiS (14–13) line emission coming from the envelope covered by our CARMA observation has a thermal origin. The remaining 75% is maser emission: 60% comes from the extended maser emission, and 40% from the compact structures. The low fraction of thermal emission is directly related to the filtering of large scales performed by the interferometer, whereas the emission from compact structures is fully recovered. Thus, most of the emission of this line observed with single-dish telescopes is thermal.

The optical depth of every contribution to the synthetic emission depends on the physical conditions prevailing throughout the envelope and on the solid angle subtended by them, which is not known in our case. Hence, even though the model reproduces the observations reasonably well (see Figures 2, 3, 8, and 9), we cannot provide meaningful sizes and opacities for the clumps. However, in our model, the optical depth adopts different values varying from nearly 0 to a maximum of ≃27.8 very close to the star for the thermally excited component of the envelope. In the case of the maser-emitting regions, the optical depth ranges from negative values close to zero, related to very negative excitation temperatures (several times −100 K) or low population inversions, to ≃−11.3 at about 7 ${R}_{\star }$ from the star, which is associated with gas volumes with strongly inverted populations or low negative excitation temperatures (only a few times −10 K).

Modeling the most redshifted channels required an increment of the width of the outer acceleration shell at ≃20 ${R}_{\star }$. Initially, we adopted a shell width of ≃0.5 ${R}_{\star }$, but the observation suggests that its width is ≃2.5 ${R}_{\star }$. This modified acceleration shell has a very limited effect on models derived from previous single-dish observations, so it could be easily overlooked (e.g., Fonfría et al. 2008, 2015, 2017; Agúndez et al. 2012).

The interferometer data sets taken with CARMA and ALMA were complemented neither with on-the-fly maps acquired with single-dish telescopes nor with total power observations. The lack of short-baseline visibilities prevented us from properly covering the central region of the uv plane, filtering a significant amount of flux. The consequence of this flux loss is twofold: (1) the observed absolute flux is a lower limit to the total absolute flux, and it is not possible to compare it with the flux determined from lower angular resolution or single-dish observations (e.g., our IRAM 30 m data), and (2) the outer parts of the observed emission distribution (corresponding in our case to the lowest-level contours) could be affected by significant sidelobes in the PSF. However, we can properly deal with these effects by applying the same constraints to the models as the CARMA array applied on the actual brightness distribution, filtering the same flux as during the observations (for a deeper explanation, see Fonfría et al. 2014). This results in realistic predictions not affected by observational issues.

The SiS (14–13) data presented in this work were observed at the same time with a lower spectral resolution SiS (14–13) data set (≃12.3 km s⁻¹; Fonfría et al. 2014). A model composed of several compact emitting regions was proposed to explain the low spectral resolution emission. The lack of spectral resolution affecting these observations hampered our effort to get accurate positions of the maser-emitting regions and their physical properties. We can now compare models from both high and low spectral resolution observations. This comparison suggests that the main emitting components derived in the current paper were already proposed to exist by Fonfría et al. (2014). Our current model significantly refines the position and size of the components of the maser-emitting structure and the gas excitation conditions, providing us, in addition, with information about the spatial distribution along the line of sight, something impossible to determine from the low spectral resolution data set.

Since a significant part of the SiS (14–13) emission has maser nature, one important question is how the pumping mechanism works in IRC+10216. Fonfría Expósito et al. (2006) suggested that the overlap between lines of SiS and the abundant species C₂H₂ or HCN in the MIR could produce the population inversion involved in the maser emission. In particular, they proposed that the SiS rovibrational line 1–0R(13) could be blended under certain conditions with the C₂H₂ line ν₅R_e(9), with rest frequency 5.9 km s⁻¹ blueshifted with respect to the SiS line frequency. This blending would pump SiS molecules initially in the rotational level J = 13 of the vibrational ground state to the rotational level J = 14 in the vibrational state v = 1. The excited SiS molecules would de-excite again to the vibrational ground state populating the rotational level J = 15, since the SiS rovibrational transitions must fulfill the selection rule ΔJ = ± 1. This mechanism would invert the populations of levels J = 13 and 15, explaining the maser emission found in the SiS pure rotational lines J = 15–14 and 14–13. The small frequency shift between the SiS and the C₂H₂ rovibrational lines can be shortened by means of the Doppler effect if the C₂H₂ line is produced in a different part of the envelope than the SiS line. This indeed occurs in a radially expanding envelope, since every gas volume moves away from any other.

5.3.1. Extended Maser Emission

5.3.2. Compact Maser-emitting Regions

Although the presence of the extended emission can be assumed as a direct consequence of the pumping mechanism working in a radially expanding gas envelope, it does not explain so clearly the presence of the derived compact maser structures (the clumps and the arc, which can also be understood as a group of very close clumps). These clumps are probably related to small regions of the envelope where the pumping mechanism is the same as in the structure that produces the extended maser emission. However, the maser amplification in them is enhanced compared to their surroundings. This can be an effect of a local variation of the gas density/SiS abundance with respect to H₂ or of the excitation conditions/overlapping degree between the rovibrational lines involved in the pumping mechanism. It is still unclear what could force these physical parameters to vary, but we venture that they might be a consequence of (1) the increase of the gas density due to spatial irregularities in the mass-loss rate or the presence of arcs, (2) the increase of the molecular abundances caused by a variation in the gas-phase chemistry related to, e.g., an increment in the dissociating radiation field coming from outside the envelope that could reach its inner shells, or (3) inhomogeneities in the dust grain distribution that reduce the dust opacity, favoring the increase of the pumping radiation in the maser-emitting region. These scenarios could be compatible with the complex structure composed of thin arcs found in the envelope of IRC+10216 (Mauron & Huggins 1999, 2000; Leão et al. 2006; Cernicharo et al. 2015; Decin et al. 2015; Quintana-Lacaci et al. 2016; Guélin et al. 2018). These structures are predicted to appear naturally in envelopes surrounding single stars by Woitke (2006) and in binary systems by Mastrodemos & Morris (1999), and/or with the peculiar formation of molecules such as H₂O, C₂H₄, C₄H₂, and carbon chains close to the star (Agúndez et al. 2010, 2017; Fonfría et al. 2017, 2018).

Our model suggests that the only maser emission coming from clumps with impact parameters below $5\,{R}_{\star }$ in every velocity channel map is actually produced in Region II. In fact, there is no need to include clumps in Region I ($r\lesssim 5\,{R}_{\star }$) to reproduce the observations (Figure 10). The arguments discussed in the previous paragraph can be used to explain the existence of the clumps in Region II, but they are also compatible with the existence of clumps closer to the star, where the line width is high, favoring the radiative amplification due to the overlaps expected to exist between lines formed in different emitting volumes. This apparent incompatibility disappears if we realize that the pumping mechanism invoked to explain the maser emission of the SiS (14–13) line is a radiative process that can be collisionally canceled in regions with a high kinetic temperature and density, as in Region I. In order to prove it, we have solved the statistical equilibrium equations for SiS at different distances from the star assuming the physical conditions in IRC+10216 derived by Fonfría et al. (2015). The intensity of the pumping radiation associated with the overlap between the SiS line 1–0R(13) and the C₂H₂ line ν₅ R_e(9) was modeled as an enhancement of the continuum intensity that connects the SiS levels v = 0 J = 13 with v = 1 J = 14, keeping untouched the intensity at the frequency of any other line. The results of this analysis are graphically included in Figure 11, where we show the populations of some rotational levels of the vibrational ground state calculated with and without the continuum intensity enhancement. We conclude from these results that the rotational excitation of SiS can be considered to be dominated by collisions inward from ≃3 ${R}_{\star }$, disabling the radiative population inversion mechanism and preventing the maser radiation emission. This effect explains the hollowness of the extended maser emission inferred from our analysis of the CARMA observations. From 3 to 5 ${R}_{\star }$ the populations depart from thermal equilibrium, as the collisions cannot counteract the effect of the pumping radiation, getting inverted eventually for the levels J = 13, 14, and 15. The position in the envelope where the population inversion is reached depends on the intensity of the pumping radiation: the stronger this intensity, the closer to the star the population inversion occurs.

Zoom In Zoom Out Reset image size

Figure 11 shows that it is in principle possible to find maser emission beyond ≃20 ${R}_{\star }$, in disagreement with the results derived from the modeling of our observations. Moreover, the population inversion in Region III might even be boosted by the decrease of the gas density owing to the gas acceleration occurring around $20\,{R}_{\star }$, which weakens the effect of collisions on the level populations. However, the incompatibility between theory and observations is actually apparent, and it can be explained by the gas kinematic properties around ≃20 ${R}_{\star }$. The turbulent velocity in this region of the envelope is as low as ≃1 km s⁻¹, and the expected acceleration of ≃3–4 km s⁻¹ (Huggins & Healy 1986; Keady et al. 1988; Skinner et al. 1999; Fonfría et al. 2008, 2015; Agúndez et al. 2012) could disable the pumping mechanism that we propose to invert the populations, which is very sensitive to frequency shifts. In this situation, the overlap in the MIR between the C₂H₂ line coming from Region II, where the bulk of the pumping radiation comes from, and the SiS line in Region III is significantly reduced, strongly hampering the population inversion of the SiS rotational levels. Additionally, we also need to have in mind that the fact that the populations of two levels are inverted does not necessarily mean that we could observe maser emission. The population inversion has to be strong enough to produce noticeable effects. Otherwise, its effects could be negligible or be confused with suprathermal emission.

As we noted in Section 3.2, not all the SiS emission clearly follows the stellar pulsation. There are shorter-scale variations in several lines (e.g., J = 5–4 and 6–5) that in principle cannot be satisfactorily explained by this mechanism. The radiation field coming from the star and the dusty component of the envelope is expected to be driven by the stellar pulsation that varies periodically with time. Consequently, the effect of this emission on the molecules in the circumstellar envelope has to be also periodical. However, short-timescale variations of the radiation field due to, for instance, changes of the stellar pulsation, which is not very well known so far, could leave noticeable tracks on the molecular emission. Moreover, we note that the envelope of IRC+10216 is a very complicated environment that contains a large number of molecules affecting each other, and unexpected effects could arise if we consider a realistic time-dependent envelope model.

The comparison of the SiS (14–13) and SiS (15–14) lines acquired with CARMA (2011 January), ALMA (2012 April), and the IRAM 30 m telescope (2004 June and monitoring program, 2015 to present) reveals that their maser structure comprises two substructures that have different time dependences. Here we assume that the maser emission of both SiS lines is produced by the same structure because the pumping mechanism involves both lines. On one hand, the absence/presence of the most blue- and redshifted maser components (A* and D*) when lines J = 14–13 (CARMA) and 15–14 (ALMA) are compared (≃1.3 yr lapse between both observations) indicates significant changes with time (see Figure 1). This is supported by the fact that these maser components are also present in the IRAM 30 m telescope observations of both lines taken by us in 2004 and in the period from 2015 to present, as well as in the single-dish monitoring of the SiS (14–13) line recently carried out by He et al. (2017). The time monitorings show the strongest maser peaks appearing and disappearing, while the rest of the line profiles stay reasonably constant (Figure 6). Hence, the formation of components A* and D* depends strongly on the pulsation phase. On the other hand, components A, B, C, and D vary much less with respect to the bulk of the line, indicating that they are mildly affected by changes in the stellar emission or the envelope evolution.

One of the most stunning results of our analysis is the existence of the phase difference of ≃0.17 (≃110 days) between the NIR light curve and that of component A* with respect to the J = 14–13 and 15–14 SiS lines (Figure 7). The obvious similarity between this dependence and the NIR light curve denotes that this component is driven by the radiation coming from the central star. The time variation of other molecular species confirms the existence of this phase difference (Pardo et al. 2018).

The observations carried out by Shenavrin et al. (2011) and used by Menten et al. (2012) were taken in bands J, H, K, L, and M in the period from 1984 to 2008, covering ≃5.5 pulsation cycles. The fit to the observed data is reasonably good, and although small departures can be noticed in a couple of cycles, they could only imply phase differences of a few percent with respect to the fitted light curve. However, IRC+10216 can only be observed during part of the year, preventing us from tracing properly the light curve (Witteborn et al. 1980; Le Bertre 1988, 1992; Dyck et al. 1991; Monnier et al. 1998; Shenavrin et al. 2011). Consequently, we know neither how the actual light curve is nor whether this shape is constant over time. Thus, the assumption of a sine-like dependence on time imposes strong constraints on the reference Julian Date (usually at maximum light), the pulsation period, and the amplitude of the brightness variation, which perhaps do not represent accurately the real pulsation process. In fact, significant variations of the pulsation period over time, usually accepted as produced by thermal pulses, have been reported for several Mira stars (Templeton et al. 2005, and references therein). Among them, T UMi is one of the best examples of a Mira star, showing such a long-term variation with a fast decrease of ≃10 days every three to four pulsation cycles. Other stars could even display sudden variations on the pulsation period superimposed on their long-term behavior. Hence, it is possible that the light curve has changed from the end of the data acquisition to date.

Nevertheless, we have to consider the possibility that the time dependence molecular excitation can produce unexpected delays in the observed molecular emission. The radiative and collisional transitions happen in timescales shorter than 20 days for most of them (${A}_{v^{\prime} \to v^{\prime\prime} ,J^{\prime} \to J^{\prime\prime} }\gtrsim 6\times {10}^{-7}$ s⁻¹ for all the lines but the rotational J = 1−0; Cernicharo 2012), and they are usually much faster. Only the timescales related to the radiative transitions J = 1–0 in every vibrational state are longer than the detected phase difference (${A}_{v\to v,J=1\to 0}\,\simeq 7\times {10}^{-8}$ s⁻¹; Cernicharo 2012). Hence, it is unlikely that this offset is a consequence of the time evolution of the populations of the SiS rovibrational levels. However, a time evolution excitation effect cannot be completely ruled out if we consider the whole envelope in the problem. In particular, IR light could take as long as ≃10–20 days to reach the external layers where Pardo et al. (2018) have observed the spectra of several molecules showing a phase lag similar to what we have found for SiS. An accurate description of the molecular emission dependence with time ought to rely on a complex nonlocal, time-dependent excitation model, something that is beyond the scope of the current work but that we will address in the future.