MOLsphere and pulsations of the Galactic Center's red supergiant GCIRS 7 from VLTI/GRAVITY

GCIRS 7, the brightest star in the Galactic central parsec, formed $6\pm2$ Myr ago together with dozens of massive stars in a disk orbiting the central black-hole. It has been argued that GCIRS 7 is a pulsating body, on the basis of photometric variability. We present the first medium-resolution ($R=500$), K-band spectro-interferometric observations of GCIRS 7, using the GRAVITY instrument with the four auxiliary telescopes of the ESO VLTI. We looked for variations using two epochs, namely 2017 and 2019. We find GCIRS 7 to be moderately resolved with a uniform-disk photospheric diameter of $\theta^*_\text{UD}=1.55 \pm 0.03$ mas ($R^*_\text{UD}=1368 \pm 26$ $R_\odot$) in the K-band continuum. The narrow-band uniform-disk diameter increases above 2.3 $\mu$m, with a clear correlation with the CO band heads in the spectrum. This correlation is aptly modeled by a hot ($T_\text{L}=2368\pm37$ K), geometrically thin molecular shell with a diameter of $\theta_\text{L}=1.74\pm0.03$ mas, as measured in 2017. The shell diameter increased ($\theta_\text{L}=1.89\pm0.03$ mas), while its temperature decreased ($T_\text{L}=2140\pm42$ K) in 2019. In contrast, the photospheric diameter $\theta^*_\text{UD}$ and the extinction up to the photosphere of GCIRS 7 ($A_{\mathrm{K}_\mathrm{S}}=3.18 \pm 0.16$) have the same value within uncertainties at the two epochs. In the context of previous interferometric and photo-spectrometric measurements, the GRAVITY data allow for an interpretation in terms of photospheric pulsations. The photospheric diameter measured in 2017 and 2019 is significantly larger than previously reported using the PIONIER instrument ($\theta_*=1.076 \pm 0.093$ mas in 2013 in the H band). The parameters of the photosphere and molecular shell of GCIRS 7 are comparable to those of other red supergiants that have previously been studied using interferometry.


Introduction
The stellar population of the central parsec of the Galaxy has been widely studied (Genzel et al. 2010;Morris et al. 2012, and references therein), where the presence of a disk of young stars is well recognized (Genzel et al. 2000(Genzel et al. , 2003Paumard et al. 2006;Lu et al. 2009; Bartko et al. 2009;Yelda et al. 2014). Most of these stars are massive O-type supergiants and Wolf-Rayet stars (Martins et al. 2007; Bartko et al. 2010;Sanchez-Bermudez et al. 2014). GCIRS 7 is one of the few evolved late-type stars (an Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under the programme IDs 098.D-0250 and 103.B-0032. GRAVITY is developed in a collaboration by the Max Planck Institute for extraterrestrial Physics, LESIA of Observatoire de Paris / Université PSL / CNRS / Sorbonne Université / Université de Paris and IPAG of Université Grenoble Alpes / CNRS, the Max Planck Institute for Astronomy, the University of Cologne, the CENTRA -Centro de Astrofisica e Gravitação, and the European Southern Observatory. Corresponding author: G. Rodríguez-Coira. M1 red supergiant or RSG, Blum et al. 1996b) and a SiO maser source (Menten et al. 1997) as well as the brightest star (in the H and K bands, with H = 9.5 and K = 6.5) of all the central parsec (Becklin & Neugebauer 1975). The works of Yusef-Zadeh & Morris (1991), and Serabyn et al. (1991) reported a cometary tail whose origin comes from GCIRS 7; more recently, Tsuboi et al. (2020) revealed the presence of an ionised shell in the core of the cometary tail, estimating the mass loss of GCIRS 7 with ALMA observations. This stellar population is permeated by the complex interstellar medium (ISM) environment and interacts with it. The central parsec is surrounded by a 2-7 pc-wide clumpy torus, the Circumnuclear Disk (CND), composed of dust and neutral gas (Becklin et al. 1982). The H ii region Sgr A West (the Minispiral; e.g., Lacy et al. 1980;Lo & Claussen 1983) consists of tidallysheared streamers and smaller patches and filaments of dust and ionised gas that orbit and penetrate the central parsec (Liszt 2003;Paumard et al. 2004;Mužić et al. 2007;Irons et al. 2012;Tsuboi et al. 2017). The volume surrounding these components inside the central cavity of the CND is not empty but, rather, Article number, page 1 of 12 arXiv:2105.09832v1 [astro-ph.SR] 20 May 2021 A&A proofs: manuscript no. main filled with hot (≈ 1.3 keV) plasma detected in X-ray (Baganoff et al. 2003). Finally, warm H 2 (with an excitation temperature of T e ≈ 2000 K) has been detected throughout the central parsec, presumably at the surface of many dusty clumps (Ciurlo et al. 2016(Ciurlo et al. , 2019. Ferrière (2012) gives an interesting overview of the ISM content of the central parsec. The large and dense clumps that form the Minispiral are believed to be infalling from the CND and beyond, but the origin of the lighter and less dense features (filaments, smaller clumps, X-ray plasma) is less clear. A fraction of them could originate in the feedback from the massive stars. G1 (also known as Sgr A*-f, Clénet et al. 2004Clénet et al. , 2005 and G2 (Gillessen et al. 2012) may well be extreme examples of such feedback clumplets (e.g., Meyer & Meyer-Hofmeister 2012;De Colle et al. 2014;Schartmann et al. 2015).
Thanks to the performance of stellar interferometers, the understanding of the structure and evolution of RSG has improved significantly. The closest ones have been widely studied, not only obtaining measurements of their sizes but also revealing single-layer atmospheres (Perrin et al. 2005(Perrin et al. , 2007Montargès et al. 2014), multi-layer atmospheres (Ohnaka et al. 2009(Ohnaka et al. , 2011(Ohnaka et al. , 2013Hadjara et al. 2019), complex structures in the photosphere (Haubois et al. 2009;Chiavassa et al. 2010;Ravi et al. 2011;O'Gorman et al. 2017;Ohnaka et al. 2017), and even the temporal evolution of the stellar surface (Ohnaka et al. 2011(Ohnaka et al. , 2013Montargès et al. 2016Montargès et al. , 2018Climent et al. 2020). Moreover, imaging of a RSG was performed in Baron et al. (2014), Monnier et al. (2014), and more recently, in Wittkowski et al. (2017a) and Climent et al. (2020). Although the sample of spatially resolved RSGs has been increasing over the last decade (Arroyo-Torres et al. 2013Wittkowski et al. 2017b), this sample is still not very large due to the shortness of the RSG phase, hence, only a small number of stars can be resolved using interferometers. When available, the study and characterization of the outer atmosphere of any new RSG and its temporal evolution would add valuable knowledge to the understanding of their mass loss processes, which have not yet been fully described from the first principles (Beasor et al. 2020). Paumard et al. (2014) compiled almost 40 years of nearinfrared photometric data on GCIRS 7 and exhibited two periods in the light curves: a short "fundamental" period, P 0 ≈ 470 days, and a long "secondary" period, P LSP ≈ 2800 days, as are often seen in RSGs. Those periods are believed to be a sign of pulsations, especially for the P 0 of the fundamental or first overtone radial mode (Yang & Jiang 2012, and references therein). Such pulsations are expected to play a major role in the mass loss of RSGs. However, they have never been confirmed on the basis of direct size measurements. GCIRS 7 has been observed via interferometry on the VLTI using AMBER in the K band and PIONIER in the H band (Pott et al. 2008;Paumard et al. 2014); however the AMBER data do not have sufficient spectral resolution and (u, v)-coverage to disentangle the stellar disk from the circumstellar environment so that only the PIONIER data provide a trustworthy uniform-disk diameter (θ UD (2013) = 1.076 ± 0.093 mas).
The GRAVITY instrument has tremendously increased the sensitivity of the VLTI (Gravity Collaboration et al. 2017), allowing us to observe GCIRS 7 at moderate spectral resolution (R = 500) in single-field mode using the four 1.8-m auxiliary telescopes (AT) at two epochs (2017 and 2019), with the goal of detecting variations in the photospheric diameter of the star and in its circumstellar environment. The data sets, the data reduction, and the calibration processes are described in Sect. 2. The methods and models used to measure the parameters of the star are described in Sect. 3. The results and their implications

Data
The log of the data is presented in Table 1 with science on-target (SCI), calibrator on-target (CAL), and sky (SKY) frames. The data were taken at two different epochs (two SCI data frames with run ID 098.D-0250(B), corresponding to the night of 18 March 2017 and seven SCI data frames with run ID 103.B-0032(F), corresponding to the nights of the 5 and 6 July 2019) with two different baseline configurations, as shown in Figure  1. The maximum baseline was 132.5 m in 2017 and 129.3 m in 2019. Turbulence in the beams was corrected with the NAOMI adaptive optics (AO) system (Gonté et al. 2016;Woillez et al. 2019) on axis for the 2019 data, while only tip-tilt stabilization using STRAP on a nearby visible star was possible in 2017. Fringe-tracking was performed on-axis at both epochs. Only the science beam combiner data were used in this study since the fringe-tracker data do not have sufficient spectral resolution for our purpose.
We chose a single target to be the calibrator for both spectroscopy and interferometry. We used this calibrator to remove the atmospheric features using the appropriate template (see Sect. 2.1.2). The calibrator is always observed at the same air mass and integration time as the science sequence (DIT = 10 s, NDIT = 30, air mass = 1.01 in 2017 and DIT = 5 s, NDIT = 30, air mass = 1.01 in 2019). We recorded the 2017 data in combined polarization mode and the 2019 data in split polarization mode, making use of the Wollaston prism. For this reason, the second run in 2019 results in two simultaneously-recorded data sets, one for each polarization (P1 and P2).

Interferometric quantities
The data were reduced with the GRAVITY pipeline (Lapeyrere et al. 2014). As the source is moderately resolved, both the closure phase signal and the differential phase are 0 • ± 3 • and we only used the squared visibilities (V 2 ) for the interferometric analysis. The pipeline was also used to produce the pho- tometric spectra. The squared visibilities were calibrated with HD160852 (θ LD = 0.148 ± 0.004 mas, V 2 = 0.997 ± 0.001, Chelli et al. 2016) in 2017 and HD161703 (θ LD = 0.38 ± 0.01 mas, V 2 = 0.983 ± 0.001, Chelli et al. 2016) in 2019, where θ LD is the limb-darkened diameter, which, in the original work, was obtained via a polynomial fitting. Both calibrators are observed at the same air mass as GCIRS 7.

Photometric spectra calibration
We used several templates from Pickles (1998) to perform the absolute flux calibration of the spectra. In order to get a robust estimate of the error bars related to the template choice, we chose templates A0V, A3V, F2V, and F5V for 2017 and K0III, K2III, and K3III for 2019. The following formula yields the calibrated flux density of GCIRS 7, measured in Wcm −2 µm −1 : where F(λ) template is the corresponding spectral template from Pickles (1998) for each epoch (four templates in 2017, three templates in 2019). For each template, a spectra F IRS7 λ is obtained. We take the average of these individual estimates as our final calibrated spectrum for each epoch, and their standard deviation as the uncertainty on this calibration.
While the 2017 data were recorded in combined polarization mode resulting in a single spectrum, the 2019 data were recorded using the Wollaston prism, resulting in two independent spectra which were calibrated individually and then co-added into a single spectrum.

Local interstellar extinction and spectrum normalization
It is possible to measure the effect of the interstellar extinction in the direction of the source using the slope of the continuum emission and considering an extinction law. Here, we consider the law derived in Fritz et al. (2011): with α = −2.11 ± 0.06, λ 0 = 2.166 µm (Brackett γ). An average value of A 0 = 2.62 ± 0.11 is given in Fritz et al. (2011) for the Galactic Center. However, it is known that extinction varies significantly throughout the central parsec, from a lineof-sight to the next and along a given line-of-sight (e.g., Ciurlo et al. 2019). We therefore decided to derive the extinction from the source itself. Starting from a gray atmosphere approximation (F λ ∝ B(λ, T )10 −0.4A λ ), an expression for A 0 can be obtained by introducing A λ from Eq. 2 and taking the derivative of the logarithm of F λ with respect to λ: where h, k, c are the Planck and Boltzmann constants and the speed of light respectively. We used the effective temperature To estimate the uncertainties, in addition to the flux error bars, we considered the deviation of the SED of a typical RSG from a blackbody. For this purpose, we used the de-reddened and absolute-calibrated spectrum of Betelgeuse from Rayner et al. (2009) to estimate the deviation. This star shares the spectral class with GCIRS 7 and was found to have a same effective temperature, which has not substantially changed in recent decades (several measurements are presented in Sec. 4.3). An estimate of the slope for the sub-band used in our work to estimate the extinction (λ = 2.1 − 2.2 µm) gives a deviation of 11% while comparing the slope of the blackbody with a linear fit of the spectrum of Betelgeuse in the continuum sub-band. This deviation in the slope is propagated by using Eq. 3, giving an error bar of 0.20. This value has been considered in addition to the error propagation of the fluxes to estimate the uncertainty.
The two values are compatible and can be averaged down to a single estimate of the extinction:Ā 0 = 3.27 ± 0.20. We can use this value to de-redden F IRS7 λ (Eq. 1) and integrate over the K S band to determine the broad-band extinction specifically for the photosphere of GCIRS 7: Then each spectrum is normalized separately using the average A 0 for each epoch. Finally, a weighted average of the normalized spectra is computed for each epoch, the weights being the inverse of the square of the error bars. The weighted average of the normalized spectra is noted as F N λ throughout this work.

Squared visibilities
The target is only moderately resolved by GRAVITY (V 2 > 0.4) and therefore the data are not sensitive to limb darkening. We modeled the visibility in the continuum (2.1-2.2 µm) with a uniform disk plus a spatially over-resolved background as in Perrin et al. (2007): where V aD (q, a, θ * UD ) is the visibility of this uniform-disk plus the spatially over-resolved background, V UD (q, θ * UD ) is the usual uniform-disk visibility, q = √ u 2 + v 2 is the spatial frequency (u and v are the spatial frequency components), θ * UD is the angular diameter of the disk, and a is the fraction of stellar flux over the total injected flux.
To explore the wavelength dependence of the parameter a, we performed a uniform-disk plus background fit for each spec- tral channel in the continuum region (2.1-2.2 µm). The results for a do not reveal any significant wavelength dependency and therefore we assumed the stellar flux and background to exhibit the same spectrum.
We estimated θ * UD and a using a Monte-Carlo Markov chain (MCMC) algorithm based on the Python package emcee (Foreman-Mackey et al. 2013). We let 100 walkers evolve for 600 steps for each file and each polarization state (for 2019) individually, using all wavelengths from 2.1 to 2.2 µm. The individual Monte-Carlo simulations (2 in 2017, 14 in 2019) yield distributions that are too far apart compared to their internal scatter (Fig. 2). This is because systematic errors resulting from the variation in the transfer function dominate over statistical errors. For this reason, we combined all the samples from these simulations in a single histogram per epoch. The 2017 combined histogram is clearly bi-modal because there are only two individual frames, while the 2019 data set is rich enough that the 14 individual Gaussian-like histograms merge into a broad Gaussian-like peak. The two polarization states in the 2019 data give similar distributions and are combined together. As our best estimate for the two parameters, we use the median of the 1D combined histograms (which is equivalent to taking the average of the median values of the individual 1D histograms). The 2019 data set is well suited to determine the uncertainties because it contains enough individual measurements for their scatter to make sense. Using the 16% and 84% percentiles of the 1D histogram yields uncertainties that combine statistical errors, instrumental stability, and degeneracy between the two parameters. The same approach on the 2017 data leads to a value for the uncertainty on θ * UD that we deem too small, being dominated by the scatter of two individual values. We therefore adopt the 2019 error bar instead also for 2017. Our final estimates for the two years are as follows: -2017: a = 0.927 ± 0.008, θ * UD = 1.547 ± 0.030 mas; -2019: a = 0.968 ± 0.006, θ * UD = 1.549 ± 0.030 mas.  It is worth noting that this uncertainty on the uniform disk diameter (30 µas) is itself quite an achievement, on the same order as the astrometric measurements performed by GRAVITY (e.g., Gravity Collaboration et al. 2020). These average best-fit models are shown in Fig. 3., compared with the best single UD fit without considering the background.
The fraction of coherent flux a means that the fibers of GRAVITY were fed with 7% incoherent light from the circumstellar background in 2017 and 3% in 2019. This decrease may be due to the improved turbulence correction offered by the AO system NAOMI compared to tip-tilt stabilization with STRAP, since the residuals of this correction directly translate into a widening of the fibre field-of-view (Perrin & Woillez 2019). We removed this background from the visibilities by dividing them by a for each epoch. The use of a single value instead of a value per spectral channel is justified as the spectral resolution is moderate and the flux ratio between the emission of the circumstellar dust and stellar flux varies only slowly across the K band.
After correcting for this background, a simple uniformdisk fit was performed for each spectral channel in the interval (2.05-2.4 µm). The results are displayed in Fig. 4, showing both the normalized spectra derived in Sect. 3.1 and the chromatic uniform-disk diameter profile as a function of wavelength. The diameter is close to constant in the continuum between 2.1 and 2.2 µm, again showing that the continuum is very well fitted by a uniform-disk plus the spatially over-resolved background. The unbiased weighted standard deviation of θ UD through the continuum (again 2.1-2.2 µm) is 0.020 mas in 2017 and 0.018 in 2019. At this level, the main source of uncertainty on θ * UD is the presence of many absorption lines in the spectrum, which form at various altitudes in the atmosphere of the star. Therefore, the limitation is fundamental: it is the definition of the photosphere itself.
The average value of θ UD above 2.29 µm is 1.66 mas in 2017 and 1.77 mas in 2019, that is, a 12-σ departure from the continuum level. This sharp increase follows the absorption features that can be seen in the spectrum, with local maxima matching the deep CO band heads. Such variations correlated with the CO optical depth point towards a molecular shell above the photosphere as previously evidenced, for instance, by Perrin et al. (2004), Perrin et al. (2005), or Hadjara et al. (2019) for other evolved stars.
In addition to the maxima in the CO bands, the uniform disk diameter increases monotonically after 2.29 µm in 2019. Such a feature has also been observed in Betelgeuse (Montargès et al. 2014) and other RSGs (Arroyo-Torres et al. 2013 and it can be attributed to the presence of H2O (Montargès et al. 2014). However, this monotonic increase in the angular size is not clear in 2017, which hints at a variation in the stellar atmosphere between the two epochs.
A lower limit on the diameter of the shell is given by the maximum uniform-disk size in the molecular band: θ S ≥ θ CO UD , which we estimated by taking the average of the three θ UD values around the CO 5-3 band head on each epoch: This result provides evidence for the need to carry out a more complex model than that with a grey atmosphere to describe the physics of the object. In the next section, we model the molecular features with a simple shell model, which will also allow us to give a proper interpretation to the time variability which we see.

Single-layer shell model
A geometrically thin molecular layer model has proven successful to reproduce the visibilities of other RSG stars (Perrin et al. 2004(Perrin et al. , 2005Montargès et al. 2014) and this is what we chose to use for this study. It has allowed for the interpretation of interferometric observations of RSG stars surrounded by a so-called MOLsphere according to the term coined by Tsuji (2000). In this model, the shell is modeled with a temperature and an optical thickness but has zero geometrical thickness.
Both the stellar photosphere and the shell are modeled with black body functions. The specific intensity at angular distance r from the center of the star as seen from the observer at wavelength λ is given by: where T * and T L are the temperatures of the photosphere and of the molecular layer, R * and R L are their angular radii (hence R * = θ * UD /2), τ λ 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 cos β = 1 − (r/R L ) 2 . A sketch of the model is displayed in Fig. 5.
The center-to-limb variation is illustrated in Fig. 6 for various optical depths. The sharp variation near r/R L = 0.6 corresponds to the edge of the photosphere, which is assumed to be a uniform disk. The increase of the projected layer optical depth with increasing r (due to the cos β factor) causes a slight limb darkening on the disk over the photosphere, and a limb brightening between the limb of the photosphere and the limb of the shell. The model fails to reproduce infinitely strong flux attenuation by the shell unless its temperature reaches 0 K. For example, the maximum attenuation reaches a factor 1.55 for a temperature of T L = 2360 K and is typically reached for optical depths of 5 or larger. The fact that at high optical depth the shell itself behaves as a photosphere is one of the shortcomings of this simple We adopted the effective temperature of Paumard et al. (2014) for the photosphere: T * = 3600 K. We also fixed the photospheric diameter to the continuum uniform-disk diameter from Sect. 3.2: R * = θ * UD /2 = 0.775 mas. In this section, we focus on the molecular features above 2.28 microns (Montargès et al. 2014), which includes CO bands but also water vapor to obtain measurements for T L and R L .
A physical modeling of the wavelength-dependent optical depth of the thin layer τ λ , while possible, would significantly add to the complexity of the model without begin necessarily accurate because of the very simplistic geometrical model. However, for a set of parameters (T * , T L , R * , R L ) and for each wavelength, the relation between F N λ and τ λ is bijective, except where the model saturates as explained above. For any quadruplet (T * , T L , R * , R L ) and for each wavelength, we determine τ λ univocally by finding the root of the following quantity: which is implemented by a minimization since Eq. 7 is always positive. When the model does not saturate, this minimum reaches zero. When it does saturate (this happens in a few spectral channels of the CO band heads), the exact (large) value of τ λ does not matter as the shell is then optically thick. A visibility model of the star and the shell V 2 shell can then be built and compared to the corrected squared visibility data V 2 c,i = V 2 i /a 2 , in which the effect of the incoherent background has been removed. This yields the following χ 2 : Article number, page 6 of 12 GRAVITY Collaboration: MOLsphere and pulsations of GCIRS 7 from VLTI/GRAVITY Fig. 7. Contour plot of the χ 2 maps (Eq. 8) considering T L , θ L , which are the two parameters of interest. Solid lines: Decrease of likelihood of 1, 3, and 5 σ relative to the hot solution minimum of each epoch (Press et al. 1992). Grey dashed line: Temperature profile of a spherically thin shell enclosing a perfect black body (Eq. 9). Triangles with error bars: Cold solutions for both epochs corresponding to the absolute minima. Crosses with error bars: Hot solutions for both epochs (corresponding to the local minima).
where M is the total number of data points considered for each epoch: M = 6 × D × n λ where D is the number of observations, 6 the number of baselines and n λ the number of spectral channels in the band of interest. T L and R L are the only two remaining free parameters of the model. The fitting process involves then two steps: -For each pair (T L , R L ), τ λ is measured from the spectrum for each wavelength by minimising Eq. 7. -By using the resulting set (T L , R L , τ λ ), the visibility squared of the model is computed and compared with the data via Eq. 8. Figure 7 displays a contour plot of this χ 2 for each epoch in the (T L /T * , θ L /θ * = R L /R * ) plane (where θ i = 2R i ). The map does not show a single clearly defined minimum, but a trough instead, thereby displaying a degeneracy between the radius and the temperature of the single layer of the shell model. For both 2017 and 2019, the trough has one global minimum at (T L , θ L ) = (838 ± 23 K, 5.36 ± 0.39 mas) with χ 2 r = 3.78 in 2017 and (T L , θ L ) = (839 ± 31 K, 5.96 ± 0.58 mas) with χ 2 r = 1.98 in 2019, which we refer to as the 'cold' solution and one local minimum at (T L , θ L ) = (2368 ± 37 K, 1.74 ± 0.03 mas) with χ 2 r = 3.93 in 2017; and (T L , θ L ) = (2140 ± 42 K, 1.89 ± 0.03 mas) with χ 2 r = 2.18 in 2019, which we refer to as the 'hot' solution. In 2017, the two solutions are incompatible at 5σ.
Despite the fact that the cold solution corresponds to the absolute minimum, it cannot be the dominating solution of the model as its temperature is below the condensation temperature of silicate dust. In the context of mass loss, CO could be detected beyond the dust condensation radius as the CO gas can be dragged by the dust set in motion by the radiation pressure of the star. However, CO is naturally a major component of the molecular shell and therefore there should be an inner region where CO must be formed before being dragged by the dust. This inner shell must present a higher column density and a stronger spectroscopic signal in the K band than the cold solution. These requirements are satisfied by the hot solution.
In addition, we plotted in Fig. 7 the position in radius and temperature for a fully absorbing and geometrically thin shell in thermal equilibrium above the photosphere: This model is more consistent with the hot solution than with the cold solution. For all these reasons, we believe the hot solution makes more physical sense than the cold solution although both are compatible with our data. However, the reduced χ 2 of the hot solution is not ideal, which we interpret as the signature that our model is still too simple to perfectly reproduce the data.
In an attempt to include these modeling errors into the statistical uncertainties, we rescaled the errors by the square root of the χ 2 of their respective hot solution.
The optical depths obtained for the hot solution of each epoch are presented Fig. 8. The error bars were obtained through Monte-Carlo error propagation in the 1σ contour of the hot solution. Only the values with τ ± ∆τ < 5 are presented as the model reaches saturation near the peaks. A summary of the parameters of the model is presented in Table 2 for the hot solution.
A visibility contrast function is introduced to better show the effect of the molecular features on the visibilities. It is the average of the band to continuum visibility ratio over a subset of baselines: This contrast function, where the spectral and spatial information are combined in a synthetic way, is an interferometric counterpart to the normalized spectrum. Figure 9 shows this function for the data of the two epochs as well as for synthetic data corresponding to the hot solution. For this plot, we chose the three partially overlapping longest baselines for each epoch. The visibility contrasts are plotted between 2.28 and 2.4 µm to show the upper part of the continuum and the molecular features. The single-layer shell model matches the data quite well except at the bottom of the band heads, where flux attenuation saturates with this simple model, as explained above.

Interstellar extinction
We compare our measurement of the extinction with previous ones. In Sec. 3.1, we measured the extinction from the  changes incurred to the slope of the spectrum in the continuum. The extinction found in Brackett γ was then converted to A K S = 3.18 ± 0.20 through integration across the K S band, quite stable between our two epochs. Others also estimated the extinction in the direction of GCIRS 7. Blum et al. (1996a) obtained A K = 3.72 ± 0.13 via near-infrared (NIR) photometry and an assumed intrinsic color, and later Blum et al. (2003) measured A K = 3.48 ± 0.09. Both values are more than 1σ higher than our measurement. This discrepancy with archival data could be due to the better spatial filtering offered by GRAVITY thanks to the mono-mode fibers and interferometric spatial resolution, which allows for most of the infrared excess to be rejected from the surrounding material, or to a change in the circumstellar contribution, subject to any variation of the fundamental parameters of the star as it is pulsating (Paumard et al. 2014). On the other hand, the interstellar medium (ISM) in the central parsec is known to be clumpy and heterogeneous, responsible for variations of more than one magnitude in the K band from one lineof-sight to the next and along a given line-of-sight (Ciurlo et al. 2016, and references therein). The motion of dusty clumps in the ISM over decade-long time scales (Gillessen et al. 2012(Gillessen et al. , 2013Ciurlo et al. 2019Ciurlo et al. , 2020 could easily explain variations of the interstellar extinction by a fraction of a K-band magnitude. A value of 3.18 is above average when compared to the extinction map of Schödel et al. (2010); but it is rather typical in view of the extinction map of Ciurlo et al. (2016) who, by using the ISM instead of the stars, was arguably able to prove deeper along the line-of-sight. Therefore, a large part of the excess extinction compared to the field average towards GCIRS 7 is attributable to ISM local to the central parsec, as is possibly of circumstellar origin in some part.

Photospheric size
We measure a photospheric diameter of θ * UD = 1.55 ± 0.03 mas, yielding R * = 1368 ± 27 R at 8.246 kpc (Gravity Collaboration et al. 2020). This value is in good agreement with the values of θ * UD = 1.5 − 2 mas of Paumard et al. (2014) on K-band data obtained with AMBER in 2008, but much larger than their Hband measurement on PIONIER 2013 data: θ * UD = 1.076 ± 0.093 mas. One possible explanation for a smaller photosphere at Hband compared to K-band is the wavelength-dependence of the opacity of negative hydrogen H − , which determines the opacity of RSGs (Gray 2005). However, this process alone can hardly account for a factor 1.44 in size between the two bands. Therefore, these measurements provide another clue that GCIRS 7 has been pulsating over the last ten years. Such variations in stellar diameter are also corroborated by photometric estimates using SINFONI (Paumard et al. 2014) and ALMA (Tsuboi et al. 2020). The θ * UD diameters from Paumard et al. (2014), Tsuboi et al. (2020) and this work are displayed in Fig. 11, together with a model assuming the pulsation periods (470 and 2620 days) obtained by Paumard et al. (2014). The phases and amplitudes of the two modes of pulsations and the average size were semimanually adjusted, returning a value for χ 2 r that is well below 1. However, given the little data available, too many solutions are possible to make it useful to quote the best-fit parameters. This shows that the data are consistent with pulsations but we are not able to provide further constraints at this stage.
We measure the same photospheric size in 2017 and 2019, with a 3σ upper limit on the difference of 0.13 mas ( 8%). This is plausibly explained by the observational gap. With the interplay of the two periods (470 and 2620 days) revealed in Paumard et al. (2014), the size of the photosphere could well have varied during the 840 days that separate the observations performed in 2017 and in 2019 and may have come back to a very similar value (Fig. 11). However, it is also possible that the size of the star did not vary between 2017 and 2019 or varied less than expected from the 2013-2017 era. The pulsations of red supergiants are irregular and intertwined with convective mechanisms. Their periods, phases, and amplitudes can vary over short timescales.

The MOLsphere of GCIRS 7 in context
Here, we compare the results of this work with the parameters of the MOLspheres of Antares, Betelgeuse and µCep as derived in various works following similar approaches as ours. The contours shown correspond to 1 sigma from their respective hot solutions of GCIRS 7 obtained in this work, which are represented by the crosses with errorbars in blue and orange. The green, red, and purple symbols correspond to previous works for Antares, Betelgeuse, and µCep respectively (Table 3). The circles correspond to models where only visibilities were used, while the crosses represent models which use both visibilities and spectra. The dotted lines connect the inner and the outer layers on the multiple shell models in a same work.
A single thin shell, similar at the one used in this work, was previously modeled in Perrin et al. (2005) for µCep via interferometry in the K band. The work of Verhoelst et al. (2006) combined spectroscopy in NIR, MIR, and FIR with interferometry in the K and L bands to characterize the extended atmosphere of Betelgeuse revealing the presence of alumina, while Ryde et al. (2006) presented new MIR spectroscopic data which could not rule out a scenario based only on a cool photosphere (T = 3250 K). For the same star, the most recent work involving a thin shell model was presented by Montargès et al. (2014), who used K-band interferometry to estimate the temperature, diameter and density of the thin shell.
Several multi-layer models have been also proposed for RSGs. Tsuji (2006) used NIR and MIR spectroscopy, together with K-band interferometry, to constrain a two-layer model (inner and outer) for µCep and Betelgeuse. A year later, Perrin et al. (2007) used N-band interferometry to model the close stellar environment of Betelgeuse with the use of another two-layer model estimating the density of the H 2 O, SiO, and Al 2 O 3 molecules in the MOLsphere. Ohnaka et al. (2009) also resolved spatially the structure of the atmosphere of Betelgeuse and characterized a MOLsphere by using two layers and K-band interferometry. This extended component was also observed by Ohnaka et al. (2011) (K-band interferometry) who found a significant variation of the visibilities and phases in the CO overtone lines between two different epochs, but no significant variation in the continuum. For the star Antares, Ohnaka et al. (2013) used a similar approach (K-band interferometry) characterizing a double layer model at two different epochs, but without any significant variations in temperature or size. A model consisting in several concentric layers was used in Hadjara et al. (2019) to characterize a sample of M stars via K-band interferometry. In the case of Antares, they modeled seven layers between 1.06 and 1.76 stellar radii.
These measurements are listed in Table 3 and displayed in Fig. 10, as well as the results of this work. The properties of GCIRS 7 are compatible in 2019 within one sigma with all the single shell models, except for Verhoelst et al. (2006), and they are also compatible with the outer shell of Antares from Ohnaka et al. (2013). On the other hand, 2017 reveals a smaller shell, close to the inner layer of Antares from Hadjara et al. (2019).
In addition, considering the multiple shell models, both the inner and the outer shell of Perrin et al. (2007) and the outer shell of Ohnaka et al. (2013) are also compatible within 1σ to the through pointing to the cold solution in 2019. This is an indication of the shortcomings of the single shell model to interpret 2019 data for GCIRS 7.

Size and expansion of the molecular shell
The ratio of the shell diameter to the photospheric diameter (1.16-1.26) is very similar to previous findings. Perrin et al. (2005) found a ratio of 1.32 for µ Cep, Montargès et al. (2014) found 1.25 for Betelgeuse. This is still compatible with the lower altitude layers of α Sco for which Hadjara et al. (2019) found ratios in the range 1.06 to 1.76, and temperatures between 2360 and 1900 K using a multi-layer model (see the next section for the discussion of the temperature). In addition, Wittkowski et al. The uniform disk diameters for the molecular band presented on Fig. 4 show globally higher values in 2019 than in 2017. This is also evidenced by the significantly lower visibility contrast in 2019 beyond 2.33 µm (Fig. 9). This is confirmed by the analysis with the single-layer model where the hot solutions found in 2017 and 2019 in Fig. 7 are different by 4σ. Our results are compatible with a layer expansion of 8% leading to a layer temperature decrease from 2017 to 2019.
However, the single shell model constitutes a simplified approach for a MOLsphere. Indeed, although a strong contribution of the shell due to a larger size or a lower temperature is found in 2019, the apparent difference in the shell between the two epochs may be a consequence of the degeneracies of the parameters. A more complex model would be needed to solve the degeneracies and, therefore, to be able to conclude about a possible shell expansion.

Column density of CO
The optical depth profiles shown in Fig. 8 are consistent in 2017 and 2019 except for the peaks where the model saturates. With these optical depth profiles, we used the method of Goorvitch (1994) to estimate the column density of CO. The method uses the first band head of CO at 2.293 µm and the temperature of the shell T L measured for each epoch. The optical depths of Fig. 8 reveal values slightly larger for 2017 than for 2019 as long as the wavelength increases, except for the first head band, which is the same at the two epochs. This is the effect of the monotonic increase of the diameter mentioned in the uniform disk approach (Fig. 4). In addition, in 2017 none the head bands except for the first one can be reproduced by the optical depths of the thin shell due to saturation, while the first two band heads are reproduced in 2019. This is an indicator that the density of the MOLsphere of GCIRS 7 might have decreased.
It is possible to test if this decrease can be due to an expansion of the molecular shell by scaling the measured density by the ratio of the surface of the shell at the two epochs. Assuming no sudden mass loss episode, if the shell has expanded since 2017 at the measured diameter (from 1.74 to 1.89 mas), the column density of CO in 2019 would be N CO = 2.27 × 10 20 mol cm −2 , which is 1σ compatible with the measured value. Therefore, our single-shell interpretation allows for a possible expansion of the molecular envelope, but the fact that the model is simplified does not allow for it to be fully confirmed without a deeper study involving an extended atmosphere approach.
Indeed, although CO is the main component of the molecular shell, other molecular species such as water have also been observed in the MOLspheres of other RSGs (Verhoelst et al. 2006;Ryde et al. 2006;Perrin et al. 2007;Montargès et al. 2014), and the decrease of the temperature of the molecular shell may also play a role in the density variation. For GCIRS 7, the molecular shell in 2017 could have been warm enough for water to dissociate (T = 2400 K), although it may still have existed at T = 2200 K in 2019.
The location of GCIRS 7 is particular as it is located at less than 1 pc from a supermassive black hole and plenty of massive young stars (Krabbe et al. 1995) and it is known that the most external layers are being blown away by the effect of their wind (Tsuboi et al. 2020;Yusef-Zadeh & Morris 1991;Serabyn et al. 1991). The question of whether or not their presence affects the observed molecular shell of GCIRS 7 could be addressed by the use of a more complex model to fully characterize the outer atmosphere of GCIRS 7, such as that of Hadjara et al. (2019). In addition, either increasing the signal to noise ratio in the K band or obtaining more interferometric data at longer baselines, as well as exploring other wavelengths (H, L bands), would bring more constraints to better characterize the outer atmosphere of GCIRS 7. Such approaches would help to understand the nature of the molecular shell observed, as well as its temporal evolution.

Conclusion
We report on the spectro-interferometry of GCIRS 7 in the K band using GRAVITY at ESO/VLTI. With the sensitivity of current interferometers, we proved that wavelength-dependent structures can be observed in evolved stars even if they are not fully resolved.
We find that GCIRS 7 presents the behavior of a typical RSG. We detect a molecular shell above the photosphere and estimate the sizes and temperatures for the photosphere and the shell as well as the column density for CO, based on optical depths constrained with a single-layer model. This is, to date, the RSG star with the smallest apparent size and the farthest for which a molecular shell has been spatially resolved from the star and characterized. We also obtained an estimate of the local interstellar extinction with the spectral data.
The extinction (A K S = 3.18 ± 0.16) and size of the photosphere (θ * UD = 1.55 ± 0.03) were the same within uncertainties at the two epochs in 2017 and 2019. However, the photospheric size must have changed from > 1.5 mas in 2008 to 1.1 ± 0.1 mas in 2013 (Paumard et al. 2014) and back to 1.55 ± 0.03 mas in 2017-2019.
The spectro-differential visibility signal demonstrate the presence of CO above the photosphere. In the context of a thin spherical shell model, the temperature (≈ 2200-2400 K) and diameter (10-20% larger than the photosphere) of this shell are in line with what has been found for other similar RSG stars. The size and temperature of the shell have significantly changed between the two epochs and are compatible with an expansion.
The column density for the molecular shell presents a value in the same line than the column density of the shells for other RSGs measured with similar methods. The model fails to reproduce all the band heads in 2017 except the first, while in 2019 the first two band heads are reproduced. This suggests that the density must have been higher in 2017. An interpretation based on a shell expansion from 2017 to 2019 is compatible with our data.
This work corresponds to a first-order description of the outer atmosphere of GCIRS 7. Overall, our results support the interpretation in terms of stellar pulsations proposed by Paumard et al. (2014) and hint at an expansion of a molecular shell. Follow-up observations over a good fraction of the ≈ 2800day period with contemporaneous H and K photometry, spectroscopic effective temperature, and interferometric size measurements, together with a more detailed multi-layer model made relevant by the spectral resolution of GRAVITY, would be needed to further confirm the pulsations and study the associated massloss processes.