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

A. Chiavassa, K. Kravchenko, M. Montargès, F. Millour, A. Matter, B.

Freytag, M. Wittkowski, V. Hocdé, P. Cruzalèbes, F. Allouche, et al.

► To cite this version:

A. Chiavassa, K. Kravchenko, M. Montargès, F. Millour, A. Matter, et al. The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE. Astronomy and Astrophysics - A&A, 2022, 658, pp.A185. 10.1051/0004-6361/202142514. hal-03583566

HAL Id: hal-03583566

https://hal.science/hal-03583566

Submitted on 21 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE*

A. Chiavassa<sup>1,2,3</sup> ◐, K. Kravchenko<sup>4</sup> ◐, M. Montargès<sup>5</sup> ◐, F. Millour<sup>1</sup>, A. Matter<sup>1</sup>, B. Freytag<sup>6</sup>, M. Wittkowski<sup>2</sup> ◐, V. Hocdé<sup>7</sup>, P. Cruzalèbes<sup>1</sup>, F. Allouche<sup>1</sup>, B. Lopez<sup>1</sup>, S. Lagarde<sup>1</sup>, R. G. Petrov<sup>1</sup>, A. Meilland<sup>1</sup>, S. Robbe-Dubois<sup>1</sup>, K.-H. Hofmann<sup>8</sup>, G. Weigelt<sup>8</sup> ◐, P. Berio<sup>1</sup>, P. Bendjoya<sup>1</sup>, F. Bettonvil<sup>9</sup>, A. Domiciano de Souza<sup>1</sup>, M. Heininger<sup>8</sup>, Th. Henning<sup>10</sup>, J. W. Isbell<sup>10</sup>, W. Jaffe<sup>9</sup>, L. Labadie<sup>11</sup>, M. Lehmitz<sup>10</sup>, K. Meisenheimer<sup>10</sup>, A. Soulain<sup>12</sup>, J. Varga<sup>9,13</sup>, J.-C. Augereau<sup>12</sup>, R. van Boekel<sup>10</sup>, L. Burtscher<sup>9</sup>, W. C. Danchi<sup>14</sup>, C. Dominik<sup>15</sup>, J. Drevon<sup>1</sup>, V. Gámez Rosas<sup>9</sup>, M. R. Hogerheijde<sup>9,15</sup>, J. Hron<sup>16</sup>, L. Klarmann<sup>10</sup>, E. Kokoulina<sup>1</sup>, E. Lagadec<sup>1</sup>, J. Leftley<sup>1</sup>, L. Mosoni<sup>17</sup>, N. Nardetto<sup>1</sup>, C. Paladini<sup>18</sup>, E. Pantin<sup>19</sup>, D. Schertl<sup>8</sup>, P. Stee<sup>1</sup>, L. Szabados<sup>13,20</sup>, R. Waters<sup>21,22</sup>, S. Wolf<sup>23</sup>, and G. Yoffe<sup>9</sup>

*(Affiliations can be found after the references)*

Received 22 October 2021 / Accepted 6 December 2021

## ABSTRACT

**Context.** VX Sgr is a cool, evolved, and luminous red star whose stellar parameters are difficult to determine, which affects its classification.

**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<sup>5</sup>BOLD and a uniform disc model to estimate the apparent diameter and interpret the stellar surface structures. Moreover, we employed the radiative transfer codes OPTIM3D and RADMC3D to compute the spectral energy distribution for the L, M, and N bands, respectively.

**Results.** MATISSE observations unveil, for the first 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 firmly conclude on the interpretation of the current data because of the difficulty in constraining the model parameters using the limited accuracy of our absolute flux 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 first time.

**Key words.** stars: atmospheres – stars: late-type – stars: individual: VX Sgr – techniques: interferometric – infrared: stars – stars: imaging

## 1. Introduction

VX Sagittarii (HD 165674, VX Sgr) is a cool and luminous evolved star. It is a semi-regular variable (732-day period, Kholopov et al. 1987) with unusually large magnitudes spanning from a spectral type M5.5 to M9.8 at the visual maximum and minimum phases, respectively (Lockwood & Wing 1982; Kiss et al. 2006). VX Sgr has been categorised as an oxygen-rich star (Speck et al. 2000) with an effective temperature ranging over three values: 2900 K (García-Hernández et al. 2007), 3150 K (Liu et al. 2017), and ~3700 K (Arroyo-Torres et al. 2013). VX Sgr manifests a strong, asymmetric, and variable mass loss with rates of ~[1–6] × 10<sup>-5</sup> M<sub>☉</sub> yr<sup>-1</sup> (Chapman & Cohen 1986; De Beck et al. 2010; Mauron & Josselin 2011; Liu et al. 2017; Gordon et al. 2018; Gail et al. 2020).

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<sub>2</sub>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 Fe-containing silicates (Gail et al. 2020). Observations in the radio domain, with H<sub>2</sub>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.

<footnote>* Based on the observations made with VLTI-ESO Paranal, Chile under the programme IDs 0103.D-0153(D, E, G). The data are available at oidb.jmmc.fr</footnote>
Using a 22 GHz maser's emission of H₂O, Xu et al. (2018) measured a distance of 1.56<sup>+0.11</sup><sub>-0.10</sub> 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<sup>+0.655</sup><sub>-0.109</sub> 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, 2018). As a consequence of these uncertainties, the radius of VX Sgr ranges between 1000 and 2000 R<sub>☉</sub> (Monnier et al. 2004; Chiavassa et al. 2010; Xu et al. 2018) as well as the luminosity between 1.95 × 10⁵ and 5.5 × 10⁵ L<sub>☉</sub> (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 obser-
vations 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 obser-
vational 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 pho-
tometric 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 prop-
erties of the stellar surface and circumstellar environment in the
wavelength bands: *L*, *M*, and *N*.

2. Interferometric observations

2.1. Short description of the instrument

Multi AperTure Interferometric and SpectroScopic Experi-
ment (MATISSE; Lopez et al. 2022) is the second-generation
four-telescope interferometric beam combiner of the Very
Large Telescope Interferometer (VLTI), which enables mil-
liarcsecond(s) angular resolution imaging in the *L* (3–4 µm),
*M* (4.5–5 µm), and *N* (8–13 µm) mid-infrared bands. It pro-
vides exquisite spatial resolution: *λ*/2*B* = 2.25 mas, and *λ*/*B* =
7.50 mas for a baseline of *B* = 140 m and for wavelengths *λ* =
3 µm and *λ* = 10 µm, respectively. Moreover, it can provide
spectrally dispersed data at low (*R* = *λ*/δ*λ* = 35, the one used
in this work), medium (*R* = 500), high (*R* = 960), and high+
(*R* = 3600) spectral resolutions in *L* and/or *M* bands, and low
(*R* = 35) or high (*R* = 300) in *N* band.

In each of the *L*, *M*, and *N* bands, MATISSE uses a multi-
axial 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 so-
called 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 back-
ground 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 photom-
etry measurements for faint targets in the *L* and *M* bands. The
second modulation is an optical path difference (OPD modula-
tion, 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, detec-
tor integrations, OPD modulation) are synchronised all together
to ensure minimum overheads when observing. Finally, a beam
modulation (or beam commutation) is applied to the observa-
tions 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.

## 2.2. 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 con-
figurations within 15 days in June, while the large configuration
could only be observed late at the end of August 2019. The cal-
ibrators 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 reduc-
tion software (matisse-drs, Millour et al. 2016). It adopts the
classical Fourier transforms scheme (Perrin 2003), with pho-
tometric 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 sum-
ming spectral densities<fnref n="1" />, bispectra<fnref n="2" />, and interspectra<fnref n="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

<footnote n="1">¹ Squared modulus of the Fourier transform of the interferograms.</footnote>

<footnote n="2">² The product of the complex coherent fluxes (i.e. the Fourier trans-
form of the interferograms) over a baseline triplet. The argument of such
quantity is the closure phase.</footnote>

<footnote n="3">³ The product of the complex coherent flux with the complex coherent flux measured at a specific wavelength (or a mean of the complex coherent flux over the spectral bandwidth). The argument of such quantity is the differential phase.</footnote>
<table><caption>**Table 1.** Summary of all the data obtained for VX Sgr and its calibrators.</caption><thead><tr><th>Date</th><th>UT <sup>(a)</sup></th><th>Seeing science <sup>(b)</sup><br>(arcsec)</th><th>tau0 science <sup>(c)</sup><br>(ms)</th><th>Calibrators <sup>(d)</sup></th><th>Seeing calibrator<br>(arcsec)</th><th>tau0 calibrator<br>(ms)</th><th>AT configuration</th></tr></thead><tbody><tr><td>2019-06-09</td><td>04:33:19</td><td>0.6</td><td>2.8</td><td>ν Oph</td><td>0.6</td><td>3.4</td><td>K0G2D0J3/compact</td></tr><tr><td>2019-06-09</td><td>05:41:38</td><td>0.4</td><td>3.7</td><td>β Sgr/ν Oph</td><td>1.2/0.7</td><td>2.1/2.2</td><td>K0G2D0J3/compact</td></tr><tr><td>2019-06-09</td><td>07:21:32</td><td>0.4</td><td>3.7</td><td>ν Oph</td><td>1.2</td><td>2.1</td><td>K0G2D0J3/compact</td></tr><tr><td>2019-06-09</td><td>08:58:50</td><td>1.6</td><td>2.0</td><td>β Sgr</td><td>1.1</td><td>2.8</td><td>K0G2D0J3/compact</td></tr><tr><td>2019-06-25</td><td>01:02:23</td><td>0.6</td><td>2.8</td><td>δ Sgr</td><td>0.9</td><td>0.9</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-25</td><td>02:41:26</td><td>1.1</td><td>0.7</td><td>δ Sgr</td><td>1.1</td><td>1.1</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-25</td><td>03:38:24</td><td>1.0</td><td>0.8</td><td>δ Sgr</td><td>1.0</td><td>0.7</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-25</td><td>04:33:20</td><td>1.0</td><td>0.8</td><td>β Sgr/δ Sgr</td><td>1.1/1.0</td><td>0.8/1.0</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-25</td><td>06:12:47</td><td>1.1</td><td>1.0</td><td>δ Sgr</td><td>1.1</td><td>0.8</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-25</td><td>07:04:58</td><td>1.2</td><td>0.8</td><td>δ Sgr</td><td>1.2</td><td>1.0</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-06-26</td><td>00:42:32</td><td>1.6</td><td>0.8</td><td>δ Sgr</td><td>1.3</td><td>0.8</td><td>A0B2D0C1/medium</td></tr><tr><td>2019-08-29</td><td>00:16:25</td><td>0.5</td><td>3.3</td><td>β Sgr/ν Oph/α Ser</td><td>0.9/0.8/1.4</td><td>1.5/1.7/1.5</td><td>A0G1J2J3/Large</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>UT time at the start of the data acquisition. <sup>(b)</sup>Seeing during the observation for the science target. <sup>(c)</sup>Coherence time in the visible spectral region of the science target. <sup>(d)</sup>Calibrator stars used. Some parameters of the calibrators are reported in Table 2.

<table><caption>**Table 2.** Parameters of the calibrators for VX Sgr extracted from Cruzalèbes et al. (2019) (left) and stellar parameters of the synthetic spectra (PHOENIX grid, Husser et al. 2013) used for the calibrators (see text): effective temperature ($T_{eff}$), surface gravity (log g), and metallicity [$Fe/H$] (right).</caption><thead><tr><th>Name</th><th>Spectral type</th><th>LDD <sup>(a)</sup> diameter<br>(mas)</th><th>T<sub>eff</sub><br>(K)</th><th>log g<br>(cgs)</th><th>[Fe/H]</th></tr></thead><tbody><tr><td>HD 163917 (ν Oph)</td><td>G9IIIa</td><td>2.789 ± 0.005</td><td>5100</td><td>2.5</td><td>+0.5</td></tr><tr><td>HD 188603 (β Sgr)</td><td>K2.5Ib</td><td>2.895 ± 0.013</td><td>4200</td><td>1.0</td><td>+0.5</td></tr><tr><td>HD 168454 (δ Sgr)</td><td>K2.5IIIa</td><td>5.874 ± 0.026</td><td>4600</td><td>1.5</td><td>+0.5</td></tr><tr><td>HD 140573 (α Ser)</td><td>K2IIIb</td><td>4.770 ± 0.013</td><td>4200</td><td>3.0</td><td>0.0</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Limb-Darkened Disc (LDD).

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

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** UV coverage. The spectral coverage contains wavelengths between 2.8 and 13.0 µm. The resolution of the interferometer corresponds to a beam of a 132 × 132 m telescope, 3.1 × 3.1 mas and 7.8 × 7.8 mas at 4.0 and 10.0 µm, respectively.</figcaption></figure>

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.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Representative example of visibility amplitudes (**top panels**) and closure phases (**bottom panels**) for the **L** and **M** bands (**left column**) and the **N** band (**right column**). The different colours indicate different baselines for visibilities and triplets for closures phases. The spectral region between 4.27 and 4.44 µm suffers from noise due to the telluric lines in this range. The night shown in this plot is 25 June 2019, as in Fig. 3.</figcaption></figure>

### 3. 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<fnref n="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_{ij} = \frac{B_{ij}}{\lambda}$), $B_{ij}$ being the projected baseline between telescopes *i* and *j*, and $\lambda$ the wavelength of observation. On the other hand, phases are computed as the FT phase at $f_{ij}$. The closure phases are computed as the sum of the three phases at the three observation spatial frequencies ($f_{ij}$, $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 ($\chi^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

scaled together (thanks to the use of the “hyperparameter” ($\mu^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⁻¹ 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 $\mu = 5 \times 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).

<footnote n="4">⁴ https://software.llnl.gov/yorick-doc/</footnote>

<footnote n="5">⁵ The hyperparameter is used to adjust the weight of the constraints set by the measurements relatively to the ones set by the priors.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Black curve: representative example of calibrated MATISSE flux of VX Sgr (see Sect. 2.2) in *L* and *M* bands (**top**) and *N* band (**low**) overplotted to the telluric spectrum at the observing site (red) computed with the ESO Skycalc web application and sampled to the same resolving power of the MATISSE observations. The night shown in this plot is 25 June 2019.</figcaption></figure>

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 $\lambda$ ~ 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.

# 4. 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.

## 4.1. 3D simulations of AGB and RSG stars

We used the radiation-hydrodynamics (RHD) simulations of stellar surface convection computed with the CO<sup>5</sup>BOLD (Freytag et al. 2012) code. The code solves the coupled non-linear 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_{\text{eff}}$ and $L_{\star}$, 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.

## 4.2. 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 $\mu$ (where $\mu$ = cos($\theta$), with $\theta$ 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Reconstructed images in *L* and *M* bands. They have 256 × 256 pixels with 0.39 mas pixel⁻¹ and are convolved with the interferometric beam (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⁻⁴. 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.</figcaption></figure>

<table><caption>**Table 3.** Parameters of the RHD simulations used in this work.</caption><thead><tr><th>Simulation</th><th>Stellar type</th><th><i>M</i><sub>★</sub> (<i>M</i><sub>☉</sub>)</th><th><i>L</i><sub>★</sub> (<i>L</i><sub>☉</sub>)</th><th><i>R</i><sub>★</sub> (<i>R</i><sub>☉</sub>)</th><th><i>T</i><sub>eff</sub> (K)</th><th>log <i>g</i> (cgs)</th><th><i>t</i><sub>avg</sub> (yr)</th><th>Grid (grid points)</th><th><i>x</i><sub>box</sub> (<i>R</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td>st29gm06n001 <sup>(a)</sup></td><td>AGB</td><td>1</td><td>6956.3 ± 547.4</td><td>350.43 ± 11.3</td><td>2814.7 ± 68.9</td><td>-0.65 ± 0.03</td><td>25.35</td><td>281<sup>3</sup></td><td>1381</td></tr><tr><td>st35gm04n38 <sup>(b)</sup></td><td>RSG</td><td>5</td><td>41517.3 ± 1074.4</td><td>582.03 ± 4.7</td><td>3414.2 ± 16.8</td><td>-0.40 ± 0.01</td><td>11.46</td><td>401<sup>3</sup></td><td>1631</td></tr></tbody></table>

**Notes.** The first column shows the simulation name, and the following five columns show the stellar parameters such as the mass (*M*<sub>★</sub>), the average luminosity (*L*<sub>★</sub>), radius (*R*<sub>★</sub>), effective temperature (*T*<sub>eff</sub>), and surface gravity (log *g*). The different quantities are averaged over spherical shells and epochs (7th column, *t*<sub>avg</sub>). Errors are one standard-deviation fluctuations with respect to the time average (see Chiavassa et al. 2011a, 2009). The solar metallicity is assumed. <sup>(a)</sup>Simulation from the AGB grid of Freytag et al. (2017). <sup>(b)</sup>Simulation of the RSG star presented in detail in Kravchenko et al. (2019).

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 (Θ<sub>AGB</sub> = 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₂O and CO molecules as dominant absorber in the outer layers.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Reconstructed images in the *N* band. They have 512 × 512 pixels with 0.59 mas pixel⁻¹ and are convolved with the interferometric beam (Fig. 1). The peak intensity is normalised between [0, 1] and the corresponding noise maps are shown in Fig. A.2: the highest noise value is ~3 × 10⁻⁵. The contour lines correspond to [0.20, 0.30, 0.40, 0.50, 0.70, 0.95]. East is to the left, north is up. As reported in Sect. 2.2, the data quality of the images at wavelengths between 9.3 and 9.9 µm, is compromised by the telluric spectrum at the observing site.</figcaption></figure>

Figure 8 displays all the measured angular diameters of 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<fnref n="6" /> that, in turn, make the stellar size larger.

<table><caption>**Table 4.** Apparent diameters (in mas) in the *L*, *M*, and *N* bands obtained from a uniform disc ($\Theta_{UD}$) model fitted to the MATISSE squared visibilities and fitting the synthetic RSG ($\Theta_{RSG}$) and AGB ($\Theta_{AGB}$) visibilities (computed from the RHD simulations of Table 3).</caption><thead><tr><th>Wavelength (µm)</th><th>Θ<sub>UD</sub><br>(mas)</th><th>Θ<sub>RSG</sub><br>(mas)</th><th>Θ<sub>AGB</sub><br>(mas)</th></tr></thead><tbody><tr><td colspan="4" style="text-align:center;">L and M bands</td></tr><tr><td>2.92 ↓</td><td>26.16 ± 0.17</td><td>26.83 ± 0.29</td><td>20.29 ± 0.22</td></tr><tr><td>2.99</td><td>26.21 ± 0.14</td><td>26.77 ± 0.25</td><td>20.13 ± 0.18</td></tr><tr><td>3.05</td><td>26.38 ± 0.14</td><td>27.10 ± 0.26</td><td>20.86 ± 0.19</td></tr><tr><td>3.11</td><td>25.44 ± 0.15</td><td>26.30 ± 0.29</td><td>20.07 ± 0.20</td></tr><tr><td>3.17</td><td>24.09 ± 0.16</td><td>24.07 ± 0.24</td><td>18.99 ± 0.15</td></tr><tr><td>3.23</td><td>21.94 ± 0.12</td><td>22.80 ± 0.18</td><td>18.82 ± 0.13</td></tr><tr><td>3.29</td><td>21.36 ± 0.11</td><td>22.35 ± 0.22</td><td>18.27 ± 0.14</td></tr><tr><td>3.35 ↑</td><td>20.55 ± 0.12</td><td>21.22 ± 0.20</td><td>18.31 ± 0.14</td></tr><tr><td>3.41</td><td>20.32 ± 0.13</td><td>21.27 ± 0.21</td><td>18.61 ± 0.16</td></tr><tr><td>3.47</td><td>20.00 ± 0.13</td><td>21.09 ± 0.22</td><td>19.35 ± 0.17</td></tr><tr><td>3.54</td><td>19.71 ± 0.13</td><td>21.05 ± 0.23</td><td>19.84 ± 0.19</td></tr><tr><td>3.60</td><td>19.63 ± 0.14</td><td>20.99 ± 0.24</td><td>20.14 ± 0.21</td></tr><tr><td>3.66</td><td>19.53 ± 0.14</td><td>21.18 ± 0.25</td><td>20.40 ± 0.28</td></tr><tr><td>3.72</td><td>19.95 ± 0.16</td><td>21.32 ± 0.29</td><td>20.25 ± 0.30</td></tr><tr><td>3.78</td><td>20.28 ± 0.18</td><td>21.52 ± 0.33</td><td>21.50 ± 0.35</td></tr><tr><td>3.84</td><td>20.60 ± 0.20</td><td>21.83 ± 0.36</td><td>22.55 ± 0.39</td></tr><tr><td>3.90</td><td>21.01 ± 0.21</td><td>21.91 ± 0.38</td><td>22.48 ± 0.40</td></tr><tr><td>3.96</td><td>21.67 ± 0.22</td><td>22.92 ± 0.39</td><td>23.96 ± 0.42</td></tr><tr><td>4.02 ↓</td><td>22.73 ± 0.20</td><td>23.77 ± 0.35</td><td>23.72 ± 0.39</td></tr><tr><td>4.08</td><td>23.87 ± 0.21</td><td>24.25 ± 0.32</td><td>22.67 ± 0.32</td></tr><tr><td>4.14</td><td>25.29 ± 0.22</td><td>25.26 ± 0.35</td><td>23.10 ± 0.31</td></tr><tr><td>4.20</td><td>27.02 ± 0.25</td><td>26.75 ± 0.40</td><td>24.33 ± 0.34</td></tr><tr><td>4.27</td><td>28.80 ± 0.29</td><td>28.36 ± 0.48</td><td>25.14 ± 0.41</td></tr><tr><td>4.33</td><td>31.84 ± 0.35</td><td>31.11 ± 0.58</td><td>28.74 ± 0.54</td></tr><tr><td>4.39</td><td>34.10 ± 0.37</td><td>31.64 ± 0.59</td><td>28.54 ± 0.54</td></tr><tr><td>4.45</td><td>38.55 ± 0.48</td><td>33.50 ± 0.66</td><td>28.79 ± 0.51</td></tr><tr><td>4.51</td><td>38.06 ± 0.50</td><td>32.36 ± 0.64</td><td>28.43 ± 0.52</td></tr><tr><td>4.57</td><td>38.65 ± 0.54</td><td>32.21 ± 0.67</td><td>27.71 ± 0.54</td></tr><tr><td>4.63</td><td>41.38 ± 0.68</td><td>33.85 ± 0.77</td><td>28.14 ± 0.57</td></tr><tr><td>4.69 ↑</td><td>44.15 ± 0.98</td><td>28.88 ± 0.56</td><td>25.73 ± 0.55</td></tr><tr><td colspan="4" style="text-align:center;">N band</td></tr><tr><td>8.10</td><td>133.66 ± 5.00</td><td>—</td><td>91.87 ± 6.55</td></tr><tr><td>8.40</td><td>149.59 ± 3.84</td><td>—</td><td>103.64 ± 5.19</td></tr><tr><td>8.70</td><td>164.17 ± 3.89</td><td>—</td><td>116.32 ± 5.63</td></tr><tr><td>9.00</td><td>179.36 ± 3.26</td><td>—</td><td>128.73 ± 4.78</td></tr><tr><td>9.30 ↓</td><td>191.71 ± 2.83</td><td>—</td><td>143.71 ± 4.66</td></tr><tr><td>9.59</td><td>202.89 ± 1.73</td><td>—</td><td>153.91 ± 3.54</td></tr><tr><td>9.89 ↑</td><td>210.49 ± 1.91</td><td>—</td><td>157.80 ± 3.73</td></tr><tr><td>10.19</td><td>215.67 ± 2.36</td><td>—</td><td>171.12 ± 4.89</td></tr><tr><td>10.49</td><td>221.74 ± 2.85</td><td>—</td><td>174.10 ± 5.62</td></tr><tr><td>10.79</td><td>226.60 ± 4.02</td><td>—</td><td>173.66 ± 6.84</td></tr><tr><td>11.09</td><td>222.91 ± 5.00</td><td>—</td><td>173.46 ± 11.07</td></tr><tr><td>11.39</td><td>230.87 ± 5.03</td><td>—</td><td>178.95 ± 9.12</td></tr><tr><td>11.68</td><td>229.73 ± 4.68</td><td>—</td><td>176.73 ± 8.59</td></tr><tr><td>11.98</td><td>221.32 ± 5.28</td><td>—</td><td>173.06 ± 10.96</td></tr><tr><td>12.28</td><td>216.82 ± 5.53</td><td>—</td><td>161.03 ± 10.03</td></tr><tr><td>12.58</td><td>221.34 ± 6.35</td><td>—</td><td>161.99 ± 10.10</td></tr><tr><td>12.88</td><td>223.97 ± 7.87</td><td>—</td><td>169.31 ± 12.88</td></tr></tbody></table>

**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).

<footnote n="6">⁶ As formulated for the first time by Tsuji (2000), it is a spherical region that is optically thin for particular molecular lines.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Top row: colour composite images in the *H* and *K* bands (observation taken in 2008 with the AMBER instrument, with 21 spectral channels combined, Chiavassa et al. 2010), *L* and *M* (31 spectral channels from Fig. 4 combined), and *N* band (17 spectral channels from Fig. 5 combined). Bottom row, left panel: colour composite images for the *H* − *K* (pink scale) and *L* − *M* (blue scale) bands. Bottom row, right panel: all the bands, including *N* (green scale), are combined. The photospheric structures (*H* to *M* bands) are located in the centre of the *N* band maps.</figcaption></figure>

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₂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.

4.3. 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⁻¹ 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<sub>rms</sub>/⟨I⟩ for each snap-
shot, where *I* is the intensity. This is performed only for the
pixels within the stellar surface with a cut at 70% of the max-
imum of the ring-averaged intensity profile<sup>7</sup>. 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 χ<sup>2</sup><sub>snapshot</sub> =

|C<sub>observations</sub> - C<sub>3D</sub>|² / (σ²<sub>C<sub>observations</sub></sub> + σ²<sub>C<sub>3D</sub></sub>) , where C<sub>observations</sub> and C<sub>3D</sub> are the intensity con-
trasts for the reconstructed image (observations) and the syn-
thetic one (3D), while σ<sub>C<sub>observations</sub></sub> is the standard deviation of
the observed intensity contrast. Afterwards, we computed the
reduced χ² as χ² = (1/N - 1) Σ χ²<sub>snapshot</sub> where N is the total number

<footnote n="7">⁷ A similar method based on the stellar radial cut is used in Chiavassa et al. (2020) and Montargès et al. (2018).</footnote>
<figure><img src="image_8.png" /></figure>

**Fig. 7.** Intensity profiles (normalised to their maximum) as a function of the wavelength in *L* and *M* bands for a particular snapshot of the RSG (*top panel*) and the AGB simulations (*bottom panel*) of 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 $\bar{\chi}^2$.

In the *L* and *M* bands, the best matching RSG snapshot has $\bar{\chi}^2=1.12$, and the one with the worst 29.80. For the AGB's best snapshot, $\bar{\chi}^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 $\bar{\chi}^2$ values are very close. In the end, our 3D RHD simulations return a good agreement with the observations in terms

<figure><img src="image_9.png" /></figure>

**Fig. 8.** Angular diameters measured (Table 4) with the RHD simulations of Table 3 and a UD model in the *L* and *M* bands (*top panel*) and the *N* band (*bottom panel*). The green horizontal lines highlight the principal chemical species.

of contrast of surface structures, meaning that they are adequate for interpreting the observed photospheric inhomogeneities as convection-related surface structures.

## 5. 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.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Intensity contrast in the case of the AGB (top) and the RSG (bottom) simulations as a function of the wavelength in the L and M bands. Each coloured point indicates a different snapshot and the vertical bar on the right gives its number in the simulated temporal coverage. The black line displays the contrast of reconstructed images of Fig. 4 and the red line the best matching snapshot across the wavelength bins.</figcaption></figure>

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 ($\rho$) gradient ($\rho_{max} - \rho_{min}$)/$\rho_{max}$ > 0.05, $\rho_{min}$ and $\rho_{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 $\alpha$ 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₂Al₂SiO₇, Mutschke et al. 1998), corundum (Al₂O₃ Zeidler et al. 2013), and olivine (MgFeSiO₄, 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<fnref n="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.5–5] µm, with a size distribution of $n(a) \propto 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_\star$ to match the reconstructed images in Fig. 5. The dust mass-loss rate is set at $2 \times 10^{-7} M_\odot \text{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.

## 6. 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⁵BOLD, to calculate synthetic intensity maps in the same observed spectral channels of the MATISSE observations. As a first step, we measured VX Sgr monochromatic 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₂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). As a consequence, it makes the levitation of material possible that

<footnote n="8">⁸ https://www.astro.uni-jena.de/Laboratory/OCDB/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** *Top row:* reconstructed image (left) and synthetic intensity maps computed for the best matching RHD simulations of Table 3 at 4.57 µm. The intensity is normalised between [0, 1]. The images, 256 × 256 pixels with 0.39 mas pixel⁻¹, 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 $\bar{\chi}^2 = 1.12$, while the AGB one is $\bar{\chi}^2 = 1.98$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** MATISSE calibrated fluxes for the three observing nights (Table 1; see Sect. 2.2 for details for the details about calibration).</figcaption></figure>

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₂Al₂SiO₇), corundum (Al₂O₃), and olivine (MgFeSiO₄). 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
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Top panel: SEDs computed with RADMC3D compared to MATISSE calibrated flux from Fig. 11. The black curve shows the model with all three dust species included (melilite, corundum, and olivine), while the different colours (red, orange, and light blue) show one species at a time. Bottom panel: SEDs computed with RADMC3D for olivine only and different grain size ranges.</figcaption></figure>

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 circum-
stellar 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, multi-
wavelength simultaneous observations will be key to probing the
stellar parameters of VX Sgr and characterising its photometric
variability.

Acknowledgements. This work was granted access to the HPC resources of Observatoire de la Côte d'Azur - Mésocentre SIGAMM. This research has also benefited from the help of SUV (http://www.jmmc.fr/suv.htm), the VLTI user support service of the Jean-Marie Mariotti Center (http://www.jmmc.fr). Moreover, this work was supported by the "Programme National de Physique

Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. BF acknowl-
edges funding from the European Research Council (ERC) under the European
Union's Horizon 2020 research and innovation programme Grant agreement no.
883867, project EXWINGS) and the Swedish Research Council (Vetenskap-
srådet, grant number 2019-04059). The computations of 3D models were enabled
by resources provided by the Swedish National Infrastructure for Computing
(SNIC) at UPPMAX. V.H. is supported by the National Science Center, Poland,
Sonata BIS project 2018/30/E/ST9/00598. This research made use of IPython,
Numpy, Matplotlib, SciPy, and Astropy (available at http://www.astropy.
org/), a community-developed core Python package for Astronomy (Astropy
Collaboration 2013).

References

Arroyo-Torres, B., Wittkowski, M., Marcaide, J. M., & Hauschildt, P. H. 2013, A&A, 554, A76

Arroyo-Torres, B., Wittkowski, M., Chiavassa, A., et al. 2015, A&A, 575, A50

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615

Chapman, J. M., & Cohen, R. J. 1986, MNRAS, 220, 513

Chen, X., Shen, Z.-Q., & Xu, Y. 2007, Chinese J. Astron. Astrophys., 7, 531

Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351

Chiavassa, A., Lacour, S., Millour, F., et al. 2010, A&A, 511, A51

Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011a, A&A, 535, A22

Chiavassa, A., Pasquato, E., Jorissen, A., et al. 2011b, A&A, 528, A120

Chiavassa, A., Freytag, B., & Schultheis, M. 2018, A&A, 617, L1

Chiavassa, A., Kravchenko, K., Millour, F., et al. 2020, A&A, 640, A23

Coudé du Foresto, V., Ridgway, S., & Mariotti, J. M. 1997, A&AS, 121, 379

Cruzalèbes, P., Petrov, R. G., Robbe-Dubois, S., et al. 2019, MNRAS, 490, 3158

De Beck, E., Decin, L., de Koter, A., et al. 2010, A&A, 523, A18

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Dullemond, C. P. 2012, RADMC-3D: A multi-purpose radiative transfer tool, astrophysics Source Code Library

Duvert, G., Young, J., & Hummel, C. A. 2017, A&A, 597, A8

Freytag, B., Steffen, M., Ludwig, H.-G., et al. 2012, J. Comput. Phys., 231, 919

Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137

Gail, H.-P., Tamanai, A., Pucci, A., & Dohmen, R. 2020, A&A, 644, A139

García-Hernández, D. A., García-Lario, P., Plez, B., et al. 2007, A&A, 462, 711

Gordon, M. S., Humphreys, R. M., Jones, T. J., et al. 2018, AJ, 155, 212

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Höfner, S., & Olofsson, H. 2018, A&ARv, 26, 1

Humphreys, R. M. 1974, ApJ, 188, 75

Humphreys, R. M., Strecker, D. W., & Ney, E. P. 1972, ApJ, 172, 75

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Jaeger, C., Mutschke, H., Begemann, B., Dorschner, J., & Henning, T. 1994, A&A, 292, 641

Kholopov, P. N., Samus', N. N., Kazarovets, E. V., & Kireeva, N. N. 1987, Inform. Bull. Variable Stars, 3058, 1

Kiss, L. L., Szabó, G. M., & Bedding, T. R. 2006, MNRAS, 372, 1721

Kravchenko, K., Chiavassa, A., Van Eck, S., et al. 2019, A&A, 632, A28

Lançon, A., Hauschildt, P. H., Ladjal, D., & Mouhcine, M. 2007, A&A, 468, 205

Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973

Liu, J., Jiang, B. W., Li, A., & Gao, J. 2017, MNRAS, 466, 1963

Lockwood, G. W., & Wing, R. F. 1982, MNRAS, 198, 385

Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202141785

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Mauron, N., & Josselin, E. 2011, A&A, 526, A156

Millour, F., Berio, P., Heininger, M., et al. 2016, SPIE Conf. Ser., 9907, 990723

Monnier, J. D., Millan-Gabet, R., Tuthill, P. G., et al. 2004, ApJ, 605, 436

Montargès, M., Norris, R., Chiavassa, A., et al. 2018, A&A, 614, A12

Montargès, M., Cannon, E., Lagadec, E., et al. 2021, Nature, 594, 365

Mutschke, H., Begemann, B., Dorschner, J., et al. 1998, A&A, 333, 188

Paladini, C., Klotz, D., Sacuto, S., et al. 2017, A&A, 600, A136

Pauls, T. A., Young, J. S., Cotton, W. D., & Monnier, J. D. 2005, PASP, 117, 1255
Perrin, G. 2003, A&A, 400, 1173

Petrov, R. G., Malbet, F., Weigelt, G., et al. 2007, A&A, 464, 1

Reid, M. J., Muhleman, D. O., Moran, J. M., Johnston, K. J., & Schwartz, P. R. 1977, ApJ, 214, 60

Scicluna, P., Kemper, F., Siebenmorgen, R., et al. 2020, MNRAS, 494, 3200

Speck, A. K., Barlow, M. J., Sylvester, R. J., & Hofmeister, A. M. 2000, A&AS, 146, 437

Su, J. B., Shen, Z. Q., Chen, X., & Jiang, D. R. 2018, ApJ, 853, 42

Tabernero, H. M., Dorda, R., Negueruela, I., & Marfil, E. 2021, A&A, 646, A98

Tatulli, E., Millour, F., Chelli, A., et al. 2007, A&A, 464, 29
Thiébaut, E. 2008, SPIE Conf. Ser., 7013, 701311

Tremblay, P. E., Ludwig, H. G., Freytag, B., Steffen, M., & Caffau, E. 2013, A&A, 557, A7

Tsuji, T. 2000, ApJ, 538, 801

van Leeuwen, F. 2007, A&A, 474, 653

Verhoelst, T., van der Zypen, N., Hony, S., et al. 2009, A&A, 498, 127

Wittkowski, M., Chiavassa, A., Freytag, B., et al. 2016, A&A, 587, A12

Wittkowski, M., Bladh, S., Chiavassa, A., et al. 2019, The Messenger, 178, 34

Xu, S., Zhang, B., Reid, M. J., et al. 2018, ApJ, 859, 14

Yoon, D.-H., Cho, S.-H., Yun, Y., et al. 2018, Nat. Commun., 9, 2534

Zeidler, S., Posch, T., & Mutschke, H. 2013, A&A, 553, A81

¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lagrange, CS 34229, Nice, France
e-mail: andrea.chiavassa@oca.eu

² European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

³ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany

⁴ Max Planck Institute for extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany

5 LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France

⁶ Theoretical Astrophysics, Department of Physics and Astronomy at Uppsala University, Regementsvägen 1, Box 516, 75120 Uppsala, Sweden

⁷ Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, Bartycka 18, 00-716 Warszawa, Poland

⁸ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

⁹ Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

¹⁰ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

¹¹ I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937, Köln, Germany

¹² Univ. Grenoble Alpes, CNRS, IPAG, 38000, Grenoble, France

¹³ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary

¹⁴ NASA Goddard Space Flight Center, Astrophysics Division, Greenbelt, MD 20771, USA

¹⁵ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090 GE Amsterdam, The Netherlands

¹⁶ Department of Astrophysics, University of Vienna, 1180 Vienna, Türkenschanzstrasse 17, Austria

¹⁷ Zselic Park of Stars, 064/2 hrsz., 7477 Zselickisfalud, Hungary

¹⁸ European Southern Observatory, Alonso de Cordova 3107, Casilla 19, Santiago 19001, Chile

¹⁹ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France

²⁰ CSFK Lendület Near-Field Cosmology Research Group, Budapest, Hungary

²¹ Institute for Mathematics, Astrophysics and Particle Physics, Radboud University, PO Box 9010, MC 62, 6500 GL Nijmegen, the Netherlands

²² SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands

²³ Institute of Theoretical Physics and Astrophysics, Kiel University, Leibnizstr. 15, 24118 Kiel, Germany
Appendix A: Intensity RMS maps of reconstructed images

<figure><img src="image_14.png" /><figcaption>Fig. A.1. Intensity RMS of the reconstructed images in the *L* and *M* bands from Fig. 4. 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. The z-axis scale is the same for all the panels.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. A.2.** Same as Fig. A.1, but for the N band. The telluric spectrum dominates between 9.3 and 9.9 µm.</figcaption></figure>