# Disk Evolution Study through Imaging of Nearby Young Stars (DESTINYS): A Panchromatic View of DO Tau's Complex Kilo-astronomical-unit Environment

Jane Huang, Christian Ginski, Myriam Benisty, Bin Ren, Alexander J. Bohn, Élodie Choquet, Karin I. Öberg, Álvaro Ribas, Jaehan Bae, Edwin A. Bergin, et al.

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

## ► To cite this version:

Jane Huang, Christian Ginski, Myriam Benisty, Bin Ren, Alexander J. Bohn, et al., Disk Evolution Study through Imaging of Nearby Young Stars (DESTINYS): A Panchromatic View of DO Tau's Complex Kilo-astronomical-unit Environment. The Astrophysical Journal, 2022, 930, 10.3847/1538-4357/ac63ba. insu-03690860

HAL Id: insu-03690860

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

Submitted on 20 Jun 2022

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

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

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

# Disk Evolution Study through Imaging of Nearby Young Stars (DESTINYS): A Panchromatic View of DO Tau's Complex Kilo-astronomical-unit Environment

Jane Huang<sup>1,23</sup>*, Christian Ginski<sup>2,3</sup>*, Myriam Benisty<sup>4,5</sup>*, Bin Ren (任彬)<sup>6</sup>*, Alexander J. Bohn<sup>3</sup>*, Élodie Choquet<sup>7</sup>*, Karin I. Öberg<sup>8</sup>*, Álvaro Ribas<sup>9</sup>*, Jaehan Bae<sup>10</sup>*, Edwin A. Bergin<sup>1</sup>*, Til Birnstiel<sup>11,12</sup>*, Yann Boehler<sup>5</sup>*, Stefano Facchini<sup>13,14</sup>*, Daniel Harsono<sup>15</sup>*, Michiel Hogerheijde<sup>2,3</sup>*, Feng Long<sup>8</sup>*, Carlo F. Manara<sup>14</sup>*, François Ménard<sup>5</sup>*, Paola Pinilla<sup>16,17</sup>*, Christophe Pinte<sup>5,18</sup>*, Christian Rab<sup>11,19</sup>*, Jonathan P. Williams<sup>20</sup>*, and Alice Zurlo<sup>7,21,22</sup>†

## Abstract

While protoplanetary disks are often treated as isolated systems in planet formation models, observations increasingly suggest that vigorous interactions between Class II disks and their environments are not rare. DO Tau is a T Tauri star that has previously been hypothesized to have undergone a close encounter with the HV Tau system. As part of the DESTINYS ESO Large Programme, we present new Very Large Telescope (VLT)/SPHERE polarimetric observations of DO Tau and combine them with archival Hubble Space Telescope (HST) scattered-light images and Atacama Large Millimeter/submillimeter Array (ALMA) observations of CO isotopologues and CS to map a network of complex structures. The SPHERE and ALMA observations show that the circumstellar disk is connected to arms extending out to several hundred astronomical units. HST and ALMA also reveal stream-like structures northeast of DO Tau, some of which are at least several thousand astronomical units long. These streams appear not to be gravitationally bound to DO Tau, and comparisons with previous Herschel far-IR observations suggest that the streams are part of a bridge-like structure connecting DO Tau and HV Tau. We also detect a fainter redshifted counterpart to a previously known blueshifted CO outflow. While some of DO Tau's complex structures could be attributed to a recent disk-disk encounter, they might be explained alternatively by interactions with remnant material from the star formation process. These panchromatic observations of DO Tau highlight the need to contextualize the evolution of Class II disks by examining processes occurring over a wide range of size scales.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Radio interferometry (1346); Polarimetry (1278); Planet formation (1241)*

*Supporting material: figure set*

## 1. Introduction

Planet formation has often been modeled as taking place in isolated protoplanetary disks (e.g., Lambrechts &

Johansen 2012; Mordasini et al. 2012; Zhang et al. 2018). However, observations are increasingly challenging the general applicability of such an assumption. Gaps and rings have been detected in the millimeter continuum emission of Class I disks, potentially indicative of planet-disk interactions in these young embedded systems (e.g., Sheehan & Eisner 2018; Alves et al. 2020; Segura-Cox et al. 2020). Meanwhile, molecular and scattered-light observations have revealed large-scale streams, spirals, and tails associated with the ostensibly more evolved Class II protoplanetary disks, suggestive of ongoing infall or

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

<footnote n="1">¹Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA; jnhuang@umich.edu</footnote>

<footnote n="2">²Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands</footnote>

<footnote n="3">³Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386, Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="6">⁶Department of Astronomy, California Institute of Technology, MC 249-17, 1200 East California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷Aix Marseille University, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="8">⁸Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="9">⁹European Southern Observatory, Alonso de Córdova 3017, Vitacura, Santiago, Chile</footnote>

<footnote n="10">¹⁰Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="11">¹¹University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="12">¹²Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany</footnote>

<footnote n="13">¹³Università degli Studi di Milano, via Giovanni Celoria 16, I-20133 Milano, Italy</footnote>

<footnote n="14">¹⁴European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="15">¹⁵Institute of Astronomy, Department of Physics, National Tsing Hua University, Hsinchu, Taiwan</footnote>

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

<footnote n="17">¹⁷Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK</footnote>

<footnote n="18">¹⁸School of Physics and Astronomy, Monash University, Clayton, Vic 3800, Australia</footnote>

<footnote n="19">¹⁹Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, D-85748 Garching, Germany</footnote>

<footnote n="20">²⁰Institute for Astronomy, University of Hawaii, Honolulu, HI 96822, USA</footnote>

<footnote n="21">²¹Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="22">²²Escuela de Ingeniería Industrial, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile
Received 2021 December 14; revised 2022 March 31; accepted 2022 April 1; published 2022 May 17</footnote>

<footnote n="23">²³NASA Hubble Fellowship Program Sagan Fellow.</footnote>
close encounters with other stars (e.g., Nakajima & Goli-
mowski 1995; Cabrit et al. 2006; Tang et al. 2012; Akiyama
et al. 2019; Huang et al. 2021). Together, these observations
indicate that many disks are perturbed by their surroundings
while planet formation is ongoing. Key theorized consequences
of these external perturbations include modification of the total
mass budget available for planet formation, disk truncation,
inducement of misalignments or instabilities, formation of disk
substructures, and stellar outbursts (e.g., Clarke & Pringle 1993;
Bae et al. 2015; Harsono et al. 2011; Manara et al. 2018; Cuello
et al. 2019; Dullemond et al. 2019; Kuffmeier et al.
2020, 2021; Kuznetsova et al. 2022). Thus, mapping disk
environments is crucial for understanding what processes
control protoplanets' formation location and timescale, mass,
orbital behavior, and survival.

High-resolution panchromatic observations are essential for
characterizing interactions between disks and their environ-
ments. Millimeter continuum emission primarily traces the
distribution of millimeter-sized dust grains and is therefore a
good probe of the bulk of the disk mass reservoir. However,
because grain growth is most effective in dense regions and
larger dust grains are more susceptible to inward radial drift
and vertical settling, the millimeter continuum does not trace
extended, more tenuous structures well. Optical and infrared
scattered light trace the distribution of submicron-sized dust
grains. In the context of disk observations, they have primarily
been used to characterize the surface layers of disks (e.g., Krist
et al. 2000; Fukagawa et al. 2004; Avenhaus et al. 2018).
However, scattered-light imaging has also proven effective for
identifying large-scale structures that do not have detected
millimeter continuum counterparts (e.g., Grady et al. 1999;
Garufi et al. 2020). Furthermore, because the strength of the
scattered-light signal and the degree of linear polarization
depend in part on the scattering angle, they can be used to infer
how structures are oriented relative to the illuminating star
(e.g., Stolker et al. 2016; Ginski et al. 2021). Although small
dust grains and gas are thought to be well coupled, molecular
line observations can still reveal material not visible in
scattered light because the detectability of structures with
scattered light is affected by disk self-shadowing and the rapid
decrease in stellar illumination with distance (e.g., Huang et al.
2020; Garufi et al. 2022a). Molecular line observations are also
critical for probing kinematics and gas mass.

The Disk Evolution Study Through Imaging of Nearby
Young Stars (DESTINYS) Large Programme has been
conducting a near-infrared polarimetric survey of young stars
with the SPHERE instrument on the Very Large Telescope
(Ginski et al. 2020, 2021). Whereas early scattered-light studies
disproportionately targeted Herbig Ae stars (see, e.g., Garufi
et al. 2018, and references therein), DESTINYS aims to
provide a comprehensive understanding of the disk structures
of lower-mass pre-main-sequence stars. While the original
primary aims of DESTINYS were to characterize dust
evolution in disks, identify small-scale substructures that may
be signs of planet-disk interactions, and search for stellar and
planetary-mass companions, DESTINYS has serendipitously
provided insights into interactions between disks and their
environments, as demonstrated by the detection of SU Aur's
spectacular infalling tail structures and misaligned disk in
Ginski et al. (2021).

Another DESTINYS target exhibiting signs of interactions
with its surroundings is DO Tau, an M0.3 T Tauri star located

138.4 ± 0.7 pc away in the Taurus star-forming region
(Herczeg & Hillenbrand 2014; Gaia Collaboration et al.
2021; Bailer-Jones et al. 2021). It hosts a Class II proto-
planetary disk associated with a jet and molecular outflow (e.g.,
Hirth et al. 1994; Kenyon & Hartmann 1995; Fernández-López
et al. 2020). Coronagraphic Hubble Space Telescope (HST)
optical and Subaru infrared images revealed an arc-like
structure extending up to a few hundred astronomical units
north of the star (Grady 2004; Itoh et al. 2008). In projection,
DO Tau is situated ~91" (12,600 au) northwest of the T Tauri
triple star system HV Tau. HV Tau A and B have a spectral
type of M1 and have a projected separation from one another of
10 au (Herbig & Bell 1988; Simon et al. 1996). HV Tau C is an
M0.5 star that hosts an edge-on protoplanetary disk and has a
projected separation of ~550 au from HV Tau AB (e.g., Woitas
& Leinert 1998; Duchêne et al. 2010). DO Tau and HV Tau are
surrounded by nebulosities visible at both optical and far-
infrared wavelengths (e.g., Struve & Straka 1962; McGroarty
& Ray 2004; Howard et al. 2013). Based on the morphology of
these nebulosities, Winter et al. (2018) hypothesized that DO
Tau and HV Tau originally constituted a quadruple system and
underwent a close encounter that stripped material from their
disks and ejected DO Tau.

In this work, we combine new SPHERE observations of DO
Tau with the Atacama Large Millimeter/submillimeter Array
(ALMA) molecular line observations and archival optical and
infrared images in order to investigate the spatial and kinematic
relationship between DO Tau's disk and its complex
surroundings. The observations and data reduction are
summarized in Section 2. The scattered-light structures are
described in Section 3, and the molecular emission structures
are described and compared to the scattered-light structures in
Section 4. The properties and potential origins of the structures
are discussed in Section 5, and the results are summarized in
Section 6.

## 2. Observations and Data Reduction

In this section, we describe the various new and archival
observations analyzed in this work. The observations are also
listed in Table 1.

### 2.1. SPHERE/IRDIS Observations

DO Tau was observed with the Very Large Telescope's
SPHERE/IRDIS instrument (Beuzit et al. 2019) using the
broadband *H* filter (central wavelength: 1.6255 µm) on 2019
November 27 (program I.D. 0104.C-0850(A), P.I.: Y. Boehler)
and 2019 December 20 (corresponding to the DESTINYS
Large Programme with I.D. 1104.C-0415(E), P.I.: C. Ginski).
Both programs used dual-beam polarimetric imaging (DPI)
with pupil tracking (de Boer et al. 2020; van Holstein et al.
2020a). In each epoch, the main observing sequence consisted
of 56 frames with the DO Tau star behind an apodized Lyot
coronagraph, which has an inner working angle of 92.5 mas
(Carbillet et al. 2011; Guerri et al. 2011). The integration time
of each frame was 64 s, for a total time of 59.7 minutes per
epoch. Each coronagraphic sequence was preceded and
followed by five flux calibration frames, in which DO Tau
was observed offset from the coronagraph and with the ND1.0
neutral density filter to prevent saturation. Each flux calibration
frame lasted for 0.84 s, for a total integration time of 8.4 s per
epoch. Weather conditions were excellent on both nights. On
<table><caption>**Table 1**
DO Tau Observations Analyzed in This Work</caption><thead><tr><th>Facility</th><th>Program</th><th>Observation Type<sup>a</sup></th><th>References</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center; font-weight: bold;">Main text</td></tr><tr><td>VLT/SPHERE</td><td>0104.C-0850(A)</td><td>H-band polarimetric and total intensity imaging</td><td>1</td></tr><tr><td>VLT/SPHERE</td><td>1104.C-0415(E) (DESTINYS)</td><td>H-band polarimetric and total intensity imaging</td><td>1</td></tr><tr><td>ALMA</td><td>2016.1.00627.S</td><td><sup>12</sup>CO, <sup>13</sup>CO, C<sup>18</sup>O, and CS line imaging</td><td>1, 2, 3, 4</td></tr><tr><td>HST/STIS</td><td>HST-GO-9136</td><td>Broadband optical imaging</td><td>5, 6</td></tr><tr><td>Apache Point Observatory</td><td>Sloan Digital Sky Survey</td><td>g-, r-, i-band imaging</td><td>7, 8</td></tr><tr><td>Herschel PACS</td><td>KPOT_bdent_1 (GASPS)</td><td>160 µm imaging</td><td>9, 10</td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Appendix only</td></tr><tr><td>HST/NICMOS</td><td>HST-GO-7418</td><td>F110W and F160W imaging</td><td>11, 12</td></tr><tr><td>Subaru/CIAO</td><td>o05146</td><td>H-band imaging</td><td>13</td></tr></tbody></table>

**Note.**

<sup>a</sup> This column describes only the observations that appear in this work, not necessarily all observations obtained by the program.

**References.** (1) This work, (2) Bergner et al. (2019), (3) Le Gal et al. (2019), (4) Pegues et al. (2020), (5) Grady (2004), (6) Ren et al. (2017), (7) York et al. (2000), (8) Finkbeiner et al. (2004), (9) Dent et al. (2013), (10) Howard et al. (2013), (11) Choquet et al. (2014), (12) Hagan et al. (2018), (13) Itoh et al. (2008).

2019 November 27, the typical seeing was 0.′′6 and the typical coherence time was 8 ms. On 2019 December 20, the typical seeing was 0.′′5 and the typical coherence time was 10 ms.

The data from each night were individually reduced with the publicly available IRDIS Data reduction for Accurate Polarimetry (IRDAP) pipeline (van Holstein et al. 2017, 2020a, 2020b), yielding polarized intensity images through polarimetric differential imaging (PDI) and total intensity images through classical angular differential imaging (cADI). After checking that results from the two nights were consistent, the reduced images were then averaged together to increase the signal-to-noise ratio (S/N) of the faint extended features.

Because ADI results in self-subtraction artifacts, we also produced total intensity images of DO Tau using reference differential imaging (RDI). DESTINYS Large Programme observations of IP Tau (taken on 2019 December 15) were used to create a point-spread function (PSF) reference library. The observational setup and weather conditions for IP Tau were similar to those for DO Tau. As an M0.6 T Tauri star located in the Taurus star-forming region (Herczeg & Hillenbrand 2014), IP Tau provides a close color match to DO Tau. However, because IP Tau also has a disk, using it as a PSF reference is only useful for recovering structures beyond an arcsecond or so from the star. PSF subtraction was performed with the Karhunen–Loève Image Projection (KLIP) algorithm described in Soummer et al. (2012). Because the standard KLIP procedure results in the overestimation of the stellar signal when extended structures are present, we applied the iterative disk feedback reference differential imaging (IDF-RDI) approach presented in D. Vaendel et al. (in preparation) and Ginski et al. (2021). In brief, the PSF-subtracted image produced with KLIP is treated as an initial estimate of the nonstellar signal. This image is subtracted from the original stack of science images of DO Tau (i.e., the non-PSF-subtracted images) in order to remove as much nonstellar signal as possible. KLIP is then used on the modified science image stack to obtain an improved estimate of the stellar signal, and PSF subtraction is performed on this modified stack. If the resulting PSF-subtracted image still shows significant nonstellar structures, additional iterations are performed until the stellar signal estimate no longer improves. For DO Tau, 10 iterations of IDF-RDI were performed for each epoch, and the results for each epoch were averaged together. Appendix A

shows a comparison of the total intensity images of DO Tau produced with the standard RDI with KLIP (equivalent to the zeroth iteration of IDF-RDI) and after 10 iterations of IDF-RDI.

## 2.2. ALMA Observations

DO Tau was observed with ALMA in two Band 6 spectral settings (henceforth referred to as the 1.3 mm and 1.1 mm settings) as part of program 2016.1.00627.S (PI: K. Öberg). Portions of DO Tau observations from this program have been published previously: C<sup>18</sup>O *J* = 2–1 and upper limits for H<sup>13</sup>CN *J* = 3–2 and C<sub>2</sub>H *N* = 3–2 in Bergner et al. (2019), H<sub>2</sub>CO *J*<sub>K<sub>a</sub>K<sub>c</sub></sub> = 3<sub>03</sub> – 2<sub>02</sub> in Pegues et al. (2020), CS *J* = 5–4 in Le Gal et al. (2019), and upper limits on DCN *J* = 3–2 and HC<sup>15</sup>N *J* = 3–2 in Bergner et al. (2020). In this article, we present the <sup>12</sup>CO and <sup>13</sup>CO *J* = 2–1 observations from this program for the first time. We also reprocess the C<sup>18</sup>O and CS observations for consistency; whereas the aforementioned works only analyzed the Keplerian disk, the present article examines the extended structures in further detail. (The extended structures are not detected in the other previously published line observations.)

Details of the observations, including the dates, array and correlator setups, and calibrators, are provided in Appendix B. Because of the narrow spectral windows used in both spectral settings, bandwidth switching was employed, which required manual calibration by ALMA/NAASC staff. The associated calibration scripts are available on the ALMA archive. Following data delivery from ALMA, we performed additional processing and imaging steps in CASA 5.6.1 (McMullin et al. 2007). For each individual execution block, a set of pseudo-continuum visibilities was produced by flagging channels with strong line emission and then averaging together the remaining channels. A preliminary continuum image was produced for each execution block using the Clark CLEAN algorithm as implemented in the tclean task. After measuring the location of the continuum peak by fitting a two-dimensional Gaussian to the disk emission, each set of visibilities was shifted to align the continuum peak with the phase center. We then compared the fluxes of different execution blocks taken with the same spectral setting. For the 1.1 mm setting, the continuum flux of the third execution block was ~10% lower than the earlier two
<table><caption>**Table 2**
ALMA Imaging Summary</caption><thead><tr><th>Transition</th><th>Rest Frequency<br>(GHz)</th><th>Synthesized Beam<br>(arcsec × arcsec (°))</th><th>Peak I<sub>ν</sub><br>(mJy beam<sup>-1</sup>)</th><th>rms Noise<sup>a</sup><br>(mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2–1</td><td>230.5380000</td><td>0.84 × 0.53 (33°7)</td><td>576.9</td><td>4.7</td></tr><tr><td><sup>13</sup>CO J = 2–1</td><td>220.3986842</td><td>0.87 × 0.56 (32°9)</td><td>201.6</td><td>4.7</td></tr><tr><td>C<sup>18</sup>O J = 2–1</td><td>219.5603541</td><td>0.88 × 0.56 (33°3)</td><td>61.9</td><td>3.9</td></tr><tr><td>CS J = 5–4</td><td>244.9355565</td><td>0.69 × 0.55 (21°3)</td><td>56.2</td><td>3.3</td></tr></tbody></table>

**Note.**

execution blocks. For the 1.3 mm setting, the continuum flux of the second execution block was ~15% lower than the first execution block. These discrepancies are comparable to ALMA's nominal absolute flux calibration accuracy in Band 6 (e.g., Andreani et al. 2016). For the 1.3 mm setting, we rescaled the visibilities of the first execution block such that the continuum flux matched that of the second execution block, and for the 1.1 mm setting, we rescaled the visibilities of the third execution block such that the continuum flux matched that of the first two. The choice of which execution block served as the reference for rescaling is arbitrary, but the interpretation of the data in this work is not affected by this choice.

For each spectral setting, the pseudo-continuum visibilities of execution blocks were imaged together with Clark CLEAN to form an initial model for self-calibration. One round of phase self-calibration and then one round of amplitude self-calibration were applied to the pseudo-continuum visibilities. The self-calibration tables derived from continuum imaging were applied to the full-resolution spectral windows containing line data. Continuum subtraction was then performed in the uv plane using the uvcontsub task. Image cubes were produced for the ¹²CO J = 2–1, ¹³CO J = 2–1, C¹⁸O J = 2–1, and CS J = 5–4 transitions using multiscale CLEAN (scales of [0, 0′′8, 1′′6, 3′′2]) and a Briggs robust parameter of 1.0. Because the line emission is extended and irregular, we used CASA's auto-multithresh automasking algorithm (Kepley et al. 2020) to draw CLEAN masks. After some initial experimentation, we set the automasking parameter values to sidelobethreshold=2.0, noisethreshold=3.5, lownoisethreshold=1.5, minbeamfrac=0.4, and negativethreshold=15.0. The channel spacing was set to 0.25 km s⁻¹ to balance between achieving adequate sensitivity in individual channels and resolving kinematic detail. A primary beam correction was applied to all CLEANed image cubes using the impbcor task. The rms of each image cube was measured within line-free channels of the non-primary-beam-corrected images. The synthesized beam size, rms, and peak intensity of each image cube are given in Table 2.

## 2.3. Archival Optical/Infrared Observations

### 2.3.1. HST STIS

DO Tau was imaged with the HST STIS instrument as part of program HST-GO-9136 (P.I.: C. Grady) on 2001 December 8 and 2003 February 12, with a total exposure time of 2268 s per epoch. The STIS CCD has a broad bandpass that extends from 2000 to 10,300 Å, with a central wavelength of 5740 Å. For each set of observations, DO Tau was placed under the coronagraphic mask at the WEDGE1.0 location, which has a

width of 1". These STIS data were first presented in Grady (2004).

The calibrated, flat-fielded science files for DO Tau were retrieved from the Mikulski Archive for Space Telescopes. Following the procedure described in Ren et al. (2017), a median filter was applied to correct bad pixels. The two rolls were rotated such that north was up and east was to the left, and then aligned using the centerRadon package (Ren et al. 2019). The two rolls were then combined, with pixels averaged in regions where the two rolls overlapped. However, if a pixel from one roll fell within a diffraction spike, the combined image only used the pixel value from the other roll. The portion of the combined image lying within 8" of DO Tau's stellar position was then replaced with the KLIP PSF-subtracted image presented in Ren et al. (2017). PSF subtraction over the entire field of view of the DO Tau observations was not feasible due to limitations in the field of view of the PSF reference library. The PSF wings are negligible outside a radius of 8" but the intensities are discontinuous at the boundary between the PSF-subtracted portion of the combined DO Tau image and the non-PSF-subtracted portion. Thus, caution should be applied in interpreting features at this boundary.

### 2.3.2. SDSS

The field around DO Tau was imaged by the Sloan Digital Sky Survey (SDSS) on 2002 December 31 (MJD 52639). The original SDSS run-camcol-field identifier for the data is 3559-5-55. Details of the SDSS observations and data reduction are provided in York et al. (2000), Finkbeiner et al. (2004), and Ahn et al. (2012). We downloaded pipeline-processed SDSS g-, r-, and i-band cutouts from NASA Skyview (McGlynn et al. 1998) and produced a color composite using the Astropy (Astropy Collaboration et al. 2013) implementation of the Lupton et al. (2004) algorithm.

### 2.3.3. Herschel PACS

DO Tau was imaged at 160 µm with the Herschel PACS instrument (Poglitsch et al. 2010) on 2011 March 30 as part of the GASPS Open Key Time Project (Dent et al. 2013; Howard et al. 2013). We retrieved the pipeline-processed Level 2.5 Unimap map (which combined data from observation IDs 1342217478, 1342217479, 1342217480, and 1342217481) from the Herschel Science Archive.

## 3. Morphology in Scattered Light

DO Tau is surrounded by a series of complex structures spanning several orders of magnitude in size scales, as shown in the SDSS, HST STIS, and SPHERE images presented in Figure 1. The SDSS optical color image shows the broad

<footnote>a With channel widths of 0.25 km s⁻¹.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) SDSS color composite image of DO Tau and HV Tau made from g-, r-, and i-band images. The blue box shows the relative size of the STIS image shown in panel (b). The white stars mark the locations of HV Tau A and B, HV Tau C, and DO Tau (only a single star is used to mark HV Tau A and B because of their close proximity). (b) HST STIS image of DO Tau, smoothed by convolving with a 2D Gaussian with a standard deviation of 0.1 (equal to two pixels). The gray curves mark the locations of the large northeastern arc and the α, β, and γ streams. The yellow star shows the location of DO Tau. Pixels falling within the coronagraphic wedges or stellar diffraction spikes are masked. The dotted white circle denotes the region where PSF subtraction was performed with KLIP. The pink box shows the relative size of the SPHERE image in part (c). The axes are marked with the angular offsets from DO Tau. (c) SPHERE H-band total intensity image of DO Tau, produced with cADI and smoothed by convolving with a 2D Gaussian with a standard deviation of 0.0368 (equal to three pixels). The solid curves mark the northern and eastern arms, the dotted white curve marks the outflow shell, and the dotted line marks the jet. The gray circle denotes the extent of the coronagraph. Versions of the STIS and SPHERE images with color bars and no annotations are presented in Appendix C.</figcaption></figure>

nebulosity surrounding DO Tau and HV Tau, which are separated by ~12,600 au. This nebula is relatively well known, having previously been listed in catalogs such as those by

Struve & Straka (1962) and Magakian (2003). Within this broad nebulosity, an arc-like structure (which we call the “large northeastern arc”) is visible to the northeast of DO Tau in the
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Left: smoothed SPHERE $Q_{\phi}$ H-band image of DO Tau. The gray circle marks the extent of the coronagraph. The solid white curves mark the eastern and northern arms. The dotted white curve marks part of an outflow shell. Offsets from the star are labeled on the axes. North is up, and east is to the left. A version of this image without annotations is presented in Figure 17. Right: an inset of the SPHERE $Q_{\phi}$ H-band image of DO Tau. Unlike the left-hand panel, no smoothing is applied to the right-hand panel in order to show more detail at smaller separations from the star. Note that the color stretch is also different from that of the left-hand panel in order to show brightness variations across the disk region more clearly.</figcaption></figure>

SDSS image but can be seen in finer detail in the HST STIS optical image. The southeastern end of the arc is located at a position angle (P.A.) east of north of ~110° and a projected separation of ~14" (~1900 au) from DO Tau, while the northwestern end is located at a P.A. of ~310° and a separation of ~27" (~3700 au) from the star. There are also at least three stream-like structures that are oriented nearly perpendicularly to the large northeastern arc. To our knowledge, these stream-like structures have not been explicitly identified in the literature; previous publications of this HST data (Grady 2004; Itoh et al. 2008; Ren et al. 2017) only presented cropped versions of the image. We label the stream-like structures as α, β, and γ (in order from north to south). The southwestern ends of these structures appear to point toward a smaller arc north of DO Tau. This smaller arc was first detected by Grady (2004) in the same set of STIS observations.

Although the STIS wedges obscure much of the material within several hundred astronomical units of DO Tau, SPHERE provides a clearer view of these structures. We identify a counterpart to the HST image's smaller arc, which we call the “northern arm.” We also more clearly identify an “eastern arm” that is only partially visible in the HST STIS image. Although the eastern and northern arms appear as though they could be a single structure in the total intensity images, we identify them as separate structures because they are kinematically distinct and because they appear to connect separately to the disk in polarized scattered light, as discussed later in this section and in Section 4. In addition, a narrow structure that is west of DO Tau and partially obscured by the wedges in the STIS image turns out to be part of a larger arc-like structure southwest of DO Tau. Based on visual inspection of the SPHERE cADI image in Figure 1, it can be approximated as an elliptical arc centered ~1"7 (~240 au) west and ~0"4 (~55 au) south of the star, with a P.A. of 80°, an angular extent of ~300°, a semimajor axis of ~1"65 (~230 au), and a semiminor axis of ~1"4 (~190 au). This geometry matches one of the blueshifted CO outflow shells identified in

Fernández-López et al. (2020). In the cADI image, the outflow shell is bisected by a linear feature that has the same orientation as the blueshifted side of the [Fe II] jet imaged by Erkal et al. (2021). Erkal et al. (2021) detected emission from the 1.53 and 1.64 µm [Fe II] lines, which fall within the wavelength range (but are not spectrally resolved by) of SPHERE's broadband H filter. Outside of the arm structures, a background point source is visible ~3"5 east of DO Tau and is discussed further in Appendix D.

While the SPHERE cADI image is dominated by self-subtraction artifacts within ~1" of the star, the $Q_{\phi}$ image clarifies the relationship between some of the arm structures and the circumstellar disk (Figure 2). The quantity $Q_{\phi}$ is defined as

$$Q_{\phi} = -Q_{\text{IPS}} \cos(2\phi) - U_{\text{IPS}} \sin(2\phi), \quad (1)$$

where $\phi$ is the azimuthal angle in the observer frame and $Q_{\text{IPS}}$ and $U_{\text{IPS}}$ are the linear polarization components of the Stokes vector following correction for instrumental polarization (e.g., de Boer et al. 2020). $Q_{\phi}$ is related to the linearly polarized intensity $PI$ by $PI = \sqrt{Q_{\phi}^2 + U_{\phi}^2}$, where $U_{\phi}$ is defined as

$$U_{\phi} = Q_{\text{IPS}} \sin(2\phi) - U_{\text{IPS}} \cos(2\phi). \quad (2)$$

A positive signal in the $Q_{\phi}$ image traces azimuthal polarization, which is expected to be dominant for relatively low-inclination disks such as DO Tau (27°.6 as measured from millimeter continuum emission in Long et al. 2019).

In projection, the northern arm appears to emerge from the disk at a P.A. of ~114° east of north. The arm curves up to a separation of ~3"3 (~460 au) northwest of DO Tau and then turns over toward the southwest, terminating at a P.A. of ~256° and a separation of ~3"7 (~510 au) from the star. (These quantities are estimated via visual inspection because there is no sharp delineation between the disk and arm structures). The eastern arm extends from the north side of the
disk and terminates at a P.A. of ~117° and a separation of ~2.2"3 (~320 au) from the star.

The circumstellar disk itself appears relatively compact in scattered light, with a radius of ~80 au as determined by visual inspection. However, the size of the disk is not straightforward to quantify due to the aforementioned extended structures. Furthermore, scattered light only provides a lower bound on the disk size because sensitivity decreases with the square of the distance from the star and many disks are also self-shadowed (e.g., Garufi et al. 2022a). DO Tau's disk appears about twice as large in scattered light compared to millimeter continuum. The former predominantly traces submicron-sized dust grains in the disk atmosphere, while the latter primarily traces millimeter-sized dust grains in the disk midplane. Long et al. (2019) find that 95% of the millimeter continuum flux is enclosed within a radius of only 36.4 au. The millimeter continuum image from Long et al. (2019) exhibits no substructures at an angular resolution of ~0.′′1 (14 au). Likewise, no gaps or rings are visible in the scattered-light image. The disk signal has a broader extent on the west side compared to the east side, which we interpret as a projection effect of viewing a flared, inclined disk surface. Some contamination from the outflow shell may also enhance the difference between the west and east sides. Because the east side appears to be foreshortened, we infer that it is inclined toward the observer. Based on the obscuration of DO Tau's jet by the eastern side of the disk, Erkal et al. (2021) likewise concluded that the east side is tilted toward the observer.

The degree of linear polarization (DoLP), defined as the ratio of the linearly polarized intensity to total intensity, provides constraints on the three-dimensional orientation of the extended structures around DO Tau. Because the extended structures around DO Tau are not detected in the $U_φ$ images, we assume that $PI ≈ Q_φ$, so that we can estimate the DoLP without the additional noise contribution from the $U_φ$ component. To calculate the DoLP, we use the total intensity image produced with IDF-RDI rather than cADI because the latter procedure generally results in significant self-subtraction of the extended structures. However, because IDF-RDI does not yield a perfect subtraction of the stellar PSF, the resulting DoLP values must still be treated with caution. Therefore, our analysis focuses on relative rather than absolute DoLP values. To increase the S/N of the DoLP map, we smoothed the $Q_φ$ and total intensity images using a Gaussian kernel with a standard deviation of one pixel (0.013") prior to division and then binned the resulting DoLP map by a factor of 2.

Panels (a) and (b) of Figure 3 show the IDF-RDI total intensity image and DoLP map of DO Tau. Within the DoLP map, the northern arm stands out most clearly. Part (c) of Figure 3 shows how the DoLP of the northern arm varies from east to west. As the arm emerges from the eastern side of the disk, the DoLP values are low (<0.2). The DoLP values generally increase from east to west across the arm, peaking at a value of ~0.5, but then decrease again near the western terminus of the arm. If we assume that the DoLP has a roughly bell-shaped dependence upon the scattering angle, with a peak at ~90° (e.g., Min et al. 2005), then the DoLP pattern of the northern arm could be explained by one of several scenarios: (1) The scattering angles increase from east to west across the arm, starting at a value <90° and ending at a value >90°. (2) The scattering angles decrease from east to west across the arm, starting at a value >90° and ending at a value <90°. (3)

Starting from a value <90°, the scattering angles increase from east to west up to 90°, then decrease again at the western end of the arm. (4) Starting from a value >90°, the scattering angles decrease from east to west down to 90°, then increase again at the western end of the arm. Given that the northern arm appears to emerge from the eastern side of the disk, which is tilted toward the observer, we expect the eastern side of the northern arm to have small scattering angles. Furthermore, the eastern side of the arm is brighter in total intensity compared to the western side, suggesting stronger forward scattering on the eastern side. This brightness difference was previously noted by Itoh et al. (2008) in Subaru/CIAO H-band observations of DO Tau and also ascribed to forward scattering. Scenarios 2 and 4 therefore appear unlikely. We favor scenario 1 over scenario 3 because we do not observe the total intensity of the arm increasing at the western terminus, as one might expect if scattering angles were decreasing from 90° and forward scattering became stronger. However, scenario 3 is not definitively ruled out because the distance from the star also affects total intensity. The eastern arm appears to have total intensity and DoLP values comparable to the eastern side of the northern arm, suggesting that they have similar (small) scattering angles. Meanwhile, the DoLP in the vicinity of the outflow shell appears to be close to zero. Because the disk has a relatively low inclination and CO kinematics show that the outflow shell is in front of the disk (Fernández-López et al. 2020), the shell is expected to have low scattering angles and therefore a low DoLP. A schematic of the proposed geometry is presented in panel (d) of Figure 3.

# 4. Molecular Line Observations

## 4.1. Overview

Integrated intensity maps and intensity-weighted velocity maps of $^{12}$CO J = 2–1, $^{13}$CO J = 2–1, C$^{18}$O J = 2–1, and CS J = 5–4 are shown in Figure 4. The velocity integration ranges are selected based on the extent over which emission above the 3$σ$ level is detected: –13.75 to 20.5 km s$^{-1}$ for $^{12}$CO, –0.75 to 9.25 km s$^{-1}$ for $^{13}$CO, 4.0 to 8.75 km s$^{-1}$ for C$^{18}$O, and 4.0 to 8.25 km s$^{-1}$ for CS. Because the spatial distribution of the emission changes substantially from channel to channel, most of the channels at any given spatial location are emission free even though the overall velocity range of the emission is large. To reduce the contribution of signal-free regions to the integrated intensity maps, we only include pixels in the image cube with intensities exceeding the 2$σ$ level. For the intensity-weighted velocity maps, a 5$σ$ pixel clip is used instead because higher moments are more sensitive to outliers. The full set of channel maps is presented in Appendix E.

The characteristic rotation velocity pattern for a Keplerian disk is visible in the velocity-weighted intensity maps of $^{13}$CO, C$^{18}$O, and CS, but $^{12}$CO is dominated by nondisk emission. The extended emission structures traced by all four of these molecules are described in more detail below.

## 4.2. Outflow Shells

The $^{12}$CO emission is dominated by a series of blueshifted, overlapping ring-like structures southwest of DO Tau. These structures were previously detected in an independent set of ALMA $^{12}$CO J = 2–1 observations by Fernández-López et al. (2020), who attributed them to a series of outflow shells. We do not detect these shells in the other three molecular tracers. A
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** (a) SPHERE total intensity image of DO Tau produced with IDF-RDI. The image is binned by a factor of 2 relative to the native IRDIS pixel resolution in order to increase the S/N of the extended structures, which are marked by white curves. The central region is masked because the presence of a disk around the PSF reference star prevents accurate estimates of the total intensity at small separations from DO Tau. (b) Degree of linear polarization map for DO Tau. Values less than 0 or greater than 1 are masked because they are unphysical. (c) DoLP values measured across the northern arm. The blue curve denotes the mean values measured within 6 pixel × 6 pixel boxes (roughly 0"/15 × 0"/15) centered on points along the white curve corresponding to the northern arm in panel (b). The shaded region shows the standard deviation measured within each of these boxes. (d) A schematic of an overhead view of the different structures traced by SPHERE around DO Tau, showing how various structures are inclined relative to the line of sight. North points out of the page, and south points into the page. The structures are not drawn to scale. The redshifted jet and outflow are drawn with dashed curves because they are not detected with SPHERE but are known to be present based on other observations (see Erkal et al. 2021 and Section 4.)</figcaption></figure>

comparison of the SPHERE H-band cADI image to the ¹²CO integrated intensity map shows the coincidence between the southwestern arc detected in scattered light and the innermost CO outflow shell (Figure 5).

Thanks to the longer integration time of our observations compared to those presented in Fernández-López et al. (2020) (~4× longer), we are able to detect a fainter, redshifted outflow. Figure 6 compares integrated intensity maps of the redshifted and blueshifted outflows. The blueshifted outflow is detected at local standard of rest (Kinematic), or LSRK, velocities ranging from roughly −13.75 to 4.0 km s⁻¹, while the redshifted outflow is detected from approximately 9.5 to

20.5 km s⁻¹. (Fernández-López et al. 2020 identify ¹²CO outflow shells at LSRK velocities between 4.4 and 8.2 km s⁻¹, but we attribute this emission instead largely to the northern and eastern arms based on comparison with the ¹³CO emission in Section 4.3.) The peak integrated intensity of the blueshifted outflow is 2.1 Jy beam⁻¹ km s⁻¹, while that of the redshifted outflow is only 0.21 Jy beam⁻¹ km s⁻¹. Likewise, Erkal et al. (2021) found that the blueshifted side of the [Fe II] jet is brighter than the redshifted side. Whereas the blueshifted outflow manifests as multiple rings of different sizes, the redshifted outflow appears to be a single large arc west of DO Tau (although the widening of the emission at the
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Left column: integrated intensity maps of ¹²CO J = 2–1, ¹³CO J = 2–1, C¹⁸O J = 2–1, and CS J = 5–4 toward DO Tau (in order from top to bottom). The synthesized beam is shown as a solid white ellipse in the lower-left corner of each panel. The white star marks the position of the millimeter continuum peak, which is located at the phase center. Offsets from the phase center (in arcseconds) are marked in the lower-left panel. North is up, and east is to the left. The dotted white ellipse marks the FWHM of the ALMA primary beam. An arcinsh stretch is applied to the color scale to highlight faint emission features. Right column: corresponding intensity-weighted velocity maps. Note that the color scale for ¹²CO is different from the other molecules to accommodate the wider velocity range over which emission is detected.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** A comparison between the SPHERE H-band cADI image of DO Tau (left) and the ¹²CO integrated intensity map (right) on the same relative scale. The dotted white curve marks the position of a scattered-light arc that coincides with the innermost outflow shell as traced by ¹²CO.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Overlaid integrated intensity maps of the high-velocity ¹²CO emission toward DO Tau. No clipping is applied to make these integrated intensity maps because the redshifted outflow emission levels are low in individual channels. The velocity integration ranges are given at the top of the figure. The blue contours, corresponding to the blueshifted emission, are drawn at 75, 150, 400, 700, 1000, 1300, and 1600 mJy beam⁻¹ km s⁻¹. The red contours, corresponding to the redshifted emission, are drawn at 50, 100, and 150 mJy beam⁻¹ km s⁻¹. The black star marks the location of DO Tau. The dotted black circle marks the FWHM of the primary beam. Offsets from the phase center (in arcseconds) are marked in the lower-left panel. North is up, and east is to the left.</figcaption></figure>

northern end suggests that it may trace two slightly offset rings). Because of the faintness of the redshifted outflow in individual channels (see Appendix E), it is ambiguous how the morphology varies as a function of velocity.

As shown in Figure 6, the two ends of the redshifted arc extend toward the FWHM of the primary beam, and the emission blends into the noise as sensitivity decreases away from the phase center. Given the presence of ring-like blueshifted outflow emission, we consider it likely that the

redshifted arc is part of a ring-like outflow shell that extends beyond ALMA’s field of view. To constrain the geometry of the redshifted outflow, we first smoothed the image with a Gaussian kernel with a standard deviation of 0.3 to increase the S/N of the faint arc and then measured the positions $x_i, y_i$ of local intensity maxima along horizontal (east-west) slices of the image spaced 0.8 apart in the vertical (north-south) direction. The northern and southern bounds of where measurements were taken were determined based on whether a distinct peak was identifiable in a given image slice. We then modeled the arc as an ellipse with the free parameters $x_0, y_0, a_0, f$, and $\theta$, where $x_0$ and $y_0$ are offsets (in arcseconds) from DO Tau, $a_0$ is the semimajor axis (also in arcseconds), $f$ is the ratio of the semiminor axis to semimajor axis, and $\theta$ is the angle (in radians) between the semimajor axis and the northern direction. The parameters $x_0$ and $y_0$ are defined such that the positive direction is toward the east and north, respectively. The log-likelihood function is written as

$$ \log \mathcal{L} = -\frac{1}{2} \sum_i \left[ \frac{d_i^2}{\sigma^2} + \log(2\pi\sigma^2) \right], \quad (3) $$

where $d_i$ is the smallest distance between point i and the model ellipse, and $\sigma$ is the standard deviation of the major axis of the synthesized beam. We calculate $d_i$ using the analytical expressions described in Zhang (1997). Flat priors were specified for all free parameters, with bounds of [0", 20"] for $x_0$, [-20", 20"] for $y_0$, [5", 20"] for $a_0$, [0, 1] for $f$, and [0, $\pi$] for $\theta$. We used emcee to explore the posteriors, with 40 walkers evolved over 2000 steps each. After discarding the first 1000 steps as burn-in, we find the following posterior medians, with the uncertainties corresponding to the 16th and 84th percentiles: $x_0 = 7.3_{-1.9}^{+2.6}$ arcseconds, $y_0 = -0.8_{-0.8}^{+1.1}$ arcseconds, $a_0 = 14.6_{-1.5}^{+2.6}$ arcseconds, $f = 0.90_{-0.08}^{+0.07}$, and $\theta = 1.2 \pm 0.9$ rad. The ellipse with the median parameter values is shown in Figure 7. The center of the ellipse is tentatively offset to the south from the axis of DO Tau’s [Fe II] jet; the former has a P. A. of 96°±7°, while the latter has a P.A. of 260° ± 3° (Erkal et al. 2021). However, Erkal et al. (2021) note that the redshifted side of the [Fe II] jet may be misaligned from the blueshifted side. The redshifted outflow also appears to be
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Model ellipse with median parameter values (in red) overlaid on an integrated intensity map of the redshifted outflow as traced by $^{12}$CO (including velocities between 9.5 and 20.5 km s$^{-1}$). The red cross marks the center of the model ellipse and associated $1\sigma$ error bars. The dotted red points denote the positions of the local intensity maxima measured along horizontal cuts across the image. The yellow star marks the position of DO Tau. The white dotted circle denotes the FWHM of the primary beam. The dashed red and blue line segments show the P.A. of DO Tau's jet as measured by Erkal et al. (2021). The shaded regions show the $1\sigma$ uncertainties in the jet orientation. The line segments are drawn longer than the observed extent of the jet in order to show the apparent offset from the center of the redshifted outflow.</figcaption></figure>

slightly misaligned relative to the axis of the blueshifted molecular outflow, which has a P.A. of 253°.6 ± 0°.1 (Fernández-López et al. 2020). Deeper, mosaicked CO observations will be required to measure the full extent of the redshifted molecular outflow in order to confirm whether it is indeed misaligned. Erkal et al. (2021) inferred that DO Tau's jet was precessing, so precession may account for misalignment of the molecular outflow.

The detected redshifted and blueshifted outflow shells have strikingly different radii. Whereas the shell radii in the blueshifted outflow range from ~200–1000 au (Fernández-López et al. 2020), the redshifted outflow shell has a much larger radius of 2000 au. We estimate the dynamical age of the redshifted outflow shell as $t_{dyn} = D_{outflow}/v_{outflow}$, where $D_{outflow}$ is the deprojected distance from DO Tau to the edge of the outflow shell and $v_{outflow}$ is the deprojected velocity of the outflow relative to DO Tau ($v_{outflow} = (v_{LSRK} - v_{sys})/\cos i$). DO Tau's systemic velocity $v_{sys}$ is 5.9 km s$^{-1}$ (Pegues et al. 2020). To compute $D_{outflow}$, we make the simplifying assumption that the outflow axis is perpendicular to the outer disk, for which the inclination has been constrained. Because the outflow misalignment does not seem to be extreme, this simplification should not significantly affect the dynamical timescale estimate. We take the disk inclination to be 27.6, based on the millimeter continuum measurement in Long et al. (2019). Emission from the redshifted outflow is detected at LSRK velocities from ~9.5 to 20.5 km s$^{-1}$ and does not have a well-defined peak along the spectral axis. We therefore use the upper and lower ends of the velocity range to

estimate that $t_{dyn}$ is between ~900 and 3500 yr, which is comparable to or older than the range of dynamical ages that Fernández-López et al. (2020) estimated for the blueshifted outflow shells (~460–1090 yr).

### 4.3. The Northern and Eastern Arms

The northern and eastern arms detected in scattered light also have counterparts in $^{12}$CO, $^{13}$CO and CS emission (Figure 8). A $^{18}$O counterpart is not visible in the moment maps, but a portion of it is tentatively visible at an LSRK velocity of 5.5 km s$^{-1}$ in the channel maps (Appendix E). In $^{12}$CO emission, the northern arm is only partially visible due to spatial filtering in channels near the systemic velocity and to contamination from other structures. The northern arm is more clearly visible in $^{13}$CO emission and appears to extend farther south compared to its scattered-light counterpart. It is also possible, though, that this apparent southward extension is contamination from another structure. In CS, only the eastern side of the northern arm is detected. Peculiarly, the portion of the arm brightest in CS emission coincides with the region where the arm becomes more tenuous in $^{13}$CO emission. Observing additional transitions of CS and $^{13}$CO to better constrain the physical conditions across the arm may yield useful insight into why CS exhibits this emission pattern.

The kinematics of the northern and eastern arms are best probed by $^{13}$CO, which is less contaminated by outflows compared to $^{12}$CO but has a higher S/N compared to CS. Figure 9 compares the SPHERE $Q_\phi$ image to the $^{13}$CO channel maps. Emission from the northern arm is detected roughly between 5.00 and 7.75 km s$^{-1}$. The eastern end of this arm emerges from the side of the Keplerian disk that is blueshifted relative to the systemic velocity (5.9 km s$^{-1}$). The behavior of the western side of the northern arm is more difficult to interpret. Blueshifted emission overlaps the western side of the arm in the channels from 5.50 to 5.75 km s$^{-1}$. Little to no emission is detected on the western side of the arm in the channels from 6.00 to 6.25 km s$^{-1}$. Then, a bright redshifted component overlaps the northern arm from 6.50 to 7.75 km s$^{-1}$. The two distinct kinematic components may indicate that what appears to be a single continuous arm in the SPHERE image may actually be multiple structures overlapping in projection, or perhaps that the arm is expanding. Meanwhile, the eastern arm is detected between 6.00 and 7.75 km s$^{-1}$, emerging from the redshifted side of the Keplerian disk.

In Section 3, we inferred from the total intensity and DoLP map that the northern arm is oriented such that the eastern side is tilted toward the observer. If this is indeed the case, then its blueshifted motion indicates that it is moving away from the disk. Meanwhile, the redshifted motion of the eastern arm indicates that it is moving opposite to the eastern side of the northern arm.

### 4.4. Stream-like Structures

We detect several stream-like structures (labeled $\alpha_{CO}$, $\beta_{CO}$, $\delta_{CO}$, and $\epsilon_{CO}$) northeast of DO Tau in $^{12}$CO emission at velocities between 6.75 and 8.75 km s$^{-1}$ (Figure 10), which is redshifted relative to DO Tau's systemic velocity of 5.9 km s$^{-1}$. We also tentatively identify a $\gamma_{CO}$ feature. Its geometry is less certain due to low sensitivity outside the primary-beam FWHM as well as confusion with the $\delta_{CO}$ and $\epsilon_{CO}$ structures. In $^{13}$CO emission, we identify counterparts to $\beta_{CO}$, $\delta_{CO}$, and $\epsilon_{CO}$. Much
<figure><img src="image_13.png" /><figcaption>**Figure 8.** A comparison of the SPHERE $Q_φ$ image of DO Tau to insets of the $^{12}$CO, $^{13}$CO, and CS integrated intensity maps. The integration range of $^{12}$CO is truncated to 4.00–8.25 km s$^{-1}$ in order to minimize contamination from the outflows. The white curves mark the locations of the northern and eastern arms. The white star denotes the position of DO Tau. Offsets from the phase center (in arcseconds) are marked on the leftmost panel. North is up, and east is to the left.</figcaption></figure>

of this emission appears at a velocity range (5.75 km s$^{-1}$ to 6.5 km s$^{-1}$), where $^{12}$CO exhibits severe spatial filtering. Thus, in the intensity-weighted velocity maps (Figure 4), the stream-like structures appear blueshifted in $^{13}$CO compared to $^{12}$CO. C$^{18}$O also exhibits extended emission that overlaps with the $\delta_{CO}$ feature.

In projection, $\alpha_{CO}$, $\beta_{CO}$, $\delta_{CO}$, and $\epsilon_{CO}$ all appear to converge on the western side of DO Tau's northern arm. The $\alpha_{CO}$ and $\beta_{CO}$ streams can be traced up to the edge of the image, with a projected separation of ~19" (~2600 au) from DO Tau. The appearance of the streams is adversely affected by foreground contamination, loss of sensitivity away from the phase center, and spatial filtering (the maximum recoverable scale of the ALMA configuration used is ~4"). However, the correspondence between the $\alpha_{CO}$ and $\beta_{CO}$ and the $\alpha$ and $\beta$ streams detected in the HST STIS image lends greater confidence to the characterization of their geometry (Figure 11). The $\gamma_{CO}$ stream appears to be slightly offset to the southeast from the $\gamma$ stream identified with HST, but observations with better uv coverage and sensitivity should be obtained to confirm whether this offset is genuine.

A comparison of the stream locations with the Herschel PACS 160 µm image (Figure 11) suggests that the $\alpha_{CO}$ and $\beta_{CO}$ (and perhaps $\gamma_{CO}$) streams form the base of an elongated structure that extends northeast of DO Tau. This elongated structure has previously been described as a tidal tail that forms part of a bridge structure connected to HV Tau (Winter et al. 2018). Because the Herschel PACS image has a relatively low spatial resolution (~11" at 160 µm), we are unable to identify individual infrared counterparts to the ALMA and HST streams.

To determine whether the stream-like structures are gravitationally bound, we compute the radius $r$ at which the line-of-sight velocity relative to DO Tau is equal to the escape velocity, $\sqrt{2GM_*/r}$. Structures outside this radius must be moving faster than the escape velocity and are thus gravitationally unbound. (Structures inside this radius may or may not be unbound depending on what the other velocity components are.) In making this calculation, we assume that the mass of the gas and dust surrounding DO Tau is small compared to the stellar mass of $0.54 \pm 0.07 M_\odot$ (Braun et al. 2021). The resulting escape

velocity radii are plotted in Figure 12. Because the streams generally lie outside these radii, we conclude that they are not gravitationally bound.

### 4.5. Mass Constraints on Structures Surrounding DO Tau

Under the assumption that $^{13}$CO is optically thin outside DO Tau's Keplerian disk, we can use it to estimate a lower bound for the gas mass of the structures around the disk. However, it should be kept in mind that the flux may be significantly underestimated due to spatial filtering. Furthermore, the ALMA field of view does not capture the full extent of these structures.

In the LTE limit, column density can be estimated from optically thin emission with

$$N = \frac{Q(T)}{g_u} \exp\left(\frac{E_u}{T}\right) \frac{4\pi}{A_{ul}hc} \int I_{\nu} dv, \quad (4)$$

where $Q$ is the partition function, $g_u$ is the upper-state degeneracy, $E_u$ is the upper-state energy, $T$ is the excitation temperature, $A_{ul}$ is the Einstein A coefficient, and $\int I_{\nu} dv$ is the integral of the intensity over the velocity axis (e.g., Goldsmith & Langer 1999). For $^{13}$CO $J=2-1$, $g_u=10$ (note that the Cologne Database for Molecular Spectroscopy includes hyperfine splitting in its partition function and degeneracy), $E_u=15.9$ K, and $A_{ul}=6.08 \times 10^{-7} \text{ s}^{-1}$ (Müller et al. 2001, 2005).

To estimate the flux, we made a $^{13}$CO integrated intensity map without clipping in order to avoid biasing measurements upward. We then extracted the flux within the highlighted region shown in Figure 13, finding a value of ~4.5 Jy km s$^{-1}$. This value excludes emission from the disk, which we estimate has a radial extent of ~250 au based on visual inspection. Although the boundary between the disk and the extended structures is not well defined, the exact choice of value does not affect the order of magnitude of the mass estimate for the extended structures because they dominate the flux due to their large emitting area. The mean integrated intensity within the region of extraction is ~11 mJy beam$^{-1}$ km s$^{-1}$. To derive a mean column density, we adopt a gas temperature of 30 K, based on temperatures that
<figure><img src="image_14.png" /><figcaption>**Figure 9.** A comparison of the northern and eastern arms as traced by polarized scattered light (black and white image) and ¹³CO emission. The ¹³CO contours are drawn at the 4σ, 6σ, and 8σ levels. North is up, and east is to the left. The synthesized beam is drawn as a white ellipse in the lower-left corner of each panel, while the LSRK velocity (km s⁻¹) is in the top-right corner. The blue curves mark the northern and eastern arms. The red curves mark the Keplerian disk.</figcaption></figure>

Winter et al. (2018) estimated from Herschel dust maps of the region around DO Tau. At this temperature, $Q = 23.4$ for ¹³CO (Müller et al. 2001, 2005). We find a mean ¹³CO

column density of ~3 × 10¹⁴ cm⁻². We assume that the composition of the gas is similar to that of local molecular clouds, i.e., the gas is predominantly molecular hydrogen,
<figure><img src="image_15.png" /><figcaption>**Figure 10.** A comparison of integrated intensity maps of stream-like structures detected in ¹²CO (left), ¹³CO (middle), and C¹⁸O (right). The velocity range of the integrated intensity maps is labeled above each panel. The color scale is saturated to make faint features more visible. Due to severe spatial filtering between 5.25 and 6.5 km s⁻¹, ¹²CO is integrated over a different velocity range from the other two isotopologues. The locations of the αCO, βCO, γCO, δCO, and ϵCO streams are marked with white curves. The gray curve on the ¹³CO image shows the position of DO Tau's northern arm. The gray star marks the position of DO Tau in each panel. The white dotted circle marks the FWHM of the primary beam.</figcaption></figure>

¹²CO:H₂ = 10⁻⁴, and ¹²CO:¹³CO = 69 (e.g., Frerking et al. 1982; Wilson 1999). We estimate a lower bound of ~4 × 10⁻⁴ M☉ for the gas mass of the extended structures. Because ¹³CO is likely optically thick within the Keplerian disk, we do not use it to estimate the disk mass. However, Kwon et al. (2015) estimated a total disk mass of 0.014 M☉ based on radiative transfer modeling of the millimeter continuum. Thus, the mass of the extended structures appears to be at least a few percent of the disk mass.

Based on Herschel far-infrared maps of dust emission, Winter et al. (2018) derived a lower bound of 10⁻² M☉ for the total mass of extended structures surrounding DO Tau and HV Tau. This is much higher than our lower bound of ~4 × 10⁻⁴ M☉, which is due at least in part to the Herschel measurements being taken over a much larger area. DO Tau and HV Tau have a projected separation of ~91", whereas the FWHM of the primary beam of the ALMA ¹³CO image is only ~29". Spatial filtering of ¹³CO is almost certainly responsible for part of the discrepancy as well. Other possibilities are that the gas-to-dust ratio or ¹²CO:H₂ ratio differs from interstellar medium (ISM) values. This could be the case if the material originated from within the protoplanetary disks because the gas-to-dust and ¹²CO:H₂ ratios in disks are both expected to decrease over time (e.g., Gorti et al. 2015; Miotello et al. 2019). Sturm et al. (2022) find evidence from [C I] observations that carbon is depleted by a factor of 17 from ISM levels in the outer regions of the DO Tau disk, which would in turn reduce the ¹²CO:H₂ ratio. Wide-field observations with better uv coverage will be necessary to robustly examine the consistency between the ¹³CO-derived mass and the Herschel-derived mass.

## 4.6. Summary of Structures Associated with DO Tau

Figure 14 presents a schematic showing the spatial relationship between the various structures around DO Tau identified in scattered-light and molecular line emission. The smallest structure characterized is the circumstellar disk itself, with a radial extent of ~80 au in scattered light. The northern and eastern arms, both connected to the disk, trace structures up to several hundred astronomical units from DO Tau. The streams and the large northeastern arc are much larger in scale, with extents of several thousand astronomical units. The scales of the molecular outflow shells range from hundreds to thousands of astronomical units.

## 5. Discussion

Here we discuss potential origins for the arms, stream-like structures, and large northeastern arc near DO Tau. We limit the discussion of the jet and molecular outflow because they have already recently been analyzed at length by Fernández-López et al. (2020) and Erkal et al. (2021).

### 5.1. Evidence for an Encounter between DO Tau and HV Tau

Winter et al. (2018) presented stellar encounter simulations reproducing the general morphology of the dust “bridge” connecting DO Tau and HV Tau, as imaged by Herschel. They inferred that DO Tau and HV Tau began as a quadruple system, but then subsequently became either unbound or only loosely bound. However, the simulation resolution was not high enough for comparison with the complex network of structures detected by HST, SPHERE, and ALMA. Arm-like structures reminiscent of DO Tau’s northern and eastern arms, though, do appear in high-resolution stellar flyby simulations of other
<figure><img src="image_16.png" /><figcaption>**Figure 11.** (a) ALMA $^{12}$CO integrated intensity map of DO Tau, annotated in white with the locations of the $\alpha_{CO}$, $\beta_{CO}$, and $\gamma_{CO}$ streams. Only velocities from 6.75 to 8.75 km s$^{-1}$ are included. The gray star marks the position of DO Tau. Axes are marked with the angular offsets with respect to the star. (b) HST STIS image of DO Tau, with the dashed gray curves denoting the locations of the $\alpha_{CO}$, $\beta_{CO}$, and $\gamma_{CO}$ streams. The solid gray curves mark the $\alpha$, $\beta$, and $\gamma$ features identified within the HST image from Section 3. The yellow star marks the location of DO Tau. The bright end of the color scale is saturated in order to show the contrast between the large-scale structures more clearly. (c) Herschel 160 µm image showing the nebulosity surrounding HV Tau and DO Tau. The blue curves correspond to the white annotations drawn on the ALMA image denoting the $\alpha_{CO}$, $\beta_{CO}$, and $\gamma_{CO}$ streams. The pink contours correspond to intensity levels of [5, 10, 15, 20, 25] mJy pixel$^{-1}$. The red stars mark the locations of HV Tau and DO Tau. (Only one star is used to represent HV Tau A and B due to their close proximity. HV Tau C is to the northeast of HV Tau A and B, as shown in Figure 1).</figcaption></figure>

systems (e.g., Clarke & Pringle 1993; Pfalzner et al. 2008; Dai et al. 2015; Cuello et al. 2019).

The DO Tau system also bears notable similarities to the RW Aur system, which has been hypothesized to have undergone a recent stellar encounter (e.g., Dai et al. 2015; Rodriguez et al. 2018). Like DO Tau, RW Aur features multiple arm-like structures in CO emission (Cabrit et al. 2006; Rodriguez et al. 2018). Both systems have also been noted for their unusually strong photometric variability on timescales of several weeks or longer (Rodriguez et al. 2016; Rigon et al. 2017). Among other explanations, a misaligned disk has been hypothesized to be responsible for RW Aur's dimming events (Facchini et al. 2016; Koutoulaki et al. 2019). A similar link may exist for DO Tau, which has been hypothesized to have a misaligned inner disk based on evidence of a precessing jet (Erkal et al. 2021). Hydrodynamical simulations have demonstrated that stellar encounters can lead to disk warps and misalignments, albeit with rapid realignment timescales (e.g., Clarke & Pringle 1993; Cuello et al. 2019; Nealon et al. 2020). In general, disk

misalignments have been more commonly attributed to the presence of a misaligned embedded companion, which simulations suggest can yield more extreme and longer-lasting disk warps compared to stellar encounters (e.g., Facchini et al. 2018; Zhu 2019; Nealon et al. 2020). Nevertheless, the evidence for misalignments in both the RW Aur A and DO Tau disks motivates further examination of the long-term impact of external perturbations on the evolution of the inner disk. Characterizing the mechanisms that can lead to disk misalignment is key for determining the extent to which processes that occur during the gas-rich protoplanetary disk phase may be responsible for the misaligned orbits observed in mature planetary systems (e.g., Batygin 2012; Lai 2014; Spalding & Batygin 2016).

One characteristic that sets DO Tau apart from most other systems that have been hypothesized to have undergone a stellar encounter (e.g., RW Aur, AS 205, UX Tau, FU Ori) is that DO Tau's current projected separation from HV Tau is ~12,600 au, whereas the current projected separations between
<figure><img src="image_17.png" /><figcaption>**Figure 12.** $^{12}$CO J = 2–1 channel maps annotated in white with the locations of the α$_{CO}$, β$_{CO}$, γ$_{CO}$, δ$_{CO}$, and ε$_{CO}$ streams. The γ$_{CO}$ feature is marked with a dashed rather than a solid line because its identification is less certain. The dotted pink circle denotes the radius at which the line-of-sight velocity relative to DO Tau is equal to the escape velocity. The LSRK velocity (km s$^{-1}$) appears in the top right of each panel. The synthesized beam is shown as a solid white ellipse in the lower-left corner. A gray cross marks the position of the millimeter continuum peak, which is located at the phase center. Offsets from the phase center (in arcseconds) are marked in the lower-left panel. The top end of the color bar is saturated in order to show the faint extended features in more detail. North is up, and east is to the left.</figcaption></figure>

the other hypothesized stellar encounter pairs are typically only on
the order of a few hundred astronomical units (e.g., Forgan &
Rice 2010; Kurtovic et al. 2018; Zapata et al. 2020; Ménard et al.
2020; Borchert et al. 2022). Dong et al. (2022) recently identified
a more extreme case in which Z CMa appears to have undergone
an encounter with a source now separated in projection by ~4700
au, but this is still considerably smaller than the projected
separation between DO Tau and HV Tau. A difficulty with the

stellar encounter scenario for DO Tau and HV Tau is that in
general, simulations suggest that arm-like structures produced by
stellar encounters should dissipate within several thousand years
(e.g., Cuello et al. 2019). However, a kinematic analysis by
Winter et al. (2018) indicates that an encounter between DO Tau
and HV Tau would have had to take place ~0.1 Myr ago to match
their current separation. (The simulated bridge structure between
DO Tau and HV Tau in Winter et al. 2018 starts to dissipate
<figure><img src="image_18.png" /><figcaption>**Figure 13.** $^{13}$CO integrated intensity map with the colored region between the two solid pink ellipses showing where the flux from the extended structures is measured. The star marks the location of DO Tau. The white dotted circle marks the FWHM of the primary beam.</figcaption></figure>

before the stars reach their present-day separation, but this
behavior was due to numerical limitations of the SPH simulation
setup).

Winter et al. (2018) raised the caveat that the proper motions of
DO Tau and HV Tau AB suggest that rather than moving away
from each other (as one would expect if DO Tau had been ejected
after an encounter), the systems may be moving toward one
another. The Gaia EDR3 (Gaia Collaboration et al. 2021) proper
motion of DO Tau is $\mu_{\alpha}^* = 6.27 \pm 0.04$ mas yr$^{-1}$,
$\mu_{\delta} = -21.05 \pm 0.03$ mas yr$^{-1}$ and the proper motion of HV
Tau AB is $\mu_{\alpha}^* = 4.78 \pm 0.06$ mas yr$^{-1}$, $\mu_{\delta} = -21.4 \pm 0.05$ mas
yr$^{-1}$. Proper motion information is not available for HV Tau C.
However, Winter et al. (2018) also point out that because HV Tau
AB is not resolved by Gaia, its derived proper motion does not
necessarily correspond to that of its center of mass and could
therefore give an inaccurate impression of its velocity relative to
DO Tau (see also Clarke 2020).

Stellar ejections have been hypothesized as a way to power
molecular outflows (e.g., Bally et al. 2011). Based on the dynamical
ages of the outflow shells estimated in this work and Fernández-
López et al. (2020), though, the launch of these outflow shells
appears to be a more recent event than the hypothesized encounter
between DO Tau and HV Tau, which Winter et al. (2018)
estimated took place 0.1 Myr ago. Nevertheless, it is worth
considering how the after-effects of the stellar encounter might have
affected the behavior of the outflow, which is noticeably
asymmetric. One possible cause for jet and outflow asymmetries
is environmental heterogeneity (e.g., Hirth et al. 1994; Arce &
Goodman 2002). DO Tau's local environment is quite complex,
which may be due in part to the disruptive influence of a stellar
encounter. However, molecular clouds are also intrinsically
heterogeneous, so it is not necessarily straightforward to isolate
the influence of an encounter, as discussed below.

## 5.2. Evidence for Remnant Envelope Material

DO Tau has traditionally been identified as a Class II
system based on its SED (e.g., Kenyon & Hartmann 1995;

Luhman et al. 2010). Because of the disk’s bright gas and dust
emission, it has frequently been included in studies aimed at
characterizing trends in Class II disk properties (e.g., Najita
et al. 2013; Howard et al. 2013; Williams & Best 2014; Long
et al. 2019; Bergner et al. 2019). Stellar age estimates for DO
Tau have spanned a wide range, from ~0.4 to 6 Myr (e.g.,
Andrews et al. 2018; Long et al. 2019), but the system exhibits
several attributes more characteristic of the younger side of the
estimated age range. DO Tau’s large northeastern arc resembles
the arc-like reflection nebulae that have been observed around
the partially embedded HL Tau and DG Tau systems (e.g.,
Nakajima & Golimowski 1995; Stapelfeldt et al. 1995). In
addition, CO emission maps of HL Tau and DG Tau have
revealed arcs and streams on scales of hundreds to thousands of
astronomical units, reminiscent of those detected in the CO
maps of DO Tau (e.g., Yen et al. 2017; Güdel et al. 2018). DO
Tau’s strong CO outflow is also suggestive of youth, although
molecular outflows have on rare occasions been detected in
other Class II systems (e.g., Frank et al. 2014; Louvet et al.
2018). Another potential indication of relative youth is DO
Tau’s spectral energy distribution (Figure 15), which exhibits
higher mid-IR and far-IR emission than most other Class II
systems in Taurus. (While a value of $A_v = 3$ from Andrews
et al. 2013 was assumed in order to deredden DO Tau’s SED,
the discrepancy with other Taurus Class II sources is still
present if we instead adopt $A_v = 0.75$ from Herczeg &
Hillenbrand 2014.) This excess may be due to the presence of
residual envelope material and/or to DO Tau’s disk being
comparatively less settled (e.g., Wilking et al. 1989; Dullemond
& Dominik 2004). Both characteristics are often (though not
exclusively) associated with younger systems.

While the simulations in Winter et al. (2018) could
qualitatively reproduce the large bridge structure connecting
DO Tau and HV Tau solely through the ejection of disk
material following a close encounter, the morphological
similarities between DO Tau and sources known to be partially
embedded raise the question of whether some or perhaps even
most of the material surrounding DO Tau actually originated
from outside the disk. The formation of bridge-like structures
between stars is not uniquely an outcome of tidal stripping by
stellar encounters. Bridge-like structures appear in simulations
of binary formation through turbulent fragmentation, although
separations exceeding $10^4$ au are uncommon past 100 kyr (e.g.,
Offner et al. 2016; Lee et al. 2019). Simulations by Kuffmeier
et al. (2019) demonstrate that bridge-like structures can form
between protostars in multiple systems due to compression of
their natal filament. However, given that bridge structures that
form in this manner only survive for a few tens of thousands of
years in the simulations, this mechanism does not provide a
fully satisfactory explanation for the bridge structure connect-
ing DO Tau and HV Tau.

The possibilities that DO Tau underwent a stellar encounter
and that the system is partially embedded are not mutually
exclusive. Stellar interactions are expected to be more probable
during the embedded phase because protostars often form in
close groups and because the presence of an envelope increases
the cross section of a system (e.g., Gutermuth et al. 2005;
Bate 2018). The large-scale structures around DO Tau could
therefore represent some combination of disk material stripped
by a stellar encounter and disturbed remnant envelope material.

At any rate, the observations of DO Tau provide further
motivation to reexamine the environments of other disks
<figure><img src="image_19.png" /><figcaption>**Figure 14.** A schematic of the structures surrounding DO Tau, based on SPHERE and HST scattered-light observations and ALMA molecular line observations presented in this work, as well as jet observations from Erkal et al. (2021). The yellow star marks the position of DO Tau. All other structures are labeled within the schematic. Note that the CO streams likely extend further eastward than depicted in this schematic, but we are limited by ALMA's and HST's field of view.</figcaption></figure>

traditionally considered to be Class II disks to determine if they likewise are (partially) embedded or interacting with cloud material. Star formation simulations from Padoan et al. (2014) suggest that Class II sources experience infall from nearby cloud material even up to ages of at least several million years. In this case, DO Tau's complex large-scale structures would not necessarily be indicative of youth. A growing number of Class II disks have been identified as potentially experiencing “late infall,” at ages of roughly a Myr or more (e.g., Tang et al. 2012; Dullemond et al. 2019; Ginski et al. 2020; Huang et al. 2021; Mesa et al. 2022). Ongoing accretion of external material during the Class II phase has been suggested as a way to increase the budget of material available to form planets within a given disk (Manara et al. 2018; Dullemond et al. 2019).

### 5.3. Avenues for Clarifying the Origins of DO Tau's Extended Structures

Determining the extent to which the material surrounding DO Tau is primordial is key for understanding how disks evolve. Under the assumption that all the material surrounding DO Tau and HV Tau was stripped from their disks, the simulations in Winter et al. (2018) indicated that the stellar

encounter caused the disks to lose half their mass. However, if some of the external material is primordial, then the impact of stellar encounters on disk mass may not be so extreme. Furthermore, if DO Tau is still interacting with its envelope, it may not be appropriate to group it with envelope-free Class II disks in population studies.

One limitation in our analysis is that some of the structures detected by ALMA extend beyond its field of view, so we do not know their complete morphology. While the low-resolution Herschel dust maps suggest that the streams extend toward HV Tau, mosaicked ALMA observations would be needed to check if they are connected. Meanwhile, the kinematics of the large northeastern arc are unknown because it does not have a clear counterpart in the molecular emission maps, which could be due to spatial filtering. Mosaicked maps with better uv coverage could be used to investigate whether the structure is gravitationally bound to DO Tau.

High-resolution hydrodynamical simulations of stellar encounters with disks both with and without envelopes could help to clarify where DO Tau's extended structures might have originated. Attempting to reproduce the morphology and kinematics of DO Tau's structures in detail would be a very
<figure><img src="image_20.png" /><figcaption>**Figure 15.** A comparison between the spectral energy distribution of DO Tau and other Class II disks in Taurus. The blue points (with $1\sigma$ error bars) denote dereddened photometric measurements presented in Andrews et al. (2013) and Ribas et al. (2017). An $A_v$ value of 3 is adopted based on Andrews et al. (2013). The red curve corresponds to the Spitzer/IRS spectrum retrieved from the Cornell Atlas of Spitzer/IRS Sources (Lebouteiller et al. 2011, 2015). The dashed black curve shows the estimated stellar photosphere. The solid black curve corresponds to the median dereddened Taurus Class II SED from Ribas et al. (2017), while the gray shaded area corresponds to the interquartile range.</figcaption></figure>

computationally formidable task, though, because of the wide range of size scales involved and the large space of initial conditions that would need to be explored. Simulations spanning disk to sub-parsec scales have generally focused on the first couple hundred thousand years of a protostellar system's life (e.g., Kuffmeier et al. 2017; Bate 2018; Kuffmeier et al. 2018; Lebreuilly et al. 2021), but observations of complex large-scale structures associated with more evolved systems are highlighting the need to extend these types of simulations to later stages.

Another potential avenue to distinguish between material ejected from the disk and remnant material from the star formation process is to compare the composition of the gas around DO Tau to that of (partially) embedded systems, isolated Class II disks, and Class II stellar flyby candidate systems. While spatially resolved molecular observations of Class I and II disks have been steadily growing (e.g., Pegues et al. 2020; Bergner et al. 2020; van't Hoff et al. 2020), published spatially resolved molecular line observations of stellar flyby candidates have largely been limited to CO (e.g., Kurtovic et al. 2018; Rodriguez et al. 2018; Zapata et al. 2020). Further molecular line observations of the latter type of object, in conjunction with astrochemical modeling, will be needed to determine whether material stripped from the disk has a chemical signature distinct from envelope material. For example, CS emission in the embedded DG Tau system exhibits arm-like structures reminiscent of those detected around DO Tau (Garufi et al. 2022b). However, it remains to be seen whether other stellar flyby candidate systems have similar CS emission arms. Other species worth targeting in DO Tau and other stellar flyby candidates include SO and $SO_2$, which often appear to be associated with accretion shocks in embedded systems but are rarely detected in envelope-free Class II systems (e.g., Sakai et al. 2014; Artur de la Villarmois et al. 2019; Garufi et al. 2022b). It will also be useful to improve constraints on DO Tau's C/O ratio, which is expected to increase between the Class I and II stages (e.g., Cleeves et al. 2018; Miotello et al. 2019). Results have thus far been mixed;

Sturm et al. (2022) derived a relatively high C/O ratio of ~1 based on spatially unresolved [C I] observations but found that the $C_2H$ upper limit from Bergner et al. (2019) was better explained by an ISM-like ratio of 0.47. Obtaining deeper and higher resolutions of [C I] and $C_2H$, or estimating the C/O ratio from SO/CS, can help to resolve this ambiguity.

While DO Tau is not itself an FU Ori object, McGroarty & Ray (2004) suggested that the presence of nearby Herbig-Haro objects may be evidence of a past FU Ori-like outburst. The complex system of extended structures around DO Tau is also reminiscent of scattered-light and molecular emission images of FU Ori systems (Liu et al. 2016; Zurlo et al. 2017). The mechanism responsible for triggering the strong outbursts in FU Ori-type objects (and potentially sculpting complex large-scale structures around them) is still debated; hypotheses include gravitational instability, accretion from surrounding cloud material, binary interactions, and stellar encounters (e.g., Bonnell & Bastien 1992; Vorobyov & Basu 2005; Forgan & Rice 2010; Audard et al. 2014; Dullemond et al. 2019; Borchert et al. 2022; Dong et al. 2022). Further observational and theoretical investigation of FU Ori-type systems may provide insight into DO Tau's history as well.

## 6. Summary

We analyzed new and archival scattered-light and molecular line observations to map the extent, relative orientation, and kinematics of the series of complex structures surrounding DO Tau. Our findings are as follows:

1. We detect DO Tau's circumstellar disk in polarized scattered light for the first time. Like the millimeter continuum emission presented in Long et al. (2019), the disk does not exhibit any clear substructure in scattered light.

2. The SPHERE scattered-light and ALMA molecular line observations show that the disk is connected to two arm-like structures, one on the eastern side and one on the northern side. They extend up to several hundred astronomical units away from the disk.

3. The HST scattered light and ALMA CO emission trace multiple stream-like structures extending northeast of DO Tau over scales of at least a couple thousand astronomical units. In projection, these streams mostly appear to converge at DO Tau's northern arm. The streams are redshifted with respect to DO Tau's systemic velocity and are not gravitationally bound to DO Tau. Based on their orientation, the streams appear to be part of the bridge-like structure connecting DO Tau and HV Tau as seen in Herschel far-infrared images.

4. We also detect a faint redshifted counterpart to the blueshifted outflow previously detected in CO emission. The redshifted outflow appears to be misaligned with respect to the jet (as traced by [Fe II]), but wider-field CO mapping will be necessary to confirm the geometry of the redshifted outflow.

5. While some of DO Tau's complex structures are compatible with a previously hypothesized stellar encounter, the system also shows signposts of still being partially embedded. Future high-resolution simulations of protostellar evolution up to the Class II stage and observational characterization of the composition of the
gas in and around DO Tau will be key for clarifying the
origins of these structures.

This case study of DO Tau highlights the utility of spatially
resolved panchromatic observations for probing the environ-
ments of protoplanetary disks. For DO Tau, no single tracer
reveals all of the associated large-scale structures. Similar joint
analyses of scattered-light and molecular line observations of
other systems will be key for understanding how the disk
population is affected by interactions with their environments.

This work is based on observations collected at the European
Southern Observatory under ESO program(s) 0104.C-0850(A)
and 1104.C-0415(E), and based on observations made with the
NASA/ESA Hubble Space Telescope, obtained from the data
archive at the Space Telescope Science Institute. STScI is
operated by the Association of Universities for Research in
Astronomy, Inc. under NASA contract NAS 5-26555. This
paper makes use of ALMA data ADS/JAO.ALMA#2016.1.
00627.S.

ALMA is a partnership of ESO (representing its member
states), NSF (USA) and NINS (Japan), together with NRC
(Canada) and NSC and ASIAA (Taiwan), 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 Asso-
ciated Universities, Inc. Funding for SDSS-III has been
provided by the Alfred P. Sloan Foundation, the Participating
Institutions, the National Science Foundation, and the U.S.
Department of Energy Office of Science. The SDSS-III website
is http://www.sdss3.org/. SDSS-III is managed by the
Astrophysical Research Consortium for the Participating
Institutions of the SDSS-III Collaboration including the
University of Arizona, the Brazilian Participation Group,
Brookhaven National Laboratory, Carnegie Mellon University,
University of Florida, the French Participation Group, the
German Participation Group, Harvard University, the Instituto
de Astrofísica de Canarias, the Michigan State/Notre Dame/
JINA Participation Group, Johns Hopkins University, Lawr-
ence Berkeley National Laboratory, Max Planck Institute for
Astrophysics, Max Planck Institute for Extraterrestrial Physics,
New Mexico State University, New York University, Ohio
State University, Pennsylvania State University, University of
Portsmouth, Princeton University, the Spanish Participation
Group, University of Tokyo, University of Utah, Vanderbilt
University, University of Virginia, University of Washington,
and Yale University. This research has made use of data
reprocessed as part of the ALICE program, which was
supported by NASA through grants HST-AR-12652 (PI: R.
Soummer), HST-GO-11136 (PI: D. Golimowski), HST-GO-
13855 (PI: E. Choquet), HST-GO-13331 (PI: L. Pueyo), and
STScI Director's Discretionary Research funds, and was
conducted at STScI which is operated by AURA under NASA
contract NAS5-26555. This work has made use of the SPHERE
Data Centre, jointly operated by OSUG/IPAG (Grenoble),
PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange
(Nice), Observatoire de Paris/LESIA (Paris), and Observatoire
de Lyon/CRAL, and is supported by a grant from Labex

OSUG@2020 (Investissements d'avenirANR10 LABX56). We
acknowledge the use of NASA's SkyView facility (http://
skyview.gsfc.nasa.gov) located at NASA Goddard Space Flight
Center. Herschel is an ESA space observatory with science
instruments provided by European-led Principal Investigator
consortia and with important participation from NASA. This
research is based in part on data collected at the Subaru
Telescope, which is operated by the National Astronomical
Observatory of Japan. We are honored and grateful for the
opportunity to use observations from Maunakea, which has
cultural, historical, and natural significance in Hawaii. This
research has made use of NASA's Astrophysics Data System.
We thank Sean Andrews, Sarah Sadavoy, and Andrew Winter
for useful discussions. We also thank the referee for comments
improving the clarity of the manuscript.

Support for J.H. was provided by NASA through NASA
Hubble Fellowship grant \#HST-HF2-51460.001-A awarded
by the Space Telescope Science Institute, which is operated by
the Association of Universities for Research in Astronomy,
Inc., for NASA, under contract NAS5-26555. M.B. acknowl-
edges funding from the European Research Council (ERC)
under the European Union's Horizon 2020 research and
innovation program (grant PROTOPLANETS No.
101002188). T.B. acknowledges funding from the European
Research Council (ERC) under the European Union's Horizon
2020 research and innovation program under grant agreement
No 714769 and funding by the Deutsche Forschungsge-
meinschaft (DFG, German Research Foundation) under grants
361140270, 325594231, and Germany's Excellence Strategy-
EXC-2094-390783311. D.H. is supported by CICA through a
grant and grant No. 110J035319 from the Ministry of Education
of Taiwan. A.Z. acknowledges support from the FONDECYT
Iniciación en investigación project number 11190837. C.H.R.
is grateful for support from the Max Planck Society. This
project has received funding from the European Union's
Horizon 2020 research and innovation program under the
Marie Sklodowska-Curie grant agreement No 823823
(DUSTBUSTERS).

Facilities: ALMA, VLT:Melipal, HST (NICMOS, STIS).

Software: analysisUtils (https://casaguides.nrao.edu/
index.php/Analysis_Utility), AstroPy (Astropy Collabora-
tion et al. 2013), bezier (Hermes 2017), CASA (McMullin
et al. 2007), centerRadon (Ren et al. 2019), cmasher (van
der Velden 2020), IRDAP (van Holstein et al. 2020b),
matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020),
SkyView (McGlynn et al. 1998), Tiny Tim (Krist et al. 2011).

**Appendix A**

Comparison of Standard RDI and IDF-RDI Images of
DO Tau

Figure 16 compares total intensity images of DO Tau produced through standard RDI with KLIP and through IDF-RDI. While the extended structures around DO Tau appear similar in the two images, IDF-RDI mitigates the over-subtraction artifacts visible in the RDI image.
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Left: total intensity image of DO Tau produced with standard RDI (equivalent to the zeroth iteration of IDF-RDI). Right: total intensity image of DO Tau after 10 iterations of IDF-RDI.</figcaption></figure>

## Appendix B
### Details of ALMA Observational Setup

Details of the ALMA observations of DO Tau are provided in Tables 3, 4, and 5.

<table><caption>**Table 3**
Summary of ALMA Observational Setup</caption><thead><tr><th>Date</th><th>Antennas</th><th>Baselines (m)</th><th>Time on Source (minutes)</th><th>Flux Calibrator</th><th>Bandpass Calibrator</th><th>Phase Calibrator</th></tr></thead><tbody><tr><td colspan="7">1.3 mm Setting</td></tr><tr><td>2016 Dec 1</td><td>44</td><td>15-704</td><td>10</td><td>J0423-0120</td><td>J0237+2848</td><td>J0426+2327</td></tr><tr><td>2016 Dec 1</td><td>44</td><td>15-704</td><td>10</td><td>J0423-0120</td><td>J0510+1800</td><td>J0426+2327</td></tr><tr><td colspan="7">1.1 mm Setting</td></tr><tr><td>2016 Dec 1</td><td>44</td><td>15-704</td><td>9</td><td>J0510+1800</td><td>J0510+1800</td><td>J0426+2327</td></tr><tr><td>2016 Dec 2</td><td>45</td><td>15-704</td><td>9</td><td>J0510+1800</td><td>J0510+1800</td><td>J0426+2327</td></tr><tr><td>2016 Dec 3</td><td>44</td><td>15-704</td><td>9</td><td>J0423-0120</td><td>J0237+2848</td><td>J0426+2327</td></tr></tbody></table>
<table><caption>**Table 4**
1.3 mm Spectral Setting</caption><thead><tr><th>Center Frequency<sup>a</sup><br>(GHz)</th><th>Bandwidth<br>(MHz)</th><th>Native Channel Spacing<br>(MHz)</th></tr></thead><tbody><tr><td>217.227</td><td>58.6</td><td>0.061</td></tr><tr><td>218.210</td><td>58.6</td><td>0.061</td></tr><tr><td>218.313</td><td>58.6</td><td>0.061</td></tr><tr><td>218.721</td><td>58.6</td><td>0.061</td></tr><tr><td>219.548</td><td>117.2</td><td>0.061</td></tr><tr><td>220.387</td><td>117.2</td><td>0.061</td></tr><tr><td>230.525</td><td>117.2</td><td>0.061</td></tr><tr><td>231.309</td><td>117.2</td><td>0.061</td></tr><tr><td>232.406</td><td>937.5</td><td>0.244</td></tr></tbody></table>

**Note.**

<table><caption>**Table 5**
1.1 mm Spectral Setting</caption><thead><tr><th>Center Frequency<sup>a</sup><br>(GHz)</th><th>Bandwidth<br>(MHz)</th><th>Native Channel Spacing<br>(MHz)</th></tr></thead><tbody><tr><td>242.903</td><td>937.5</td><td>0.244</td></tr><tr><td>244.212</td><td>117.2</td><td>0.061</td></tr><tr><td>244.925</td><td>117.2</td><td>0.061</td></tr><tr><td>258.146</td><td>117.2</td><td>0.061</td></tr><tr><td>259.001</td><td>117.2</td><td>0.061</td></tr><tr><td>260.507</td><td>117.2</td><td>0.061</td></tr><tr><td>262.031</td><td>117.2</td><td>0.061</td></tr></tbody></table>

**Note.**

<footnote><sup>a</sup> Topocentric frequency.</footnote>

<footnote><sup>a</sup> Topocentric frequency.</footnote>
## Appendix C
Nonannotated versions of Scattered-light Images

Nonannotated versions of the HST and SPHERE images presented in the main text are shown in Figure 17.

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

<figure><img src="image_23.png" /><figcaption>Figure 17. Nonannotated versions of the HST STIS and SPHERE images presented in Figures 1 and 2.</figcaption></figure>

## Appendix D
Redetection of a Point Source in the Field of DO Tau

Itoh et al. (2008) reported the detection of a point source east of DO Tau in H-band observations taken with the Subaru CIAO instrument on 2005 November 12. We redetect this point source in both epochs of the SPHERE H-band observations as well as in archival NICMOS observations in the F110W and F160W filters.

Together, these observations provide a 21 yr baseline to examine whether the point source is bound to DO Tau.

### D.1. *Archival Data Reduction*

Because Itoh et al. (2008) did not report a position angle for the point source in their Subaru H-band observations from program o05146 (PI: M. Tamura), we retrieved the raw data from Subaru's
<figure><img src="image_24.png" /><figcaption>**Figure 18.** A gallery of images showing the point source detected in the vicinity of DO Tau. Each panel is labeled with the instrument, filter, and date of observation. The white arrow points to the location of the point source. The axes of the leftmost image are labeled with the offset in arcseconds with respect to DO Tau. The black circles in the NICMOS images and gray circles in the CIAO and SPHERE images denote the extent of the coronagraphic mask. An azimuthally averaged radial profile has been subtracted from the CIAO image to make the point source stand out more clearly. The two SPHERE cADI images have been convolved with a Gaussian with a standard deviation of 3 pixels (0.0368).</figcaption></figure>

SMOKA archive for reanalysis. Details of the observations are provided in Itoh et al. (2008). We performed dark subtraction, filtered bad pixels, recentered the central star, and median-combined the individual exposures.

To obtain a longer baseline to analyze the motion of the point source, we checked the HST archive for earlier observations. DO Tau was imaged by HST NICMOS with the NIC2 coronagraph on 1998 December 02 in the F110W (central wavelength: 1.1 µm) and F160W (central wavelength: 1.6 µm) filters as part of program HST-GO-7418 (P. I.: D. Padgett). The exposure time in each filter was 512 s. The observations were processed with the ALICE pipeline as described in Choquet et al. (2014) and Hagan et al. (2018). The sizes of the DO Tau images included in the ALICE public data release<fnref n="24" /> are truncated to 6" × 6"; however, because the point source reported by Itoh et al. (2008) is outside the field of view of the public version of the images, we produced new 10" × 10" versions. In brief, the input science images and PSF reference library images come from the database of NICMOS observations recalibrated by the LAPLACE project (Schneider et al. 2010). PSF subtraction for DO Tau was performed with KLIP. We used a PSF library of 18 frames and subtracted two KL modes for the F110W image, while we used a PSF library of 100 frames and subtracted one KL mode for the F160W image.

The NICMOS, CIAO, and SPHERE images of the point source are shown in Figure 18. The point source is not detected in the HST STIS optical image, for which Itoh et al. (2008) estimated a detection limit of 28 mag.

## D.2. Astrometric Analysis

For the NICMOS observations, the astrometry was performed following the methodology described in Choquet et al. (2015). For each filter, a synthetic NICMOS PSF was generated using the Tiny Tim software (Krist et al. 2011). The synthetic PSF was used as a matched filter template to find the position and flux values that maximize the correlation with the data. The 1σ astrometric uncertainty was taken to be half the size of the NIC2 pixel. The S/N was estimated by measuring the standard deviation in an annulus around the point source.

<figure><img src="image_25.png" /><figcaption>**Figure 19.** Measurements of the P.A. of the point source in the field of DO Tau, shown as a function of time. The gray ribbon shows the expected values for a stationary background source. The dashed black lines show the expected position angle as a function of time if the point source had a face-on circular orbit (i.e., the orientation at which the change in P.A. would be maximized). The blue dashed-dotted line marks the P.A. measured for the most recent observation.</figcaption></figure>

However, the background estimates may be skewed by the extended structures surrounding DO Tau.

In the Subaru observations, the point source noted by Itoh et al. (2008) is visible near the bright diffraction spike to the east of the central star. To extract the position of the point source in detector coordinates, a 2D Gaussian was fitted at the source position, which yields an uncertainty in the P.A. of the source relative to the DO Tau of 0.33. However, this uncertainty does not account for systematic offsets from a dedicated true north astrometric calibration. Following the discussion in Ginski et al. (2014), we incorporate an additional systematic offset of ~1° into the uncertainties that we report.

The astrometric extraction of the SPHERE H-band data was performed at the SPHERE Data Center (Delorme et al. 2017; Galicher et al. 2018) by applying the TLOCI reduction algorithm (Marois et al. 2010; Galicher et al. 2018) to the intensity data and injecting negative point sources to null the source signal. While the astrometric calibration of SPHERE is

<footnote n="24">²⁴ https://archive.stsci.edu/doi/resolve/resolve.html?doi=10.17909/T9W89V</footnote>
<table><caption>**Table 6**
Point Source Properties</caption><thead><tr><th>Date</th><th>Instrument</th><th>Filter</th><th>Separation (′′)</th><th>P.A. (deg.)</th><th>S/N</th></tr></thead><tbody><tr><td>1998 Dec 2</td><td>NICMOS</td><td>F110W</td><td>3.51 ± 0.04</td><td>90.7 ± 0.6</td><td>2.5</td></tr><tr><td>1998 Dec 2</td><td>NICMOS</td><td>F160W</td><td>3.51 ± 0.04</td><td>90.8 ± 0.6</td><td>5.6</td></tr><tr><td>2005<br>Nov 12</td><td>CIAO</td><td>H</td><td>3.449 ± 0.020</td><td>87.85<sup>+0.33</sup><sub>-1.11</sub></td><td>6.9</td></tr><tr><td>2019<br>Dec 20</td><td>SPHERE</td><td>BB_H</td><td>3.418 ± 0.008</td><td>86.0 ± 0.2</td><td>61.9</td></tr></tbody></table>

very stable (Maire et al. 2016), we still used a dedicated astrometric calibration epoch with the cluster 47 Tuc, which was acquired within three weeks of the observations of DO Tau. The initial TLOCI extraction yielded an uncertainty in the P.A. of 0.1°. We then accounted for the uncertainty in the true north alignment of the detector, which is also 0.1°.

The astrometric measurements for the different epochs are given in Table 6. The position angles of the point source are plotted in Figure 19 and compared to the expected P.A. of a bound source on a face-on circular orbit (i.e., the orientation at which the change in P.A. would be maximized) and of a stationary (distant) background source. The P.A. of the latter is calculated given a proper motion for DO Tau of $\mu_{\alpha}^* = 6.27 \pm 0.04$ mas yr$^{-1}$, $\mu_{\delta} = -21.05 \pm 0.03$ mas yr$^{-1}$ (Gaia Collaboration et al. 2021). The change in P.A. of the observed point source is much larger than expected for a bound source but smaller than for a stationary background source. Based on the SPHERE and NICMOS observations (and excluding the more uncertain Subaru observations), we

estimate that the proper motion difference between the point source and DO Tau is $\Delta\mu_{\alpha}^* = -4.7 \pm 2.0$ mas yr$^{-1}$, $\Delta\mu_{\delta} = 13.6 \pm 2.4$ mas yr$^{-1}$. To examine whether the object could be a background source that is still within Taurus, we retrieved Gaia DR2 proper motions of the Taurus members listed in Luhman (2018) and computed their dispersions, which were 3.5 mas yr$^{-1}$ in the R.A. direction and 4.0 mas yr$^{-1}$ in decl. Because the relative proper motion of the point source deviates by more than $3\sigma$ in the decl. direction from the dispersion in Taurus, we conclude that the point source is most likely not a member of Taurus, but not distant enough from the observer to appear stationary.

## Appendix E
Channel Maps of Molecular Line Observations

Channel maps of $^{12}$CO J = 2–1, $^{13}$CO J = 2–1, C$^{18}$O J = 2–1, and CS J = 5–4 toward DO Tau are provided in Figure 20.
<figure><img src="image_26.png" /><figcaption>**Figure 20.** Channel maps of $^{12}$CO J = 2–1 emission toward DO Tau, part 1. The LSRK velocity (km s$^{-1}$) appears in the top right of each panel. The synthesized beam is shown as a solid white ellipse in the lower-left corner. A gray cross marks the position of the millimeter continuum peak, which is located at the phase center. Offsets from the phase center (in arcseconds) are marked in the lower-left panel. North is up, and east is to the left. An arcinsh stretch is applied to the color scale to highlight faint emission features.</figcaption></figure>

(The complete figure set (8 images) is available.)
## ORCID iDs

Jane Huang https://orcid.org/0000-0001-6947-6072

Christian Ginski https://orcid.org/0000-0002-4438-1971

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Bin Ren (任彬) https://orcid.org/0000-0003-1698-9696

Alexander J. Bohn https://orcid.org/0000-0003-1401-9952

Élodie Choquet https://orcid.org/0000-0002-9173-0740

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

Álvaro Ribas https://orcid.org/0000-0003-3133-3580

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Til Birnstiel https://orcid.org/0000-0002-1899-8783

Yann Boehler https://orcid.org/0000-0002-8692-8744

Stefano Facchini https://orcid.org/0000-0003-4689-2684

Daniel Harsono https://orcid.org/0000-0001-6307-4195

Michiel Hogerheijde https://orcid.org/0000-0001-5217-537X

Feng Long https://orcid.org/0000-0002-7607-719X

Carlo F. Manara https://orcid.org/0000-0003-3562-262X

François Ménard https://orcid.org/0000-0002-1637-7393

Paola Pinilla https://orcid.org/0000-0001-8764-1780

Christophe Pinte https://orcid.org/0000-0001-5907-5179

Christian Rab https://orcid.org/0000-0003-1817-6576

Jonathan P. Williams https://orcid.org/0000-0001-5058-695X

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

## References

Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21

Akiyama, E., Vorobyov, E. I., Baobabu Liu, H., et al. 2019, AJ, 157, 165

Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, ApJL, 904, L6

Andreani, P., Carpenter, J., Diaz-Trigo, M., et al. 2016, ALMA Cycle 4 Proposer's Guide and Capabilities, ALMA Doc. 4.2 v1.0, https://arc.iram.fr/documents/cycle4/ALMA_04_proposer_guide.pdf

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157

Arce, H. G., & Goodman, A. A. 2002, ApJ, 575, 911

Artur de la Villarmois, E., Jørgensen, J. K., Kristensen, L. E., et al. 2019, A&A, 626, A71

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

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

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Bae, J., Hartmann, L., & Zhu, Z. 2015, ApJ, 805, 15

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

Bally, J., Cunningham, N. J., Moeckel, N., et al. 2011, ApJ, 727, 113

Bate, M. R. 2018, MNRAS, 475, 5618

Batygin, K. 2012, Natur, 491, 418

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020, ApJ, 898, 97

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Bonnell, I., & Bastien, P. 1992, ApJL, 401, L31

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

Braun, T. A. M., Yen, H.-W., Koch, P. M., et al. 2021, ApJ, 908, 46

Cabrit, S., Pety, J., Pesenti, N., & Dougados, C. 2006, A&A, 452, 897

Carbillet, M., Bendjoya, P., Abe, L., et al. 2011, ExA, 30, 39

Choquet, É., Pueyo, L., Hagan, J. B., et al. 2014, Proc. SPIE, 9143, 914357

Choquet, É., Pueyo, L., Soummer, R., et al. 2015, Proc. SPIE, 9605, 96051P

Clarke, C. J. 2020, MNRAS, 491, L72

Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155

Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

Deo, F., Facchini, S., Clarke, C. J., & Haworth, T. J. 2015, MNRAS, 449, 1996 de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63

Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proc. of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé et al., 347

Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477

Dong, R., Liu, H. B., Cuello, N., et al. 2022, NatAs, 6, 331

Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112

Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075

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

Erkal, J., Dougados, C., Coffey, D., et al. 2021, A&A, 650, A46

Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459

Facchini, S., Manara, C. F., Schneider, P. C., et al. 2016, A&A, 596, A38

Fernández-López, M., Zapata, L. A., Rodríguez, L. F., et al. 2020, AJ, 159, 171

Finkbeiner, D. P., Padmanabhan, N., Schlegel, D. J., et al. 2004, AJ, 128, 2577

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

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 451

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Fukagawa, M., Hayashi, M., Tamura, M., et al. 2004, ApJL, 605, L53

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92

Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A62

Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94

Garufi, A., Dominik, C., Ginski, C., et al. 2022a, A&A, 658, A137

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

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

Ginski, C., Ménard, F., Rab, C., et al. 2020, A&A, 642, A119

Ginski, C., Schmidt, T. O. B., Mugrauer, M., et al. 2014, MNRAS, 444, 2280

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29

Grady, C. A. 2004, in ASP Conf. Ser., 321, Extrasolar Planets: Today and Tomorrow, ed. J. Beaulieu, A. Lecavelier Des Etangs, & C. Terquem (San Francisco, CA: ASP), 244

Grady, C. A., Woodgate, B., Bruhweiler, F. C., et al. 1999, ApJL, 523, L151

Güdel, M., Eibensteiner, C., Dionatos, O., et al. 2018, A&A, 620, L1

Guerri, G., Daban, J.-B., Robbe-Dubois, S., et al. 2011, ExA, 30, 59

Gutermuth, R. A., Megeath, S. T., Pipher, J. L., et al. 2005, ApJ, 632, 397

Hagan, J. B., Choquet, É., Soummer, R., & Vigan, A. 2018, AJ, 155, 179

Harsono, D., Alexander, R. D., & Levin, Y. 2011, MNRAS, 413, 423

Herbig, G. H., & Bell, K. R. 1988, in Third Catalog of Emission-Line Stars of the Orion Population, ed. G. H. Herbig & K. R. Bell (Santa Cruz, CA: Lick Observatory), 90

Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97

Hermes, D. 2017, JOSS, 2, 267

Hirth, G. A., Mundt, R., Solf, J., & Ray, T. P. 1994, ApJL, 427, L99

Howard, C. D., Sandell, G., Vacca, W. D., et al. 2013, ApJ, 776, 21

Huang, J., Andrews, S. M., Öberg, K. I., et al. 2020, ApJ, 898, 140

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, ApJS, 257, 19

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

Itoh, Y., Hayashi, M., Tamura, M., et al. 2008, PASJ, 60, 223

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

Kepley, A. A., Tsutsumi, T., Brogan, C. L., et al. 2020, PASP, 132, 024505

Koutoulaki, M., Facchini, S., Manara, C. F., et al. 2019, A&A, 625, A49

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Krist, J. E., Stapelfeldt, K. R., Ménard, F., Padgett, D. L., & Burrows, C. J. 2000, ApJ, 538, 793

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

Kuffmeier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161

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

Kuffmeier, M., Goicovic, F. G., & Dullemond, C. P. 2020, A&A, 633, A3

Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7

Kurtovic, N. T., Pérez, L. M., Benisty, M., et al. 2018, ApJ, 869, L44

Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92

Kwon, W., Looney, L. W., Mundy, L. G., & Welch, W. J. 2015, ApJ, 808, 102

Lai, D. 2014, MNRAS, 440, 3532

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Le Gal, R., Öberg, K. I., Loomis, R. A., Pegues, J., & Bergner, J. B. 2019, ApJ, 876, 72

Lebouteiller, V., Barry, D. J., Goes, C., et al. 2015, ApJS, 218, 21

Lebouteiller, V., Barry, D. J., Spoon, H. W. W., et al. 2011, ApJS, 196, 8

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJL, 917, L10

Lee, A. T., Offner, S. S. R., Kratter, K. M., Smullen, R. A., & Li, P. S. 2019, ApJ, 887, 232

Liu, H. B., Takami, M., Kudo, T., et al. 2016, SciA, 2, e1500875

Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49
Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120

Luhman, K. L. 2018, AJ, 156, 271

Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111

Lupton, R., Blanton, M. R., Fekete, G., et al. 2004, PASP, 116, 133

Magakian, T. Y. 2003, A&A, 399, 141

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

McGlynn, T., Scollick, K., & White, N. 1998, in IAU Symp. 179, New Horizons from Multi-Wavelength Sky Surveys, ed. B. J. McLean et al. (Dordrecht: Kluwer), 465

McGroarty, F., & Ray, T. P. 2004, A&A, 420, 975

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

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

Mesa, D., Ginski, C., Gratton, R., et al. 2022, A&A, 658, A63

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

Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69

Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, JMoSt, 742, 215

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49

Najita, J. R., Carr, J. S., Pontoppidan, K. M., et al. 2013, ApJ, 766, 134

Nakajima, T., & Golimowski, D. A. 1995, AJ, 109, 1181

Nealon, R., Cuello, N., & Alexander, R. 2020, MNRAS, 491, 4108

Offner, S. S. R., Dunham, M. M., Lee, K. I., Arce, H. G., & Fielding, D. B. 2016, ApJL, 827, L11

Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32

Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2020, ApJ, 890, 142

Pfalzner, S., Tackenberg, J., & Steinhausen, M. 2008, A&A, 487, L45

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Ren, B., Pueyo, L., Perrin, M. D., Debes, J. H., & Choquet, É. 2017, Proc. SPIE, 10400, 1040021

Ren, B., Wang, J. J., & Pueyo, L. 2019, centerRadon: Center determination code in stellar images, Astrophysics Source Code Library, ascl.: 1906.021

Ribas, Á., Espaillat, C. C., Macías, E., et al. 2017, ApJ, 849, 63

Rigon, L., Scholz, A., Anderson, D., & West, R. 2017, MNRAS, 465, 3889

Rodriguez, J. E., Loomis, R., Cabrit, S., et al. 2018, ApJ, 859, 150

Rodriguez, J. E., Reed, P. A., Siverd, R. J., et al. 2016, AJ, 151, 29

Sakai, N., Sakai, T., Hirota, T., et al. 2014, Natur, 507, 78

Schneider, G., Silverstone, M. D., Stobie, E., Rhee, J. H., & Hines, D. C. 2010, in Hubble after SM4. Preparing JWST, ed. S. Deustua & C. Oliveira, 15

Segura-Cox, D. M., Schmiedeke, A., Pineda, J. E., et al. 2020, Natur, 586, 228

Sheehan, P. D., & Eisner, J. A. 2018, ApJ, 857, 18

Simon, M., Holfeltz, S. T., & Taff, L. G. 1996, ApJ, 469, 890

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Spalding, C., & Batygin, K. 2016, ApJ, 830, 5

Stapelfeldt, K. R., Burrows, C. J., Krist, J. E., et al. 1995, ApJ, 449, 888

Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70

Struve, O., & Straka, W. C. 1962, PASP, 74, 474

Sturm, J. A., McClure, M. K., Harsono, D., et al. 2022, arXiv:2201.04089

Tang, Y. W., Guilloteau, S., Piétu, V., et al. 2012, A&A, 547, A84

van der Velden, E. 2020, JOSS, 5, 2004

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

van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020b, IRDAP: SPHERE-IRDIS polarimetric data reduction pipeline, Astrophysics Source Code Library, ascl:2004.015

van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, Proc. SPIE, 10400, 1040015

van't Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

Vorobyov, E. I., & Basu, S. 2005, ApJL, 633, L137

Wilking, B. A., Lada, C. J., & Young, E. T. 1989, ApJ, 340, 823

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Wilson, T. L. 1999, RPPh, 62, 143

Winter, A. J., Booth, R. A., & Clarke, C. J. 2018, MNRAS, 479, 5522

Woitas, J., & Leinert, C. 1998, A&A, 338, 122

Yen, H.-W., Takakuwa, S., Chu, Y.-H., et al. 2017, A&A, 608, A134

York, D. G., Adelman, J., Anderson, J. E. J., et al. 2000, AJ, 120, 1579

Zapata, L. A., Rodríguez, L. F., Fernández-López, M., et al. 2020, ApJ, 896, 132

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhang, Z. 1997, Image and Vision Comput., 15, 59

Zhu, Z. 2019, MNRAS, 483, 4221

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