Locating dust and molecules in the inner circumstellar environment of R Sculptoris with MATISSE

Context. Asymptotic giant branch (AGB) stars are one of the main sources of dust production in the Galaxy. However, it is not yet clear what this process looks like and where the dust happens to be condensing in the circumstellar environment. Aims. By characterizing the location of the dust and the molecules in the close environment of an AGB star, we aim to achieve a better understanding the history of the dust formation process. Methods. We observed the carbon star R Scl with the thermal-infrared VLTI/MATISSE instrument in L - and N -bands. The high angular resolution of the VLTI observations (as small as 4.4mas in the L -band and 15mas in the N -band with ATs), combined with a large u v -plane coverage allowed us to use image reconstruction methods. To constrain the dust and molecules’ location, we used two different methods: one using MIRA image reconstruction algorithm and the second using the 1D code RHAPSODY . Results. We found evidence of C 2 H 2 and HCN molecules between 1 and 3.4 R (cid:63) which is much closer to the star than the location of the dust (between 3.8 and 17.0 R (cid:63) ). We also estimated a mass-loss rate of 1 . 2 ± 0 . 4 × 10 − 6 M (cid:12) yr − 1 . In the meantime, we conﬁrmed the previously published characteristics of a thin dust shell, composed of amorphous carbon (amC) and silicon carbide (SiC). However, no clear SiC feature has been detected in the MATISSE visibilities. This might be caused by molecular absorption that can affect the shape of the SiC band at 11.3 µ m. Conclusions. The appearance of the molecular shells is in good agreement with predictions from dynamical atmosphere models. For the ﬁrst time, we co-located dust and molecules in the environment of an AGB star. We conﬁrm that the molecules are located closer to the star than the dust. The MIRA images unveil the presence of a clumpy environment in the fuzzy emission region beyond 4.0 R (cid:63) . Furthermore, with the available dynamic range and angular resolution, we did not detect the presence of a binary companion. To solve this problem, additional observations combining MATISSE and SAM-VISIR instrument should enable this detection in future studies.


Introduction
R Scl is a bright infrared source known to be a semi-regular pulsating C-rich asymptotic giant branch (AGB) star (Srb) with a pulsation period of about 370 days (Samus et al. 2009) and a mass-loss rate (published in the literature) ranging from 2 × 10 −10 M ⊙ yr −1  up to 1.6 × 10 −6 M ⊙ yr −1 (De Beck et al. 2010). Carbon stars are spectrally classified in R or N-type stars supplementing the classical hotter O, B, A, F, G, K, and M spectral types. While the N-type stars are colder than R stars, they are marked with a number identifying the effective temperature range of the star's location. Barnbaum et al. (1996) classified R Scl as a C-N5-type star in the revised MK classification of Keenan (1993), with C/O≈ 1.4 (Hron et al. 1998). Its distance is estimated to be 360±50 pc by Maercker et al. (2018), and measured as 440±30 pc (DR2) and 393 ± 12 pc (DR3) by Gaia Collaboration (2018). These two Gaia distances are marginally compatible within 1.5 σ with the values from Maercker et al. (2018) and with each other; therefore, we retained the conservative value of Maercker et al. (2018) (360±50 pc) for the purposes of this work.
R Scl is one of just a few detached shell objects (Olofsson et al. 1996) and ALMA observations revealed a spiral structure and the presence of clumpy and non-centrosymmetric structures inside the shell. Maercker et al. (2012) suggested it could be explained by the presence of a so far unseen companion of 0.25 M ⊙ hidden at 60 au ( 34 R ⋆ ) in the dusty shell. VLTI/PIONIER images do not show that companion, but instead show an extended photosphere with a dominant spot that is very likely to be of a convective origin (Wittkowski et al. 2017). Polarimetric observations confirm the clumpy structure of the dust forming region (Yudin & Evans 2002). However, such clumps were not detected in the MIDI observations of Sacuto et al. (2011). Such a non-detection was very likely due to the limitations in terms of uv-coverage and visibility sampling of MIDI, which is not optimized for the detection of large-scale asymmetric structures (Paladini et al. 2012(Paladini et al. , 2017. This paper is organized as follows. In Section 2, we present new VLTI/MATISSE observations in L-and N-bands of the AGB star R Scl, as well as the data reduction processes. In Section 3, we present the photometric data of this star and we derived a hydrostatic model to describe the central source. The resulting best-fit model is used as input for the DUSTY modeling used to interpret both the photometry and the interferometric MATISSE data. In Sections 4 and 5, we present the very first chromatic images in the L (3-4 µm) and N (8-13 µm) midinfrared (mid-IR) bands of the star's photosphere and its immediate vicinity. In addition, to locate the molecules and the dust in the circumstellar shell, we modeled the measured visibility data with a combination of ring-shaped layers within the close environment of the star and discuss our results. Finally, in Section 6, we provide our conclusions.

Observations and data processing
The Multi Aperture Interferometric SpectroScopic Experiment (MATISSE) is a four-telescope interferometric beam combiner that disperses the coherent light from four telescopes of the Very Large Telescope Interferometer (VLTI) in the L (3-4 µm), M (4.5-5 µm), and N (8-13 µm) mid-infrared bands (Lopez et al. 2022). It can combine either the four 1.8 m telescopes (auxiliary telescopes, AT), or the four 8 m telescopes (unit telescopes, UT). MATISSE delivers coherent fluxes, closures phases, differential phases, and photometry which can be used with the coherent flux to derive the relevant visibilities.
The four ATs can be moved along rails on the VLTI site, allowing for the baselines (separation vector between two telescopes with coordinates (u,v)) to more densely populate the (u,v)-plane than when only 4 fixed telescopes are used. This makes MATISSE an imaging instrument that offers an exquisite spatial resolution: reaching, in the L-band, up to λ L = 3 µm, λ L /B = 4.4 mas, and in the N-band up to λ N = 10 µm, λ N /B = 15 mas for a baseline of B = 140 m. The reliability of the image reconstruction process, the signal-to-noise-ratio (S/N), and the dynamic range of the image depend on the completeness of the (u,v)-plane coverage.
The observations of R Scl were conducted using the ATs of the VLTI during the commissioning run of MATISSE in December 2018. Table A.1 presents the logbook of the observations with the used AT configurations. Since on December 5 and 6, the observing nights were of poor quality characterized by technical glitches, only a fraction of the data from these two nights are used here. Fig. 1: SED of R Scl (log-log scales). Orange spectra: MA-TISSE data-sets in L-and N-bands with their associated error bars. Blue spectra: the ISO/SWS spectrum. Red points: photometric data from the literature. The solid black line is the best COMARCS+DUSTY model (see Section 3). Spectral features are labeled with arrows. Insets zoom onto the MATISSE spectrum (linear scales). SiC is marked as boldface as it is solid state, unlike the other gaseous species.
The data are processed using the MATISSE data reduction software drsmat version 1.5.0 (Millour et al. 2016), which adopts a classical Fourier transforms scheme (Perrin 2003), with photometric calibration in the case of the L-band (V. Coudé du Foresto et al. 1997). Chopping correction (subtracting sky frames from target frames) and optical path difference (OPD) demodulation procedures are applied before summing up the spectral densities, bispectra, and interspectra, leading to estimates of the spectro-interferometric observables: spectral distributions in squared visibility, closure phase, and differential phase, in a similar way as was done with the AMBER instrument (Tatulli et al. 2007). This leads to four reduced OIFITS2 files per pointed tar- Fig. 2: Plot of the MATISSE squared calibrated visibilities (grey dots) at 3.10 µm around the C 2 H 2 +HCN absorption feature (panel a), 3.52 µm in the L-band pseudo-continuum (panel b), 11.30 µm around the SiC feature (panel c), and 12.00 µm in the N-band pseudo-continuum (panel d). We over-plotted the data with several models: Hankel profile (blue), DUSTY model (orange), and MIRA mean 1D radial profile (red) and its dispersion (solid red interval). get (one for each beam commutation), which contains as many observables as exposures on the target.
In the next step, an interferometric calibration procedure is applied to the data: the targets used as calibrators (CAL) are corrected from their expected observables given the predicted value of their angular diameters, and each science target (SCI) is corrected from this transfer function estimate, owing to a calibrated OIFITS2 file for each CAL-SCI pair 1 . It is important to notice that calibrators have been chosen with no IR excess (see Cruzalèbes et al. (2019) for a description of the IR excess-free criterion) to avoid bias on the infrared measurements and in the calibration process. The drsmat has been developed in C with the ESO-specific library CPL, and is interfaced with the ESO recipes environment esorex. A graphical user interface can be used in the form of python tools developed by the consortium 2 to handle the data and recipes, and display the results.

Matching hydrostatic photosphere and 1D dust envelope model with the SED
Following the approach of Sacuto et al. (2011), we performed an estimation of the spectral energy distribution (SED) of R Scl at the pulsation phase of the MATISSE observations. Indeed, the SED is formed by selected photometric data points listed in Table A.3 and by the ISO/SWS spectra convolved at the same spectral resolution as MATISSE (R ∼ 30) and chosen at a pulsation phase that is as close as possible to that of the MATISSE observations (φ ∼ 0.95, see Figure 1).

Stellar pulsation
To estimate the pulsation phase of the MATISSE data, we made a periodogram analysis of the AAVSO V-band photometry. For this purpose, we downloaded data from between JD 2411972.3 and 2458893.5, and we selected only the best-quality data (based on the AAVSO data flag system). The periodogram not shown in this article is produced using the Lomb-Scargle method, an efficient algorithm for detecting and characterizing periodicity in unevenly-sampled time-data (VanderPlas 2018). The primary peak was detected at P = 372 ± 5 days which is consistent with the period 376 days from Wittkowski et al. (2017). We also detected a secondary peak at 15000 days (≈ 41 years), which was previously unknown. We ensured this was not an alias of the shorter period by checking that the secondary peak remains after we removed the primary pulsation from the data. Furthermore, the 41-year pulsation period is consistent with long sec-ondary periods detected in several AGB stars (Wood et al. 1999;Soszyński et al. 2021). Henceforth, we refer to the "pulsation" as the first 372-day period and "variability" as the long-term variability of 41 years. The zero-phase time is set at 2456245.5 days (=2012-11-14). The MATISSE observations were carried out close to maximum luminosity.

Photometric data and dereddening
The V, R, and I magnitudes reported in Table A.3 are taken from the catalog of Johnson et al. (1966). The errors on the flux are the results of errors in the magnitude, the color index, and the extinction. The J, H, K, L band magnitudes at the same pulsation phase as MATISSE data are taken from Whitelock et al. (2006). The photometric data are corrected from reddening using the visual extinction equation of Bagnulo et al. (1998). To retrieve the extinction at each fiducial wavelength, we assumed a totalto-selective extinction ratio in the visible of R = 3.1 (Fitzpatrick 1999a). We notice that the dereddening effect induced by the interstellar medium on the observed photometric data points is small.

Stellar contribution: COMARCS model fitting
In this section, the stellar contribution will be described using hydrostatic models. The detailed fit of the spectral energy distribution which would require more advanced models such as dynamic ones (Höfner et al. 2005) is beyond the scope of the current paper. Our aim here is to derive the stellar parameters and get a hydrostatic model that will be used in the following section as central source for the radiative transfer equation solver DUSTY 3 in order to model the total emission of R Scl. We fit the photometric data with the grid of precomputed low-resolution synthetic spectra COMARCS (Aringer et al. 2016(Aringer et al. , 2019. The spectra are calculated starting from 1D hydrostatic model atmospheres that take into account atomic and molecular absorption directly computed using the opacity generation code from the Copenhagen Opacities for Model Atmospheres (COMA) program ). The fit is limited to the wavelength region where the dust contribution is negligible, namely, between the V and I-bands (after dereddening). In such a wavelength domain, we could expect strong absorption from the circumstellar envelope (CSE) which is not considered in the models. Such dust extinction could, in principle, impact the effective temperature estimation of the COMARCS model. However, the bias described above is not significant for our purpose, as previously discussed in Brunner et al. (2018). Only a subset of 2600 COMARCS spectra matching metallicity, effective temperature, and stellar mass was selected from the original grid of 11248 models available in the official COMARCS repository 4 (Aringer et al. 2016). After computing χ 2 on the spectra subset, we selected one of the models within 3σ of the minimum χ 2 . Using Equation 2 from Cruzalèbes et al. (2013) we estimated the Rosseland radius associated to the model chosen in the previous step. Our result R ⋆ = 5.3 ± 0.6 mas of the star is consistent with Cruzalèbes et al. (2013) stellar limb-darkened radius of R ⋆,C = 5.03 ± 0.03 mas, based on hydrostatic MARCS model, and it is also consistent with the Rosseland radius R ⋆,W = 4.45 ± 0.15 mas of Wittkowski et al. (2017), obtained using dynamic atmosphere models. Even using different models, the latter two estimations are consistent within 1.5 σ with our Rosseland radius, meaning that the differences are not statistically significant. Finally, we also tested the COMARCS model against a distance of 393 pc (Gaia Collaboration 2018), following the same steps than described previously, and we did not notice any significant change in the SED and the stellar contribution.

Dust contribution: DUSTY model fitting
To describe the circumstellar envelope of the star, we use the 1D radiative transfer code DUSTY (Ivezić & Elitzur 1997). We set: (1) the central source as our best-fitting COMARCS model; (2) the dust grain-size distribution following the classical MRN (Mathis Rumpl Nordsieck) distribution (Mathis et al. 1977). For the latter, we assumed a minimum and maximum grain size a min = 0.005 µm and a max = 0.25 µm, respectively, which are largely adopted for size distribution among stellar envelopes and interstellar dust; (3) the inner shell temperature was set to 1200 K, which aptly represents the sublimation-condensation temperature of dust in particular AmC (Gail & Sedlmayr 1999). Such a temperature is affected, in practice, by unknowns related to dust grain optical and physical properties. As a consequence, we refer to the DUSTY radius where the dust grains condensate as the "DUSTY inner radius;" (4) then, we arbitrarily set the outer radius of the envelope to 1000 R in consistently to Sacuto et al. (2011); (5) finally, we assumed the density distribution driven by the pressure on dust grains. Indeed, the Radiatively Driven Winds Analytic (RDWA) approximation (Elitzur & Ivezić 2001) was used for the dust density profile. The approach is appropriate for AGB stars (in most cases) and according to the DUSTY manual, it offers the advantage of a much shorter computing time. Then, fitting the model on both observed SED and the MA-TISSE visibilities, we finetuned the amount of silicone carbide (SiC) and amorphous carbon (amC), as well as the optical depth at 1 µm. Thus, we estimated the inner shell radius and the massloss rate. We assume that the main error sources in the DUSTY fitting process come from the stellar temperature and luminosity. We used the lower and higher temperatures and luminosities of the COMARCS grid of models based on the error bars computed in Sect. 3.3. We repeated the process described in Sect. 3.4 and deduced the confidence interval of the DUSTY parameters. The summary of the DUSTY parameters are given in Table 1. Table 1 summarizes the parameters of the composite COMARCS+DUSTY best-fit model. The fractional abundance of the relevant grains in the dust shell resulting from this fit is a mixture of amorphous carbon (AmC, 88 ± 11%) and silicon carbide (SiC, 12 ± 11%), with the optical depth at 1µm τ λ=1µm = 0.19 ± 0.05. Our results are consistent with the findings of Sacuto et al. (2011), who used a grid of DUSTY models and least-squares fitting minimization to find a mixture of 90 ± 10% AmC and 10 ± 10% of SiC and an optical depth at 1 µm of τ λ=1µm = 0.18 ± 0.05. Then, we also estimated a mass-loss rate ofṀ = (1.2 ± 0.4) × 10 −6 M ⊙ yr −1 , using a dust grain bulk density ρ s = 1.85 g.cm −3 (Rouleau & Martin 1991) and a gas-to-dust-ratio of r gd = 590 (Schöier et al. 2005). Our estimation of the mass-loss rate is also consistent with that of De Beck et al. (2010), namely, 1.6 × 10 −6 M ⊙ yr −1 . In Fig. 1, we display the selected photometry data sets and the composite model of the SED. The MATISSE and the ISO/SWS spectra clearly show signatures possibly associated to acetylene (C 2 H 2 ) and hydrogen cyanide (HCN) molecules at 3.1 µm as well as the solid state SiC at 11.3 µm (Yang et al. 2004). We suspect the vertical-axis offset of the MATISSE spectrum with respect to the ISO observations and COMARCS+DUSTY models is due to calibration issues in terms of the MATISSE absolute flux or a field-of-view effect (Paladini et al. 2017). We note that COMARCS + DUSTY models cannot reproduce properly the CO + C 3 absorption feature at 5.2 µm. Part of this discrepancy can be linked to a well known problem about the opacity data and equilibrium constants involved in the C 3 band chemical description (Aringer et al. 2019).

Results of the composite model fitting
Since DUSTY does not handle molecules, and knowing the coexistence of molecules and dust, we expected to see strong discrepancies between the DUSTY model and the MATISSE data in the L-band, where strong molecular features are located. Although dust emission usually dominates the N-band, molecular contributions are also expected there, namely: CS and possibly SiO at 8 µm; HCN+C 2 H 2 molecular bands between 11 and 16 µm (Chubb et al. 2020) which can affect the shape of the SiC band at 11.3 µm.
The resulting fit to the SED (black curve in Fig. 1) provides results that are in reasonable agreement with the above expectations: a satisfactory fit to the N-band visibilities from MATISSE is found (orange curves in Fig. 2 c and d), but the L-band visibilities are not well matched with this star+dust model alone (Fig. 2  a and b). We therefore need to refine this model, introducing additional ingredients for the L-band.

R Scl as seen by MATISSE
To locate the dusty and molecular regions, we used two independent methods: the first one based on the MIRA image reconstruction algorithm and the second one a new approach that is presented in this paper and based on an intensity radial profile reconstruction method using the Hankel transform. The Hankel method has the advantage of requiring less parameters than image reconstruction, leading to fewer solution non-uniqueness problems, a better convergence, and a better dynamic range, given that the source is centro-symmetric at approximately all wavelengths.

MIRA image reconstruction
In Appendix. A.1, we show the (u,v)-plane coverage obtained during the observing run. The coverage is very dense, allowing us to reconstruct high-fidelity monochromatic images of R Scl. For that purpose, we fed the calibrated visibilities and closure phases into the MIRA software of Thiébaut (2008). Additional scripts written by Millour et al. (2012) allowed us to slice the wavelengths, and improve the image quality by low-frequencies filling.
The images were reconstructed with 128 pixel size, using total variation regularization, a hyperparameter value of 10,000 (as suggested in Renard et al. 2011), a random start image for the first wavelength, and the median image of all previous wavelengths for the subsequent ones. In rows c and d of Fig. 3, we show a subset of the reconstructed images.
Since the reconstructed interferometric image has a finer resolution than the diffraction-limited resolution, we also show in Fig. A.3 the reconstructed images convolved by an interferometer clean beam with a Gaussian profile. Furthermore, our MIRA image reconstruction effectively reproduces the closure phase of MATISSE data at all wavelengths -even at 3.11 µm, as shown in Appendix A.2, where strong asymmetries are observed.

Hankel Profile reconstruction
The Hankel profile method uses a large number of concentric uniform narrow rings of same width as that given by the angular resolution of the instrument. A detailed description of the RHAPSODY tool is given in Appendix B. Although it is built on simple models, RHAPSODY is not a simple model-fitting tool; rather, it is a radial profile reconstruction tool based on a Bayesian approach.
In this work, we reconstructed monochromatic radial profiles from visibilities using this technique. Each radial profiles is generated varying several parameters, including the number of rings, N step , the number of iterations in the minimization routine, the initial radial profile, and the value of the hyperparameter µ that we varied between µ = 1 and µ = 10 10 , in increments of the power of ten. The L-curve (χ 2 as a function of µ) provides us a guideline to select the optimal value using both regularization methods: 1) the quadratic smoothness f smooth , which tends to minimize the intensity variations between consecutive radii; and 2) the total variation f var which tends to reconstruct uniform intensities with sharp edges. Then we obtain the best radial profile reconstructions (shown in Figs. 3 and 4) with µ = 10 4 and a median value of the minimum reduced χ 2 of 3.3, using the total variation regularization method.
Images of the reconstructed Hankel profiles at selected wavelengths (i.e., pseudo-continuum (3.52 µm) and in C 2 H 2 and/or HCN absorption bands (3.10 µm, > 12 µm), as well as the SiC Fig. 3: Reconstructed Hankel distribution from radial profiles (row a), the same with the identified features highlighted (row b), images reconstructed with MIRA (row c), and the same with features highlighted (row d). Each panel shows the reconstructed image at the wavelength shown at its left top corner, covering the L-band (3-4 µm) and the N-band (8-13 µm). In these panels, the blue dashed circle in the center shows the calculated Rosseland radius of the stellar photosphere, the green circles represent the inner and outer boundary of the distinct molecular shell (seen here only at 3.52 µm and only clearly identifiable in the Hankel reconstruction), the blue circle represents the inner boundary of the dust envelope predicted by DUSTY, while the white circle at each bottom left corner shows the theoretical angular resolution of the interferometer. The white arrow shows the north-west arc structure referenced in Sect. 4.3.3. signature (11.30 µm)) are shown in Fig. 3. The estimated Rosseland diameter is marked as a dashed circle and the typical angular resolution of the interferometer for the fiducial wavelength, computed as 0.5 λ/B (Monnier 2003), is marked as a white circle in the bottom-left corner.
In Fig. 4, we show the obtained spectro-radial maps of the Hankel profile in the L-and N-bands. Profile reconstruction artifacts, coming from the (u,v)-plane sampling voids, are clearly identifiable as their distance from the star depends on the wavelength across the two L-and N-bands (traced as dashed white lines in Fig. 4). Spatially well-constrained structures show no dependence of the distance on wavelength (d = Const[λ]). We refer to Appendix C for more details on the fidelity of the profiles. The two grey vertical bands from 3.05 to 3.30 µm and from 3.75 to 4.00 µm cover the spectral ranges where the centro-symmetric Hankel profile is not able to properly reproduce the asymmetric Fig. 4: Spectra of the intensity profiles in L-and N-bands. Upper panels correspond to the observed MATISSE spectrum normalized by the median and then divided by a black body spectra at T=2700 K to underline the observed emission and absorption features. Lower panels are the spectro-radial maps in L-and Nband obtained by plotting the best intensity Hankel profile normalized at one for each observed wavelength. The faint red structures highlighted with the inclined dashed-white lines are reconstruction artifacts (see Sect. 4.2 for details). The dashed and solid blue horizontal lines show the position of the Rosseland radius and the DUSTY inner radius respectively. The solid green horizontal lines delimit the extension of a hot distinct molecular layer (see Sect. 5.3). The grey vertical bands cover the spectral ranges where the centro-symmetric Hankel profile is not able to properly reproduce the asymmetric shape of R Scl revealed by non-zero closure phases.
shape of R Scl, revealed by non-zero closure phases and by the reconstructed images (see the following section). The Rosseland radius is represented as a dashed blue line, the DUSTY inner radius in solid blue line, and a new, so far unexpected emission discussed in Section 5.3, which is between 12 and 18 mas (∼ 2.3-3.4 R ⋆ ) within the two solid green lines.
As seen in Fig. 4, the extent of the inmost circumstellar shell is not identical at all wavelengths but this variation is not coming from any instrumental or data reduction issue. The shell is significantly larger between 2.9 µm and 3.3 µm, and also 3.6 µm and 4.0 µm, and except for those specific bandwidths, it steadily increases towards longer wavelengths all the way up to 12 µm. Furthermore, this is in good agreement with Paladini et al. (2009): from dynamical model, the pseudo continuum radius is larger than the Rosseland radius because of the atmospheric extension and the pseudo-continuous molecular opacity. Notably, even at 3.525 µm, where we are expecting to find only the pseudo-continuum (Paladini et al. 2009), the radial profile extent is slightly larger than the value of the Rosseland radius. Therefore, our observations are qualitatively supportive of the synthetic observations derived on the basis of dynamic models described in Paladini et al. (2009).
The two reconstruction processes (profile and image) are different and, therefore, they produce different results that must be interpreted adequately: the Hankel profiles are computed using the visibilities across all position angles altogether, reconstructing only the centro-symmetric information of the object, whereas the MIRA images are computed using the visibilities and the closure phases as a function of position angles and they bear the asymmetry information. As a consequence, it is not surprising to find asymmetric structures in the MIRA images while not seeing these asymmetries in the Hankel profiles. Therefore, we should not go too far in the comparison of images coming from both methods: what can be compared are the radial positions of features, but not their relative intensities (except for azimuthally integrated profiles of MIRA images compared to the Hankel profiles shown in Fig. 6 and described in the following subsection). Furthermore, since it is difficult to estimate the true angular resolution and the shape of the PSF, the pixel size used to reconstruct the images is not meant to be the exact angular resolution of the image for both methods.

Description of the R Scl environment in the mid-IR
In this section, we describe the observations, while the interpretations can be found in Sect. 5. In Fig. 5, we provide a sketch of the region of the circumstellar envelope close to the star R Scl based on the Hankel profiles and the MIRA images, which we describe in more detail in the following sections. Such a circumstellar envelope is mainly composed of: i) the central component from 0 to 5 mas from the center of the star; ii) an inmost layer from 5 to 10 mas (i.e., ∼ 1-2 R ⋆ ); iii) a distinct layer from 12 to 18 mas (i.e., ∼ 2.3-3.4 R ⋆ ); iv) a fuzzy emission in the stellar outskirts at 20 to 90 mas (i.e., ∼ 3.8-17.0 R ⋆ ).

Central component (0 -5 mas)
In the first region between 0 and 5 mas from the center of the star (Figs. 3 and 4), we see a nearly round source dominating the flux in the pseudo-continuum (3.52 µm and 10 µm). Figure 6 shows spectra extracted from the Hankel profiles and the MIRA images in different regions around the star. These separated spectra allow us to better understand the nature and the composition of the different layer around the star. In Fig. 6, we see that the spectrum of the central source is very similar to the total one: an absorption feature at 3.05 µm and two bumps in the N-band at 8.7 µm and around 11.1 µm. This central part is clearly asymmetric in the 3.1 µm C 2 H 2 /HCN absorption band in the MIRA images.

Inmost layer (5 -10 mas)
A layer beyond the estimated Rosseland radius of R Scl is seen as a faint environment between 5 and 7 mas (∼ 1.0-1.3 R ⋆ ), assuming that the border of the photosphere is at 5 mas from the center of the star, as also seen in the images in Fig. 3. The Nband images also exhibit a larger central object compared to the L-band central star, meaning that we see the contributions from the star and the inmost layer, but we cannot disentangle them due to the lower angular resolution at these longer wavelengths.
We also observe a similar structure extending up to 10 mas (1.8 R ⋆ ) and attached to the central component, visible in the MIRA images at 3.1 µm and 3.9 µm as an irregular extended layer attached to the central component (Fig. 3) and as a ring around the central object in the Hankel profiles. In the spectra of Fig. 6, this is seen as two emissions bumps at 3.05 µm and 3.9 µm.

Distinct layer (12 -18 mas)
A layer that is apparently distinct from the previous one (i.e., continuum at 3.52 µm) is detected between 12 and 18 mas, below the DUSTY inner radius of 24.5 mas. Such layer is seen in the Lband in the MIRA image reconstruction (Figs. 3 and A.4) as a fuzzy arc to the north-west of the star and a few blobs at the same distance from the photosphere. The RHAPSODY intensity profile reconstruction shows also such layer as a continuous ring. The nature of such layer is further discussed in Sect.5.3.

Fuzzy emission in the outer skirts (beyond 20 mas)
A fuzzy emission above 20 mas (marked as a solid blue line in Figs. 3 and 4) is detected with MIRA and RHAPSODY. This emission is present both in L-and N-band. In the RHAPSODY images, the dynamic range of this emission is so low that only flux artifacts are reconstructed (see Appendix C). In the MIRA images, the (u,v)-coverage does not allow us to meaningfully reconstruct the shape of the brightness distribution. As a consequence, the fuzzy emission appears as blobs or artifacts in the north-eastern and south-western directions.

Discussion
In this paper, we describe 2D images that are projections of the actual R Scl 3D environment. Consequently, the distances mentioned in the following sections refer by default to distances projected onto the plane of the sky.

Central part: Star veiled by molecules
The data show a structure between 0-5 mas which can be interpreted as the star's photosphere. At the wavelengths where molecular bands are expected, the photosphere is no longer visible due to the presence of molecular shells in front of it. In the Fig. 6: MATISSE spectrum (black line) together with spectra extracted from circumstellar regions located at 0-5 mas (pink), 5-10 mas (green), 12-18 mas (blue), and 20-90 mas (golden) for both the Hankel (solid lines) and MIRA reconstructed images (dashed lines).
L-band at 3.05 µm, we can see a spectral feature which, according to, for instance, Aringer et al. (2019), could be attributed to C 2 H 2 + HCN. The feature detected at 3.9 µm is associated to C 2 H 2 .
The bump near 11.3 µm in the N-band (see Figs. 1 and 6) resembles the SiC signature. Two hypothesis are able to explain such feature: i) The signature is not from dust but rather from a pseudo-continuum between adjacent molecular features; ii) The signature is produced by an optically thick layer of SiC dust located along the line of sight. In the case where the layer temperature is lower than the one of the background (i.e., the stellar surface), an optically thin layer scenario would produce an absorption (Rutten 2000); hence, this would not be compatible with the emission observed. The limited angular resolution and the dynamic range of the MATISSE N-band data do not allow us to correctly probe the SiC emission region, which provides key information in making the choice of one scenario over the other.

Inmost hot molecular layer (5-10 mas)
The inmost layer has clear spectral features (seen in the green curve of Fig. 6), appearing as two emission bumps at 3.05 µm and 3.9 µm, namely, the signatures of C 2 H 2 and HCN. The hypothesis of the existence of an absorption band between these bumps is unlikely because spectral features at 3.05 µm and 3.9 µm should therefore also be in absorption. In addition, the data angular resolution in the L-band is high enough to avoid any contamination from the stellar emission. Furthermore, the existence of these emission bands is consistent with the lack of a strong continuum emission coming from the back hemisphere of the envelope, in the ring between 5 and 10 mas. On the other hand, C 2 H 2 and HCN are already formed in the stellar atmosphere at a temperature between 2500 K and 3000 K (Eriksson et al. 1984); thus, we suggest that this layer is attached to the star and could be part of the extended stellar atmosphere.
The inmost hot molecular layer is also visible at all other wavelengths in the L-band as a faint environment between 5 and 7 mas (∼ 1.0-1.3 R ⋆ ), as seen in the images in Fig. 3. It may also be seen in the N-band, as the images show a larger central object with respect to the photosphere of the star.

Possibility of a distinct hot molecular layer (12-18 mas)
In L-band panels of Fig. 3, the MIRA images indicate the presence of irregular emissions between 12 and 18 mas while RHAPSODY images (as also seen in Fig. 4) shows clearly the presence of a layer in emission. In contrast, in the N-band, neither MIRA nor RHAPSODY show any trace of such layer, but there is still a noticeable enlargement with regard to the angular diameter of the photosphere. The non-detection of such layer might be induced by the limited dynamic range and N-band angular resolution of the data (see Appendix C). The composition of this layer can then be of two different natures: either formed by dust or filled with molecules.

Dust hypothesis
The dust hypothesis, based on the assumption the 12-18 mas layer (2.3-3.4 R ⋆ ) is formed by dust, offers three possibilities: i) Dust grains in the layer are only made of amorphous carbon. According to our estimation using Equation 1 in Ivezić & Elitzur (1996), such layer should radiate at a mean temperature of 1500 K. This temperature is in accordance with the condensation temperature of amC in dynamic models (Nowotny et al. 2010). Wittkowski et al. (2017) also confirmed the presence of amC in the line of sight of their R Scl PIONIER images; ii) Dust grains in the layer are made of SiC. If we assume, as the temperature for the SiC dust condensation, T e f f = 1350 K from Yasuda & Kozasa (2012), then such a scenario would be discarded because the estimated temperature of the layer in this work is higher than the latter. On the other hand, we should keep in mind that the condensation temperature of SiC is not well known. Yasuda & Kozasa (2012) used select data on the condensation of Si x C y onto SiC dust grains that are affected by large errors. Gobrecht et al. (2017) analyzed part of the chemical reactions used by Yasuda & Kozasa (2012) and found significant differences. Agúndez et al. (2020) gave a condensation temperature for SiC of ≈ 1300-1400 K, while Men'shchikov et al. (2001) provided a temperature of 2000 K. Given these uncertainties, the scenario cannot be fully discarded; iii) Dust grains in the layer are made of titanium carbide (TiC). TiC also condenses at those temperatures (Agúndez et al. 2020), however such a chemical species is not included in this paper since there is no clear spectral signature of TiC in the MATISSE spectral range (von Helden et al. 2000).

Molecular hypothesis
In the case of a molecular layer, it might be formed by C 2 H 2 + HCN molecules. Indeed, in the ISO/SWS spectra of Fig. 1, we observe a minimum around 10 µm in N-band, which could be the pseudo-continuum between C 2 H 2 + HCN featuresat 7.5 µm [6.5 µm-9.0 µm] and 13.0 µm [11.0 µm-16.0 µm] Chubb et al. 2020. Such a possible pseudo-continuum would be bordered at long wavelengths by the SiC band. Furthermore, we should also notice that weaker molecular bands as SiH 4 , NH 3 , or C 2 H 4 might also contribute to the shape of the SED at those wavelengths. On the other hand, similar layers were detected around AGB stars in previous works: HD 187796 (χ Cyg) was observed by Perrin et al. (2004) with a molecular layer detected between 1.8-1.9 R ⋆ , Zhao-Geisler et al. (2012) revealed a molecular layer made of C 2 H 2 and HCN at about 2 R ⋆ around the carbon rich star V Hya. Similar sizes can be found in oxygen-rich AGB stars (BK Vir: 1.2-4.5 R ⋆ , SW Vir: 1.2-3.0 R ⋆ in Hadjara et al. 2019).
The molecule detection at this distance from the star echoes the results of Wittkowski et al. (2017) which modeled AMBER and PIONIER data of R Scl using COMA showing evidence for molecular contributions by CO, CN, C 2 , C 2 H 2 , HCN, and C 3 close to the star. It is worthwhile to note in this respect that PAHs and very small grains seem to show broad emission plateaus in the N-band which would produce a similar spectral signature (Peeters et al. 2017). Such particles belong to the formation process of amC and would thus be present in the dust condensation region.
Following the reasoning presented in Sect. 5.3.1 and Sect. 5.3.2, it is very likely that the 12-18 mas layers are composed by C 2 H 2 , HCN, and (possibly) amC dust.

Comparison with PIONIER data
In Figure 7, we compare our VLTI/MATISSE images with those from VLTI/PIONIER from Wittkowski et al. (2017) 5 .
The PIONIER 1.76 µm observations give us information mainly on the stellar surface and its extended atmosphere, being outside of the 1.53 µm C 2 H 2 +HCN molecular absorption band. We need to keep in mind that C 2 bands are also present in the observing bandwidth of PIONIER and could affect the images of the continuum emission. The 1.76 µm VLTI/PIONER data were observed in 2014, corresponding to a pulsation phase of ϕ = 0.78, and the VLTI/MATISSE data were observed in 2018 at a pulsation phase of ϕ = 0.96. The two epochs at four years apart are expected to show an evolution: the long-term variability of the star had time to completely reshape the region of the envelope close to the inmost layers of the vicinity of the star. The comparison shows similarities: the size of the central source (0-5 mas) with PIONIER is akin to our central source at 3.51 µm (in the pseudo-continuum).
It also exhibits some differences: asymmetries are not located at the same place between the two epochs. Two hypotheses for this variation: a significant evolution of the stellar and circumstellar environment shape or wavelength-dependent opacities related to the presence of molecules and possibly dust in front of the star. Then, this comparison has underlined the fact that it is possible to observe stellar surface variability at a timescale of few years for this kind of object. This means that in the case of a high observation frequency with MATISSE or another imaging interferometer, we should be able to estimate the characteristic timescale and sizes of patterns caused by convection onto the stellar surface.

Conclusion
In this paper, for the first time, we present our processing of simultaneous L-band and N-band observational data of the AGB star R Scl obtained with the VLTI/MATISSE instrument. Thanks to the good (u,v)-plane coverage of the MATISSE observations, we derived a consistent MIRA image and Hankel intensity profile reconstruction. The use of both independent and complementary methods allows us to locate the dust and molecules around R Scl. In the images (and intensity profiles), we find the signature of a sporadic mass-loss event maybe responsible for the formation of a distinct layer. Based on the MIRA images, we also detected some asymmetric structures confirming the presence of clumps  in the envelope of R Scl. We did not find any clear signature of SiC forming regions around the star. We present the main steps of our work as follows: 1. Data obtained with the VLTI/MATISSE instrument during the seven-days of the long commissioning session from 2018-12-03 to 2018-12-15 was used in this work. 2. Using the COMARCS stellar atmospheric model and the dust radiative transfer code DUSTY, we confirm the previously published characteristics of a thin dust shell, mainly composed of amorphous carbon with a small amount of silicon carbide (10%), leading to the location of the inner boundary estimated by DUSTY at about 4.6 R ⋆ from the central star (i.e., 24.4 mas) and a mass-loss rate of 1.2 ± 0.4 × 10 −6 M ⊙ yr −1 . 3. In order to locate the dust and molecule regions, we compared two different methods: one based on the MIRA image reconstruction and an other one based on a brand-new approach which consists of reconstructing the Hankel intensity profile of the circumstellar environment using the code RHAPSODY 6 . This method allowed us to simultaneously reproduce the spectral energy distribution and the MATISSE visibilities. 4. These complementary and independent approaches allow us to get a better view of both the composition and the geometry of the close environment of R Scl. We locate the dust emission region between 20 and 90 mas (∼ 3.8-17.0 R ⋆ ), as well as the presence of molecular shell around 5-10 mas (∼ 1-2 R ⋆ ) and a distinct one around 12-18 mas (∼ 2.3-3.4 R ⋆ ). This distinct layer could be made of both: dust and molecules. To reconstruct the formation history of such layer would require a monitoring of the source. 6 Link to RHAPSODY: https://github.com/jdrevon/RHAPSODY Our results are in qualitative agreement with the predictions of dynamic models (Paladini et al. 2009) and a detailed comparison with such models is the subject of further study.
Acknowledgements. MATISSE has been built by a consortium composed of French (INSU-CNRS in Paris and OCA in Nice), German (MPIA, MPIfR and University of Kiel), Dutch (NOVA and University of Leiden), and Austrian (University of Vienna) institutes. It was defined, funded and built in close collaboration with ESO. The Conseil Départemental des Alpes-Maritimes in France, the Konkoly Observatory and Cologne University have also provided resources to manufacture the instrument. The results are based on public data released from the MATISSE commissioning observations at the VLT Interferometer under Programme IDs 60.A-9257(E). We have great remembrance of the contributions from our two deceased colleagues, Olivier Chesneau and Michel Dugué to the MATISSE instrument. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. This work was supported by the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES, as well as from the Programme National de Physique Stellaire (PNPS) from CNRS. It has also been supported by the French government through the UCA-JEDI Investments in the Future project managed by the National research Agency (ANR) with the reference number ANR-15-IDEX-01 and the Hungarian governent through NKFIH grant K132406. Vincent Hocdé is supported by the National Science Center, Poland, Sonata BIS project 2018/30/E/ST9/00598. The authors made extensive use of the Jean-Marie Mariotti Center (JMMC) tools to prepare this paper. We can cite among others ASPRO 7 used to prepare observations, and OIdb 8 to publish the reduced data. The authors are very grateful to the anonymous referee for his/her careful review of the manuscript, resulting in insightful recommendations and constructive suggestions leading to the great improvement of its content and presentation. We are also very grateful to Orlagh Creevey for her precious information regarding the Gaia distances of R Scl, and Bernhard Aringer for the fruitful discussions about the COMARCS models. Appendix A.3: (u,v)-plane coverage  Notes: φ is the pulsation phase and A λ is the extinction value given with relative errors of 25%, computed thanks to the tabulated extinction values E(B-V) from Fitzpatrick (1999b). Magnitude is the observed magnitude provided by the catalog and the last column Flux is the dereddened flux computed using the A λ values.