© The Authors 2025

1 Introduction

Dust and gas inside protoplanetary disks are the basic material, which will be reprocessed at later evolutionary stages, to form terrestrial and gaseous planets. One of the intermediary stages in star formation is the Class II young stellar objects (YSOs) (e.g., Lada 1987; Williams & Cieza 2011). These are variable in the optical and near-infrared (NIR) due to eruptive events (outbursts) followed by changes in the rates at which they accrete mass from their protoplanetary disks. Typical mass-accretion rates range from 10⁻⁸ to 10⁻⁴ M_⊙ yr⁻¹. These stars can be classified in different groups based on the strength and duration of said variability (see, e.g., Fischer et al. 2023 for a recent review on the subject matter). Prominent sub-groups are the FUor and the EXor YSOs, so named after their archetypes, FU Ori and EX Lup, respectively. Stars from these sub-groups are typically on the lower-mass end (≤2 M_⊙). Their eruptions can be in the range of 1–6 mag in the optical, and such outbursts last for several days to years (EXors), or even decades (FUors). On the other hand, the intermediate-mass YSOs (2–12 M_⊙) also show variability in the form of (out)bursts, although these are typically of lower amplitude; that is, ≤2 mag in the optical. The class of Herbig Ae/Be eruptive stars belongs in this category (Herbig 1960). In all of the above-mentioned cases, the eruptions are irregular.

There are currently two predominant theories as to what causes eruptions. On the one hand, they may result from physical changes – instabilities and/or fragmentations – within the protoplanetary disk itself (Vorobyov & Basu 2010; Vorobyov et al. 2024; Nayakshin et al. 2024). On the other hand, the orbital motion of a companion or a flyby may perturb the disk (Bonnell & Bastien 1992; Skliarevskii & Vorobyov 2023). Regardless of the cause, the eruptions occur at the innermost regions of the protoplanetary disks within a distance ≤1 au from the protostar. In particular, these protoplanetary disks are thought to have a defined structure: (1) a hot and/or warm inner part emitting in the infrared within ≤10 au of the protostar, and (2) a cooler exterior emitting at far-IR and submillimeter wavelengths (e.g., Dullemond & Monnier 2010, and references therein). Their composition is a mixture of gas, dust, and ices. During an eruption, infalling gas within a few stellar radii becomes optically thick and is heated, which makes it visible in the optical and NIR. This is a major difference between FUor disks and those of quiescent protostars, since the latter have an optically thin gaseous disk and their circumstellar emission is dominated by their dusty component. Herbig Ae/Be star protoplanetary disks are overall more massive and much larger in size compared to typical FUor and T Tauri disks. They are divided in two main groups depending on whether they are flared (Group I) or flat (Group II), while they can also have several other characteristics such as line emission, disk gaps, NIR and/or mid-infrared (MIR) excess, and processed dust material such as crystalline silicates or PAHs (Meeus et al. 2001; Maaskant et al. 2013; Menu et al. 2015). Such small spatial scales of a few astronomical units can only be resolved by high-angular-resolution instruments, such as the Atacama Large Millimeter/submillimeter Array (ALMA) for radii of ≥10 au, and the Very Large Telescope Interferometer (VLTI) for shorter radii closer to the protostars. Multiwavelength observations are most suited to painting a clear picture of the processes in these disks.

Z CMa is a member of the Canis Majoris R1 star forming region (e.g., Shevchenko et al. 1999) with a probable age of the order of 5 Myr (Gregorio-Hetem et al. 2009). Hereafter, we adopt the distance estimate of Dong et al. (2022) (1125 ± 30 pc), which derived a median distance based on Gaia data of nearby CMa R1 cluster members. Z CMa is a binary with a separation of ~0.1″ at a position angle of PA ~130° (Koresko et al. 1991; Barth et al. 1994). The northwest component (Z CMa NW) is a Herbig star (HBe star), while the southeastern companion (Z CMa SE) is an FUor. Constraining their individual masses has been a matter of debate. van den Ancker et al. (2004) estimated masses at ~16 M_⊙ for the HBe star and ~3 M_⊙ for the FUor. Both values are slightly higher than the assumed range for low- and intermediate-mass stars. On the other hand, Fairlamb et al. (2015) placed the mass of the HBe slightly lower at 11 ± 1.7 M_⊙. The system is highly variable in the optical with several shortterm (~100 days) bursts and long-term (~1000 days) outbursts that raise its brightness by δV ~ 3 mag (Fig. 1). This EXor-like variability is attributed to the Herbig star, while it is possible that the FUor companion has returned to quiescence (Sicilia-Aguilar et al. 2020 and references therein).

Near-infrared scattered-light imaging at large spatial scales (i.e., >1″) revealed a streamer extending south-southwest from the binary system (Canovas et al. 2012, 2015; Liu et al. 2016). Dong et al. (2022) argued that a potential flyby intruder is located at the edge of this streamer; however, Zurlo et al. (2024) showed that the streamer is in fact part of a large-scale outflow. Further adding to the system’s complexity, jets emanate from both stars. Whelan et al. (2010) estimated the jet position angles at ~235° and ~245° for the FUor and HBe jet, respectively, although the uncertainty in the measurements is somewhat large due to the clumpiness seen at sub-arcsecond scales. This clumpiness is evident in both jets, and it could be explained by jet wiggles potentially instigated by additional companions in narrow orbits around the known stars (Whelan et al. 2010; Antoniucci et al. 2016). On the other hand, a disruption induced by a stellar flyby could also explain the jet wiggles.

Near-infrared interferometric observations revealed two very compact sources with a diameter of ≤5 au, and these correspond to the hot, innermost disk regions (Monnier et al. 2005; Millan-Gabet et al. 2006; Benisty et al. 2010). Observations of the system in the MIR have been limited by the spatial resolution of single-dish telescopes since the binary is seen as an oblate source (Malbet et al. 1993; Persi et al. 2007), although the Herbig disk is the brightest of the two (Monnier et al. 2009). On the other hand, the binary was resolved in the MIR with the MID-infrared Interferometric instrument (MIDI) on the VLTI, but the interpretation of the results was model-dependent (Ratzka 2005; Varga et al. 2018).

Protoplanetary disks often exhibit substructures such as spiral arms, arcs, rings, and gaps, to name but a few. It is widely believed that such structures are (in)direct evidence of disk–planet interactions (Benisty et al. 2023). For the case of Herbig disks, it is believed that the formation of gaps might be related to their evolutionary stage (Maaskant et al. 2013; Menu et al. 2015). The substructures have been observed at different spatial scales and wavelengths, and comparisons have shown that in several cases, the cooler outer disks observed in the (sub)millimeter continuum can be misaligned with structures observed in the NIR (e.g., Andrews et al. 2018; Bohn et al. 2022). Such misalignments might be more profound in Herbig disks (e.g., Kraus et al. 2020), but this could be an observational bias since these disks are usually larger in size.

We revisited Z CMa in the MIR with the Multi-AperTure mid-Infrared SpectroScopic Experiment (MATISSE; Lopez et al. 2022), the new imaging instrument of the VLTI. Our goal was to image both disks and look for substructures within each, as well as for potential misalignments with respect to structures observed in the (sub)millimeter and NIR regimes. Furthermore, we explored the effect EXor-like (out)bursts can have on the inner disk region of the Herbig star. We also measured the sizes of the dusty, inner regions of each disk and, when possible, explored their mineralogical composition. This work focuses on the imaging results and is the first in a series of papers on the MATISSE data, which will be followed by a spatio-kinematic analysis of the gaseous inner regions and a study on the system’s variability.

The observations are described in Section 2, and the results of our analysis are shown in Section 3. Our interpretation is discussed in the penultimate section (Sect. 4), which is followed by our concluding remarks. Supplementary material can be found in the appendices.

Fig. 1 V-band light curve of Z CMa with corresponding epochs of MATISSE observations (orange lines; dashed and solid for the quiescent and the outburst phase, respectively) as in Table B.1. Photometric data from the Archive of the American Association of Variable Star Observers (AAVSO; green circles), the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL) Optical Monitoring Camera (OMC; magenta circles), and the Kamogata/Kiso/Kyoto Wide-field Survey of Variable Star Observers League in Japan (KWS; blue circles Maehara 2014). The imaging data were obtained during the outburst from late 2022 to early 2023. Also shown for reference are the epochs of archival ALMA and JVLA observations (purple lines; Takami et al. 2019).

2 Observations

2.1 Archival photometry and the variability of Z CMa

A current visual and V-band light curve of Z CMa is shown in Fig. 1. The Herbig star experienced yet another outburst in late 2022 that lasted longer than 100 days. It appears to have remained in outburst until early December 2023, with a probable short-term dimming in October, although we note that there is a broad gap in temporal coverage from May to September 2023 when Z CMa is not observable. As described below, this last outburst serendipitously occurred while our observations were obtained (orange lines). Earlier submillimeter and radio interferometric observations by ALMA and JVLA (Takami et al. 2019) are also indicated for reference (purple lines; Fig. 1).

Figure 2 shows the spectral energy distribution (SED) of Z CMa during the outburst. The SED has been dereddened for A_V = 2.5 (Appendix A), and it encompasses emission from both sources in Z CMa; however, the Herbig star and its disk would be the main contributors during the outburst. The optical photometry (BVRI; black points) represents the median values from the Kamogata/Kiso/Kyoto Wide-field Survey of Variable Star Observers League in Japan (KWS, Maehara 2014) during the outburst. Our SpeX/IRTF and MATISSE spectra (see below) are shown in red and blue, respectively. Also included in black are archival data from Galaxy Evolution Explorer (GALEX, Beitia-Antero & Gómez de Castro 2016), Midsource Space Experiment (MSX, Egan et al. 2003), AKARI (Ishihara et al. 2010), Infrared Astronomical Satellite (IRAS, Helou & Walker 1988), James Clerk Maxwell Telescope (JCMT, Weintraub et al. 1991; Dent et al. 1998; Di Francesco et al. 2008), Swedish-ESO Submillimetre Telescope (SEST, Reipurth et al. 1993), and Institute for Radio Astronomy in Millimeter Range 30-m telescope (IRAM30, Altenhoff et al. 1994). A Kurucz model (8500 K; green) and a blackbody (100 K; orange) have been fit for reference. The far-infrared to (sub)millimeter tail of the system’s SED has a spectral index of α =3.2 ± 0.2.

Fig. 2 Dereddened SED of Z CMa during the 2023 outburst. Our own SpeX/IRTF and MATISSE spectra are shown in red and blue, respectively. For a description of the archival photometry (black points), we refer the reader to Sect. 2.1. A Kurucz model (green) and a blackbody model (orange) are shown for reference.

2.2 MATISSE/VLTI

We obtained high-angular-resolution interferometric data with MATISSE/VLTI in the L (2.9–4.1 μm), M (4.5–4.9 μm), and N (8–13 μm) bands in “imaging” mode from 29 November, 2022 until 25 March, 2023 (ESO Program ID: 0110.C-4209(A), PI: F. Lykou), while Z CMa was in outburst (Fig. 1). In the imaging mode, all three array configurations (“small”, “medium”, and “large”) of the 1.8-m Auxiliary Telescopes (ATs) were utilized in order to obtain the best achievable uv-plane coverage (Table B.1), and this is shown in Fig. 3. Each hour of observation includes data acquisition for both science target and a calibrator in all three bands.

The MATISSE instrumental field of view (FOV) is that of the beam entering the instrument through a pinhole, and for the case of an AT with a given diameter, D (1.8 m), it is 1.5λ/D in the LM band, and 2λ/D in the N band (Lopez et al. 2022). Therefore, the FOV at 3.5 μm and 10.5 μm becomes 0.6″ and 1.8″, respectively. As such, both stars in Z CMa were observed simultaneously by MATISSE and it was not possible to define individual pointings.

The “GRAVITY for MATISSE” (GRA4MAT) mode (Woillez et al. 2024) was used for all observations, where the GRAVITY instrument – the VLTI instrument operating in the K band (2–2.4 μm) – is utilized as a fringe tracker to allow for better stabilization of the fringes of MATISSE, improving signal-to-noise (S/N) and calibration. Several spectral setups were chosen for these observations. For the imaging runs, we opted for the medium spectral resolution in the L band (MED-L; R ~ 500) and the low spectral resolution in the N band (LOW-N; R ~ 30).

Data reduction and analysis were performed with the MATISSE DRS pipeline version 1.7.5 and with the dedicated Python library mat_tools¹ for flux calibration and visualization. The basic concepts of the MATISSE data reduction have been introduced extensively in several other works (e.g., Varga et al. 2021; Gámez Rosas et al. 2022); therefore, we refrain from repeating them here. Only one calibrator was used for both spectral bands for the imaging observations. Those were α CMa (A1V, diameter 6.09 mas) for the small array configuration and HD 48217 (K5III, diameter 2.60 mas) for the medium and large arrays. An assessment of the data quality and caveats is given in Appendix B.

Our data are complemented with guaranteed time observations (GTO) obtained on 23 January, 2021 and 15 November, 2021 during quiescence (Fig. 1) with the medium array configuration (ESO Program ID: 0106.C-0501(B) and 0108.C-0385(D); PI: B. Lopez). The GTO runs (Table B.1) contained experimental observations in the L band with both medium and low spectral resolutions (LOW-L: R ~ 30), and both the low and the high spectral resolution in the N band (HIGH-N: R ~ 220). The first GTO data were obtained in a CAL_N-SCI-CAL_L sequence, whereby a different calibrator was used for each spectral band, namely α CMa (A1V, diameter 6.09 mas) as CAL_N and HD 58972 as CAL_L (K3III, diameter 3.39 mas). It is worth noting that since MATISSE observations were obtained simultaneously in both bands, the science target’s data can be reduced with both calibrators to check consistency. The second GTO observations were obtained with one calibrator for both bands (HD 39853; K4III, diameter 2.23 mas).

Fig. 3 MATISSE uv-plane coverage of “imaging” mode observations of Z CMa during its outburst phase. The uv coordinates are colored per respective VLTI configuration array as in Table B.1; that is, the small array included the A0-B2-D0-C1 quadruple, the medium array included the K0-G2-D0-J3 quadruple, and the large included the A0-G1-J2-K0 quadruple, which was interchanged with the A0-G1-J2-J3 to maximize uv coverage.

2.3 SpeX/IRTF

We obtained long-slit spectra with the SpeX spectrograph (Rayner et al. 2003) on the Infrared Telescope Facility (IRTF) on Mauna Kea, Hawaii, on 04 November, 2023 (program ID: 2023B039, PI: P. Ábrahám) while the system was still in outburst. We opted for the Echelle grating in short-wavelength mode (SXD), which covers the NIR wavelength region (0.7–2.55 μm). The achieved spectral resolution was R ~ 750. The slit width of 0.8″ covers both stars, and the spatial resolution is insufficient to recover individual stellar spectra.

HIP 24555 (A0V) was chosen as the telluric standard. Although we did not obtain simultaneous photometry, we found that the telluric standard was sufficient for photometric calibration. Data reduction, corrections for telluric absorption, and flux calibration were performed using the SPEXTOOL package (Vacca et al. 2003; Cushing et al. 2004).

The SpeX spectrum is rich in emission lines (red; Fig. 2), whose analysis will be the basis of a forthcoming paper. Here, we show that the photometric calibration of MATISSE agrees with SpeX, and we use both to build Z CMa’s SED in Fig. 2.

3 Results

The interferometric data – that is, the correlated fluxes (Figs. B.1–B.3), visibilities (Figs. B.4–B.6), and closure phases (Figs. B.7-B.9) – show a sinusoidal modulation in several baselines in all three bands, indicating that the binary system is resolved. The amplitude of the sinusoidal modulation appears shallower in the M-band correlated fluxes and visibilities than in the L and N bands. In the following section, we provide initial estimates of the binary’s separation and position angle through geometric model fitting of the visibilities and closure phases during quiescence. Image reconstruction results are shown in Sect. 3.2. The N-band correlated fluxes obtained with the small array show a deep absorption feature which we explore further in Sect. 3.3.

We note the presence of several emission lines in the L- and M-band correlated fluxes from the small array observations (Figs. B.1 and B.2), which, as mentioned in the caveats in Appendix B, belong to the hydrogen lines of the calibrator Sirius. Although these could be subtracted by averaging the emission around the line, we refrained from doing so here to indicate the adverse effect of this calibrator.

3.1 Geometric model fitting

Although this work focuses on the imaging results with MATISSE, we performed initial geometric model fits to the visibilities and closure phases of the GTO LOW-L and HIGH-N data obtained during quiescence (id: 1, Table B.1) to derive approximate size estimates for the two sources, as well as relative astrometry, and compared these to literature values and the imaging results.

We employed LITPro (Tallon-Bosc et al. 2008), a geometric modeling tool that does not take into account the chromaticity of the data. For simplicity, we followed literature studies and fit centrosymmetric distributions (Ratzka 2005; Monnier et al. 2005; Millan-Gabet et al. 2006; Benisty et al. 2010) such as a uniform disk, a Gaussian distribution, and a binary in the form of two point sources, where the size of each source, the flux ratio, and the position of the companion are free parameters, while the Herbig star’s position is fixed at the origin (0,0). The results of the LITPro fitting process are shown in Table 1. These include the angular sizes of each disk (θ_HBe and θ_FUor), their relative flux ratio Δf (HBe to FUor), the relative position of the FUor companion with respect to the Herbig star (Δα, Δδ), and the goodness-of-fit ${\chi }_{\text{red\hspace{0.17em}}}^2$.

Many of the fits are poor as the flux contributions and the size of the FUor were often highly correlated. This is more prominent when fitting two point sources. Since LITPro fits the entire spectral range of the MATISSE L and N bands, and does not take into account either the chromaticity of each source or the probable temperature gradient in each disk, comparisons to the results from Sect. 4 and Table 3 could be tenuous. Attempts to improve such fitting in LITPro by introducing background emission and/or incorporating a blackbody temperature for each disk failed to produce reliable results since the flux-ratio parameters were highly correlated and/or degenerate. Nevertheless, the geometric model fitting results can be interpreted as showing that the size of the HBe disk was larger than that of the FUor, irrespective of the geometric model during quiescence in 2021. For consistency, we derived the average relative astrometry from the three methods in the L band, since the spatial distribution of the dust continuum is expected to affect the estimates in the N band. As such, from Table 1 we calculate ⟨Δα⟩ = 77.17 ± 0.03 mas and ⟨Δδ⟩ = −81.96 ± 0.09 mas for the L band (see Table E.1). This result is included in our analysis in Sect. 4.3.

From the above, we argue that the simple geometric fitting is convenient for obtaining initial parametric results on the source sizes and good estimates of the relative astrometry; however, it is not always reliable since it does not take into account the brightness distribution and intricate geometries (e.g., asymmetries) of each source. Since this work is focused on imaging data with MATISSE, we expect the reconstructed images to shed more light on the spatial distribution of both disks.

3.2 Image reconstruction

The image reconstruction concentrates on the MATISSE data obtained during the outburst phase (Fig. 1). As mentioned in Sect. 2, the MATISSE FOV is 0.6″ and 1.8″ at 3.5 μm and 10.5 μm, respectively. For consistency, we opted to use a common FOV and pixel scale for all reconstruction methods. Using the full FOV of MATISSE did not improve any of the initial image reconstruction tests, as it only increases the amount of residual noise in the maps while increasing computational time. Therefore, in all further image reconstruction the FOV was 400 × 400 mas with a pixel scale of 1 mas/pixel in the LM band and 3 mas/pixel in the N band. These pixel scales are used to achieve “super-resolution” (i.e., Nyquist sampling for each image), as that is estimated at the longest baseline obtained.

Pseudo-continuum images were reconstructed in all bands, that is, for short bandwidths where no spectral lines are observed. These were at 3.4–3.6 and 4.6–4.8 μm in the L and M bands, and 8.4–8.6, 9.1–9.3, 10.3–10.5, and 11.4–11.6 μm in the N band. Henceforth, the pseudo-continuum images are mentioned as the 3.5, 4.7, 8.5, 9.2, 10.4, and 11.5 μm images, respectively. All image reconstructions are smoothed with a Gaussian of similar FWHM to the image resolution, that is, 1 mas in the LM and 3 mas in the N bands.

We used the online tool OImaging² to test two image-reconstruction methods (BSMEM and MiRA) for consistency and finally opted for the MiRA algorithm (see Appendix C). All array configurations with ids=3 to 23 and data quality ranked as A and B were used (see Table B.1). The reconstructed images are shown in Fig. 4 and are ordered by increasing wavelength from left to right and top to bottom. All images are normalized to unity and shown in a logarithmic stretch with a minimum value of 10⁻². This level is similar to the arbitrary flux-level minimum set during image reconstruction³. Any structures seen in the maps below that limit are residual artifacts due to incomplete uv-plane coverage (see Appendix C; Planquart et al. 2024 offers an extensive treatise on MATISSE image reconstruction). The theoretical MATISSE beam at each wavelength is shown in the bottom left corner of each panel.

The primary HBe is the obvious brighter and larger source in the middle of each panel in Fig. 4. The FUor companion is located as expected toward the southeast. All pseudo-continuum reconstructed images in Fig. 4 are shown at the same FOV to reveal the changes in the angular size of the sources. This is more prominent for the HBe primary. On the other hand, the FUor’s size remained nearly constant at all wavelengths.

The separation of the Z CMa binary was calculated from the reconstructed images in Fig. 4. Each position was estimated by finding the location of the brightest peaks in each map (in all bands) with the DAOFIND algorithm of the photutils package⁴ and then averaging over all maps since the observations were obtained simultaneously in all bands. Here, we did not take into account the temporal spread of the observations over four months (Table B.1), since this is a wide binary with a rather large period. The positional errors are the standard deviations of the measured locations. Our MATISSE measurements result in a separation 117.88 ± 0.73 mas and at a position angle of 139.16° ± 0.29° east of north. This epoch will be included in Sect. 4.3.

Apart from the relative astrometry, the DAOFIND algorithm also calculated fluxes for each star, that is, the central bright region in each image. Since the pseudo-continuum reconstructed images are not calibrated with respect to absolute photometry, and considering that this is a nontrivial matter in image reconstruction, these fluxes are only relative. As such, we estimated that the relative-flux ratios Δf at the nearest decimal of the HBe with respect to the FUor, and these can be found in Table 2. Here, the uncertainties reflect the range of values for apertures of 1–3 pixels wide.

It appears that in the N band the FUor companion is more than ten times fainter than the HBe primary. As a reminder, the MATISSE imaging data were obtained while the HBe star was in outburst. The brightness ratio between the FUor and HBe star has been studied before. In the infrared and up to the L′ band, Hinkley et al. (2013) and Bonnefoy et al. (2017) have shown that the flux ratio of the HBe with respect to the FUor increases at longer wavelengths, especially while the primary is in outburst. It is worth mentioning that the binary cannot be resolved with conventional AO-assisted imaging beyond the L band. Ratzka (2005) estimated an HBe/FUor flux ratio ≈${6.3}_{-1.5}^{+2.8}$ in the N band by fitting the visibilities of MIDI data obtained while the system experienced a short-term burst in late 2004. In that respect, our results are consistent with literature studies of the (out)burst phase. Moreover, Hinkley et al. (2013) reported a continuous decay in the flux of the FUor companion with respect to historical values. Considering that their observations were obtained approximately 15 years ago (at the time of writing), it is worth revisiting the system to explore whether the FUor continues its dimming and if it has returned to quiescence.

Figure 5 shows a comparison of the protoplanetary disk sizes between MATISSE and ALMA (Dong et al. 2022). The underlying reconstructed image (in blue) corresponds to the 3.5 μm pseudo-continuum, while the overlaid gray contours correspond to the 9.5 μm image reconstruction. Both have been normalized to the maximum peak brightness (here, the Herbig star), and the contour levels are at 5%, 10%, 50%, and 90% of the peak. This comparison shows that the majority of the MIR flux arises from the same region at 3.5 and 9.5 μm for the Herbig disk (size≤15 mas), while the disk may extend up to ~40 mas in size (upper limit) at lower flux levels. The FUor disk’s size remains somewhat constant in both wavelengths (size ≪15 mas), although it is much fainter compared to the HBe at 9.5 μm (5% contour level). Also overlaid in Fig. 5 are two ellipsoids (blue dashed lines), which correspond to the maximum deconvolved disk sizes (Gaussian FWHM; see Table 3) from the ALMA 1.3-mm continuum in Dong et al. (2022). For reference, the ALMA beam was 78 × 46 mas and PA 66° and the size of both sources prior to deconvolution were similar to the beam size, while Dong et al. (2022) caution that “extended dust emission may affect their fitting process”. The latter could explain the discrepancy in disk inclination between 1.3-mm and 9-mm data (see Table 3). The optically thin and cooler outer regions of protoplanetary disks that emit in the (sub)millimeter regime are expected to be significantly larger than the optically thick inner regions that are emitting in the NIR and MIR (see, e.g., Dullemond & Monnier 2010).

Fig. 4 MiRA image reconstructions in L, M, and N bands during the 2023 outburst. From left to right and top to bottom, these correspond to pseudo-continuum reconstructed images at 3.5, 4.7, 8.5, 9.2, 10.4, and 11.5 μm. These are smoothed with a Gaussian beam with a FWHM similar to the image resolution, normalized to unity, and shown in a logarithmic stretch with a minimum value of 10−2 to suppress image reconstruction residuals below that level. Such residuals are the result of insufficient uv coverage (Fig. 3). The SE component (FUor) becomes much fainter than the NW component (HBe) at longer wavelengths. The beam size is shown as a black circle at the bottom left corner of each panel.

Fig. 5 Comparison between ALMA and MATISSE. The normalized pseudo-continuum reconstructed image at 3.5 μm is shown in blue. Overlaid are contours (gray lines) at 5, 10, 50, and 90% of the peak of the normalized 9.5 μm reconstructed image. These are compared to the ALMA 1.3-mm continuum disk sizes (blue dashed-line ellipsoids) of Dong et al. (2022).

3.3 Radial variation of the silicate feature and optical depth analysis

The MATISSE N-band correlated fluxes⁵ are shown in Fig. B.3. Each panel corresponds to a separate array configuration to those in Table B.1 for the data obtained during the 2022/3 outburst. As a reminder, these are the so-called imaging data not to be confused with the GTO data obtained during quiescence.

For the medium and the two large arrays (Table B.1) in Fig. B.3, the N-band correlated spectra are dominated by three factors: (a) the sinusoidal modulation of the binary, (b) the effects of the terrestrial atmosphere (atmospheric seeing, wind speed, and coherence time), and (c) the rather noisy, low S/N for certain baselines. For the second case, the effect is increased for fainter calibrators such as the one used in these arrays. For the third case, these baselines are usually the combination of AT2-AT3 telescopes, that is, the shortest baselines for each array, and could be related to the known, occasionally poor performance of AT3. Considering these uncertainties, we excluded the correlated spectra obtained with these arrays from our analysis of the silicate feature.

The small-array correlated spectra (Fig. B.3) show the silicate feature in absorption. The strength (flux level) of the feature is decreasing with increasing baseline length, that is, from 10.4 m to 33.9 m (or else at spatial scales from approximately 32 to 104 mas). We note that despite our concerns about using a large calibrator, such as Sirius (see Appendix B), the overall shape and strength of the feature are consistent with MIDI/VLTI observations obtained with different calibrators but with similar baselines (Varga et al. 2018). Therefore, the radial variation of the strength of the silicate absorption feature with regard to baseline length is real.

The binary can be resolved even with the small array data in the N band, although the disks are marginally resolved (size <10 mas) and no other spatial information can be derived for them. Therefore, it is possible that the N-band correlated spectra also contain the sinusoidal modulation of the binary’s signal, which could slightly shift the minimum of the silicate feature. However, this appears to be stationary taking into account its proximity to the 9.4–9.9 μm ozone feature that affects that part of the spectrum.

This radial variation of the correlated spectra may be related to changes in the optical depth of dusty material along the line of sight; that is, the correlated fluxes are attenuated by a factor of ∝ e^−τ_λ (e.g., Lykou et al. 2024). The interstellar extinction should also be taken into account (e.g., Quanz et al. 2007) before considering estimating the circumstellar extinction, and to that end we dereddened each correlated spectrum using A_V = 2.5 (see Appendix A). Following that, we explored this attenuation using the methods of van Boekel et al. (2003) and Kemper et al. (2004), whereby we fit a linear continuum F_c between 8 and 12.9 μm for each correlated spectrum F_corr, and we extracted the optical depth as τ_λ = −ln(F_corr/F_c). The calculated optical depths per epoch and per baseline are shown in Fig. D.1. We note that τ_λ < 0.3 overall. We discuss our analysis further in Sect. 4.2.

4 Discussion

4.1 Protoplanetary disks

Earlier attempts to image the Z CMa system in the MIR with AO-assisted telescopes failed to resolve the binary. Malbet et al. (1993) suggested a circumbinary disk with a size of $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}38\pm 0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}06$ and a PA 161° ± 8° in the L′ (3.74 μm) and M (4.75 μm) bands. On the other hand, Monnier et al. (2009) could not fully resolve the binary in the N band, but they only found an extension toward the FUor location. They measured a source size of (68 ± 2) × (41 ± 6) mas at the presumed location of the Herbig star with a PA at 137° ±5° and an over-resolved (halo) component at (2 ± 1)% of the total flux. The MATISSE images do not reveal any circumbinary warm dust emission.

Table 3 lists the protoplanetary disk sizes corresponding to a 2D Gaussian fit to each disk in the MATISSE image reconstructions for the brightest regions (that is ≥50% of the peak). Based on their apparent separation on the sky at less than 120 mas, the radius of the disks ought to be truncated to less than half the separation; that is, <60 mas, or <70 au, at our adopted distance. Our findings for each individual disk are discussed further below.

Apart from this work, there exist only a handful of interferometric studies of the Z CMa system from the last 20 years. Most of them observed the system while the Herbig star was in quiescence (Monnier et al. 2005; Millan-Gabet et al. 2006; Dong et al. 2022), while two studies observed it during an (out)burst (Monnier et al. 2009; Benisty et al. 2010). The source size estimates from the literature are also listed in Table 3. Also listed are the protoplanetary disks’ orientations where available. For the Dong et al. (2022) deconvolved 1.3 and 9 mm continuum disk sizes, we calculated the inclination as i = arccos (a_minor/a_major), and we note the large uncertainties based on error propagation.

4.1.1 The FUor companion

In the early 1990s, the FUor companion was found to be the brightest of the pair in the visual (Barth et al. 1994). More recent studies suggest that the FUor is slowly dimming both in the visual and in the infrared (Hinkley et al. 2013; Bonnefoy et al. 2017; Sicilia-Aguilar et al. 2020). The diameter of its hot, innermost region was found to be less than 4 mas in the K band (Millan-Gabet et al. 2006), which translates to ≤2.2 au radius at our adopted distance of 1125 pc.

Our analysis of the MATISSE data shows that the MIR emitting regions of the FUor disk is resolved in all bands. At shorter wavelengths, that is, in the L and M bands, the majority of the emission (≥50% peak) rises from a region with a size (i.e., diameter) of ≤4 mas, while at lower levels (1% peak) it extends to <9 mas. In the N band, the brightest disk region extends up to ~15 mas (≥50% peak) and <20 mas at lower levels. When compared to the geometric model fitting from the quiescent epoch data, the angular sizes at lower flux levels agree relatively well (Sect. 3.1). This is not surprising since it is the Herbig star that went into outburst, and the size of the FUor disk should not have been affected.

Overall, the innermost disk region emitting in the MIR has a radius of <5 au, which is consistent with the K-band size mentioned above, confirming that the inner disk is compact. Similar results have been found for other FUor systems (Bourdarot et al. 2023; Lykou et al. 2022, 2024). The cooler, outer parts of the FUor disk as measured by ALMA, even though marginally resolved at 1.3 mm, are indeed wider than those of MATISSE with a radius of ≤30 au.

Our estimates for the disk orientation from the MATISSE image reconstructions have large uncertainties (Table 3). The mean position angle of the disk is 152.5°±12.7° and the mean inclination is 33.5°±5.7°. The uncertainties are the standard deviations of the individual estimates. Given the extremely large uncertainties in the disk PAs and inclinations calculated from the millimeter continuum maps of Dong et al. (2022), we cannot confirm a misalignment between the warm inner disk and its cooler exterior. On the other hand, the disk PA agrees rather well with the orientation of the FUor’s jet; that is, approximately 235° (Whelan et al. 2010; Antoniucci et al. 2016).

4.1.2 The Herbig primary

The HBe star is responsible for the EXor-like (out)bursts of the last 20 years (Fig. 1, and Sicilia-Aguilar et al. 2020). Previous interferometric studies in the NIR (K-band), which were obtained during quiescent and outburst phases, have found that the diameter of its hot, inner disk region is smaller than ≤4 mas (Monnier et al. 2005; Millan-Gabet et al. 2006; Benisty et al. 2010). In fact, the diameter appears to have remained relatively constant irrespective of the variability phase. MATISSE can marginally resolve that inner disk region in the MIR. Mid-infrared observations with MIDI/VLTI could resolve the pair of stars – but no spatial information could be extracted for each disk – and a binary model fit suggested the Herbig star was brighter with a flux ratio of 0.16 ± 0.05 at ~10 μm (Ratzka 2005). However, an additional extended symmetric component with a θ_FWHM ~ 17 mas had to be included in the models to reproduce the observed visibilities. That component could represent the regions probed with MATISSE well. Varga et al. (2018) modeled the N-band MIDI data of the HBe disk with a temperature gradient model (q ~ 0.87) and a disk’s inner rim at the dust sublimation radius estimated at ~3.9 au for a luminosity⁶ of ~ 3200 L_⊙ and a dust sublimation temperature of T_sub = 1500 K. The disk’s half-flux radius was at ~10.3 au. At their adopted distance of 1.05 kpc, the dust sublimation diameter from the temperature gradient model is ~7.5 mas, which is almost double that of the K-band diameter estimates for a similar distance (Monnier et al. 2005; Benisty et al. 2010). The N-band half-flux diameter is ~20 mas.

The MATISSE reconstructed images show that the majority of the emission (≥50% peak) arises from a region with a diameter <5 mas in the L band and <15 mas in the N band. This would translate to disk radii <3 and <9 au, respectively. The region seen in the L band could correspond to dust continuum from the inner rim of the disk. Although the true luminosity of the HBe is not known, we can assume 5000 L_⊙ as a lower limit based on the SED fitting (Fig. 2). The typical sublimation temperature for silicate grains is 1500 K. Therefore, the dust sublimation radius at the given luminosity becomes R_sub ≈ 5 au. This radius can extend as far as 7 au when using the accretion luminosity calculated by Fairlamb et al. (2015), which was a factor of two larger (10 000 L_⊙). It is therefore possible that the L-band emission arises from inside the theoretical dust sublimation temperature. On the other hand, the N-band region where dust emission dominates is more consistent with the latter case when assuming a higher luminosity.

If the discrepancy among the region detected at 2 μm, the region detected in the L band, the dust sublimation radius from the temperature-gradient model, and the above-mentioned theoretical R_sub estimates is real, then one could argue for a potential cavity or gap in the disk and whether that could be detected in the MIR (e.g., Maaskant et al. 2013; Menu et al. 2015). However, the dereddened SED of the 2023 outburst shows NIR excess, even for higher values of A_V, and it is relatively flat with a MIR excess (Fig. 2). These suggest a continuous disk. Unlike in Hinkley et al. (2013), which provided photometry for each star, the SpeX spectra and the averaged archival photometry used in this work originate from the entire system. However, we presume that the optical and NIR fluxes originate from the Herbig star that is in outburst. Overall, the current SED agrees well with the scenario that Z CMa belongs to Group Ib of flared disks (Meeus et al. 2001).

At lower levels, that is, down to 1% of the peak, the diameter of the HBe disk is ≤ 38 mas at 10.4 μm, and ≤48 mas at 11.5 μm. These can be compared to the estimates of Monnier et al. (2009) at 10.7 μm with a source size 68 × 41 mas at the presumed location of the Herbig star. Although that dataset was also obtained at an epoch when the Herbig star was in outburst, it is more likely that the difference in sizes is related to the different observing methods⁷ and not the stellar variability. The abovementioned extent of the disk at those lower flux levels beyond 9 μm would suggest that the warm inner HBe disk is of similar size or slightly larger that the cooler, optically thin dust emission detected at longer wavelengths by ALMA and JVLA (Dong et al. 2022). This discrepancy might be explained if the MIR-emitting disk regions were significantly brighter and therefore appear larger during the 2023 outburst, as opposed to the cooler outer region that was imaged during quiescence (ALMA+JVLA; Fig. 1). However, the MIR spectrum of Z CMa has not changed much overall within the last 20 years (e.g., Polomski et al. 2005; Varga et al. 2018) within the 20% flux uncertainties of MATISSE. Since this size discrepancy would violate the typical model for protoplanetary disks (Dullemond & Monnier 2010), we believe that the millimeter size estimates are either under-estimated if the disks were marginally resolved by ALMA and JVLA, or the millimeter-emitting region is more optically thick than expected. The latter might be the case if the spectral index of the HBe disk is α ≤ 2 (Liu 2019). The discrepancies between the MIR- and millimeter-emitting regions, as well as that of the inner disk rim and the dust sublimation temperature, could be explored with a radiative transfer model for the disk.

Normalized radial profiles were created from each reconstructed map in Fig. 4 within a 60 mas radius of the peake-mission pixel, which we have to presume is the unresolved protostar. These radial profiles are shown in Fig. 6. We set a threshold at 1% of the peak as applied to the image reconstructions. At a first glance, the radial profiles indicate a continuous disk whose brightness diminishes toward larger radii, as would be expected from a temperature gradient. This occurs inside a maximum radius of 50 au, congruent with the expected truncation radius from the binary’s separation. The majority of the MIR emission originates within a 5 mas radius, that is, <6 au (vertical dotted line). This is the location of the first plateau in the 3.5 μm reconstructed map. The inner disk region emitting at shorter wavelengths (3.5 and 4.7 μm) is contained within a radius of <10 mas at its farthest extent (otherwise, ≤11 au). The outer-disk radius can extend up to 36 mas in the N band (i.e., ≤40 au). The N-band profiles show a distinct shoulder at radii of approximately 8–13 mas; that is, at the farthest extent of the inner-most region. Regarding the Herbig star’s variability, it is difficult to ascertain real changes in the disk sizes between quiescence and outburst based on the simplicity of the geometric model fitting (Sect. 3.1, Tables 1).

The MATISSE images suggest an average PA at 156.5°± 12.1° and inclination at 26.6°± 3.8° for the HBe disk (Table 3). As is the case with its FUor companion, this orientation agrees well with the derived PA for the HBe jet (245°; Whelan et al. 2010; Antoniucci et al. 2016). Moreover, the orientations of both disks agree well with each other within the uncertainties. This could suggest that the disks have not been perturbed by a tertiary companion. Once again, we do not find clear evidence for a misalignment between the warm inner HBe disk detected by MATISSE and its cooler outer counterpart viewed by ALMA and JVLA due to the large uncertainties in determining the orientation of the latter.

Based on all of the above, we do not find any clear evidence for substructures such as gaps in the HBe disk, and we do not find any evidence of misalignments either. The current observations might not have had sufficient resolution to reveal such structures, but perhaps this could be revisited in the future with the new 200-m baseline of the VLTI, while the same can be argued for future ALMA observations that could explore spatially resolved kinematics for each disk. Nevertheless, such observations will still be limited by the large distance to Z CMa.

Fig. 6 Normalized radial profiles of HBe disk from reconstructed images in Fig. 4. The horizontal dotted line demarcates the relative flux threshold set for the image reconstructions. The vertical dashed line demarcates the location of the first plateau at 3.5 μm and, consequently, the region from which the majority of the emission originates at all wavelengths.

4.2 Dust distribution

Near-infrared polarimetric imaging studies of Z CMa (Canovas et al. 2012; Hinkley et al. 2013; Canovas et al. 2015; Liu et al. 2016) suggest that the Herbig star is enshrouded in a dusty cocoon. This could potentially be a portion of material spanning hundreds of astronomical units within the overall outflow (e.g., Zurlo et al. 2024). Considering the wavelengths probed, this is due to dust scattering from small-sized grains distributed uniformly over the system, since none of these studies could resolve the binary behind coronagraphic masks. Canovas et al. (2012) modeled the polarization seen in this extended (size >1″) cocoon with a spherical shell around the Herbig star extending to 70 au, with an inner hole at a radius of 20 au, which in turn is embedded in a large-scale, infalling molecular envelope. This morphology could allow for the emanation of a jet, a bipolar outflow, and/or disk winds from the HBe star (Poetzel et al. 1989; Benisty et al. 2010).

At our adopted distance, the cocoon’s dimensions could be translated to angular sizes of ~124 mas for the shell and ~35 mas for the inner hole. Both scales are feasible for MATISSE, and here we explore whether the structures detected in the MIR could be related to the alleged cocoon seen in scattered light through our analysis of the radial variation of the silicate feature (Sect. 3.3).

Despite the limited uv coverage of MATISSE’s small array data (see Fig. 3), we explored the distribution of dust by mapping the optical depth at three different wavelengths, namely 8.5, 10, and 11.5 μm. The optical depth values were extracted at these wavelengths for each baseline or for each uv point (Fig. D.1). We therefore re-imaged the optical depth distribution with respect to each uv point. These uv maps are shown in Fig. 7, and from top to bottom they correspond to 8.5, 10, and 11.5 μm, while the strength of the optical depth is indicated in each individual color bar.

These uv maps show that the dust is not distributed uniformly over the source, but there is a preferential distribution: smaller τ_λ at large spatial scales (i.e., toward the center of each map), and larger τ_λ at small spatial scales (i.e., toward the outer regions of each map). Therefore, this nonuniform distribution is found on spatial scales of approximately 30–100 mas or within an area with a size of 35–110 au. Conversely, the separation of the binary is >100 au, while, as we showed earlier, the FUor companion is >10x fainter than the Herbig star in the MIR.

The nonuniform distribution is more clearly illustrated in Fig. 8, where the uv coverage of Fig. 7 has been converted to angular radii (λ/4B_λ), and these were in turn translated to physical radii for our adopted distance. The optical depth for each wavelength is color-coded with respect to the parameter η, which we define as: $$\eta =|\cos (B_{\mathrm{P}\mathrm{A}}-\theta )|,$$(1)

where B_PA is the baseline PA and θ in the PA of the binary’s orbit estimated from MATISSE data (Sect. 3.2). For baselines parallel to the binary’s orbit, η=1, while for those perpendicular to it, η=0. However, the optical depth does not seem to be dependent on the baseline PA.

As shown in the top panel of Fig. 8, the optical depth at 8.5 μm is relatively flat, with τ_{8.5 μm} <0.1. At longer wavelengths (middle and bottom panels), the optical depth varies with respective physical radius. That is, the highest optical depth resides at radii of <40 au. This is directly comparable to the extent of the HBe disk at 10.4 μm (gray dotted lines) located at approximately 7 and 35 au (upper limits) at 50% and 1% of the normalized peak, respectively. Therefore, we suggest that the variation of the silicate feature is related to the Herbig star’s protoplanetary disk and possibly material around it, while optically thin material lies beyond 40 au and could extend up to the presumed truncation radius of ≈ 65 au.

Previous studies (Cohen 1980; Cohen & Witteborn 1985; Polomski et al. 2005; Schütz et al. 2005; Do-Duy et al. 2020) have explored the silicate absorption feature using MIR long-slit spectroscopy and spectro-photometry. Their modeling of the feature with different dust mixtures, primarily silicates, and assigning black-body temperatures for the absorbing material, suggest optical depths τ_λ ≈ 0.2–0.5 and dust temperatures of 300–500 K. As such, our estimates for the optical depth agree with literature values, although it is worth noting that they were constrained to radii of ≥q20 au, and a moderate interstellar extinction correction has already been applied here. The optical depth would be much higher on smaller spatial scales near the protostar.

Several of the studies mentioned above suggested that the absorption feature originates from the FUor companion. The MATISSE images have shown that the FUor companion is indeed >10× fainter than the HBe primary. For the silicate feature to appear in emission, the total flux spectrum ought to be dereddened with A_V > 7 mag. From the reddening law of Gordon et al. (2023) (see Appendix A), the extinction at 10 μm becomes 0.57 mag and therefore the optical depth τ_{10 μm} is 0.52. This is somewhat higher than what is found from our analysis above, considering that the correlated spectra had already been dereddened for extinction A_V = 2.5 mag. If the interstellar extinction is not taken into account, our analysis shows that the maximum τ_λ is 0.35 on these large spatial scales.

As mentioned above, we argued that the absorption feature originates from the disk of the HBe star and it is contained within a radius of <40 au. This may contradict the scenario of a carved-out dusty cocoon (e.g., Canovas et al. 2012), while the blueshifted jet emission (Poetzel et al. 1989; Stelzer et al. 2009; Whelan et al. 2010; Antoniucci et al. 2016; Liimets et al. 2023) should originate from a much smaller region (<10 au). Furthermore, we note the lack of any large-scale, circumbinary warm dust emission in the reconstructed maps.

Evidently, the variability of the silicate feature, as suggested by the MATISSE results, cannot be fully explored with simple analytical models as used in literature studies, but it should be treated with radiative transfer simulations to study the disk’s geometry and the dust’s chemical composition, surface density, and grain size. We aim to explore this computationally expensive task in future work.

Fig. 7 Optical-depth uv maps of Z CMa. These correspond, from top to bottom, to 8.5, 10, and 11.5 μm, and the uv coverage is that of the small array (see also Fig. 3). The optical depth varies per wavelength, as indicated by each individual color bar, as well as per spatial scale probed by MATISSE (see Sects. 3.3 and 4.2).

Fig. 8 Optical depth variation at different physical radii. Three wavelengths are explored in the panels from top to bottom at 8.5, 10, and 11.5 μm, respectively. The radius is a 1/2 spatial scale probed by the interferometer (see Sect. 4.2) converted to physical units at the adopted distance. The optical depth is color-coded with regard to our parameter, η (see Sect. 4.2). The vertical gray dotted lines in the middle panel mark the upper limit of the radial extent of the HBe disk in the 10.4 μm image (Fig. 4) at 50% and 1% of the normalized peak.

Fig. 9 Best-fit orbital solution from orbitize using relative astrometry for Z CMa (Table E.1). Left panel shows one orbit (purple line, color-coded per year as shown in the color bar) with the relative astrometry of the FUor companion (orange points) and the location of the HBe star (asterisk). Right panels show the solution of the same orbit (gray lines) against the angular separation (top) and PA (bottom) that is also included in the left panel.

4.3 Binary orbit

Using the MATISSE astrometry and archival data (Table E.1), and based on an earlier analysis by Antoniucci et al. (2016) and Bonnefoy et al. (2017), we explored the binary’s orbit. Figure 9 shows the relative position of the FUor (orange points) with respect to the HBe star (black star), which is set at the origin (0,0). These relative astrometric measurements are essential in deriving the orbital elements and and estimating a mass ratio of Z CMa, assuming a fixed mass for one of binary components. It is worth noting that there have been no previous observational attempts to estimate masses for each component; for example, by measuring the radial velocities of each star through high-spectral-resolution studies. Since this is a wide binary, we expect that it would take several decades to collect radial velocity data. Hence, relative astrometry can offer the best solution here.

We used the tool orbitize⁸ by Blunt et al. (2020) to estimate the orbital elements with the Markov chain Monte Carlo (MCMC) solution. We opted for 200 walkers with 20 temperatures and 50 000 steps per walker after allowing a burn-in phase of 10 000 steps. The initial parameters required for this model were the relative astrometry, the source’s parallax, and the total mass of the binary. As stated earlier, we adopted the distance estimate of Dong et al. (2022), which we converted to a parallax (0.85 ± 0.05 mas). This parameter is fixed as a prior to minimize degeneracies. We also tested models where the parallax is a free parameter, and it was found to range within the estimates of Dong et al. (2022). Since the mass ratio of the system is unknown, we opted for a range of 10 ± 5 M_⊙ for the total mass, which should include the entire range of possible masses for each star. Figure 9 shows the best-fit orbital solution from orbitize. The median values of the orbital elements are the semi-major axis, a (${245.6}_{-42.9}^{+35.7}$ au), the eccentricity, e (${0.172}_{-0.057}^{+0.058}$), the inclination, i (${65.9}_{-4.4}^{+2.9}$ deg), the argument of the periastron of the secondary’s orbit, ω (${142.1}_{-30.5}^{+9.2}$ deg), and the longitude of ascending node, Ω (${200.1}_{-1.8}^{+1.7}$ deg). The system’s total mass, M_total, is estimated at ${16.4}_{-2.3}^{+2.1}$ M_⊙. All given errors represent the 68% confidence intervals from the MCMC fitting (Fig. E.1). The orbital period is estimated at ${950}_{-256}^{+218}$ years. The parameters ω and Ω cannot be constrained further in the absence of radial-velocity information (e.g., Nowak et al. 2024).

Our fit provides a conservative estimate of the system’s total mass. Assuming that the FUor companion is expected to be <2 M_⊙, the mass of the primary HBe would range from 12–16 M_⊙, and thus at the most massive end of its sequence. This could make Z CMa an early precursor of a supernova explosion.

There exist indirect estimates of the Herbig star’s mass; however, several other literature studies have inadvertently used wrong distances⁹ and/or cumulative photometric measurements in the SED analysis that did not take into account the source’s variability (Alonso-Albi et al. 2015; Vioque et al. 2018; Marcos-Arenal et al. 2021; Guzmán-Díaz et al. 2021; Stapper et al. 2022; Grant et al. 2023; Vioque et al. 2023; Thomas et al. 2023). van den Ancker et al. (2004) argued for a 16 M_⊙ mass for the Herbig star and approximately 3 M_⊙ for the FUor. More recently, Fairlamb et al. (2015) estimated a mass of 11.0 ± 1.7 M_⊙ for the HBe star by fitting theoretical spectra to their XSHOOTER results and for their adopted distance of 1050 ± 210 pc, which is closer to the lower end of our mass estimate for the HBe.

4.4 Nondetection of a tertiary companion

The detection of the streamer in NIR polarimetric imaging from the last decade (Millan-Gabet & Monnier 2002; Canovas et al. 2015; Liu et al. 2016) had been attributed to a stellar fly-by with a probable identification of the intruder at 1.3-mm by Dong et al. (2022). More recently, Zurlo et al. (2024) showed that the streamer feature is in fact part of the extended outflow of the system, while they disproved the presence of said intruder (otherwise a nondetection). The authors identified another background source at approximately $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}75$ northeast of the Z CMa binary as a tertiary companion. They claimed that the new source, evident as scattered light in the H band, as well as in the 1.3-mm continuum (Dong et al. 2022), was marginally detected in the L′ band. However, one could argue that the last detection is somewhat similar to spurious noise from PSF subtraction or that the source lies inside the jet of Z CMa.

The current relative astrometry does not allow for an evaluation of multiplicity for Z CMa beyond the known FUor companion. Moreover, this potential tertiary companion falls outside the MATISSE FOV in the L band, while in image reconstructions with the broader FOV of the N band, no other source is identified fainter than the FUor companion, the latter being more than ten times fainter than the primary Herbig star in the N band.

5 Conclusions

We present our MIR interferometric imaging results from MATISSE/VLTI of the young eruptive stars binary Z CMa, which were obtained during the 2023 EXor-like major outburst of the Herbig primary. This is the first part in a series of papers from these MATISSE observations.

We provide mono-chromatic images at several wavelengths essentially corresponding to continuum emission in the L, M, and N bands. The binary is resolved in all bands, and extended circumstellar emission is evident around each star. From the images we calculated relative brightness ratios in all three bands and find that the FUor is much fainter than the Herbig star, and more than 10x fainter in the N band in particular.

Optical depth maps indicate a nonuniform distribution of silicate dust on large spatial scales, or else within a region of 40–130 au in diameter, where the optical depth increases monotonically toward smaller radii. This distribution is more probably associated with the bright Herbig star and its protoplanetary disk and should be explored further with dust radiative transfer simulations. The dusty cocoon surrounding the Herbig star that was alleged by earlier NIR polarimetric imaging studies may be transparent in the MIR since emission is clearly seen from the protoplanetary disks of both binary components. We do not find any evidence for substructures, such as gaps, in the MATISSE images of the HBe disk, but this could be due to the limited angular resolution taking into account the object’s distance.

The MIR-emitting region of the FUor companion’s disk is rather small in all spectral bands compared to that of the primary, at an approximately constant angular size of ≤15 mas, which translates to a physical radius of <9 au. This is congruent with recent MATISSE results on FUor disks, which allude that these are rather compact, especially when the system is returning to quiescence. The HBe disk extends to radii of ≤8 au in the L band and ≪35 au in the N band. At longer wavelengths and lower flux scales, the disk is marginally smaller than the 1.3-mm optically thin emission in ALMA studies (Dong et al. 2022). The disk position angles at 152.5° ± 12.7° and 156.5° ± 12.1° and inclinations at 33.5° ± 5.7° and 26.6° ± 3.8° for the FUor and HBe, respectively, suggest that the disks are aligned, and they agree with the orientations extrapolated from the jets (Whelan et al. 2010; Antoniucci et al. 2016).

From the MATISSE reconstructed images, we calculated a binary separation of 117.88 ± 0.73 mas at a position angle of 139.16° ± 0.29°. This is complementary to approximately 20 years of relative astrometry data. Our orbital fitting provided estimates for the total mass of the system M_total = ${16.4}_{-2.3}^{+2.1}$ M_⊙ and a period approximately 950 years. The fit orbital inclination (i ~ 66°) is somewhat different to the calculated inclinations, within the large uncertainties, of the disks from 1.3-mm and 9-mm images (Dong et al. 2022), as well as those estimated from the MATISSE images. Due to the large uncertainties, we cannot ascertain whether there is a misalignment between the millimeter continuum and MIR wavelengths. We expect that the orbital solutions can be constrained better with future astrometric observations (e.g., Gaia DR4).

Z CMa is a highly complex system with an evolving Herbig star, and an FUor companion that may have returned to quiescence. Future MIR observations of the system during the Herbig star’s quiescence with MATISSE/VLTI and METIS/ELT may reveal more regarding extended circumbinary emission. We encourage new ALMA observations at higher angular resolution to investigate the discrepancy between MATISSE and ALMA disk sizes. In future works as part of this series, we intend to employ radiative transfer methods to investigate the geometric structure and chemical composition of the HBe disk, and exploit the spectro-interferometric MATISSE data.

Acknowledgements

We would like to thank the anonymous referee for their feedback that helped improve this manuscript. This work received funding from the Hungarian NKFIH OTKA project no. K-132406 and K-147380, and the NKFIH Excellence grant TKP2021-NKTA-64. This work was also supported by the NKFIH NKKP grant ADVANCED 149943. Project no. 149943 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the NKKP ADVANCED funding scheme. FCSM received financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Starting Grant “Chemtrip”, grant agreement No. 949278). Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes: 0106.C-0501(B), 0108.C-0385(D), and 0110.C-4209(A). Spectroscopic observations were obtained as Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration This research has made use of the Jean-Marie Mariotti Center OImaging service (Available at http://www.jmmc.fr/oimaging) part of the European Commission’s FP7 Capacities programme (Grant Agreement Number 312430), and LITpro service co-developed by CRAL, IPAG and LAGRANGE (LITpro software available at http://www.jmmc.fr/litpro). We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. Based on data from the OMC Archive at CAB (INTA-CSIC), preprocessed by ISDC and further processed by the OMC Team at CAB. The OMC Archive is part of the Spanish Virtual Observatory project. Both are funded by MCIN/AEI/10.13039/501100011033 through grants PID2020-112949GB-I00 and PID2019-107061GB-C61, respectively.

Appendix A Interstellar extinction

The current literature includes several estimates for the interstellar extinction, A_V, toward Z CMa (Table A.1). The analysis of Bayestar19 (Green et al. 2019) and STILISM (Lallement et al. 2019) indicates the presence of a foreground cloud within 1 kpc. This could explain the extremely high extinction values estimated by IRAS and Planck satellites, whose wide beams would have engulfed not only the foreground cloud but also dust from the surrounding CMa R1 star forming region. The large uncertainties by Connelley & Reipurth (2018) are a result of their fitting the spectrum of Z CMa to that of FU Orionis. The protostars themselves are dusty, therefore we cannot exclude that the remaining estimates in Table A.1 may be affected by circumstellar and/or circumbinary extinction as well. However, the contribution for each protoplanetary disk is unclear, since both are visible in the infrared and (sub)millimeter wavelengths. The large uncertainties by Connelley & Reipurth (2018) are a result of their fitting the spectrum of Z CMa to that of FU Orionis, while Hinkley et al. (2013) attempted to fit the NIR photometry of each star with reddened blackbody spectra to estimate circumstellar extinction.

In this work, we adopt an interstellar extinction of A_V = 2.5 mag. This is the average of all values where A_V ≤ 5 mag in Table A.1, and it is in agreement with the extinction used in Sicilia-Aguilar et al. (2020). When necessary, we used the extinction law of Gordon et al. (2023) to deredden spectrophotometry with R_V = 3.1, which for our adopted A_V, results in extinction of A_{10 μm}/A_V ≃ 0.08 or else A_{10 μm} ≃ 0.2 mag for the silicate feature.

Appendix B MATISSE data quality

A summary of the MATISSE/VLTI data is shown in Table B.1. Observational blocks that were not successful and/or aborted due to adverse atmospheric conditions and/or technical faults at the telescope are also indicated as ‘X’, and these datasets are discarded from this work. The GTO observations from 2021 were obtained in MED-L, LOW-L, and HIGH-N (id=1) and LOWL and LOW-N (id=2) spectral resolutions, while the imaging observations were all obtained in MED-L and LOW-N (ids=3–23). Only the LOW-L and LOW-N data from id=1 were used in this work.

GRA4MAT and LM band. In the majority of cases, chopped photometry failed in the LM band even during the best observing conditions (i.e., atmospheric coherence times t_coh ≥ 3 ms) and even with the bright calibrator (α CMa). It is unclear why this is the case for GRA4MAT even with a bright science target like Z CMa (F_{3.5 μm} = 40 Jy). Therefore, only the non-chopped LM-band (and subsequently M) data is used in this work. This issue does not affect the N-band data. To validate the photometric accuracy of the non-chopped data from the 1.8-m ATs in all bands, we compared the averaged MATISSE spectra against the archival M and N band photometry from the Midsource Space Experiment (MSX) and AKARI space telescopes, as well as our own SpeX/IRTF spectroscopy. We found that all are in good agreement (Fig. 2). We also note that Z CMa saturates (NEO)WISE photometry and therefore these measurements are unreliable.

GRA4MAT and N band. The standard MATISSE pipeline calibration (ver.1.7.5 and later) failed for the case of the N-band data with GRA4MAT. We suspect that this was caused by the optical path difference (OPD) selected by GRAVITY for GRA4MAT (co-phasing) in the N-band. To correct for this, we reduced the N-band data with an earlier version of the pipeline (ver. 1.6) that is used in the standalone MATISSE mode. This issue does not affect the LM-band data. Furthermore, GRA4MAT does not record simultaneous photometry for the N band. This can be recovered instead from the interferometric data by dividing the correlated fluxes with the visibilities. The quality of this data has been verified to be in good agreement against archival MIDI spectra (Varga et al. 2018) and archival photometry (Fig. 2).

Rotational flip in the closure phases. The phase sign ambiguity for the N-band data that was mentioned by Gámez Rosas et al. (2022) was not resolved in the MATISSE DRS pipelines used for this work (ver. 1.7.5 and earlier). This was confirmed in our data, since the orientation of the Z CMa binary is well known (see Sect. E). Therefore, all our models and image reconstruction results in the N band have been rotated by 180° to correct for this issue.¹⁰

Choosing calibrators. Sirius (α CMa) was selected as a calibrator for our small array observations (ids: 7–11, 20, and 21 in Table B.1) due to its high relative brightness for the ATs. Sirius is a binary system, but a wide one with a current separation ≫ 2″ and therefore outside the ATs FOV. The white dwarf companion (Sirius B) is not known to emit in the MIR. Unfortunately, due to its proximity to the Sun, the A-type main sequence star (Sirius A) can be resolved with the small array (⊘ ≈ 6.14 mas; Bourgés et al. 2014). Moreover, the hydrogen spectral signatures of Sirius are clearly detected in the calibrated data (single-dish spectra and correlated fluxes). Although the hydrogen lines in Sirius A’s spectrum are in absorption, since the calibrated data are processed by dividing the response of the science target with that of the calibrator, they appear here as emission lines (cf. Fig. B.1 and B.2). Hence, when necessary in our analysis, we either subtract the hydrogen lines from the small array data by replacing them with the average continuum or we select a continuum region that does not include these lines. Image reconstruction is not affected, since this is based only on the closure phases and visibilities where we do not see any influence by Sirius’ spectral response.

Fig. B.1 MATISSE L-band correlated fluxes. These are color-coded w.r.t. baseline (length, PA) as listed within each panel. From left to right and top to bottom each panel corresponds to sequential observations with the array configurations listed in Table B.1, that is ids=3, 5-13, 15-16, and 18-23.

Fig. B.2 As in Fig. B.1 for the M-band correlated fluxes.

Fig. B.3 As in Fig. B.1 for the N-band correlated fluxes.

Fig. B.4 MATISSE L-band squared visibilities from the imaging data. The color-coding and sequential order are as in Fig. B.1.

Fig. B.5 As in Fig. B.4 for the M-band squared visibilities.

Fig. B.6 As in Fig. B.4 for the N-band squared visibilities.

Fig. B.7 MATISSE L-band closure phases from the imaging data. These are color-coded per baseline triangle as shown in each panel. Each panel is sequentially ordered as in Fig. B.1.

Fig. B.8 As in Fig. B.7 for the M-band closure phases.

Fig. B.9 As in Fig. B.7 for the N-band closure phases. Here, the phase convention is flipped w.r.t. to the LM data.

Appendix C Image reconstruction methods

Unlike radio interferometry where there exist robust image reconstruction algorithms, which take advantage on the well-calibrated information of the visibilities’ phase, the sparse uv coverage and absence of absolute (visibility) phase measurements is a well-known problem that affects image reconstruction in infrared interferometry. Furthermore, there exist multiple image reconstruction techniques to tackle such data. It is thus quite common practice to compare images produced from at least two image reconstruction algorithms. Here, we showcase the differences found between BSMEM ver.2.3 (Buscher 1994) and MiRA ver.2.3.2 (Thiébaut 2008) algorithms on OImaging for the same set of MATISSE data, which will explain our final choice on selecting an algorithm throughout this work (i.e., MiRA). The comparison is made for the L (3.5 μm) and N (8.5 μm) band data.

Figure C.1 shows image reconstructions with BSMEM (left panels) and MiRA (right panels) for continuum emission in the L band (top row; 3.4–3.6 μm; namely 3.5 μm) and the N band (bottom row; 8–9 μm; namely 8.5 μm). The spectral bandwidths were chosen to resemble ‘pseudo-continuum’ ranges without any spectral lines. The FOV is 400×400 mas with a pixel scale 1 mas at 3.5 μm and 3 mas at 8.5 μm. North is oriented up and east left, while the theoretical MATISSE beam is shown as a white circle in the left panels. All images have been convolved with a Gaussian beam with FWHM similar to the image resolution, and they are normalized to the peak intensity, and scaled in a logarithmic stretch to a minimum level of 10⁻³ as indicated in the colorbar. This enhances emission from the faint FUor as well as from the residuals.

In OImaging the MiRA algorithm forces the brightest component (here the HBe star) in the middle of the map, while BSMEM applies a random shift in the map. Both methods produce image reconstruction residuals due to the incomplete uv coverage of the observations (Planquart et al. 2024, provides a thorough treatise on MATISSE image reconstructions). In both cases this is manifested by multiple blobs throughout the map, while the ‘floral pattern’ of the uv coverage is seen in the MiRA reconstructions. Although BSMEM can accurately predict the location and orientation of the binary’s components, the individual components have rather irregular patterns reminiscent of a patchwork (at least for the HBe star) in the non-convolved images due to the ‘super-resolution’ effect. Similarly, MiRA presents folded-like artifacts around the HBe disk at 3.5 μm, but it is overall consistent in all bands. Overall, both methods can predict the location of the binary companion and a similar average size for both sources, therefore we trust that our image reconstruction is accurate taking into account the limited uv-coverage of our data.

Fig. C.1 Comparison between BSMEM (left) and MiRA (right) reconstructions at 3.5μm (top row) and 8.5μm (bottom row). Each reconstruction image is normalized to the peak intensity, and scaled in a logarithmic stretch to a minimum level of 10−3 as indicated in the colorbar. The FOV is 400 × 400 mas with north up and east toward left. The theoretical MATISSE beam is shown as a white circle in the lower-left corner of the left panels. A Gaussian smoothing with a FWHM equal to the beam has been applied throughout to enhance the FUor’s disk.

Appendix D Optical depth variability at large scales

The optical depths calculated for each individual, small-array baseline as in Sect. 3.3 are shown in Fig. D.1. The spatial scales probed correspond to angular sizes of approximately 30 to 100 mas, which can be converted to physical scales of approximately 34-115 au. The binary’s separation is < 120 mas, or else < 135 au and thus larger than the spatial scales probed by the small array. At the time the MATISSE observations were obtained, the pointing accuracy was not optimized to focus on any of the two stars. The imaging results show that the FUor disk is rather small and faint compared to the HBe disk, while there is no apparent extended or circumbinary warm dust emission in the MIR (e.g., Monnier et al. 2009). Here, we argue that the spatial scales probed by the small array should correspond to the outer regions of the Herbig star’s protoplanetary disk.

Fig. D.1 Optical depth (τ) per baseline. The details of each baseline name, projected length and position angle, are stated within each panel. Each row corresponds to a different small-array observation and their ids (blue) can be found in Table B.1.

Appendix E Relative astrometry

Table E.1 includes the relative astrometry of the FUor star w.r.t. the HBe star from literature data, as well as our own derivation from our MATISSE observations, and the literature (Leinert & Haas 1987; Koresko et al. 1991; Haas et al. 1993; Malbet et al. 1993; Barth et al. 1994; Thiebaut et al. 1995; Fischer et al. 1998; Millan-Gabet & Monnier 2002; Ratzka 2005; Szeifert et al. 2010; Canovas et al. 2015; Antoniucci et al. 2016; Bonnefoy et al. 2017; Takami et al. 2019; Dong et al. 2022; Thomas et al. 2023). For the case of Szeifert et al. (2010), we provide the correct epochs of the NACO observations as extracted from the ESO archive. Following the analysis of Ratzka (2005) on the 2004 MIDI/VLTI data, who modeled the binary as two uniform disk sources, we employed the same method to estimate the orbit from all archival MIDI/VLTI observations (Varga et al. 2018). We supplement Table E.1 with our analysis of archival, raw, acquisition data in the NIR and MIR from NIRC2/Keck (program IDs: C216N2, H205N2L, N112), and NIRI/Gemini (program ID: GN-2010B-Q-74). For the case of NIRI/Gemini and NIRC2/Keck observations, the binary is clearly resolved even in raw acquisition data. Here, we repeat the same method as for the MATISSE images. The uncertainties are the standard deviations of the fitted separations and position angles. The analysis of the relative astrometry data is presented in Sect. 4.3. Figure E.1 is the MCMC fitting corner plot from orbitize.

Fig. E.1 Corner plot of the MCMC fitting from orbitize for the model presented in Section 4.3.

References

All Tables

All Figures

Fig. 1 V-band light curve of Z CMa with corresponding epochs of MATISSE observations (orange lines; dashed and solid for the quiescent and the outburst phase, respectively) as in Table B.1. Photometric data from the Archive of the American Association of Variable Star Observers (AAVSO; green circles), the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL) Optical Monitoring Camera (OMC; magenta circles), and the Kamogata/Kiso/Kyoto Wide-field Survey of Variable Star Observers League in Japan (KWS; blue circles Maehara 2014). The imaging data were obtained during the outburst from late 2022 to early 2023. Also shown for reference are the epochs of archival ALMA and JVLA observations (purple lines; Takami et al. 2019).

In the text

	Fig. 2 Dereddened SED of Z CMa during the 2023 outburst. Our own SpeX/IRTF and MATISSE spectra are shown in red and blue, respectively. For a description of the archival photometry (black points), we refer the reader to Sect. 2.1. A Kurucz model (green) and a blackbody model (orange) are shown for reference.
In the text

Fig. 3 MATISSE uv-plane coverage of “imaging” mode observations of Z CMa during its outburst phase. The uv coordinates are colored per respective VLTI configuration array as in Table B.1; that is, the small array included the A0-B2-D0-C1 quadruple, the medium array included the K0-G2-D0-J3 quadruple, and the large included the A0-G1-J2-K0 quadruple, which was interchanged with the A0-G1-J2-J3 to maximize uv coverage.

In the text

Fig. 4 MiRA image reconstructions in L, M, and N bands during the 2023 outburst. From left to right and top to bottom, these correspond to pseudo-continuum reconstructed images at 3.5, 4.7, 8.5, 9.2, 10.4, and 11.5 μm. These are smoothed with a Gaussian beam with a FWHM similar to the image resolution, normalized to unity, and shown in a logarithmic stretch with a minimum value of 10−2 to suppress image reconstruction residuals below that level. Such residuals are the result of insufficient uv coverage (Fig. 3). The SE component (FUor) becomes much fainter than the NW component (HBe) at longer wavelengths. The beam size is shown as a black circle at the bottom left corner of each panel.

In the text

	Fig. 5 Comparison between ALMA and MATISSE. The normalized pseudo-continuum reconstructed image at 3.5 μm is shown in blue. Overlaid are contours (gray lines) at 5, 10, 50, and 90% of the peak of the normalized 9.5 μm reconstructed image. These are compared to the ALMA 1.3-mm continuum disk sizes (blue dashed-line ellipsoids) of Dong et al. (2022).
In the text

	Fig. 6 Normalized radial profiles of HBe disk from reconstructed images in Fig. 4. The horizontal dotted line demarcates the relative flux threshold set for the image reconstructions. The vertical dashed line demarcates the location of the first plateau at 3.5 μm and, consequently, the region from which the majority of the emission originates at all wavelengths.
In the text

	Fig. 7 Optical-depth uv maps of Z CMa. These correspond, from top to bottom, to 8.5, 10, and 11.5 μm, and the uv coverage is that of the small array (see also Fig. 3). The optical depth varies per wavelength, as indicated by each individual color bar, as well as per spatial scale probed by MATISSE (see Sects. 3.3 and 4.2).
In the text

Fig. 8 Optical depth variation at different physical radii. Three wavelengths are explored in the panels from top to bottom at 8.5, 10, and 11.5 μm, respectively. The radius is a 1/2 spatial scale probed by the interferometer (see Sect. 4.2) converted to physical units at the adopted distance. The optical depth is color-coded with regard to our parameter, η (see Sect. 4.2). The vertical gray dotted lines in the middle panel mark the upper limit of the radial extent of the HBe disk in the 10.4 μm image (Fig. 4) at 50% and 1% of the normalized peak.

In the text

Fig. 9 Best-fit orbital solution from orbitize using relative astrometry for Z CMa (Table E.1). Left panel shows one orbit (purple line, color-coded per year as shown in the color bar) with the relative astrometry of the FUor companion (orange points) and the location of the HBe star (asterisk). Right panels show the solution of the same orbit (gray lines) against the angular separation (top) and PA (bottom) that is also included in the left panel.

In the text

	Fig. B.1 MATISSE L-band correlated fluxes. These are color-coded w.r.t. baseline (length, PA) as listed within each panel. From left to right and top to bottom each panel corresponds to sequential observations with the array configurations listed in Table B.1, that is ids=3, 5-13, 15-16, and 18-23.
In the text

	Fig. B.2 As in Fig. B.1 for the M-band correlated fluxes.
In the text

	Fig. B.3 As in Fig. B.1 for the N-band correlated fluxes.
In the text

	Fig. B.4 MATISSE L-band squared visibilities from the imaging data. The color-coding and sequential order are as in Fig. B.1.
In the text

	Fig. B.5 As in Fig. B.4 for the M-band squared visibilities.
In the text

	Fig. B.6 As in Fig. B.4 for the N-band squared visibilities.
In the text

	Fig. B.7 MATISSE L-band closure phases from the imaging data. These are color-coded per baseline triangle as shown in each panel. Each panel is sequentially ordered as in Fig. B.1.
In the text

	Fig. B.8 As in Fig. B.7 for the M-band closure phases.
In the text

	Fig. B.9 As in Fig. B.7 for the N-band closure phases. Here, the phase convention is flipped w.r.t. to the LM data.
In the text

Fig. C.1 Comparison between BSMEM (left) and MiRA (right) reconstructions at 3.5μm (top row) and 8.5μm (bottom row). Each reconstruction image is normalized to the peak intensity, and scaled in a logarithmic stretch to a minimum level of 10−3 as indicated in the colorbar. The FOV is 400 × 400 mas with north up and east toward left. The theoretical MATISSE beam is shown as a white circle in the lower-left corner of the left panels. A Gaussian smoothing with a FWHM equal to the beam has been applied throughout to enhance the FUor’s disk.

In the text

	Fig. D.1 Optical depth (τ) per baseline. The details of each baseline name, projected length and position angle, are stated within each panel. Each row corresponds to a different small-array observation and their ids (blue) can be found in Table B.1.
In the text

	Fig. E.1 Corner plot of the MCMC fitting from orbitize for the model presented in Section 4.3.
In the text