11institutetext: Nicolaus Copernicus Astronomical Center, Polish Academy of
Sciences, Rabiańska 8, 87-100 Toruń, Poland 22institutetext: Université Côte
d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, France
33institutetext: European Southern Observatory, Karl-Schwarzschild-Str. 2,
85748 Garching bei Munchen, Germany 44institutetext: European Southern
Observatory, Alonso de Cordoba 3107, Vitacura, Santiago, Chile
55institutetext: Astronomy Department, University of Michigan, Ann Arbor, MI
48109, USA 66institutetext: Astrophysics Group, Department of Physics &
Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK
77institutetext: Institut de Planetologie et d’Astrophysique de Grenoble,
Grenoble 38058, France 88institutetext: The CHARA Array of Georgia State
University, Mount Wilson Observatory, Mount Wilson, CA 91203, USA
99institutetext: NASA Ames Research Center, Moffett Field, CA 94035, USA

# Reconstructing the mid-infrared environment in the stellar merger remnant
V838 Monocerotis

Muhammad Zain Mobeen 11 Tomasz Kamiński 11 Alexis Matter 22 Markus Wittkowski
33 John D. Monnier 55 Stefan Kraus 66 Jean-Baptiste Le Bouquin 77 Narsireddy
Anugu 88 Theo Ten Brummelaar 88 Claire L. Davies 66 Jacob Ennis 55 Tyler
Gardner 66 Aaron Labdon 44 Cyprien Lanthermann 88 Gail H. Schaefer 88 Benjamin
R. Setterholm 55 Nour Ibrahim 55 Steve B. Howell 99

###### Abstract

Context. V838 Mon is a stellar merger remnant that erupted in 2002 in a
luminous red novae event. Although it is well studied in the optical, near
infrared and submillimeter regimes, its structure in the mid infrared
wavelengths remains elusive. Over the past two decades only a handful of mid-
infrared interferometric studies have been carried out and they suggest the
presence of an elongated structure at multiple wavelengths. However, given the
limited nature of these observations, the true morphology of the source could
not be determined.

Aims. By performing image reconstruction using observations taken at the VLTI
and CHARA, we aim to map out the circumstellar environment in V838 Mon.

Methods. We observed V838 Mon with the MATISSE ($LMN$ bands) and GRAVITY ($K$
band) instruments at the VLTI and also the MIRCX/MYSTIC ($HK$ bands)
instruments at the CHARA array. We geometrically modelled the squared
visibilities and the closure phases in each of the bands to obtain constraints
on physical parameters. Furthermore, we constructed high resolution images of
V838 Mon in the $HK$ bands using the MIRA and SQUEEZE algorithms to study the
immediate surroundings of the star. Lastly we also modelled the spectral
features seen in the $K$ and $M$ bands at various temperatures.

Results. The image reconstructions show a bipolar structure that surrounds the
central star in the post merger remnant. In the $K$ band the super resolved
images show an extended structure (uniform disk diameter $\sim 1.94$ mas) with
a clumpy morphology that is aligned along a north-west position angle (PA) of
–$40\degr$. Whereas in the $H$ band the extended structure (uniform disk
diameter $\sim 1.18$ mas) lies roughly along the same PA. However the northern
lobe is slightly misaligned with respect to the southern lobe which results in
the closure phase deviations.

Conclusions. The VLTI and CHARA imaging results show that V838 Mon is
surrounded by features that resemble jets that are intrinsically asymmetric.
This is further confirmed by the closure phase modelling. Further observations
with VLTI can help to determine whether this structure shows any variation
over time and also if such bi-polar structures are commonly formed in other
stellar merger remnants.

###### Key Words.:

instrumentation: interferometers – techniques: interferometric – stars:
individual: V838 Monocerotis – circumstellar matter

##  1 Introduction

At the start of 2002 V838 Monocerotis erupted in a luminous red nova event
(Munari et al. 2002b; Tylenda 2005) and in a few weeks brightened by almost
two orders of magnitude, finally reaching a peak luminosity of
$10^{6}L_{\sun}$ (Tylenda 2005; Sparks et al. 2008; Bond et al. 2003). The
event is thought to have been the result of a stellar merger. According to the
scenario proposed in Tylenda & Soker (2006), an 8 $M_{\sun}$ B-type main
sequence star coalesced with a 0.4 $M_{\sun}$ young stellar object. The
outburst was soon followed by a gradual decrease in temperature, and its
spectra soon evolved to that of a late M-type supergiant (Evans et al. 2003;
Loebman et al. 2015). Spectra taken in the 2000s revealed the presence of
various molecules in V838 Mon, including water and transition-metal oxides
(Banerjee & Ashok 2002; Kamiński et al. 2009). Dust was also observed to be
produced in the post merger environment (Wisniewski et al. 2008; Kamiński et
al. 2021). Additionally, a B-type companion was observed in the vicinity of
the central merger remnant, which thus suggests that the merger had taken
place in a hierarchical triple system (Munari et al. 2002a; Kamiński et al.
2021). The companion was obscured by dust formed in the aftermath of the 2002
eruption (Tylenda et al. 2009). V838 Mon is the best studied luminous red nova
in the Milky Way, although, many others have also been found within the Galaxy
as well as elsewhere in the Local Group (Pastorello et al. 2019).

As the merger remnant in V838 Mon is enshrouded by dust, it is therefore an
ideal target for mid-infrared interferometric studies. The first of these
studies was conducted by Lane et al. (2005) in which they observed V838 Mon
using the Palomar Testbed Interferometer (PTI). By modelling the squared
visibilities in the $K$-band at 2.2 $\mu$m they were able to measure the size
of the merger remnant of $1.83\pm 0.06$ mas. There were also hints of
asymmetries in the object, but due to scarce measurements these could not be
confirmed. Chesneau et al. (2014) followed up these measurements between 2011
and 2014, using the Very Large Telescope Interferometer (VLTI) instruments:
Astronomical Multi-BEam combineR (AMBER; Petrov et al. 2007) in $H$ and $K$
bands, and the MID-infrared Interferometric instrument (MIDI; Leinert et al.
2003) in $N$ band. Fitting uniform disk models to the AMBER measurements gave
an angular diameter of $1.15\pm 0.2$ mas, which –according to the authors –
indicates that the photosphere in V838 Mon had contracted by about 40% over
the course of a decade. Also, their modeling of the AMBER data suggests that
an extended component was present in the system, with a lower limit on the
full width at half-maximum (FWHM) of $\sim$20 mas. Modeling of the MIDI
measurements seems to point towards the presence of a dusty elongated
structure whose major axis varies as a function of wavelength between 25 and
70 mas in $N$ band. Submillimeter observations obtained with the Atacama Large
Millimeter/sub millimeter Array (ALMA) in continuum revealed the presence of a
flattened structure with a FWHM of 17.6$\times$17.6 mas surrounding V838 Mon
(Kamiński et al. 2021). Recent $L$ band measurements by Mobeen et al. (2021)
also seem to paint a similar picture. Mobeen et al. (2021) geometrically
modelled the squared visibilities and closure phases in the $L$-band, obtained
using the Multi AperTure mid-Infrared SpectroScopic Experiment instrument
(MATISSE) at the VLTI in 2020. They found that the structure in the $L$-band
is well represented by an elliptical disk tilted at an angle of –40$\degr$.
Furthermore, the closure phases showed small but non-zero deviations, which
suggest the presence of asymmetries in the system. The interferometric
measurements span across the wavebands (from 2.2 $\mu$m to 1.3 mm) and trace a
dusty structure oriented roughly along the same direction, with PA in the
range –10° (MIDI) to –50° (ALMA). This might indicate either a single
overarching structure in the post-merger remnant, or multiple similarly
aligned structures. Simulations of stellar merger events also suggest the
presence of a disk like structure in the post-merger remnant, which is thought
to be a reservoir for the pre-merger binary angular momentum (e.g. Webbink
1976; Pejcha et al. 2017).

V838 Mon serves as an excellent source to advance our understanding of the
post-merger environment decades after the luminous red nova event. Thus, it
provides us with crucial insights into the physical processes at play in these
merger events and their final products in the long term. In this paper, we
analyze and interpret recent interferometric observations obtained with a
variety of instruments that span many near to mid infrared wavelengths.

The format of the paper is as follows. In Sect. 2 we present all of our VLTI
and CHARA observations and outline the main steps of the data reduction. We
also analyze and interpret recent optical speckle interferometric observations
obtained at 562 nm and 832 nm. In Sect. 3 we mainly present the results of
geometrically modelling the interferometric observables (squared visibilities
and closure phases) observed with the MATISSE and GRAVITY instruments at VLTI
and with MIRCX/MYSTIC at CHARA. Sect 4 centers around our image reconstruction
attempts for the VLTI and CHARA datasets using two distinct image
reconstruction algorithms. The modelling and imaging results are discussed in
depth in Sect. 5 finally followed by Sect. 6 in which we present the main
conclusions of this study.

##  2 Observations and data reduction

V838 Mon was observed with the VLTI located at Paranal observatory in Chile,
in 2021 and 2022. Observations were carried out using the 1.8 m Auxiliary
Telescopes (ATs) and two instruments, MATISSE (Lopez et al. 2022) and GRAVITY
(GRAVITY Collaboration et al. 2017). MATISSE is a four telescope beam combiner
which covers the $L$ (2.8–4.2 $\mu$m), $M$ (4.5–5 $\mu$m) and $N$ (8–13
$\mu$m) bands, while GRAVITY combines light in the $K$ (1.9–2.4 $\mu$m) band.
For both instruments, we intended to get 18 observing blocks (OBs) to perform
image reconstruction. However, we were only able to obtain three OBs for
MATISSE and fifteen OBs for GRAVITY. The technical details of our VLTI
observations are presented in Tables 1–2. In the case of MATISSE, only the
large and small configurations were employed, while for GRAVITY, observations
using all three configurations were carried out. In particular, intermediate
configurations were also used to sample better the UV plane. The UV coverages
obtained for the MATISSE and GRAVITY observations are shown in Figs. 1 & 2,
respectively.

Figure 1: UV tracks for MATISSE observing runs in 2021–2022.

MATISSE observations were carried out in low spectral resolution (R$\sim 30$)
using the GRA4MAT mode, in which the GRAVITY fringe tracker is used to
stabilize fringes for MATISSE (Lopez et al. 2022). Each MATISSE observation
for V838 Mon consisted of a CAL-SCI-CAL observing sequence in which two
calibrator stars, 20 Mon (spectral type K0 III) and HD 52666 (spectral type K5
III) were observed to calibrate the $LMN$ bands. The source and calibrator
fluxes in the $LMN$ bands are given in Table 4.

The GRAVITY observations were carried out in medium spectral resolution
(R$\sim$500) using its single-field mode in which the light is split equally
between the fringe tracker channel and the science channel. We adopted a CAL-
SCI-CAL sequence for each GRAVITY run. However, one calibrator was observed
here, HD 54990 (spectral type K0/1 III).

Figure 2: UV coverage for GRAVITY observing runs.

Interferometric observations of V838 Mon were also obtained at the Center for
High Angular Resolution Astronomy (CHARA) Array on UT 2022 Mar 2, 3, and 9.
The CHARA Array is operated by Georgia State University and is located at
Mount Wilson Observatory in southern California. The array consists of six
1-meter telescopes arranged in a Y-configuration with baselines ranging in
length from 34 m to 331 m (ten Brummelaar et al. 2005). We combined the light
from 4 to 5 of the telescopes using the Michigan InfraRed Combiner-eXeter
(MIRC-X) beam combiner in the $H$-band using the low spectral resolution
(R=50) prism that disperses light over 8 spectral channels (Anugu et al.
2020). On the last night, additional simultaneous observations were obtained
using the Michigan Young STar Imager (MYSTIC) combiner (Setterholm et al.
2022) in the $K$-band in low spectral resolution (R=49) mode (Monnier et al.
2018). The CHARA MIRC-X and MYSTIC observation log is given in Table 3. On
each night, we alternated between observations of V838 Mon and calibrator
stars to calibrate the interferometric transfer function. The $H$ band fluxes
are presented in Table 6, while the uniform disk diameters of the calibrators
in the $HK$ bands were adopted from Bourgés et al. (2014) and Bourgés et al.
(2017) and are listed in Table 7.

Figure 3: UV coverage for CHARA observing runs in the $HK$ bands.

###  2.1 MATISSE data reduction

The $LMN$-band MATISSE data were reduced using version 1.7.5 of the ESO data
reduction pipeline. The final products were OIFITS files (version 2)
comprising uncalibrated interferometric observables, including six dispersed
squared visibilities and three independent closure phase measurements per
exposure. The observations of V838 Mon were then calibrated using our
measurements of the calibrator stars listed in Table 4, whose diameters and
fluxes were taken from the Mid-infrared stellar Diameters and Fluxes
compilation Catalog111MDFC is available through VizieR service at
https://vizier.u-strasbg.fr/viz-bin/VizieR (MDFC; Cruzalebes et al. 2019, see
Table 2). The calibration of squared visibilities was performed by dividing
the raw squared visibilities of V838 Mon by those of the calibrators,
corrected for their diameter (which is called the interferometric transfer
function). The calibrated $LMN$ spectra of V838 Mon were obtained by
multiplying the ratio between the target and calibrator raw fluxes measured by
MATISSE at each wavelength by a model of the absolute flux of the calibrator.
This model was taken from the PHOENIX stellar spectra grid (Husser et al.
2013).

Upon completing the data reduction and inspecting the resultant products, we
found that the $N$ band squared visibilities were of poor quality. This was
due to the fact that the average $N$-band total flux of V838 Mon ($\sim$15 Jy)
lies very close to the photometric sensitivity limit of MATISSE in the $N$
band. Furthermore, due to the increased thermal background effects beyond 11
$\mu$m, the quality of the photometric data became even worse. This rendered
the $N$-band photometry, and thus the absolute visibilities, unusable, and
they were discarded. However, we were able to make use of the correlated
fluxes ($\sim$5 Jy), which turn out to be slightly greater than the
sensitivity limit of the GRA4MAT mode in $N$-band. Those observables allowed
us to model geometrically the source in the $N$ band and put constraints on
the extent of the structure in this band.

###  2.2 GRAVITY data reduction

GRAVITY data reduction was performed with the ESO GRAVITY data reduction
pipeline version 1.4.1 This was carried out in the ESO Reflex workflow
environment. For the calibration of squared visibilities we used the
visibility calibration workflow in Reflex. Furthermore, the reduced $K$ band
spectra were flux calibrated similarly to the $L$ band spectra obtained from
MATISSE. The flux calibration routine from the MATISSE consortium was used for
this purpose and, just like in the case of the $M$ band, the $K$ band spectra
were also similarly normalized and later compared to synthetic spectra.

###  2.3 CHARA data reduction

CHARA data was reduced and calibrated by the support astronomers at the CHARA
array. The data were reduced using the standard MIRC-X and MYSTIC pipeline
(version 1.3.5) written in
python222https://gitlab.chara.gsu.edu/lebouquj/mircx_pipeline.git and
described by Anugu et al. (2020). The pipeline produces calibrated visibility
amplitudes for each pair of telescopes and closure phases for each combination
of three telescopes. To assess the data quality, the calibrators were checked
against each other on each night. The calibrators showed no evidence for
binarity based on a visual inspection of the data, and the diameters derived
from the measured visibilities were consistent with the expected values within
uncertainties. The calibrated OIFITS files will be available in the Optical
Interferometry Database333https://oidb.jmmc.fr/index.html and the CHARA Data
Archive444https://www.chara.gsu.edu/observers/database.

###  2.4 Gemini South observation

V838 Mon was also observed twice with the Gemini South 8-m telescope using the
high-resolution Zorro instrument (Scott et al. 2021; Howell & Furlan 2022).
Zorro provides simultaneous speckle imaging in two optical bands, here
centered at 562 nm and 832 nm. V838 Mon was observed on 25 Feb 2021 UT and,
about one year later, on 20 March 2022 UT. The speckle imaging on each night
consisted of a number of 60 ms frames taken in a row, the February 2021
observation taking 12,000 frames and the March 2022 observation taking 16,000
frames. The February observations occurred during a night of average seeing
(0$\aas@@fstack{\prime\prime}$6) while the longer March 2022 observations
occurred during good seeing (0$\aas@@fstack{\prime\prime}$45). The data were
reduced using a standard speckle imaging pipeline with reduced output data
products including reconstructed images and 5$\sigma$ contrast limits (Howell
et al. 2011). The two sets of observations agreed well with the March 2022
data, providing a higher S/N result.

Fig. 4 shows the resultant contrast curves for both filters and our 832 nm
reconstructed image from March 2022. Fig. 5 shows a similar plot but for
February 2021 along with image reconstructions at 832 nm and at 562 nm. There
are no close (¡1$\aas@@fstack{\prime\prime}$2) stellar companions detected
within the angular and contrast limits achieved. However, as seen in detail in
Fig. 6, the image at 832 nm is extended beyond just a point source with a
slight elongation in the north-south direction and to the east by about
0$\aas@@fstack{\prime\prime}$03\. Similarly, at 562 nm, we see also see a
noticeable elongation in the north-western direction, as shown in Fig. 7. The
562 nm elongation is in very good agreement with what we observe in the mid-
infrared bands (see Fig. 30). While the elongation varies significantly at 832
nm, we note that the orientation (north-south) of this particular structure is
remarkably similar to the PA ($\sim$3°) we obtain by fitting the $N$ band
visibility amplitudes with an elliptical disk model. Only with future
extensive observations at Gemini South will it be possible to reliably
constrain the orientation of the elongation of V838 Mon in the $V$ and $I$
bands.

Figure 4: Contrast (in magnitudes) curves for both filters from March 2022. The speckle image reconstruction at 832 nm is shown in the top right corner. Figure 5: Same as 5 but for February 2021. The speckle image reconstruction at 562 nm is also shown. Figure 6: Speckle image reconstruction from March 2022 at 832 nm. A noticeable north-south elongation is visible. The image field of view (FoV) is 0$\aas@@fstack{\prime\prime}$85 by 0$\aas@@fstack{\prime\prime}$55\. North is up, east is left. Figure 7: Speckle image reconstruction from February 2021 at 562 nm. A prominent north-west elongation is visible. The image FoV is 0$\aas@@fstack{\prime\prime}$21 by 0$\aas@@fstack{\prime\prime}$14\. North is up, east is left. Table 1: VLTI/GRAVITY observation log. Date | Configuration | Seeing
---|---|---
01/04/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$42
02/03/2022 | A0-B2-C1-D0 | 0$\aas@@fstack{\prime\prime}$86
02/03/2022 | A0-B2-C1-D0 | 0$\aas@@fstack{\prime\prime}$82
02/04/2022 | A0-B2-C1-D0 | 0$\aas@@fstack{\prime\prime}$90
02/04/2022 | A0-B2-C1-D0 | 1$\aas@@fstack{\prime\prime}$00
02/04/2022 | A0-B2-C1-D0 | 0$\aas@@fstack{\prime\prime}$62
02/05/2022 | A0-D0-G2-J3 | 0$\aas@@fstack{\prime\prime}$57
02/06/2022 | D0-G2-J3-K0 | 0$\aas@@fstack{\prime\prime}$50
02/06/2022 | D0-G2-J3-K0 | 1$\aas@@fstack{\prime\prime}$04
02/10/2022 | D0-G2-J3-K0 | 0$\aas@@fstack{\prime\prime}$56
02/11/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$69
02/13/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$66
03/11/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$55
03/25/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$41
03/25/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$40
Table 2: VLTI/MATISSE observation log. Date | Configuration | Seeing
---|---|---
10/30/2021 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$41
03/02/2022 | A0-G1-J2-J3 | 0$\aas@@fstack{\prime\prime}$45
03/31/2022 | A0-B2-C1-D0 | 0$\aas@@fstack{\prime\prime}$77
Table 3: CHARA observation log. Date | Combiner | Mode | Configuration
---|---|---|---
03/02/2022 | MIRC-X | H-Prism50 | S1-S2-E2-W1-W2
03/03/2022 | MIRC-X | H-Prism50 | S1-S2-E2-W1-W2
03/09/2022 | MIRC-X | H-Prism50 | S1-S2-E2-W1-W2
03/09/2022 | MYSTIC | K-Prism49 | S1-S2-E2-W1-W2
Table 4: Total $LMN$-band fluxes of V838 Mon and the calibrators.

Object |  $L$ band |  $M$ band |  $N$ band | Reference
---|---|---|---|---
| [Jy] | [Jy] | [Jy] |
V838 Mon | 5 | 4 | 30 | 1
HD 52666 | 132 | 67 | 18 | 2
20 Mon | 30.3 | 19.2 | 4.31 | 2

(1) Kamiński et al. (2021); (2) (MDFC; Cruzalebes et al. 2019)

Table 5: Total $K$-band fluxes of V838 Mon and the calibrator.

Object |  $K$ band | Reference
---|---|---
| [mag] |
V838 Mon | 5.08 | 1
HD 54990 | 3.77 | 2

(1) Kamiński et al. (2021); (2) (MDFC; Cruzalebes et al. 2019)

Table 6: Total $H$-band fluxes of V838 Mon and the calibrator.

Object |  $H$ band | Reference
---|---|---
| [mag] |
V838 Mon | 5.86 | 1
HD 54990 | 3.89 | 2
HD 61039 | 5.54 | 2
HD 59230 | 4.16 | 2

(1) Kamiński et al. (2021); (2) (MDFC; Cruzalebes et al. 2019)

Table 7: $H$-band calibrator uniform disk diameters.

Calibrator | UDD, $H$ | UDD, $K$ | eUDD
---|---|---|---
| [mas] | [mas] | [mas]
HD 54990 | 0.863 | 0.867 | 0.079
HD 61039 | 0.377 | 0.379 | 0.009
HD 59230 | 0.749 | 0.753 | 0.072
555Uniform disk diameter (UDD) in the $HK$ bands, taken from the JMMC Stellar
Diameters Catalog (Bourgés et al. 2014). eUDD is the errorbar on the diameter
measurement.

##  3 Geometrical modelling

The resultant data products for the above-mentioned VLTI and CHARA instruments
are the interferometric observables, i.e. squared visibilities ($V^{2}$) and
the closure phases. The squared visibilities represent the fringe contrast,
thus an object is said to be completely resolved in the case that the value
for the squared visibility is zero, and completely unresolved in the case the
value is one. The squared visibilities can be used to constrain the size of
the source. Closure phases are the sum of the individual phase measurements by
telescopes within a particular triangular configuration in the array. This
results in the atmospheric phase cancelling out, thereby only leaving the sum
of intrinsic phases of the object (Jennison 1958). Closure phases are a probe
for asymmetries, so deviations from values of 0° or 180° would indicate some
deviation from centro-symmetry of the source.

###  3.1 $L$-band geometrical modelling (MATISSE)

Previous geometrical modelling and preliminary imaging results in the $L$ band
seem to suggest the presence of an elongated and tilted structure that
enshrouds the central merger remnant (Mobeen et al. 2021). The previous VLTI
observations in 2020 were obtained using only the large configuration. In the
current study, V838 Mon was observed with the small (maximal baseline is
$\approx 30$ m) and large (minimal baseline $\approx$ 140 m) configurations in
the $L$-band. A quick look at the squared visibilities as a function of
spatial frequency suggests that the source is mostly resolved at the longer
baselines (see Fig. 8). However, even with the most extended array
configuration, V838 Mon is not completely resolved as the squared visibilities
drop to a minimum of about 0.3. At the shortest baselines, V838 Mon is mostly
unresolved, with the squared visibilities being at around 0.9. As a function
of wavelength, though, it is clear that the squared visibilities do not show
much variation, which suggests that no major chromatic effects are in play. At
the very edge of the wavebands, the data was much noisier. We excluded it from
our analysis and restricted ourselves to the wavelength range of 3–4 $\mu$m in
the $L$ band.

A qualitative look at the closure phases (see Fig. 8) indicates that they are
mostly close to zero, with some very small deviations of about a few degrees
(maximum $\sim$3°). This suggests a slight asymmetry in the system. Again,
similar to the squared visibilities we do not see much variation in the
closure phases as a function of wavelength, which suggests that the shape of
the object is not changing over the covered wavelength range.

To interpret the MATISSE data, we employed geometrical modelling, because the
UV coverage of the MATISSE 2021/2022 data alone is insufficient for
comprehensive imaging. We used the modelling software LITPRO provided by Jean-
Marie Mariotti Center666https://www.jmmc.fr/english/tools/data-
analysis/oifits-explorer/ (JMMC), which allows the user to fit simple
geometrical models such as disks, point sources and ellipses to the observed
visibilities and closure phases. LITPRO uses a $\chi^{2}$ minimization scheme
to compute the model parameters and their error bars. For the purpose of
fitting the V838 Mon data, we tried four models: a uniform circular disk (CD),
an elliptical disk (ED), a circular Gaussian, and an elliptical Gaussian. With
a reduced $\chi^{2}$ of 6, we find the elliptical models to better represent
the data. The plots of the modelled and observed visibilities are shown in
Fig. 8. Additionally, we are able to obtain precise estimates for the size and
PA for the various models. These are listed in the Table 8. The obtained
parameters for size and orientation are in good agreement with what previous
MATISSE observations had indicated (Mobeen et al. 2021). In Mobeen et al.
(2021), the elliptical disk model yielded a semi-major axis PA of $-40\degr\pm
6\degr$, an angular diameter of $3.28\pm 0.18$ mas, and a stretch ratio of
$1.40\pm 0.1$. In this study, we get values of
–41$\aas@@fstack{\circ}$34$\pm$0$\aas@@fstack{\circ}$49, $3.69\pm 0.02$ mas,
and $1.56\pm 0.01$, respectively. It is evident that these parameters have
remained largely the same between Jan/Mar 2020 and Oct 2021/Mar 2022. This
means that the feature in the $L$ band is non-transient on a timescale of
years. The absence of any change in PA also suggests that no dynamical
variations have occurred in the post merger environment. Thus, it would seem
that the disk-like feature is stable and long-lived in the $L$ band. This
allowed us to combine the two datasets and attempt imaging (see Sect. 4.1).

Figure 8: Observed (blue) and best fits for the $L$ band squared visibilities
and closure phases. The models ED and CD refer to elliptical disk and circular
disk, respectively.

###  3.2 $N$-band geometrical modelling (MATISSE)

As stated previously, despite the generally poor quality of the $N$ band data
due to thermal background and sensitivity limits, we were still able to obtain
the correlated fluxes for V838 Mon for two of the three nights. Being
measurements of the fringe amplitude, correlated fluxes are a proxy for fringe
contrast and allow the determination of the source size. We modelled the $N$
band visibility amplitudes similar to the geometrical modelling we did for the
squared visibilities in other bands. Given the scarce number of data points,
we used simple geometries to model the observables. In the $N$ band, our goal
was to obtain a size measurement for the structure and also constraints on its
ellipticity. We used two models, a uniform disk and an elliptical disk. The
uniform disk model fit yields an angular diameter of $87.73\pm 0.01$ mas,
while the elliptical disk model yields an angular diameter of $48.70\pm 2.40$
mas, a PA of 5$\aas@@fstack{\circ}$03$\pm$0$\aas@@fstack{\circ}$70, and a
stretch ratio of $5.68\pm 0.31$. These fits are shown in Fig. 9. If we
consider the size estimate obtained from the elliptical disk model (semi-major
axis of $\sim 300$ mas), then it would seem that the most extended structure
surrounding V838 Mon exists in the $N$ band. The PA points to a north-south
elongation in the structure, which is drastically different from that seen in
other bands. This means that the $N$ band structure is a distinct feature that
may not be dynamically related to structures seen in the $HKL$ bands.

Figure 9: Observed (blue) and best fits for the $N$ band visibility
amplitudes. The models ED and CD refer to elliptical disk and circular disk,
respectively.

###  3.3 $M$-band geometrical modelling (MATISSE)

For the first time, we were able to observe V838 Mon in the $M$ band using
MATISSE. However, due to the scarce data in this band, we did not attempt any
image reconstruction. As can be seen in Fig. 10, the shape of the squared
visibilities as a function of spatial frequency shows a sinusoidal-like
modulation. It is prominent across all spatial frequencies, but becomes more
noticeable at longer baselines. The amplitude variation lies in the range
0.3–0.4, and the frequency of these modulations increases at lower spatial
frequencies. The $M$-band closure phases are mostly close to zero. However, at
larger spatial frequencies they tend to deviate by a few degrees (maximum
$\sim 5\degr$), thereby hinting towards asymmetries similar to what we see in
the $L$-band.

The visibilities vary also as a function of wavelength and show a broad
depression near 4.7 $\mu$m. This is illustrated in Fig. 11 for visibilities
averaged over all baselines. A minimum is reached at 4.7 $\mu$m which is close
to the fundamental band of CO (Fig. 12). The presence of the feature explains
the sinusoidal variations seen in Fig. 10. Closure phases do not seem to show
much wavelength dependence throughout this particular waveband. Our flux
calibrated spectra, however, readily show an absorption feature at around 4.6
$\mu$m, which we identify as the fundamental band of CO. In Fig. 12, we
compare the entire calibrated $LM$ spectrum to a simulation of CO $\Delta
v=\pm 1$ absorption from an optically thin slab at thermal equilibrium and at
a gas temperature of 30 K. The match in wavelength is not perfect, but can be
explained by the complex pseudo-continuum baseline on both sides of the
feature. The observed feature may not be pure absorption. Also, telluric
absorption of water makes this part of the band particularly difficult to
observe and calibrate. Our identification of the CO band is strengthened by
observations of the same band in V838 Mon 2005 at a much higher spectral
resolution (Geballe et al. 2007), and its presence in other red novae, for
instance in V4332 Sgr (Banerjee et al. 2004).

Although the $L$-band spectrum shows clear signatures of circumstellar
molecular features, limited by the small number of measured visibilities, we
are currently able to make basic size measurement only by combining data in
the entire observed band. We used a wide number of achromatic models, which
include a circular disk (CD), circular Gaussian (CG), elliptical disk (ED),
elliptical Gaussian (EG) and a two point source model (2P). The fits to these
models and their details are given in Table 8, while the models of the
visibilities and closure phases are also shown in Fig. 10. None of the models
are able to reproduce exactly the modulation shown by the visibilities. Disk-
like models with a size of 4–6 mas fit the data much better than a double
source. The derived size most likely represents the dominant circumstellar
component.

Figure 10: Observations (blue) and best fit models for the $M$ band squared
visibilities and closure phases. Models ED, CD and 2P refer to elliptical
disk, circular disk, and two point sources, respectively. Figure 11: Observed
squared visibilities as a function of wavelength near the 4.7 $\mu$m dip. Data
for all baselines were averaged. Figure 12: Observed $LM$-band spectrum (solid
line) is compared to a simple simulation (dashed line) of absorption at the CO
fundamental band near 4.6 $\mu$m and at a gas temperature of 30 K. The
simulated spectrum was smoothed to the spectral resolution of MATISSE.

###  3.4 $K$-band geometrical modelling (GRAVITY and MYSTIC)

The GRAVITY $K$ band observations for V838 Mon are the most extensive of all
the VLTI observations that we have obtained so far. Observations on all three
array configurations, including the intermediate configurations, were
successfully executed, which resulted in good sampling of the UV plane and
allowed us to perform image reconstruction (see Sec. 4.2). We obtained 15 out
of the requested 18 observing runs. Results are shown in Fig. 13. The source
is well resolved at the longest baselines, as the squared visibilities fall to
a minimum of 0.1. The longer baselines also reveal smaller peak-like features,
which point to a complex morphology in the $K$ band that cannot be explained
by just a simple disk or a Gaussian model (shown more clearly in Fig. 14). The
squared visibilities do not change considerably with wavelength. This suggests
that the shape of the structure is not wavelength dependent throughout the
entire $K$ band. We do see, however, subtle spectral effects near 2 and 2.3
$\mu$m which we assign to water and CO absorption (see Figs. 16 and 17). The
closure phases seem to scatter around zero by a few degrees and reach a
maximum of $\sim$5°). The root-mean-square (RMS) scatter is 2° which suggests
that the observed deviations in the closure phases can be attributed mostly to
observational noise (within 4$\sigma$) and maybe some very small asymmetries
present in V838 Mon. They also do not seem to vary with wavelength, thus
indicating a non-variable morphology in the $K$ band.

Although the $K$ band data were used extensively for imaging, we nevertheless
thought it prudent to also model the visibilities directly using simple
geometrical models in order to obtain constraints on the size of the
structure. The best fit parameters for all models are shown in Table 8. The
size of the source lies in the range 1.2–2.6 mas, depending on the choice of
model. The value of the PA in the $K$ band is similar to measurements in the
other bands.

Figure 13: $K$ band squared visibilities and closure phases measured with
GRAVITY (gray). Best-fit geometrical models for elliptical disk (ED) and
circular disk (CD), are also shown (both have zero closure phases). Figure 14:
Zoomed-in section on some of the ‘bump’ features in the squared visibilities
from Fig. 13.

We also were able to obtain $K$ band spectra with a medium resolution
(R$\sim$500). The spectrum was calibrated using the star HD 54990. This
resulted in the removal of the major telluric feature caused by H2O near 2
$\mu$m, leaving us with absorption bands in the range 2.3 – 2.35 $\mu$m and
less intense absorption features near 2 $\mu$m. We identified these as CO and
H2O features, respectively. To ensure that this identification was accurate,
we modelled the absorption at various temperatures using pgopher (Western
2017) and spectroscopic data from the EXOMOL database (Tennyson et al. 2020).
The models represent a simple slab of gas in local thermodynamic equilibrium
(i.e., characterized by a single temperature) and optically thin absorption.
Observed spectra are compared to the simple simulations in Figs. 16, 17, and
18. After experimenting with a wide variety of temperatures for both species,
we find that a model with a CO temperature of 3000 K and an H2O temperature of
1000 K best represent the data, but due to the poor spectral resolution these
are only very rough constraints. The absorption may arise in the photosphere
or very close to it in the circumstellar gas (cf. Geballe et al. 2007). We use
geometrical modelling to determine the sizes of the source when the
visibilities are limited to spectral regions where molecular absorption is
strongest. For CO, we restrict the GRAVITY data to 2.29–2.37 $\mu$m, while for
water 1.97–2.07 $\mu$m was selected. By fitting circular disk models to each
of these datasets, we obtain angular diameters of $2.1\pm 0.1$ mas and $2.0\pm
0.1$ mas, respectively. With the uncertainties, these are consistent with the
size obtained for the entire band.

We also obtained a single night of MYSTIC $K$ band observations and
geometrically modelled the visibilities with a uniform disk and an elliptical
disk model. The results are shown in Fig. 15. The uniform disk model yields an
angular diameter of $1.57\pm 0.01$, mas while the elliptical disk model yields
an angular diameter of $1.72\pm 0.01$ mas, a PA of $-41.6\pm 1.5$°, and a
stretch ratio of $1.28\pm 0.01$. Due to the longer CHARA baselines, V838 Mon
is much better resolved, which yields a smaller angular diameter than that
obtained with just the GRAVITY observations. Additionally, we also combined
both the GRAVITY and CHARA squared visibilities and performed a combined
uniform disk fit to the dataset, which gave an angular diameter estimate of
$1.94\pm 0.01$ mas. This value is essentially the same as the one obtained
from fitting only the GRAVITY squared visibilities, which suggests that
combining them with the MYSTIC squared visibilities does not result in a
significant change in the size estimate for V838 Mon.

The mean of the MYSTIC closure phases is $\sim$2° with an RMS of $\sim$11°,
which strongly hints towards the presence of asymmetries. Unlike the GRAVITY
dataset, the MYSTIC observations lasted only a single night because of which
we could only geometrically model the closure phases. The results of this
modelling are presented in Sect. 4.5.

Figure 15: $K$ band squared visibilities as measured with MYSTIC (blue). Best-
fit geometrical models for elliptical disk (ED) and circular disk (CD) are
also shown. Figure 16: Synthetic spectra generated for H2O at 1000 K and for
CO at 3000 K are compared to the observed $K$ band spectrum from GRAVITY over
the entire observed band. The synthetic spectra are arbitrary scaled. The
simulation does not include high-opacity effects, which may result in the poor
fit to the observed saturated water band near 2 $\mu$m. The synthetic spectra
were convolved to the spectral resolution of the GRAVITY resolution. Figure
17: Synthetic spectra of water generated at various temperatures are compared
to the observed spectrum (blue). Figure 18: Synthetic spectra generated for
the first-overtone CO bands at various temperatures. The solid blue line is
the observed spectrum.

###  3.5 $H$-band geometrical modelling

CHARA MIRC-X observations allowed us to observe V838 Mon in the $H$ band as
well. Previous measurements in the $H$ band were taken in 2012 using the AMBER
instrument, and prior to that in 2004 (Lane et al. 2005) using the Palomar
Testbed Interferometer (PTI). These recent CHARA measurements serve as a
direct follow-up to the previous $H$ band observations. The squared
visibilities, shown in Fig. 19, suggest that the source is resolved at the
longest CHARA baselines ($\lesssim$331 m), as the visibilities fall to a
minimum of 0.1. The closure phases split into two groups. One is scattered
around null phase, and the other, more populous group, is centered at phases
$\sim$6°, thereby hinting towards deviations from centro-symmetry in some
parts of the remnant. This pattern is observed at all spatial frequencies and
requires a more complex source structure, which we discuss in Sect. 4.3.

Figure 19: Observations (dark blue) and best fit models for the $H$ band
squared visibilities and closure phases. Circular (CD) and elliptical (ED)
disk models are shown.

##  4 Image reconstruction

###  4.1 $L$ band image reconstruction

For the purpose of image reconstruction, we combined the 2021/2022 MATISSE $L$
band data set with the one from 2020 used by Mobeen et al. (2021), since no
significant variability occurred in the source over the course of the year.
This is further confirmed by our modelling of the recent visibilities and
phases presented in Sect. 3.1. Since both of the datasets, except for one,
consist of observations made with the large configuration, any reconstructed
image will be sensitive to compact features in V838 Mon. The image was
constructed using the image reconstruction algorithm Multi-aperture Image
Reconstruction Algorithm MIRA (Thiébaut 2008) which uses a likelihood cost
function minimization method. For the $L$ band image, we set a pixel size of
0.1 mas, which is smaller than the theoretical best resolution
($\lambda_{L}/2B_{\rm max}\approx 2$ mas). This was done as an attempt to
reveal any finer substructure in V838 Mon at this band. A FoV of 5 mas was
chosen, taking into consideration the size of the $L$ band structure inferred
from geometrical modelling (see Table 8). The regularization scheme we
employed was the hyperbolic regularization with a regularization parameter of
$10^{6}$. These regularization settings resulted in images with the lowest
$\chi^{2}$ values.

The resulting image displayed in Fig. 20 shows an asymmetrical, elliptical
feature. The geometric orientation is along the PA that we obtained from our
modelling presented in the prior sections.

Figure 20: Image reconstruction of V838 Mon done in the $L$ band using the
combined 2020–2022 MATISSE data. Flux is in a logarithmic scale in arbitrary
units.

###  4.2 $K$ band image reconstructions

Given the comprehensive nature of our GRAVITY observations, we were able to
perform detailed image reconstruction of V838 Mon from the measured
interferometric observables in the $K$ band. For this purpose, we used two
separate image reconstruction algorithms, SQUEEZE (Baron et al. 2010) and MIRA
(Thiébaut 2008). Both these algorithms attempt to solve for the best image
given the data and a model of the image. They employ different methodologies
and minimization techniques from one another. Thus, by using multiple
algorithms to reconstruct an image and comparing them to one another, one can
distinguish between authentic features present in a source and artifacts that
could have resulted from the different algorithms. For the purpose of our
image reconstruction, we limited our FoV to 10 mas and set the pixel size to
0.1 mas, thus making the image 100 by 100 pixels. This choice for FoV and
pixel size was adopted in order to super resolve any potential sub-mas
features. With MIRA, just like in the $L$ band, we used the hyperbolic
regularization and found that the $\chi^{2}$ was minimized with a
regularization parameter of $10^{6}$. While with the SQUEEZE algorithm, we
used total variation and found the $\chi^{2}$ to be minimized with a parameter
value of $10^{3}$. Total variation with a high enough regularization parameter
value is identical to the hyperbolic regularization scheme. The values for the
regularization parameters were determined by reconstructing images for a wide
range of parameter values and by noting the value at which no further
reduction in $\chi^{2}$ takes place. This parameter is significant since it
sets the weight for the image model against the data when the minimization
scheme is implemented. We combined all the GRAVITY OIFITS files and merged
them into a single file by using OI-
Tools777https://github.com/fabienbaron/OITOOLS.jl/tree/main/demos. This
yielded a single OIFITS file that was then used as the input file in the image
algorithms. The number of iterations was set at 500. Finally, a randomly
generated image was used as the starting image for all of our reconstruction
attempts. The reconstruction results with both MIRA and SQUEEZE are shown in
Fig. 21 and 22 respectively. We performed image reconstruction using a
different strategy as well, whereby the combined GRAVITY data OIFITS file was
binned into smaller wavelength increments of 0.05 $\mu$m. The image in the
first wavelength bin in this process was constructed with a random starting
image. However, in subsequent reconstructions for the other bins, the starting
image used would be the image obtained for the preceding bin. Using this chain
imaging method helped us to significantly decrease the value of the reduced
$\chi^{2}$, i.e. from $\approx 150$ to $\approx 8$ in the final image. Thus,
the final image is at the endpoint of the waveband. The result of this
strategy is displayed in Fig. 23.

Figure 21: Reconstructed $K$ band image of V838 Mon obtained using the MIRA
imaging algorithm. Flux level is logarithmic, and the units are arbitrary.
Figure 22: Reconstructed $K$ band image of V838 Mon obtained using the SQUEEZE
imaging algorithm.  Figure 23: Reconstructed $K$ band image of V838 Mon
obtained using the chain imaging method with the MIRA imaging algorithm.

The resulting images (Figs. 21 and 22) from both the algorithms show a similar
morphology in the $K$ band. At 0.1 mas, multiple features are immediately
noticeable. These include the clump-like features (labelled as A, B, D, E, and
G in Fig. 21) that are roughly distributed along an axis at the PA value
($\sim$–40°) which we have obtained from geometrical modelling and which is
similar to the PA found for the disk-like feature found in V838 Mon in other
bands. Furthermore, we also see two narrow linear features (labelled C and F
in Fig. 21) stretching along the northeast direction, which appear to be
nearly perpendicular to the axis defining the brighter clumps. When we look at
the reconstructed image from SQUEEZE, we see that despite being generally
similar to its MIRA counterpart, the structure in the SQUEEZE image is a lot
more continuous, somewhat resembling a clumpy disjointed ring seen at an
intermediate inclination. The features that correspond to regions A, B, D, E,
and G in Fig. 21, are present in Fig. 22, however, features C and F are not
clearly distinguishable in the SQUEEZE image. This could mean that the linear
features are imaging artifacts.

Our highest-resolution (0.1 mas) $K$-band images reveal a very clumpy
morphology that surrounds the central elliptical structure (see Fig. 21). The
major clumps seem to almost trace out a ring that has a size of about 5 mas,
while the inner source spans 2 mas. It is worth noting that at the adopted
pixel size of 0.1 mas, the source is already super-resolved when compared to
the actual achievable resolution with the longest baseline at the VLTI (i.e.
$\sim$1 mas). Thus, the finer details such as the clumps and linear features
could be simply image artefacts that arise as a consequence of using too high
a resolution. To further explore the effects of super-resolution, we performed
image reconstruction using larger pixel sizes of 0.5 and 0.8 mas. These
results can be seen in Figs. 31–34. The resulting less super-resolved images
show only a single extended feature, whereas the inner feature seen in the
more resolved images ceases to be distinguishable. Furthermore, the clumps
which were prominent in the 0.1 mas resolution image also seem to have
vanished, however, the extended structure in the 0.5 mas and 0.8 mas images
seems to suggest that the circumstellar environment is asymmetrical and with a
morphology that could be described as bipolar. The orientation of this bipolar
structure is along the same axis defined by the PA from the geometrical model
fitting. It is likely that the outer structure is indeed an example of
outflows that are expected to be produced in post-merger environments, similar
to what is seen in the sub-mm regime in the stellar merger remnant V4332 Sgr
(Kamiński et al. 2018).

###  4.3 $H$ band image reconstruction

Similar to how we analyzed the $K$ band data, we employed the same imaging
algorithms and attempted image reconstruction for the CHARA observations in
the $H$ band. The image FoV was set to 10 mas, while a pixel size of 0.1 mas
was chosen and the starting image was randomly generated. The value of the
regularization was set to 1000 for both MIRA and SQUEEZE imaging. The
resulting images (Figs. 24 and 25) show an elongated structure that is
oriented at a PA of –40°. The physical extent of this feature is about 2 mas,
thus allowing us to probe the innermost parts of the circumstellar
environment. Despite the fact that the images provide a very good fit to the
data (i.e. reduced $\chi^{2}=1$) some artifacts are also present. These might
have occurred due to the fact that we were only able to use 4–5 telescopes on
any given night and not the full six telescope configuration of CHARA.

It would seem that in the $H$ band too we are seeing a similar bipolar feature
but at much smaller spatial scales, placing it closer to the central merger
remnant. The reconstructed images in the $H$ band show a structure that is
similar in morphology to its counterparts in other wavelength regimes. The
images indicate the presence of asymmetrical structures in the south-east and
north-west regions of V838 Mon. The asymmetry is most noticeable in the
northern part of the bipolar structure, which seems to be misaligned with
respect to the southern feature. A crucial difference in the $H$ band
structure is the relatively high closure phase deviations and the apparent
asymmetrical shape of the northern lobe visible in the image reconstructions.
The orientation of the CHARA feature is in good agreement with what is seen in
the $L$, $M$, and $K$ bands. This can be seen in the sketch of the mid-
infrared structures in Fig. 30.

Figure 24: Reconstructed $H$ band image of V838 Mon obtained using the MIRA
imaging algorithm. Figure 25: Reconstructed $H$ band image of V838 Mon
obtained using the SQUEEZE imaging algorithm.

###  4.4 $H$ band closure phase modelling

As mentioned previously, we were unable to explain $H$ band closure phases
with simple geometrical models (circular and elliptical disks, see Sect. 3.5).
Given the magnitude of these closure phase deviations, it is apparent that the
circumstellar structure in V838 Mon possesses significant asymmetries. This is
further corroborated by our imaging results presented in Sect. 4.3. Similar to
the $K$ band images, the $H$ band ones seem to show a bipolar structure which
is roughly oriented along the same PA as all the other structures that encase
the merger remnant. Given our imaging results in the $H$ band, we again
attempt to reproduce the closure phase deviations, this time using a more
complicated multi-component model that better represents the structure as
revealed by the $H$ band image. A sketch of this model is shown in Fig. 26.
The model includes a central 1 mas disk that represents the central star. In
addition to that, we also include two elliptical disks (“clumps”). The semi-
major axis of the lower ellipse is oriented roughly along the PA that the
structures in other bands follow, while that of the upper ellipse is nearly
horizontal. The modelled and observed closure phases are plotted as a function
of spatial frequency in Fig. 27. We see that the more complex model inspired
by the image reconstructions is able to explain the closure phase deviations
quite satisfactorily.

Figure 26: Model we used for explaining the $H$ band closure phases. The red
central component represents the star, while the blue ellipses may represent
bipolar outflows hinted by the $H$ band image reconstructions (cf. Figs. 24
and 25). About 95% of the model flux is contained within the central
component.  Figure 27: Simulated and observed $H$ band closure phases as a
function of spatial frequency. The simulated closure phases are represented by
red triangles, while the observations are represented by the blue circles.
Simulations correspond to the structure shown in Fig. 26 and match the
observations better than the simpler models described in Sect. 3.5.

###  4.5 MYSTIC $K$ band closure phase modelling

As stated before, the $K$ band closure phases as observed by MYSTIC do show
significant non-zero deviations. Therefore, we decided to model these closure
phases geometrically to see if the observations can be well explained with an
asymmetric intensity distribution. We used the interferometric modelling
software PMOIRED (Mérand 2022) for this purpose. Using PMOIRED, we
experimented with various combinations of geometrical building blocks (disks,
Gaussians, etc.) in order to reproduce the closure phases as accurately as
possible. Ultimately, we found the best fit model ($\chi^{2}\\!\sim$1.5) which
consists of an elongated primary disk with a size of 1 mas and a smaller (0.5
mas) less intense secondary disk that is separated from the primary by
$\sim$1.7 mas. The geometry of this model is shown in Fig. 28.

Figure 28: Image for the best-fit geometric model comprising two components in
the $K$ band. The primary component is the elliptical disk with a size of
$\sim$1 mas and elongated at a PA of $-40\degr$. The secondary component is a
smaller circular disk with a diameter of 0.5 mas and separated from the
primary by $\sim$1.7 mas. The bulk of the model flux ($\sim$ 99%) is contained
in the primary component.

The closure phases produced by the model and their fit to the observations are
shown in Fig. 29. We can see from the figure that the non-zero closure phases
are explained quite satisfactorily by our model at all baselines. The
orientation of the primary disk seems to agree well with the general
orientation of V838 Mon as seen in the $HK$ band image reconstructions.
Furthermore, the PA of the smaller disk component with respect to the primary
also closely resembles the general orientation of V838 Mon in the near to mid
infrared bands. The modelling results prove that the closure phase deviations
are well explained with the addition of an off-centre “clump”.

Figure 29: MYSTIC $K$-band closure phases (top) for the observed (points) and
modelled (lines) visibilities representing the model in Fig. 28. The bottom
panel shows the respective residuals. Figure 30: A sketch of the best fitting
elliptical models for V838 Mon in the $HKL$ bands as observed by MIRC-X,
GRAVITY and MATISSE.

##  5 Discussion

###  5.1 Jets

Lane et al. (2005) suggested in their study that V838 Mon could potentially
possess an asymmetric structure, however due to their limited observations and
absence of closure phase measurements they were unable to prove this
conclusively. In Mobeen et al. (2021), the $L$ band closure phases seemed to
be hinting towards minor asymmetries, as indicated by the small but non-zero
deviations. However, we interpreted them with great caution given their small
magnitude and the limited data. In this study, we have conclusively shown the
existence of these asymmetries via geometrical modelling and image
reconstruction in the $HK$ bands. The GRAVITY and CHARA images (see Figs. 21
and 24) show an inhomogeneous bipolar structure. This inhomogeneous structure
is more pronounced in the $H$ band image, as we can see that the northern lobe
is oriented at a different PA than its southern counterpart. The asymmetric
lobes model is further corroborated by $H$ band closure phases that are
satisfactorily explained when we assume such a morphology for our geometrical
model (see Figs. 26 and 27).

These outflows resemble jets, which are strongly believed to be present in
Intermediate Luminosity Optical Transients (ILOT) such as Luminous Red Novae
(LRN) and are thought to play an important role by providing the extra energy
needed to unbind the common envelope. Soker (2016) shows that energy from the
jets is transferred to the primary envelope during the grazing envelope
evolution stage, which means that the jets could be the mechanism responsible
for the outbursts seen in LRNs including V838 Mon. In a subsequent study,
(Soker 2020) it was shown that the luminosity and total energy output of V838
Mon could be adequately explained if jets present in the system collide with a
slow moving shell of material. The author concluded that this was a plausible
scenario for V838 Mon and V4332 Sgr. In the case of V4332 Sgr it was shown
using sub millimetre observations carried out at ALMA (Kamiński et al. 2018)
that the object does indeed possess bipolar outflows consisting of molecular
material. This means that jets were a main component in the LRN event in V4332
Sgr. If jets are the reason for the LRN in V838 Mon, then that necessitates
that we should be able to observe jet like features in the source. Since we
see a clear bipolar structure in the image reconstructions for V838 Mon
presented in Sects. 4.2 and 4.3 this means that jets played a role in the
energy output of the outburst in 2002 which resulted in V838 Mon brightening
by several orders of magnitude. Soker (2023) shows that only jets can
satisfactorily explain the rapid rise in luminosity that is seen in ILOTS such
as V838 Mon and V1309 Sco just prior to the outbursts. As of now, no jets have
been observed in V1309 Sco, but the clear detection of jets in V838 Mon, V4332
Sgr and CK Vulpeculae suggests that such bipolar features are perhaps
universal amongst many other stellar merger remnants. Recently the Blue Ring
Nebula was found to possess bipolar conical outflows which led Hoadley et al.
(2020) to conclude that the object is a result of a stellar merger. Our
interferometric observations of V838 Mon seem to be giving credence to the
idea that jets are an intrinsic feature of LRNs and persist even in the post-
merger environment. In V838 Mon it seems that the jets are only prominent in
the mid-infrared $HK$ bands since the ALMA maps show the merger ejecta to be
spherically symmetric at a spatial scale of $\sim$ 25 mas and any jet-like
feature seems to be absent (Kamiński et al. 2021). This contrasts with the
ALMA maps for V4332 Sgr which show a very clear signature of a bipolar outflow
(Kamiński et al. 2018). With extensive future VLTI observations of V4332 Sgr
we will be able to determine the extent to which the milliarcsecond scale
environment differs from the large scale structures as revealed by ALMA. It
could be that the nature of every LRN jet is unique and is determined by a
number of other factors such as the mass and stellar type of the progenitor
stars, the surrounding medium and duration of the outburst, etc.

###  5.2 Clumps

The super resolved GRAVITY $K$ band image (see Fig. 21) points to the possible
presence of clumps in the immediate circumstellar vicinity in V838 Mon. The
GRAVITY closure phases alone do not seem to give us much insight, as the
deviations are much too small (see lower panel in Fig. 13). However, since the
CHARA baselines are longer by a factor of about three, the closure phases
measured using MYSTIC show much higher deviations from $0\degr$ (discussed in
Sect. 4.5). These MYSTIC closure phases, as evidenced by the modelling,
suggest the presence of a small clump separated from the central star by $\sim
1.7$ mas. The size of the clump is 0.5 mas, which makes it too small to be
detected by the VLTI but large enough to appear in our CHARA observations.

The clump lies approximately along the PA traced by the bipolar structure in
V838 Mon, in fact its orientation seems to be quite similar to that of the
northern lobe seen in the $H$ band. This clump is the first ever detection of
a sub milliarcsecond feature in V838 Mon. The exact nature of the clump,
whether stable or transient as of now, remains elusive. ALMA maps of V838 Mon
by Kamiński et al. (2021) also showed the presence of prominent clumps in the
wind component of the merger remnant in the submillimeter regime. The
formation of dust clumps is a prevalent phenomenon amongst red supergiants
(RSGs) and asymptotic giant branch stars. VX Sagitarri which is a variable
supergiant star has been found to possess an inhomogeneous atmosphere which
contains spots (Chiavassa et al. 2010). It was shown by Chiavassa et al.
(2009) using radiative hydrodynamical simulations that the interferometric
observations for the RSG $\alpha$ Ori could be explained using limb darkening
models, which suggested the presence of asymmetric surface features. In
addition to this, Kamiński (2019) showed using ALMA observations that multiple
dust clumps had formed in the RSG VY Canis Majoris (VY CMa). A dust clump is
also thought to have formed in the (RSG) Betelgeuse during the Great Dimming
event of 2019–2020 (Montargès et al. 2021). Follow-up VLTI/MATISSE
observations in 2020 were obtained by Cannon et al. (2023) which showed a
complex closure phase signal in Betelgeuse. The authors constructed multiple
radiative transfer models, such as those consisting of clumps and spots.
Although they could not reproduce the closure phase signal exactly due to
their use of only rudimentary models, they concluded that either clumps or
other surface features were the main source of the closure phase deviations.
In the case of V838 Mon, the closure phases are not as high as those in
Betelgeuse observed by Cannon et al. (2023), but, as shown in Sect. 4.5, the
MYSTIC closure phases can only be explained by the smaller clump (see Fig.
28).

Since the spectra for V838 Mon resembles that of a late M-type supergiant,
(Evans et al. 2003; Loebman et al. 2015) we expect the circumstellar
environment in V838 Mon to be somewhat similar to other evolved stars with
clumpy surroundings. This particular clump that we have been able to constrain
with MYSTIC could in fact be a smaller feature within the lobe/jet like
structure that is observed at longer wavelengths. This would mean that the
lobes themselves could be clumpy and inhomogeneous. These clumps could be the
result of the jet interacting with the surrounding less dense material in the
V838 Mon during its common envelope phase. Soker (2019) propose that this
interaction between the jet and the envelope results in the formation of jet-
inflated low-density bubbles. As a consequence, the surrounding medium is
prone to clumping due to Rayleigh-Taylor instabilities. This would result in
the jet departing from a pure axisymmetric morphology and become variable in
direction and brightness. It could be the case that the observed $H$ band
clump was the result of such jet-inflated bubbles.

If the clump is a transient feature, then multiple follow-up observations with
MIRCX/MYSTIC at CHARA could help to trace out its trajectory with respect to
the central merger product. These observations could also help determine the
time evolution of the size and brightness of the clump, which can help to
determine whether the jet-inflated bubble formation scenario is plausible.

###  5.3 Clump ejection scenario

The outburst in 2002 showed a complex multi-peaked light curve that was
explained to be the result of multiple shell ejection events in V838 Mon.
These short duration bursts of energy are thought to coincide with the
ejection of the shells and material outflows (Tylenda 2005; Soker & Tylenda
2007). Recent $V$ band photometry of V838 Mon also shows that the source has
been gradually increasing in brightness over the past decade, especially in
the $UBV$ bands (Liimets et al. 2023). The photometry also shows amplitude
variations of about $\sim$0.4, which occur with an approximate period of
$\sim$300 days 888https://asas-sn.osu.edu.

If the observed clump is a mass ejection event or part of the bipolar outflows
seen in the $H$ band images (see Fig. 24), then it could very likely be the
case that the formation of such features correlates with the small brightness
increases that are seen in the light curve modulations for V838 Mon. If true,
then this would mean that the post-merger circumstellar environment is quite
dynamic and is evolving over time scales of about a year. This is in stark
contrast to the $L$ band structure that appears to be stable and non-transient
over timescales of years. Future extensive observations of V838 Mon with
MYSTIC/MIRC-X could help to determine if the observed amplitude variations
correlate with the appearance of ’clumps’, which can help to constrain the
source of the brightness fluctuations observed in V838 Mon.

###  5.4 Contraction controversy

Previously, it was claimed that V838 Mon was gradually shrinking due to the
geometrical thinning of the surrounding ejecta (Chesneau et al. 2014). Our
recent size measurement in the $K$ band shows that contraction has apparently
”ceased” and the size of V838 Mon at 2 $\mu$m is now greater than that
measured by Lane et al. (2005) just a few years after the outburst. This would
imply that V838 Mon instead has now started expanding since the previous size
measurement by Chesneau et al. (2014). Previous interferometric observations
by Lane et al. (2005) and Chesneau et al. (2014) were limited in scope due to
a wide variety of factors.

For example, the observations by Lane et al. (2005) were done using the PTI
which had only two baselines with a maximum baseline length of 85 m. As a
result of this, V838 Mon was not properly resolved, as is evidenced by the
fact that the lowest squared visibility measured by Lane et al. (2005) was
about 0.6. Similarly, the VLTI-AMBER squared visibilities at 2 $\mu$m fitted
by Chesneau et al. (2014) reached a minimum of about 0.2, which is a
significant improvement over the previous measurement, but a small fraction of
the source flux remained unresolved. This also explains why the size
measurement obtained by Chesneau et al. (2014) is smaller than that estimated
by Lane et al. (2005). Their fit yielded a size of 1.15 $\pm$ 0.20 mas, which
is about $40\%$ less than the value measured by Lane et al. (2005) (1.83 $\pm$
0.06 mas), which seemed to have indicated that a contraction had taken place
over the course of a decade.

Our GRAVITY and MYSTIC $K$ band observations on the other hand were able to
better resolve V838 Mon due to the longer baselines (maximum baselines of
$\sim$ 130 m and $\sim$ 300 m on the VLTI and CHARA respectively) and greater
instrument sensitivity. The smallest observed value for $V^{2}$ with VLTI was
about 0.1 while with CHARA it was $\sim$ zero. These size estimates for V838
Mon are based on observations that almost fully resolve V838 Mon, which was
not the case in the above-mentioned previous studies. For example the uniform
disk diameter measured with CHARA-MYSTIC as stated earlier is $1.57\pm 0.01$
mas while with VLTI-GRAVITY it is $1.939\pm 0.001$ mas (see Table 8), since
the CHARA baselines are longer than the VLTI ones by about a factor of three
we conclude that the CHARA size estimate is indeed the most accurate. If the
initial size estimate by Lane et al. (2005) is assumed to be accurate, then
our results imply that indeed V838 Mon has begun to expand. However, as we
have discussed above, the 2005 measurements were limited by baseline length,
and therefore it is highly likely that the angular diameter for V838 Mon was
overestimated. This would explain the apparent contraction that was believed
to be occurring in V838 Mon. Future observations with CHARA can help to
determine the actual evolution of the size of V838 Mon by comparing them to
the previous measurement from 2013 and the ones presented in this study.

###  5.5 Polarization in V838 Mon

The Gemini speckle image reconstructions in the $V$ and $I$ band speckle
images show a morphology that is consistent with the interferometric imaging
and modelling results in the mid-infrared. Both images (see Figs. 6 and 7)
seem to suggest that the innermost vicinity in V838 Mon is quite noticeably
bright at 562 nm and 832 nm. If this emission results from scattered light,
then it is likely to produce a net polarization. It was shown by Wisniewski et
al. (2003b) that a low degree of polarization (1$\%$) was present in 2003 in
the $V$ band, with a measured PA of $150\degr$. The authors suggested that
this polarization signal was pointing to asymmetries in the surroundings of
V838 Mon. A follow-up measurement, six months later, showed that the PA of the
polarisation vector flipped by $90\degr$ (Wisniewski et al. 2003a). This
change over the course of six months seemed to suggest the dynamic, evolving
nature of the dusty environment in V838 Mon. The Gemini observations reveal a
complex circumstellar medium in V838 Mon that might still be dynamic and could
show signs of varying intrinsic polarization and PA, especially in the $V$
band. The above-mentioned observations by Wisniewski et al. (2003b) and
Wisniewski et al. (2003a) were carried out over two decades ago, and since the
dramatic change in PA occurred over a period of only six months, we expect the
environment in V838 Mon to have significantly changed since then. Future
polarimetric measurements could help to determine if linear polarization in
V838 Mon has increased significantly in the last two decades. As we discuss
earlier in the section, it could be possible that the observed jet like
features in the $K$ band could be transient phenomenon, in which case it would
be important to study how the PA and degree of polarization change alongside
the jets. We also note that the clumping scenario that we have discussed
before could also be probed with further polarimetric observations. For
instance, if the clump like feature that we see in our MYSTIC $K$ band
modelling results is indeed a mass ejection event, then it could cause
significant changes in the degree of polarization and the orientation of the
PA as it propagates through the dusty circumstellar environment.

These polarimetric observations can be combined alongside our interferometric
results in this paper for radiative transfer modelling of the dusty
environment in V838 Mon. Synthetic images at the $HK$ bands can be produced
from these simulations, which can then be compared to the reconstructed images
that we presented in this study.

##  6 Conclusions

In this study, we were able to perform a multi-wavelength interferometric
analysis of the stellar merger remnant V838 Mon in the $HKLM$ bands using the
VLTI instruments MATISSE and GRAVITY as well as the CHARA instrument
MIRC-X/MYSTIC. In the $HK$ bands we imaged the circumstellar environment in
the immediate vicinity of the central merger remnant, while in the other bands
we resorted to geometrical modelling to obtain constraints on the geometry and
orientation of the post-merger environment. We were also able to compare the
mid-infrared image reconstructions with $I$ and $V$ band speckle imaging from
Gemini South. We also modelled the spectral features that were present in the
$K$ and $M$ bands and were able to put constrains on the temperatures. Our
main findings are as follows.

  * •

Geometrical modelling suggests that the $L$ band structure has not varied at
all since previous observations presented in Mobeen et al. (2021).

  * •

Image reconstruction in the $HK$ bands reveals a bipolar morphology that
resembles jets. These ’jets’ are present in both bands, while in the $H$ band
they appear to be slightly asymmetric as well.

  * •

The super resolved GRAVITY $K$ band images at 0.1 mas reveal clumpy outflows
that surround the inner feature.

  * •

The MYSTIC $K$ band closure phases are well explained by a small (0.5 mas),
off-centre circular feature that we identify as a potential clump.

  * •

The orientation of the extended circumstellar feature (i.e. the outflows)
tends to be along the same general direction i.e. north-west, with a PA that
varies from –30$\degr$ to –50$\degr$ and an ellipticity in the range 1.1–1.6,
depending on the wavelength, which is in agreement with prior studies
(Chesneau et al. 2014; Kamiński et al. 2021). The $N$ and $I$ band structure
however is oriented along the north-south direction

  * •

The $K$ band CO and water features are best modelled at temperatures of 2000 –
3000 K, thus hinting at their photospheric origin. We were also able to
measure the sizes of the CO and water emitting regions, which are 2.1 mas and
1.985 mas respectively. In the $M$ band, we observe CO for the first time at
4.6 $\mu$m for the first time using interferometry.

###### Acknowledgements.

T. K. and M. Z. M. acknowledge funding from grant no 2018/30/E/ST9/00398 from
the Polish National Science Center. Based on observations made with ESO
telescopes at Paranal observatory under program IDs 0104.D-0101(C) and
0108.D-0628(D). This research has benefited from the help of SUV, the VLTI
user support service of the Jean-Marie Mariotti Center
(http://www.jmmc.fr/suv.htm). This research has also made use of the JMMC’s
Searchcal, LITpro, OIFitsExplorer, and Aspro services. This work is based upon
observations obtained with the Georgia State University Center for High
Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA
Array is supported by the National Science Foundation under Grant No.
AST-1636624 and AST-2034336. Institutional support has been provided from the
GSU College of Arts and Sciences and the GSU Office of the Vice President for
Research and Economic Development. Time at the CHARA Array was granted through
the NOIRLab community access program (NOIRLab PropID: 2022A-426176; PI: M.
Mobeen). SK acknowledges funding for MIRC-X from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (Starting Grant No. 639889 and Consolidated Grant No. 101003096).
JDM acknowledges funding for the development of MIRC-X (NASA-XRP NNX16AD43G,
NSF-AST 1909165) and MYSTIC (NSF-ATI 1506540, NSF-AST 1909165). Some of the
observations in this paper made use of the High-Resolution Imaging instrument
Zorro and were obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105.
Zorro was funded by the NASA Exoplanet Exploration Program and built at the
NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and
Emmett Quigley. Zorro was mounted on the Gemini South telescope of the
international Gemini Observatory, a program of NSF’s OIR Lab, which is managed
by the Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation. on behalf of the
Gemini partnership: the National Science Foundation (United States), National
Research Council (Canada), Agencia Nacional de Investigación y Desarrollo
(Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina),
Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and
Korea Astronomy and Space Science Institute (Republic of Korea).

## References

  * Anugu et al. (2020) Anugu, N., Le Bouquin, J.-B., Monnier, J. D., et al. 2020, AJ, 160, 158
  * Banerjee & Ashok (2002) Banerjee, D. P. K. & Ashok, N. M. 2002, A&A, 395, 161
  * Banerjee et al. (2004) Banerjee, D. P. K., Varricatt, W. P., & Ashok, N. M. 2004, ApJ, 615, L53
  * Baron et al. (2010) Baron, F., Monnier, J. D., & Kloppenborg, B. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7734, Optical and Infrared Interferometry II, ed. W. C. Danchi, F. Delplancke, & J. K. Rajagopal, 77342I
  * Bond et al. (2003) Bond, H. E., Henden, A., Levay, Z. G., et al. 2003, Nature, 422, 405
  * Bourgés et al. (2014) Bourgés, L., Lafrasse, S., Mella, G., et al. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay, 223
  * Bourgés et al. (2017) Bourgés, L., Mella, G., Lafrasse, S., et al. 2017, VizieR Online Data Catalog, II/346
  * Cannon et al. (2023) Cannon, E., Montargès, M., de Koter, A., et al. 2023, arXiv e-prints, arXiv:2303.08892
  * Chesneau et al. (2014) Chesneau, O., Millour, F., De Marco, O., et al. 2014, A&A, 569, L3
  * Chiavassa et al. (2010) Chiavassa, A., Lacour, S., Millour, F., et al. 2010, A&A, 511, A51
  * Chiavassa et al. (2009) Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351
  * Cruzalebes et al. (2019) Cruzalebes, P., Petrov, R. G. P., Robbe-Dubois, S., et al. 2019, VizieR Online Data Catalog, II/361
  * Evans et al. (2003) Evans, A., Geballe, T. R., Rushton, M. T., et al. 2003, MNRAS, 343, 1054
  * Geballe et al. (2007) Geballe, T. R., Rushton, M. T., Eyres, S. P. S., et al. 2007, A&A, 467, 269
  * GRAVITY Collaboration et al. (2017) GRAVITY Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A, 602, A94
  * Hoadley et al. (2020) Hoadley, K., Martin, D. C., Metzger, B. D., et al. 2020, Nature, 587, 387
  * Howell et al. (2011) Howell, S. B., Everett, M. E., Sherry, W., Horch, E., & Ciardi, D. R. 2011, AJ, 142, 19
  * Howell & Furlan (2022) Howell, S. B. & Furlan, E. 2022, Frontiers in Astronomy and Space Sciences, 9, 871163
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Jennison (1958) Jennison, R. C. 1958, MNRAS, 118, 276
  * Kamiński (2019) Kamiński, T. 2019, A&A, 627, A114
  * Kamiński et al. (2009) Kamiński, T., Schmidt, M., Tylenda, R., Konacki, M., & Gromadzki, M. 2009, ApJS, 182, 33
  * Kamiński et al. (2018) Kamiński, T., Steffen, W., Tylenda, R., et al. 2018, A&A, 617, A129
  * Kamiński et al. (2021) Kamiński, T., Tylenda, R., Kiljan, A., et al. 2021, A&A, 655, A32
  * Lane et al. (2005) Lane, B. F., Retter, A., Thompson, R. R., & Eisner, J. A. 2005, ApJ, 622, L137
  * Leinert et al. (2003) Leinert, C., Graser, U., Przygodda, F., et al. 2003, Ap&SS, 286, 73
  * Liimets et al. (2023) Liimets, T., Kolka, I., Kraus, M., et al. 2023, A&A, 670, A13
  * Loebman et al. (2015) Loebman, S. R., Wisniewski, J. P., Schmidt, S. J., et al. 2015, AJ, 149, 17
  * Lopez et al. (2022) Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192
  * Mérand (2022) Mérand, A. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12183, Optical and Infrared Interferometry and Imaging VIII, ed. A. Mérand, S. Sallum, & J. Sanchez-Bermudez, 121831N
  * Mobeen et al. (2021) Mobeen, M. Z., Kamiński, T., Matter, A., Wittkowski, M., & Paladini, C. 2021, A&A, 655, A100
  * Monnier et al. (2018) Monnier, J. D., Le Bouquin, J.-B., Anugu, N., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10701, Optical and Infrared Interferometry and Imaging VI, ed. M. J. Creech-Eakman, P. G. Tuthill, & A. Mérand, 1070122
  * Montargès et al. (2021) Montargès, M., Cannon, E., Lagadec, E., et al. 2021, Nature, 594, 365
  * Munari et al. (2002a) Munari, U., Desidera, S., & Henden, A. 2002a, IAU Circ., 8005, 2
  * Munari et al. (2002b) Munari, U., Henden, A., Kiyota, S., et al. 2002b, A&A, 389, L51
  * Pastorello et al. (2019) Pastorello, A., Mason, E., Taubenberger, S., et al. 2019, A&A, 630, A75
  * Pejcha et al. (2017) Pejcha, O., Metzger, B. D., Tyles, J. G., & Tomida, K. 2017, ApJ, 850, 59
  * Petrov et al. (2007) Petrov, R. G., Malbet, F., Weigelt, G., et al. 2007, A&A, 464, 1
  * Scott et al. (2021) Scott, N. J., Howell, S. B., Gnilka, C. L., et al. 2021, Frontiers in Astronomy and Space Sciences, 8, 138
  * Setterholm et al. (2022) Setterholm, B. R., Monnier, J. D., Le Bouquin, J.-B., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12183, Optical and Infrared Interferometry and Imaging VIII, ed. A. Mérand, S. Sallum, & J. Sanchez-Bermudez, 121830B
  * Soker (2016) Soker, N. 2016, New Astronomy, 47, 16
  * Soker (2019) Soker, N. 2019, MNRAS, 483, 5020
  * Soker (2020) Soker, N. 2020, ApJ, 893, 20
  * Soker (2023) Soker, N. 2023, arXiv e-prints, arXiv:2306.07702
  * Soker & Tylenda (2007) Soker, N. & Tylenda, R. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 363, The Nature of V838 Mon and its Light Echo, ed. R. L. M. Corradi & U. Munari, 280
  * Sparks et al. (2008) Sparks, W. B., Bond, H. E., Cracraft, M., et al. 2008, AJ, 135, 605
  * ten Brummelaar et al. (2005) ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453
  * Tennyson et al. (2020) Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2020, J. Quant. Spec. Radiat. Transf., 255, 107228
  * Thiébaut (2008) Thiébaut, E. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7013, Optical and Infrared Interferometry, ed. M. Schöller, W. C. Danchi, & F. Delplancke, 70131I
  * Tylenda (2005) Tylenda, R. 2005, A&A, 436, 1009
  * Tylenda et al. (2009) Tylenda, R., Kamiński, T., & Schmidt, M. 2009, A&A, 503, 899
  * Tylenda & Soker (2006) Tylenda, R. & Soker, N. 2006, A&A, 451, 223
  * Webbink (1976) Webbink, R. F. 1976, ApJ, 209, 829
  * Western (2017) Western, C. M. 2017, Journal of Quantitative Spectroscopy and Radiative Transfer, 186, 221, satellite Remote Sensing and Spectroscopy: Joint ACE-Odin Meeting, October 2015
  * Wisniewski et al. (2003a) Wisniewski, J. P., Bjorkman, K. S., & Magalhães, A. M. 2003a, ApJ, 598, L43
  * Wisniewski et al. (2008) Wisniewski, J. P., Clampin, M., Bjorkman, K. S., & Barry, R. K. 2008, ApJ, 683, L171
  * Wisniewski et al. (2003b) Wisniewski, J. P., Morrison, N. D., Bjorkman, K. S., et al. 2003b, ApJ, 588, 486

##  Appendix A Images at lower resolutions

Figure 31: $H$ band image for V838 Mon with MIRA algorithm at a resolution of
0.2 mas. Figure 32: $H$ band image for V838 Mon with MIRA algorithm at a
resolution of 0.5 mas.  Figure 33: $K$ band image for V838 Mon with MIRA
algorithm at a resolution of 0.5 mas.  Figure 34: $K$ band image for V838 Mon
with MIRA algorithm at a resolution of 0.8 mas.  Table 8: Fitted models and
their parameters.

Model | Band | $\chi^{2}_{\rm r}$ | Size | $PA$ | Stretch | Flux |
---|---|---|---|---|---|---|---
|  |  | [mas] | [°] | ratio | [%] |
CD | $L$ | 10.34 | 3.00 $\pm$ 0.01 | - | - | $100^{a}$ |
ED | $L$ | 4.22 | 3.69 $\pm$ 0.02 | –41.34 $\pm$ 0.49 | 1.56 $\pm$ 0.01 | $100^{a}$ |
CG | $L$ | 9.50 | 1.82 $\pm$ 0.01 | - | - | $100^{a}$ |
EG | $L$ | 4.02 | 1.44 $\pm$ 0.01 | –41.56 $\pm$ 0.51 | 1.54 $\pm$ 0.2 | $100^{a}$ |
CD | $M$ | 21.94 | 5.86 $\pm$ 0.08 | - | - | $100^{a}$ |
EG | $M$ | 18.75 | 4.35 $\pm$ 0.09 | –53 $\pm$ 2.87 | 1.42 $\pm$ 0.03 | $100^{a}$ |
CD | $H$ | 16.92 | 1.180 $\pm$ 0.004 | - | - | $100^{a}$ |
ED | $H$ | 16.66 | 1.20 $\pm$ 0.01 | –30.69 $\pm$ 3.91 | 1.10 $\pm$ 0.01 | $100^{a}$ |
CD | $K$ | 341.5 | 1.939 $\pm$ 0.001 | - | - | $100^{a}$ |
ED | $K$ | 195.5 | 2.619 $\pm$ 0.002 | –47.09 $\pm$ 0.06 | 0.615 $\pm$ 0.01 | $100^{a}$ |
CG | $K$ | 328.1 | 1.179 $\pm$ 0.001 | - | - | $100^{a}$ |
999The models used to represent the $L$, $M$, $H$ and $K$ band data taken with
MATISSE, GRAVITY and MIRC-X. The parameters are the size (angular diameter or
FWHM), the PA of the major axis, the stretch ratio and the flux level. The
models used are uniform circular or elliptical disk (CD or ED, respectively),
a circular or elliptical Gaussian (CG or EG, respectively). The stretch ratio
is between that of the major and minor axis size. $a$$a$footnotetext: Flux
value was fixed from the start. The size measurements in the $K$ band were
obtained using GRAVITY observations only.