VLTI-MATISSE L- and N-band aperture-synthesis imaging of the unclassified B[e] star FS Canis Majoris

K. -H. Hofmann, A. Bensberg, D. Schertl, G. Weigelt, S. Wolf, A. Meilland, F.

Millour, L. B. F. M. Waters, S. Kraus, K. Ohnaka, et al.

► To cite this version:

K. -H. Hofmann, A. Bensberg, D. Schertl, G. Weigelt, S. Wolf, et al., VLTI-MATISSE L- and N-band aperture-synthesis imaging of the unclassified B[e] star FS Canis Majoris. Astronomy and Astrophysics - A&A, 2022, 658, pp.A81. 10.1051/0004-6361/202141601. insu-03656928

HAL Id: insu-03656928

https://insu.hal.science/insu-03656928

Submitted on 2 May 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# VLTI-MATISSE L- and N-band aperture-synthesis imaging of the unclassified B[e] star FS Canis Majoris*

K.-H. Hofmann¹, A. Bensberg², D. Schertl¹, G. Weigelt¹, S. Wolf², A. Meilland³, F. Millour³, L. B. F. M. Waters⁴,⁵,
S. Kraus⁶, K. Ohnaka⁷, B. Lopez³, R. G. Petrov³, S. Lagarde³, Ph. Berio³, F. Allouche³, S. Robbe-Dubois³,
W. Jaffe⁸, Th. Henning⁹, C. Paladini¹⁰, M. Schöller¹¹, A. Mérand¹¹, A. Glindemann¹¹, U. Beckmann¹, M. Heininger¹,
F. Bettonvil¹², G. Zins¹⁰, J. Woillez¹¹, P. Bristow¹¹, P. Stee³, F. Vakili³, R. van Boekel⁹, M. R. Hogerheijde⁸,¹³,
C. Dominik¹³, J.-C. Augereau¹⁴, A. Matter³, J. Hron¹⁵, E. Pantin¹⁶, Th. Rivinius¹⁰, W.-J. de Wit¹⁰, J. Varga⁸,¹⁷,
L. Klarmann⁹, K. Meisenheimer⁹, V. Gámez Rosas⁸, L. Burtscher⁸, J. Leftley³, J. W. Isbell⁹, G. Yoffe⁹,
E. Kokoulina³, W. C. Danchi³,¹⁸, P. Cruzalèbes³, A. Domiciano de Souza³, J. Drevon³, V. Hocdé³, A. Kreplin⁶,
L. Labadie¹⁹, C. Connot¹, E. Nußbaum¹, M. Lehmitz⁹, P. Antonelli³, U. Graser⁹, and C. Leinert⁹

*(Affiliations can be found after the references)*

Received 21 June 2021 / Accepted 14 November 2021

## ABSTRACT

**Context.** FS Canis Majoris (FS CMa, HD 45677) is an unclassified B[e] star surrounded by an inclined dust disk. The evolutionary stage of FS CMa is still debated. Perpendicular to the circumstellar disk, a bipolar outflow was detected. Infrared aperture-synthesis imaging provides us with a unique opportunity to study the disk structure.

**Aims.** Our aim is to study the intensity distribution of the disk of FS CMa in the mid-infrared *L* and *N* bands.

**Methods.** We performed aperture-synthesis imaging of FS CMa with the MATISSE instrument (Multi AperTure mid-Infrared SpectroScopic Experiment) in the low spectral resolution mode to obtain images in the *L* and *N* bands. We computed radiative transfer models that reproduce the *L*- and *N*-band intensity distributions of the resolved disks.

**Results.** We present *L*- and *N*-band aperture-synthesis images of FS CMa reconstructed in the wavelength bands of 3.4–3.8 and 8.6–9.0 µm. In the *L*-band image, the inner rim region of an inclined circumstellar disk and the central object can be seen with a spatial resolution of 2.7 milliarcsec (mas). An inner disk cavity with an angular diameter of ~6 × 12 mas is resolved. The *L*-band disk consists of a bright northwestern (NW) disk region and a much fainter southeastern (SE) region. The images suggest that we are looking at the bright inner wall of the NW disk rim, which is on the far side of the disk. In the *N* band, only the bright NW disk region is seen. In addition to deriving the inclination and the inner disk radius, fitting the reconstructed brightness distributions via radiative transfer modelling allows one to constrain the innermost disk structure, in particular the shape of the inner disk rim.

**Key words.** techniques: interferometric – techniques: image processing – circumstellar matter – stars: emission-line, Be – stars: imaging – stars: individual: FS CMa

## 1. Introduction

B[e] stars are emission line stars with low-excitation permitted emission lines (e.g., Fe II), forbidden lines of [Fe II] and [O I], and an infrared excess from circumstellar dust (see Zickgraf 1998 and Lamers et al. 1998). This B[e] star definition has the disadvantage that B[e] stars are not a homogeneous group of stars, but it consists of different types of stars of different evolutionary stages. FS CMa (HD 45677, MWC 142) belongs to the subgroup of unclassified B[e] stars. Its distance is 620±41 pc (Vioque et al. 2020). The evolutionary stage of FS CMa is still debated. It is not known whether FS CMa is a young star, a star near the main sequence, or a more evolved star (Zickgraf 1998; Laws et al. 2020 and references). Vioque et al. (2020) conclude that FS CMa together with several other famous B[e] stars are most likely not pre-main-sequence stars (see list of these stars in Appendix A of their paper).

FS CMa is famous for its large, inclined circumstellar disk, which was studied with several infrared interferometers (Monnier et al. 2006; Menu et al. 2015; Lazareff et al. 2017; Gravity Collaboration 2019). Kluska et al. (2014, 2020) reported the first *H*-band aperture-synthesis images of the inclined, asymmetric disk of FS CMa. Laws et al. (2020) resolved a bipolar outflow perpendicular to the disk plane. Aperture-synthesis imaging with the new MATISSE interferometry instrument (Multi AperTure mid-Infrared SpectroScopic Experiment) at the Very Large Telescope Interferometer (VLTI) of the European Southern Observatory provides us with the unique opportunity to image the disk of FS CMa with unprecedented spatial resolution in the *L* and *N* bands.

In this paper, we present the first milliarcsecond-resolution mid-infrared *L*- and *N*-band aperture-synthesis imaging of FS CMa and radiative transfer modeling for the interpretation of the images. In Sect. 2, we describe the observations and data reduction. In Sect. 3, we present the mid-infrared *L*- and *N*-band aperture-synthesis images. In Sect. 4, we show our modelling of the disk structure of FS CMa based on the reconstructed high-resolution images. The results are summarized in Sect. 5.

<footnote>* Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere within the commissioning of the VLTI-MATISSE instrument (ID 60.A-9257(E)).</footnote>
## 2. Observations and data reduction

FS CMa was observed with the four Auxiliary Telescopes (ATs) of the VLTI (Schöller 2007) and the mid-infrared MATISSE instrument (Lopez et al. 2014, 2022; Matter et al. 2016; Petrov et al. 2018) during the MATISSE commissioning run between December 2 and 15, 2018, as summarized in Table A.1. The observations were carried out in the mid-infrared *L* and *N* bands with low spectral resolution of *R* ~ 30. The observations were performed with seven different AT quadruplet configurations with projected baseline lengths ranging between 10.5 and 137.0 m. Figure A.1 shows the *uv* plane coverage of all observations of the FS CMa observing run. The spectrally dispersed *L*- and *N*-band interferograms were recorded simultaneously. The detector integration time (DIT) was 75 ms in the *L* band and 20 ms in the *N* band. The calibrator stars are listed in Table A.2. The *L*-band interferograms were obtained in the science-photometry (SciPhot) mode, that is, interferograms and photometric data were recorded simultaneously. On the other hand, the *N*-band interferograms were obtained in the high-sensitivity (HighSense) mode, that is, interferograms and photometric data were recorded sequentially.

Each *L*- and *N*-band MATISSE observation usually consists of at least four interferogram data sets corresponding to four different configurations of the Beam Commuting Device (BCD, see Millour et al. 2008 and Lopez et al. 2014). These four data sets are recorded without chopping. Additionally, each observation consists of eight more data sets recorded with chopping. In the *N* band, these eight chopping data sets contain the photometric data from the four telescopes (HighSense mode) recorded with two different configurations of the BCDs. During the recording of these eight chopping data sets, additional *L*-band interferograms with photometric data (SciPhot mode) are recorded with the same two BCD configurations.

The data obtained during all nights of this run, except one night, were used for the subsequent data reduction and image reconstruction. From the four nights with the AT configuration A0-B2-J2-C1, the data of the night 2018-12-05 with the worst seeing conditions were not used because a reliable calibration was not possible (see Table A.1 for more details).

### 2.1. Data reduction of the *L*-band data

The *L*- and *N*-band data were reduced with the MATISSE Data Reduction Software (DRS) of pipeline version 1.5.3. After the data reduction, (a) four sets of interferometric data (closure phases and visibilities) are obtained from the data sets corresponding to all four BCD configurations without chopping, and (b) two sets of raw interferometric data from the eight chopping data sets recorded with two different BCD configurations. The transfer functions of the raw interferometric data of the target were calibrated using the calibrator stars listed in Table A.2. In Appendix B.1, we show examples of observed closure phases and visibilities. The estimation of the visibility and closure phase errors is also described in Appendix B.1.

### 2.2. Data reduction of the *N*-band data

In order to reduce the noise in the data, the *N*-band interferograms were spectrally binned, that is, neighboring spectral channels were averaged before calculating the average closure phases and visibilities. In Appendix B.2, we show calibrated FS CMa *N*-band closures phases and visibilities without spectral binning and with three different spectral binnings of 9, 11, and

17 spectral channels. Spectral binning of 7 (e.g., seven neighboring spectral channels in the spectrally dispersed interferograms are averaged) does not reduce the spectral resolution of the data because of the spectral oversampling of the dispersed interferograms recorded in the MATISSE low spectral resolution mode. Data processing without spectral binning (i.e., binning 1) is just shown for illustration in Appendix B.2, but is usually not used because of the severe noise. The transfer functions of the derived four interferometric data sets were calibrated with the calibrator stars listed in Table A.2. The estimation of the visibility and closure phase errors is described in Appendix B.2.

## 3. Aperture-synthesis image reconstruction of the dust disk of FS CMa

The MATISSE image reconstruction package IRBis (Hofmann et al. 2014) was applied to the calibrated interferometric data of FS CMa. This software package is part of the MATISSE data reduction pipeline<fnref n="1" />. For the reconstruction process, the two built-in minimization engines, the six regularization functions, and the three cost functions were applied (details are described in Hofmann et al. 2014 and Appendix D). The start image of each image reconstruction run was a circular Gaussian fit to the measured visibilities of the target. We used this Gaussian as a prior for the image reconstruction. From the reconstructions that were calculated, the best ones were selected using the image quality measure *qrec*. The quality measure *qrec* is calculated from the reduced $\chi^2$ values of the closure phases and visibilities, and a so-called residual ratio (see Hofmann et al. 2014 for more details). Because of the inhomogeneous *uv* plane coverage, we weighted the *uv* point measures with the power of 0.5 or zero of the inverse *uv* plane density function (see Table 1; power zero means equal weights for all *uv* points).

### 3.1. *L*-band images of FS CMa

The field-of-view (FOV) and pixel grid used for image reconstruction in the *L* band were 80×80 mas and 64×64 pixels, respectively. The image reconstruction was perfomed in the wavelength range 3.4–3.8 µm. All six calibrated interferometric data sets (four corresponding to the four BCD configurations without chopping, and two to the two BCD configurations with chopping) of each observation were used in the image reconstruction process, as described in Sect. 2.1. Figure 1 (top) shows the obtained *L*-band reconstruction. The reconstruction parameters and quality measures are listed in Table 1. The reconstructed image has a spatial resolution of ~2.7 mas (it was not convolved to reduce the resolution). In Fig. C.1, the measured closure phases and visibilities of an example observation and the ones derived from the reconstructed image (Fig. 1 top) are displayed. This comparison is a useful method to investigate whether the closure phases and visibilities of the observations are in acceptable agreement with the closure phases and visibilities of the reconstruction.

The *L*-band image shows a central object and a ring-like, elliptical circumstellar structure that looks like an inclined disk consisting of a bright NW disk region and a much fainter SE region. The first aperture-synthesis image of FS CMa was reported by Kluska et al. (2014, 2020), who imaged the disk of FS CMa in the *H* band. In this *H*-band image, the NW disk region and the central object can be seen, whereas the faint SE *L*-band region is not visible. The inner disk cavity in the *L*-band

<footnote n="1">¹ https://www.eso.org/sci/software/pipelines/</footnote>
<table><caption>**Table 1.** Image reconstruction parameters and quality of the reconstructed FS CMa images.</caption><thead><tr><th>Wavelength <i>r</i><sup>(a)</sup></th><th>Bin<sup>(b)</sup></th><th>CP<sup>(c)</sup></th><th><i>V</i><sup>2</sup><sup>(d)</sup></th><th>CP err.<sup>(e)</sup></th><th><i>V</i><sup>2</sup> snr<sup>(f)</sup></th><th>FOV<sup>(g)</sup></th><th>npix<sup>(h)</sup></th><th>Cost<sup>(i)</sup></th><th>Reg.<sup>(j)</sup></th><th>Opt.<sup>(k)</sup></th><th>Weight<sup>(l)</sup></th><th>Mask<sup>(m)</sup></th><th>CP res.<sup>(n)</sup></th><th><i>V</i><sup>2</sup>res<sup>(n)</sup></th></tr></thead><tbody><tr><td>3.4–3.8 µm</td><td>1</td><td>11 355</td><td>17 033</td><td>2.9°</td><td>18.4</td><td>80</td><td>64</td><td>3</td><td>6</td><td>2</td><td>0.5</td><td>30–42</td><td>2.9°</td><td>0.009</td></tr><tr><td>8.6–9.0 µm</td><td>11</td><td>2084</td><td>3079</td><td>10.0°</td><td>11.5</td><td>360</td><td>128</td><td>2</td><td>3</td><td>1</td><td>0.0</td><td>160–190</td><td>7.5°</td><td>0.021</td></tr></tbody></table>

**Notes.** (a)Wavelength range used for image reconstruction. (b)Number of spectral channels binned (spectral binning; see Sect. 2.2 for details). (c)Number of closure phases. (d)Number of visibilities. (e)Mean closure phase error of the data. (f)Mean signal-to-noise ratio of the visibilities. (g)FOV of the reconstruction in mas. (h)Pixel grid npix × npix used in the image reconstruction. (i)Number of the cost function as defined in Hofmann et al. (2014) and the Appendix D. (j)Regularization function number as defined in the same paper. (k)Optimization routine number as defined in the Appendix D. (l)Power of the uv density weight (see Hofmann et al. 2014). (m)Radii used for the image space mask in mas. (n)Closure phase (CP) and squared visibility ($V^2$) residuals derived from the measured closure phases and visibilities, and from the closure phases and visibilities calculated from the reconstructed images.

image has a size of ~6 × 12 mas. The asymmetric intensity distribution suggests that the NW region is bright because we are looking at the bright, inner wall of the NW disk rim. Therefore, this intensity distribution suggests that the NW rim is on the far side.

Reconstructed interferometric infrared images can consist of both real target structures and reconstruction noise caused by gaps in the uv plane coverage and calibration errors. To estimate this noise in our reconstructed images, we computed error maps of the reconstructed images presented in Figs. E.1 and E.2. For this goal, we first performed computer simulations, in which we used the same uv coverage as the one of the reported MATISSE observations. As computer simulation targets, we chose the model images shown in Figs. E.1a and E.2a (which are similar to the model images shown in Figs. 6 and 7, but with the central object intensity increased to the central object intensity of the reconstructed image). The computed visibilities and closure phases were degraded by simulated noise corresponding to the noise in the MATISSE observations. Figures E.1b and E.2b show the obtained reconstructed images of the computer simulations. The error maps shown in Figs. E.1c and E.2c are computed as the absolute value of the difference between the computer simulation targets (Figs. E.1a and E.2a) and the reconstructions of the computer simulations (Figs. E.1b and E.2b). The L-band error map presented in Fig. E.1c suggests that the clumpiness in the disk structure is partially degraded by the aforementioned reconstruction noise. This is consistent with the fact that there are differences between the positions of the brightest clumps in Kluska’s H-band image (Kluska et al. 2014, 2020) and the positions of the brightest clumps in the MATISSE image. Furthermore, the computer simulation clearly shows that the elongation of the reconstructed central object is at least partially a bridge-like reconstruction artifact caused by the neighborhood of the bright disk rim. This bridge effect is even partially seen in the original computer object in Fig. E.1a. However, a close suspected companion star (Baines et al. 2006) or a complex circumstellar environment may also partially contribute to the elongated image of the central object.

## 3.2. N-band images of FS CMa

For the N-band image reconstruction, the FOV and pixel grid used were 360 × 360 mas and 128 × 128 pixel, respectively. The image reconstruction was performed in the wavelength range of 8.6–9.0 µm, where the N-band data have their highest quality. All four calibrated interferometric data sets (corresponding to the four BCD configurations without chopping) of each observation were used for image reconstruction, as described in Sect. 2.2.

Figure 1 (bottom) shows the obtained N-band reconstruction. The reconstruction parameters and quality measures are listed in Table 1. The reconstruction has a spatial resolution of ~6.6 mas. Figure C.3 shows the comparison between the measured closure phases and visibilities of a representative measurement and the ones derived from the reconstruction (Fig. 1 bottom).

# 4. Radiative transfer modeling of the FS CMa images

The reconstructed images presented in Sect. 3 provide the basis to constrain the disk structure of FS CMa. In Sect. 4.1, the applied model and radiative transfer method are briefly outlined. In Sect. 4.2, the fitting process for the L-band image as well as the comparison with the N-band image are described. For this, we are using a well-studied basic disk model. Subsequently, we extend this model motivated by the remaining differences between the simulated and reconstructed images.

## 4.1. Model description

**Basic disk model 1.** We start with a basic model of a circumstellar disk that has been successfully used for fitting of observations of several protoplanetary disks (e.g., Wolf et al. 2003; Sauter et al. 2009; Brunngräber et al. 2016). The stellar parameters are fixed using values derived by Vioque et al. (2018). These and the disk inclination determined by Kluska et al. (2020) are summarized in Table 2.

In addition to the model referred to above, an exponential decay of the density in the outer regions based on considerations by Lynden-Bell & Pringle (1974), Kenyon & Hartmann (1995), and Hartmann et al. (1998) is applied. The density distribution is described by

$$ \rho_g(r, z) = \frac{\Sigma_g(r)}{\sqrt{2\pi} h_g(r)} \exp\left[-\frac{1}{2}\left(\frac{z}{h_g(r)}\right)^2\right]. \quad (1) $$

Here (r, z) denote cylindrical coordinates, $\Sigma_g$ the vertically integrated surface density, and $h_g$ the scale height, which can be written as

$$ \Sigma_g(r) = \sqrt{2\pi} \rho_0 h_{\text{ref}} \left(\frac{r}{R_0}\right)^{-\gamma} \exp\left[-\left(\frac{r}{R_{\text{trunc}}}\right)^{2-\gamma}\right] \quad (2) $$

and

$$ h_g(r) = h_{\text{ref}} \left(\frac{r}{R_0}\right)^{\beta}, \quad (3) $$
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** L- and N-band MATISSE aperture-synthesis images of FS CMa: top: Reconstructed L-band image (3.4–3.8 µm). The inner dark disk cavity in the L-band image has an angular diameter of ~6 × 12 mas (i.e., ~3.7 × 7.4 au for a distance of 620 pc). Bottom: Reconstructed N-band image (8.6–9.0 µm). The FOV is 64 mas (39.7 au for a distance of 620 pc). North is up and east to the left. The inserted circles represent the 2.7 and 6.6 mas diameters of the L- and N-band resolution beams. See Sect. 3 for more details.</figcaption></figure>

with the density scaling parameter $\rho_0$ determined by the total disk mass, $R_0$ the radius of the reference scale height $h_{ref}$, and $R_{trunc}$ the truncation radius at which the exponential surface-density drop-off sets in. The parameters $\beta$ and $\gamma$ determine the disk flaring and radial density profile.

**Dust properties.** Given the uncertain nature of FS CMa – a pre-main-sequence star or an evolved star – we first tentatively adopt the mixture of astronomical silicate and graphite. If FS CMa is a pre-main-sequence star, the dust chemical composition of its disk should be similar to that of the interstellar

<table><caption>**Table 2.** Selected stellar and disk parameters derived from earlier observations of FS CMa.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Distance (pc) <sup>(1)</sup></td><td>d</td><td>620<sup>+41</sup><sub>-33</sub></td></tr><tr><td>Source temperature (K) <sup>(1)</sup></td><td>T<sub>eff</sub></td><td>16 500<sup>+3000</sup><sub>-750</sub></td></tr><tr><td>Source luminosity (L<sub>☉</sub>) <sup>(1)</sup></td><td>log L<sub>★</sub></td><td>2.88<sup>+0.32</sup><sub>-0.17</sub></td></tr><tr><td>Inclination (°) <sup>(2)</sup></td><td>i</td><td>45 ± 8</td></tr><tr><td>Position angle (°) <sup>(2)</sup></td><td>PA</td><td>76 ± 11</td></tr></tbody></table>

**References.** <sup>(1)</sup> Vioque et al. (2018); <sup>(2)</sup> Kluska et al. (2020).

medium. For this reason, we first use a dust mixture consisting of 62.5% astronomical silicate and 37.5% graphite, which follows the $\frac{1}{3}-\frac{2}{3}$ relation for parallel and perpendicular orientations by Draine & Malhotra (1993). On the other hand, if FS CMa is an evolved star, its disk should be oxygen-rich, consisting of silicate alone without graphite. We describe models with this dust chemistry later in Sect. 4.5. The grain size distribution is described by the typical power law $n(s) \propto s^{-q}$ of Mathis et al. (1977) with $q = 3.5$ and grain sizes between 5 nm and 250 nm. Since this also includes small dust grains that are unusual for stars with high uv emission, we make sure that all discussed models in the following sections have temperatures below the respective sublimation temperature. The impact of larger dust grains is also discussed later in Sect. 4.5. To derive the required optical properties, we apply the wavelength-dependent refractive indices by Draine & Lee (1984), Laor & Draine (1993), and Weingartner & Draine (2001).

**Radiative transfer.** To analyze the reconstructed images from the observations, we calculate corresponding model images based on the previously described model. For this purpose, we apply the publicly available 3D Monte-Carlo radiative transfer code POLARIS (Reissl et al. 2016)$^2$. POLARIS calculates the temperature distribution of the disk using the immediate correction method of Bjorkman & Wood (2001) together with the instant reemission method of Lucy (1999). A raytracing algorithm is then used to calculate the thermal emission of the dust grains based on the temperature distribution. Finally, the scattered stellar light and scattered reemitted radiation are added using an anisotropic Monte-Carlo scattering method.

## 4.2. Fitting process

Due to the higher spatial resolution of the reconstructed L-band image compared to the N-band image, the former has been chosen to guide the fitting. In order to constrain the model parameters, we first analyze their impact on the inner disk structures visible in the reconstructed images. The inner radius $R_{in}$, disk flaring $\beta$, scale height $h_{100}$ at 100 au, and inclination $i$ turned out to be the most crucial parameters for the observed appearance of the innermost disk region. In contrast to this, the outer radius $R_{out}$, truncation radius $R_{trunc}$, radial density parameter $\gamma$, and total dust mass $M_{dust}$ have hardly any impact and thus cannot be constrained. Therefore, we use parameter values which have been either constrained by former observations of the large-scale disk structure of FS CMa or can be considered as typical for this class of objects.

<footnote>$^2$ https://www1.astrophysik.uni-kiel.de/~polaris/</footnote>
<figure><img src="image_3.png" /></figure>

<table><caption>**Fig. 2.** Comparison of the observed L-band image with the image of the basic disk model 1: (left) intensity distribution of the L-band model disk image, (middle) overlay of the observed L-band image with a contour plot of the model image, (right) difference image (observed image minus model image). Contour lines show 50% (solid) and 10% (dashed) levels of the maximum flux. The best-fit parameter values are compiled in the table below. The central star is masked in reconstructed and model images to focus on the difference caused by the disk structure. The flux ratio of central object and disk is discussed in Sect. 4.5.</caption><thead><tr><td>Parameters of model 1</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Inner radius (au)</td><td>R<sub>in</sub></td><td>5</td></tr><tr><td>Outer radius (au)</td><td>R<sub>out</sub></td><td>1000</td></tr><tr><td>Truncation radius (au)</td><td>R<sub>trunc</sub></td><td>200</td></tr><tr><td>Total dust mass (M<sub>☉</sub>)</td><td>M<sub>dust</sub></td><td>10<sup>-3</sup></td></tr><tr><td>Disk flaring exponent</td><td>β</td><td>1.2</td></tr><tr><td>Radial density exponent</td><td>γ</td><td>1.1</td></tr><tr><td>Scale height (au)</td><td>h<sub>100</sub></td><td>20</td></tr><tr><td>Inclination (°)</td><td>i</td><td>42</td></tr><tr><td>Position angle (°)</td><td>PA</td><td>20</td></tr></tbody></table>

To quantify the goodness of the fit, we (1) calculate images of the difference between the reconstructed and simulated images. In addition, we (2) consider the absolute mean difference by averaging the absolute value of the difference images, as well as (3) the difference between observed and modeled visibilities and closure phases. Relative flux units are used in the reconstructed images and in the simulated images. In the next step, we mask the extended image of the central star in the reconstructed L-band image and in the simulated image to first focus on reproducing the brightness distribution of the disk. However, unmasked images of the simulation and reconstruction are shown later in Figs. 6 and 9. The masked simulated images are shifted within a range of a few pixels to find the best match with the reconstructed image (i.e., the minimum of the absolute mean difference). In a second step (Sect. 4.6), we compute the difference between model closure phases and observed closure phases and model visibilities with observed visibilities (called residuals below). The comparison of the total flux of the observed and simulated images is discussed in Sect. 4.4.

**Inner radius and inclination.** The appearance of the inner rim of the disk as seen in the L-band image is mainly determined by the inner radius and disk inclination. While the inner radius directly affects the size of the inner disk cavity, the inclination determines the eccentricity of the elliptical model image (see Fig. 2, left), with potential shadowing by outer disk regions and – in interplay with the anisotropic scattering function of the dust grains – the relative contributions of forward and back-scattered radiation. Both parameters are fit simultaneously.

We constrain the inner radius $R_{in}$ by varying the inner radius between 3 au and 7 au and calculating the absolute mean difference between the reconstructed and model images by averaging the absolute value of the difference images (see Fig. 2 right)

to find the best-fit image. To constrain the inclination, we perform a series of simulations for inclinations between 35° and 50°. These are again compared to the reconstructed image by calculating the absolute mean difference between the reconstructed and model images. Based on these simulations, we deduce first best-fit values for the inclination of $42 \pm 3°$ and the inner radius $R_{in} = 5 \pm 0.5$ au, respectively. We note, that the temperature of the dust of the rim does not exceed 1200 K. The inner radius is therefore not perfectly corresponding to the sublimation temperature of the dust composition, which is 1300–1500 K for silicate (see, e.g., Pollack et al. 1994) and 2500 K for graphite (see, e.g., Kobayashi et al. 2011).

**Scale height and flaring parameter.** While the inner radius and inclination are directly related to the inner ellipse seen in the L-band model, the scale height $h_{100}$ and the flaring parameter $\beta$ determine the radial extent of the bright inner rim. We conduct simulations with scale heights between 15 au and 30 au in steps of 5 au and compare the resulting radial extent of the inner wall with that of the reconstructed L-band image. Constraining the flaring parameter $\beta$ is less straightforward because it affects various features in the image. Here, we have to compare the overall shape of the simulated images with the reconstructed images. To cover a sufficiently wide range of flarings, we run simulations with a flaring parameter $\beta$ between 0.5 and 1.5 in steps of 0.1. By increasing the flaring we increase the flux on the far side of the disk. We find that at a certain degree of flaring the disk starts to cover parts of the inner region. Considering this and excluding unrealistically high and low scale heights, we constrain the flaring parameter to 1.2 with a best-fit scale height at 100 au of around 20 au.

Figure 2 presents a comparison of the best-fit basic model 1 with the observed image: (from left to right) the intensity
<figure><img src="image_4.png" /><figcaption>Fig. 3. Vertical cuts through the density distributions of the different shapes of the modified inner rim and unmodified (yellow) inner rim for comparison ($R_\epsilon$ fixed at $1.4 R_{in} = 7$ au). The shape parameter $\epsilon$ is increasing from left to right and top to bottom (0.25, 0.5, 0.7, 0.9, 1.0). The upper left image shows the unmodified inner rim described in Eq. (1) for reference.</figcaption></figure>

distribution of the simulated model 1 image, an overlay of the observed L-band image with a contour plot of the model image, as well as the corresponding difference image (the central star is masked in the figure).

It can easily be seen that the inner elliptical gap and the extent of the bright region are well fit. However, the difference image also shows that the brightness distribution does not fit the characteristic asymmetry of the brightness distribution of the disk rim.

### 4.3. Modified curved inner rim (models 2 and 3)

While the basic model 1 described above allows us to roughly reproduce the shape of the image of the disk as well as the radial extent of the bright region, it is not capable of reproducing the asymmetry of the ring-like brightness distribution, as the difference image in Fig. 2 shows. Thus, the shape of the inner rim has to be modified toward a lower vertical rim extent in our model. Besides the motivation derived from our observational finding, this modification is also in agreement with previous studies based on hydrodynamical simulations and analytical considerations for the inner rim of circumstellar disks (e.g., Isella & Natta 2005; McClure et al. 2013; Flock et al. 2016). To match the brightness distribution along the inner disk rim and to take the theoretical studies into account, we expand our model by a simple parameterization of differently curved shapes of the inner rim described by the modified scale height:

$$h_g(r) = h_0(R_\epsilon) \cdot \left( \frac{r - R_{in}}{R_\epsilon - R_{in}} \right)^{\epsilon}, \quad (4)$$

where $R_\epsilon$ is the radius up to which the rim is modified, $h_0(R_\epsilon)$ is the unmodified scale height of $R_\epsilon$, and $\epsilon$ is a parameter to describe the different shapes of the inner rim. For illustration, selected examples for the different inner rim shapes are shown in Fig. 3. We now fix the best-fit parameters described in the section above and determine the best-fit of the brightness distribution of the inner rim by sampling a parameter space of (1.2–2.4) $R_{in}$ for

<figure><img src="image_5.png" /><figcaption>Fig. 4. Absolute mean difference between simulated (model 3) and reconstructed L-band images for different inner rim parameters $\epsilon$ and $R_\epsilon$.</figcaption></figure>

$R_\epsilon$ in steps of 0.2 $R_{in}$ and 0.1–1.0 in steps of 0.1 for $\epsilon$, respectively. The resulting map of the absolute mean difference for each model can be found in Fig. 4. In the following, we shortly discuss the impact of the two rim parameters $\epsilon$ and $R_\epsilon$.

**Inner rim parameter $\epsilon$.** The parameter $\epsilon$ allows us to modify the shape and therefore illumination of the inner rim of the model. As the case $\epsilon = 0.1$ describes a vertical rim that extends up to $R_\epsilon$ with a constant height, increasing $R_\epsilon$ results in a worse fit (see Fig. 4). Increasing $\epsilon$ to 0.5 changes the structure of the inner rim toward a parabolic shape. This configuration does not only result in the required radially more efficient illumination of the inner disk region, but also allows the observer to see more of the inner region on the near side of the disk. Besides the qualitatively better match of the observed inner rim structure, the absolute mean difference between reconstructed and simulated images is decreasing. Values of $\epsilon$ up to 1 are describing the transition from a parabolic shape to a peaked, wedge-like inner rim, leading to an extension of the bright region to the outer parts.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison of the observed *L*-band image with the image of model 2: (left) best-fit *L*-band model disk image with modified inner rim, (middle) overlay of contour plot of model image and reconstructed *L*-band image, and (right) difference image. Contour lines show 50% (solid) and 10% (dashed) levels of the maximum flux. As outlined in Sect. 4.2, the given flux unit is relative. Rim parameter: $R_\epsilon = 1.8 R_{in} = 9$ au, $\epsilon = 0.7$. The central object is masked in reconstructed and model image to focus on the difference caused by the disk.</figcaption></figure>

The corresponding simulated images show a broad, extended ring that is inconsistent with the reconstructed images.

**Radius of modified inner rim R<sub>ε</sub>.** The radius of the modified inner rim describes the point up to which the basic model is modified. It is therefore directly related to the shape of the simulated image for different values of the parameter ε. In the case of low values of ε, for example, a vertical rim with constant height, an increase of R<sub>ε</sub> results in an increasing mismatch between the simulated and reconstructed images. A similar trend can be found for high values of ε, which result in a bright extended ring, as described above, that increases even more when increasing R<sub>ε</sub>. For shapes of the inner rim between parabolic and tapered, we see a decrease of the absolute difference up to a value of R<sub>ε</sub> = 1.8 R<sub>in</sub>. This corresponds to an extension of the bright region as well as a blurring of the bright rim. Values of R<sub>ε</sub> > 1.8 R<sub>in</sub> will not improve the model because the modified rim is extending seamlessly into the flared disk in this case, which means changes of R<sub>ε</sub> in this regime will not influence the shape of the inner rim. Values of R<sub>ε</sub> below 1.2 R<sub>in</sub> can also be neglected, because the modification of the rim would be too small to be noticed.

According to the absolute mean differences measured for the different parameter sets, we find best-fit values of $R_\epsilon = 1.8 R_{in}$ and $\epsilon = 0.7$ (model 2; Fig. 5). We note that the model with $\epsilon = 0.6$ is deviating less than 1% from the best-fit case and is even best-fitting for larger values of $R_\epsilon$. The simulated image, the contour plot as well as the plot of the differences for the best-fit case are shown in Fig. 5. It can be seen that the difference image now has smaller values. We note that the difference plot of observed image minus model image still shows two regions with higher differences than on average. In these regions, the flux in the simulated image is higher than the measured flux on the near side.

To make sure that the modified rim is not affecting the best-fit parameters derived in Sect. 4.2, we repeat the calculations for the entire parameter grid. Comparing the differences confirms the best-fit parameters depicted in the table of Fig. 2 except for the inclination for which we now find a best-fit value of 51°. The corresponding simulated image of this model 3, the contour plot as well as the plot of the differences are shown in Fig. 6.

As a further check, we simulate the corresponding *N*-band image of our model and the respective difference images. The best-fit parameter values in this case are $R_\epsilon = 1.8 R_{in}$ and $\epsilon = 0.9$. In contrast to the best-fit of the reconstructed *L*-band image, the simulated *N*-band images with $\epsilon = 0.7$ or 0.8 are of comparable goodness, as their respective absolute mean differences are also within the measured uncertainties and are deviating by only 3%.

The simulated image, contour plot and the plot of the differences between simulated and observed images of the *N*-band best-fit parameter values are shown in Fig. 7. The best-fit parameter values derived on the basis of the reconstructed *L*-band image fits the *N* band equally well. As in the case of the *L* band, the differences are on the order of the uncertainties of the reconstructed image. The trends found for the parameters describing the shape of the inner rim as depicted in Fig. 4 for the *L*-band image can also be confirmed with the *N*-band simulations. Thus, we consider the best-fit model 3 is consistent with the reconstructed *L*- and *N*-band images.

### 4.4. Total flux

We measure a total *L*-band model flux of 35 Jy in the model FOV of 80 mas, while we derive an observed flux of 47 Jy (estimated uncertainty about ±50%) from the FS CMa interferograms in their larger 0.6 arcsec FOV of the *L*-band interferograms. The observed flux was derived by comparing the integrated flux of the interferograms of the target and the calibrator stars. In the *N*-band, we obtain a total model flux of 110 Jy in the model FOV of 360 mas and derived an observed flux of 213 Jy (±50%) from the FS CMa *N*-band interferograms in their FOV of 2.3 arcsec. In the model, we are thus covering about 3/4 of the observed flux in the *L* band, but only 1/2 of the flux in the *N* band. The differences may be caused by the flux from regions outside the model FOV.

We note that the thermal model emission always outweighs the scattered contribution by a factor of 10 in the *L* band and a factor of 300 in the *N* band, respectively. A detailed image of the ratio between thermal and scattered emission for both bands can be found in Fig. 8.

### 4.5. Model image with pure astrosilicate and a dust envelope with constant density (model 4)

Model 3 presented in Fig. 6 in the Sect. 4.3 is able to reproduce the asymmetric structure of the inclined FS CMa disk. However, the brightness of the central object in the observed *L*-band image and the model image disagree. In the observed *L*-band image, the integrated *L*-band flux of the central object is about (4.5 ± 0.9)% of the total image flux, whereas in the model, the *L*-band flux from the central object is only 0.61%.

Even though the central star shows photometric variability of more than a magnitude within years as discussed for example in de Winter & van den Ancker (1997), the stellar parameters have been well constrained by multiple previous studies
<figure><img src="image_7.png" /><figcaption>**Fig. 8.** Ratio of thermal to scattered emission of the *L*-band (top) and *N*-band (bottom) model. Direct emission of the star is included in the scattered contribution. The ratio of the total thermal emission to total scattered and direct stellar emission is 10 in case of the *L*-band image and 380 in the case of the *N*-band image.</figcaption></figure>

modeled flux fraction of central object (0.61–1.1% of the total flux of central object plus dust disk) is lower than in the reconstructed image (where it is 4.5%). This disagreement may be caused by a compact gaseous disk surrounding the central star (Sect. 4.7). Because of this gaseous disk, the central object of FS CMa seems to be much more complicated than just a normal star. This is probably the reason, why our model, which does not include a gaseous disk, is not able to reproduce the correct flux ratio of central object and dust disk. Because of this uncertainty, our main goal is only to model the shape of the dust disk. Modeling of the entire complex system of dust disk and central object is out of the scope of this paper.

In this complex situation, it became clear that the just applied modeling criterion of using the difference images of model and reconstructed images has the advantage that we can find the best-fit disk image by ignoring the reconstructed image of the central object. For this goal, we first masked the central object in previous sections before we calculated the average absolute value of the difference images to obtain an integral measure of the observation-model difference. In the next Sect. 4.6, we use another important criterion to find the best-fit image, namely the differences between modeled and observed closure phases and visibilities.

### 4.6. Comparison of the closure phases and visibilities derived from the model intensity distribution to the measured closure phases and visibilities

Instead of calculating the differences between reconstructed and the modeled images to constrain model parameters (as in the

previous sections), we now compute the average absolute value of the difference between the observed and model visibilities and closure phases called “model-observation residuals”.

The model-observation residuals of the closure phases and the visibilities for the *L*-band model 3 image (Fig. 6) are 43° and 0.040, respectively. These big residuals are not surprizing because we know already from the previous sections that the model flux ratio of the central object and the disk disagrees with the observations.

These *L*-band model-observation residuals are very big in comparison to another important type of residuals, the “reconstruction-observation” residuals, which are the residuals between the visibilities and closure phases of the reconstructed image (Fig. 1) and the observed visibilities and closure phases. The reconstruction-observation residuals are a measure of the quality of the reconstructed image. The reconstruction-observation residuals for the *L*-band image in Fig. 1 are only 2.9° and 0.009 for closure phase and visibilities, respectively (see Table 1).

One of the reasons for the big aforementioned model-observation residuals (closure phase: 43°, visibility: 0.040) is the small flux fraction of the central object (only 0.61–1.1% of the total flux) of the models 3 and 4 (Figs. 6 and 9). Other reasons are most likely additional disk asymmetries or clumpiness, which are not included in the modeling. This strong dependence of the model-observation residuals on the too small flux fraction of the central object is problematic if we want to use the model-observation closure phase and visibility residuals to constrain the parameters of the model disk alone because these residuals depend on both the required dust disk parameters and the wrong flux ratio of central object and dust disk. We can constrain the disk parameters only if we avoid the influence of the wrong flux ratio of central model star and disk by compensating the wrong flux ratio, which can approximately be performed by applying a suitable flux scaling factor to the central model star. This scaling factor is only needed as a tool to evaluate the model-observation residuals of the disk and is not used in radiative transfer modeling. For these reasons, we first compute the dependence of the model-observation residuals on the scaling factor of the model central star. The results are shown in Figs. 10 and 11 for our two models 3 and 4 presented in Figs. 6 and 9.

Figure 10 shows the closure phase and visibility residuals of the model 3 in Fig. 6 (i.e., model-observation residuals) as a function of different scaling factors for the brightness of the model central star. The minimum closure phase residuals of 19.9° are achieved with a scaling factor of 9.4. Figure 11 shows the residuals of model 4 (modified model image displayed in Fig. 9) as a function of the scaling factors. In this case, the minimum closure phase residuals of 24.7° are obtained with a scaling factor of 3.9, whereas the closure phase residuals in the case of no scaling are 48.2°. The remaining closure phase residuals of ~25° are most likely caused by additional asymmetries in the dust disk, which are not taken into account in the models. This discrepancy is also shown in Fig. C.2, where the measured interferometric data are plotted together with the interferometric data derived from the two models 3 and 4, but modified with the above model star scaling factors 9.4 and 3.9, respectively.

In addition to the above model-observation residual method, we used the $\chi^2$ method to determine the average deviation of the model visibilities and closure phases from the observed values. In Appendix F, we show that the $\chi^2$ method provides similar scaling factors of 9.0 and 3.7 (instead of the above 9.4 and 3.9) for the two models 3 and 4, respectively.
<figure><img src="image_8.png" /></figure>

**Fig. 9.** Comparison of the observed *L*-band image with the image of model 4: best-fit *L*-band model image with pure astrosilicate and a dust envelope with constant density and $\tau = 0.1$ (left), contour plot over reconstructed *L*-band image (middle) and difference image (right). Contour lines show 50% (solid) and 10% (dashed) levels of the maximum flux. As outlined in Sect. 4.2, the given flux unit is relative. Rim parameter: $R_{\epsilon} = 1.8 R_{\text{in}} = 9$ au, $\epsilon = 0.7$. Inclination: $i = 51^{\circ}$.

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

**Fig. 10.** Closure phase and visibility residuals of the best-fit *L*-band model 3 shown in Fig. 6 as function of the central star scaling factor.

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

**Fig. 12.** Closure phase and visibility residuals of the modified best-fit *N*-band model shown in Fig. 7 as a function of the central star scaling factor.

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

**Fig. 11.** Closure phase and visibility residuals of the modified best-fit *L*-band model 4 shown in Fig. 9 as function of the central star scaling factor.

The model star scaling factors derived from the closure phase residuals are approximately consistent with the modeling results obtained when we used the absolute mean differences between the observed and modeled images as criterion. In this case, the model images have a flux fraction in the central object of only 0.61% (model 3) or 1.1% (model 4) of the total flux, whereas the flux fraction of the central object is (4.5 ± 0.9)% in the observed *L*-band image, corresponding to ratio of 4.5/0.61 = 7.4 (model 3) and 4.5/1.1 = 4.1 (model 4). These model-reconstruction deviations are similar to the model central star

scaling factors of 9.4 (model 3) and 3.9 (model 4) derived from the closure phase residuals.

Figure 12 shows the model-observation residuals of the *N*-band model 3 image displayed in Fig. 7 as a function of the different scaling factors of the central star. The minimum closure phase residuals of 12.19° are obtained with a central star scaling factor of 9.2, whereas the closure phase residuals in the case of no scaling are 12.23°. The difference between scaling and no scaling of the central star is very small because the integrated stellar intensity is only 0.03% of the total flux of the best-fit *N*-band model. In Fig. C.4, the measured interferometric data are plotted together with the interferometric data derived from the *N*-band reconstruction shown in Fig. 1 and derived from the model displayed in Fig. 7 (model 3), but modified with the above scaling factor 9.2.

We have shown by three different methods that the model flux ratio of the central star and the total flux (of disk and central object) is about 3–9 times too small. In Sect. 4.7, we discuss a likely cause for this discrepancy. As we discussed in previous sections, our goal is only to constrain model parameters of the dust disk. The modeling of the complex star-disk system is out of scope of this paper because of the complex central object discussed in Sect. 4.7.

Therefore, we focus now on the application of the model-observation closure phase residuals and the aforementioned $\chi^2$ results to constrain the same disk shape parameters $R_{\epsilon}$ and $\epsilon$,
<figure><img src="image_12.png" /><figcaption>**Fig. 13.** Closure phase residuals of L-band models (different versions of model 3) with different inner rim parameters $\epsilon$ and $R_\epsilon$ as a function of the central star scaling factor. The best-fit models are those two with $\epsilon = 0.6$ and $R_\epsilon = 1.8$ and $2.2 R_{in}$ at the central star scaling factor 8.4.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 14.** Model-observation closure phase residuals of L-band models (different versions of model 3) for different inner rim parameters $\epsilon$ and $R_\epsilon$ at the best-fit scaling factor 8.4 of the central star.</figcaption></figure>

which were constrained in Sect. 4.3 using the presented difference between the reconstructed image and the model image. A third required free modeling parameter is the aforementioned scaling factor of the brightness of the central star because of the derived flux ratio discrepancy.

Figure 13 shows the model-observation closure phase residuals of several versions of model 3 (i.e., $R_\epsilon$—$\epsilon$ pairs with $R_\epsilon = 1.8, 2.0, 2.2$ and $\epsilon = 0.5, 0.6, 0.7$) as a function of the central star scaling factor. The best-fit models, with the smallest closure phase residuals of $16.5^\circ$, have a central star scaling factor of $8.4, \epsilon = 0.6$, and $R_\epsilon = 1.8$ and $2.2$. Furthermore, the $\chi^2$ method presented in Appendix F provides the same disk parameters $\epsilon$ and $R_\epsilon$ as derived with the residual method applied in this section. Figures 14 and 15 show the same type of $R_\epsilon$—$\epsilon$ plot as in Fig. 4, but now using the model-observation closure phase residuals or $\chi^2$ results, instead of the absolute mean differences between simulated and reconstructed images. The results of the three fits (in Figs. 4, 14, and 15) are similar. They differ in the inner rim parameter value $\epsilon$ only, 0.6 in Figs. 14 and 15 instead of 0.7 in Fig. 4. In all models presented in this study, the dust disk temperature is below the sublimation temperature. Figures 10 and 11 suggest that model 3 can reproduce the disk structure slightly better than model 4 because the closure phase residuals of model

3 at minimum scaling factor are smaller (~20°; Fig. 10) than the closure phase residuals of model 4 (Fig. 11; ~25°).

### 4.7. Influence of an additional gaseous B[e] star disk on the model brightness ratio between the central object and the dust disk of FS CMa

In the models discussed above, the model image of the central object is always too faint with respect to the dust disk. If we assume that the central star is surrounded by an additional compact (i.e., unresolved) gaseous disk, the model brightness ratio between central object (= star plus gaseous disk) and dust disk would change, as discussed below. The theory of gaseous envelopes and disks of Be stars is discussed in Lamers & Waters (1984), Waters (1986), Lee et al. (1991), Carciofi (2011), Rivinius et al. (2013), Vieira et al. (2015, 2016). The gaseous disks of several Be and B[e] stars were studied with high spatial resolution using optical and infrared long-baseline interferometers (e.g., Thom et al. 1986; Mourard et al. 1989; Quirrenbach et al. 1993, 1997; Stee et al. 1995; Berio et al. 1999; Tycner et al. 2006; Domiciciano de Souza et al. 2007; Meilland et al. 2009; Millour et al. 2011; Touhami et al. 2011; Kraus et al. 2012).

Below we discuss that the infrared continuum flux from a gaseous disk can be several times brighter than the stellar flux, and the brightness of the gaseous disk depends on the wavelength and the viewing angle (e.g., it is brightest for face-on viewing direction; Vieira et al. 2015, 2016). Because of this viewing angle dependence, a gaseous disk can change the brightness ratio between the central object (= star plus gaseous disk) and the dust disk of FS CMa. A gaseous disk can much increase the flux from the central object seen by the observer (inclination angle ~51°). However, the illumination of the dust disk by the gaseous disk may be smaller if the gaseous and dust disk are aligned.

The dependence of the intensity ratio of total flux to stellar flux on the inclination angle, density, and wavelength can be calculated using the DISCO code<fnref n="3" /> based on the Vieira et al. (2015, 2016) model for classical Be star gaseous disks. The stellar parameters used are given in Table 2. The value of the inclination angle is assumed to be 51°, as derived during the modeling process. The parameters of the gaseous disk model are given in Table 3.

Figure 15 shows the dependence of the model intensity ratio of total flux (i.e., stellar plus gaseous disk flux) to the stellar flux on wavelength for different gas densities at the base of the gaseous disk and for the inclination angle of 51° of FS CMa. Figure 16 shows the intensity ratio of the total flux (i.e., stellar plus gaseous disk flux) to the stellar flux as a function of the wavelength for different inclination angles for a gas density at the base of the gaseous disk of $10^{-7} \text{ kg m}^{-3}$. For example, at L-band wavelengths, the total flux is about 3–4 times larger than the stellar flux for the inclination range of 0–51° (0° corresponds to face-on) of the gaseous disk in contrast to larger inclination angles. This inclination angle dependence of the brightness of such a gaseous disk would have a strong influence on the FS CMa model image if FS CMa is surrounded by a bright gaseous disk and the gaseous disk is aligned with the dust disk. In this case, the gaseous disk is seen edge-on (i.e., low brightness) from viewing locations within the dust disk, whereas we see a bright 51° inclined disk. Therefore, such a gaseous disk may explain the difference between the modeled and observed flux ratio of the

<footnote n="3">³ https://amhra.jmmc.fr/</footnote>
<table><caption>**Table 3.** Parameters of the modeled gaseous disk.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Disk outer radius (R<sub>☉</sub>)</td><td>50</td></tr><tr><td>Disk temperature (K)</td><td>9900</td></tr><tr><td>Density at the disk base (kg m<sup>-3</sup>)</td><td>0.5–2.5 × 10<sup>-7</sup></td></tr><tr><td>Density power law exponent</td><td>-3.5</td></tr><tr><td>Power law coefficient of disk flaring</td><td>1.5</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 15.** Wavelength dependence of the intensity ratio of the total flux (stellar plus gaseous disk) to the stellar flux for the inclination angle of 51° for different gas densities at the base of the gaseous disk.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 16.** Wavelength dependence of the intensity ratio of total flux (stellar plus gaseous disk) to the stellar flux for different inclination angles with a gas density at the base of the gaseous disk of $10^{-7}$ kg m$^{-3}$.</figcaption></figure>

central object and the dust disk because a gaseous disk is not taken into account in our radiative transfer modeling.

## 5. Summary and conclusions

FS CMa belongs to the subgroup of unclassified B[e] stars. It is not known whether FS CMa is a young Herbig B[e] or a star in a later evolutionary stage. Vioque et al. (2020) concluded that FS CMa is not a young star. FS CMa is famous for its large, inclined circumstellar disk and spectacular bipolar outflow perpendicular to the disk plane.

We performed aperture-synthesis imaging of FS CMa with the new VLTI-MATISSE interferometry instrument and carried out radiative transfer modeling to interpret the images. The following results were obtained:

- We reconstructed the first aperture-synthesis L- and N-band images of the stellar disk and the central object. Image reconstruction was performed with the IRBis method. The L- and N-band images have resolutions of 2.7 mas and 6.6 mas, respectively;

- In the L-band image, the inner rim region of the circumstellar dust disk and the central star can be seen. The size of the resolved inner dark disk cavity is ~6×12 mas. The resolved L-band disk consists of a bright NW disk region and a much fainter SE one. The images suggest that we are looking at the bright inner wall of the NW disk rim, while the SE rim is partially self-shadowed and appears therefore dimmer in the image. Therefore, the large brightness of NW rim suggests that the NW rim is on the far side. The N-band image shows only the bright NW disk region;

- The reconstructed images provide the basis for a detailed analysis of the disk structure, and thus the spatial distribution of the dust in the innermost disk region. Performing radiative transfer simulations, we derived an inner disk radius of 5 au and a disk inclination of 51°. Most importantly, the good quality of the reconstructed images, in particular the high-resolution L-band image, allowed us to constrain the shape of the vertical structure of the inner few au of the inner disk. Using a parameterized disk model, the shape of the rim is described by two parameters characterizing a) the curvature of the disk rim and b) the radial range beyond which the density profile is described by the global disk model. In particular, we find a rather flat shape of the disk between the inner disk radius and 9 au. Regardless of the specific evolutionary state of FS CMa, this finding provides direct constraints for further studies of the physical mechanisms resulting in this particular shape, that is, the interaction of the stellar radiation with the gas and dust disk surrounding the central star. However, a detailed analysis taking the specifics of FS CMa into account is beyond the scope of this study;

- While the relative brightness distribution can be well fitted, the flux ratio between the central object and the net flux of the disk cannot be matched on the basis of a star-to-dust disk model alone. Depending on the specific model setup, characterized by the chemical composition and size distribution of the dust, even taking a hypothetical envelope with a low optical depth as a source of possible contamination into account, the latter results in a star-to-disk flux ratio which is too low. All model images have a flux fraction in the central object of only ~0.61–1.1% of the total flux, whereas the flux fraction of the central object is ~4.5% in the observed L-band image. To explain this discrepancy, we suggest the possible additional contribution of a compact gaseous B[e] star disk in the vicinity of the star. Taking the dependence of the brightness on the viewing angle into account, we find that this gaseous disk may be sufficiently bright to explain the observed flux ratio. More detailed studies of FS CMa are needed in the future to improve our understanding of its evolutionary stage, suspected companion, brightness variations, and physical properties of its suggested gaseous disk.

*Acknowledgements.* MATISSE was designed, funded and built in close collaboration with ESO, by a consortium composed of institutes in France (J.-L. Lagrange Laboratory – INSU-CNRS – Côte d’Azur Observatory – University of Côte d’Azur), Germany (MPIA, MPIfR and University of Kiel), the Netherlands (NOVA and University of Leiden), and Austria (University of Vienna). We thank all ESO colleagues for the excellent collaboration. This work has been supported by the French government through the UCAJEDI Investments in the Future project managed by the National research Agency (ANR) with
the reference number ANR-15-IDEX-01. The Konkoly Observatory and Cologne
University have also provided some support in the manufacture of the instru-
ment. K.O. acknowledges the support of the Agencia Nacional de Investigación
y Desarrollo (ANID) through the FONDECYT Regular grant 1180066 and
1210652. S.K. and A.K. acknowledge support from an STFC Consolidated Grant
(ST/V000721/1) and ERC Starting Grant (Grant Agreement No. 639889). The
research of J.V. and M.H. is supported by NOVA, the Netherlands Research
School for Astronomy. T.H. acknowledges support from the European Research
Council under the Horizon 2020 Framework Program via the ERC Advanced
Grant Origins 83 24 28. A.G. acknowledges support from the European Research
Council (ERC) under the European Union's Horizon 2020 research and inno-
vation programme under grant agreement No 695099 (project CepBin). P.A.
acknowledges support from the Hungarian NKFIH OTKA grant K132406, and
from the European Research Council (ERC) under the European Union's Hori-
zon 2020 research and innovation programme under grant agreement No 716155
(SACCRED). J.S.B. acknowledges the support received from the UNAM PAPIIT
project IA 101220 and from the CONACyT project 263975. AD thanks sup-
port from FAPESP #2019/02029-2. F.N. acknowledges FAPESP for support
through proc. 2017/18191-8. This research has made use of the services of
the ESO Science Archive Facility. This publication makes use of the SIM-
BAD database operated at CDS, Strasbourg, France. This research has made
use of the Jean-Marie Mariotti Center (JMMC) – MOIO AMHRA service at
https://amhra.jmmc.fr/. We thank the referee for helpful suggestions.

References

Baines, D., Oudmaijer, R. D., Porter, J. M., & Pozzo, M. 2006, MNRAS, 367, 737

Berio, P., Stee, P., Vakili, F., et al. 1999, A&A, 345, 203

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

Bourges, L., Mella, G., Lafrasse, S., et al. 2017, VizieR Online Data Catalog, II/346

Brunngräber, R., Wolf, S., Ratzka, T., & Ober, F. 2016, A&A, 585, A100

Byrd, R. H., L. P., & Nocedal, J. 1995, SIAM J. Sci. Stat. Comput., 16, 1190

Carciofi, A. C. 2011, in Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 272, 325

Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23

de Winter, D., & van den Ancker, M. E. 1997, A&AS, 121, 275

Domiciano de Souza, A., Driebe, T., Chesneau, O., et al. 2007, A&A, 464, 81

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Draine, B. T., & Malhotra, S. 1993, ApJ, 414, 632

Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144

Gravity Collaboration (Perraut, K., et al.) 2019, A&A, 632, A53

Hager, W. W., & Zhang, H. 2005, SIAM J. Optim., 16, 170

Hager, W. W., & Zhang, H. 2006, SIAM J. Optim., 17, 526

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Hofmann, K.-H., Weigelt, G., & Schertl, D. 2014, A&A, 565, A48

Isella, A., & Natta, A. 2005, A&A, 438, 899

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Kluska, J., Malbet, F., Berger, J. P., et al. 2014, A&A, 564, A80

Kluska, J., Berger, J. P., Malbet, F., et al. 2020, A&A, 636, A116

Kobayashi, H., Kimura, H., Watanabe, S. i., Yamamoto, T., & Müller, S. 2011, Earth, Planets Space, 63, 1067

Kraus, S., Monnier, J. D., Che, X., et al. 2012, ApJ, 744, 19

Lamers, H. J. G. L. M., & Waters, L. B. F. M. 1984, A&A, 136, 37

Lamers, H. J. G. L. M., Zickgraf, F. J., de Winter, D., Houziaux, L., & Zorec, J. 1998, A Revised Classification Scheme for B[e] Stars, eds. A. M. Hubert, & C. Jaschek, 233, 277

Laor, A., & Draine, B. T. 1993, ApJ, 402, 441

Laws, A. S. E., Harries, T. J., Setterholm, B. R., et al. 2020, ApJ, 888, 7

Lazareff, B., Berger, J. P., Kluska, J., et al. 2017, A&A, 599, A85

Lee, U., Osaki, Y., & Saio, H. 1991, MNRAS, 250, 432

Lopez, B., Lagarde, S., Jaffe, W., et al. 2014, The Messenger, 157, 5

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

Lucy, L. B. 1999, A&A, 344, 282

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

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

Matter, A., Lopez, B., Antonelli, P., et al. 2016, SPIE Conf. Ser., 9907, 99070A

McClure, M. K., D'Alessio, P., Calvet, N., et al. 2013, ApJ, 775, 114

Meilland, A., Stee, P., Chesneau, O., & Jones, C. 2009, A&A, 505, 687

Menu, J., van Boekel, R., Henning, T., et al. 2015, A&A, 581, A107

Millour, F., Petrov, R. G., Vannier, M., & Kraus, S. 2008, SPIE Conf. Ser., 7013, 70131G

Millour, F., Meilland, A., Chesneau, O., et al. 2011, A&A, 526, A107

Monnier, J. D., Berger, J. P., Millan-Gabet, R., et al. 2006, ApJ, 647, 444

Mourard, D., Bosc, I., Labeyrie, A., Koechlin, L., & Saha, S. 1989, Nature, 342, 520

Petrov, R. G., Allouche, F., Berio, P., et al. 2018, SPIE Conf. Ser., 10701, 1070109

Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615

Quirrenbach, A., Hummel, C. A., Buscher, D. F., et al. 1993, ApJ, 416, L25

Quirrenbach, A., Bjorkman, K. S., Bjorkman, J. E., et al. 1997, ApJ, 479, 477

Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87

Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69

Sauter, J., Wolf, S., Launhardt, R., et al. 2009, A&A, 505, 1167

Schöller, M. 2007, New Astron. Rev., 51, 628

Stee, P., de Araujo, F. X., Vakili, F., et al. 1995, A&A, 300, 219

Thom, C., Granes, P., & Vakili, F. 1986, A&A, 165, L13

Touhami, Y., Gies, D. R., & Schaefer, G. H. 2011, ApJ, 729, 17

Tycner, C., Gilbreath, G. C., Zavala, R. T., et al. 2006, AJ, 131, 2710

Vieira, R. G., Carciofi, A. C., & Bjorkman, J. E. 2015, MNRAS, 454, 2107

Vieira, R. G., Carciofi, A. C., & Bjorkman, J. E. 2016, in ASP Conf. Ser., 506, Bright Emissaires: Be Stars as Messengers of Star-Disk Physics, eds. T. A. A. Sigut, & C. E. Jones, 135

Vioque, M., Oudmaijer, R. D., Baines, D., Mendigutía, I., & Pérez-Martínez, R. 2018, A&A, 620, A128

Vioque, M., Oudmaijer, R. D., Schreiner, M., et al. 2020, A&A, 638, A21

Waters, L. B. F. M. 1986, A&A, 162, 121

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Wolf, S., Padgett, D. L., & Stapelfeldt, K. R. 2003, ApJ, 588, 373

Zhu, C., B. R. H., & Nocedal, J. 1997, ACM Trans. Math. Softw., 23

Zickgraf, F. J. 1998, Current Definition of B[e] Stars, eds. A. M. Hubert, & C. Jaschek, 233, 1

1 Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany
e-mail: khh@mpi.fr-bonn.mpg.de

2 University of Kiel, Institute of Theoretical Physics and Astro-physics, Leibnizstrasse 15, 24181 Kiel, Germany

3 Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Boulevard de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

4 Institute for Mathematics, Astrophysics and Particle Physics, Radboud University, PO Box 9010, MC 62, 6500 GL Nijmegen, The Netherlands

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

6 School of Physics, Astrophysics Group, University of Exeter, Stocker Road, Exeter EX4 4QL, UK

7 Departamento de Ciencias Físicas, Facultad de Ciencias Exactas, Universidad Andrés Bello, Fernández Concha 700, Las Condes, Santiago, Chile

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

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

10 European Southern Observatory, Casilla 19001, Santiago 19, Chile

11 European Southern Observatory, Karl Schwarzschild Strasse 2, 85748 Garching, Germany

12 ASTRON, Dwingeloo, Netherlands

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

14 Univ. Grenoble Alpes, CNRS, IPAG, 38000, Grenoble, France

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

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

17 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Konkoly Thege Miklósút 15-17, 1121 Budapest, Hungary

18 Astrophysics Science Division, NASA/GSFC, Greenbelt, MD 20771, USA

19 I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937, Köln, Germany
<figure><img src="image_16.png" /><figcaption>Fig. A.1. uv coverage of the FS CMa observations.</figcaption></figure>

# Appendix A: Observations

The VLTI-MATISSE observations of FS CMa are summerized in Table A.1. The uv coverage is shown in Fig. A.1. The calibrator stars used in this project are listed in Table A.2.

# Appendix B: Estimation of the visibility and closure phase errors and selection of the calibrated interferometric data

## Appendix B.1: L-band

Figure B.1 shows examples of our obtained calibrated visibilities and closure phases (CPs). All six data sets of an observation discussed in Sect. 2.1 are displayed. To estimate the visibility and CP errors, we calculated the standard deviations of the visibilities and CPs of the six calibrated interferometric data sets of each observation (see Sect. 2.1) for each spectral channel. This is possible, because the six calibrated interferometric data sets of each observation are recorded within about 15 min, and therefore are covering nearly the same projected baseline vectors. The standard deviations were used as the common errors of the six calibrated interferometric data sets of each observation.

To obtain a reliable data set for image reconstruction, we eliminated all CPs with errors >60° and all visibilities with a signal-to-noise ratio <0.5. In a second selection, the reconstructed four CPs of each single observation were tested if they are affected by systematic errors. This test is called CP test and is discussed below in more detail. The four CPs measured simultaneously with a 4-telescope array are not independent, that is, each one of the four CPs can be calculated by a linear combination of the three other CPs (see, for example, Chael et al. 2018). The distance between each measured CP and the corresponding one derived from the other three simultaneously measured CPs was calculated. This distance is a linear combination of the four CP deviations (e.g., the deviation between the measured CP and its true value). In many cases this distance is small, that is, there is a good agreement between the measured CPs and the CPs derived from the three other CPs. But there may be some observations where this distance is large. The reason

for this could be, for example, a nonoptimal noise bias compensation in the average bispectrum. Figure B.2 shows an example observation, where the measured CPs and the CPs derived from the other three ones are very similar with respect to their noise. In Fig. B.3 an observation is shown with a strong mismatch between the measured CPs and the derived CPs. In a second selection, we eliminated those observations where the mean deviation between the measured CPs and the derived CPs are larger than ~8°.

The feasibility of the aforementioned closure phase test was investigated by a computer simulation. The theoretical target was a 10 Jy source of similar structure, size, and sky position as FS CMa. The bispectrum of that target was simulated as if it were derived from data observed with the large AT configuration A0-G1-J2-K0 in L-band at low spectral resolution. The bispectrum was simulated with the noise bias terms caused by the noise of a 10 Jy source and strong sky background, but without noise in the bispectrum. Figure B.4 shows such a closure phase test for the simulated bispectrum with different levels of noise bias compensation (100% corresponds to perfect compensation). Because of the residual noise bias, the reconstructed CPs systematically deviate from their theoretical values. The deviations between the measured CPs and the theoretical ones is zero in the case of a perfect subtraction of the noise bias term. Figure B.4 shows that the distance between the measured CP and the same one derived from the three other CPs, as discussed above (e.g., the CP test, and denoted by $|CP - CP_{calc}|$ in Fig. B.4) is comparable to the distances between the measured and theoretical (true) CPs (denoted by $|CP_{theo} - CP_{meas}|$ in Fig. B.4). This experiment shows that systematic errors of the CP can roughly be determined by such a CP test.

## Appendix B.2: N-band

The visibility and CP errors for the calibrated interferometric N-band data, obtained with spectral binning of 11, were derived in the same way as in the L band: the standard deviations derived from the four calibrated interferometric data sets (one for each of the four BCD configurations) of each observation (see Sect. 2.2) were calculated for each spectral channel. The derived standard deviations were used as the common errors of the four interferometric data sets of each observation.

The data set used for image reconstruction was selected in the same way as in the case of the L-band data. CPs with errors > 60° and visibilities ($V^2$) with signal-to-noise ratio < 0.5 were not used for the image reconstruction. In a second selection, the four CPs of each single observation were tested if they are affected by systematic errors. This was done with the CP test described in Appendix B.1 for the L-band data. In the second selection, we eliminated all observations where the mean deviation between the measured CPs and the derived ones are larger than ~1.7× the mean CP error of the corresponding observation.

Figure B.5 shows examples of calibrated visibilities and CPs of a typical observation without and with different spectral binnings. The data without spectral binning (red lines) show very strong noise, which is strongly reduced by spectral binning (blue, green, and brown lines corresponding to spectral binnings of 9, 11, and 17, respectively).
2018-12-09T08.16.26-07.53_28averaged

<figure><img src="image_17.png" /><figcaption>**Fig. B.1.** Examples of calibrated visibilities and CPs of FS CMa L-band observation: all six interferometric data sets of an observation are shown (see Sect. 2.1 for more details).</figcaption></figure>
<table><caption>**Table A.1.** Summary of the VLTI-MATISSE observations of FS CMa.</caption><thead><tr><th>date</th><th>n<sup>a</sup></th><th>UT start<sup>b</sup></th><th>UT end<sup>c</sup></th><th>array</th><th>JD</th><th>seeing<sup>d</sup></th><th>τ<sub>0</sub><sup>e</sup></th></tr></thead><tbody><tr><td>2018-12-02</td><td>3</td><td>06:25</td><td>08:18</td><td>A0 B2 J2 C1</td><td>2458454.277</td><td>1.10</td><td>3.47</td></tr><tr><td>2018-12-03</td><td>2</td><td>04:00</td><td>05:13</td><td>A0 B2 J2 C1</td><td>2458455.178</td><td>0.49</td><td>4.73</td></tr><tr><td>2018-12-05</td><td>2</td><td>07:14</td><td>07:36</td><td>A0 B2 J2 C1</td><td>2458457.308</td><td>1.17</td><td>2.44</td></tr><tr><td>2018-12-06</td><td>5</td><td>04:00</td><td>08:33</td><td>A0 B2 J2 C1</td><td>2458458.333</td><td>0.53</td><td>4.07</td></tr><tr><td>2018-12-07</td><td>5</td><td>02:44</td><td>07:28</td><td>A0 B2 D0 C1</td><td>2458459.170</td><td>1.25</td><td>2.28</td></tr><tr><td>2018-12-08</td><td>4</td><td>04:33</td><td>08:46</td><td>A0 B2 D0 C1</td><td>2458460.296</td><td>0.61</td><td>3.76</td></tr><tr><td>2018-12-09</td><td>4</td><td>06:06</td><td>09:03</td><td>K0 B2 D0 J3</td><td>2458461.314</td><td>0.46</td><td>4.41</td></tr><tr><td>2018-12-10</td><td>6</td><td>04:48</td><td>08:40</td><td>K0 G2 D0 J3</td><td>2458462.219</td><td>0.54</td><td>8.12</td></tr><tr><td>2018-12-11</td><td>2</td><td>07:23</td><td>08:24</td><td>K0 G2 D0 J3</td><td>2458463.313</td><td>0.60</td><td>12.66</td></tr><tr><td>2018-12-12</td><td>3</td><td>05:07</td><td>08:20</td><td>A0 G1 D0 J3</td><td>2458464.290</td><td>0.71</td><td>4.82</td></tr><tr><td>2018-12-13</td><td>4</td><td>04:58</td><td>08:11</td><td>A0 G1 J2 J3</td><td>2458465.254</td><td>0.87</td><td>5.63</td></tr><tr><td>2018-12-14</td><td>2</td><td>05:17</td><td>06:19</td><td>A0 G1 J2 J3</td><td>2458466.226</td><td>0.74</td><td>4.87</td></tr><tr><td>2018-12-15</td><td>2</td><td>04:52</td><td>07:45</td><td>A0 G1 J2 K0</td><td>2458467.332</td><td>0.53</td><td>5.49</td></tr></tbody></table>

**Notes.** <sup>a</sup> Number of data sets recorded at different hour angles per night. <sup>b</sup> UT time at the start of the first data set of the night. <sup>c</sup> UT time at the start of the last data set of the night. <sup>d</sup> Average visible seeing ("" during the night). <sup>e</sup> Average coherence time $\tau_0$ (ms) in the visible during the night.

<figure><img src="image_18.png" /><figcaption>**Fig. B.2.** CP test: good agreement between the measured CPs and the CPs calculated by the other three CPs.</figcaption></figure>

<table><caption>**Table A.2.** Interferometric calibrator stars (JMMC Stellar Diameters Catalog - JSDC, Version 2; Bourges et al. 2017).</caption><thead><tr><th>name</th><th>spectral type</th><th>UDD-L<sup>a</sup></th><th>UDD-N<sup>b</sup></th></tr></thead><tbody><tr><td>7 Cet</td><td>M1III</td><td>5.50 ±0.55 mas</td><td>5.52 ±0.55 mas</td></tr><tr><td>ε Lep</td><td>K4III</td><td>5.87 ±0.56 mas</td><td>5.92 ±0.56 mas</td></tr><tr><td>V407 Pup</td><td>M2Ia/ab</td><td>2.60 ±0.25 mas</td><td>2.60 ±0.25 mas</td></tr><tr><td>β Col</td><td>K1III CN</td><td>3.61 ±0.41 mas</td><td>3.64 ±0.41 mas</td></tr><tr><td>π Eri</td><td>M1III</td><td>5.28 ±0.47 mas</td><td>5.30 ±0.47 mas</td></tr><tr><td>θ CMa</td><td>K3/4III</td><td>4.00 ±0.39 mas</td><td>4.03 ±0.39 mas</td></tr></tbody></table>

**Notes.** <sup>a</sup> L-band uniform disk diameter. <sup>b</sup> N-band uniform disk diameter.

## Appendix C: Comparison of measured CPs and visibilities with the same ones derived from the reconstructions and from the models

## Appendix D: Image reconstruction

The current version of the MATISSE image reconstruction package IRBis (mat_cal_imarec) consists of two optimization routines: a) a large-scale, bound-constrained nonlinear optimization routine called ASA_CG (engine 1; Hager & Zhang (2005), Hager & Zhang (2006)), and b) the new built-in optimization routine, a limited memory algorithm for bound-constrained optimization called L-BFGS-B (engine 2; Byrd & Nocedal (1995), Zhu & Nocedal (1997)).

The version of IRBis described in Hofmann et al. (2014) contains two cost functions: $Q_1[o_k(\mathbf{x})]$ (cost fct. 1) corresponding to the $\chi^2$ function of the bispectrum, and $Q_2[o_k(\mathbf{x})]$ (cost fct. 2) cor-
<figure><img src="image_19.png" /><figcaption>Fig. B.3. CP test: strong mismatch between the measured CPs and the CPs computed by the other three CPs.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. B.4. Computer simulation of the CP test, which shows the dependence of all CP distances on the percentage of noise bias compensation. The distance between the measured CP and the one derived from the three other CPs (|CP - CPcalc|; red line) is comparable to the deviations between the four measured and theoretical CP (|CPtheo - CPmeas|; blue lines; see text in Appendix B.1 for details).</figcaption></figure>

responding to the $\chi^2$ function of the bispectrum phasors plus the bispectrum moduli. $o_k(x)$ denotes the current iterated image and $x$ is a two-dimensional vector in image space. The actual version of IRBis additionally contains a third cost function $Q_3[o_k(x)]$ (cost fct. 3), which is the $\chi^2$ function of the bispectrum phasors

plus the squared visibilities, given by

$$Q_3[o_k(\mathbf{x})] := \\ = \int_{\mathbf{f}_1, \mathbf{f}_2 \in \mathcal{M}} \frac{w_d(\mathbf{f}_1, \mathbf{f}_2)}{\sigma_{\beta}^2(\mathbf{f}_1, \mathbf{f}_2)} \cdot |\gamma_0 \exp[i\beta_k(\mathbf{f}_1, \mathbf{f}_2)] - \\ - \exp[i\beta(\mathbf{f}_1, \mathbf{f}_2)]|^2 d\mathbf{f}_1 d\mathbf{f}_2 + \\ + f_0 \cdot \int_{\mathbf{f} \in \mathcal{M}_p} \frac{w_d^p(\mathbf{f})}{\sigma_p^2(\mathbf{f})} \cdot |\gamma_0 |O_k(\mathbf{f})|^2 - |O(\mathbf{f})|^2|^2 d\mathbf{f}. \quad (D.1)$$

Here, $\beta(\mathbf{f}_1, \mathbf{f}_2)$ and $\beta_k(\mathbf{f}_1, \mathbf{f}_2)$ denote the measured CPs and the CPs of the current iterated image, respectively. The functions $\sigma_{\beta}(\mathbf{f}_1, \mathbf{f}_2)$ and $\sigma_p(\mathbf{f})$ are the errors of the bispectrum phasors and the squared visibilities, respectively. $O_k(\mathbf{f})$ denotes the Fourier transform of the iterated image $o_k(\mathbf{x})$. $\mathbf{f}_1, \mathbf{f}_2$, and $\mathbf{f}$ are two-dimensional spatial frequency vectors. The functions $w_d(\mathbf{f}_1, \mathbf{f}_2)$ and $w_d^p(\mathbf{f})$ are the uv point density weights of the bispectrum and the squared visibilities $|O(\mathbf{f})|^2$, respectively. The factor $\gamma_0$ scales the bispectrum phasor and the squared visibilities of the iterated image to minimize the value of $Q_3[o_k(\mathbf{x})]$. The variable $f_0$ is the relative weight between the CP term and the squared visibility term in $Q_3[o_k(\mathbf{x})]$.

## Appendix E: Estimation of the error maps of the reconstructed images

In this section, we describe the building of the error maps of the reconstructed L- and N-band images shown in Fig. 1. For this goal, we performed computer simulations. In order to cover the noise introduced by the gaps in the uv coverage and by the noise in the measured data, we a) used the same uv coverage as the one of the MATISSE observations, and b) degraded the computed visibilities and CPs by simulated noise corresponding to the noise in the MATISSE observations, respectively. As computer simulation targets, we have chosen the models similar to those shown in Figs. 6 and 7, but with the central star intensity
2018-12-09T08_16_26-07_53_28 Documents average

<figure><img src="image_21.png" /><figcaption>**Fig. B.5.** Example of calibrated visibilities and CPs of one *N*-band observation without and with spectral binning: red line denotes data without spectral binning, and the blue, green and brown lines correspond to data with spectral binnings of 9, 11 and 17, respectively. Data processing without spectral binning (i.e., binning 1) is just shown for illustration, but is usually not used because of the severe noise (see Sect. 2.2 for details).</figcaption></figure>
2018-12-09T08_16_26-07_53_28averaged

<figure><img src="image_22.png" /><figcaption>**Fig. C.1.** Comparison of measured visibilities and CPs in *L* band with those derived from the reconstructed image: (red) measured data, and (blue) data derived from the reconstruction shown in Fig. 1 top. The data are averaged over all six data sets of each observation.</figcaption></figure>
2018-12-09T08_16_26-07_53_28averaged

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

**Fig. C.2.** Comparison of various types of observed and model visibilities and CPs in the *L* band: (red) measured data, (blue) data derived from the reconstruction shown in Fig 1 top, (green) data derived from the best-fit model 3 (Fig. 6, but with the central star scaled by factor 9.4), and (light brown) data from the best-fit model 4 (Fig. 9, but by scaling of the central star with factor 3.9). For more details see Sects. 4.6 and 4.7. The data are averaged over all six data sets of each observation.
2018-12-09T08.16.26-07.53.28averaged

<figure><img src="image_24.png" /><figcaption>**Fig. C.3.** Comparison of measured visibilities and CPs in *N* band (spectral binning 11) with those derived from the reconstructed image: (red) measured data, and (blue) data derived from the reconstruction shown in Fig. 1 bottom. The data are averaged over all four BCD configurations of each observation.</figcaption></figure>
2018-12-09T08_16_26-07_53_28averaged

<figure><img src="image_25.png" /><figcaption>**Fig. C.4.** Comparison of various types of observed and model visibilities and CPs in the *N* band (spectral binning 11): (red) measured data, (blue) data derived from the reconstruction shown in Fig. 1 bottom, and (green) data derived from the best-fit *N*-band model 3 (Fig. 7, but with the central star scaled by factor 9.2). The data are averaged over all four BCD configurations of each observation.</figcaption></figure>
increased to the central star intensity of the reconstructed images in Fig. 1. The image reconstructions were performed with IRBis and the same parameters were used as for the reconstructed images shown in Fig. 1 and listed in Table 1. The error map of each of the reconstructed images (Fig. 1) is the absolute value of the difference between the computer simulation reconstruction and the computer simulation target. Before calculating the difference, the reconstruction was a) scaled to the total intensity and b) shifted to the center of the computer simulation target. To assign the reconstructions in Fig. 1 to the above produced error maps, each of the two reconstructions was a) normalized to the same total intensity as the computer simulation target, and b) shifted to the center of the computer simulation target. Figures E.1 and E.2 present the computer simulation experiments for deriving the error maps of the L- and N-band reconstructions in Fig. 1, respectively. In all four images shown in Fig. E.1 (and E.2), the same intensity is coded with the same color. Figures E.1a and E.2a are the L- and N-band computer simulation targets, which are convolved with point spread functions corresponding to sizes of 2.7 mas and 6.6 mas, respectively. Figures E.1b and E.2b present the reconstructed images obtained from the simulated data. Figures E.1c and E.2c show the error maps for the L- and N-band reconstructions (Fig. 1), derived from the computer simulation targets (Figs. E.1a and E.2a) and the computer simulation reconstructions (Figs. E.1b and E.2b). Finally, Figures E.1d and E.2d are the reconstructions shown in Fig. 1.

## Appendix F: Comparison of the measured CPs and visibilities with the ones derived from the model intensity distribution using the $\chi^2$ method

The model-observation $\chi^2$ values (reduced $\chi^2$) of the CPs and visibilities of the L-band model image displayed in Fig. 6 are 738 and 1163, respectively. These big $\chi^2$ values are mainly due to the fact that the model flux ratio of the central object and the disk disagree with the observations as discussed in Sect. 4.5. These model-observation $\chi^2$ values are also very big in comparison to the reconstruction-observation $\chi^2$ values between the visibilities and CPs of the reconstructed L-band image (Fig. 1 top) and the observed visibilities and CPs. The CP $\chi^2$ value is 2.88 and the visibility $\chi^2$ value is 3.97. The reconstruction-observation $\chi^2$ values of the reconstructed N-band image (Fig. 1 bottom) are 1.06 for the CPs and 4.94 for the visibilities.

Figure F.1 shows the CP and visibility $\chi^2$ values of model 3 (Fig. 6) as a function of different scaling factors of the model central star. The minimum CP $\chi^2$ value of 178 is achieved with a scaling factor of 9.0, which is similar to 9.4 obtained with the residual method (see Fig. 10).

Figure F.2 shows $\chi^2$ values of the model 4 (Fig. 9) as a function of the scaling factors. In this case the minimum CP $\chi^2$ value of 369 is obtained with a scaling factor of 3.7, which is similar to 3.9 derived with the residual method (see Fig. 11).

Figure F.3 presents the CP and visibility $\chi^2$ values of the N-band model image displayed in Fig. 7 as a function of the different scaling factors of the model central star. The minimum CP $\chi^2$ value of 10.6 is obtained with a scaling factor of 27, whereas the CP $\chi^2$ value in the case of no scaling is 11.0. The difference between scaling and no scaling of the central model star is very small because of its very low flux (i.e., 0.03% of the total flux of the best-fit N-band model). The very low flux of the central star is responsible also for the different optimal scaling factors obtained with CP residuals (optimal scaling factor is equal to 9.2) and CP $\chi^2$ values (optimal scaling factor is 27 in this case).

<figure><img src="image_26.png" /><figcaption>**Fig. E.1.** Estimation of the error map of the L-band reconstruction (Fig. 1 top) from synthetic data: a) theoretical object used for generating the synthetic data (CPs and visibilities with the same noise and uv coverage as in the measured data (see text for more details), b) image reconstructed from the synthetic data, c) derived error map, d) reconstruction shown in Fig. 1 (top). The FOV is 64 mas. North is up and east to the left.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. E.2.** Estimation of the error map of the *N*-band reconstruction (Fig. 1 bottom) from synthetic data: a) theoretical object, b) image reconstructed from the synthetic data, c) derived error map, d) reconstruction shown in Fig. 1 (bottom). The FOV is 64 mas. North is up and east to the left.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. F.1.** CP and visibility $\chi^2$ values of the best-fit *L*-band model 3 shown in Fig. 6 as a function of the central star scaling factor.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. F.2.** CP and visibility $\chi^2$ values of the modified best-fit *L*-band model 4 shown in Fig. 9 as a function of the central star scaling factor.</figcaption></figure>

Figure F.4 shows the model-observation CP $\chi^2$ values for several versions of model 3 (i.e., $R_\epsilon-\epsilon$ pairs with $R_\epsilon = 1.8, 2.0, 2.2$ and $\epsilon = 0.5, 0.6, 0.7$) as a function of the central star scaling factor. The best-fit models, with the smallest CP $\chi^2$ values of 110 and 111, have a central star scaling factor of 7.8, which is similar to 8.4 obtained with the residual method (see Fig. 13), and $\epsilon = 0.6$ and $R_\epsilon = 1.8$ and 2.2, which is identical to the results of the residual method used in Sect. 4.6. For this central star scaling factor of 7.8, Fig. F.5 shows the same type of $R_\epsilon-\epsilon$ plot as in Fig. 14, but using the model-observation CP $\chi^2$ values instead of the CP residuals. The results of the two fits provide, in both cases, best-fit models with $\epsilon = 0.6$ and $R_\epsilon = 1.8$ and 2.2.
<figure><img src="image_30.png" /></figure>

**Fig. F.3.** CP and visibility $\chi^2$ values of the modified best-fit *N*-band model shown in Fig. 7 as a function of the central star scaling factor.

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

**Fig. F.4.** CP $\chi^2$ values of *L*-band models (different versions of model 3) with different inner rim parameters $\epsilon$ and $R_\epsilon$ as function of the central star scaling factor. The best-fit models are those two with $\epsilon = 0.6$ and $R_\epsilon = 1.8$ and $2.2 R_{\text{in}}$ at the central star scaling factor 7.8.

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

**Fig. F.5.** Model-observation CP $\chi^2$ values of *L*-band models (different versions of model 3) for different inner rim parameters $\epsilon$ and $R_\epsilon$ at the best-fit scaling factor 7.8 of the central star.