ELT-METIS imaging simulations for disks and envelopes associated with FU Ori-type objects

In Paper I, we developed a semi-analytic method to calculate extended MIR emission for existing FUors with an order-of-magnitude accuracy. The calculations are made using (1) the observed distribution of the polarized intensity (PI) in the $H$-band ($\lambda$=1.65 µm); (2) the observed SEDs at ultraviolet (UV) to MIR wavelengths; and (3) assumed optical properties for dust grains in a disk or an envelope, which produces extended infrared emission. Our method allows us to easily calculate the emission distributions for various cases. These cases are specifically: (A) when the central radiation source at UV to NIR wavelengths ($\lambda<3$ µm) is a flat compact self-luminous disk (e.g., Hartmann & Kenyon, 1996; Zhu et al., 2008) or a star (e.g., Herbig et al., 2003; Elbakyan et al., 2019), (B) when the infrared extended emission is associated with a disk or an envelope, and (C) with different dust models for light scattering and thermal radiation from the extended disk or envelope. This semi-analytic method is complementary to full radiative transfer simulations, which offer more accurate calculations but only with specific dynamical models and significant computational time.

We used some assumptions and simplification for the geometry of the extended disks/envelopes to derive the MIR images. This makes our calculations less accurate, but still with accuracies sufficient for our purpose, that is, determining whether the extended MIR emission is observable using HCI. In Paper I, we performed comparisons with some numerical simulations, and these suggest that our semi-analytic method yields intensity distributions with an accuracy within a factor of 2. However, this method does not include radiation heating in the inner disk edge or adiabatic heating, which potentially enhances thermal emission at $\lambda\sim 10$ $\micron$ (see also Section 4 for future work).

As in Paper I, we calculated the MIR emission for two FUors, FU Ori and V1735 Cyg. We tabulate the key parameters for these objects in Table 2. Figure 1 shows the images in $H$-band observed using the Subaru-HiCIAO instrument. For each object, we set the intensity within 0$\farcs$2 of the central source to be zero as we were not able to make reliable measurement due to the central source being significantly brighter than the extended emission. We note that this software mask is significantly larger than the inner working angle of the instrument tabulated in Table 1. In each panel, we also show a few contours to indicate approximate emission distributions. These contours will be used in Section 3.1 to perform comparisons with the point-spread functions (PSFs) of the central source and their residuals for subtraction for the METIS observations.

For (2) described above, we used the SEDs described in Paper I. In practice, V1735 Cyg cannot be observed from the ELT site due to its high declination. However, we will still use this target to investigate the detection of MIR emission associated with FUors in a more general context.

The infrared extended emission ($r>10$ au), either observed or to be observed, must be due to dust scattering and thermal dust radiation in the disk or envelope (Section 1; see also Whitney et al., 2003a, b, 2004, 2013; Audard et al., 2014). At UV to NIR wavelengths ($\lambda<3$ µm), dust grains in the extended disk or envelope are illuminated and heated by radiation from the central source, which is either a disk or a star as described above. At MIR wavelengths ($\lambda>3$ µm), the observed SEDs suggest that the radiation from the central source is dominated by the former (Liu et al., 2016, Paper I).

We approximated an extended disk with an optically thin atmosphere and an optically thin interior. For the disk, we then calculated emission via single scattering (the central source$\rightarrow$scattering in the disk atmosphere$\rightarrow$the observer), double scattering (the central source$\rightarrow$the first scattering in the disk atmosphere$\rightarrow$the second scattering in the disk interior$\rightarrow$the observer), and thermal emission from the disk atmosphere and the interior. The light with more than two scatterings, which is not included in our calculations, would enhance the total intensity by only $\lesssim$5 %. This estimate is based on the fact that light even via double scattering contributes to the total intensity only up to $\sim$20 % as described below. For the extended envelope, we calculated the scattered and thermal emission assuming that the envelope is optically thin and therefore double scattering and self absorption are negligible.

For the optical properties of dust grains, we used three models as for Paper I. These dust models were originally developed for the interstellar medium without ice coating (‘Dust1’), for a molecular cloud with ice coating (‘Dust2’), and for the surface of the disk associated with HH 30 (‘Dust3’), respectively and were used in Whitney et al. (2003a, b, 2004, 2013). Each dust model uses a number of homogeneous spherical particles with “astronomical silicates” (Draine & Lee, 1984) and graphite, with the size distributions adjusted to reproduce some observations. Whitney et al. (2003a, b) calculated the optical parameters (dust opacities, scattering albedos, forward throwing parameters, and polarization for scattered light) for various wavelengths for these dust models based on the Mie theory and the geometrical optics algorithm (Wood et al., 2002). See Paper I for more details on these dust models and how to calculate the distributions of total intensity (Stokes $I$) at MIR wavelengths from polarized intensity distributions in the $H$-band.

The dust opacities for these models are tabulated in Table 3. For this paper, we used the Dust1 and Dust2 models for the extended disks as well as the extended envelopes, as the Dust3 model cannot explain the optical properties of all of the known circumstellar disks associated with young stars.

The calculated distributions of total intensity for the three MIR wavelengths ($\lambda$=3.8, 4.8, and 11.3 µm) are dominated by the single scattering process, as was the case for those executed at similar wavelengths in Paper I. Figure 2 shows the calculated intensities for the single scattered light normalized by the following case: when both the bright central emission and the extended emission are from a disk with ‘Dust1’. The scattered emission is brighter for the following cases: (1) observations at short wavelengths; (2) when the extended emission is due to a disk rather than an envelope; (3) when the central radiation source is also a disk, not a star at any wavelength; and (4) for ‘Dust2’ rather than ‘Dust1’, and ‘Dust3’ rather than ‘Dust2’. Double scattering in the disk enhances the emission by up to $\sim$20 % only. The thermal emission is responsible for up to $\sim$10 % at $\lambda$=11.3 µm in the regions close to the central source, but is negligible for all the other cases.

When the central radiation source at NIR wavelengths is a star, the calculated intensities are also affected by the assumed typical grazing angle of the disk surface or the envelope $\overline{\gamma}$ ($\sim\overline{z/r}$, where $z$ and $r$ are the height and radius, respectively). In contrast, the intensity for single scattering is independent of $\overline{\gamma}$ if the central radiation source is a disk at all wavelengths. Therefore, the calculated intensity including all the radiation processes is nearly independent of $\overline{\gamma}$ for the latter case. The intensity derived using a star is smaller than that using a disk by approximately a factor of $\overline{\gamma}$. For this paper, we used $\overline{\gamma}$=0.1 as a lower limit to yield self-consistent calculations. See Paper I for details about the effects and limitations of $\overline{\gamma}$.

As described in Section 2.1, the bright central source is a self-luminous compact disk at the wavelengths of the HCI observations. Although this is nearly a point source at the given angular resolution, the emission from the outer disk radii leaks out from the coronagraphic mask, and therefore enhances the speckle noise, as demonstrated in later sections.

For this work, we used the following disk models: (1) conventional optically-thick accretion disk models (Pringle, 1981); and (2) a Gaussian distribution with a HWHM of 1, 2, and 4 au. Interferometric observations of disks at NIR and MIR wavelengths have been powerful in constraining their spatial distributions, but have not yet been able to determine detailed distributions as described below. For (1), the intensity distribution depends on the wavelengths of the observations and the mass accretion rate (and therefore the resultant luminosity). This model successfully explained the observed SEDs and MIR interferometric visibilities ($\lambda\sim 10$ µm) for FU Ori (Quanz et al., 2006). This study was corroborated by Labdon et al. (2021), who executed multi-band NIR interferometry ($\lambda$ = 1.2-2.2 µm) for the same star. The models in (2) were used for the following studies.  Liu et al. (2019); Lykou et al. (2022) attributed their NIR-to-MIR interferometric observations ($\lambda$ = 2-3.5 µm) of FU Ori to a Gaussian disk with a HWHM of 0.5-3 au. Quanz et al. (2006) demonstrated that the MIR interferometric visibilities ($\lambda\sim 10$ µm) measured for FU Ori are consistent with a Gaussian disk with a HWHM of $\sim$4 au as well as the conventional disk model.

Figure 3 shows the radial intensity distribution for the above disk models in three HEEPS bands. The distribution of emission for the conventional disks was smaller than any of the Gaussian distributions we used. The conventional disks show a larger spatial extent at longer wavelengths because of contributions from the cooler outer regions. Their spatial extent is slightly larger for FU Ori than V1735 Cyg due to a larger disk luminosity.  For all the cases, most of the flux from the central disk is distributed within the diffraction core of the telescope.  As a result, their spatial distributions are almost identical to those of a point source without a coronagraph.

The inclination angle of the disk+envelope system is not known for our target objects (Paper I). For our simulation, we use these compact disks with a face-on view, which yields the maximum leakage from the coronagraphic mask and therefore conservative detection limits.

To investigate the signal-to-noise ratio for the extended MIR emission, we used HEEPS v1.0.0¹¹1https://github.com/vortex-exoplanet/HEEPS (Delacroix et al., 2022a), an open-source python-based software for HCI (Delacroix et al., 2022b). HEEPS calculates the speckle noise associated with bright sources in the given detector format through the following steps: (1) derivation of a temporal series of single-conjugate adaptive optics (SCAO) residual phase screens, including the predicted instantaneous pointing errors of the observations; (2) propagation of the SCAO residual phase screens through the individual optical components of the instrument, using one of the vortex coronagraphy modes (CVC or RAVC); (3) accumulation of the instantaneous coronagraphic PSFs to produce a mock observing sequence in pupil-stabilised mode, including the METIS radiometric budget for the considered target star in the considered filter; and (4) computation of the post-processed contrast. For each filter, these calculations are made monochromatically for the given representative wavelength.

Using HEEPS, we first calculated the following sequences (cubes) of the PSFs for the default field of view (FoV) of HEEPS of 2.2$\times$2.2 arcseconds:

(I) On-axis PSF cubes for the central sources described in Section 2.2. We calculated the angular distribution of the emission, and derived the PSF cubes for an on-source integration of 30 minutes for each wavelength, coronagraph (CVC or RAVC), and disk model. Each cube consists of 6000 images for the given integration and time interval.
HEEPS has a built-in function to execute such simulations for slightly extended emission such as these disks, approximating them using multiple point sources. However, this function requires significant computational time proportional to the number of the point sources. To overcome this problem, we developed Monte-Carlo simulators to eject photon packets from the individual locations of the above disk models. We added the coordinate of each photon packet to the pointing error for the simulation for each image of the cube provided with a 300-msec interval. This method allowed us to calculate a PSF cube for each case with a computational time of approximately 45 minutes using an 8-12 CPU server. The use of 6000 photon packets for a 30-minute integration reduced the statistical error for the intensity distribution of the PSF to within 2 %.

(II) Off-axis PSF cubes used to convolve extended MIR emission derived in Section 2.1. We also used the off-axis PSF at $\lambda$=3.8 µm to study the impact of a binary companion (Section 3.3).

(III) On-axis PSF cubes for reference stars (point sources) to subtract the bright central emission from the target exposures (see Appendix A for necessity). For these cubes, we used different sets of phase screens from (I)(II) in order to reproduce the observations at different times. We used a star 3 times as bright as each target object, and executed simulations for a 30-minute integration which yielded 6000 frames. In Section 4, we will briefly discuss the use of stars with different brightness levels at varying integrations.
As described below, we derived median PSFs using the PSF cubes and subtracted them from the individual science frames. In the future, we may use PSF models with the Locally Optimized Combination of Images (LOCI; Lafrenière et al., 2007) and Principal Component Analysis (PCA; Soummer et al., 2012; Amara & Quanz, 2012; Juillard et al., 2024) that have previously been successful when the central science target is a star.

To calculate speckle noise, we used the SCAO residual phase screens derived with a 300 msec sampling for the following conditions, as in Carlomagno et al. (2020); Delacroix et al. (2022a). We split the currently available 12000 phase screens into two to obtain the cubes for (I)(II) and (III), respectively. The mock observations were made for median atmospheric conditions at the ELT site, at an assumed declination of $-$5° and with a $K$-band ($\lambda$=2.2 µm) magnitude of 5 for SCAO corrections. The actual target $K$ magnitudes of 5-7 for FU Ori and V1735 Cyg yield similar performances as our mock observations (Feldt et al., 2024). The assumed declination yields a rotation of the parallactic angle of about 20^∘ during the given integration when the target was crossing the meridian. For the thermal background, we used constant values estimated for the individual filters tabulated in Table 1. We approximated that the PSF subtraction described above does not yield any residual patterns for the thermal background except for the photon noise.

Using the cubes (I)-(III), one would execute the following simulations for extended emission in order to approximate the actual observations in the sky: (A1) apply field rotations to the extended emission for the mock observing sequence; (A2) convolve the images of the extended emission using the off-axis PSFs calculated for the individual time sequences; (A3) add the on-axis PSF for the central source; (A4) subtract the stacked reference PSF; (A5) de-rotate the images to match their sky coordinates; and (A6) stack the cube to create the final image. However, this sequence also requires significant computational time due to the number of combinations of the bands, coronagraphic modes, central sources, and extended emission in our study. In practice, the process (A5) returns the location of the extended emission to its original position before applying the field rotations (A1). Therefore, we alternatively executed the following processes, which yielded identical results but with a significantly smaller computational time: (B1) subtracted the stacked reference PSF from each image in the on-axis PSF cube; (B2) derotated the images in the on-axis and off-axis PSFs, respectively, to match their sky coordinates; (B3) stacked the on-axis and off-axis PSF cubes, respectively; and (B4) convolved the extended emission using the stacked off-axis PSF and added it to the stacked on-axis PSF.

Before the PSF subtraction, we stacked the reference PSF without field rotation, and scaled it using the fluxes of the target PSF (without the extended emission) and the reference PSF measured within a radius of 0$\farcs$5. This scaling process yielded the cases when the ideal PSF subtraction was achieved. In contrast, the images obtained through actual observations of the targets would include both the on-axis PSF and extended emission, and therefore require complicated optimization of the PSF subtraction to avoid over-subtraction. This optimization is beyond the scope of this paper, that is, to investigate whether the target emission is actually detectable over speckle and photon noise.

The photon noise was added to the final images as follows. Before the process (B1), we first obtained the image of the photon counts for the stacked reference PSF with the thermal background, derived the image of the Poisson noise, and added it to the stacked reference PSF. Secondly, we obtained the image of the photon counts for the target frames with the thermal background, through (B2)-(B4), but without PSF subtraction. We then derived the distribution of the Poisson noise and added it to the image derived through (B1)-(B4). This method yields the same level of photon noise as we add Poisson noise to each image with 300-msec sampling, but with significantly fewer random numbers and therefore total time for calculation.

We then converted the units of the extended emission to W m^-2 µm^-1 arcsec^-1. In addition to the images with extended emission, we also created images using the same PSFs but without extended emission. We used these images to estimate the detection limits of the extended emission in later sections.

The read noise of the detector, which is not included in our calculations, is marginal or negligible compared with the total noise if we select a sufficiently long exposure. Under these conditions, the results presented in the rest of the paper are independent of the actual exposure time. Such exposures saturate the central source in limited cases, but do not significantly affect the observations of the extended emission in which we are interested (Appendix B). To minimize the overheads for the detector readouts, and therefore the total time for the observations, one would select an actual exposure time to be as long as possible without the detector saturation significantly hampering the image alignment for the image stacking process.

Figure 4 shows how the two coronagraphs work to reduce the flux from the bright central disk in various cases, that is, for different disk models and three bands, and when the target is FU Ori or V1735 Cyg. These are shown as the flux ratios of the on-axis (when the emission is centered on the coronagraphic mask) to the off-axis observations (when the emission is out of the coronagraphic mask) measured within $r$=5$\lambda$/D of the central source (0$\farcs$12, 0$\farcs$15, and 0$\farcs$37 at $\lambda$=3.8, 4.8, and 11.3 µm, respectively; where $D$ is the diameter of the telescope). This region for each wavelength would cover 90–95 % of the total flux without a coronagraph. The derived ratios, which differ from those using the conventional method with the central position only (“peak rejection rate”; e.g., Carlomagno et al., 2020), allow us to evaluate the flux reduction independent of different shapes of PSFs.

In each panel of Figure 4, we also show the flux ratios for a point source for reference. The horizontal axis of each panel is organized from left to right for small to large spatial extent of the bright central source (Section 2.2; Figure 3). Speckle noise and the pointing errors for the AO were included but photon noise is not.

In general, a large spatial extent for the bright central source yielded a larger flux leakage from the mask. However, the differences were marginal between the following cases: (1) the central source is a point source and with the conventional accretion disk for the observations at $\lambda$=3.8 and 4.8 µm, and (2) all the cases for the observations at $\lambda$=11.3 µm. The latter trend is attributed to the relatively large diffraction pattern at the wavelength of the observations (Figure 3). The coronagraph worked better for V1735 Cyg than FU Ori, in particular for the Gaussian disk models, because of the smaller angular scales at a larger distance. At $\lambda$=3.8 and 4.8 µm, for which both of the CVC and RAVC coronagraphs are available, RAVC yielded a better capability to reduce the flux of the central source by a factor of up to $\sim$5.

Figure 5 shows the on-axis PSFs for FU Ori with the CVC coronagraph for three bands and different central disks. To show the instrumental PSFs, we stacked the PSF cubes without field rotation. At $\lambda$=3.8 and 4.8 µm, the PSF for each central source consists of a bright core with six diffraction spikes at 60^∘ intervals, and a ring with a diameter of about 1 arcsecond. The six diffraction spikes are due to the shadow of the support structure that hold the secondary mirror of ELT. The ring corresponds to the control radius of the SCAO with the ELT-M4 deformable mirror (0.8, 1.0, and 2.5 arcsec for $\lambda$=3.8, 4.8, and 11.3 $\micron$, respectively).

We find marginal differences in the PSFs between the smallest and largest disks shown in the left and right ends of the figure. When the central source is larger, the core and the spikes are brighter while the brightness of the ring remains the same. The brightest PSFs at $\lambda$=3.8 and 4.8 µm are also associated with six more faint diffraction spikes between the six bright spikes. The PSFs at $\lambda$=11.3 µm resemble the central part of the PSFs at $\lambda$=3.8 and 4.8 µm. The PSFs for the RAVC coronagraph also show the same trends described above (see Appendix C).

Figure 6 shows the same PSFs but after subtracting the reference PSF and applying the field rotation before stacking (B2 and B3 in Section 2.3) in order to show how the residual of PSF subtraction can affect the observations of the extended emission. The residual pattern of the PSF subtraction does not show significant differences between different central disks but does in the absolute level of the speckle noise associated with the spike patterns, and a negative ring occurs due to the flux scaling of the reference PSF. The use of the RAVC coronagraph yields a similar trend (Appendix C). For both cases, the spike patterns become significantly less noticeable because of (1) the PSF subtraction; and (2) the field derotation of the individual frames before stacking. For FU Ori and V1735 Cyg, the major features of our interest in the extended emission are located within the negative ring pattern described above.

Figure 7 shows the 5-$\sigma$ detection limits of the extended emission for various cases of the CVC observations as a function of angular distance from the central disk. The calculation for each dot was made by measuring the root mean square of the speckle and thermal noise in a 21$\times$21-pixel box (approximately 0$\farcs$1$\times$0$\farcs$1; see Table 1) between the residual of the bright spikes after the PSF subtraction. Before the measurements, we convolved the image using a 2-D Gaussian with a FWHM of 30 mas to improve the detection limit without significantly degrading the images for the target extended emission (Section 3.2).

At longer wavelengths ($\lambda$=4.8 and 11.3 µm), the constant noise level in the outer radii is due to photon noise. Any other curves, whose spatial variations are due to the speckle noise, show that the detection limits are smaller for larger distances at the inner radii, but these increase at the outer radii due to speckle noise associated with the negative rings (Figure 6).

In Figure 8, we compare the detection limits for the CVC and RAVC coronagraphs using the smallest and largest central disks. CVC yields better detection limits than RAVC. Therefore, we use CVC for the rest of the paper in order to investigate the detection of the target extended emission for various cases.

Figures 9 and 10 show PSF-subtracted images for four combinations of the extended emission (a disk or an envelope), the central radiation source at NIR wavelengths (a disk or a star), and dust grain models (Dust1, 2, and 3). These combinations are tabulated in Table 5 as Cases 1-4. Cases 1 and 2 are the brightest cases where we use the extended disk and the envelope, respectively. Cases 3 and 4 are the same but for the faintest cases (see Figure 2). For all cases, we used the conventional accretion disk as the bright central source in the MIR. For reference, we also measured the signal-to-noise at the peak positions indicated in the figures. These values are tabulated in Table 4.

In general, the observations at shorter wavelengths yielded better signal-to-noise ratios due to the brighter nature of the target extended emission (Figure 2) and significantly fainter thermal background (Table 1). In contrast, the target extended emission can be observed only for the brightest cases at $\lambda$=11.3 µm.

The use of the Gaussian disks instead of the conventional accretion disk yields similar images to those in Figures 9 and 10, (but with comparable or lower signal-to-noise ratios) for the following reasons: (1) the major emission features shown in these figures are located within the negative ring pattern caused by the PSF subtraction; and (2) these central compact disks yield similar noise pattern but the absolute intensity scale within the ring (see Section 3.1). In Figure 11, we plot the brightness of the target extended emission and 5-$\sigma$ detection limits at the emission peaks indicated in Figures 9 and 10. If the bright central source in the MIR is the conventional accretion disk, the observations at $\lambda$=3.8 µm will allow us to detect extended emission over 5-$\sigma$ for many cases (8-10 out of 12 for each of FU Ori and V1735 Cyg). At $\lambda$=4.8 µm, the detections will be limited to 4-9 cases for each target. At $\lambda$=11.3 µm, the extended emission will be detected only for the brightest cases. These are: when the central radiation source at NIR wavelengths is a disk, and the extended emission is due to a disk with Dust3.

The detection rates are worse for larger central disks. If it is a Gaussian disk with a HWHM of 4 au, the extended emission will be detected only at $\lambda$=3.8 and 4.8 µm, if the central radiation source at NIR wavelengths is a disk, and the extended emission is associated with an extended disk.

FU Ori is known to be associated with a binary companion (FU Ori S) at a separation of 0$\farcs$5 (Wang et al., 2004; Reipurth & Aspin, 2004). Reipurth & Aspin (2004) measured $L^{\prime}$-band ($\lambda=3.8$ µm) magnitudes of FU Ori and FU Ori S of 4.2 and 8.1, respectively. This companion was only marginally visible in the images in Section 3.2 because the imaging polarimetry technique used for the original $H$-band images significantly reduced its flux (see Figure 1).

Figure 12 shows the simulated images with the binary companion and the extended emission at $\lambda$=3.8 µm with the CVC coronagraph and for Cases 1-3. As in Figure 9, the bright central source at the observing wavelength is assumed to be the conventional disk.  As shown in Figure 9, these are the best cases for increasing the signal-to-noise ratio for the extended emission in terms of the observing wavelength (see Figure 9) and the bright central source in the MIR (Figure 11). Figure 12 still shows the bright arm in the east. However, these images suggest that the presence of a binary companion can significantly hamper observations of the extended emission if they spatially overlap.

This problem could be resolved if we use a PSF reference for the binary companion as well as the bright central emission. We need different integrations for these PSFs than those for the bright central source because the companion and the bright emission are located off-axis and on-axis of the coronagraph optics, respectively.