© The Authors 2023

1 Introduction

When low-to-intermediate-mass stars reach the Asymptotic Giant Branch (AGB) they are characterized by an inert C and O nucleus surrounded by He and H burning shells (Höfner 2008). The early AGB phase begins with He-shell burning and subsequently evolves to a thermal pulse AGB phase where the H-shell burns. At the beginning of this phase, the atmosphere of the star is dominated by O-rich molecules (such as CO, TiO, SiO, H₂O, and VO; Wittkowski & Paladini 2014), but eventually (and depending on the mass of the star) the third dredge-up takes place, transporting the C formed in the He burning shell to the atmosphere. Depending on the efficiency of this dredge-up, the star will be of M-type (C/O ~ 0.5), C-type (C/O> 1), or an intermediate star (0.5 < C/O < 1; Höfner & Olofsson 2018).

Models of dust formation in AGBs suggest that pulsations are not sufficient to form the observed outflows and mass-loss rates. However, it is believed that there is a momentum transfer due to propagating shock-waves in the atmosphere that lift gas creating density-enhanced layers around the star where the nucleation of dust can occur (Jeong et al. 2003).

The models suggest that the mechanisms for driving dust differ depending on the dust composition. C-type AGBs form amorphous carbon grains with large opacities that can be accelerated by radiation pressure (Gautschy-Loidl et al. 2004; Sacuto et al. 2011; Rau et al. 2015, 2017). For M-type AGB stars, the high chemical complexity around them makes identifying the dust grains driving the dusty winds more complex. Theoretical models of the formation of dust-driven winds in M-type AGBs suggest that the distance(s) from the photosphere to the innermost molecular layers (located at a few stellar radii) determine how far the shock waves can levitate the gas. This is strongly coupled with the condensation region that the dust of a given chemistry must be formed to produce the observed winds (Nowotny et al. 2005; Bladh & Höfner 2012; Liljegren et al. 2016). Therefore, high-angular-resolution observations of the innermost molecular layers are particularly important to constrain the gas-dust coupling and the dust chemistry (see, e.g., Norris et al. 2012).

During the last decade, due to the enhanced spectral and spatial resolution made available, near-infrared interferometry has become an important technique to map the innermost structure of AGBs’ wind in order to link the observed morphology with the mass-loss models. Most of the analyses based on interferometric data consist of model fitting the observables directly (squared visibilities -V²- and closure phases -CPs-) with parametric models of the star’s atmosphere and/or more sophisticated radiative transfer and hydrodynamic simulations with codes like CODEX (Ireland et al. 2008, 2011) or CO5BOLD (Freytag et al. 2002, 2003). Some examples of these studies include Ireland et al. (2004), Wittkowski et al. (2008, 2011, 2016), and Haubois et al. (2015). However, due to the complexity of the environment, it is still difficult to fully reproduce the asymmetries traced by the CPs of the interferometric data with those models. Therefore, complementary techniques, such as aperture-synthesis image reconstruction can be used to better depict the inner circumstellar layers of AGB stars. Examples of interferometric imaging applied to AGB stars include the analysis presented by Wittkowski et al. (2017) on the C-rich AGB star R Sculptoris using the H–band (λ₀ ~ 1.76 μm) beam combiner PIONIER (Le Bouquin et al. 2011) at the Very Large Telescope Interferometer (VLTI). In that study, the reconstructed images are consistent with a dynamic atmosphere and a stellar disk dominated by giant convection cells, which has a strong impact on the distribution of the molecular layers and the formation of dust at a few stellar radii.

A further example includes interferometric imaging applied to the long period variable (LPV) Mira using the Interferometric Optical Telescope Array (IOTA) in the H-band, as reported by Perrin et al. (2020). In that study, a large localized absorbing patch is observed in the images. This patch is interpreted as a local dust formation zone, and it is consistent with the combined effect of pulsation and convection in the mass-loss process of AGBs (see also Paladini et al. 2018). These examples highlight the importance of interferometric imaging to confront the observed morphologies with theoretical models of dust-driven winds and mass-loss processes in AGBs. However, there is still limited observational information, in particular for the case of M-type AGBs.

In this work, we aim to characterize the innermost environment of the M-type Mira star R Car using infrared interferometry in order to link its morphology with models of the formation of dust-driven winds. For that purpose, we present new K-band (λ₀ = 2.2 μm) interferometric data obtained in 2018 with the instrument GRAVITY (Abuter et al. 2017) at the VLTI. These data were analyzed to depict the innermost structure of R Car at two phases within one pulsation cycle. Our analyses include: (i) parametric model fitting of the squared visibilities (V²); (ii) image reconstruction of both the pseudo-continuum and across the first (λ₀ ~ 2.29 μm) and second (λ₀ ~ 2.32 μm) CO band heads; and (iii) an estimation of the CO layer(s) physical parameters such as temperature, size, and optical thickness by using a single-layer model. This paper is structured as follows: observations and data reduction are presented in Sect. 2; in Sect. 3, we describe the results of our geometrical model fitting, our image reconstruction procedure, and our single-layer model. Our results are discussed in Sect. 4. Finally, we present our conclusions in Sect. 5.

2 Observations and data reduction

2.1 Properties of R Car

Miras are AGB stars with large-amplitude variations over long periods (10² days), for which asymmetries in the circumstellar material are quite frequent (Cruzalèbes et al. 2015; Ragland et al. 2006; Lacour et al. 2009). R Car is a bright, near-infrared Mira star (K–band: −1.23 mag) located at 182 ± 16 pc from Earth (Gaia Collaboration 2021). This source exhibits variations in the V band from +3.9 to +10.5mag. Figure 1 shows a visual light curve of R Car based on data obtained from the American Association of Variable Star Observers (AAVSO) database¹.

The photosphere of the star was observed in 2014 as part of the Interferometric Imaging Contest (Monnier et al. 2014) using the H-band four-beam combiner interferometer PIONIER attached to the VLTI. The reconstructed image using these data revealed bright spots on its surface and (potential) inner layers of circumstellar material and a stellar disk with a radius of 5 mas.

2.2 GRAVITY observations and interferometric data reduction

R Car was observed in January and February, 2018 for a total of eight hours, with the beam-combiner GRAVITY. The data were collected using six baselines with the Auxiliary Telescopes (ATs) configuration A0-B2-C1-D0, which provides minimum and maximum baselines of 7.67 m (equivalent to an angular resolution of ∆θ ~ 59 mas in this wavelength range) and 32 m (∆θ ~ 15mas), respectively. The data cover a wavelength range between 2.0 μm and 2.4 μm with a spectral resolution of R ~ 4000. The complete log of observations is reported in Table A.1.

We reduced and calibrated our data using the latest release of the GRAVITY pipeline (version: 1.5.0). Detailed descriptions of the recipes are available in the ESO pipeline manual² and in (Lapeyrere et al. 2014). The same standard data reduction procedures were applied to the scientific target and the calibrator star (HD 80404; which has an angular size of 1.8 mas) to correct for the DARK, FLAT, and BAD calibrations, as well as to extract the interferometric observables. The calibrated V² of the target were obtained following the standard procedure described in the ESO Gravity pipeline manual, which consists of dividing the raw observables by the ones of the calibrator star. Calibrated CPs were obtained by subtracting the raw CPs of the calibrator from the ones of the source. Due to the observed changes in the observables at similar baselines, the data were split into two epochs (January and February, with equivalent pulsation phases of 0.56 and 0.66, respectively). The calibrated R Car data files can be downloaded directly from the Optical Interferometric Data Base (OiDB) archive hosted at the Jean-Marie Mariotti Center³.

Fig. 1 Visual light curve of R Car based on AAVSO database (blue points). The red line represents the best fit of a sine function to the curve which results in a period of P = 314 ± 1 days. The dashed vertical lines and the colored diamonds show the epochs of our 2018 GRAVITY observations (see the labels on the plot).

2.3 GRAVITY spectrum calibration

GRAVITY data provide us with K–band spectra (R ~ 4000) of the target and calibrator. For the posterior analysis of the CO band heads observed in the R Car spectrum, we calibrated and normalized it to remove the telluric contribution from the Earth’s atmosphere using the following procedure. First, we averaged the R Car and calibrator spectra to obtain the high signal-to-noise spectra and , respectively. To remove the telluric lines from the averaged R Car spectrum, we used the following expression:

(1)

where was calculated by averaging 15 theoretical spectra from the theoretical spectra web server⁴ (Allard et al. 2011) of our calibrator spectral type A7 Ib. Second, we flattened and normalized (in the range of 2.21–2.34 μm) using a third-order polynomial function. Finally, we applied a wavelength correction (equivalent to a Doppler velocity of 17 km s⁻¹) for the local stan-dardofrest(LSR), taking into accountthe propermotion and the distance of the source. The normalized R Car spectrum of each epoch is presented in Fig. 2.

Fig. 2 R Car calibrated spectra and V2. Top panel shows the normalized calibrated spectra of the GRAVITY data. The blue spectrum corresponds to the January epoch, while the red one corresponds to the February one (see the labels on the plot). The middle and bottom panels show plots of the V2 versus wavelength for the two epochs observed.

3 Analysis and results

3.1 Geometrical model

To obtain the angular size of R Car across the K band, we applied a geometrical model of a uniform disk (UD) to the V² data. The visibility function of this model is given by the following equation (Berger & Segransan 2007):

(2)

where , u and υ are the spatial frequencies sampled by the interferometric observations, J₁ is the first-order Bessel function, and Θ_UD and F_r are the angular diameter of the uniform disk profile and a scaling factor that accounts for the over-resolved flux in the observations, respectively. Prior to modeling the pseudo-continuum data, we averaged 200 spectral bins into one single measurement. This allowed us to obtain eight spectral channels across the 2.0–2.4 μm range. Considering the calibration error in the averaged quantities, the final precision of the V² is and the precision of the CPS σ^CP is ~ 1 deg. Figure A.1 shows the calibrated V² and CPs as well as the u−υ coverage of the two epochs of observation for the eight spectral channels probed.

To fit the data, we used a Markov chain Monte Carlo (MCMC) algorithm based on the Python package emcee (Foreman-Mackey et al. 2013). We let 250 walkers evolve for 150 steps using the data of each spectral bin independently. Table 1 shows the best fit of Θ_UD and F_r for each one of the observed epochs. Figure 3 shows the V² data and their corresponding best-fit models. Figure 4 shows the evolution of the best-fit Θ_UD versus wavelength of the two epochs of observation. From this figure, we observe two significant changes in the pseudo-continuum structure of R Car.

First, the size of the target changes with wavelength, decreasing as we move toward the center of the bandpass (λ ~ 2.23 μm) and then increasing again. This trend is expected due to the presence of molecules such as H₂O at ~ 2 μm and CO at ~ 2.29–2.33 μm lying above the stellar surface (see, e.g., Wittkowski et al. 2008). The smallest Θ_UD value (at ~2.23 μm) corresponds to the closest pseudo-continuum-emitting region due to the photosphere. This size is the best estimate of the stellar disk diameter in the K band. Ireland et al. (2004) derived a wavelength-dependent Gaussian full width at half maximum (FWHM) between 15 and 20 mas at 700–920 nm. Our best-fit UD diameter at ~2.23 μm indicates that the star appears to be larger near the V band than in the K band. This size difference is due to the presence of dust layers obscuring the stellar-disk in the V band, but that are transparent in the K band, allowing us to have a closer look at the continuum-emitting region.

Second, the disk diameter is larger during the January epoch by up to 12% compared with the February epoch. From the light curve reported in Fig. 1, it can be observed that the magnitude of the star increased in the February epoch, when the diameter of the star is smaller. This behavior comes with an expected increment of the temperature during February, which makes the star look brighter at the cost of it having a smaller radius (see, e.g., Wittkowski et al. 2012, 2018; Haubois et al. 2015). A rough estimate of the temperature difference (T_J/T_F) between the observed epochs could be obtained using the Steffan–Boltzmann law and the magnitude difference (|∆V_mag| = 1.8) from Fig. 1, resulting in T_January/T_February ~ 0.62.

Fig. 3 V2 of the pseudo-continuum versus spatial frequencies of both epochs. The lines indicate the best-fit UD models. Each color corresponds to a different wavelength (see the label on the plot). We highlight that the spatial frequencies increase in opposite directions depending on the epoch. This allowed us to have a better visual comparison of the trends present in the V2 data between the two epochs.

Fig. 4 Best-fit ΘUD as a function of wavelength for the two analyzed epochs. Each color corresponds to a different epoch (see the labels on the plot).

3.2 Spherically thin layer: MOLsphere model

The GRAVITY spectra show prominent CO band heads (see Fig. 2). The data show drops in the V² of the CO 2-0 (λ₀ ~2.2946 μm) and CO 3-1 (λ₀ ~2.3240 μm) vibro-rotational transitions compared with the pseudo-continuum, which indicates that the CO line-emitting region is more extended (see, e.g., Wittkowski et al. 2012). The innermost CO layers are important to trace the dust nucleation region since most of the carbon or oxygen present in the atmosphere is concentrated in the stable bound of this molecule (Höfner & Olofsson 2018).

To characterize the regions where those transitions originate, we employed a single-layer model (called MOLsphere) that has proven to be successful in explaining the absorption of the CO molecule that appears in the visibilities of evolved stars (Perrin et al. 2004, 2005; Montargès et al. 2014; Rodríguez-Coira et al. 2021). This model consists of a stellar disk with a compact layer around it. The star is modeled by a stellar surface of radius R_* that emits as a black body at a temperature of T_*. It is surrounded by a compact spherical layer with a radius of R_L that absorbs the radiation emitted by the star and re-emits it like a black body. The MOLsphere is characterized by its temperature, T_L, radius, R_L, and its optical depth, τ_λ. The region between the stellar photosphere and the layer is assumed to be empty (see Fig. 5 for a sketch of the model). The analytical expression of the model is given by

(3)

where ; T_* and T_L are the temperatures of the photosphere and of the CO layer, respectively; R_* and R_L are the angular radius of the star and the layer, respectively; τ_λ is the optical depth of the molecular layer at wavelength λ; B_λ(T) is the Planck function (at wavelength λ and temperature T); and β is the angle between the radius vector and the line-of-sight so that .

For the fitting, we adopted a T_* = 2800 K and R_* = 5 mas (Monnier et al. 2014; McDonald et al. 2012). We used the disk estimation reported in the H band because it is considerably less affected by molecular contribution, in contrast with the K band where molecules such as CO, H₂O, and OH, among others (see, e.g., Wittkowski et al. 2018; Paladini 2011), are present over the entire band. The unknown parameters are then T_L, R_L, and τ_L. We performed a two-step procedure to estimate these parameters. As a first step, for each pair of T_L and R_L, we performed a least-squares minimization to find the best-fit value of τ_L that reproduces the spectrum . The next step consisted of using an MCMC method based on the Python library emcee to estimate the best combinations of T_L, R_L, and τ_L that reproduce the V².

To account for the over-resolved flux when modeling the V² data, we added a scaling factor of Fr. This value was estimated simultaneously with the other parameters in the model. Table 2 shows the corresponding results for the good quality spectral channels across the two CO band heads. Figure 6 shows the R Car’s spectra with the corresponding flux values obtained with the single-layer model. Figures 7 and 8 show the best fit to the V² with the parameters obtained from the model. From the bestfit models, we derive a temperature range between 1300 and 1500 K, which is consistent with the expected temperature values to produce the CO vibro-rotational transitions observed in the K band (Geballe et al. 2007). The values of the optical depth, τ_L, support that the layer is optically thick across all the spectral channels, with higher τ_L toward the centers of the band heads. In Sect. 4, we discuss the implications of these results in the context of the wind-driving dust candidates that can coexist with the CO.

Fig. 5 Illustration of single-layer model. The yellow disk represents the star, and the orange ring represents the layer. The parameters in the image are as described in Sect. 3.2. This image was adapted from Rodríguez-Coira et al. (2021).

3.3 Image reconstruction

3.3.1 Regularized reconstructed images

The geometrical model described in Sect. 3.1 and the MOL-sphere model in Sect. 3.2 show that the source is symmetric. However, the CPs observed in our data differ considerably from zero or 180 degrees, which indicates that the source is asymmetric (see Fig. A.1). So, in order to depict the complex nature of the source across the K band and at the CO band heads, we reconstructed aperture-synthesis images using the software BSMEM (Baron & Young 2008). This code uses an iterative regularized minimization algorithm of the following form:

(4)

where f_data(x) is a measurement of the log-likelihood and it characterizes the discrepancy between the image model (at a given iteration) and the observables; f_prior(x) contains the prior information included in the reconstruction, and the hyperparameter μ is used to adjust the relative weight of the constraints between the measurements and the ones set by the priors (Renard et al. 2011). BSMEM uses Entropy (Gull & Skilling 1984) as regularizer in f_prior (x).

As initial setup for our reconstruction, we created images of 128 × 128 pixels with a pixel scale of 0.469 mas. The best-fit sizes of the UD models were used as starting point in the reconstruction process. We performed the reconstructions by taking all the V² and CPs on each individual channel separately. Since our data are constrained by a relatively sparse u−υ coverage, this could affect the quality of the reconstructed images. Therefore, to properly distinguish intrinsic astrophysical features from image reconstruction artifacts (Haubois et al. 2015; Sanchez-Bermudez et al. 2018) and to provide an estimate on the quality of our reconstructed images, we employed a bootstrapping method (Zoubir & Iskander 2004).

We built new data samples from the original one. The samples were created by keeping the original number of data points unaltered, but allowing random sampling with replacement. This created data sets with points with higher weights than in the original data set and some others with zero weights. This algorithm also produces slight changes in the u−υ plane, which allows us to trace the impact of the u−υ coverage on the reconstructed images, without loosing the statistical significance of the original number of data points. According to Babu & Singh (1983), the statistical moments such as mean, variance, and standard deviation of the bootstrapped samples are good approximations to the statistical moments of the original data sets.

We employed this method by creating 50 different sampled data sets from the original observables, and CP_* = CP₁,…, CP₅₀, using the Python scikit learn package (Pedregosa et al. 2011). We performed image reconstructions on each sample to obtain a set of bootstrapped images: I_* = I₁…, I₅₀. All the reconstructions from the samples converged, and the best-fit image reproduces the trend in the sampled observables (V² and CPs) quite well. We calculated the mean and standard deviation per pixel along the 50 bootstrapped images to obtain a mean image I_mean and a standard deviation map I_σ. We used I_mean/I_σ to estimate the signal-to-noise ratio (S/N) of the structures present in the mean image. We discarded all the features with an S/N < 3 as they are related to reconstruction artifacts. Figures 9–11 show the I_mean images for the pseudo-continuum and CO band head channels recovered. All the observed structures in the I_mean images lie within the three S/N contours. Figures B.1–B.3 show the best-fit synthetic observables obtained from the I_mean images to the original data points. As can be seen, the trends are well-reproduced at all angular scales and wavelengths.

From a visual inspection of our best reconstructed images, I_mean (Figs. 9–11), we note that all the pseudo-continuum images are asymmetric and show structural differences between the observed epochs. In our January images, we observe a spot toward the south-east; while in February, a more elongated (and centered) structure is observed. Except the change in size across the band pass, we do not see significant structural changes with wavelength.

Fig. 6 Best fit to spectra with the single-layer model for the two CO band heads and epochs. The gray dots correspond to the data points and the solid-lines correspond to the best-fit model (see the labels on the plot).

3.3.2 Principal component analysis to understand the structure of R Car

To have a more precise characterization of the asymmetries in the reconstructed images of the CO band heads, we used the principal component analysis (PCA) described by Medeiros et al. (2018). Those authors demonstrated that the visibilities of the principal components are equal to those of the visibilities. This method is useful for tracing the changes across a set of images that have the greatest effect on the visibility (amplitude and closure phase) profile. In our case, we estimate the most significant structural changes of the observed asymmetric structures across wavelength for each of the observed epochs. The following procedure was applied to the ensemble of wavelength-dependent images per epoch to extract their principal components.

For the PCA analysis, we normalized each data set by subtracting the corresponding mean image and dividing by their standard deviation image across the wavelength range. To perform the PCA analysis, we used the CASSINI-PCA⁵ package. This software allowed us to compute the covariance matrix of the data set and to transform it into the space of the principal components. This allowed us to determine the eigenvectors (or eigenimages in our case) and their corresponding eigenvalues. Since we only have seven images per data set, and the possible number of components must be smaller than or equal to the number of images in the data set, we decided to only keep the first four principal components. From our tests, we observe that those components explain, at least 93% of the variance in the data sets.

The maps that correspond to the first four principal components of the CO band heads are reported in Fig. 12. For both CO data sets, we observe that the first two components trace the variance associated with the general changes in brightness and size observed in the images of the CO band heads across wavelengths. However, components 3–4 show considerably more asymmetric changes in the external structure of the target (for contours between 50% and 10% of the peak). The observed varying structures in components 3–4 are evenly distributed across the different position angles in the structure of the source. We suspect that the variance explained by these components is associated with the more asymmetric structure in the CO images. This could also be connected with the possibility of having a clumpy CO distribution of material in the envelope around R Car. As an example of the contribution of each component to the asymmetric structure of the target, in Fig. 13 we include the CPs obtained at a wavelength of 2.2946 μm from the recovered images with cumulative components between 1 and 4. As can be observed, the recovered images improve the recovery of the asymmetries traced by the CPs as we add more components to them. This is an indication that the variance observed in the eigenimages is tracing a more changing environment in the CO distribution across the sampled wavelengths. Unfortunately, the limited resolution of our observations is not enough to properly resolve these asymmetries completely.

Fig. 7 Best fit to V2 from the best single-layer models across the first CO band head for our two epochs. The colored dots and the shaded regions correspond to the best-fit models and their corresponding error-bars. Each color correspond to a different wavelength (see labels on the plot). The data are shown with black dots. We highlight that the spatial frequencies increase in opposite directions depending on the epoch. This allowed us to have a better visual comparison of the trends present in the V2 data between the two epochs. Also, for better visualization, we displaced the V2 vertically.

Fig. 8 Best fit to V2 from the best single-layer models across the second CO band head for our two epochs. The panels are as described in Fig. 7.

Fig. 9 Best pseudo-continuum reconstructed images for the (a) January and (b) February epochs. The white ellipse located in the right corner of the image at 2.0160 μm corresponds to the mean synthesized primary beam. The blue contours represent 10, 30, 50, 70 90, 95, 97, and 99% of the intensity’s peak.

Fig. 10 Best reconstructed images per spectral channel for the (a) January and (b) February epochs across the first CO band head. The white ellipse located at the right corner of the image at 2.2941 μm and the blue contours are as described in Fig. 9.

3.4 Pure-line CO images

Our CO reconstructed images contain a contribution associated with the underlying stellar component. To remove such contribution, we first centered the images across the CO band heads in a common reference position. For this purpose, we used the differential phases (DPs) in the GRAVITY data. DPs measure the photocenter displacement of an emission line with respect to the continuum for a given baseline length and orientation. We only restricted the pure-line CO images to the first band head, since the DPs of the second band head have much larger uncertainties. To compute the global 2D photocenter displacement vector at a given spectral channel, we used the following relation (Kraus 2012):

(5)

where ϕ_i is the GRAVITY differential phase measured on baseline i, is the corresponding baseline vector [B_x, B_y], and λ is the central wavelength of the working spectral channel. By solving the matrix represented by Eq. (5), we were able to find the global [p_x, p_y] displacements of the emission line. We solved the system of equations employing the least-squares algorithm implemented in the Python function numpy.linalg.lstsq. This function inverts the matrix formed by each one of the observed differential phases per baseline per channel to find the solution of the astrometric displacements. To obtain an estimation of the precision (i.e., 1σ error-bars) of the reported displacements, we applied the same matrix inversion process to 100 different samples obtained with random points using a Gaussian distribution with a standard deviation equivalent to the error bars of the differential phases. The final astrometric displacements are equivalent to the mean of the 100 different samples, while their error bars are equivalent to the standard deviation of the displacements from the samples. We obtained the displacements for the different spectral channels across the CO band heads, as well as for a few continuum channels. Figure 14 shows the first CO band head offsets for our two epochs. We can observe that the wavelengths corresponding to the continuum emission lie close to the reference position [0,0], and the ones of the band head show displacements within 1 mas. It is interesting to highlight that the CO centroid position of the two epochs is different, confirming the variability of the structure of the source between them. Once the centroid position of the images was computed, we used them to extract the pure-line CO images using the following method. First, CPs are shift-invariant quantities, it means that images reconstructed with this observable could be placed at different positions in the pixel grid and, still, reproduce the CPs in the data. Therefore, to account for this effect, we computed the continuum and CO images, and we performed a sub-pixel shift to the corresponding centroid position obtained from the astrometric displacements [p_x, p_y] of the differential phases. For the subsequent analysis, we took the continuum images at λ_cont = 2.228 μm. Second, once all the centroids of the images matched with the astrometric displacements, we subtracted the continuum image from each one of the CO images across the different spectral channels. Since the total flux in the reconstructed images is normalized to unity, before performing the subtraction we calibrated the total flux on each image by scaling them with their corresponding normalized flux scale from the spectra in Fig. 2.

Figure 15 shows the CO images with the continuum contribution subtracted. We masked the central region to better visualise the CO-emitting region. From a visual inspection of the images, we can identify some features. First, there are different structures between wavelengths and epochs, suggesting the inhomogeneous nature of CO. The images also show bright and faint regions at different position angles. We relate them with zones of lower (faint regions) and higher (bright regions) densities and temperatures.

Fig. 11 Best reconstructed images per spectral channel for the (a) January and (b) February epochs across the second CO band head. The white ellipse located in the right corner of the image at 2.3229 μm and the blue contours are as described in Fig. 9.

Fig. 12 First four principal components for the January and February epochs for the first CO (plots a and b) and the second CO band heads (plots e and f). The relative scale of the structures in the eigenimages is displayed with a color bar at the bottom of each panel. Only the central 20 mas of the eigenimages are shown in each panel. The white contours correspond to the mean image across wavelength (per given data set) and they represent 10, 30, 50, 70, 90, 95, 97, and 99% of the intensity’s peak.

Fig. 13 Representative example of CPs versus spatial frequencies of the data set that corresponds to the first CO band head at 2.2946 μm (gray squares). The colored dots show the CPs from the recovered images obtained with different numbers of principal components (see label on the plot).

Fig. 14 Two-dimensional photocenter offsets in mas. The black dots correspond to the photocenter position of the continuum emission and the labeled colored dots to the photo-centers of the CO band head. The error bars of the reported displacements are smaller than the used symbols.

4 Discussion

4.1 CO column densities from the single-layer model

The single-layer model described in Sect. 3.2 allowed us to obtain information about the CO layer’s physical conditions (e.g., opacity, temperature, and sizes). Using the optical thickness and temperature derived for the layer, we estimated the corresponding column density in units of cm⁻² using the following relation (Savage & Sembach 1991):

(6)

where τ_L is the optical thickness obtained from our single-layer model, ∆v is the width of the band head in km s⁻¹, f_ij is the absorption oscillator strength that depends on the temperature, λ_ij is the wavelength of the transition, c is the speed of light, e⁻ is the electron charge and m_e its mass, and the constant term π e²/m_e c² = 8.8523×10⁻¹³ cm mol⁻¹. We calculated the column density at the center of the first CO band head, where the optical thickness has its largest value, by setting λ_ij = 2.2946 μm. Then, using Eq. (3) from Goorvitch (1994), we calculated f_ij for the vibro-rotational transition V_0->2, J_20->21 (V and J are the vibrational and rotational states, respectively). We obtained their corresponding Einstein coefficients, A, from Chandra et al. (1996). Finally, we calculated ∆υ by fitting an inverse Lorentzian profile across the first CO band head for the two epochs. Since ∆υ is expressed in km s⁻¹, ∆υ = c∆λ/λ₀.

The obtained column densities for both epochs lie in the N_CO ~ 9.19×10¹⁸–1.01 × 10¹⁹ cm⁻² range. This result has two implications. On one side, these values are consistent with the range of column densities reported by authors including Rodríguez-Coira et al. (2021) and Tsuji (2006) for other evolved stars. Furthermore, these values, on average, are high enough to shield the dust particles allowing them to grow to micron-sized scales, where the stellar radiation could act efficiently to drive the observed winds in this kind of object. On the other side, our current single-layer model returns large uncertainties in the estimated parameters. Therefore, we cannot confirm if CO formation happened between both epochs.

To obtain a direct comparison between the size of the photosphere (R_*) and the size of the CO layer (R_L), we calculated the ratio R_L/R_*. We take the R* value reported by Monnier et al. (2014) of 5 mas as a reference. From R_L/R_*, we find a mean extension of the layer R_L ~ 2.4 R_*. To verify the consistency of our results, we compared them with other O-rich Miras for which the single-layer model has been applied. Table 3 shows the values of the layer temperature and R_L/R_* for the Mira stars R Leo, χ Cyg, and T Cep, as well as for our R Car estimates. From the values reported, we confirm that the parameters derived for R Car are consistent with other estimates reported in the literature for similar type of objects.

4.2 Observational restrictions to the wind driving candidates for M-type Mira stars

Given the high chemical complexity around M-type AGB stars, it is not easy to determine the chemical composition of dust-driven winds. Even though observations tell us about the existence of different dust species in the circumstellar environment (Molster et al. 2010), not all are efficient in absorbing or scattering the star’s radiation to trigger a stellar wind. A condition for a dust grain to be considered a wind driver is that it must be formed closer to the stellar surface within reach of shock waves produced by the star. The shortest distance from the star at which a dust grain can be formed depends on its chemical and optical properties (Bladh & Höfner 2012).

As mentioned in Sect. 1, the observed CO vibro-rotational transitions (with ∆V = 2) originate in the deep photospheric layers where dust is forming. Therefore, with the CO layer parameters, obtained from our single-layer model and the pure-line CO images, we can set constraints on the types of dust grains that are considered wind-driving candidates and that can coexist with similar physical conditions to the CO innermost layers detected with our GRAVITY observations. It is important to remember that for winds to be triggered, the star pulsations must levitate the gas to a certain distance (defined as levitation distance, R_ı) where dust can condense (known as condensation distance, R_c). Bladh & Höfner (2012) define the levitation distance with the following expression:

(7)

where u₀ is the initial velocity of the gas at a distance, R₀, and u_esc = (2GM_*/R_*)^1/2 is the escape velocity. The levitation radius can be extracted directly from the size of the layer from our single-layer model or from the pure-line CO-reconstructed images. The single-layer model provides us an average estimate of R_ı = 2.43 R_*, since the model only considers a symmetric structure. However, the reconstructed images provide us a range of R₁ between 2.25 and 2.72 R_* due to the observed asymmetries. Assuming u_esc = 44 km s⁻¹, M = 0.87 M_⊙ (Takeuti et al. 2013), T_eff = 2800 K (McDonald et al. 2012), and R_* = 167 R_⊙ (Monnier et al. 2014) for Eq. (7) and that the CO band heads that we observed are formed in a region from the stellar surface at R₀ = 2 R_*, we computed an initial velocity of the CO of u₀ = 13 km s⁻¹ considering the single-layer model and u₀ between 10 and 16 km s⁻¹ using the CO reconstructed images. These values are consistent with the ones reported in the literature on the order of u₀ ~ 10 km s⁻¹ (Nowotny et al. 2010).

To determine the dust condensation radius, R_c, (Eq. (8)) we have to consider a Planckian radiation field and a thermal equilibrium for the determination of the dust grain temperature:

(8)

where T_c is the dust’s condensation temperature, T_* is the star’s effective temperature, and p is the dust absorption coefficient. We estimated the condensation radii for several types of dust grains assuming the condensation coefficients and temperatures reported by Bladh & Höfner (2012) and references therein. These values are reported in a p versus T_c diagram. Considering that the levitation radius has to be larger or equal to the condensation radius, we obtained curves of possible p and T_c values with the R_l estimates from our single-layer model and pure-line CO reconstructed images.

Figures 16 and 17 show the T_c versus p, diagram. In both figures, the red dots correspond to a variety of dust grains with different condensation temperatures and values of p. The solidblack line in Fig. 16 displays the values of T_c and p for the R_l value obtained from the single-layer model. The colored solid lines in Fig. 17 correspond to the R_l values extracted from the pure-line CO images for both observed epochs at different orientations: 0°, 90°, 180°, and 270°. Since we have seven images per epoch, we adopted the median value per image data set at each orientation.

Those figures show that our CO layer radius is consistent with the condensation radius for composites Mg₂SiO₄ and MgSiO₃. However, we find layer temperatures around 1300 ± 100 K (blue and green shaded regions on Fig. 16), which are slightly higher than the expected condensation temperature for those two dust composites (T_c ~ 1100 K). This temperature discrepancy could be explained by the fact that the single-layer model considers a homogeneous and symmetric distribution of the CO layer. However, the reconstructed images clearly show that the CO distribution is neither homogeneous nor symmetric as shown in Fig. 17 where, depending on the position angle, either one or both composites lie(s) below the estimated levitation radius. Therefore, regions with lower temperatures and larger levitation radii could exist around R Car. This supports the hypothesis that dust formation could occur in an anisotropic way over the CO shell extension.

Magnesium-iron silicates (olivines and pyroxenes) with nano grain sizes (~0.1 μas) have already been detected in the circumstellar environments of M-type AGB stars (Norris et al. 2012; Ohnaka et al. 2017), and they appear to be the main thing responsible for the dust-driven winds observed in these stars. While our GRAVITY observations and models suggest the possible coexistence of those dust types at the inner gaseous (CO) shells, we require additional interferometric observations in the mid-infrared (potentially with the beam-combiner MATISSE-VLTI) to confirm this scenario. Very recent models of dusty winds in M-type AGBs suggest Fe enrichment of the wind as it moves far away from the photosphere (Höfner et al. 2022). Those models suggest an inner core of Mg₂SiO₄ (or of a similar composite) surrounded by more outer layers of MgFeSiO₄. New MATISSE-VLTI interferometric observations could help us to confirm this scenario in the proximity of R Car in the near future. By confirming this scenario, mid-infrared interferometric observations could shed new light on the problem of a lack of the 10 μm silicate feature in the spectra of M-type AGB stars.

Fig. 15 Pure-line CO images of first CO band head for the (a) January and (b) February epochs. The region that corresponds to the subtracted continuum contribution has been masked for better visualization of the CO region. The white ellipse located in the lower right corner of the image at 2.2941 μm corresponds to the mean synthesized primary beam.

Fig. 16 Curve of constant levitation radius (Rl = 2.43 R*) represented with a solid-black line. The plot shows the condensation temperature, Tc, as a function of the dust absorption coefficient, p. The red circles indicate the positions of a variety of dust grains. The blue and green shaded regions correspond to the CO temperature ranges estimated from the single-layer model (see labels on the plot).

Fig. 17 Panels similar to Fig. 16 but showing the median condensation radii obtained directly from the pure-line CO reconstructed images at different position angles (see labels on the panels) for (a) the January and (b) February epochs.

4.3 The role of CO asymmetries on mass-loss estimates

The mass-loss rate is a crucial parameter for understanding the evolution of M-type AGB stars. Bladh et al. (2019) presented an extensive grid of DARWIN dynamic atmosphere and wind models to establish the mass-loss rates, wind velocities, and grain sizes for M-type AGB stars. These models assume that the principal mechanism for wind triggering is the radiation scattering on Mg₂SiO₄ composites. This grid of models produce mass-loss rates on the order of 1 × 10⁻⁶ M_⊙ yr⁻¹ for stellar parameters similar to R Car (e.g., M_* =0.75 M_⊙; log[L_*/L_⊙] = 3.55; P = 283 days).

More importantly, those models show that, for the mass (0.75–1 M_⊙) and luminosity (log[L_*/L_⊙] = 3.7) ranges of R Car, large values of the pulsation amplitude in the models and high chemical enrichment develop stable winds based on the Mg₂SiO₄ composites. However, there is a region (when ∆u_p = 2 km s⁻¹ and n_d/n_H = 1×10¹⁵) of low chemical enrichment and small pulsation amplitudes where those DARWIN models do not develop stable dust-driven winds. While those models are important to demonstrate the dynamical evolution of dust-driven winds with chemical properties based on magnesium-iron silicates, they do not incorporate the importance of asymmetries on the development of the wind.

Liljegren et al. (2018) studied the effects of complex con-vective structures on wind production and mass loss in M-type AGB stars. The authors used the CO5BOLD code (Freytag et al. 2003) to simulate the convective surface and pulsations ofAGBs, using the results as boundary conditions in their DARWIN models. Those models show that gas is not levitated uniformly above the star’s surface (only about 70% or less of the stellar surface is covered), causing dust not to be formed in a spherical layer around the star. Instead, it forms clumps. This appears to have a direct implication in an overestimation of the mass-loss rate determined with homogeneous models. In this context, our work shows direct evidence of the asymmetric and variable gas distribution at the innermost layers of M-type AGB stars. A specific DARWIN model of R Car is beyond the scope of this paper and it is left for future studies. However, by contrasting the structures observed in our reconstructed images and future dedicated CO5BOLD/DARWIN models, we could obtain better constraints on the mass-loss evolution and the role of magnesium dust composites on the development of stable winds.

5 Conclusions

In this work, we reported an analysis of the GRAVITY-VLTI (K band) interferometric data of the M-type Mira star, R Car. Our primary scientific goal is to provide new observational evidence on the mass-loss mechanisms associated with M-type AGB Mira stars. In particular, we were interested in linking the structure of the gas around R Car with the origin of the dusty winds observed in M-type Mira stars. For this purpose, we characterize the morphology of the target within ~3 stellar radii.

To obtain a first-order estimate of the size of the source, we applied a geometrical model to the V² across the K band for our two epochs of observation with GRAVITY-VLTI, which are separated by ~30 days that correspond to ~ 10% of the full pulsation period. We find that the size of the target changes with wavelength due to the presence of molecules including H₂O and CO at ~2 μm and ~2.29–2.33 μm, respectively. The best estimate of the stellar disk diameter in the K band was derived at a wavelength of ~2.23 μm. Our parametric model allowed us to obtain a diameter of 16.67±0.05 mas (3.03 au) in January, 2018 and 14.84 ± 0.06 mas (2.70 au) in February, 2018.

We used a parametric single-layer model, which allowed us to have a physical estimate of the size, temperature, and optical thickness of the CO layer. Using the estimated optical thickness and temperature at the minimum peak of the first CO band head, we derived a column density of N_CO ~ 9×10¹⁸–1×10¹⁹ cm⁻². This estimate is in agreement with the column densities reported for other evolved stars, and it is high enough for allowing coalescence dust formation. The precision of our geometrical single-layer model is not enough to confirm the formation of CO. For this, higher angular resolution measurements in combination with more complex models of the CO layer would allow us to better constrain the sizes, temperatures, and optical thicknesses of the gas, which would result in better estimations of the column densities.

From the temperature and size of the CO layer, we set constraints on the dust candidates that could coexist with the gas. We find that the composites Mg₂SiO₃ and MgSiO₄ could coexist at the same distances as the ones we obtained for the CO innermost layers (~2–2.4R_*). However, the condensation temperatures for these dust particles are lower than the temperatures obtained from the CO layer model. We propose that the temperature discrepancy is due to the nature of the single-layer model, which considers a homogeneous and symmetric distribution of the CO layer, giving mean values for T_L, R_L‚ and τ_L. Our hypothesis is supported by the varying interferometric CPs and the posterior image reconstruction, which show that the CO distribution is neither homogeneous nor symmetric.

To complement the information obtained from our symmetric parametric models, we reconstructed interferometric images from the GRAVITY data. As we expected from the variation of the CPs, we observed several asymmetries in the images across the spectral ranges and for both epochs of observation. The presence of those asymmetries were confirmed (above 3 σ detection level) using bootstrap, PCA, and spectro-astrometric analyses. We propose that these asymmetries trace a complex, perhaps clumpy structure of the CO molecule distribution, which also changes within a period of 30 days or less.

Recent models of the convection and pulsation mechanisms of M-type AGB stars support that an asymmetric levitation of the gas at the innermost regions of these stars has a direct impact on the estimation of mass-loss rates and the formation of stable dust-driven winds. A dedicated characterization of the asymmetries observed in the CO distribution of R Car, together with new GRAVITY-VLTI data (with the longest baselines available at the VLTI) could provide unique observational evidence to test this scenario. Furthermore, complementary mid-infrared interferometric data with MATISSE-VLTI could help us to test the role of magnesium-based dust composites in the formation of dust-driven winds in M-type Mira stars and their subsequent chemical evolution.

Acknowledgements

The authors of this work thank the anonymous referee for his/her valuable comments to improve our manuscript. A.R.-G. acknowledges the support received through the PhD scholarship (No. 760678) granted by the Mexican Council of Science CONACyT. J.S.-B. acknowledges the support received from the UNAM PAPIIT projects IA 101220 and IA 105023; and from the CONACyT “Ciencia de Frontera” project CF-2019/263975. R.S. acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). R.S. acknowledges financial support from national project PGC2018-095049-B-C21 (MCIU/AEI/FEDER, UE). This project is based on observations collected at the European Southern Observatory under ESO programme 0100.D-0835.

Appendix A Gravity observation list and observables

Fig. A.1 GRAVITY calibrated interferometric observables. The left plot shows the V2 and CPs plotted versus spatial frequencies for the data sets collected and the right one shows the corresponding u-v coverage. Panel (a) corresponds to the January epoch, while panel (b) corresponds to the February one. The colors in the plots show the different wavelengths (see the labels on the plot).

Appendix B Best-fit geometrical models

In this section, we show the best-fit images to the observed V² and CPs across the pseudo-continuum, the first and second CO band heads for our January and February epochs (see Sect. 3.3.1). In all figures, the red dots in the main panels correspond to the synthetic V² and CPs extracted from the raw reconstructed images, while the data are shown with gray dots in the main panels. The lower panels show the corresponding residuals (in terms of the number of standard deviations) coming from the comparison between the data and the best-fit reconstructed images.

Fig. B.1 Best-fit to observed V2 and CPs from the best reconstructed images across the pseudo-continuum for the (a) January and (b) February epochs. The red dots in the main panel correspond to the synthetic V2 and CPs extracted from the raw reconstructed images. The data are shown with gray dots in the main panels. The lower panels show the residuals (in terms of the number of standard deviations) coming from the comparison between the data and the best-fit reconstructed images.

Fig. B.2 Fit to observed V2 and CPs from the best reconstructed images across the first CO band head for the a) January and b) February epochs. The panels are as described in Figure B.1.

Fig. B.3 Fit to observed V2 and CPs from the best reconstructed images across the second CO band head for the a) January and b) February epochs. The panels are as described in Figure B.1.

References

All Tables

All Figures

	Fig. 1 Visual light curve of R Car based on AAVSO database (blue points). The red line represents the best fit of a sine function to the curve which results in a period of P = 314 ± 1 days. The dashed vertical lines and the colored diamonds show the epochs of our 2018 GRAVITY observations (see the labels on the plot).
In the text

	Fig. 2 R Car calibrated spectra and V2. Top panel shows the normalized calibrated spectra of the GRAVITY data. The blue spectrum corresponds to the January epoch, while the red one corresponds to the February one (see the labels on the plot). The middle and bottom panels show plots of the V2 versus wavelength for the two epochs observed.
In the text

Fig. 3 V2 of the pseudo-continuum versus spatial frequencies of both epochs. The lines indicate the best-fit UD models. Each color corresponds to a different wavelength (see the label on the plot). We highlight that the spatial frequencies increase in opposite directions depending on the epoch. This allowed us to have a better visual comparison of the trends present in the V2 data between the two epochs.

In the text

	Fig. 4 Best-fit ΘUD as a function of wavelength for the two analyzed epochs. Each color corresponds to a different epoch (see the labels on the plot).
In the text

	Fig. 5 Illustration of single-layer model. The yellow disk represents the star, and the orange ring represents the layer. The parameters in the image are as described in Sect. 3.2. This image was adapted from Rodríguez-Coira et al. (2021).
In the text

	Fig. 6 Best fit to spectra with the single-layer model for the two CO band heads and epochs. The gray dots correspond to the data points and the solid-lines correspond to the best-fit model (see the labels on the plot).
In the text

Fig. 7 Best fit to V2 from the best single-layer models across the first CO band head for our two epochs. The colored dots and the shaded regions correspond to the best-fit models and their corresponding error-bars. Each color correspond to a different wavelength (see labels on the plot). The data are shown with black dots. We highlight that the spatial frequencies increase in opposite directions depending on the epoch. This allowed us to have a better visual comparison of the trends present in the V2 data between the two epochs. Also, for better visualization, we displaced the V2 vertically.

In the text

	Fig. 8 Best fit to V2 from the best single-layer models across the second CO band head for our two epochs. The panels are as described in Fig. 7.
In the text

	Fig. 9 Best pseudo-continuum reconstructed images for the (a) January and (b) February epochs. The white ellipse located in the right corner of the image at 2.0160 μm corresponds to the mean synthesized primary beam. The blue contours represent 10, 30, 50, 70 90, 95, 97, and 99% of the intensity’s peak.
In the text

	Fig. 10 Best reconstructed images per spectral channel for the (a) January and (b) February epochs across the first CO band head. The white ellipse located at the right corner of the image at 2.2941 μm and the blue contours are as described in Fig. 9.
In the text

	Fig. 11 Best reconstructed images per spectral channel for the (a) January and (b) February epochs across the second CO band head. The white ellipse located in the right corner of the image at 2.3229 μm and the blue contours are as described in Fig. 9.
In the text

Fig. 12 First four principal components for the January and February epochs for the first CO (plots a and b) and the second CO band heads (plots e and f). The relative scale of the structures in the eigenimages is displayed with a color bar at the bottom of each panel. Only the central 20 mas of the eigenimages are shown in each panel. The white contours correspond to the mean image across wavelength (per given data set) and they represent 10, 30, 50, 70, 90, 95, 97, and 99% of the intensity’s peak.

In the text

	Fig. 13 Representative example of CPs versus spatial frequencies of the data set that corresponds to the first CO band head at 2.2946 μm (gray squares). The colored dots show the CPs from the recovered images obtained with different numbers of principal components (see label on the plot).
In the text

	Fig. 14 Two-dimensional photocenter offsets in mas. The black dots correspond to the photocenter position of the continuum emission and the labeled colored dots to the photo-centers of the CO band head. The error bars of the reported displacements are smaller than the used symbols.
In the text

	Fig. 15 Pure-line CO images of first CO band head for the (a) January and (b) February epochs. The region that corresponds to the subtracted continuum contribution has been masked for better visualization of the CO region. The white ellipse located in the lower right corner of the image at 2.2941 μm corresponds to the mean synthesized primary beam.
In the text

Fig. 16 Curve of constant levitation radius (Rl = 2.43 R*) represented with a solid-black line. The plot shows the condensation temperature, Tc, as a function of the dust absorption coefficient, p. The red circles indicate the positions of a variety of dust grains. The blue and green shaded regions correspond to the CO temperature ranges estimated from the single-layer model (see labels on the plot).

In the text

	Fig. 17 Panels similar to Fig. 16 but showing the median condensation radii obtained directly from the pure-line CO reconstructed images at different position angles (see labels on the panels) for (a) the January and (b) February epochs.
In the text

	Fig. A.1 GRAVITY calibrated interferometric observables. The left plot shows the V2 and CPs plotted versus spatial frequencies for the data sets collected and the right one shows the corresponding u-v coverage. Panel (a) corresponds to the January epoch, while panel (b) corresponds to the February one. The colors in the plots show the different wavelengths (see the labels on the plot).
In the text

Fig. B.1 Best-fit to observed V2 and CPs from the best reconstructed images across the pseudo-continuum for the (a) January and (b) February epochs. The red dots in the main panel correspond to the synthetic V2 and CPs extracted from the raw reconstructed images. The data are shown with gray dots in the main panels. The lower panels show the residuals (in terms of the number of standard deviations) coming from the comparison between the data and the best-fit reconstructed images.

In the text

	Fig. B.2 Fit to observed V2 and CPs from the best reconstructed images across the first CO band head for the a) January and b) February epochs. The panels are as described in Figure B.1.
In the text

	Fig. B.3 Fit to observed V2 and CPs from the best reconstructed images across the second CO band head for the a) January and b) February epochs. The panels are as described in Figure B.1.
In the text