A Multi-Wavelength Study of the Dynamic Environment surrounding the FUor V960 Mon

V960 Mon was observed by ALMA in band 6 (2016.1.00209.S, PI: M.Takami) on 20 April 2017 (compact configuration) and on 27 July 2017 (extended configuration). The observation included in total five spectral windows: two of a bandwidth of 1.875 GHz dedicated to the continuum at central frequencies of 216.877 GHz and 232.178 GHz; and three spectral windows of a bandwidth of 0.059 GHz to cover the $J=2-1$ transitions of ¹²CO, ¹³CO and C¹⁸O centered at 230.514 GHz, 220.375 GHz and 219.537 GHz, and channel widths of 15.259 kHz, 30.158 kHz, and 30.158 kHz, respectively.

The band 6 continuum data of a reduced field of view (FOV) around V960 Mon were presented in Kóspál et al. (2021, ${\rm FOV}\sim 4^{"}$) and further analyzed in Weber et al. (2023, ${\rm FOV}\sim 2^{"}$). Here, we used the band 6 continuum reduction with natural weighting presented in Weber et al. (2023), where self-calibration was performed on the compact and extended configurations separately before combining them. For details about the data reduction, we refer the reader to that article. For the spectral windows targeting the CO molecular lines, we subtracted the continuum from the lines and cleaned each channel with the CASA task tclean. We predefined cleaning masks by only permitting signal above three times the channel’s rms-noise, but where appropriate, we added user-defined masks to avoid losing the faint signal that reappears over several channels. We further inspected the continuum channels for the detection of molecular lines (see Appendix A), and identified SiO, N₂D⁺, DCO⁺ and DCN. Because each molecular emission region has a unique spatial extent, we examined each case individually to identify baselines that produced sinusoidal artifacts from extended emission contamination and flagged them accordingly.

ALMA band 3 and band 4 data (both 2019.1.01144.S, PI: H.B. Liu) were taken on the 23 August 2021 and 22 August 2021, respectively. We reduced these data using natural weights. Both observations included four spectral windows dedicated to the continuum of a bandwidth of 1.875 GHz. They were centered at 86.00 GHz, 87.86 GHz, 98.19 GHz and 100.00 GHz for band 3 and 146.00 GHz, 147.86 GHz, 158.19 GHz and 160.00 GHz for band 4. Because the continuum emission shows a low signal-to-noise ratio (SNR) in bands 3 and 4, and no significant improvement was expected, we did not perform any self-calibration. We list the resulting dimensions of the clean beam and the SNR of each band in Table 1.

V960 Mon was observed on the night of 24 January 2021 using the Multi Unit Spectroscopic Explorer (MUSE) at the VLT (0106.C-0510, PI: F. Cruz-Saenz de Miera) with a total integration time of 25.6 min. The observing conditions were excellent (seeing: 0.7”, a coherence time of 5 ms, a wind speed of 6 m s^-1 at 30 m altitude). The spectrograph divided the light into 3721 channels of $\Delta\lambda=0.125\,$nm, ranging from $470\,$nm to $935\,$nm. To process the data, we used the mpdaf.obj-package²²2https://mpdaf.readthedocs.io/en/latest/index.html (Bacon et al., 2016, version 3.6). To generate an image with sufficient signal, we stacked the spectral channels. Our analysis of the data revealed that the signals of interest, particularly those emanating from the large-scale environment of the FUor object, are predominantly found in the longer wavelength channels. Therefore, in our stacking process, we selectively included the 1000 channels with the longest wavelengths (810–935 nm). We then binned these 1000 channels into 50 representative channels. In each representative channel, we subtracted the average of columns and rows to remove the vertical and horizontal stripes seen on the detector. A spectral analysis of the line emission will be carried out in a later study.

The $H$-band VLT/SPHERE image presented in Weber et al. (2023), focused exclusively on the polarized component of the observed light. However, the same observation (098.C-0422, PI: L. Cieza) also provides the total intensity (i.e., the sum of both beams on the detector). To obtain the total intensity, we summed all the aligned object frames. We produced the total intensity using the IRDAP reduction pipeline³³3https://irdap.readthedocs.io/en/latest/ (van Holstein et al., 2020, version 1.3.5).

V960 Mon was observed by ALMA in band 6, band 4 and band 3. In a previous study, Weber et al. (2023) reported a fragmenting spiral arm in the band 6 continuum. Their analysis focused on unresolved clumps detected at roughly the 5 $\sigma_{\rm rms}$ level. We do not find any band 3 or band 4 counterparts to these clumps. This non-detection is consistent with thermal emission from the clumps, which would yield fluxes of $\lesssim 2\,\sigma_{\rm rms}$ in band 4 and $\lesssim 1\,\sigma_{\rm rms}$ in band 3 for the rms-noise values listed in Table 1.

When inspecting the large-scale environment of the FUor, several bright sources appear in the FOV as shown in the band 6 image in the central panel of Fig. 1. Following Kóspál et al. (2021), we labeled these sources V960 Mon for the outbursting star at the center of the image, V960 Mon N for the northern object, V960 Mon E for the eastern object, V960 Mon SE for the southeastern object and V960 Mon S for the southern object. The sources, their projected position relative to V960 Mon, their continuum flux in the different bands, the corresponding estimated dust masses and spectral indices are listed in Table 2. To convert flux into mass, we assumed a dust temperature of 20 K and an opacity of $\kappa_{\nu}=1\,{\rm cm}^{2}\,{\rm g}^{-1}\,(\nu/100\,{\rm GHz})$.

The spectral windows of the ALMA band 6 observation cover several molecular transition lines. For line detection, we extracted spectra directly from the visibilities of each spectral window, identifying signals corresponding to rotational transitions of ¹²CO, ¹³CO, C¹⁸O, SiO, DCO⁺, N₂D⁺ and DCN. This approach is effective for detecting the brightest emission lines, as visibilities inherently integrate spatial information including noise. To achieve more robust detections and isolate emission from individual sources, we subsequently performed the analysis on the reconstructed data cubes for each line. These molecules, the detected rotational transitions, their reference frequency, the velocity resolution of the observation in the respective spectral window, and the upper energy level of the transitions are listed in Table 3.

The left panel of Fig. 6 displays the image obtained by stacking the MUSE channels as described in Sec. 2.2. Across the entire MUSE detector field of view ($65\farcs 6\times 67\farcs 6$), we detected 24 point sources. Within the area covered by the ALMA band 6 primary beam (${\rm FWHM}\sim 26^{\prime\prime}$), only V960 Mon and V960 Mon N among the band 6 continuum sources have a clear MUSE counterpart. In addition, we find a faint point source in the northeast of V960 Mon, which we label V960 Mon NE.

In the vicinity of the primary object (at separations $\lesssim 3^{\prime\prime}$), V960 Mon’s luminosity dominates at these shorter wavelengths. However, the MUSE image reveals intriguing large-scale structures. Notably, the source identified as V960 Mon SE appears as the nucleus of a diffuse clump or cloudlet, intriguingly connected to the primary system via a structure of intertwined gas. Remarkably, this represents the first detection of material associated with V960 Mon SE at wavelengths below 4.5 µm, seen in reflection rather than in emission. Previously, V960 Mon SE was reported as undetected at these shorter wavelengths (Kóspál et al., 2015). The observed decrease in brightness with increasing distance from V960 Mon suggests that the structure is primarily illuminated by the central FUor object, rather than being self-luminous at these wavelengths.

Meanwhile, the enigmatic object to the east, V960 Mon E, shows no visible counterpart in the MUSE image. One plausible explanation is that material closer to the primary star casts a shadow over V960 Mon E, preventing direct illumination from the FUor object due to the observed geometry.

To the west of V960 Mon, the image reveals what appears to be the illuminated interior of a surrounding envelope. This observation corroborates the subtle emission noted in the ¹³CO moment-0 map, as shown in the central column of Fig. 2. While the molecular lines trace emission, the MUSE counterpart likely shows light reflected from mircon-sized grains entrailed in the structure. The observed structure is present in a large wavelength range, excluding the possibility of being associated with distinct molecular line emission in the optical and NIR.

Finally, the MUSE image in Fig. 6 underscores the value of stacked MUSE data in probing the interactions between young stars and their large-scale environments.

Polarized intensity ($I_{\rm pol}$) is a measure for imaging circumstellar environments because initially unpolarized light becomes polarized after scattering off dust grains. Consequently, observations of polarized light trace the surfaces of dusty environments. However, this technique has several limitations. First, by design, polarized light observations do not capture light directly emitted from companions (which is generally unpolarized), nor do they capture light scattered from unresolved dust structures around the light source, such as that expected from a circumstellar or circumplanetary disk. However, if this disk is sufficiently inclined or very asymmetric, the observed light from the central unresolved beam could be polarized. Therefore, this method is not suitable for direct imaging of companions and can only detect them if their emission becomes net-polarized within the telescope’s point spread function along the light path (see e.g., van Holstein et al., 2021, for a discussion). Furthermore, the assumption that scattered light contributes to the polarized intensity relies on single scattering of photons. If the majority of stellar photons in a certain region undergo multiple scatterings before escaping the system toward the observer, the emergent light becomes unpolarized due to the summation over randomly oriented polarization states.

Consequently, we analyzed the total intensity in which photons from companions and multiple scattering events still contribute, even if at a contrast that is difficult to detect. The central panel of Fig. 6 displays the VLT/SPHERE $H$-band total intensity image, where the flux is scaled by the squared distance from the center to emphasize large-scale features. A bright elongated structure appears to the south of V960 Mon. Additionally, in the right panel of Fig. 6, we show the polarized light $H$-band image (presented in Weber et al., 2023) in black and white. This comparison shows that the elongated total intensity feature disappears approximately at the distances to where the spiral structures observed in scattered light extend. Notably, we do not detect any polarization from the $H$-band counterpart, even after applying the same $r^{2}$-scaling. Here, we discuss the reason for why this might be. Even if the assumption of single-scattering holds true, the efficiency of polarization in a scattering event still depends on the scattering angle $\theta_{\rm scat}$. In particular, polarization is predicted to be suppressed for backward scattering ($\theta_{\rm scat}>90^{\circ}$) (e.g., Min et al., 2005). Observations of multiple disks confirm a decrease in polarization efficiency at large scattering angles (e.g., Ginski et al., 2023).

Consequently, the lack of polarized light from the elongated feature in V960 Mon indicates a low polarization efficiency, likely due to backward scattering of the primary’s light. This, in turn, implies that the elongated feature is situated behind the star from our viewpoint. In fact, the feature disappears closer to the star, where the spiral structure is visible in scattered light, a result consistent with the spirals’ optically thick region blocking the line of sight.

Finally, we also note that although V960 Mon E lies within the SPHERE field of view, as in the MUSE image, it is detected in neither the total intensity nor the polarized intensity at $H$-band. We conclude that V960 Mon E is deeply embedded, causing its flux to be heavily attenuated at short wavelengths.

Assuming the Gaia-measured distance to V960 Mon is accurate, the VLT/SPHERE image in scattered light (Weber et al., 2023) shows material present around V960 Mon up to a separation of $\gtrsim 2000\,$au. How such a vast environment can be sustained is an intriguing question. As suggested by the detection of the elongated feature (see Fig. 6, Sec. 3.4), one possibility is that V960 Mon is currently being fed by asymmetric infall from the neighboring SE source, separated by about 8$\farcs$4. While infalling streamers onto circumstellar environments have been increasingly reported in recent years (e.g., Ginski et al., 2021; Garufi et al., 2022; Cacciapuoti et al., 2024), this could mark an interesting case where the associated reservoir of the streamer’s origin is observed. From the optical and NIR evidence alone, it is not discernible whether the feature connecting V960 Mon SE and V960 Mon represents infall onto the FUor. A similar feature can also be observed during a flyby in which the gravitational interaction between the two objects creates an over-dense bridge (e.g. Cuello et al., 2020), or from colliding flows in dense star formation regions (e.g. Kuffmeier et al., 2019).

Guided by the detection in the optical and NIR, we investigated the ALMA band 6 cube for molecular counterparts. From a channel-by-channel analysis, we found faint ¹³CO signal between 21.17 km s^-1 and 22.36 km s^-1 co-located with the MUSE elongated feature. This implies that the line of sight velocity in the connecting gas is blue-shifted by $2-3$ km s^-1 with respect to V960 Mon. We combine the channels in this velocity range of the cube using bettermoments, masking out signal below 2$\,\sigma_{\rm rms}$ in each of the 14 channels. The resulting moment-0 map limited to this velocity range has an rms-noise of $\sim 3$ Jy beam^-1 m s^-1.

Fig. 7 shows the MUSE image superimposed by the emission contours from these combined ¹³CO channels at 6 Jy beam^-1 m s^-1 (corresponding to $\sim 2\sigma_{\rm rms}$). In addition to the emission indicating the connection between V960 Mon and V960 Mon SE, there is faint ¹³CO emission at these velocities west of the primary, which could trace the inner wall of the surrounding envelope.

Having concluded in Sec. 3.4 that the elongated structure is likely back-grounded with respect to V960 Mon, the blue-shifted velocity of the gas in the feature indicates that it is moving towards the outbursting star V960 Mon. Hence, the combined analysis of optical, NIR and ¹³CO data gives evidence that the observed feature indeed represents infall onto the primary.

Kuffmeier et al. (2018) demonstrated that in a circumstellar environment the infall of material can locally reduce the Toomre parameter, leading to a localized GI and the formation of spiral structures in the disk around the star. However, the simulation results presented there focused on an area much closer to the star. The very extended structure around V960 Mon raises the question of whether the spirals are part of a gravitationally unstable disk or if they arise from infalling material that becomes twisted as it moves inward. Testing a connection between the spirals and disk instability requires measuring underlying Keplerian rotation, evidence that enhanced ¹³CO observations could provide. Moreover, although Kuffmeier et al. (2018) showed that localized GI induced by infall can generate a spiral structure, no fragmentation was observed in their simulations. However, Longarini et al. (2023) suggest that in certain scenarios, fragmentation occurs only in the dust phase, a finding consistent with the ALMA continuum observation and not at odds with Kuffmeier et al. (2018), where dust was not implemented. If the spiral structure arises from non-Keplerian motion, it remains an open question under what conditions the infalling material might fragment.

GI in the outer parts of a protoplanetary disk can generate significant angular momentum transport. This may lead to an increased gas flow towards the inner disk region where the switching on and off of the magnetorotational instability may trigger outbursts (Armitage et al., 2001; Zhu et al., 2009; Stamatellos et al., 2011). If the disk gas is replenished via infall, outbursts could also be initiated due to GI-enhanced accretion (Kuffmeier et al., 2018) or due to clump formation by fragmentation and subsequent accretion (Vorobyov & Basu, 2015).

Originally identified by Kóspál et al. (2021), the continuum source V960 Mon E was re-detected in the band 3 and band 4 continuum data, yet remains undetected in the optical and NIR. This implies that V960 Mon E is deeply embedded in dust that is optically thick at these wavelengths. As discussed in Section 3.2, a substantial portion of the ¹²CO emission lies in the foreground of V960 Mon and appears spatially anchored near V960 Mon E. Fig. 8 shows the ¹²CO moment-1 map where regions of negative moment-0 are masked. The velocity field is centered on V960 Mon E’s velocity, with blue- and red-shifted lobes on opposite sides of its continuum emission. The inset of Fig. 8 highlights the spatial alignment between the SiO jet and the blue-shifted cone of ¹²CO emission. This alignment supports the interpretation of the observed ¹²CO feature as part of an outflow driven by V960 Mon E. From ¹³CO we measure a central line of sight velocity of approximately $23.19\pm 0.03$ km s^-1(see Appendix B). We note that this value is subject to large systematic errors from contamination within the region that is not included in the presented uncertainty, so the resulting velocity estimate should be treated with caution. This velocity is blue-shifted by about $1$ km s^-1 with respect to V960 Mon. This measurement sets as a lower limit for their relative motion and suggests that V960 Mon E is above the escape velocity of the system, hence probably not on a gravitationally bound orbit around V960 Mon. Still, the close spatial and spectral proximity between V960 Mon E to V960 Mon suggests a dynamical association, further supported by a tentative connecting structure, faintly suggested in both C¹⁸O and N₂D⁺ moment-0 maps.

These observations raise the question of the dynamic history of the object driving the ¹²CO outflow. The detection of C¹⁸O, DCO⁺, and N₂D⁺ is commonly interpreted as indicative of a protostellar envelope (Tychoniec et al., 2021). Thus, V960 Mon E may represent a protostar forming from the same molecular reservoir that V960 Mon is accreting from (similar to scenarios described by e.g. Kuffmeier et al., 2019), or it could be a stellar flyby. In fact, the object V960 Mon SE presents another viable flyby candidate. The line of sight velocity of V960 Mon SE is measured to be 23.55 km s^-1 from ¹²CO, so blue-shifted by 0.7 km s^-1 with respect to V960 Mon and likely above the escape velocity of the system. Flybys have been shown to produce spiral features in circumstellar dust and gas and connecting structures between the flyby components (e.g., Cuello et al., 2020), suggesting that such a scenario may be at play here. Interestingly, recent work by Fu et al. (2024) discussed the possibility that such a connecting feature induced by a flyby may gravitationally fragment to form future free-floating giant planets.

In terms of stellar accretion behavior, several studies pointed out that flybys may generate periods of enhanced accretion rates on the interacting stars, causing FUor events (e.g., Forgan & Rice, 2010; Vorobyov et al., 2021; Borchert et al., 2022). However, the simulated flybys produced sharply peaked bursts that are very unlike the observed extended classic FUor outbursts (cf. section 5.4 in Nayakshin & Elbakyan, 2024, for a discussion). To determine whether V960 Mon E or V960 Mon SE are bound to or unbound from the FUor, more accurate line of sight velocity measurements are needed. In addition, a deeper examination of the surrounding environment may reveal connecting structures that can inform on their dynamical relationship.

Finally, we want to discuss the scenario that some of the observed mm-sources represent compact objects ejected from the circumstellar environment of V960 Mon due to a stellar intruder, and the connection of this with the outbursting state of V960 Mon.

Recent work indicated that disk Thermal Instability (Bell & Lin, 1994) in conjunction with thermal evaporation of massive planets formed by GI in the outer parts of protoplanetary disks presents an attractive scenario for FU Ori itself (Nayakshin et al., 2023), and potentially for several other observed FUor objects (Nayakshin & Elbakyan, 2024; Nayakshin et al., 2024). The key to this scenario is a massive protoplanetary disk that is able to both hatch the fragments at $\sim 50$ au and to also drive them inward into the inner $\sim 0.1$ au rapidly enough. We note that GI fragment interactions with other fragments were shown to occasionally speed up the process of fragment migration from the outer into the inner disk (Cha & Nayakshin, 2011; Vorobyov & Elbakyan, 2018). The mechanics of this process is simple. The two bodies exchange angular momentum, with one of them gaining and the other losing angular momentum in a close interaction, the former thus gets onto an eccentric orbit with a much smaller peri-center distance, whereas the latter gets kicked to a larger orbit or ejected. The more massive the interaction partner, the more pronouced is the dynamic effect. In this regard, Vorobyov et al. (2017) have shown that close passages of stellar intruders may unbind GI-fragments from the disk, scattering them towards the star or ejecting them from the system. Importantly, these fragments may be massive enough to be in the high-mass brown dwarf and even low-mass star regime (cf. their Table 1). Previous 2D and 3D GI disk simulations showed prolific formation of such fragments (e.g. Vorobyov & Basu, 2006, 2010; Cha & Nayakshin, 2011), but were unable to resolve the innermost region where FUor outbursts originate. Nayakshin et al. (2023) showed that there in the inner disk Thermal Instability creates such a hot environment that these fragments may evaporate in a steady-state manner, providing an explanation for the bewilderingly long duration of FUor outbursts.

There is a potential connection of this scenario for FUor outbursts to the very disturbed and surprisingly large-scale gaseous environment of V960 Mon. There are several candidates for an intruder star in the immediate environment of V960 Mon. Caratti o Garatti et al. (2015) reported a stellar object at a separation of about $227\,$mas from the primary, confirmed as a bright source in SPHERE $H$-band in Weber et al. (2023). Further, VLT/ERIS observations in $L^{\prime}$-band ($\lambda_{\rm obs}=3.8\,$µm) detected a candidate companion at a distance of $0\farcs 9$ to the southwest of V960 Mon, roughly co-located with the fragmented spiral arm observed in the band 6 continuum (Dasgupta et al. in prep.). A potential intruder star may carry with it a part of the disk of the primary; and the ejected objects also tend to be surrounded by a gaseous disk. As a result, shortly after passage, the system looks very disturbed, with one or more gaseous filaments connecting the primary with the intruder and the ejectees (cf. Figs. 2 and 3 in Vorobyov et al., 2017). The structure naturally has a scale much larger than that of a typical protoplanetary disk, extending well outside a thousand au.

Finally, even if the interaction between fragment and intruder failed to be sufficiently strong, flybys may still encourage faster GI-fragment destruction via more rapid disk migration. As 3D simulations of disk flybys show (Forgan & Rice, 2010; Vorobyov et al., 2021; Borchert et al., 2022), the inner $\sim 10$ au of disks get sufficiently perturbed by the interaction to lead to higher accretion rates onto the primary star. Although this effect itself is unlikely to produce a steady FUor-outburst-type light curve, as we argued above, the increased gas accretion rate in the inner $\sim 10$ au may promote much more rapid fragment migration, which could then evaporate in the inner disk as in the scenario of Nayakshin et al. (2023).

Thus, an interaction with an intruder could simultaneously explain the formation of the southeastern elongated feature that is connecting V960 Mon SE and V960 Mon, and V960 Mon’s FUor outburst. This scenario could be further tested by future observations of the system. In particular, constraints on the grain size in the southeastern elongated structure could help to differentiate between an infalling filament and the flyby “bridge” scenarios, as the grains would be larger in the latter case.