The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE

Context. VX Sgr is a cool, evolved, and luminous red star whose stellar parameters are difﬁcult to determine, which affects its classiﬁcation. Aims. We aim to spatially resolve the photospheric extent as well as the circumstellar environment. Methods. We used interferometric observations obtained with the MATISSE instrument in the L (3–4 µ m), M (4.5–5 µ m), and N (8– 13 µ m) bands. We reconstructed monochromatic images using the MIRA software. We used 3D radiation-hydrodynamics simulations carried out with CO 5 BOLD and a uniform disc model to estimate the apparent diameter and interpret the stellar surface structures. Moreover, we employed the radiative transfer codes O PTIM 3D and R ADMC 3D to compute the spectral energy distribution for the L , M , and N bands, respectively. Results. MATISSE observations unveil, for the ﬁrst time, the morphology of VX Sgr across the L , M , and N bands. The reconstructed images show a complex morphology with brighter areas whose characteristics depend on the wavelength probed. We measured the angular diameter as a function of the wavelength and showed that the photospheric extent in the L and M bands depends on the opacity through the atmosphere. In addition to this, we also concluded that the observed photospheric inhomogeneities can be interpreted as convection-related surface structures. The comparison in the N band yielded a qualitative agreement between the N -band spectrum and simple dust radiative transfer simulations. However, it is not possible to ﬁrmly conclude on the interpretation of the current data because of the difﬁculty in constraing the model parameters using the limited accuracy of our absolute ﬂux calibration. Conclusions. MATISSE observations and the derived reconstructed images unveil the appearance of VX Sgr’s stellar surface and circumstellar environment across a very large spectral domain for the ﬁrst time.

These mass ejections cause a pronounced extended atmosphere in the K band (Monnier et al. 2004;Chiavassa et al. 2010) where interferometric observations revealed evidence of departure from circular symmetry, which could be caused by the opacity (H 2 O and CO) throughout the atmosphere (Chiavassa et al. 2010). In addition to this, its circumstellar environment is detected with SPHERE (Scicluna et al. 2020), where the authors reveal the presence of ejecta confined to non-spherical or clumpy configurations. Moreover, this environment is expected to be characterised by the presence of alumina, metallic iron, and Fecontaining silicates (Gail et al. 2020). Observations in the radio domain, with H 2 O and SiO masers (Yoon et al. 2018) or OH (Reid et al. 1977), showed larger expansion layers with respect to the typical semi-regular variables. Radio observations also made it possible to determine precise distances, even if the measures do not entirely overlap among the different works. Using the proper motions of the SiO maser at 43 GHz, Su et al. (2018) reported a distance of 1.10 ± 0.11 kpc, and Chen et al. (2007) found 1.57 ± 0.27 kpc.
Using a 22 GHz maser's emission of H 2 O, Xu et al. (2018) measured a distance of 1.56 +0. 11 −0.10 kpc. Other authors estimated the distance with different methods: Humphreys et al. (1972) placed VX Sgr in the vicinity of the Sgr OB1 cluster at 1.7 kpc;Humphreys (1974) calculated a photometric distance of 0.8 kpc; and Hipparcos parallaxes gave a value of 0.262 +0.655 −0.109 kpc (van Leeuwen 2007). It should be noted that the latter two measurements could be highly contaminated by convection-related stellar surface structures (Chiavassa et al. 2011b. As a consequence of these uncertainties, the radius of VX Sgr ranges between 1000 and 2000 R (Monnier et al. 2004;Chiavassa et al. 2010;Xu et al. 2018) as well as the luminosity between 1.95 × 10 5 and 5.5 × 10 5 L ( Levesque et al. 2005;Chiavassa et al. 2010;Xu et al. 2018), depending on the distance assumed.
Due to the limitations on precisely determining the stellar parameters, the classification of VX Sgr is still under debate. On one side, if one assumes the distances obtained in the radio, the star appears to be in agreement to what is expected for standard red supergiant (RSG) stars (Arroyo-Torres et al. 2015); but, on the other hand, Tabernero et al. (2021) provided new insights into its luminosity, evolutionary stage, and its pulsation period, which seem to point towards some sort of an extreme AGB star.
In this context, long-baseline infrared interferometric observations play an important role because they allow us to explore and characterise the geometrical extent of the stellar photosphere as well as its circumstellar environment. However, from an observational point of view, RSGs may bear similarities with AGBs (e.g. Arroyo-Torres et al. 2015), even though they largely differ in terms of mass, effective temperature, surface gravity, and photometric variability (Levesque et al. 2005;Höfner & Olofsson 2018).
In this work, we report observations of VX Sgr obtained with the MATISSE instrument at VLTI and aim to explore the properties of the stellar surface and circumstellar environment in the wavelength bands: L, M, and N.
In each of the L, M, and N bands, MATISSE uses a multiaxial all-in-one beam combination, that is the four separated telescope beams are focussed onto the detector by a common lens and produce six dispersed fringe patterns that are mixed into one single focal spot (an airy disc), representing the socalled interferometric channel. In the L and M bands, the total incoming flux is split between the interferometric channel and the four photometric channels (SiPhot mode). This mode allows us to simultaneously perform the source photometry and thus to properly calibrate the visibilities. In the N band, the source photometry is obtained after the fringes recording by closing alternatively the MATISSE shutters to feed the four photometric channels. That allows us to improve the sensitivity by sending 100% of the total incoming flux in the different channels (high sens mode).
MATISSE operates in the mid-infrared, where thermal background emission is significant. As a consequence, two additional modulation steps have been added to calibrate the data: the first one is chopping modulation (observing the target and an empty area of the sky alternatively). It is a mandatory step to extract all source photometries in N band while allowing reliable photometry measurements for faint targets in the L and M bands. The second modulation is an optical path difference (OPD) modulation, aimed at removing the low-frequency peak contamination, containing the contribution from the thermal background, to the fringe peaks (in the Fourier transform of the interferograms). The different modes of the instrument (telescope chopping, detector integrations, OPD modulation) are synchronised all together to ensure minimum overheads when observing. Finally, a beam modulation (or beam commutation) is applied to the observations by swapping two-by-two beams in a four-step sequence. This last feature is aimed at removing detection artifacts in the phases measured by the instrument (differential phases and closure phases). These artifacts are instrument phase residuals located in the optical train between the BCD and the detector (Lopez et al. 2022). The BCD (Beam Commuting Device) is a module used to calibrate closure and differential phases.

Data acquisition and reduction
We observed VX Sgr at its visible photometric minimum over three nights (Table 1) with the MATISSE instrument with the aim of reconstructing images across the different bands. The data collection spanned 2.5 months: the compact and medium configurations within 15 days in June, while the large configuration could only be observed late at the end of August 2019. The calibrators used are reported in Table 2, and the UV coverage is displayed in Fig. 1. We used the low spectral resolution mode (R = 35) in the L, M, and N bands.
For the data reduction, we applied the MATISSE data reduction software (matisse-drs, Millour et al. 2016). It adopts the classical Fourier transforms scheme (Perrin 2003), with photometric calibration in the L and N bands (Coudé du Foresto et al. 1997). Chopping correction (subtracting sky frames from target frames) and OPD demodulation are applied before summing spectral densities 1 , bispectra 2 , and interspectra 3 , leading to estimates of the spectrointerferometric observables; spectrally dispersed in squared visibility, closure phase, and differential phase, similarly to what was previously done with the AMBER instrument (Petrov et al. 2007;Tatulli et al. 2007). A MATISSE observation sequence is composed of twelve one-minute-long exposures with simultaneous interferometric and photometric data in the L and M bands and four interferometric exposures followed by eight photometric exposures in the N band. The first four exposures are taken without chopping while the following eight exposures are chopped. Each exposure is taken in one of the four configurations of the BCD. Overall, this leads to four  Table 2. reduced OIFITS2 (Duvert et al. 2017) files per pointed target (one for each BCD configuration) containing, per exposure, six dispersed squared visibilities, differential visibilities phases, and three independent dispersed closure phases (out of four in total). The interferometric calibration is then applied to the data: the targets used as calibrators (CAL) are corrected from their expected observables given the predicted value of their angular diameters, and each science target (SCI) is corrected from the transfer function thus estimated. The calibrated OIFITS2 files of the four BCD configurations are merged to produce a single calibrated OIFITS2 file for each CAL-SCI pair. Figure 2 shows a representative example of the L, M, and N bands' calibrated squared visibilities and closure phases taken with the small AT configuration (A0-B2-D0-C1) on 15 June 2019. The data taken at low spatial frequencies are of good quality, except for the ranges of the atmospheric absorption, and they already show a spatially resolved emission in all the bands, with strong brightness asymmetries in the N band indicated by the significant non-null closure phases. Figure 3 displays the calibrated MATISSE spectrum of VX Sgr for one particular night (25 June 2019) overplotted on the telluric spectrum at the observing site. For each observation block of VX Sgr, the absolute flux calibration was performed by dividing the raw spectrum of VX Sgr by the raw spectrum of the associated calibrator and then multiplying this ratio by the calibrator synthetic spectrum. For the night shown in the figure, the synthetic spectrum of the calibrator δ Sgr (Table 2) comes from the PHOENIX grid (ACES-AGSS-COND, Husser et al. 2013). According to Fig. 3, the data quality is good, while the telluric line dominates in specific spectral regions: lower than ∼3.3 µm, between 4.0 and 4.7 µm, and between 9.3 and 9.9 µm. In the following, we include all spectral regions in our work, but we highlight the telluric spectrum at the observing site.  Closure phase (degree)

The observed stellar surface of VX Sgr
We performed image reconstruction for the observed data using the MIRA software package developed by Thiébaut (2008). This package is written in the scientific language yorick 4 . It uses interferometric data in the form of OIFITS2 files (Pauls et al. 2005). The reconstructed image is compared to the visibility and closure phase data by means of Fourier transforms (FT): a FT is first computed on the image, and visibilities are computed as the FT amplitude at the observation spatial frequencies ( f i j = B i j λ ), B i j being the projected baseline between telescopes i and j, and λ the wavelength of observation. On the other hand, phases are computed as the FT phase at f i j . The closure phases are computed as the sum of the three phases at the three observation spatial frequencies ( f i j , f jk , f ki ) corresponding to a telescopes triangle. The pixels values are modulated in intensity to minimise the so-called "cost function", which is the sum of a regularisation term plus data-related terms (χ 2 ). The data terms enforce the agreement of the model image with the measured data (visibilities and closure phases). The regularisation term is a distance minimisation between the reconstructed image and an expected image, called the prior. The data term and the regularisation terms must be 4 https://software.llnl.gov/yorick-doc/ scaled together (thanks to the use of the "hyperparameter" (µ 5 ) to keep a similar magnitude in the characterise process.
We reconstructed images using the observed visibilities and closure phases at each spectral channel individually and assuming a resolution of 0.39 and 0.59 mas pixel −1 and a field of view (FOV) of 256 × 256 pixels and 512 × 512 pixels in the L, M, and N bands, respectively. These represent (i) a total variation regularisation, (ii) a circular Gaussian prior image with a full width at half maximum (FWHM) equal to 0.5 × FOV, (iii) a hyperparameter value of µ = 5 × 10 4 . For each observation block, the absolute flux calibration was performed by dividing the raw spectrum of the science target by the raw spectrum of the associated calibrator. Figure 4 displays the reconstructed images in the L and M bands, while Fig. 5 shows the same in the N band. Figure 6 shows colour composite images across different wavelength ranges and illustrates the much larger structures in the N band compared to shorter wavelengths. The corresponding intensity RMS (i.e. noise maps) for all the images is reported in Figs. A. 1 and A.2. In the L and M bands, the images show a complex morphology with a central brighter area, whose position depends on the wavelength probed. The photocentre deviates from the geometrical centre, as already pointed out for another evolved star using the same reconstruction method (Chiavassa et al. 2020). In addition to this, the stellar apparent diameter is larger at shorter wavelengths, smaller between ∼3.5 and ∼4.0 µm, before getting larger again for longer wavelengths. In the N band, at λ ∼ 8.40 µm where the brightness is about 10% of the central region. Then the contrast increases with wavelengths to ∼30% at ∼11 µm and reaches a maximum of almost 40% at ∼13 µm.

Interpreting the observations in the L and M bands
In this section, we present the simulations and methodology used to make a comparison to the observations.

3D simulations of AGB and RSG stars
We used the radiation-hydrodynamics (RHD) simulations of stellar surface convection computed with the CO 5 BOLD (Freytag et al. 2012) code. The code solves the coupled nonlinear equations of compressible hydrodynamics and non-local radiative energy transfer in the presence of a fixed external spherically symmetric gravitational field on a 3D Cartesian grid. No artificially-induced pulsations are added to the simulations (e.g. by a piston) but they are self-excited. Molecular opacities are taken into account, but radiation transport is treated in a grey approximation, ignoring radiation pressure and dust opacities. Dynamical pressure allows the density to drop much more slowly as a function of the distance than expected for a hydrostatic atmosphere (for more details see Freytag et al. 2017).
Since the uncertainty on the stellar properties of VX Sgr is large (Sect. 1), we used two different simulations: (i) for the AGB stellar type, we used the one from Freytag et al. (2017) with the highest T eff and L , and (ii) for the RSG we used the one presented in detail in Kravchenko et al. (2019). The averaged stellar parameters of the simulations used in this work are reported in Table 3. The temporal variability, used to increase the number of possible matching images, is ensured by the 100 snapshots (5.75 yr, with temporal timestep of ∼20 days) for the AGB simulation and 180 (10.35 yr, with temporal timestep of ∼20 days) for the RSG one. These snapshots are extracted from the relaxed simulated time sequence reported in Table 3.
Afterwards, we employed OPTIM3D (Chiavassa et al. 2009) to compute synthetic intensity maps in the same observed spectral channels as MATISSE and using as input the RHD simulations snapshots. This code takes into account the Doppler shifts caused by the convective motions. The radiative transfer is computed in detail using pre-tabulated extinction coefficients per unit mass, as for the hydrostatic model atmosphere code MARCS (Gustafsson et al. 2008), as functions of temperature, density, and wavelength for the solar composition (Asplund et al. 2009). The temperature and density distributions are optimised to cover the values encountered in the outer layers of the RHD simulations.
The MATISSE observations cover the L, M, and N bands. While the L and M bands are mostly characterised by molecular features, the N-band spectrum may display several dust features (e.g. Paladini et al. 2017), which in turn affect the general shape of the star and its circumstellar environment. On the theoretical side, there are no source terms or dedicated opacities for dust in RHD simulations (Freytag et al. 2017). In the end, we used the RSG and AGB models to compute synthetic images in the L and M bands, and, only for a testing purpose, the AGB model for the N band.

Apparent diameter across wavelengths
We measured the stellar diameters using uniform disk (UD) model and the 3D RHD simulations of Table 3. As a first step, we computed the azimuthally average intensity profiles for all the synthetic maps (i.e. snapshots) of the simulations and for each observed spectral channels. They were constructed using rings regularly spaced in µ (where µ = cos(θ), with θ being the angle between the line of sight and the radial direction). Afterwards, we calculated the temporal averages. Figure 7 shows examples of snapshots for the RSG and AGB simulations across the L and M bands. These profiles are an effective proxy to investigate the overall shape and extension of their atmospheres. The AGB simulation displays flatter and more extended atmospheres as already shown in Wittkowski et al. (2016). Moreover, these profiles are also useful to point out the wavelength dependence of the synthetic images: shorter and medium wavelengths (dark and intermediate green) returns smaller radii while longer wavelengths (light green) larger sizes. This aspect is noticeable in the reconstructed images of Fig. 4.
Eventually, we used the temporally averaged (over the simulation time, Table 3) intensity profiles to obtain the synthetic visibility amplitudes using the Hankel transform and fitted the visibility data separately for each observed spectral channel (this  (Fig. 1). The peak intensity is normalised between [0, 1], and the corresponding noise maps are shown in Fig. A.1. The highest noise value is ∼1.5 × 10 −4 . The contour lines correspond to [0.10, 0.30, 0.70, 0.95, 0.97, 0.99]. East is to the left, and north is up. As reported in Sect. 2.2, the data quality of the images at wavelengths lower than ∼3.3 µm and between 4.0 and 4.7 µm is compromised by the telluric spectrum at the observing site. Notes. The first column shows the simulation name, and the following five columns show the stellar parameters such as the mass (M ), the average luminosity (L ), radius (R ), effective temperature (T eff ), and surface gravity (log g). The different quantities are averaged over spherical shells and epochs (7th column, t avg ). Errors are one standard-deviation fluctuations with respect to the time average (see Chiavassa et al. 2011aChiavassa et al. , 2009). The solar metallicity is assumed. procedure has already been used, e.g. in Chiavassa et al. 2020).
The wavelength-dependent apparent diameters of VX Sgr are reported in Table 4. As a comparison, Chiavassa et al. (2010) measured a photospheric diameter (corresponding to the reference radius at 1.04 µm) of Θ = 8.82 ± 0.50 mas, which is at least two times smaller than the smallest diameter measured in the L band (Θ AGB = 18.27 ± 0.14 mas at 3.29 µm). However, in Chiavassa et al. (2010) the size of VX Sgr was largely increasing towards the K band up to 2.50 µm, where the apparent diameter was close to ∼20 mas. This is very similar to what is found in this work for several spectral channels. In the K band, the authors explained that this increase in size was due to the presence of H 2 O and CO molecules as dominant absorber in the outer layers.   Table 4 with highlighted in green the molecular and dust features that principally contribute to the observed flux. In the L and M bands, the apparent diameter observed depends on the opacity through the atmosphere. Theoretical work on AGB stellar atmospheres shows that, at low optical depth where molecules form, the dynamical pressure dominates over the gas pressure by factor of 5 to 10 (Freytag et al. 2017). This makes possible the levitation of dense material into cool layers and allows the formation of irregular and non-spherical MOLsphere 6 that, in turn, make the stellar size larger. Table 4. Apparent diameters (in mas) in the L, M, and N bands obtained from a uniform disc (Θ UD ) model fitted to the MATISSE squared visibilities and fitting the synthetic RSG (Θ RSG ) and AGB (Θ AGB ) visibilities (computed from from the RHD simulations of Table 3 Notes. The arrows in the first column and the horizontal dashed lines indicate the wavelength interval where the telluric spectrum at the observing site is weak (Fig. 3).  Figure 8 (top panel) displays that the UD, RSG, and AGB diameters show similar behaviours with larger sizes at short wavelengths (where CO and OH are present), then smaller ones between ∼3.25 and ∼4.00 µm, before increasing at longer wavelengths where SiO and CO dominate. In addition to these molecules, H 2 O also contributes across the full spectral range. However, abrupt discontinuities are visible for RSG and UD diameters at (i) ∼3.15 µm and (ii) ∼4.27 µm. The strong variation in UD and RSG at ∼3.15 µm indicates that these models are not adapted to describe the observed data since they underestimate the increase in size with the wavelength. This is not surprising because Arroyo-Torres et al. (2015) already pointed out that in 3D RSG simulations the dynamical pressure is not sufficient to enlarge the stellar atmosphere to the observed sizes in the K band. On the other hand, the very extended atmospheres of the 3D AGB simulation shows a more regular increase of the radius from shorter to longer wavelengths, denoting that this simulation is more adapted to represent the data. It should be noted that the observed data in the spectral range between 4.27 and 4.44 µm are partially contaminated by the presence of noise (Fig. 2) that makes the determination of the visibility fit less precise, and, as a consequence, there is an abrupt change in the diameter values at precisely 4.27 µm.
In the N band, Fig. 8 (bottom panel) shows the increase of the angular diameter with the wavelength as already visible in Figs. 5 and 6. However, this result is purely indicative since either the UD or the AGB simulation lack dedicated dust opacities, where the dust features prevail.

Stellar surface convection in the L and M bands
In this section, we compare the stellar surface structure observed in the aperture synthesis imaging in the L and M bands to the RHD synthetic maps. The main aim is to focus the analysis on the central part of the stellar surface to point out evidence for the presence of convection-related structures.
We proceeded as follows: (i) We rebinned the 3D synthetic maps of both RHD simulations of Table 3 to the resolution of 0.39 mas pixel −1 so that they effectively possess the same FOV and resolution as the reconstructed images (Fig. 4). (ii) We convolved them to the resolution of the interferometric beam ( Fig. 1). (iii) We used the definition of Tremblay et al. (2013) to estimate the intensity contrast as δI rms / I for each snapshot, where I is the intensity. This is performed only for the pixels within the stellar surface with a cut at 70% of the maximum of the ring-averaged intensity profile 7 . This procedure allows us to avoid outer areas dominated by the limb effect. (IV) We finally performed, for each snapshot and at a given wavelength, a minimisation procedure based on the χ 2 snapshot = |C observations −C 3D | 2 σ 2 C observations , where C observations and C 3D are the intensity contrasts for the reconstructed image (observations) and the synthetic one (3D), while σ C observations is the standard deviation of the observed intensity contrast. Afterwards, we computed the reducedχ 2 asχ 2 = 1 N−1 χ 2 snapshot where N is the total number  Table 3. As a reference, few wavelengths are indicated in top right corner. of wavelength channels (N = 30 in the L and M bands at low spectral resolution). In this way, our comparison is focused to matching simultaneously all the observed spectral channels with the same snapshot. This approach aims at having an additional wavelength dependent constraint for the simulation (i.e. the same snapshot should match all the wavelength channels). In the end, the best-matching snapshot is the one with the smallestχ 2 .
In the L and M bands, the best matching RSG snapshot has χ 2 =1.12, and the has worst 29.80. For the AGB's best snapshot, χ 2 =1.98 (worst 580.85). Figure 9 shows the measured intensity contrast for the reconstructed and synthetic images. Figure 10 displays the best-matching snapshots for a particular wavelength channel (4.51 µm) for the two considered RHD simulations of Table 3. In the L and M bands and restricting the comparison to the inner region only, the synthetic images of the AGB simulation are considerably wavelength dependent and, as a consequence, their morphology and photospheric extent noticeably change the apparent angular diameter (Fig. 4). In terms of surface intensity contrast, their agreement with the reconstructed images is slightly worse in the case of the RSG simulation but theχ 2 values are very close. In the end, our 3D RHD simulations return a good agreement with the observations in terms  (Table 4) with the RHD simulations of Table 3  of contrast of surface structures, meaning that they are adequate for interpreting the observed photospheric inhomogeneities as convection-related surface structures.

N-band and tentative analysis with dust radiative transfer modelling
We present in this section a qualitative analysis of the data obtained in the N band. First, we used the MATISSE calibrated flux as described in Sect. 2.2. Figure 11 displays the MATISSE calibrated flux. The calibrators for the night 2019-06-10 (see Table 1) are too faint to allow a proper calibration and the data are extremely noisy. Both nights of 25 June 2019 and 28 June 2019 show similar shapes of the spectrum, but with lower quality for 2019-08-28. However, the mean absolute flux levels between these nights are different. Although there may be concerns about imperfect background subtraction and slight differences in air mass between the source and the calibrator (which are not corrected in the flux calibration process), the main cause of this discrepancy is the PHOENIX stellar synthetic spectra used for the calibrators. In our analysis, these spectra have been extrapolated from the near-IR to the N band and, for the night of 2019-06-25, the theoretical spectrum of the calibrator underestimates the measured SED by a factor of two, while for the night 2019-08-28, there is an overestimation of the averaged flux measured. This partially explains the fact that the spectrum collected on 2019-06-25 is weaker than the one from 2019-08-28. In conclusion, the accuracy of our flux calibration process does not allow us to obtain a consistent absolute calibration in terms of flux level between the different nights. For this reason, we focussed our analysis on interpreting only the overall shape of the spectrum. For this purpose, the nights of both 2019-06-25 and 2019-08-28 are used.
To model the observations in this wavelength range, where the dust contribution is significant, we used version 2.0 of the publicly available radiative transfer code RADMC3D (Dullemond 2012). The code needs to specify a radiation source, and then, from a given input dust density distribution and dust opacities, it first computes the dust temperature distribution (dust continuum radiative transfer, Bjorkman & Wood 2001). Then, at every chosen wavelength, a scattering Monte Carlo run is performed to compute the scattering source function, which is then followed by a ray-tracing to compute the SEDs and the images. We used a spherical grid centred on the star initially with 20 grid points in each direction, with an adaptive mesh refinement scheme: each cell is refined if the dust density (ρ) gradient (ρ max − ρ min )/ρ max > 0.05, ρ min and ρ max are the minimum and maximum density values of 50 random points taken into each cell, respectively. For the stellar flux, we used the 1D hydrostatic synthetic spectra (Lançon et al. 2007) for a star with effective temperature of 3700 K and surface gravity equal to 0.0, as was proposed by Montargès et al. (2021), who analysed the observations of the RSG α Ori. Beside the fact that Arroyo-Torres et al. (2013) found and effective temperature of 3700 K for VX Sgr, we also tried other 1D spectra with effective temperature of 2900 K and 3400 K with worse results than the one at 3700 K.
To model the apparent extension beyond the stellar radius, we used an envelope surrounding the star with three dust species: melilite (Ca 2 Al 2 SiO 7 , Mutschke et al. 1998), corundum (Al 2 O 3 Zeidler et al. 2013), and olivine (MgFeSiO 4 , Jaeger et al. 1994;Dorschner et al. 1995). The choice of these species is based on the precedent findings for RSG stars (Verhoelst et al. 2009). The optical constants were obtained from the Jena database 8 . We set the stellar angular diameter to 8.82 mas (Chiavassa et al. 2010, the near infrared value derived by Chiavassa et al. 2010) to scale the synthetic stellar flux. To account for the extension of the object across the N band ( Fig. 5 and Table 4), we computed SEDs with a flux contribution coming from a diameter of 130 mas. Figure 12 (top) displays an example where we set a grain size, a, ranging between [0.  µm, with a size distribution of n(a) ∝ a −3.5 (Mathis et al. 1977), for all the three species. We set the inner edge of the dust density distribution at 3.5 R to match the reconstructed images in Fig. 5. The dust mass-loss rate is set at 2 × 10 −7 M yr −1 for each of the three species. Moreover, Fig. 12 (bottom) shows several RADMC3D for the olivine with different grain size ranges. The olivine flux decreases when reducing the grain size range and this behaviour is similar for all the other species. The general trend of the MATISSE data, for wavelength longer than ∼8.5 µm, is in qualitative agreement with the presence of olivine and melilite, while corundum does not seem to be visible. At shorter wavelengths, a further contamination of other molecules not included here (e.g. SiO) may cause the discrepancy. Unfortunately, it is not possible to firmly conclude on the interpretation of the current data because of the difficulty in constraining the parameters of RADMC3D models using the limited accuracy of the absolute flux calibration.

Conclusions
Our MATISSE observations unveil the photospheric extent for the first time, as well as the circumstellar environment of VX Sgr across the L, M, and N bands. We reconstructed monochromatic images using the MIRA software. They show a complex morphology with brighter areas, whose characteristics depend on the wavelength probed. In the N band, the central object is still visible at shorter wavelengths, but then the circumstellar brightness grows at longer wavelengths from 10 up to 40% of the central region.
We employed 3D RHD simulations of AGB and RSG stars, carried out with CO 5 BOLD, to calculate synthetic intensity maps in the same observed spectral channels of the MATISSE observations. As a first step, we measured VX Sgr monchromatic angular diameters using a UD model fit and the azimuthally averaged intensity profiles from AGB and RSG simulations applied to the MATISSE interferometric data. The apparent diameter measured in the L and M bands depends on the opacity through the atmosphere: VX Sgr features larger sizes at short wavelengths (where CO and OH are present), then it becomes smaller between ∼3.25 and ∼4.00 µm, before increasing at longer wavelengths where SiO and CO dominate. On top of these molecules, H 2 O also contributes across the full spectral range. The dynamical pressure dominates over the gas pressure at low optical depth, where molecules form (Freytag et al. 2017 Table 3 at 4.57 µm. The intensity is normalised between [0, 1]. The images, 256 × 256 pixels with 0.39 mas pixel −1 , were convolved with the interferometric beam (Fig. 1). The contour lines are the same as in Fig. 4. East is left, while north is up. Bottom row: synthetic intensity maps from the above snapshots non-convolved with the interferometric beam. The peak intensity is normalised between [0, 1]. The RSG snapshotχ 2 = 1.12, while the AGB one isχ 2 = 1.98.  in turn allows the formation of an irregular and non-spherical MOLsphere. This extended photosphere is what is observed by MATISSE. In this context, the comparison with the L and M images denoted also that the photospheric extent predicted in 3D RSG simulation is not sufficient, as already shown in the K band (Arroyo-Torres et al. 2015). The spectrum in the N band may be characterised by the signature of several dust features (Paladini et al. 2017) with a consequence on the brightness and on the general morphology of the star and its circumstellar environment. The angular diameter increases with wavelength. However, our comparison is purely indicative because, where the dust features prevail, the AGB simulation we used does not include source terms or dedicated opacities for dust.
We also compared the stellar surface structure observed in the reconstructed images in the L and M bands to RHD synthetic maps. Our 3D RHD simulations return a good agreement with the observations in terms of contrast of surface structures. We conclude that the photospheric inhomogeneities observed in the reconstructed images could be interpreted as convection-related surface structures.
To reproduce the observations in the N band, we used the publicly available radiative transfer code RADMC3D to calculate the SED of an envelope surrounding the star with three dust species: melilite (Ca 2 Al 2 SiO 7 ), corundum (Al 2 O 3 ), and olivine (MgFeSiO 4 ). The general trend of MATISSE data is in qualitative agreement with the presence of melilite and olivine but not corundum. Nevertheless, it is not possible to firmly conclude on the interpretation of the current data because of the difficulty in A&A 658, A185 (2022) contrasting the model parameters using the limited accuracy of our absolute flux calibration.
Even if the MATISSE observations are extremely useful to characterise and resolve the stellar surface as well as its circumstellar environment, they are not suitable to firmly conclude on the nature of VX Sgr. The reason is twofold: on one hand, the spectral resolving power used here (R = 35) is not sufficient to observe the spectral lines in detail, which could provide crucial information on the velocity fields and atmospheric stratification of the observed object (Chiavassa et al. 2011a); on the other hand, RSG and AGB stars display similarities in interferometric observables (e.g. Chiavassa et al. 2010;Arroyo-Torres et al. 2015;Wittkowski et al. 2019), making their interpretation difficult. Time-dependent, spatially and spectroscopically resolved, multiwavelength simultaneous observations will be key to probing the stellar parameters of VX Sgr and characterising its photometric variability.