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

# A Multiwavelength Study of the Dynamic Environment Surrounding the FUor V960 Mon

Philipp Weber<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Silvio Ulloa<fnref n="2" />,<fnref n="4" /> ID, Sebastián Pérez<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, James Miley<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Lucas Cieza<fnref n="2" />,<fnref n="5" /> ID, Sergei Nayakshin<fnref n="6" /> ID, Alice Zurlo<fnref n="2" />,<fnref n="5" />,<fnref n="7" /> ID, Hauyu Baobab Liu<fnref n="8" />,<fnref n="9" /> ID, Fernando Cruz-Sáenz de Miera<fnref n="10" /> ID, Antonio Hales<fnref n="11" /> ID, Antonio Garufi<fnref n="12" /> ID, Dimitris Stamatellos<fnref n="13" /> ID, Ágnes Kóspál<fnref n="14" />,<fnref n="15" />,<fnref n="16" /> ID, and Viviana Guzmán<fnref n="2" />,<fnref n="17" /> ID

## Abstract

The evolution of young stars and planet-forming environments is intrinsically linked to their nascent surroundings. This is particularly evident for FU Orionis (FUor) objects—a class of young protostars known for dramatic outbursts resulting in significant increases in brightness. We present a case study of V960 Mon, an FUor that has recently been found to show signs of a fragmenting spiral arm, potentially connected to planet formation. Our study explores the large-scale environment ($10^3–10^4$ au) and incorporates Atacama Large Millimeter/submillimeter Array band 3, band 4, and band 6 continuum data, molecular emissions from $^{12}$CO, $^{13}$CO, $^{C^{18}O}$, SiO, DCO$^+$, $N_2D^+$, and DCN, alongside optical and near-infrared observations from Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) and VLT Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE). We map a region of 20′′ across where we find tantalizing emissions that provide a unique view of a young group of protostars, including the discovery of a class-0 protostar to the east of the FUor. The $^{12}$CO and SiO tracers suggest that this object is at the base of an outflow, potentially impacting the surrounding medium. The MUSE and SPHERE observations indicate the presence of an elongated feature towards a prominent source to the southeast that may represent interaction between V960 Mon and its surroundings. Moreover, the $^{C^{18}O}$ emission overlaps with the clumps of the detected fragmenting spiral arm. These findings provide the strongest evidence to date for a connection between infalling material, fragmentation, and the intensity outburst of a protostar. Our case study highlights the complex interactions between young stars and their surroundings that drive the evolution of the planet forming environment.

Unified Astronomy Thesaurus concepts: FU Orionis stars (553); Observational astronomy (1145); Planet formation (1241); Submillimeter astronomy (1647); Young stellar objects (1834)

## 1. Introduction

FU Orionis (FUor) events represent a transient phase in early stellar evolution, characterized by sudden and dramatic increases in luminosity. These luminous outbursts, typically elevating the star's brightness by several magnitudes, are believed to result from episodic accretion events, during which a significant amount of material from the circumstellar disk is rapidly dumped onto the star (L. Hartmann & S. J. Kenyon 1996; M. Audard et al. 2014). Such events not only alter the star's appearance but also have profound implications for the evolution of the surrounding

protoplanetary disk. It has been suggested that these episodic accretion events could play a pivotal role in the early stages of planet formation.

First, the occurrence of an accretion burst has been suggested to be directly linked to planet-forming processes such as gravitational instability (GI; P. J. Armitage et al. 2001; Z. Zhu et al. 2009; D. Stamatellos et al. 2011).<fnref n="18" /> This association makes FUor objects prime targets for searching for signs of planet formation via GI. An important question in this context is how the high disk mass required to trigger GI can be sustained.

Second, the feedback from the enhanced accretion rates and consequent disk heating during FUor outbursts may influence the coagulation and composition of dust grains, potentially

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Departamento de Física, Universidad de Santiago de Chile, Avenida Victor Jara 3659, Santiago, Chile; philippweber@gmail.com</footnote>

<footnote n="2">² Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Chile</footnote>

<footnote n="3">³ Center for Interdisciplinary Research in Astrophysics and Space Exploration (CIRAS), Universidad de Santiago de Chile, Chile</footnote>

<footnote n="4">⁴ Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="5">⁵ Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Avenida Ejercito 441, Santiago, Chile</footnote>

<footnote n="6">⁶ School of Physics and Astronomy, University of Leicester, Leicester, LE1 7RH, UK</footnote>

<footnote n="7">⁷ Escuela de Ingeniería Industrial, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Avenida Ejercito 441, Santiago, Chile</footnote>

<footnote n="8">⁸ Department of Physics, National Sun Yat-Sen University, No. 70, Lien-Hai Road, Kaohsiung City 80424, Taiwan, R.O.C.</footnote>

<footnote n="9">⁹ Center of Astronomy and Gravitation, National Taiwan Normal University, Taipei 116, Taiwan</footnote>

<footnote n="10">¹⁰ Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UT3-PS, OMP, CNRS, France</footnote>

<footnote n="11">¹¹ 8 National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA</footnote>

<footnote n="12">¹² INAF—Istituto di Radioastronomia, Via Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="13">¹³ Jeremiah Horrocks Institute for Mathematics, Physics and Astronomy, University of Central Lancashire, Preston, PR1 2HE, UK</footnote>

<footnote n="14">¹⁴ Konkoly Observatory, HUN-REN Research Centre for Astronomy and Earth Sciences, MTA Centre of Excellence, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="15">¹⁵ Institute of Physics and Astronomy, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="16">¹⁶ Max-PLanck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="17">¹⁷ Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
Received 2025 February 25; revised 2025 April 1; accepted 2025 April 5; published 2025 May 27</footnote>

<footnote n="18">¹⁸ For a broader discussion on the origin of FUor events, see the review by W. J. Fischer et al. (2023).</footnote>
accelerating the formation of planetesimals (e.g., Z. Zhu et al. 2009). As such, understanding FUor phenomena could provide valuable insights into the dynamic processes underpinning the birth and early evolution of planetary systems.

Building on these theoretical foundations, observational studies have significantly advanced our understanding of FUor phenomena and their implications for planet formation. Early observational efforts primarily focused on documenting the outburst characteristics, such as the rapid rise in brightness and spectral evolution (G. H. Herbig 1977; L. Hartmann & S. J. Kenyon 1985). More recently, high-resolution observations in the (sub)millimeter and polarimetric imaging in the near-infrared (NIR) have provided more information on the FUors' surroundings. Some of these observations have revealed intricate details of the circumstellar environment during and after FUor outbursts, including changes in disk structure, temperature profiles, and dust grain properties (L. A. Cieza et al. 2016; A. Zurlo et al. 2017; M. Takami et al. 2018; A. Zurlo et al. 2024). The rapid temperature rise around FUor objects also reveals their chemical composition by evaporating disk compounds that would otherwise remain in solid form. This has facilitated the study of disks' chemistry, in particular the detection and study of complex organic molecules (J.-E. Lee et al. 2019, J.-H. Jeong et al. 2025, F. Cruz-Sáenz de Miera et al. 2025) and water vapor (J. J. Tobin et al. 2023).

Recent numerical simulations support the notion that the origin of an FUor event can trace back to the extensive environment of the young star. For example, M. Kuffmeier et al. (2018) found that the material that falls from the molecular cloud onto the circumstellar disk could locally induce GI. Meanwhile, GI facilitates the rapid transport of angular momentum and a subsequent accretion burst. Similarly, C. P. Dullemond et al. (2019) described a scenario in which material is captured from an adjacent cloudlet, triggering an accretion burst.

In a broader context, the structure and composition of the disk have been found to depend on the interaction with the environment, which can affect the processes of stellar evolution and planet formation (J. E. Pineda et al. 2023; A. J. Winter et al. 2024).

V960 Mon has received significant interest in the astronomical community due to its characteristic outburst properties and potential insights into early stellar evolution. First identified in its outburst state in 2014 (L. Hillenbrand 2014; H. Maehara et al. 2014), V960 Mon presents an exceptional opportunity to study the mechanisms and consequences of FUor-type eruptions in real time. According to Gaia DR3, the system is located at a distance of $2189 \pm 281$ pc (Gaia Collaboration et al. 2023). We employ this value with caution because a high RUWE indicates possible systematic uncertainties, but it remains our distance estimate throughout this study. At a distance of $\gtrsim 2$ kpc—which prevents resolving substructures around individual objects—this system offers an exceptional view into the accretion dynamics of young stellar systems and their interaction with the protostellar envelope.

V960 Mon poses a spectacular circumstellar environment, which putatively revealed the first evidence of GI-fragments of planetary masses. P. Weber et al. (2023) presented a polarimetric image in the NIR obtained with the Very Large Telescope's (VLT's) SPHERE instrument which shows characteristic large-scale spiral arms. A complementing Atacama

Large Millimeter/submillimeter Array (ALMA) observation detected that one of those spiral arms seems to be undergoing fragmentation, indicated by several clumps of millimeter-emission that align along a spiral arm. Here, we further study the case of V960 Mon, turning our attention to the environment at an even larger scale (~$10^4$ au), presenting continuum and molecular line emission in the millimeter regime and scattered light in the optical and NIR.

In Section 2, we describe the details of the observations and data reduction, we present the main results in Section 3, and discuss their implications in Section 4. In Section 5, we summarize our findings.

## 2. Observations and Data Reduction

### 2.1. ALMA

V960 Mon was observed by ALMA in band 6 (2016.1.00209. S, PI: M.Takami) on 2017 April 20 (compact configuration) and on 2017 July 27 (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 $^{12}$CO, $^{13}$CO and $^{18}$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, FOV ~ 4") and further analyzed in P. Weber et al. (2023, FOV ~ 2"). Here, we used the band 6 continuum reduction with natural weighting presented in P. 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_2D^+$, 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 2021 August 23 and 2021 August 22, 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.
<table><caption>**Table 1**
List of Archival ALMA Observations</caption><thead><tr><td>ALMA Band</td><td>Beam</td><td>σ<sub>rms</sub></td><td>Project-ID</td></tr></thead><tbody><tr><td>Band 6</td><td>0''20 × 0''14</td><td>28</td><td>2016.1.00209.S</td></tr><tr><td>Band 4</td><td>0''07 × 0''05</td><td>25</td><td>2019.1.01144.S</td></tr><tr><td>Band 3</td><td>0''11 × 0''08</td><td>23</td><td>2019.1.01144.S</td></tr></tbody></table>

Note. The measured rms-noise, σ<sub>rms</sub>, is given in µJy beam<sup>-1</sup>.

## 2.2. VLT/MUSE

V960 Mon was observed on the night of 2021 January 24 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 minutes. The observing conditions were excellent (seeing: 0.7, a coherence time of 5 ms, a wind speed of 6 m s<sup>-1</sup> at 30 m altitude). The spectrograph divided the light into 3721 channels of Δλ = 0.125 nm, ranging from 470 to 935 nm. To process the data, we used the `mpdaf.obj-package`<fnref n="19" /> (R. 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.

## 2.3. VLT/SPHERE Total Intensity

The H-band VLT/SPHERE image presented in P. 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<fnref n="20" /> (R. G. van Holstein et al. 2020, version 1.3.5).

# 3. Results

## 3.1. ALMA Continuum Emission

V960 Mon was observed by ALMA in band 6, band 4, and band 3. In a previous study, P. 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σ<sub>rms</sub> level. We do not find any band 3 or band 4 counterparts to these clumps. This nondetection is consistent with thermal emission from the clumps, which would yield fluxes of ≲2σ<sub>rms</sub> in band 4 and ≲1σ<sub>rms</sub> 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 Figure 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 κ<sub>ν</sub> = 1 cm<sup>2</sup> g<sup>-1</sup> (ν/100 GHz).

*Sources and Limits.* V960 Mon N has been detected several times (see e.g., Á. Kóspál et al. 2015) because it is also present in the optical and NIR. We detected it in all three analyzed continuum bands.

V960 Mon E was previously reported from the same ALMA band 6 data set in Á. Kóspál et al. (2021), but was not confirmed in any other observations. Here, we found that it also shows significant emission in ALMA bands 3 and 4, as shown in the left panel of Figure 1, where the contours show the emission at 5σ<sub>rms</sub> for the different bands. The peak emission reaches 7σ<sub>rms</sub>, 14σ<sub>rms</sub>, and 40σ<sub>rms</sub> in bands 3, 4, and 6, corresponding to 0.161 mJy beam<sup>-1</sup>, 0.35 mJy beam<sup>-1</sup>, and 1.12 mJy beam<sup>-1</sup>.

The source V960 Mon SE is even more luminous. It was first reported in Á. Kóspál et al. (2015), where the authors showed that the compact emission does not appear for observational wavelengths λ<sub>obs</sub> ≲ 4.5 µm but already dominates the flux of the system at λ<sub>obs</sub> = 70 µm. The right panel of Figure 1 shows the 5σ<sub>rms</sub> contours for the three different ALMA bands, centered on V960 Mon SE. In contrast to V960 Mon E, this source appears resolved in all three observations, at least along the north–south axis. This clearly shows that the emission is elongated along this axis, ruling out emission from a single point source. To estimate the structure’s spatial extent, we fitted a two-dimensional Gaussian to the band 4 data. The fit yielded full width at half-maximum (FWHM) values for a major axis of 0''30 ± 0''02 and a minor axis of 0''085 ± 0''01. If this continuum emission stems from a disk-like configuration, the corresponding disk radius would be ~280 au, with an inclination of i<sub>SE</sub> ≥ 74° ± 2°, and a position angle of -0.3 ± 1°.

Finally, we detected the unresolved source V960 Mon S in band 6 and band 4, while band 3 does not yield a significant detection.

*Spectral Index.* To investigate the nature of the detected emission, we calculated the spectral indices for the sources in the FOV, for different band combinations, following the formula:

$$ \alpha_{i/j} = \frac{\log F_i/F_j}{\log \nu_i/\nu_j}, \qquad (1) $$

where F is the observed flux from an object in two distinct bands i and j, and ν is the corresponding representative frequency. We calculated the representative frequencies by averaging the respective values of the continuum spectral windows and find ν<sub>3</sub> = 93 GHz, ν<sub>4</sub> = 153 GHz and ν<sub>6</sub> = 224.5 GHz. Table 2 lists the spectral indices calculated from different bands for the five most prominent flux sources in the band 6 FOV. The uncertainties given for the spectral indices in Table 2 stem from error propagation of the uncertainties in the flux measurement due to the rms-noise. We add a note of caution—while including band 6 in the spectral index calculation is advantageous—thanks to its higher S/N and greater

<footnote n="19">¹⁹ https://mpdaf.readthedocs.io/en/latest/index.html</footnote>

<footnote n="20">²⁰ https://irdap.readthedocs.io/en/latest/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Continuum emission from the region around V960 Mon, in units of mJy beam⁻¹. Center: A large-scale view of the band 6 continuum around V960 Mon. Prominent sources of millimeter-emission are marked by arrows and labeled. Left: Detected emission named V960 Mon E. The image itself shows the band 6 data, white, red, and blue contours are at 5σrms for band 6, band 4, and band 3, respectively. The ellipses in the bottom left corner show the restored beams for the respective color. Right: Same as left panel but for V960 Mon SE.</figcaption></figure>

<table><caption>**Table 2**
List of millimeter-sources Visible in ALMA Data</caption><thead><tr><th>Name</th><th><i>R</i><br>(")</th><th>P.A.<br>(°)</th><th><i>F</i><sub>3</sub><br>(mJy)</th><th><i>F</i><sub>4</sub><br>(mJy)</th><th><i>F</i><sub>6</sub><br>(mJy)</th><th><i>M</i><sub>3</sub><br>(M<sub>⊕</sub>)</th><th><i>M</i><sub>4</sub><br>(M<sub>⊕</sub>)</th><th><i>M</i><sub>6</sub><br>(M<sub>⊕</sub>)</th><th>α<sub>6/4</sub></th><th>α<sub>4/3</sub></th><th>α<sub>6/3</sub></th></tr></thead><tbody><tr><td>V960 Mon</td><td>...</td><td>...</td><td>0.045</td><td>0.36</td><td>0.94</td><td>40</td><td>77</td><td>71</td><td>1.9 ± 0.2</td><td>2.3 ± 0.5</td><td>2.1 ± 0.3</td></tr><tr><td>N</td><td>5.54</td><td>14.7</td><td>0.17</td><td>0.42</td><td>1.09</td><td>151</td><td>90</td><td>82</td><td>3.1 ± 0.2</td><td>1.7 ± 0.4</td><td>2.3 ± 0.2</td></tr><tr><td>E</td><td>5.32</td><td>96.6</td><td>0.12</td><td>0.42</td><td>1.66</td><td>107</td><td>92</td><td>125</td><td>2.8 ± 0.2</td><td>1.6 ± 0.2</td><td>2.1 ± 0.1</td></tr><tr><td>SE</td><td>8.46</td><td>152.4</td><td>0.62</td><td>2.49</td><td>6.41</td><td>558</td><td>540</td><td>482</td><td>1.7 ± 0.4</td><td>2.7 ± 0.5</td><td>2.1 ± 0.4</td></tr><tr><td>S</td><td>9.05</td><td>197.7</td><td>...</td><td>0.20</td><td>0.82</td><td>...</td><td>42</td><td>62</td><td>2.1 ± 0.4</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** The different columns list the distance to the primary, *R*, and the corresponding position angle P.A., as well as the integrated fluxes measured for band 3, band 4, and band 6, and the corresponding dust mass estimates, assuming optically thin emission at a temperature of 20 K and a dust opacity of $κ_ν = 1 \text{ cm}^2 \text{ g}^{-1}$ (ν/100 GHz). The spectral index $α$ is defined in Equation (1). The estimated errors arise from the propagated rms-noise of each band listed in Table 1.

frequency separation, it may also introduce a systematic error. The band 6 observations were carried out approximately four years before the band 4 and band 3 data. Since V960 Mon (and possibly other detected sources) appears variable, a variability on annual timescales would significantly compromise the spectral index analysis. Band 3 and band 4 data were taken on subsequent days.

In source S, the only measurement of the spectral index, $α_{6/4}$, is consistent with 2. In the remaining four sources, the values of $α_{6/3}$ are also consistent with 2. These sources are likely optically thick at 225 GHz with optical depth $τ_{225 \text{ GHz}} \gtrsim 5$, similar to the bright Class-II disks in the Taurus–Auriga region (C.-Y. Chung et al. 2024), which is expected by L. Delussu et al. (2024). We briefly discuss the better-detected sources—V960 Mon, N, E, and SE—in the following.

**V960 Mon, Source SE.** The values of $α_{6/4}$ are comparable to or slightly smaller than 2. In the latter case, the maximum dust grain sizes, $a_{max}$, may be close to ~100 µm, causing the spectra to appear anomalously reddened due to dust scattering (H. B. Liu 2019). The values of $α_{4/3}$ are larger than $α_{6/4}$. This may imply that the optical depths of the dominant (sub) millimeter emitters become ~1–3 at the central frequency of our ALMA band 4 observations, ~153 GHz, although uncertainties remain depending on the assumed dust temperatures,

$a_{max}$, and dust opacity tables (see Section 5.2.1 of C.-Y. Chung et al. 2024).

**Source N, E.** The values of $α_{4/3}$ are likely somewhat smaller than 2, while the values of $α_{6/4}$ are close to 2.5–3. These are similar to what was observed in the archetypal FUor, FU Ori (H. B. Liu et al. 2019, 2021), which were also resolved in many class 0/I young stellar objects (YSOs; J. I.-H. Li et al. 2017). There are two probable interpretations: (1) Dust emission is very optically thick at >93 GHz frequencies, and the value of $a_{max}$ in the dominant (sub)millimeter emitters is a few times 100 µm (depending on the assumed dust composition and porosity). In such cases, the dust albedo can increase around ~100 GHz and then decrease at higher frequencies (see Figure 3 in H. B. Liu 2019). This leads to anomalously reddened and bluened spectra at ~100 GHz and $\gtrsim 200$ GHz, respectively (for an example, see the SED model for DM Tau in H. B. Liu et al. 2024). (2) Dust emission is very optically thick, and there is a temperature gradient in the line of sight and the warmer emitters are obscured by the colder ones. This configuration may be realized for edge-on sources (R. Galván-Madrid et al. 2018), or for sources that are dominantly heated by dissipative processes in the gaseous disks (e.g., viscous or shock heating, etc.; H. B. Liu et al. 2019, 2021; J. Zamponi et al. 2021). These two interpretations are not mutually exclusive. To discern the values of $a_{max}$, future spatially resolved observations at more frequency bands, or the (sub)millimeter linear polarization
<table><caption>**Table 3**

Detected Molecular Transitions</caption><thead><tr><th>Molecule</th><th>J</th><th>ν</th><th>ν<sub>res</sub></th><th>E<sub>up</sub></th><th>μ<sup>2</sup>S</th><th>log<sub>10</sub>(A<sub>ij</sub>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO</td><td>2–1</td><td>230.538</td><td>0.04</td><td>16.6</td><td>0.024</td><td>−6.2</td></tr><tr><td><sup>13</sup>CO</td><td>2–1</td><td>220.399</td><td>0.08</td><td>15.9</td><td>0.024</td><td>−6.2</td></tr><tr><td>C<sup>18</sup>O</td><td>2–1</td><td>219.560</td><td>0.08</td><td>15.8</td><td>0.024</td><td>−6.2</td></tr><tr><td>DCN</td><td>3–2</td><td>217.238</td><td>0.67</td><td>20.9</td><td>80.5</td><td>−3.3</td></tr><tr><td>DCO<sup>+</sup></td><td>3–2</td><td>216.113</td><td>0.68</td><td>20.7</td><td>142.0</td><td>−2.6</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td>3–2</td><td>231.322</td><td>0.63</td><td>22.2</td><td>312.1</td><td>−3.15</td></tr><tr><td>SiO</td><td>5–4</td><td>217.105</td><td>0.67</td><td>31.2</td><td>48.0</td><td>−3.3</td></tr></tbody></table>

Note. Rest frequency is in GHz, velocity resolution in km s⁻¹, upper energy level in K, and line strength (μ²S) in Debye². Aᵢⱼ is the Einstein coefficient for the respective transition.

observations to probe signatures of dust self-scattering, are necessary. The observed, lower α₄/₃ in these two sources may also be partly explained by the contribution of optically thin free-free emission. Nevertheless, previous JVLA surveys toward the Corona Australis, Ophiuchus, and Orion molecular clouds showed that free-free emission in low-mass YSOs (without including proplyds around OB stars) is commonly not bright enough to contribute significantly at >90 GHz frequencies (H. B. Liu et al. 2014; A. Coutens et al. 2019; L. Tychoniec et al. 2021).

## 3.2. Molecular Line Emission

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 because 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.

### 3.2.1. CO Isotopologues

We measured the rms-noise per respective channel width of σ<sub>rms</sub>(¹²CO) = 9.3 mJy beam⁻¹, σ<sub>rms</sub>(¹³CO) = 7.0 mJy beam⁻¹, and σ<sub>rms</sub>(C¹⁸O) = 6.6 mJy beam⁻¹. We collapsed the CO data cubes and calculated the moment-0 (integrated intensity) and moment-1 (intensity-weighted mean velocity) maps using the bettermoments package (R. Teague & D. Foreman-Mackey 2018). For both moments, we applied sigma clipping below 2σ<sub>rms</sub> and smoothed the threshold mask over one beam size to avoid steep edges in the mask. Figure 2 shows the resulting moment-0 maps of the CO isotopes at different spatial scales and moment-1 maps for ¹²CO and ¹³CO. From these data, we inferred the line-of-sight velocities listed in Table 4 and described in Appendix B.

¹²CO. The left column of Figure 2 presents the ¹²CO isotopologue, where significant foreground emission obscures the region surrounding V960 Mon. We estimated the rms-noise of the moment-0 map to be σ<sub>rms</sub>(¹²CO) = 20 Jy beam⁻¹ m s⁻¹. Examination of the ¹²CO kinematics in the left bottom panel

reveals blue- and redshifted lobes spectrally centered around a velocity of approximately 22.3 km s⁻¹, and spatially centered on the millimeter-continuum source designated as V960 Mon E. This pattern suggests a projected outflow velocity of about 3.5 km s⁻¹, originating from the unclassified object V960 Mon E. We further discuss this feature in Section 3.2.2 and Section 4.

When integrating over the entire spectral axis, multiple emission components blend along the line of sight. Figure 3 shows the emission of ¹²CO recorded at 23.12 km s⁻¹ (we stacked channels from 22.92 km s⁻¹ to 23.32 km s⁻¹ to increase the SNR). At these velocities, the signal shows gas around the primary, while also including what seems to be interaction with slightly blueshifted gas reservoirs to the northwest and south. For comparison, the black contours show the polarized light of the SPHERE H-band image (P. Weber et al. 2023), which delineate a spiral structure that is possibly related to the large-scale interaction.

¹³CO. The central panels of Figure 2 show the moments of the optically thinner ¹³CO emission. We estimated the noise of the moment-0 map to be σ<sub>rms</sub>(¹³CO) = 7 Jy beam⁻¹ m s⁻¹. Within the FOV, we identified three principal regions of ¹³CO emission: in the vicinity of the primary source V960 Mon, near V960 Mon SE, and around V960 Mon S. Notably, all three regions coincide with band 6 continuum emission at levels above 10σ<sub>rms</sub>. The bottom central panel of Figure 2 shows the moment-1 map of the ¹³CO line, where a blue/red-pattern centered on the continuum emission of V960 Mon SE becomes apparent, suggesting rotation around a massive compact source. Around V960 Mon, we found redshifted contributions toward the northeast, while the clumps in the southwest are colocated with blueshifted line emission.

C¹⁸O. The right column of Figure 2 displays the moment-0 map of detected C¹⁸O emission around V960 Mon. The measured noise of the map is σ<sub>rms</sub>(C¹⁸O) = 3 Jy beam⁻¹ m s⁻¹. The top right panel shows the region including V960 Mon and the E source. Emission from the V960 Mon E is detected with a maximum of 23.2 Jy beam⁻¹ m s⁻¹. The bottom right panel focuses on the area closer to V960 Mon, coinciding with the fragmented continuum spiral presented in P. Weber et al. (2023). The C¹⁸O emission at the clumps' location reaches a maximum value of 25.0 Jy beam⁻¹ m s⁻¹. The detection of C¹⁸O suggests a substantial gas content associated with the clumps, indicating that the fragmentation is not confined solely to the dust phase.

### 3.2.2. Further Species

In Figure 4, we show the tracers, SiO, DCO⁺, N₂D⁺, and DCN, that we detected in the spectral windows that were dedicated to the continuum emission. To produce the moment-0 map, we used bettermoments without any sigma clipping or smoothing. The following paragraphs discuss the distribution and significance of each tracer. The deuterated species are thought to be a product of deuterium fractionation that distributes deuterium from HD into other species (H. Nomura et al. 2023). The deuterium fraction is expected to be significantly enhanced in cold areas (T ≲ 30 K), the detection of deuterium can therefore inform about the local temperature. Nevertheless, the different deuterated species are detected at very different locations in the image, necessitating investigation of each molecule individually. As seen from Table 3, DCO⁺, N₂D⁺, and DCN all have similar upper energy states (E<sub>up</sub> ≈ 20 K). This implies that differences in the distribution of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** CO($J=2-1$) line emission in Jy beam$^{-1}$ m s$^{-1}$ for the three different isotopes $^{12}$CO (left), $^{13}$CO (center), and $^{18}$O (right) at different spatial scales. The top left and center panels show the moment-0 maps for $^{12}$CO and $^{13}$CO, respectively, with the bottom panels showing the corresponding moment-1 map. The crosses indicate the locations of the band 6 continuum emission. The bottom center panel includes insets with black contours at $5\sigma$ and $10\sigma$ of the band 6 continuum emission. The color scale of the insets ranges from the displayed velocity to $\pm 1$ km s$^{-1}$. The right column shows the $^{18}$O moment-0 map at different scales, with the lower panel focusing on the area around V960 Mon. The white contours show the band 6 continuum at $5\sigma$ and the white ellipse in the bottom left corner marks the clean beam of the observation.</figcaption></figure>

<table><caption>**Table 4**
Line-of-sight Velocities for different Objects and the Employed Molecule (see Appendix B)</caption><thead><tr><th>Object</th><th>v<sub>los</sub><br>(km s<sup>-1</sup>)</th><th>Molecule</th></tr></thead><tbody><tr><td>V960 Mon</td><td>24.25</td><td><sup>12</sup>CO</td></tr><tr><td>V960 Mon E</td><td>23.19</td><td><sup>13</sup>CO</td></tr><tr><td>V960 Mon SE</td><td>23.55</td><td><sup>12</sup>CO</td></tr><tr><td>V960 Mon S</td><td>26.38</td><td><sup>12</sup>CO</td></tr></tbody></table>

their emission arise primarily due to different local abundances and line strengths, rather than excitation temperatures.

SiO. The left panel of Figure 4 displays the moment-0 map of Silicon monoxide (SiO) with an rms-noise of $\sigma_{rms}(\text{SiO}) \approx 8.8$ Jy beam$^{-1}$ m s$^{-1}$. The image shows a strong and collimated emission feature attached to V960 Mon E. The peak emission of the SiO feature is shifted by $-14.8$ km s$^{-1}$ with respect to the restframe, so about $-37.1$ km s$^{-1}$ relative to the systematic velocity of V960 Mon E ($v_E = 22.3$ km s$^{-1}$). Such collimated SiO emission is typically linked to a high velocity jet shocking with the surrounding medium. This is characteristic for actively accreting class-0 stars (e.g., C. Codella et al. 2007). While the SiO feature is spatially very confined, it is widely spread over different velocity channels, with a typical FWHM line width of about $10$ km s$^{-1}$ (from

$-30$ km s$^{-1}$ to $-40$ km s$^{-1}$ with respect to V960 Mon E). These broad emission lines are typical of shocked regions (e.g., R. Bachiller & M. Pérez Gutiérrez 1997). For Figure 5, we rotated the jet by $-25.5$ to align it horizontally. In the three upper panels of Figure 5, we show the peak intensity, moment-1 and moment-2 (intensity-weighted velocity dispersion, $\sigma_v$), overlayed with the peak intensity contours at 10, 20, and $40\sigma_{rms}$. The white contours display that the width of the jet is roughly constant with separation. We fit a Gaussian function to the jet width and retrieve a FWHM of $0.21$, similar to the band 6 beam size. For the two lower panels of Figure 5, we extract the maximum values (minimum, for moment-1) along the vertical direction in the top panels. To reduce the spiky structure of the velocity profiles, we smoothed the distribution with a Gaussian kernel. The intensity profile highlights that the emission is distributed in several knots along the jet. Similar behavior has been detected in SiO jets for several other class-0 objects (e.g., C. Codella et al. 2007; L. Podio et al. 2016; L. Y. Wang et al. 2019). Here, we measure a periodic reoccurrence of the knots with a distance of about $0.55$ between consecutive maxima and a clear detection of at least seven peaks. In the bottom panel of Figure 5 we show that the relative line-of-sight velocity to V960 Mon E becomes more blueshifted with distance, ranging from about $-30$ km s$^{-1}$ at the base of the jet, to about $-40$ km s$^{-1}$ at distances larger than $1''$. However, this effect may not represent a genuine
<figure><img src="image_5.png" /><figcaption>**Figure 3.** $^{12}$CO emission at 23.12 km s$^{-1}$in Jy beam$^{-1}$ m s$^{-1}$. The black contours delineate the spiral structure observed in NIR polarized light (P. Weber et al. 2023). The white plus signs indicate the sources of millimeter-emission in the FOV.</figcaption></figure>

acceleration but could instead reflect jet rotation within 1". By contrast, the velocity dispersion is greatest near to V960 Mon E ($\sigma_v \sim 7$ km s$^{-1}$) and declines with distance. To estimate the time interval between jet expulsions from the angular separations measured in the image plane, we assume that the jet's velocity in the plane of the sky is comparable to its line-of-sight velocity ($\sim 40$ km s$^{-1}$) and that all the knots have this same constant velocity. From this, we find that the knots correspond to expulsion events roughly every 150 yr, consistent with the periodicity found for molecular jets in other systems (see C.-F. Lee 2020, for a review). This variability implies that V960 Mon E is episodically accreting.

DCO$^{+}$. We detected the deuterated formylium ion (DCO$^{+}$) in the vicinity of V960 Mon E. The moment-0 map presented in the center left panel of Figure 4 has an rms-noise of $\sigma_{rms}(DCO^{+}) \approx 3.0$ Jy beam$^{-1}$ m s$^{-1}$. Notably, the DCO$^{+}$ emission is spatially offset from the continuum source by 0.29, corresponding to a projected distance of approximately 644 au. The presence of DCO$^{+}$ requires some amount of gas-phase CO and low temperatures. The asymmetry of the DCO$^{+}$ emission near V960 Mon E might result from low SNR or indicate an intrinsically nonuniform environment, possibly shaped by outflows from V960 Mon E itself and by the complex gravitational interactions with the surrounding system. We further detect some DCO$^{+}$ emission from around V960 Mon S, which is shown in Appendix A.

$N_2D^+$. The center right panel of Figure 4 shows the moment-0 map of $N_2D^+$ with an rms-noise of $\sigma_{rms}(SiO) \approx 5.4$ Jy beam$^{-1}$ m s$^{-1}$. The ion $N_2D^+$ is destroyed by reactions with CO, inhibiting its detection in CO-rich regions (e.g., P. Caselli et al. 2002). Hence, $N_2D^+$ emission is strongly anticorrelated with CO. Such conditions typically arise in cold environments where CO freeze-out onto dust grains prevails,

supporting both the presence of $N_2D^+$ and its deuterium enrichment. In contrast to the more spatially localized DCO$^+$ emission, the detected $N_2D^+$ distribution is relatively diffuse, appearing as sparse emission to the northeast of V960 Mon E and in the region between V960 Mon E and V960 Mon. The emission is only marginally above $3\sigma_{rms}$ in four channels and primarily appears in the channel centered at 23.72 km s$^{-1}$ (see Appendix A). Notably, the detected $N_2D^+$ signal locally overlaps with faint continuum band 6 emission, suggesting the presence of a potential bridging structure. Such filamentary connections have been discussed in the context of multiple stellar-formation scenarios (e.g., M. Kuffmeier et al. 2019; C. Gieser et al. 2024) and in cases involving stellar flybys (e.g., N. Cuello et al. 2020; F. Ménard et al. 2020).

DCN. The right panel of Figure 4 shows the moment-0 map of DCN with an rms-noise of $\sigma_{rms}(DCN) \approx 4.1$ Jy beam$^{-1}$ m s$^{-1}$. The DCN emission is colocated with portions of the fragmented continuum spiral arm around V960 Mon. The detection is based on four adjacent spectral channels between 23.03 and 25.05 km s$^{-1}$, with the maximum emission at 23.70 km s$^{-1}$ (see Appendix A). In the right panel of Figure 4, we overlay the moment-0 contours of the $C^{18}$O transition on the DCN emission. This comparison highlights a clear anticorrelation: DCN is more prominent in the outer regions of the spiral, whereas it is depleted toward the inner clumps, where $C^{18}$O is strongest. The presence of DCN in the outer clumps implies a substantial gas reservoir at these locations, and the $C^{18}$O detection farther inward suggests that gas extends along the entire fragmented continuum spiral. Furthermore, the detection of DCN-alongside the nondetections of DCO$^+$ and $N_2D^+$-implies that this region has recently evolved from a very cold environment to one that is slightly warmer. The ionized deuterated species such as DCO$^+$ and $N_2D^+$ are more easily depleted than DCN once the temperature rises (e.g., S. D. Rodgers & T. J. Millar 1996; Y. Aikawa et al. 2005; H. B. Liu et al. 2015). Therefore, the observed DCN emission, combined with the lack of DCO$^+$ and $N_2D^+$, supports a scenario in which the local chemistry is transitioning from a cold to a mildly warmer state. This is consistent with the recent FUor event, which is expected to heat the surrounding region through irradiation. However, a temperature increase could also be connected to contraction of the related clumps.

As stated in L. Tychoniec et al. (2021), the protostellar envelope is well traced by $C^{18}$O, DCO$^+$, and $N_2D^+$, all of which are detected around V960 Mon E. The different locations and distributions at which we detect the signal indicate warmer and colder regions, characterized by the freeze-out of CO. Additionally, it is noteworthy that we do not detect any significant $^{13}$CO or $C^{18}$O emission above the noise level centered on the outbursting star V960 Mon. This might be due to large-scale contamination in the image that is filtering out the emission from the primary. The deuterated lines originate in smaller regions, making them less affected by this. The northern companion, V960 Mon N, shows no trace of any molecular tracer, nor does it show any connections to other objects in the FOV. This raises the question if V960 Mon N is truly associated with V960 Mon, or whether it is a foreground or background star.

### 3.3. MUSE Imaging

The left panel of Figure 6 displays the image obtained by stacking the MUSE channels as described in Section 2.2. Across the entire MUSE detector FOV (65.6 × 67.6), we
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Moment-0 maps of molecular tracers in the system of V960 Mon. The units are mJy beam⁻¹ m s⁻¹. The panels show the area between V960 Mon E and V960 Mon (there is no relevant detection outside this area). In all panels, the white contours trace the band 6 continuum flux at 5 and 10 σ<sub>rms</sub>. In the inset of the right panel, the blue contours mark the C<sup>18</sup>O at 5σ<sub>rms</sub>. In both the insets of the center left and right panels, the orange ellipse in the top left corner shows the clean beam of the observation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Analysis of the SiO jet rotated by 25°.5. The top three panels show the distributions of the flux at peak emission, the moment-1 and moment-2, respectively. They share the same colorbar of linear scale with minimum to maximum ranges of [0 to 34] mJy beam⁻¹, [−45 to −26] km s⁻¹ and [0 to 7] km s⁻¹, respectively. In all three panels, the contours trace *I*<sub>peak</sub> at 10 and 20σ<sub>rms</sub>. For a quantitative analysis, we plot the maximum intensity, the minimum velocity relative to the systematic velocity of V960 Mon E and the maximum velocity dispersion, σ<sub>v</sub>, in the lower two panels.</figcaption></figure>

detected 24 point sources. Within the area covered by the ALMA band 6 primary beam (FWHM ~ 26"), 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 ≲3"), 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 Figure 2. While the molecular lines trace emission, the MUSE counterpart likely shows light reflected from micron-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 Figure 6 underscores the value of stacked MUSE data in probing the interactions between young stars and their large-scale environments.

### 3.4. SPHERE Total Intensity

Polarized intensity (*I*<sub>pol</sub>) 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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panel: VLT/MUSE image of V960 Mon and the region around it. Center panel: VLT/SPHERE *H*-band total intensity scaled with $r^2$ to highlight faint signal at large scales. Note the reduced FOV due to the smaller detector of SPHERE. Right panel: SPHERE total intensity image overlaid with SPHERE polarized intensity image (grayscale, see P. Weber et al. 2023).</figcaption></figure>

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., R. 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 Figure 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 Figure 6, we show the polarized light *H*-band image (presented in P. 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_{scat}$. In particular, polarization is predicted to be suppressed for backward scattering ($\theta_{scat} > 90^\circ$) (e.g., M. Min et al. 2005). Observations of multiple disks confirm a decrease in polarization efficiency at large scattering angles (e.g., C. 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 FOV, 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.

## 4. Discussion

The different observations of V960 Mon reveal a dynamic system around an FUor star, including several companions, spiral structures, clumps, and a connection to the large-scale environment. The concurrence of these effects paints a complicated picture, where it remains unclear which process causes which effect. This section considers several scenarios that connect the presented observations with concepts of star and planet formation. In particular, we will discuss the possibilities of large-scale infall perturbing the circumstellar environment (Section 4.1), a possible stellar flyby (Section 4.2), and the ejection of compact objects from around V960 Mon by a stellar intruder to the circumstellar environment (Section 4.3).

### 4.1. Fragmentation and Outburst Induced by Infall?

Assuming the Gaia-measured distance to V960 Mon is accurate, the VLT/SPHERE image in scattered light (P. 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 Figure 6, Section 3.4), one possibility is that V960 Mon is currently being fed by asymmetric infall from the neighboring SE source, separated by about 8"/4. While infalling streamers onto circumstellar environments have been increasingly reported in recent years (e.g., C. Ginski et al. 2021; A. Garufi et al. 2022; L. 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 or not. 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., N. Cuello et al. 2020), or from colliding flows in dense star formation regions (e.g., M. Kuffmeier et al. 2019).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** VLT/MUSE image (same as left panel of Figure 6). The blue contours show the ¹³CO emission integrated between 21.17 and 22.36 km s⁻¹ (blueshifted by ~2–3 km s⁻¹ toward V960 Mon), where emission from the elongated feature can partly be separated from the contamination of the cloud at 2σrms.</figcaption></figure>

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⁻¹ and 22.36 km s⁻¹ colocated with the MUSE elongated feature. This implies that the line-of-sight velocity in the connecting gas is blueshifted by 2–3 km s⁻¹ with respect to V960 Mon. We combine the channels in this velocity range of the cube using bettermoments, masking out signal below 2σrms in each of the 14 channels. The resulting moment-0 map limited to this velocity range has an rms-noise of ~3 Jy beam⁻¹ m s⁻¹. Figure 7 shows the MUSE image superimposed by the emission contours from these combined ¹³CO channels at 6 Jy beam⁻¹ m s⁻¹ (corresponding to ~2σ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 Section 3.4 that the elongated structure is likely back-grounded with respect to V960 Mon, the blueshifted velocity of the gas in the feature indicates that it is moving toward 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.

M. 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 M. Kuffmeier et al. (2018) showed that localized GI induced by infall can generate a spiral structure, no fragmentation was observed in their simulations. However, C. 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 M. 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 toward the inner disk region where the switching on and off of the magnetorotational instability may trigger outbursts (P. J. Armitage et al. 2001; Z. Zhu et al. 2009; D. Stamatellos et al. 2011). If the disk gas is replenished via infall, outbursts could also be initiated due to GI-enhanced accretion (M. Kuffmeier et al. 2018) or due to clump formation by fragmentation and subsequent accretion (E. I. Vorobyov & S. Basu 2015).

## 4.2. A Possible Stellar Flyby?

Originally identified by Á. Kóspál et al. (2021), the continuum source V960 Mon E was redetected 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. Figure 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 redshifted lobes on opposite sides of its continuum emission. The inset of Figure 8 highlights the spatial alignment between the SiO jet and the blueshifted 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 ± 0.03 km s⁻¹ (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 blueshifted by about 1 km s⁻¹ 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 (L. 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., M. 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** $^{12}$CO moment-1 map, limited to where the $^{12}$CO moment-0 is positive. The velocity map is centered on the systematic velocity inferred for V960 Mon E, $v_E = 23.19$ km s$^{-1}$. The zoom-in is centered on the enigmatic object V960 Mon E, with the black contours showing the band 6 continuum emission and the orange contour showing the SiO line emission, highlighting the alignment with the blueshifted cone.</figcaption></figure>

velocity of V960 Mon SE is measured to be 23.55 km s$^{-1}$ from $^{12}$CO, so blueshifted 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., N. Cuello et al. 2020), suggesting that such a scenario may be at play here. Interestingly, recent work by Z. Fu et al. (2025) 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., D. Forgan & K. Rice 2010; E. I. Vorobyov et al. 2021; E. M. A. Borchert et al. 2022). However, the simulated flybys produced sharply peaked bursts that are very unlike the observed extended classic FUor outbursts (see Section 5.4 in S. Nayakshin & V. 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.

### 4.3. Ejection from Circumstellar Environment?

Finally, we want to discuss the scenario that some of the observed millimeter-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 (K. R. Bell & D. N. C. 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 (S. Nayakshin et al. 2023), and potentially for several other observed FUor objects (S. Nayakshin & V. Elbakyan 2024; S. Nayakshin et al. 2024). The key to this scenario is a massive protoplanetary disk that is able to both hatch the fragments at ~50 au and to also drive them inward into the inner ~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 (S.-H. Cha & S. Nayakshin 2011; E. I. Vorobyov & V. G. Elbakyan 2018). The mechanics of this process are 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 pronounced is the dynamic effect. In this regard, E. I. Vorobyov et al. (2017) have shown that close passages of stellar intruders may unbind GI-fragments from the disk, scattering them toward 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 (see their Table 1). Previous 2D and 3D GI disk simulations showed prolific formation of such fragments (e.g., E. I. Vorobyov & S. Basu 2006, 2010; S.-H. Cha & S. Nayakshin 2011), but were unable to resolve the innermost region where FUor outbursts originate. S. Nayakshin et al. (2023) showed that 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. A. Caratti & 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 P. Weber et al. (2023). Further, VLT/ERIS observations in L'-band ($\lambda_{obs}$ = 3.8 µm) detected a candidate companion at a distance of 0.99 to the southwest of V960 Mon, roughly colocated with the fragmented spiral arm observed in the band 6 continuum (A. Dasgupta et al. 2025, in preparation). 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 (see Figures 2 and 3 in E. I. 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 (D. Forgan & K. Rice 2010; E. I. Vorobyov et al. 2021; E. M. A. Borchert et al. 2022), the inner ~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 ~10 au may promote much more rapid fragment migration, which could then evaporate in the inner disk, as in the scenario of S. 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 because the grains would be larger in the latter case.

## 5. Conclusions

We presented ALMA, VLT/MUSE, and VLT/SPHERE observations of the outbursting star V960 Mon. Our analysis places the 2014 FUor outburst (H. Maehara et al. 2014) and the fragmenting spiral arm reported by P. Weber et al. (2023) into the broader environmental context:

1. We identified multiple compact objects emitting at millimeter wavelengths and found signs of direct interaction with V960 Mon for at least two of them, V960 Mon SE and V960 Mon E.

2. From the millimeter spectral index, we infer that the ALMA continuum signal of those compact objects traces thermal emission emanating from the center of dense gas cores.

3. The fragmenting spiral arm is not detected in band 4 or band 3, also consistent with thermal emission at given sensitivities.

4. V960 Mon E shows characteristics of a class-0 protostar and drives an outflow composed of a slower component traced by $^{12}$CO and a faster, more collimated jet identified through SiO emission.

5. V960 Mon SE is the most prominent in millimeter continuum emission and we detect rotation in $^{12}$CO and $^{13}$CO line data, indicating a compact object at its center.

6. MUSE and SPHERE imaging show a connection between V960 Mon SE and V960 Mon. We infer that the structure is situated behind V960 Mon, and from counterparts in $^{13}$CO channels we find that it is blueshifted—the material is moving toward V960 Mon, likely representing infall onto the system.

7. This directly raises the question how this infall is connected to the outbursting state of the FUor V960 Mon, and the large-scale spirals observed in the NIR and the fragmenting spiral arm observed by ALMA (P. Weber et al. 2023).

These multiwavelength observations suggest that V960 Mon is in a unique and transient evolutionary phase. The system displays dynamic disruption across multiple scales-ranging from eruptive accretion onto the stellar surface, to gravitationally unstable spiral arms in dust and gas, and extending to large-scale (≥10,000 au) gravitational interactions with neighboring protostars. This raises key questions about how processes at different scales interconnect and whether the observed scenario is representative of the general star and planet formation paradigm.

## Acknowledgments

We thank the referee Ruobing Dong for a constructive report that helped to improve the quality of this work. P.W. acknowledges support from FONDECYT grant 3220399. This work was funded by ANID—Millennium Science Initiative Program—Center Codes NCN2021_080 and NCN2024_001. S.P. acknowledges support from FONDECYT Regular grant 1231663. J.M. acknowledges support from FONDECYT de Postdoctorado 2024 #3240612. H.B.L. is supported by the National Science and Technology Council (NSTC) of Taiwan (grant Nos. 111-2112-M-110-022-MY3, 113-2112-M-110-022-MY3). The research activities described in this paper were carried out with contribution of the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4—Education and Research, Component 2—From Research to Business (M4C2), Investment Line 3.1—Strengthening and creation of Research Infrastructures, Project IR0000034—"STILES— Strengthening the Italian Leadership in ELT and SKA." 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). This work was also supported by the NKFIH NKKP grant ADVANCED 149943 and the NKFIH excellence grant TKP2021-NKTA-64. 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. This work made use of the Puelche cluster hosted at CIRAS/USACH. Based on observations collected at the European Southern Observatory under ESO programmes 098.C-0422(B) and 0106.C-0510(A). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00209.S and ADS/JAO.ALMA#2019.1.01144.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

**Facilities:** ALMA, VLT:Melipal, VLT: Yepun.

**Software:** This work has made use of the IRDAP-pipeline (R. G. van Holstein et al. 2020) for the processing of SPHERE/IRDIS data, the CASA software (CASA Team et al. 2022) for the processing of the ALMA data and the `mpdaf.obj-package` (R. Bacon et al. 2016) for creating the MUSE cube. We used `bettermoments` (R. Teague & D. Foreman-Mackey 2018) for creating moment-0 and moment-1 maps for the molecular line emission in the ALMA data. We used `IPython` (F. Perez & B. E. Granger 2007), `NumPy` (S. van der Walt et al. 2011) and `Matplotlib` (J. D. Hunter 2007) for data analysis and creating figures.

## Appendix A: ALMA Spectrum

We examined the band 6 continuum channels to identify any notable molecular emission. We focused on the continuum
<figure><img src="image_11.png" /><figcaption>**Figure 9.** ALMA band 6 spectrum for four sources in the FOV. The vertical lines mark the rest emission frequencies of detected molecular transitions.</figcaption></figure>

sources in the FOV, placing a circular aperture mask of radius 1′′0 on each source. By integrating the flux in each channel, we derived the spectra shown in Figure 9 (before continuum subtraction).

The spectrum of V960 Mon is dominated by DCN. The right panel of Figure 4 shows that this contribution stems from the fragmenting spiral arm inside the aperture centered on V960 Mon. In V960 Mon E, three lines (DCO$^{+}$, SiO, and N$_{2}$D$^{+}$) are clearly present. V960 Mon SE displays no significant detections, although strong negative values near the DCO$^{+}$ line are attributed to contaminated baselines. For V960 Mon S, DCO$^{+}$ is the only pronounced feature.

Figure 10 illustrates the individual channels containing notable flux for DCN, DCO$^{+}$, and N$_{2}$D$^{+}$. The FOVs vary to

highlight different sources. Both DCN and DCO$^{+}$ channels reach an rms-noise level of $\sigma_{rms} = 1.8$ mJy beam$^{-1}$, while N$_{2}$D$^{+}$ has $\sigma_{rms} = 2.5$ mJy beam$^{-1}$. In the top row of Figure 10, DCN emission exceeding $3\sigma_{rms}$ appears only around the band 6 continuum clumps, outlined by white contours. In the central row (DCO$^{+}$), emission is strongest around V960 Mon E and V960 Mon S, with tentative flux in the area between V960 Mon E and V960 Mon that appears to connect to the fragmented continuum spiral arm. This interpretation is further supported by N$_{2}$D$^{+}$ in the bottom row, which shows a similar alignment in those channels. In both cases (DCO$^{+}$ and N$_{2}$D$^{+}$), the emission from the area between V960 Mon and V960 Mon E is strongest for the channel that lies between the two objects systematic line-of-sight velocities.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** We show the intensity in units of mJy beam⁻¹ maps for the channels containing significant signal for the molecules DCN (top row), DCO⁺ (center) and N₂D⁺ (bottom row). For each of those molecular lines we found contributions limited to four successive channels. In each image, the white contours show the band 6 continuum flux at 5 and 10σ<sub>rms</sub>. The millimeter-sources in the FOV are marked with white labels.</figcaption></figure>

**Appendix B**
**Line-of-sight Velocities**

We aim to measure the line-of-sight velocities for the compact objects in the FOV from the Doppler shifts of CO molecular lines. The contamination and low SNR of the optically thinner molecules complicate precise velocity measurements around the continuum sources. We estimate the systematic velocity of V960 Mon, V960 Mon E, V960 Mon SE, and V960 Mon S. V960 Mon N does not show any associated molecular line emission, so is not included in this analysis.

For V960 Mon, the spectrum was extracted using a 1" aperture on the position of the star. We calculate the rms-noise of the spectrum from the aperture-averaged data. The four panels of Figure 11 show the spectra in dependence on the systematic velocity for four different objects. The shaded area depicts the rms-noise between −σ<sub>rms</sub> and +σ<sub>rms</sub>, for V960 Mon, V960 Mon SE, and V960 Mon S corresponding to <sup>12</sup>CO, for V960 Mon E corresponding to <sup>13</sup>CO. For V960 Mon,

for all three species, we fit a Gaussian to the wings of the emission profile. The center of these Gaussians is very similar for the different isotopologues, with $\bar{v}_{12} = (24.25 \pm 0.03) \text{ km s}^{-1}$, $\bar{v}_{13} = (24.16 \pm 0.03) \text{ km s}^{-1}$, and $\bar{v}_{18} = (24.12 \pm 0.09) \text{ km s}^{-1}$, for <sup>12</sup>CO, <sup>13</sup>CO, and C<sup>18</sup>O, respectively. We choose the value of the <sup>12</sup>CO measurement as a representative value, due to the higher SNR and the better velocity resolution. The double-peak feature for <sup>13</sup>CO and C<sup>18</sup>O probably indicates rotation around the primary. However, we cannot rule out that the diminution at the systematic velocity is due to self-absorption of the molecular line.

For V960 Mon E, the <sup>12</sup>CO emission is strongly attenuated. The SNR for <sup>13</sup>CO and C<sup>18</sup>O is low, but we detect a peak at the object's location (using a mask radius of 1/200 mas) that is spectrally roughly consistent between the two isotopologues. The fit to the data points is shown in the second panel of Figure 11, with values of $\bar{v}_{13}^E = (23.19 \pm 0.03) \text{ km s}^{-1}$ and $\bar{v}_{18}^E = (23.34 \pm 0.04) \text{ km s}^{-1}$. We note that, especially at higher velocities, the CO that is comoving with the primary seems to contaminate the spectrum.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** CO spectra from 1" around the four sources V960 Mon, E, SE, and S, from the tracers $^{12}$CO, $^{13}$CO, and C$^{18}$O (where applicable). The circles show the integrated flux in each channel, the lines correspond to Gaussian fits. We mark the selected representative line-of-sight velocity by a vertical line. The vertical shaded area corresponds to the expected rms-noise of the spectrum. This approximated by the value for the aperture-averaged flux from the $^{12}$CO, except for V960 Mon E where we use the value inferred from $^{13}$CO.</figcaption></figure>

The third panel of Figure 11 shows the $^{12}$CO and $^{13}$CO spectrum for V960 Mon SE. Here, we limit the fit to the area in which we detect significant band 6 continuum emission $I_6 \ge 0.22$ mJy beam$^{-1}$. We do this because we detect rotation centered on the continuum in this area. From visual inspection of the $^{13}$CO moment-1 map (shown in the inset of the lower central panel of Figure 2) that shows rotation around V960 Mon SE, we expect the central velocity of this rotation between 23 and 24 km s$^{-1}$. In addition, the emission around V960 Mon SE seems to be strongly attenuated at these velocities. We, therefore, fit the Gaussian only to the wings, which seem to be less affected by the contamination. Same as for V960 Mon, we fit a Gaussian to the wings of the emission from around V960 Mon SE and retrieve velocity centers of $\bar{v}_{12}^{SE} = (23.55 \pm 0.07)$ km s$^{-1}$ and $\bar{v}_{13}^{SE} = (23.84 \pm 0.14)$ km s$^{-1}$.

Finally, the right-most channel of Figure 11 shows the spectrum from V960 Mon S using a mask of 0"/5. The $^{13}$CO emission shows a low SNR, the $^{12}$CO emission shows two peaks, one blueshifted and one redshifted with respect to V960 Mon. From inspecting the channels, we find that the blueshifted peak at around 18 km s$^{-1}$ is very extended over the entire southeastern area, and thus probably due to cloud contamination. The redshifted peak, on the other hand, represents concentrated emission from around the continuum source. We fit a Gaussian to the peak and measure a velocity of 26.38 km s$^{-1}$.

## ORCID iDs

Philipp Weber https://orcid.org/0000-0002-3354-6654

Silvio Ulloa https://orcid.org/0000-0002-9966-9403

Sebastián Pérez https://orcid.org/0000-0003-2953-755X

James Miley https://orcid.org/0000-0002-1575-680X

Lucas Cieza https://orcid.org/0000-0002-2828-1153

Sergei Nayakshin https://orcid.org/0000-0002-6166-2206

Alice Zurlo https://orcid.org/0000-0002-5903-8316

Hauyu Baobab Liu https://orcid.org/0000-0003-2300-2626

Fernando Cruz-Sáenz de Miera https://orcid.org/0000-0002-4283-2185

Antonio Hales https://orcid.org/0000-0001-5073-2849

Antonio Garufi https://orcid.org/0000-0002-4266-0643

Dimitris Stamatellos https://orcid.org/0000-0002-4502-8344

Ágnes Kóspál https://orcid.org/0000-0001-7157-6275

Viviana Guzmán https://orcid.org/0000-0003-4784-3040

## References

Aikawa, Y., Herbst, E., Roberts, H., & Caselli, P. 2005, ApJ, 620, 330

Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705

Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. of Arizona Press), 387

Bachiller, R., & Pérez Gutiérrez, M. 1997, ApJL, 487, L93

Bacon, R., Piqueras, L., Conseil, S., Richard, J., & Shepherd, M. 2016, MPDAF: MUSE Python Data Analysis Framework, Astrophysics Source Code Library, ascl:1611.003

Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987

Borchert, E. M. A., Price, D. J., Pinte, C., & Cuello, N. 2022, MNRAS, 517, 4436

Cacciapuoti, L., Macias, E., Gupta, A., et al. 2024, A&A, 682, A61

Caratti o Garatti, A., Garcia Lopez, R., Ray, T. P., et al. 2015, ApJL, 806, L4

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344

Cha, S.-H., & Nayakshin, S. 2011, MNRAS, 415, 3319

Chung, C.-Y., Andrews, S. M., Gurwell, M. A., et al. 2024, ApJS, 273, 29

Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Natur, 535, 258

Codella, C., Cabrit, S., Gueth, F., et al. 2007, A&A, 462, L53
Coutens, A., Liu, H. B., Jiménez-Serra, I., et al. 2019, A&A, 631, A58

Cruz-Sáenz de Miera, F., Coutens, A., Kóspál, Á., et al. 2025, A&A, 696, A18

Cuello, N., Louvet, F., Mentiplay, D., et al. 2020, MNRAS, 491, 504

Delussu, L., Birnstiel, T., Miotello, A., et al. 2024, A&A, 688, A81

Dullemond, C. P., Küffmeier, M., Goicovic, F., et al. 2019, A&A, 628, A20

Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 355

Forgan, D., & Rice, K. 2010, MNRAS, 402, 1349

Fu, Z., Deng, H., Lin, D. N. C., & Mayer, L. 2025, SciA, 11, eadu6058

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

Galván-Madrid, R., Liu, H. B., Izquierdo, A. F., et al. 2018, ApJ, 868, 39

Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104

Gieser, C., Pineda, J. E., Segura-Cox, D. M., et al. 2024, A&A, 692, A55

Ginski, C., Tazaki, R., Dominik, C., & Stolker, T. 2023, ApJ, 953, 92

Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJL, 908, L25

Hartmann, L., & Kenyon, S. J. 1985, ApJ, 299, 462

Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

Herbig, G. H. 1977, ApJ, 217, 693

Hillenbrand, L. 2014, ATel, 6797, 1

Hunter, J. D. 2007, CSE, 9, 90

Jeong, J.-H., Lee, J.-E., Lee, S., et al. 2025, ApJS, 276, 49

Kóspál, Á., Ábrahám, P., Moór, A., et al. 2015, ApJL, 801, L5

Kóspál, Á., Cruz-Sáenz de Miera, F., White, J. A., et al. 2021, ApJS, 256, 30

Kuffmeier, M., Calcutt, H., & Kristensen, L. E. 2019, A&A, 628, A112

Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642

Lee, C.-F. 2020, A&ARv, 28, 1

Lee, J.-E., Lee, S., Baek, G., et al. 2019, NatAs, 3, 314

Li, J. I.-H., Liu, H. B., Hasegawa, Y., & Hirano, N. 2017, ApJ, 840, 72

Liu, H. B. 2019, ApJL, 877, L22

Liu, H. B., Galván-Madrid, R., Jiménez-Serra, I., et al. 2015, ApJ, 804, 37

Liu, H. B., Galván-Madrid, R., Forbrich, J., et al. 2014, ApJ, 780, 155

Liu, H. B., Mérand, A., Green, J. D., et al. 2019, ApJ, 884, 97

Liu, H. B., Muto, T., Konishi, M., et al. 2024, A&A, 685, A18

Liu, H. B., Tsai, A.-L., Chen, W. P., et al. 2021, ApJ, 923, 270

Longarini, C., Armitage, P. J., Lodato, G., Price, D. J., & Ceppi, S. 2023, MNRAS, 522, 6217

Maehara, H., Kojima, T., & Fujii, M. 2014, ATel, 6770, 1

Ménard, F., Cuello, N., Ginski, C., et al. 2020, A&A, 639, L1

Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909

Nayakshin, S., Cruz Sáenz de Miera, F., Kóspál, Á., et al. 2024, MNRAS, 530, 1749

Nayakshin, S., & Elbakyan, V. 2024, MNRAS, 528, 2182

Nayakshin, S., Owen, J. E., & Elbakyan, V. 2023, MNRAS, 523, 385

Nomura, H., Furuya, K., Cordiner, M. A., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 1075

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 233

Podio, L., Codella, C., Gueth, F., et al. 2016, A&A, 593, L4

Rodgers, S. D., & Millar, T. J. 1996, MNRAS, 280, 1046

Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32

Takami, M., Fu, G., Liu, H. B., et al. 2018, ApJ, 864, 20

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Tobin, J. J., van't Hoff, M. L. R., Leemker, M., et al. 2023, Natur, 615, 227

Tychoniec, L., van Dishoeck, E. F., van't Hoff, M. L. R., et al. 2021, A&A, 655, A65

van der Walt, S., Chris Colbert, S., Varoquaux, G., et al. 2011, CSE, 13, 22

van Holstein, R. G., Girard, J. H., & de Boer, J. 2020, A&A, 633, A64

van Holstein, R. G., Stolker, T., Jensen-Clem, R., et al. 2021, A&A, 647, A21

Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956

Vorobyov, E. I., & Basu, S. 2010, ApJL, 714, L133

Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115

Vorobyov, E. I., & Elbakyan, V. G. 2018, A&A, 618, A7

Vorobyov, E. I., Elbakyan, V. G., Liu, H. B., & Takami, M. 2021, A&A, 647, A44

Vorobyov, E. I., Steinrueck, M. E., Elbakyan, V., & Guedel, M. 2017, A&A, 608, A107

Wang, L. Y., Shang, H., & Chiang, T. Y. 2019, ApJ, 874, 31

Weber, P., Pérez, S., Zurlo, A., et al. 2023, ApJL, 952, L17

Winter, A. J., Benisty, M., & Andrews, S. M. 2024, ApJL, 972, L9

Zamponi, J., Maureira, M. J., Zhao, B., et al. 2021, MNRAS, 508, 2583

Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620

Zurlo, A., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 465, 834

Zurlo, A., Weber, P., Pérez, S., et al. 2024, A&A, 686, A309