The effect of winds on atmospheric layers of red supergiants I. Modelling for interferometric observations

Context. Red supergiants (RSGs) are evolved massive stars in a stage preceding core-collapse supernova. The physical processes that trigger mass loss in their atmospheres are still not fully understood, and remain one of the key questions in stellar astrophysics. Based on observations of α Ori, a new semi-empirical method to add a wind to hydrostatic model atmospheres of RSGs was recently developed. This method can reproduce many of the static molecular shell (or ’MOLsphere’) spectral features. Aims. We use this method of adding a semi-empirical wind to a MARCS model atmosphere to compute synthetic observables, comparing the model to spatially resolved interferometric observations. We present a case study to model published interferometric data of HD 95687 and V602 Car obtained with the AMBER instrument at the Very Large Telescope Interferometer (VLTI). Methods. We compute model intensities with respect to the line of sight angle ( µ ) for di ﬀ erent mass-loss rates, spectra and visibilities using the radiative transfer code T urbospectrum . We can convolve the models to match the di ﬀ erent spectral resolutions of the VLTI instruments, studying a wavelength range of 1 . 8 − 5 µ m corresponding to the K , L and M -bands for GRAVITY and MATISSE data. The model spectra and squared visibility amplitudes are compared with the published VLTI / AMBER data. Results. The synthetic visibilities reproduce observed drops in the CO, SiO, and water layers that are not shown in visibilities based on MARCS models alone. For the case studies, we ﬁnd that adding a wind on to the MARCS model with simple radiative equilibrium dramatically improves the agreement with the squared visibility amplitudes as well as the spectra, the ﬁt being even better when applying a steeper density proﬁle than predicted from previous studies. Our results reproduce observed extended atmospheres up to several stellar radii. Conclusions. This paper shows the potential of our model to describe extended atmospheres in RSGs. It can reproduce the shapes of the spectra and visibilities with better accuracy in the CO and water lines than previous models. The method can be extended to other wavelength bands for both spectroscopic and interferometric observations. We provide temperature and density stratiﬁcations that succeed for the ﬁrst time in reproducing observed interferometric properties of red supergiant atmospheres.


Introduction
Stellar winds impact the lives of massive stars and can change their evolutionary path in the Hertzprung-Russell diagram (HRD).These mass-loss events become important as the star leaves the main sequence phase (Chiosi & Maeder 1986).For this reason, one of the key stages for mass-loss is the red supergiant (RSG) phase, where the massive star is in a stage preceding core-collapse supernova.These mass-loss events occur in the extended atmospheres of RSGs, whose extensions go up to several stellar radii.Beyond that, the temperature is low enough to produce a dusty shell.However, the observed extensions are not at all reproduced by current dynamic model atmospheres which include pulsation and convection (Arroyo-Torres et al. 2015).
As a consequence, the mechanism that triggers mass-loss in the extended atmospheres of RSGs is still poorly understood.This is not the case for their low-and intermediate-mass counterparts (Miras), whose mass-loss processes can be explained by pulsation and dust driven winds alone (Wood 1979;Bowen 1988;Höfner & Olofsson 2018).There have been several at-tempts to explain the mechanism of stellar winds in RSGs (e.g., Josselin & Plez 2007;Chiosi & Maeder 1986;Kudritzki & Puls 2000), but there is still no consensus.Recent work by Kee et al. (2021) studied the effect of turbulent atmospheric pressure in initiating and determining the mass-loss rates of RSGs, finding promising results.However, further work is needed to unambiguously determine the dynamical processes that trigger massive stellar wind events in spatially extended atmospheres.
Hence, there have been few studies that explore the mass loss effect in cool massive stars, all of them focusing on their spectra.Most explore the mid-and far-IR region, where the dust component is dominant (e.g., Groenewegen et al. 2009;Beasor & Davies 2018;Decin et al. 2006;De Beck et al. 2010).Therefore, the models rely solely on dust modelling, such as DUSTY (Ivezic & Elitzur 1997;Ivezic et al. 1999) or RADMC3D (Dullemond et al. 2012).Recently, Davies & Plez (2021) explored the extension of the atmospheres close to the stellar surface at radii smaller than the inner dust shells in the optical and near-IR.Adding the influence of a stellar wind in the MARCS model atmospheres (Gustafsson et al. 2008), Davies & Plez (2021) ex-panded the atmosphere up to several stellar radii.Their results naturally explained the presence of mid-IR excess, as well as the mismatch between temperatures derived from the optical and the IR (Levesque et al. 2005;Davies et al. 2013;González-Torà et al. 2021).They also reproduced many of the features obtained by addition of a static molecular shell (or 'MOLsphere').In short, the work by Davies & Plez (2021) opened a new window to explore the mass-loss rates of cool massive stars.
So far these models have been constrained by comparison to stellar spectra only.The spectral computation shows the flux integrated over the stellar disk and misses the spatially resolved information.Therefore, if we want a detailed way to study spatially extended stellar atmospheres, we need to use interferometric data.Interferometry uses an array of telescopes to increase the angular resolution of the observations.It is a very powerful tool to study the topography of extended atmospheres in detail, and has been used widely both for RSGs (e.g., Arroyo-Torres et al. 2013;Wittkowski et al. 2012;Climent et al. 2020;Chiavassa et al. 2022) and Miras (e.g., Wittkowski et al. 2018;Kravchenko et al. 2020).As a consequence, interferometry represents a stronger test for models.
In this work, we employ the approach by Davies & Plez (2021) to extend atmospheres for both spectral and interferometric data.We compute the synthetic visibilities of the models for different mass-loss rates ( Ṁ) in the wavelength range of the Very Large Telescope Interferometer (VLTI) instruments (1.8−5 µm), which goes from near-IR to mid-IR.We explore the robustness of the model, and present a case study for VLTI/AMBER (Petrov et al. 2007) data for the RSGs HD 95687 and V602 Car from the sample by Arroyo-Torres et al. (2015).We have chosen these two targets because they are examples of two stars with different masses and therefore different Ṁ, as well as two distinct luminosities.
This paper is organised as follows: Section 2 describes the model used for the present study.The theoretical results from our spatially resolved model atmosphere are presented in Section 3, followed by the case study for HD 95687 and V602 Car in Section 4. Lastly, we conclude in Section 5.

Methods
In absence of models that self-consistently explain winds of RSGs, we add a stellar wind with a constant Ṁ to an initial MARCS model (Gustafsson et al. 2008), following the method by Davies & Plez (2021).These models are then used to calculate both the synthetic spectra and the squared visibility amplitudes (|V| 2 ).This is described in detail in the following sections.

Models
We started with a MARCS model atmosphere.This code assumes local thermodynamic equilibrium (LTE), hydrostatic equilibrium, and spherical symmetry.We defined a radius grid for the model, allowing to contain a more extended stratification up to ∼ 8.5 R , where R is defined as the radius where the Rosseland opacity τ Ross = 2/3.Moreover, for simplicity, we assumed that: -The wind is in LTE.A discussion to this assumption can be found in Davies & Plez (2021).-The model is 1D, so we assume spherical symmetry.
To determine the outermost density, we use the mass continuity expression, Ṁ = 4πr 2 ρ(r)v(r) (1) where ρ and v are the density and velocity as a function of the stellar radial coordinate r, respectively.The wind density ρ wind (r) has the shape proposed by Harper et al. (2001): where R max is the arbitrary outer-most radius of the model, in our case 8.5 R .The β and γ parameters define the smoothness of the extended wind region and were initially set in the semi-empirical 1D model of α Ori by Harper et al. (2001): β Harp = −1.10 and γ Harp = 0.45.In Figure 1 we show what happens when changing the γ and β parameters that define the density profile.The variations of β mostly influence the smoothness of the density profile close to the stellar surface, while the variations of γ influence the full density profile to upper or lower values.We will discuss the implications on spectra and interferometric visibility values in Section 3.2.1.
The velocity profile is found assuming a fiducial wind limit of v ∞ = 25 ± 5 km/s, that is the value matched to Richards & Yates (1998); van Loon et al. (2005); Beasor & Davies (2017), and Equation 1.We assume no velocity gradient since the acceleration region is shallow and v ∞ is due to turbulent motions (see Davies & Plez 2021).The model is sensitive to the density ρ, meaning that Ṁ and v are degenerate with one another.
For the temperature profile we first used simple radiative transfer equilibrium (R.E.), defined as: where T(τ Ross =2/3) and R are the temperature and radius at the bottom of the photosphere, and T (r wind ) and r wind are the temperature and radius of the wind extension, respectively.This will result in a smoothly decreasing temperature profile for the extended atmosphere.
In addition, Davies & Plez (2021) defined a different temperature profile in their semi-empirical 1D model of α Ori by Harper et al. (2001), based on spatially-resolved radio continuum data.The main characteristic of this profile is a temperature inversion in the chromosphere of the star, that peaks at ∼ 1.4 R , and decreases again.ALMA and VLA observations of the RSGs Antares and Betelgeuse by Lim et al. (1998);O'Gorman et al. (2017O'Gorman et al. ( , 2020) ) confirm the presence of such a lukewarm chromospheric temperature inversion, peaking at a radius of 1.3 − 1.5 R with a peak temperature of ∼3800 K.However, Lim et al. (1998) pointed out that optical and ultraviolet chromospheric signatures required higher temperatures of ∼5000 K at similar radii (Uitenbroek et al. 1996).Conversely, modelling of spectroscopic and interferometric data of the CO MOLsphere derived gas temperatures of only 2000 K at 1.2 − 1.4 R (Ohnaka et al. 2013).O'Gorman et al. ( 2020) suggested these components co-exist in different structures at similar radii in an inhomogeneous atmosphere, and are spatially unresolved by current measurements.Observations at different wavelengths may then be sensitive to different such structures.
Following Davies & Plez (2021), we include in our model setup either a temperature profile in R.E., which may be more relevant for observations of the near-IR MOLsphere, or a temperature profile with a chromospheric temperature inversion which may be more relevant for chromospheric signatures in the optical or UV or for radio continuum observations.
Once the density, temperature and velocity profiles are defined, we re-sample the model to a constant logarithmic optical depth sampling ∆ log(τ), and we use 0.01 < ∆ log(τ) < 0.05.The reasons for this re-sampling are explained in Davies & Plez (2021): if the grid is too finely sampled, rounding errors can occur leading to numerical difficulties.On the other hand, if the sampling is too coarse, the τ λ = 2/3 surface is poorly resolved for strong absorption lines.
Finally, we define the outer boundary of the model where the local temperature is < 800 K, which is reached at ∼ 8.5 R .Below this temperature our code is unable to reliably converge the molecular equilibrium.In addition some species would be depleted to dust grains.
Figure 2 shows the density, temperature, velocity and Rosseland opacity profiles for the example of an extended model with log L/L = 4.8, T eff = 3500 K, log g = 0.0, M = 15 M and Ṁ = 10 −5.5 M /yr.

Computation of model intensities
We computed both the spectra and the intensity profiles with respect to µ, where µ = cos θ, with θ being the angle between the radial direction and the emergent ray, and cos θ = 1 corresponding to the intensity at the centre of the disk.For this we used the radiative transfer code turbospectrum v19.1 (Plez 2012).Setting a wavelength range from 1.8 µm to 5.0 µm with a step of 0.1Å, to explore the spectral range of the following instruments at the VLTI: -GRAVITY (Gravity Collaboration et al. 2017) for the Kband 1.8 − 2.5 µm.-MATISSE (Lopez et al. 2022) for the L (3.2 < λ < 3.9 µm) and M-bands (4.5 < λ < 5 µm).We did not use the N-band (8 < λ < 13 µm) because it is dominated by dust emission.
For the spectral synthesis, we included a list of atomic and molecular data.Chemical equilibrium is solved for 92 atoms and their first 2 ions, including Fe, Ca, Si and Ti, and molecular data for CO, TiO, H 2 O, OH, CN and SiO is included, among about 600 species.
The spectra and visibilities can be convolved to any spectral resolution used by GRAVITY and MATISSE, or instruments at other interferometers.In this work, we show as an example the results convolved to match the HIGH spectral resolution of MA-TISSE: R = 10001 .

Computation of model interferometric visibilities
To compute the visibility from the intensity profile, we used the following Hankel transform as in Davis et al. (2000), where V Model is the visibility of our model, S λ the instrument sensitivity curve, I µ λ the computed intensities with respect to µ from Turbospectrum, J 0 is the 0th order Bessel function, θ Model is the angular diameter of the outermost layer of the model, and B is the baseline of the observation.The V Model is then normalised with respect to the total flux.To estimate θ model , we used the relation with the Rosseland angular diameter θ Ross , found in Davis et al. (2000) and Wittkowski et al. (2004), where R(τ Ross = 2/3) is the radius of the star at τ Ross = 2/3, defined as the photospheric layer, and R max is the outer most radius of our model.
We use the definition in Wittkowski et al. (2017) to scale the final visibility of the model as where A allows for the attribution of a fraction of the flux to an over-resolved circumstellar component (Arroyo-Torres et al. 2013), and V Model (θ Ross ) is the model visibility computed using Equation 4 with an associated Rosseland angular diameter θ Ross from Equation 5.

Base model
We compute the spectra, intensities and |V| 2 for a base model of  2001) and the wind limit v ∞ = 25 km/s.The temperature profile is initially set to simple R.E., as we are interested in the near-IR K-band MOLsphere (cf.Section 2.1).However, we also look at the effect of a chromospheric temperature inversion in Section 3.2.2.We use mass-loss rates of Ṁ = 10 −4 , 10 −5 , 10 −6 and 10 −7 M /yr, and a simple MARCS model without any wind.As an example, we simulate a star with θ Ross = 3 mas, a baseline of B = 60 m and without any additional over-resolved component, i.e.A = 1.
We observe an extension in all cases except for the MARCS model without the addition of a wind.The CO lines and the SiO lines seem to have the most prominent presence throughout the extended atmosphere (highest intensity compared to water or the continuum).-The visibility spectra are indicative of extended layers of water vapour (centred at 2.0 µm) for high Ṁ ( Ṁ = 10 −5 M /yr), while the flux spectra are less sensitive.This water features are present in the observations of RSGs (e.g., Arroyo-Torres et al. 2015).
-Checking the |V| 2 , we are not able to reproduce some atomic lines in the 2.10 µm < λ < 2.30 µm region for mass-loss rates Ṁ < 10 −4 M /yr, which are the most sensitive to the stratification very close to the stellar surface (Kravchenko et al. 2020).

Density profile
So far we have assumed the density parameters (Equation 2  Furthermore, in Figure 8 we show the spectra in the optical TiO region (λ = 0.5 − 0.75 µm) to see the effect of changing the β parameter in the TiO lines.Again, changing γ parameter produces a very similar plot.A discussion concerning the general effect of this semi-empirical model on the TiO bands can be found in Davies & Plez (2021).Briefly, an increase of the massloss rate will cause the TiO absorption lines to deepen, shifting the star to later spectral types (e.g., for a zero-wind model of spectral type M0, if we apply Ṁ = 10 −6 , Ṁ = 10 −5.5 and Ṁ = 10 −5 M /yr in our model the star will be classified as M1, M2 and >M5, respectively).Changing the density profile parameters with a fixed Ṁ will affect the TiO bands, also shifting the stellar classification slightly to later spectral types as we deepen the TiO bands (Figure 8).
On the other hand, the TiO bands may be more sensitive to a higher chromospheric temperature component than the molecular layers in the near-IR, which may cause the TiO lines to be less deep (cf.Section 2.1).When comparing both temperature profiles in Figures 9 and 10, we see the main difference in the CO lines: in the K-band region λ = 2.29 − 2.7 µm the CO is in emission when using the temperature inversion profile (as also predicted by O'Gorman et al. 2020), while for R.E. it remains in absorption even for very high Ṁ such as Ṁ = 10 −4 M /yr.This difference can be seen in the lower panels of Figures 9 and 10, where the R.E.shows less extension in the CO region as compared with the temperature inversion.For lower mass-loss rates ( Ṁ = 10 −6 M /yr, Figure 10), it gets harder to see the difference in both profiles.The only region that seems to make a difference is the M-band (4.5 < λ < 5 µm), where we see emission in the temperature inversion case.

Temperature profile
Observations of CO lines in the K-band generally show CO in absorption, even for an extreme case such as the RSG VY CMa (Wittkowski et al. 2012), confirming that the lower temperature based on R.E. is better suited to describe the CO MOLsphere than the higher temperature components of the chromospheric temperature inversion (cf.Section 2.1).
To sum up this Section 3, the MARCS+wind model shows significant atmospheric extension in all wavelengths compared to a simple MARCS model (Figure 3).Such an extension has been observed, but it has so far not been reproduced by current models (Arroyo-Torres et al. 2013).As we increase the Ṁ, for a R.E.temperature profile the CO, SiO and water remain relatively unchanged in the spectra (Figure 4), while for the |V| 2 we G. González-Torà et al.: The effect of winds on atmospheric layers of red supergiants  see a larger extension in all cases (Figure 6).The R.E.seems to better reproduce the spectra than the temperature inversion as we do not observe the CO in emission in our case studies (see Section 4) nor other previously published data (e.g., Wittkowski et al. 2012;Arroyo-Torres et al. 2013, 2015).Changing the β and γ parameters in Equation 2 deepens the water features in the |V| 2 .Ṁ and masses.In addition, the data is readily available, and these are two well studied RSGs whose fundamental parameters (e.g., T eff , log g, θ Ross , log L/L ) are well known.The data were taken using the AMBER medium-resolution mode (R ∼ 1500) in the K − 2.1 µm and K − 2.3 µm bands.
The parameters used for each initial MARCS model are shown in Table 1, following Arroyo-Torres et al. (2015).HD 95687 is characterised by a smaller luminosity, mass and radius than V602 Car.In addition, HD 95687 shows a weaker atmospheric extension than V602 Car.
To estimate both θ Ross and A, we compute for each studied model the |V| 2 as a function of their spatial frequency B/λ o where λ o corresponds to the continuum region 2.23 − 2.27 µm.
We have compared the model and data |V| 2 , and found the best fitting θ Ross and A by means of a χ 2 minimisation.Both models with a Ṁ = 10 −4.0 M /yr.We can see that the main differences are in the regions λ = 2.35 − 3.0 µm and λ > 4.0 µm due to CO, where for the Harper+01 case the molecular lines appear more strongly in emission, and the visibilities are decreased due to a larger apparent diameter of the star.
For the temperature profile we use R.E., since the temperature inversion would either show depleted CO lines (λ = 2.29 − 2.7 µm) for the spectra, which do not match with the observations, or not enough extension for the |V| 2 .Therefore it is not possible to find a model with the temperature inversion profile that fits both spectra and |V| 2 simultaneously.
We estimate for our best fit models a θ Ross = 5.35 ± 0.7 mas and 2.9±0.8mas, and A = 1.0±0.14 and 1.0±0.08 for V602 Car and HD 95687, respectively.The errors in θ Ross and A are derived by the minimum values in the 68% dispersion contours of the χ 2 fit, that for 2 degrees of freedom corresponds to χ 2 < χ 2 min + 2.3 Table 1.Parameters of the MARCS models used for the analysis of each RSG.From left to right: luminosity log L/L , effective temperature T eff , surface gravity log g, metallicity [Z], microturbulence ξ, mass as in Arroyo-Torres et al. (2015).The last column shows the radius of the star in the photosphere in solar units (defined at τ Ross = 2/3).In this work we show that, when adding a wind to a MARCS model, we can now qualitatively fit the spectra and |V| 2 .This is something that current existing models are unable to do.This initial fit can be further improved by modifying the inner wind density profile.Figures 13 and 14 show both the spectra and |V| 2 of the new fit changing the density parameters in comparison with the data and the initial MARCS+wind fit for HD 95687 and V602 Car, respectively.The main difference between both MARCS+wind models can be found in the |V| 2 water region at λ = 2.29 − 2.5 µm, where the new γ and β fit better.

RSG
We notice that, although our best-fit model for HD 95687 can accurately reproduce both flux and |V| 2 , this succeeds for V602 Car to a lower extent: the flux is well reproduced in Figure 14, but the |V| 2 is still missing some extension, especially in the region λ = 2.29 − 2.5 µm.As mentioned, this region not only includes CO but also the presence of water.A possible explanation for this mismatch could be that for increasing Ṁ, the models still fail to reproduce the extension of the water or CO layers.Another possibility is that since our model neglects velocity gradients, it underestimates the equivalent widths of lines.Broader and stronger lines would help to increase the apparent stellar extension at those wavelengths.

Summary and conclusion
We present 1D modelling for the extended atmospheres of RSGs based on simple R.E. and chromospheric temperature inversion by Harper et al. (2001), and compute synthetic flux spectra and synthetic interferometric visibility spectra.When comparing our models to a simple MARCS or PHOENIX models, our synthetic |V| 2 shows a stronger atmospheric extension and can fit for the first time the observed extension in the case studies.
Regarding the temperature profile, we find that the R.E.reproduces the spectra better than the chromospheric temperature inversion, since we do not observe any emission in the CO bands, that are the result of models based on a temperature inversion.The possible reason that R.E.fits better than the temperature inversion, even though RSGs are known to have a chromosphere, could be: on the one hand, the presence of different spatial cells with different temperatures in the hot luke-warm chromospheres of RSGs (O'Gorman et al. 2020).On the other hand, we do not know the effect that the dust could make where T > 800 K, although we expect it to be small.Moreover, localised gaseous ejections, related to magnetic fields and surface activity were recently suggested as a major contributor to mass loss from RSGs (Humphreys & Jones 2022;Andrews et al. 2022;López Ariste et al. 2022).To explore this effect in detail, we would need to use 3D models, which is out of the scope of this paper.Our 1D modelling approach relies on an azimuthally averaged stratification, which is a good approximation for many aspects, but may not reproduce some of the observed features.
When compared to the observations, we obtain a mass-loss rate that is in accordance with typical mass-loss prescriptions (e.g., de Jager et al. 1988;Schröder & Cuntz 2005;Beasor et al. 2020).However, in order to fit both the water and CO extensions simultaneously, the density shape should be steeper close to the surface of the star than previously expected by Harper et al. (2001).
Most importantly, we are able to reproduce the |V| 2 extension of the case studies.Simple stellar atmosphere models such as MARCS do not show extension at all.However, the description very close to the stellar surface may not be optimum yet, as we are not able to reproduce some atomic lines in the 2.10 µm < λ < 2.30 µm |V| 2 region, which are the most sensitive to the stratification very close to the stellar surface.This is the first extended atmosphere model to our knowledge that can reproduce in great detail both the spectra and |V| 2 simultaneously.Therefore, we have shown the immense potential of this semi-empirical model of MARCS+wind, not only to match the spectral features without the need of dusty shells, but also the visibilities obtained by interferometric means.
In the future, we want to compare this model with data from a wider wavelength range, to see the full effects in higher wavelengths such as the L or M-bands.

Fig. 1 .
Fig. 1.Left: The different density profiles for γ = 0.45 and variations of β with β = −1.10(blue dots), −1.35 (green triangles) and −1.60 (red squares), we see that as we decrease β, the wind density near R/R − 0.99 = 10 −1 increases and gets steeper.Right: The different density profiles for β = −1.10 and the variations of γ with γ = 0.45 (blue dots), 0.25 (black triangles) and 0.05 (yellow squares).In this case, as we increase γ the slope of the profile remains the same, but the values of the wind density gradually increase.The model used with log L/L = 4.8, T eff = 3500 K, log g = 0.0, M = 15 M and Ṁ = 10 −5.5 M /yr, corresponds to the stellar parameters of HD 95687.

Fig. 2 .
Fig. 2. From top to bottom: The extended profiles for the density (blue dots), temperature (blue squares for temperature inversion by Harper et al. (2001) a.k.a "Harper", and red circles for simple radiative equilibrium), wind velocity (purple) and Rosseland optical depth (green).The extended model with log L/L = 4.8, T eff = 3500 K, log g = 0.0, M = 15 M , Ṁ = 10 −5.5 M /yr and R max = 8.5 R , corresponds to the stellar parameters of HD 95687.
R and R max = 8.5 R , corresponding to a RSG similar to HD 95687 (Arroyo-Torres et al. 2015).The density parameters in Equation 2 are β Harp = −1.10 and γ Harp = 0.45 as in Harper et al. ( Figures 4,5 and 6  show the spectra, the normalised spectra to the continuum, and squared visibility amplitudes (|V| 2 ) computed from our base model with the different Ṁ, from highest Ṁ = 10 −4 M /yr, to the simple MARCS model without wind.When comparing the spectra and |V| 2 in the Figures 4, 5 and 6, there are several things to notice:-In Figures4 and 5, the spectral signatures of CO in the wavelength range of λ = 2.29 − 2.7 µm (K-band) do not strongly depend on the Ṁ, as they remain relatively unchanged.Only at high Ṁ and high resolution spectra the low excitation lines start becoming stronger, as predicted by Tsuji (1988).-In the region λ > 4.0 µm (L and M-bands), the spectra shows the presence of SiO lines at wavelengths up to ∼ 4.3 µm, which seem to remain in absorption up to high Ṁ (∼ 10 −4 M /yr).At 4.3 µm we observe the presence of CO as we increase the Ṁ.The CO lines in the M-band are observed in emission already at low Ṁ (starting at Ṁ = 10 −6 M /yr), while the CO at the K-band remains in absorption.-In the |V| 2 (Figure 6), the most important thing to notice is that the extended molecular layers, mostly of the CO lines in λ = 2.29−2.7 µm (K-band), are seen in the extended models.This extension was not reproduced earlier with MARCS or PHOENIX (Hauschildt & Baron 1999) models alone.-The CO extension increases with increased Ṁ (in all K, L and M-bands).The atmospheric extension is best observed in the M-band, as we see a drop of |V| 2 already for a low
) from Harper et al. (2001): β Harp = −1.10 and γ Harp = 0.45. Figure 7 shows the spectra and |V| 2 for the K-band for the different β values defined in Figure 1.We did not include the variations

Fig. 4 .
Fig.4.Normalised model spectra with respect to the mean flux for Ṁ = 10 −4 (red), 10 −5 (green), 10 −6 (blue) and 10 −7 M /yr (purple).We also plot the spectra and visibilities based on the MARCS model without a wind (orange).This is the case for R.E.We have convolved the results with the spectral resolution of R = 1000.The main differences within models can be seen in the water, CO and SiO molecular bands, the corresponding wavelength regions are highlighted in light grey.

Figures 9
Figures 9 and 10 show the spectra and |V| 2 for our two temperature stratifications defined in Section 2, R.E and temperature inversion, with Ṁ = 10 −4 M /yr and Ṁ = 10 −6 M /yr, respectively.When comparing both temperature profiles in Figures 9 and 10, we see the main difference in the CO lines: in the K-band region λ = 2.29 − 2.7 µm the CO is in emission when using the temperature inversion profile (as also predicted by O'Gorman et al. 2020), while for R.E. it remains in absorption even for very high Ṁ such as Ṁ = 10 −4 M /yr.This difference can be seen in the lower panels of Figures9 and 10, where the R.E.shows less extension in the CO region as compared with the temperature inversion.For lower mass-loss rates ( Ṁ = 10 −6 M /yr, Figure10), it gets harder to see the difference in both profiles.The only region that seems to make a difference is the M-band (4.5 < λ < 5 µm), where we see emission in the temperature inversion case.Observations of CO lines in the K-band generally show CO in absorption, even for an extreme case such as the RSG VY CMa(Wittkowski et al. 2012), confirming that the lower temperature based on R.E. is better suited to describe the CO MOLsphere than the higher temperature components of the chromospheric temperature inversion (cf.Section 2.1).

Fig. 5 .
Fig. 5. Same as Figure 4 but for the flux normalised to the continuum.

Fig. 7 .
Fig. 7. Normalised model spectra (upper panel) and |V| 2 (lower panel) for a fixed γ = 0.45 and the different β = −1.10(blue), −1.35 (green) and −1.60 (red) in the K-band.The baseline used is B = 63.8 m, corresponding to the case study of HD 95687 in Section 4. We can see that as the β gets lower, the water features in λ = 2.35−2.5 µm become slightly more prominent.
Figures 11 and 12 show the MARCS model fit to the data of Arroyo-Torres et al. (2015) and our initial MARCS+wind model fit with β Harp = −1.10 and γ Harp = 0.45, compared to the data of HD 95687 and V602 Car, respectively.For our model fit, we check both the spectra and |V| 2 .We use a range of mass-loss rates of −7 < log Ṁ/M < −4 with a grid spacing of ∆ Ṁ/M = 0.25.

Fig. 9 .
Fig. 9. Normalised model spectra (upper panel) and |V| 2 (lower panel) for our base model with simple radiative equilibrium (red) and the chromospheric temperature inversion profile by Harper et al. (2001) (blue).Both models with a Ṁ = 10 −4.0 M /yr.We can see that the main differences are in the regions λ = 2.35 − 3.0 µm and λ > 4.0 µm due to CO, where for the Harper+01 case the molecular lines appear more strongly in emission, and the visibilities are decreased due to a larger apparent diameter of the star.

Fig. 10 .
Fig.10.Same as Figure9, but with a Ṁ = 10 −6.0 M /yr.In this case, the main differences are only in the λ > 4.0 µm region due to CO, where for the Harper+01 case the molecular lines appear more strongly in emission, and the visibilities are decreased due to a larger apparent diameter of the star.

Fig. 13 .
Fig. 13.Upper left: Normalised flux for the RSG HD 95687 (grey), as observed with VLT/AMBER for the K − 2.1 µm bands.Our initial best-fit model with β = −1.10 and γ = 0.45 is shown in red, corresponding to the parameters of the density profile by Harper et al. (2001).The final best-fit model with β = −1.60 and γ = 0.05 is shown in green.Upper right: Same as the upper left panel but for the K − 2.3 µm band.Lower left: Same as the upper left panel but for the |V| 2 and a baseline of B = 60.8 m.Lower right: Same as the lower right but for the K − 2.3 µm band and a baseline of B = 63.2 m.The best fit model for β and γ can reproduce the water features in λ = 2.29 − 2.5 µm with better accuracy than our initial best fit.