© The Authors 2023

1 Introduction

In late 2019 and early 2020, the red supergiant (RSG), Betelgeuse (α Ori), experienced an unprecedented dimming event (Guinan et al. 2020) with its brightness in the V band diminishing by approximately 70% at its minimum in February 2020. Multiple hypotheses regarding the physical processes responsible for this great drop in brightness have been proposed. These include dust nucleation above the surface obscuring the star (Safonov et al. 2020; Cotton et al. 2020; Levesque & Massey 2020), an ejection of dense chromospheric material from the southern hemisphere subsequently forming molecules and dust (Dupree et al. 2020), photospheric cooling as a result of stochastic convective motions and pulsations (Dharmawardena et al. 2020; Harper et al. 2020; Alexeeva et al. 2021), and a ‘molecular plume’ causing a local opacity increase (Kravchenko et al. 2021; Davies & Plez 2021).

To further explore and differentiate between these scenarios, spatially resolved observations of the stellar disk and direct surroundings of the star were needed. Observations with the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument at the Very Large Telescope (VLT) from Montargès et al. (2021) provided such data and revealed that the dimming was most apparent in the southern hemisphere of the star. Based on 3D dust radiative transfer modelling and fitting of the spectral energy distribution (SED), they conclude that their data of the Great Dimming are in agreement with both a localised photospheric cooling and the formation of dust concentrated in the line of sight towards the southern part of the star, or the (possibly additive) effects of both phenomena. This hypothesis is further strengthened by the recent study by Taniguchi et al. (2022) who come to the same conclusion using photometric monitoring of Betelgeuse from the Himawari-8 meteorological satellite.

Understanding the cause of the Great Dimming of Betelgeuse may prove vital in unravelling the mystery behind the mechanism(s) driving copious amounts of mass from RSGs (Ṁ ~ 10⁻⁷−10⁻⁵ M_⊙ yr⁻¹; e.g. De Beck et al. 2010; Beasor et al. 2020). Gustafsson & Plez (1992) and Josselin & Plez (2007) propose the mass-loss trigger is linked to surface activity, where turbulent motions, pulsations, and large convective cells up-welling from the stellar envelope lower the effective gravity, potentially allowing radiation pressure to launch material. Large convective cells have been observed, detected as hot and cool spots on the stellar surface, on Betelgeuse (e.g. Haubois et al. 2009; Chiavassa et al. 2010; Montargès et al. 2016; O’Gorman et al. 2017) using optical and radio interferometry. Kee et al. (2021) pioneered a theoretical mass-loss prescription in which the wind is driven by turbulent pressure. Such turbulence may be connected to (sub-)surface activity, specifically to pulsations and the violent surfacing of large convective cells. If indeed it is possible to link surface activity to the ejection of material from the star and the consequent formation of clouds of dust, it could bring us one step closer to understanding the (or at least one) wind-launching mechanism of RSGs. These stars are the most frequent progenitors of type II-P supernovæ (SNe). Recent observations (Moriya et al. 2018; Jacobson-Galán et al. 2022) have shown how mass loss during the RSG phase can modify the SN light curve. Close-by RSGs such as Betelgeuse are ideal candidates to search for such observational constraints due to their large angular diameters allowing for spatial resolution of the star.

In this paper, we present VLT Interferometer (VLTI) N-band observations of Betelgeuse with the Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE) instrument, obtained near the V-band brightness minimum, in February 2020. We tested the compatibility of the cool spot and dust clump scenarios with these data. To this end, we performed 3D radiative transfer modelling of the star and near surroundings. In Sect. 2 we discuss the acquisition of the observations and the data reduction process. We then present the observational analysis in Sect. 3. The setup and analysis of the dust radiative transfer models are described in Sect. 4, and Sect. 5 shows the same for our cool spot model. We discuss our findings in Sect. 6 and finish with a summary and conclusions in Sect. 7.

2 Observations and data reduction

MATISSE (Lopez et al. 2014, 2022) is a four-telescope beam combiner instrument at ESO’s VLTI, which covers the mid-infrared range from 2.8 to 13 µm. Observations of Betelgeuse were taken with the small, medium, and large configurations of the 1.8 m auxiliary telescopes (ATs) during three separate nights in February 2020. Each observation of Betelgeuse was preceded and/or followed by a spectro-interferometric calibrator observation. Sirius and Epsilon Leporis (ϵ Lep) were used for that purpose. For this study we focused on the N-band (8– 13 µm) data, which were obtained in the high spectral resolution (R ~ 220) mode of MATISSE. The log of these observations is shown in Table 1.

The data were reduced with version 1.5.5 of the MATISSE pipeline¹. In this work, we did not aim to study the accurate dust mineralogy of Betelgeuse. Hence, to optimise the signal-to-noise ratio on our data, we thus rebinned them to a lower spectral resolution (R ~ 30). This step was performed in the pipeline on the fringe data before the computation of the interferometric observables.

The pipeline outputs are in the form of OIFITS files (version 2), which contain the uncalibrated interferometric observables. Six spectrally dispersed squared visibilities, three independent closure phases (out of four in total), and four total spectra (from the four telescopes) per exposure are included.

To calibrate the observations, we used two calibrators, Sirius and ϵ Lep (see Table 1). They were taken from the JMMC Stellar Diameter Catalog (JSDC) version 2 (Chelli et al. 2016), with corresponding angular diameters of 6.16 ± 0.47 mas and 5.92 ± 0.56 mas, respectively. For each exposure, the visibility calibration was performed by dividing the raw squared visibility of the science target by the raw visibility of the calibrator corrected for its diameter. The calibrated total N-band spectrum of Betelgeuse was obtained on the first and third snapshots of the first night by multiplying the ratio between the target and calibrator raw fluxes measured by MATISSE at each wavelength, with a model of the absolute flux of the calibrator (Sirius).

The latter was taken from the PHOENIX grid (ACES-AGSS-COND; Husser et al. 2013). After that calibration step, we merged the four exposures to obtain a final calibrated OIFITS file for every snapshot.

After a first quality check, it appeared that some of the snapshots on Betelgeuse were corrupted and thus we excluded them from our modelling and analysis. For the first night, the visibilities of the four consecutive interferometric exposures of the second snapshot appeared very variable and significantly differed by more than 3σ. A fast-varying seeing, thus implying an unstable chopping and photometry measurement, is very likely the cause of it. For the second night, the first snapshot showed nearly zero visibilities and extremely noisy closure phases, both very different from the second snapshot that showed good non-zero visibilities and closure phases. The first snapshot likely suffered from problems of fringe coherencing during the observations and then residual optical path difference (OPD) correction at the data reduction stage. Finally, for the third night where the large AT array was used, it turned out that for all the baselines except for the shortest one (A0-G1), Betelgeuse was extremely resolved (V ≲ 10⁻³). That implied N-band correlated fluxes below the MATISSE N-band sensitivity limit with ATs (~5 Jy) and thus the data were unusable.

Figure 1 shows the calibrated squared visibilities and closure phases used in our modelling. The error bars contain two contributions: (1) a short-term one affecting the individual 1-min exposures, as computed by the MATISSE pipeline, and (2) the standard deviation between the merged exposures. We note that the calibration error associated with the stability of the MATISSE N-band transfer function over an observing night is estimated to be about 5% in average seeing conditions (Lopez et al. 2022), as was the case for our night (see Table 1). In our case, that was smaller than the two other error contributions. Figure 2 shows the (u, v) coverage of the small and medium configuration, and one baseline from the large configuration. The small and medium configurations probe scales of 60–300 mas and 20–80 mas, respectively. The longer baselines from the large configuration can probe even smaller scales; however, only the shortest baseline of this configuration is usable from these observations, which probes similar scales to the medium configuration. It is important to notice the sparseness of this coverage that prevents us from performing an image reconstruction.

In order to compare the flux retrieved from MATISSE to that of the star pre-dimming, we used a spectrum obtained through the Short Wave Spectrometer (SWS, de Graauw et al. 1996) instrument of the Infrared Space Observatory (ISO, Kessler et al. 1996). The version of the ISO data presented in this paper correspond to the Highly Processed Data Product (HPDP) set called ‘a uniform database of SWS 2.4–45.4 µm spectra’ by Sloan et al. (2003), available for public use in the ISO Data Archive².

Fig. 1 VLTI/MATISSE calibrated observations of Betelgeuse. Top: squared visibilities with 1-σ error bars shown in grey. Bottom: closure phases. The colour bar indicates the wavelength.

Fig. 2 (u, v) coverage of the small configuration shown in blue, the medium one in pink, and the large one in grey.

3 Data analysis

3.1 Squared visibilities and stellar diameter

From the van Cittert-Zernike theorem, complex visibilities are the Fourier transform of the objects’ brightness distribution. A value of |V| = 1 indicates the object is not resolved. An object is considered fully resolved when the first zero of its visibility function is observed, as is the case of the observations in this study. The squared visibilities for the three telescope configurations are shown in the top panel of Fig. 1 against spatial frequency, that is to say the length of the projected baseline divided by the wavelength. Going to larger spatial frequencies reveals finer details. Throughout this paper spatial frequency is expressed in arcsec⁻¹. While the wavelength range of the observations extends to 13 µm, we made a cut at 12.5 µm due to a high noise level beyond this wavelength. The dips in the visibilities of many of the baselines in the first lobe are particularly noticeable, with a relative minimum at approximately 10 µm. Given the wavelength at which these dips occur, we suggest that this could be a result of silicate dust around the star. This is explored further with dust radiative transfer modelling in Sect. 4.

We performed a per spectral channel uniform disk (UD) fit to the squared visibilities obtained from the two snapshots in the small configuration (which cover much of the first lobe), the results of which can be seen in Fig. 3. The fit is not extended towards longer baselines to avoid contamination by small-scale structures. The apparent stellar diameter appears to increase with wavelength, peaking between 10 and 11.5 µm. The presence of dust in the circumstellar environment is likely responsible for this. We expect that, in the N-band wavelength range, we should see the smallest contribution from dust for wavelengths <8.75 µm. Beyond this wavelength the contributions of dust and other components in the field of view become important; therefore, values for the disk diameter from λ > 8.75 µm should not be taken at face value as they are the result of a complex combination of the apparent angular diameter, resolved out flux, and the angular separation and flux ratio of components in the circumstellar environment. Previous interferometric observations of Betelgeuse where multi-component models have been applied to separate out the extended contribution at these longer wavelengths obtain disk diameters of ≈53–57 mas at ~11 µm (e.g. Danchi et al. 1994; Bester et al. 1996; Weiner et al. 2000; Perrin et al. 2007). It is worth noting that differences between angular diameter measurements could be caused by different spectral bandwidths, the evolution of the circumstellar environment, and, in the case of resolved dust emission, also the (u, v) coverage.

For the purposes of this study, we took the mean over the 8–8.75 µm interval to obtain a disk diameter, of 59.02 ± 0.64 mas. O’Gorman et al. (2017) retrieved a major axis diameter of 57.8 ±0.1 mas from a uniform elliptical disk fit to Atacama Large Millimeter Array (ALMA) band 7 (0.8–1.1 mm) observations. Kervella et al. (2018) also determined the diameter in this ALMA band to be 59 ±0.28 mas. VLTI/GRAVITY K-band (2.0–2.4 µm) data taken in the same month as our MATISSE observations give ${\theta }_{\star }^{\text{NIR}}=42.11\pm 0.05$ mas (Montargès et al. 2021), indicating that the star appears much larger in the N band than K band. This change in apparent diameter could be due to a close molecular layer of H₂O (Tsuji 2000; Ohnaka 2004; Montargès et al. 2014) and SiO (Decin et al. 2008; Kervella et al. 2018) around the star. This suggests that the diameter we measured here is a combination of the photo-spheric diameter and molecular layer closely surrounding the star (MOLsphere, Tsuji 2000). In reality, the opacity of these molecules varies over the N-band wavelength range, causing some variations in the apparent diameter, where SiO is most prominent under 10 µm and H₂0 is present across the entire range (see Perrin et al. 2007). While we find our N-band diameter result comparable to the disk diameter measured in the sub-millimetre, it is important to note that different opacity sources are dominant in each wavelength regime. For instance, H⁻ opacity is dominant for RSGs such as Betelgeuse at sub-millimetre wavelengths (Harper et al. 2001).

Since we do not have access to the angular diameter for wavelengths larger than 8.75 µm, except when using geometrical modelling assuming a certain type of geometry, we decided to use this constant diameter value, ${\theta }_{\star }^{\text{MIR}}$, for the photosphere. Figure 4 shows the first lobe part of the observations over-plotted with a 59.02 mas UD model; here, we can see that the model fits the general trend of the visibilities outside of the 10 µm feature, which we aim to reproduce using dust radiative transfer modelling in Sect. 4.

Fig. 3 Results of uniform disk fit as a function of wavelength, with 1-σ error bars shown in grey. The wavelength region 8–8.75 µm suffering the least from dust emission was used to fit a uniform disk diameter (see Fig. 4). The gap in the 9.2–9.5 µm region corresponds to flagged data.

Fig. 4 Squared visibilities from the small configuration are shown in blue with 1-σ error bars shown in light grey. A model of a uniform disk with a diameter of 59.02 mas is plotted in dark grey.

3.2 Closure phases

Closure phases are the phase sum over three baselines, which allows for the recovery of partial phase information that would otherwise be lost as a result of atmospheric turbulence. Deviations from 0 or +180° give information on the asymmetries of the system. The closure phases are shown in the bottom panel of Fig. 1. The small configuration, which probes scales of 60– 300 mas, shows little departure from 0° in the closure phases, indicating that we are not detecting significant asymmetries at these scales. However, the medium configuration, probing scales of approximately 20–80 mas, appears much more complex, indicating deviations from a centro-symmetric system at these scales. These deviations could be caused by clumpiness in the environment around the star and/or brightness variation across the surface of the star.

An attempt to reconstruct an image was made. However, due to the sparse (u, v) coverage of our observations (see Fig. 2), we were unable to converge on a definitive solution.

4 Radiative transfer modelling

In order to model the circumstellar dust around Betelgeuse, we used the 3D radiative transfer code RADMC3D (version 2.0, Dullemond et al. 2012). This code enabled us to input an arbitrary 3D dust density distribution, dust grain properties, and stellar parameters to create cubes of intensity maps, spanning the wavelength range of our observations. The interferometric observables, visibilities, and closure phases were then obtained by running these models through a VLTI simulator, ASPRO2 (Bourgès et al. 2013). From which point, we could directly compare our observations and models.

4.1 Parameters and assumptions

4.1.1 Stellar parameters

Levesque et al. (2005) measured a pre-dimming effective temperature, T_eff, of 3650 + 25 K. During the dimming, this temperature was measured at 3600 + 25 K (Levesque & Massey 2020), though a subsequent study by Začs & Pukītis (2021) suggests a larger temperature drop. As a drop in temperature during the dimming is expected to be localised (see e.g. Montargès et al. 2021), we opted for the higher pre-dimming temperature for our models. We note that a difference of 100 K has limited impact on the continuum SED at these mid-infrared (MIR) wavelengths. The SED of Betelgeuse was approximated with a PHOENIX stellar atmosphere model from Husser et al. (2013) with T_eff = 3700 K and surface gravity log g = 0 (Arentsen et al. 2019) scaled by the angular diameter in the K band (close to the H⁻ opacity minimum), ${\theta }_{\star }^{\text{NIR}}=42.11$ mas. However, these atmosphere models do not extend to our wavelength range and we therefore extrapolated using the Rayleigh–Jeans law. To test the validity of this atmosphere model and its extrapolation, we ran a comparison test with a MARCS stellar atmosphere model from Decin (2000) with the following parameters: T_eff = 3600 K, log g = 0.00, [Fe/H] = 0.00, and M = 15 M_⊙. We find that the visibilities produced using both models are virtually indistinguishable. This is unsurprising because the visibilities depend on the relative weights of the components and not the absolute values. We therefore opted to continue with the extrapolated PHOENIX model to remain consistent in comparing our findings to Montargès et al. (2021). The stellar diameter for the light emitting surface in the N band was set to ${\theta }_{\star }^{\text{MIR}}=59.02$ mas (see Sect. 3.1) at a distance of ${222}_{-34}^{+48}$ pc (Harper et al. 2017; Joyce et al. 2020). The MIR-disk diameter is thus ${13.1}_{-2.0}^{+2.9}$ au. It is currently not possible to include an inhomogeneous photosphere in the RADMC3D dust models. In Sect. 5, we explain how we tested the effect of a cool spot on the photosphere on our observables without the inclusion of dust.

4.1.2 Dust composition and grain size

Led by the findings of Verhoelst et al. (2009), we investigated three dust compositions: olivine (MgFeSiO₄), alumina (Al₂O₃), and the dust mixture found by these authors to best fit the SED of Betelgeuse, comprising of majorly melilite (Ca₂Al₂SiO₇) with smaller amounts of alumina and olivine (mass fractions of 0.64, 0.20, and 0.16, respectively). Optical constants were obtained from Dorschner et al. (1995), Begemann et al. (1997), and Mutschke et al. (1998) for olivine, alumina, and melilite, respectively. It was clear from early tests that alumina, as the sole composition of the radially outflowing wind, did not fit the shape of the observed squared visibilities; in particular, it does not reproduce the dip at 10 µm that can be seen most prominently at the two shortest baselines (see Fig. A.1). With this in mind, we continued on with the remaining two compositions for the rest of the study (see Sect. 4.3.1). We adopted a dust grain size range of a = 0.01−1 µm from Verhoelst et al. (2009) for our spherical grains, with a size distribution n(a) ∝ a^−3.5 (Mathis et al. 1977) as we deem this more physical than a single grain size in the ambient medium of RSGS where dust nucleation and growth are ongoing. Employing a size distribution also avoids non-physical resonances in radiative transfer modelling.

4.2 Dust distribution

Given the constraints of our data set, we did not aim to reconstruct the exact complex morphology of the inner wind of Betelgeuse. Instead we opted for three relatively simple morphological setups. The idea behind these setups is that they capture the basic components thought to be present in the ambient medium of RSGs. These are that a (roughly) spherical outflow is present through which mass is lost; however, part of the mass loss may be attributed to dusty clumps embedded in this outflow, suggesting episodic events of localised mass loss (e.g. Kervella et al. 2011; Ohnaka 2014).

With this in mind, the first setup consists of a radial outflow expanding at a constant velocity. Adopting a constant dust-togas ratio, the dust density distribution falls off as r⁻². We take the inner radius of our dust shell to be $13R_{\star }^{\text{NIR}}$ (see Sect. 3.1), as determined through SED fitting by Verhoelst et al. (2009) and a constant expansion velocity, v_exp, of 14 km s⁻¹ (De Beck et al. 2010), where we assume the gas and dust are well mixed and part of one flow, expanding at the same velocity. We then varied both the composition and the dust mass-loss rate (see Sect. 4.3.1). Kervella et al. (2011) presented VLT/VISIR MIR images of the circumstellar environment of Betelgeuse and identified a partial dust shell at an onset radius of $\approx 24R_{\star }^{\text{NIR}}$; however, these images also show the presence of dust closer in towards the star than said shell. For this reason, we chose this closer dust shell onset in order to approximate the dusty environment as a whole. For this and the following setups, we initialized a spherical grid (n_r, n_θ, n_ϕ) = (20, 20, 20), spanning radii from 8 to 400 au ($1.7-86R_{\star }^{\text{NIR}}$), and we employed four levels of adaptive mesh refinement to adequately resolve the dusty regions.

The second setup added a dusty clump positioned in front of the star in addition to the dusty envelope. The size, position, and density of the homogeneous and spherical clump are taken from Montargès et al. (2021), where they were determined through comparison of radiative transfer models to images from VLT/SPHERE-ZIMPOL. Our clump parameters are therefore as follows: r_c = 4.5 au, x_c = −1.9 au, y_c = −1.8 au, z_c = 20 au (where x and y are indicated in Fig. 5 and z is the line-of-sight position which is positive towards the observer), and the clump dust density is ρ_c = 2.0 × 10⁻¹⁸ g cm⁻³.

The third dust morphology setup we investigated removes the dust shell and places two more clumps in the close environment around the star, with identical properties to that described above. For this test, we used olivine as the dust composition for all clumps. The decision to add two clumps for this test scenario was motivated by previous SPHERE/ZIMPOL observations of the dust in the inner environment of the RSGs Betelgeuse and Antares, which both show a few patches of larger dust density within several stellar radii (Kervella et al. 2016; Cannon et al. 2021). We point out that our aim is not to try to constrain the properties and spatial configuration of an inner clumpy environment – apart from modelling challenges, we lack the (u, v) coverage to attempt this. Our aim is rather to probe the effects of this type of morphology on the visibility and closure phase signal. Schematics of all three RADMC3D models are shown in Fig. 5.

4.3 Modelling results

4.3.1 First setup: Composition and mass-loss rate

After a coarse exploration of parameter space to map the effect of the dust mass-loss rate on the visibilities, we zoomed in on the range from 1 × 10⁻¹⁰ to 1 × 10⁻⁹ M_⊙ yr⁻¹ in steps of 1 × 10⁻¹⁰ Μ_⊙ yr⁻¹for our two main chemical compositions. To assess how well our models recreate the observations, we performed a χ² estimation, simultaneously on both small configuration snapshots. We are able to provide an estimate of the confidence interval because we normalized the reduced χ² (${\chi }_v^2$) distribution so that the best model had ${\chi }_v^2=1$. From that point, we calculated the P value. All models with a P value higher than 0.05 (within the 95% confidence interval) are deemed acceptable models from which we take our parameter ranges. To account for our limited sampling rate, we performed a linear interpolation to the upper P-value points to give a better estimate on the confidence interval of the dust mass-loss rate. Considering the coarse sampling of our model grid, due to the exploratory nature of the dust geometries, using other statistical techniques is not possible. Since the main uncertainty lies in the model geometry and our (u, v) coverage, we do not deem it necessary to expand our model grid. The current estimates are meant to give an order of magnitude of the dust mass-loss rate from the MATISSE observations.

Here, we find that both an outflow composed of just olivine and one with the mixture of melilite, alumina, and olivine give us similarly good matches (see Fig. 6). We find dust mass-loss rates of $3_{-0.9}^{+0.9}\times {10}^{-10}{M}_{\odot }y{r}^{-1}$ and $4_{-0.9}^{+0.9}\times {10}^{-10}{M}_{\odot }y{r}^{-1}$ for olivine and the melilite mix, respectively. However, it is worth noting that a change of R_in or v_exp, which were kept constant in this study, would have an impact on the derived dust mass-loss rate. Since the model input here assumes a spherically symmetric outflow, we do not expect to see the closure phases depart from 0 or ±180°, which is in line with the model shown in Fig. 5. The few stray points in the closure phases are a result of the minor asymmetries caused by the pixel size (3.6 mas) of our model images inputted to ASPRO2.

4.3.2 Second setup: Impact of the dusty clump

The addition or omission of the dusty clump in front of the star, described in Sect. 4.2, appears to have little effect on the visibility shape giving almost the exact same χ², value either way. This is also evident when comparing the SEDs of our models: the addition of the dusty clump in the line of sight has a negligible impact (see Fig. 7). We return to this in Sect. 6. In line with these findings are the model images at 10 µm in Fig. 5 (second row) where both the images with and without a clump are virtually indistinguishable. Where we do clearly see the signature of the clump is in the closure phases of the medium configuration, as seen in the fourth row of Fig. 5. There, we can see departures from symmetry of our model at the same spatial scales as the observations; however, the observed closure phase (Fig. 1) signal is much more complex and unlikely the result of a singular dusty clump. We note that the temperature of the dust in this clump ranges from approximately 1200 K to 1800 K. Only a minor fraction (less than 10%) has a temperature in excess of 1500 K.

4.3.3 Third setup: Multi-clump model

The purpose behind the multi-clump model is not to reproduce the exact dust morphology in the inner environment of Betelgeuse, but to evaluate what effect multiple dust clumps would have on the visibilities and closure phases. A parameter study of possible clump configurations with the aim to match the closure phase signal is currently far beyond what is technically possible with our 3D models given the large parameter space, and the poor constraints provided by the (u, v) coverage of our data. For the multi-clump demonstration setup, we kept the dust clump in the line of sight of the star at a fixed position. We then placed two more identical dust clumps in positions within 20 au of the star (but in the plane of the sky, i.e. at z = 0). In order to find a tentative match for our visibilities, we then rotated the clumps around the star. As is apparent in the second row of Fig. 5, these additional clumps are much more visible at 10 µm than the clump in front of the star. While we did not aim to fit the visibilities with this simple test, we obtained a reasonable match to most baselines, shown in Figs. 6 and A.2, and Table. 2. This multiple clump model also creates more complex closure phases (see Fig. 5) which seem to better mimic the closure phase structure as seen in the observations.

Fig. 5 Comparison of the various models detailed in the main text with the MATISSE observables of Betelgeuse during the Great Dimming. First row: schematics of the different model setups (not to scale). Second row: intensity images (W m−2 µm−1 arcsec−2) output from RADMC3D at 10 µm for the three dust setups and one for the composite PHOENIX models representing the cool spot. The intensity is scaled with a power law with exponent 0.3 to best show the contribution of the different model elements. The image depicting the cool spot has been zoomed in on to better show the stellar disk (see legend for spatial scale used). Third row: squared visibilities from each of the models in small (blue), medium (pink), and large (grey) configurations. Bottom row: Closure phase for each of the models.

5 Cool spot model

In order to test if the MATISSE observations are consistent with a cool spot on the surface, we created a composite model following Montargès et al. (2021). We considered the area of the stellar photosphere to be π(59.02/2)² mas² and assigned an effective temperature of 3700 K (again scaling the stellar atmosphere models by ${\theta }_{\star }^{\text{NIR}}$). We then define a circular cool patch of radius 19 mas and centre (x, y) at (−2.4, −2.4) mas with respect to the stellar centre, on the surface of the star with a temperature of 3400 K. The composite image is shown in row 2 of Fig. 5, where it can be seen that the cool patch on the surface does not cause as high of a brightness contrast at 10 µm as it does in the visible: at 10 µm the specific intensity is 9.4% less in the cool spot than elsewhere on the stellar surface.

As predicted, the model visibilities now lack the 10 µm feature seen in the observations. However, the addition of a dusty envelope as described in Sect. 4.3.1, or other shapes of circumstellar dust, to the cool spot model would remedy this. The cool patch is clearly seen in the closure phases, showing a significant departure from symmetry; for more details, readers can refer to Fig. 5.

Fig. 6 Comparison of visibilities from RADMC3D models to the small configuration observations (2020 February 08(C)) split up by the baseline. The same models compared to 2020 February 08(A) are shown in Fig. A.2 and the residuals are plotted in Figs. A.3 and A.4. BL and PA denote the baseline length and position angle, respectively. For ${\chi }_v^2$ values, readers can refer to Table 2.

Fig. 7 Observed SEDs of Betelgeuse compared with the RADMC3D modelling. Left: spectral energy distribution from ISO-SWS01 taken pre-dimming (1997 October 08, see Sect. 2 for details) shown in blue and MATISSE shown in black (error bars shown in grey). Right: comparison of RADMC3D dust shell models, one with a clump positioned in front of the star and one without.

6 Discussion

6.1 Mass-loss rate and gas-to-dust ratio

From our spherical wind models, we obtained dust mass-loss rates in the range of (2.1–4.9) × 10⁻¹⁰ Μ_⊙ yr⁻¹ which agree well with those determined by Verhoelst et al. (2009). De Beck et al. (2010) derived a total mass-loss rate for Betelgeuse of 2.1 ×10⁻⁷ Μ_⊙ yr⁻¹from CO line profiles; however, Huggins et al. (1994) showed that CO is under-associated in the circumstellar envelope of Betelgeuse, suggesting this rate could be a lower limit. Hence, combining this mass-loss rate with our models suggests a minimum gas-to-dust ratio between 430 and 1000 for a homogeneous wind. A similarly high gas-to-dust ratio for Betelgeuse (≥550) and for the RSG Antares (≥600) was found by Justtanont et al. (1999). This suggests the dust formation in the outflow is not particularly efficient. It is possibly only efficient in clumps.

6.2 Detecting dust that potentially formed during the Great Dimming

The addition of the singular clump in front of the stellar disk with similar properties as derived in Montargès et al. (2021) does not significantly change the visibility signal. The dust clump also has a negligible contribution to the SED (right panel, Fig. 7) at these relatively long wavelengths in line with SED observations prior and during the dimming event (left panel, Fig. 7). Further investigation shows that the optical depth of the dust clump is such that the extinction (including scattering and absorption) caused by the clump almost exactly compensates for the dust emission; for more details, readers can refer to Figs. 8 and 9. This is not the case, however, if we use a different set of clump parameters, such as a clump with a higher density, for example. In which case, the near infrared (NIR) flux in the SED is reduced by the clumps’ presence in the line of sight. We note that if this clump was not in the line of sight, such as the other two clumps in our multi-clump setup, the IR emission does impact the SED.

Dharmawardena et al. (2020) used a spherical shell around the star to model the dust that would have caused the dimming, showing that the extra dust emission should have been detected in the SED at these wavelengths, using this as part of the reason to exclude dust as a possible cause of the Great Dimming. Harper et al. (2020) used a similar reasoning – a lack of an increase in emission compared to previous infrared photometry – as an argument against the dust hypothesis. Our modelling, however, shows that the localised presence of the newly formed dust in the line of sight of the star would go undetected at these wavelengths using these observational methods.

Alexeeva et al. (2021) conclude from TiO line fitting to spec-troscopic observations in the NIR at wavelengths shorter than 1 µm that only a cool photospheric patch could reproduce their observations. Montargès et al. (2021), however, found that a cool spot model produces a reduction of the flux at 1.6 µm that is significantly larger than what is observed during the dimming event. Further, Taniguchi et al. (2022) detected an enhancement in the optical depth at 10 µm during the dimming event, indicating the new formation of dust in the line of sight of the star.

Fig. 8 Optical depth through the centre of the dust clump, composed of olivine, as a function of wavelength as computed by RADMC3D. Beware that the clump is not covering the entire star.

Fig. 9 SED (blue dashed line, computed from the RADMC3D images outputted at different wavelengths) of the star with no dust. The orange line shows the SED of the star after each image was multiplied by the transmittance map, T = e−τ where τ is the optical depth, caused by the dust clump in the line of sight. Adding the dust clump emission to these new images results in the green line, which almost overlaps with the blue dashed line.

6.3 A potential clump, or potential clumps

While we did not fit the closure phases, it is clear from both the dust modelling (Sect. 4) and cool spot modelling (Sect. 5) that dust clumps (patches with local density enhancements) and surface features are the keys to reproduce the complex features seen in the MATISSE observations. This is unsurprising given the existing observations and theory of RSGs’ surfaces and winds that point to variations on spatial scales comparable to or smaller than the stellar diameter (see Sect. 1). While the addition of the single clump of dust, described in Sect. 4.2, in the line of sight of the star does cause some signal in the closure phase, it cannot fully explain the complexity of the observed closure phases. The situation is somewhat improved for our cool spot model. Therefore, we cannot, at present, conclusively rule out either of these scenarios, or indeed the presence of both simultaneously as concluded by Montargès et al. (2021).

Our pilot ‘multiple clump’ scenario shows that one may reproduce the visibilities and produce more intricate closure phase signals using clumps of dust in the inner circumstellar environment. While this model aims to show a proof of concept, ample evidence supports the presence of clumps; therefore, it seems more physically realistic than our other configurations, including the outflow models starting at $13R_{\star }^{\text{NIR}}$. The presence of clumps is supported by previous spatially resolved observations. One such observation by Kervella et al. (2016) using VLT/SPHERE/ZIMPOL shows a clumpy polarisation signal within 3 R_★, which indicates the presence of a patchy dusty environment. The patchy nature of the inner environment appears to show similarities to further out regions in the wind, with observations with VLT/VISIR by Kervella et al. (2011) identifying large dust clumps in the outer wind. Betelgeuse is not the only RSG to show these clumpy features. Also for Antares, whose clumpy wind has also been observed with VLT/SPHERE/ZIMPOL, VLT/VISIR, and MIRLIN at KECK II (Cannon et al. 2021; Ohnaka 2014; and Marsh et al. 2001, respectively) these inhomogeneities have been detected. NOEMA observations of µ Cep show that mass lost through clumps accounts for ≥25% of the RSG stars total mass loss (Montargès et al. 2019). Humphreys et al. (2021) also find dust clumps around the RSG VY CMa. Using data from the Hubble Space Telescope, they calculated the outward motion of these clumps allowing their ejection time to be calculated. The authors find that the ejection times of some of these clumps correspond with minima in the light curve of the star.

7 Conclusions

We obtained MATISSE N-band observations in February 2020 during the Great Dimming of Betelgeuse such as to probe the thermal emission of dust in the immediate surroundings of the star. Ideally, such observations may help to distinguish between the hypotheses proposed to explain this extraordinary decrease in visual brightness. To this end, we modelled the visibilities observed in the small VLTI configuration and examined the closure phases in the small and medium configurations.

From our parametric modelling of the VLTI/MATISSE data, we determined the uniform disk diameter of the star between 8 and 8.75 µm to be 59.02 ± 0.64 mas or a radius ${1409}_{-229}^{+319}$ R_⊙ at the adopted distance of ${222}_{-34}^{+48}$ pc. With the stellar size constrained, our three adopted dust models match the visibility data well: a spherical wind, a spherical wind with a dust clump in the line of sight, and three clumps placed around the star of which one is crossing the line of sight. This implies that these visibility data are not sufficiently sensitive to the spatial distribution and composition of dust in the field of view of MATISSE to accurately map the dust in the star’s environment. Better (u, v) coverage is needed to obtain a more complete image of the dust spatial morphology. Fitting a spherical homogeneous wind to the visibilities gives dust mass-loss rates of (2.1–4.9) × 10⁻¹⁰ Μ_⊙ yr⁻¹, suggesting a minimum gas-to-dust ratio between 430 and 1000 compared to the gas mass-loss estimate by De Beck et al. (2010), contrary to much lower values for outflows of asymptotic giant branch stars (where the canonical value is 100–200). This high gas-to-dust ratio indicates that dust formation may not be efficient in the wind or perhaps only efficient in clumps. Our models also exclude the possibility that the wind is dominated by Al₂O₃. However, this does not rule out that Al₂O₃ could be dominant in dust clumps in the inner wind (see e.g. Perrin et al. 2007; Cannon et al. 2021). The complexity of the closure phases from MATISSE suggests major asymmetries in the field of view (~1 arcsec); these could be caused by an asymmetric stellar disk, large-scale surface features, or dust clumps in the wind, or possibly a combination of all three. While our current data set, in principle, allows for a homogeneous wind (to reproduce visibilities) with stellar surface features (to account for the non-zero closure phases), there is observational evidence from VLT/SPHERE (Kervella et al. 2016) and VLT/VISIR (Kervella et al. 2011) that clearly shows an inhomogeneous dust distribution around the star.

In terms of the Great Dimming of Betelgeuse, we find that both models, cool spot and dust clump, or a combination of these models, are compatible with the observations. In particular, we find no inconsistencies with the modelling results presented by Montargès et al. (2021). We note that the dust clump from Montargès et al. (2021), positioned in the line of sight of the star, would be undetectable in the SED and visibilities in the mid-infrared as the extinction caused by the dust is directly compensated for by the dust emission. From our models, we can see that the presence of such a clump would only be detectable in the closure phases from the medium configuration (which probes scales of 20–80 mas). This clearly shows that in order to understand the nature of the Great Dimming, high angular resolution is mandatory to distinguish the photosphere from the circumstellar environment.

Future VLTI/MATISSE observations of RSGs and cool evolved stars, with good (u, v) coverage, will allow us to observe this close region where dust nucleation takes place in more detail. Mapping the dust across multiple epochs combined with ALMA observations of the gas morphology and kinematics would provide further insight into the connection of the geometry in the dust-forming region and the mass-loss properties and mechanism(s) of RSGs.

Acknowledgements

The authors acknowledge funding from the KU Leuven C1 grant MAESTRO C16/17/007. J.S.B. acknowledges the support received from the UNAM PAPIIT project IA 101220 and from the CONACYT project CF-2019-263975. M.M. acknowledges funding form the Programme Paris Region fellowship supported by the Région Ile-de-France. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant agreement no. 945298. 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 made use of the Jean-Marie Mariotti Center Aspro service (http://www.jmmc.fr/aspro). This research has also made use of the Jean-Marie Mariotti Center JSDC catalogue (http://www.jmmc.fr/catalogue_jsdc.htm).

Appendix A Additional comparison plots of RADMC3D models and MATISSE visibilities

Fig. A.1 Comparison of visibilities from RADMC3D of a radial outflow represented by a spherical dust shell composed of alumina to the small configuration observations (2020 February 08(C)) split up by the baseline. The baseline length and position angle are denoted by BL and PA, respectively.

Fig. A.2 Same as Fig. 6, but for snapshot 2020 February 08(A).

Fig. A.3 Residual of the squared visibilities of the small configuration (2020 February 08(C)) with respect to the RADMC3D models.

Fig. A.4 Same as Fig. A.3, but for snapshot 2020 February 08(A).

References

All Tables

All Figures

	Fig. 1 VLTI/MATISSE calibrated observations of Betelgeuse. Top: squared visibilities with 1-σ error bars shown in grey. Bottom: closure phases. The colour bar indicates the wavelength.
In the text

	Fig. 2 (u, v) coverage of the small configuration shown in blue, the medium one in pink, and the large one in grey.
In the text

	Fig. 3 Results of uniform disk fit as a function of wavelength, with 1-σ error bars shown in grey. The wavelength region 8–8.75 µm suffering the least from dust emission was used to fit a uniform disk diameter (see Fig. 4). The gap in the 9.2–9.5 µm region corresponds to flagged data.
In the text

	Fig. 4 Squared visibilities from the small configuration are shown in blue with 1-σ error bars shown in light grey. A model of a uniform disk with a diameter of 59.02 mas is plotted in dark grey.
In the text

Fig. 5 Comparison of the various models detailed in the main text with the MATISSE observables of Betelgeuse during the Great Dimming. First row: schematics of the different model setups (not to scale). Second row: intensity images (W m−2 µm−1 arcsec−2) output from RADMC3D at 10 µm for the three dust setups and one for the composite PHOENIX models representing the cool spot. The intensity is scaled with a power law with exponent 0.3 to best show the contribution of the different model elements. The image depicting the cool spot has been zoomed in on to better show the stellar disk (see legend for spatial scale used). Third row: squared visibilities from each of the models in small (blue), medium (pink), and large (grey) configurations. Bottom row: Closure phase for each of the models.

In the text

	Fig. 6 Comparison of visibilities from RADMC3D models to the small configuration observations (2020 February 08(C)) split up by the baseline. The same models compared to 2020 February 08(A) are shown in Fig. A.2 and the residuals are plotted in Figs. A.3 and A.4. BL and PA denote the baseline length and position angle, respectively. For ${\chi }_v^2$ values, readers can refer to Table 2.
In the text

	Fig. 7 Observed SEDs of Betelgeuse compared with the RADMC3D modelling. Left: spectral energy distribution from ISO-SWS01 taken pre-dimming (1997 October 08, see Sect. 2 for details) shown in blue and MATISSE shown in black (error bars shown in grey). Right: comparison of RADMC3D dust shell models, one with a clump positioned in front of the star and one without.
In the text

	Fig. 8 Optical depth through the centre of the dust clump, composed of olivine, as a function of wavelength as computed by RADMC3D. Beware that the clump is not covering the entire star.
In the text

Fig. 9 SED (blue dashed line, computed from the RADMC3D images outputted at different wavelengths) of the star with no dust. The orange line shows the SED of the star after each image was multiplied by the transmittance map, T = e−τ where τ is the optical depth, caused by the dust clump in the line of sight. Adding the dust clump emission to these new images results in the green line, which almost overlaps with the blue dashed line.

In the text

	Fig. A.1 Comparison of visibilities from RADMC3D of a radial outflow represented by a spherical dust shell composed of alumina to the small configuration observations (2020 February 08(C)) split up by the baseline. The baseline length and position angle are denoted by BL and PA, respectively.
In the text

	Fig. A.2 Same as Fig. 6, but for snapshot 2020 February 08(A).
In the text

	Fig. A.3 Residual of the squared visibilities of the small configuration (2020 February 08(C)) with respect to the RADMC3D models.
In the text

	Fig. A.4 Same as Fig. A.3, but for snapshot 2020 February 08(A).
In the text