# Misaligned Rotations of the Envelope, Outflow, and Disks in the Multiple Protostellar System of VLA 1623-2417: FAUST. III

Satoshi Ohashi, Claudio Codella, Nami Sakai, Claire J. Chandler, Cecilia Ceccarelli, Felipe Alves, Davide Fedele, Tomoyuki Hanawa, Aurora Durán, Cécile Favre, et al.

► To cite this version:

Satoshi Ohashi, Claudio Codella, Nami Sakai, Claire J. Chandler, Cecilia Ceccarelli, et al. Misaligned Rotations of the Envelope, Outflow, and Disks in the Multiple Protostellar System of VLA 1623-2417: FAUST. III. The Astrophysical Journal, 2022, 927, 10.3847/1538-4357/ac4cae. insu-03672047

HAL Id: insu-03672047

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

Submitted on 20 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# Misaligned Rotations of the Envelope, Outflow, and Disks in the Multiple Protostellar System of VLA 1623-2417: FAUST. III

Satoshi Ohashi¹, Claudio Codella²,³, Nami Sakai¹, Claire J. Chandler⁴, Cecilia Ceccarelli³, Felipe Alves⁵, Davide Fedele², Tomoyuki Hanawa⁶, Aurora Durán⁷, Cécile Favre³, Ana López-Sepulcre³,⁸, Laurent Loinard⁷,⁹, Seyma Mercimek²,¹⁰, Nadia M. Murillo¹, Linda Podio², Yichen Zhang¹, Yuri Aikawa¹¹, Nadia Balucani¹², Eleonora Bianchi³, Mathilde Bouvier³, Gemma Busquet³, Paola Caselli⁵, Emmanuel Caux¹³, Steven Charnley¹⁴, Spandan Choudhury⁵, Nicolas Cuello³, Marta De Simone³, Francois Dulieu¹⁵, Lucy Evans²,¹³, Siyi Feng¹⁶, Francesco Fontani²,⁵, Logan Francis¹⁷,¹⁸, Tetsuya Hama¹⁹,²⁰, Eric Herbst²¹, Shingo Hirano¹¹, Tomoya Hirota²², Muneaki Imai²³, Andrea Isella²⁴, Izaskun Jímenez-Serra²⁵, Doug Johnstone¹⁷,¹⁸, Claudine Kahane³, Romane Le Gal²⁶, Bertrand Lefloch³, Luke T. Maud²⁷, Maria Jose Maureira⁵, Francois Menard³, Anna Miotello²⁷, George Moellenbrock⁴, Shoji Mori¹¹, Riouhei Nakatani¹, Hideko Nomura²², Yasuhiro Oba²⁸, Ross O'Donoghue²⁹, Yuki Okoda²³, Juan Ospina-Zamudio³, Yoko Oya²³,³⁰, Jaime Pineda⁵, Albert Rimola³¹, Takeshi Sakai³², Dominique Segura-Cox⁵, Yancy Shirley³³, Brian Svoboda⁴,³⁴, Vianney Taquet², Leonardo Testi²,²⁷, Charlotte Vastel¹³, Serena Viti²⁹, Naoki Watanabe²⁸, Yoshimasa Watanabe³⁵, Arezu Witzel³, Ci Xue²¹, Bo Zhao⁵, and Satoshi Yamamoto²³,³⁰

¹ RIKEN Cluster for Pioneering Research, 2-1, Hirosawa, Wako-shi, Saitama 351-0198, Japan; satoshi.ohashi@riken.jp

² INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy

³ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁴ National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA

⁵ Center for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik (MPE), Gießenbachstr. 1, D-85741 Garching, Germany

⁶ Center for Frontier Science, Chiba University, I-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan

⁷ Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, A.P. 3-72 (Xangari), 8701, Morelia, Mexico

⁸ Institut de Radioastronomie Millimétrique, F-38406 Saint-Martin d'Hères, France

⁹ Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad Universitaria, A.P. 70-264, Cuidad de México 04510, Mexico

¹⁰ Università degli Studi di Firenze, Dipartimento di Fisica e Astronomia, via G. Sansone 1, I-50019 Sesto Fiorentino, Italy

¹¹ Department of Astronomy, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

¹² Department of Chemistry, Biology, and Biotechnology, The University of Perugia, Via Elce di Sotto 8, I-06123 Perugia, Italy

¹³ IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France

¹⁴ Astrochemistry Laboratory, Code 691, NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA

¹⁵ CY Cergy Paris Université, Sorbonne Université, Observatoire de Paris, PSL University, CNRS, LERMA, F-95000, Cergy, France

¹⁶ Department of Astronomy Xiamen University, Xiamen, 361005, China

¹⁷ NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada

¹⁸ Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 5C2, Canada

¹⁹ Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

²⁰ Department of Basic Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

²¹ Department of Chemistry, University of Virginia, McCormick Road, P.O. Box 400319, Charlottesville, VA 22904, USA

²² National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka-shi, Tokyo 181-8588, Japan

²³ Department of Physics, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

²⁴ Department of Physics and Astronomy, Rice University, 6100 Main Street, MS-108, Houston, TX 77005, USA

²⁵ Centro de Astrobiología (CSIC-INTA), Ctra. de Torrejón a Ajalvir, km 4, E-28850, Torrejón de Ardoz, Spain

²⁶ IRAP, Université de Toulouse, CNRS, UPS, CNES, F-31400 Toulouse, France

²⁷ European Southern Observatory, Karl-Schwarzschild Str. 2, D-85748 Garching bei München, Germany

²⁸ Institute of Low Temperature Science, Hokkaido University, N19W8, Kita-ku, Sapporo, Hokkaido 060-0819, Japan

²⁹ Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

³⁰ Research Center for the Early Universe, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

³¹ Departament de Química, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Spain

³² Graduate School of Informatics and Engineering, The University of Electro-Communications, Chofu, Tokyo 182-8585, Japan

³³ Steward Observatory, 933 N. Cherry Ave., Tucson, AZ 85721, USA

³⁴ Jansky Fellow of the National Radio Astronomy Observatory

³⁵ Materials Science and Engineering, College of Engineering, Shibaura Institute of Technology, 3-7-5 Toyosu, Koto-ku, Tokyo 135-8548, Japan

Received 2021 July 30; revised 2021 December 23; accepted 2022 January 13; published 2022 March 4

## Abstract

We report a study of the low-mass Class 0 multiple system VLA 1623AB in the Ophiuchus star-forming region, using H¹³CO⁺ (J = 3-2), CS (J = 5-4), and CCH (N = 3-2) lines as part of the ALMA Large Program FAUST. The analysis of the velocity fields revealed the rotation motion in the envelope and the velocity gradients in the outflows (about 2000 au down to 50 au). We further investigated the rotation of the circumbinary VLA 1623A disk, as well as the VLA 1623B disk. We found that the minor axis of the circumbinary disk of VLA 1623A is misaligned by about 12° with respect to the large-scale outflow and the rotation axis of the envelope. In contrast,

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
the minor axis of the circumbinary disk is parallel to the large-scale magnetic field according to previous dust polarization observations, suggesting that the misalignment may be caused by the different directions of the envelope rotation and the magnetic field. If the velocity gradient of the outflow is caused by rotation, the outflow has a constant angular momentum and the launching radius is estimated to be 5–16 au, although it cannot be ruled out that the velocity gradient is driven by entrainments of the two high-velocity outflows. Furthermore, we detected for the first time a velocity gradient associated with rotation toward the VLA 16293B disk. The velocity gradient is opposite to the one from the large-scale envelope, outflow, and circumbinary disk. The origin of its opposite gradient is also discussed.

*Unified Astronomy Thesaurus concepts: Star formation* (1569); *Stellar jets* (1607); *Protostars* (1302); *Circumstellar disks* (235); *Young stellar objects* (1834); *Interstellar medium* (847); *Interstellar molecules* (849); *Astronomy chemistry* (75)

# 1. Introduction

Angular momentum plays important roles in star formation processes (e.g., Bodenheimer 1995). Star formation proceeds via gravitational collapse of a dense core with a typical size of ~0.1 pc (Shu 1977). Because of the angular momentum of the dense core, the accretion material cannot directly fall onto the central protostar, which results in the formation of an accretion disk. The accretion disk evolves into a protoplanetary disk as the envelope dissipates. Planets are expected to form in such a protoplanetary disk.

During gravitational collapse, collimated supersonic (~100 km s⁻¹) jets and wider, slower (~10 km s⁻¹) outflows remove the angular momentum from the disk, allowing the accretion to the protostar (e.g., Snell et al. 1980; Frank et al. 2014). Recent ALMA high-resolution observations have succeeded in measuring the outflow rotation and have estimated the launching radii of the outflows (e.g., Bjerkeli et al. 2016; Alves et al. 2017; Hirota et al. 2017; Lee et al. 2018; Zhang et al. 2018). These observations suggest that the low-velocity outflows launch from a radius of about 2–25 au in the disk, although in some cases they could be launched from larger distances beyond the disk edge and coinciding with the location of the centrifugal barrier (see, e.g., BHB07-11; Alves et al. 2017). A relatively large launching radius for outflows seems to be consistent with a disk wind model (Konigl & Pudritz 2000) rather than an X wind model (Shu et al. 2000). However, it is still unclear how these outflows are related to the angular momentum of a parent core and disk formation.

The magnetic field is also crucial to launching the outflow, in addition to the rotation motion. The standard picture of protostar formation shows that the *B*-field couples with the gas component during the collapse of the protostellar core and becomes parallel to the outflow direction. In contrast, dust polarization observations have shown that the magnetic field (*B*-field) in protostellar cores does not correlate with the outflow but is randomly oriented (Hull et al. 2014; Hull & Zhang 2019). The angular momentum is also not aligned with the *B*-field (Galametz et al. 2020).

By taking into account the misalignment of the *B*-field and the angular momentum, ideal and nonideal magnetohydrodynamic (MHD) simulations have been performed (e.g., Matsumoto & Tomisaka 2004; Machida et al. 2006; Hennebelle & Ciardi 2009; Joos et al. 2012; Zhao et al. 2018; Hirano & Machida 2019). Hirano & Machida (2019) showed that the misalignment of the *B*-field and angular momentum causes different directions of an outflow and jet driven from different scales of a circumstellar disk. Furthermore, the circumstellar disk can be warped, as revealed by the recent discoveries of the warped protostellar disks

of HH 212 (Codella et al. 2014), L1527 (Sakai et al. 2019), and L1489 (Sai et al. 2020).

To understand the star formation processes, in particular the role of angular momentum, we investigate the rotation motion on various scales from the envelope to the disks in the Very Large Array (VLA) 1623–2417 star-forming system. VLA 1623–2417 (hereafter VLA 1623) is a deeply embedded protostellar system in Ophiuchus A (André et al. 1990) at a distance of 131 ± 1 pc (Gagné et al. 2018).

VLA 1623 is considered a prototypical Class 0 source (André et al. 1993). A prominent and complex outflow is observed in CO and H₂ emission (André et al. 1990; Dent et al. 1995; Yu & Chernin 1997; Caratti o Garatti et al. 2006). Observations at centimeter wavelengths resolved four compact sources, VLA 1623A, VLA 1623B (Knot A), VLA 1623W (Knot B), and Knot C (Bontemps & André 1997). The nature of VLA 1623A is proposed as the driving source of the outflow and an asymmetric jet composed of VLA 1623B and VLA 1623W (Bontemps & André 1997; Maury et al. 2012; Hara et al. 2021), or as protostars (Looney et al. 2000; Murillo & Lai 2013). The presence of multiple outflows was highlighted in early work (Dent et al. 1995; Caratti o Garatti et al. 2006), suggesting that the second scenario was more likely. Based on the spectral energy distribution (SED) for each component, including 0.42 mm (ALMA Band 9), Herschel Space Observatory PACS fluxes and Spitzer IRS spectra, all three sources are confirmed to be protostars, with VLA 1623A and VLA 1623B being Class 0, while VLA 1623W is Class I (Murillo et al. 2018).

Recent high angular resolution dust continuum observations with ALMA revealed the complicated structure of VLA 1623, including the binarity of VLA 1623A (Harris et al. 2018). Figure 1 shows the 0.87 mm dust continuum image obtained from the ALMA archive. Details of the observations are described in Harris et al. (2018). The VLA 1623A binary, with components A1 and A2 separated by ~28 au, is associated with a circumbinary disk. It is unclear whether the prominent outflow from VLA 1623A is driven by one or both of the sources in the binary. VLA 1623B is a very bright source showing a disk-like structure, with a projected separation from VLA 1623A of ~150 au to the west. A compact but strongly collimated outflow driven by VLA 1623B was detected with ALMA Cycle 0 observations (Santangelo et al. 2015). VLA 1623B is characterized by presenting a variable H₂O maser (Furuya et al. 2003) and SO emission (Murillo et al. 2013; Hsieh et al. 2020), but lacking emission in CO isotopologues or DCO⁺ (Murillo et al. 2013, 2015). Modeling of the C¹⁸O emission around VLA 1623A found a rotationally supported disk with $r_{disk} \sim 200$ au, showing Keplerian rotation out to at least 160 au and an unresolved central mass of ~0.2 M☉
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Images of the 0.87 mm dust continuum emission observed toward the VLA 1623 region with ALMA. The data are taken from the ALMA archive, and details of the observation are described by Harris et al. (2018). The Class 0 source, VLA 1623A, is resolved into two sources, VLA 1623A1 and VLA 1623A2. The beam size is shown in the lower left corner of the right panel. The contours are 0.93, 2.7, 10, and 50 mJy beam⁻¹ (1σ = 0.12 mJy beam⁻¹) to indicate the circumbinary disk around VLA 1623A and circumstellar disk around VLA 1623B.</figcaption></figure>

(Murillo et al. 2013). Subsequently, Hsieh et al. (2020) derived a total binary mass of 0.3–0.5 $M_\odot$ using better spatial resolution ALMA data. Given the large size of the VLA 1623A disk and no apparent distortions at the edge (e.g., disk truncation or heating), the separation between VLA 1623A and VLA 1623B is likely larger than 200 au. Due to these interesting features of the early stage and its multiplicity, VLA 1623 was selected as a target of the ALMA Large Program FAUST (Fifty AU STudy of the chemistry in the disk/envelope system of solar-like protostars<fnref n="36" />). This program aims at revealing the physical and chemical structure of 13 nearby protostars ($d = 137–235$ pc), at scales from a few 1000 au down to 50 au, by observing various molecular lines.

In this paper, we present a subset of the VLA 1623 results from the FAUST project based on the analysis of the velocity maps of the H¹³CO⁺ ($J=3–2$), CS ($J=5–4$), and CCH ($N=3–2$) emission. These molecular lines trace the envelope, outflow, and disk structures of VLA 1623. In Section 2 we summarize our observations. Section 3 shows the results and describes the structures of the envelope, outflow, circumbinary disk, and circumstellar disk. We investigate the rotation motion of the envelope and the outflow and show the misalignment of their rotation axes on different scales in Section 4. Then, we discuss the origin of the misalignments in Section 5. Finally, we summarize the paper in Section 6.

## 2. Observations

Single-field observations for VLA 1623 were carried out between 2018 November and 2019 April as part of the ALMA Large Program FAUST. The molecular lines in Band 6 analyzed in this study are listed in Table 1. We used the 12 m array data from two different configurations (C43-5 and C43-2 for sparse and compact configurations, respectively) and the 7 m array data of the Atacama Compact Array (ACA/Morita Array), combining these visibility data in the UV plane during the imaging. In total, the baseline lengths range from

7.2 to 737.4 m. The adopted field center was ($\alpha_{2000}$, $\delta_{2000}$) = (16h26m26s.392, -24°24'30"718), which is close to the protostellar position. The back-end correlator for the molecular line observations was set to a resolution of 122 kHz (~0.17 km s⁻¹) and a bandwidth of 62.5 MHz (~87 km s⁻¹). The data were reduced in Common Astronomy Software Applications package (CASA) 5.6.1–8 (McMullin et al. 2007) using a modified version of the ALMA calibration pipeline and an additional in-house calibration routine to correct for the $T_{\text{sys}}$ and spectral line data normalization.<fnref n="37" /> Self-calibration was carried out using line-free continuum emission for each configuration. The complex gain corrections derived from the self-calibration were then applied to all channels in the data, and the continuum model derived from the self-calibration was subtracted from the data to produce continuum-subtracted line data. A self-calibration technique was also used to align both amplitudes and phases across the multiple configurations. Images were prepared by using the *tclean* task in CASA 5.6.1, where Briggs weighting with a robustness parameter of 0.5 was employed. The primary beam correction was applied to all the images presented in this paper. The maximum recoverable scale ($\theta_{\text{MRS}}$) is estimated as ~0.6λ/L<sub>min</sub>, where λ is the observing wavelength and L<sub>min</sub> is the minimum baseline. Since the maximum recoverable scales are ~20", any structures extended more than that size will be resolved out. The rms noise levels for H¹³CO⁺, CS, and CCH are 2.5, 2.3, and 2.8 mJy beam⁻¹ channel⁻¹, respectively. The achieved synthetic beams are shown in Table 1.

To understand the relationship between the molecular gas distribution and the protostars, we compare the 0.87 mm continuum image previously reported by Harris et al. (2018) with the molecular line distributions. The absolute position in the 0.87 mm observation is fixed by comparing the VLA 1623B position of the present FAUST observations at 1.3 mm with the previous 0.87 mm continuum observations. We find

<footnote n="36">³⁶ http://faust-alma.riken.jp</footnote>

<footnote n="37">³⁷ https://help.almascence.org/kb/articles/what-errors-could-originate-from-the-correlator-spectral-normalization-and-tsys-calibration</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** A three-color image of the integrated intensities of H<sup>13</sup>CO<sup>+</sup> (blue), CCH (green), and CS (red) emission toward the VLA 1623 region. The velocity range for integration is $V_{LSR} = -6$ to 10 km s<sup>-1</sup>. The white contour shows the 0.87 mm dust continuum emission presented by Harris et al. (2018) and shown in Figure 1 to indicate the disk and protostellar positions.</figcaption></figure>

<table><caption>**Table 1**
List of Observed Lines<sup>a</sup></caption><thead><tr><th>Molecule<sup>a</sup></th><th>Transition</th><th>Frequency (GHz)</th><th>S<sub>μ</sub><sup>2</sup> (D<sup>2</sup>)</th><th>E<sub>ν</sub>k<sup>-1</sup> (K)</th><th>Beam Size</th><th>rms (mJy beam<sup>-1</sup> channel<sup>-1</sup>)</th></tr></thead><tbody><tr><td>H<sup>13</sup>CO<sup>+</sup></td><td>J = 3–2</td><td>260.25534</td><td>45.6</td><td>25.0</td><td>0.′′498 × 0.′′469 (P.A. 80°)</td><td>2.5</td></tr><tr><td>CS</td><td>J = 5–4</td><td>244.93556</td><td>19.2</td><td>35.3</td><td>0.′′529 × 0.′′495 (P.A. 82°)</td><td>2.3</td></tr><tr><td>CCH</td><td>N = 3–2, J = 7/2–5/2, F = 4–3</td><td>262.00426</td><td>2.3</td><td>25.2</td><td>0.′′493 × 0.′′465 (P.A. 81°)</td><td>2.8</td></tr><tr><td>CCH</td><td>N = 3–2, J = 7/2–5/2, F = 3–2</td><td>262.00648</td><td>1.7</td><td>25.1</td><td>0.′′493 × 0.′′465 (P.A. 81°)</td><td>2.8</td></tr></tbody></table>

**Note.**

that the peak position of VLA 1623B in the 0.87 mm continuum image is offset by ~0.′′1 from that in the 1.3 mm continuum image. This is plausibly due to the proper motion of the source, which will be discussed in a forthcoming study. Even though the offset is smaller than our spatial resolution of ~0.′′5, it is corrected by shifting the 0.87 mm continuum image in order to compare it with the molecular distributions obtained in this study. The offset in VLA 1623A is hardly identified because our spatial resolution is not enough to resolve the binary protostars of A1 and A2.

### 3. Results

In this section, we show the intensity distributions of the molecular line emission. Figure 2 shows a three-color image of the integrated intensities of H<sup>13</sup>CO<sup>+</sup> (blue), CCH (green), and CS (red) emission overlaid on the 0.87 mm dust continuum emission. The velocity range for integration is set to be $V_{LSR} = -6$ to 10 km s<sup>-1</sup>. The detailed distributions of each molecular line are described in Section 3.2. Based on Figure 2, the molecular lines are classified into two groups: the dense gas

tracer H<sup>13</sup>CO<sup>+</sup> (e.g., Caselli et al. 2002; Onishi et al. 2002; Ikeda et al. 2007), and the outflow cavity tracers CS and CCH (e.g., Zhang et al. 2018) in this source. The CS and CCH molecules are expected to be abundant in the (UV-irradiated) warm region even though the critical densities ($n_{cr}$) of these molecular lines are almost similar ($n_{cr} \sim 10^7$ cm<sup>-3</sup>).

#### 3.1. Spectra of the Molecular Lines toward the Protostars

Figures 3 and 4 show the spectral line profiles toward the protostellar positions of VLA 1623A1 and VLA 1623B. Note that our observations cannot resolve the positions of components A1 and A2. The CCH molecule has hyperfine structure, with two transitions separated by 2.2 MHz, corresponding to ~2.5 km s<sup>-1</sup>. The LSR velocity scale is defined relative to the lower frequency ($N = 3 - 2$, $J = 7/2 - 5/2$, $F = 4 - 3$) CCH transition.

The most prominent feature of these spectra is recognized as absorptions at the systemic cloud velocity of $V_{LSR} = 3.8$ km s<sup>-1</sup> (Narayanan & Logan 2006). In particular, the absorption feature of the H<sup>13</sup>CO<sup>+</sup> and CCH lines is found toward each

<footnote><sup>a</sup> Line parameters are taken from CDMS (Endres et al. 2016).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** H<sup>13</sup>CO<sup>+</sup> (J = 3–2) and CS (J = 5–4) spectra toward the protostellar positions of VLA 1623A1 and VLA 1623B are shown as the black and red lines, respectively. Vertical dashed lines stand for the source systemic velocity (3.8 km s<sup>-1</sup>).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** CCH (N = 3–2, J = 7/2–5/2, F = 4–3) and (N = 3–2, J = 7/2–5/2, F = 3–2) spectra toward the protostellar positions of VLA 1623A1 and VLA 1623B are shown. Vertical dashed lines stand for the source systemic velocity (3.8 km s<sup>-1</sup>) relative to the F = 4–3 and F = 3–2 hyperfine components. V<sub>LSR</sub> is calculated by using the F = 4–3 hyperfine component.</figcaption></figure>

protostar position, suggesting that the lines are absorbed against the optically thick continuum background of the protostellar disk(s) and/or self-absorption caused by the foreground (infalling) gas. These absorption features are also observed in the channel maps of the H<sup>13</sup>CO<sup>+</sup>, CS, and CCH emission (see the following subsections). We discuss the origin of the absorption feature in Section 4.1.

When we assume that the absorption originates from the optically thick background emission, the systemic velocity in front of these protostars is determined to be 3.8 km s<sup>-1</sup>, which is consistent with that of the parent protostellar core (Narayanan & Logan 2006). Even if the absorption is due to the self-absorption or infall motion, the systemic velocity would not be much different from 3.8 km s<sup>-1</sup>.

Another interesting feature is found in the CS line toward VLA 1623B (Figure 3, bottom panel). The CS emission is detected in a wider velocity range of –6 to 12 km s<sup>-1</sup>, while the H<sup>13</sup>CO<sup>+</sup> and CCH emission is detected in a narrower velocity range of 2.5–5 km s<sup>-1</sup> (per hyperfine component). We

discuss these CS high-velocity components in terms of a rotation motion of the VLA 1623B disk in Section 4.5.

## 3.2. The Velocity-integrated Intensity Maps

In this section, we describe the distributions of the H<sup>13</sup>CO<sup>+</sup>, CS, and CCH emission from the velocity-integrated intensity maps.

### 3.2.1. The Envelope Structure Shown by the H<sup>13</sup>CO<sup>+</sup> Emission

Figure 5 shows the integrated intensity map of the H<sup>13</sup>CO<sup>+</sup> emission overlaid on the 0.87 mm dust continuum emission tracing the circumbinary disk of VLA 1623A and the protostellar disk around VLA 1623B. The H<sup>13</sup>CO<sup>+</sup> emission is widely detected and shows an elongated structure along the northeast-to-southwest direction.

The size of the H<sup>13</sup>CO<sup>+</sup> distribution is about 4000 au, corresponding to 0.02 pc, which is similar to the typical size of compact dense cores in various star-forming regions, including high-mass star-forming regions in infrared dark clouds (IRDCs) (e.g., Ohashi et al. 2016). Therefore, the distribution of the H<sup>13</sup>CO<sup>+</sup> emission appears to represent the dense core in the VLA 1623 region.

Figure 5 (right panel) also shows the zoomed-in image around the circumbinary disk and protostars. Lower integrated intensities are found toward the protostars of VLA 1623A1, VLA 1623A2, and VLA 1623B. In particular, an absorption feature (~-0.05 Jy beam<sup>-1</sup>, corresponding to ~4 K at V<sub>LSR</sub> = 3.8 km s<sup>-1</sup>) is found toward VLA 1623B as shown in Figure 3.

### 3.2.2. The Outflow Cavity Structure Shown by the CS and CCH Emission

Figure 6 shows the integrated intensity maps of the CS and CCH emission. The color image indicates the CS distribution, while the black contours indicate the CCH emission distribution. The blue contours show the 0.87 mm dust continuum emission to indicate the circumbinary disk and protostellar positions. The zoom-in view around the protostars is also shown in the right panel.

The distributions of CS and CCH are similar to each other and show “X-shaped” structures, indicating the edges of the outflow cavities. The outflow cavity structure is elongated in the southeast-to-northwest direction, which is consistent with the direction of the CO outflow (Santangelo et al. 2015; Hara et al. 2021). We discuss the outflow properties in Section 4.4.

The distributions of CS and CCH are locally different in some parts, even though the CS and CCH distributions reveal similar outflow cavity structures at large scales. The most prominent difference is seen in the disk around VLA 1623B (see right panel of Figure 6). The integrated intensity of the CS emission is peaked at this position, while that of the CCH emission shows a hole structure. A similar feature of the CCH intensity depression near the protostar is shown by Oya et al. (2017) in L483. They suggested that the depression of the CCH emission is caused by the gas-phase destruction and/or depletion onto dust grains. We suggest that the CS emission in VLA 1623B comes from the high-velocity components of a rotation motion (Figure 3; see Section 4.5 for further discussion). In contrast, the CCH emission is not detected in the high-velocity regions and is absorbed at the systemic
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Integrated intensity maps of the H<sup>13</sup>CO<sup>+</sup> emission toward the VLA 1623 region. The velocity range for integration is 1.4–6.0 km s<sup>-1</sup>. The 1σ noise level is 2.9 × 10<sup>-3</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup>. The black contours show the 0.87 mm dust continuum emission presented by Harris et al. (2018) and shown in Figure 1, to indicate the disk and protostellar positions. The beam size of the H<sup>13</sup>CO<sup>+</sup> emission is shown in the lower left corner of each panel.</figcaption></figure>

velocity ($V_{LSR} \sim 3.8$ km s<sup>-1</sup>). Note that the CS emission also shows the absorption feature around the systemic velocity.

Furthermore, other different distributions are also found in the outflow cavity structures. For example, the CS outflow cavity structure is slightly wider compared to the CCH outflow cavity in the northern part of the southeast outflow. In addition, the CS emission is detected in the surrounding regions perpendicular to the outflow direction.

To conclude, based on the reported spatial distributions, we suggest that the CS emission traces not only the outflow cavity structure but also the disk rotation and accretion flow as similar to the binary Class 0 source of IRAS 16293 A1/A2 (Maureira et al. 2020). In contrast, the CCH emission mainly shows the outflow cavity structure.

### 3.3. Analysis of the Channel Maps

To investigate the kinematic motions of the envelope, outflow, and disks, we show the channel maps of the H<sup>13</sup>CO<sup>+</sup>, CS, and CCH emission. The velocity range for the integration and channel step are set to be 0.6 km s<sup>-1</sup> in each channel map.

#### 3.3.1. H<sup>13</sup>CO<sup>+</sup> Channel Maps

Figure 7 shows channel maps of the H<sup>13</sup>CO<sup>+</sup> emission to illustrate several components showing the different structures.

In the velocity range <2.1 km s<sup>-1</sup>, the H<sup>13</sup>CO<sup>+</sup> emission shows a filamentary structure to the south of the protostars. A similar filamentary structure is also identified by Hsieh et al. (2020) in SO emission in the same velocity range. Hsieh et al. (2020) interpret the SO emission as tracing an accretion flow toward the VLA 1623A and VLA 1623B disks.

In the channel maps with $V_{LSR} = 3.3–4.5$ km s<sup>-1</sup>, the H<sup>13</sup>CO<sup>+</sup> emission is widely distributed, tracing the envelope structure. The blueshifted components are located on the northeast side, while the red components are located on the southwest side, indicating a rotation motion of the envelope.

Figure 8 shows a zoomed-in view around the protostars. In the channel maps of $V_{LSR} = 2.1–2.7$ km s<sup>-1</sup> and $V_{LSR} = 5.1–5.7$ km s<sup>-1</sup>, the H<sup>13</sup>CO<sup>+</sup> emission is detected around the circumbinary disk. In particular, the velocity components ($V_{LSR} = 2.1$ and $5.7$ km s<sup>-1</sup>) are detected inside the circumbinary disk. The velocity gradient is found to be consistent with the envelope rotation, indicating the rotation of the circumbinary disk. Therefore, we suggest that H<sup>13</sup>CO<sup>+</sup> would be a good tracer for probing the kinematics of the disk (e.g., Tobin et al. 2011).

The absorption feature is found toward the VLA 1623B disk at the velocity of 3.3–3.9 km s<sup>-1</sup>. In addition, the H<sup>13</sup>CO<sup>+</sup> emission is depressed toward the protostars of VLA 1623A in the circumbinary disk. This drop in emission would be caused by the optically thick continuum emission of the background sources of the protostellar disks and/or less abundance of the H<sup>13</sup>CO<sup>+</sup> molecule.

#### 3.3.2. CS Channel Maps

Figure 9 shows 0.6 km s<sup>-1</sup> channel maps of the CS emission, as for the H<sup>13</sup>CO<sup>+</sup> emission. As shown in the integrated intensity maps of the CS and CCH emission (Figure 6), the CS emission mainly traces the outflow cavity structure. We also indicate the directions of the CO outflows observed by Hara et al. (2021) as the blue and red arrows. We find that the outflow cavity encompasses these two high-velocity outflows.

In the channel maps centered at $V_{LSR} = 0.9$ and $1.5$ km s<sup>-1</sup>, the CS emission shows an extension perpendicular to the outflow. In particular, the CS emission is found to the north of the circumbinary disk. This would be part of the accretion flow toward the circumbinary disk.

In the channel maps from $V_{LSR} = 2.1$ to $5.1$ km s<sup>-1</sup>, the CS emission mainly shows the outflow cavity structure. In the channel maps at $V_{LSR} = 2.1$ and $2.7$ km s<sup>-1</sup>, the northern part of the outflow is detected. In contrast, the southern part of the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Integrated intensity maps of the CS emission in color scale and the CCH emission in black contours toward the VLA 1623 region. The velocity range for integration is -6 to +10 km s⁻¹. The contours of the CCH emission are 0.02, 0.05, 0.1, and 0.2 Jy beam⁻¹ km s⁻¹, where 1σ = 0.0025 Jy beam⁻¹ km s⁻¹. The blue thick contours show the 0.87 mm dust continuum emission (Harris et al. 2018) to indicate the circumbinary disk and protostellar positions. The beam size of the CS emission is shown in the lower left corner of each panel.</figcaption></figure>

outflow is detected in the channel maps at $V_{LSR} = 4.5$ and
$5.1$ km s⁻¹.

The “X-shaped” outflow cavity structure is found in the channel map of $V_{LSR} = 3.9$ km s⁻¹, which indicates that the systemic velocity of the outflow is ~3.9 km s⁻¹. Therefore, the systemic velocities of the envelope and outflow are consistent, to within the spectral resolution of the channel maps. In addition, the velocity gradient along the outflow was not identified.

No emission is found toward the VLA 1623B disk in the channel maps at $V_{LSR} = 2.1-4.5$ km s⁻¹ as in the case for the $H^{13}CO^+$ emission. In particular, the absorption is seen in the outflow cavity structure shown in the channel maps at $V_{LSR} = 3.3-5.1$ km s⁻¹, indicating that the CS emission is absorbed against the optically thick emission of the background continuum source of the VLA 1623B disk. In other words, the background dust emission is optically thick and masks the line emission. This means that the outflow is located in front of the VLA 1623B disk.

### 3.3.3. CCH Channel Maps

Figure 10 shows the channel maps of the CCH emission. The CCH molecule has hyperfine structure with the two transitions of CCH separated by 2.2 MHz (see Table 1), corresponding to ~2.5 km s⁻¹. Therefore, the same distribution is repeated in the channel maps at an interval of 2.5 km s⁻¹.

The CCH channel maps show the clear outflow cavity structure. The blueshifted emission can be seen in the northern part of the outflow cavity in $V_{LSR} = 0.3-0.9$ km s⁻¹ and $V_{LSR} = 2.7-3.3$ km s⁻¹, while the redshifted emission appears in the southern part of the outflow cavity in $V_{LSR} = 2.1$ km s⁻¹ and $V_{LSR} = 4.5-5.7$ km s⁻¹. This velocity gradient is consistent with the CS channel maps, indicating that the velocity gradient is perpendicular to the outflow lobes. However, the accretion flow or additional components surrounding the protostars

identified in the CS emission are not observed in the CCH emission. Thus, the CCH emission traces the whole outflow cavity structure, while the CS emission traces not only the outflow cavity but also the disk components and accretion flows.

## 4. Discussion

### 4.1. The Origin of the Absorption Feature toward the Protostellar Positions

We found the absorption features at ~3.8 km s⁻¹ toward the protostellar positions shown in the spectra (Figures 3 and 4) in the channel maps (Figures 7, 9, and 10). If these are caused by self-absorption, the molecular line emission needs to be highly optically thick to produce the absorption feature, irrelevant to the optical depth of the dust continuum emission. However, we found that the $H^{13}CO^+$ and CCH emission is optically thin with $\tau \sim 0.1-0.3$ if the continuum background is assumed to be optically thin because the peak brightness temperature of the emission is ~2–10 K in the protostellar positions. We assume that the excitation temperature is consistent with the dust temperature of ~40 K, which is derived by the SED fitting toward the VLA 1623A and VLA 1623B positions (Murillo et al. 2018). Note that the brightness temperatures of the 0.87 mm continuum emission reach as high as ~60 K at those protostellar positions. If we assume the brightness temperature of 60 K as the excitation temperature, the optical depths of the molecular line emission become much lower. Even if we assume an excitation temperature of 10 K, the $H^{13}CO^+$ and CCH emission is still optically thin with $\tau \sim 0.3$. Furthermore, observations with the single-dish telescope of the Caltech Submillimeter Observatory (CSO) showed that the $H^{13}CO^+$ (J = 4–3) emission is optically thin in this region (Narayanan & Logan 2006). Therefore, it is likely that the lines are absorbed against the optically thick background emission of the protostellar disks even though the self-absorption and infall
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Channel maps of the H<sup>13</sup>CO<sup>+</sup> (J = 3–2) emission, overlaid with the 0.87 mm dust continuum emission (Harris et al. 2018) as black contours. The contour is 0.93 mJy beam<sup>-1</sup>, indicating the disks of VLA 1623A, B, and W. The channel width and step are 0.6 km s<sup>-1</sup>, and the channel centroid velocities are labeled at the upper center of each panel. The 1σ noise level for H<sup>13</sup>CO<sup>+</sup> is 1.2 × 10<sup>-3</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup> in each panel. The beam size is shown in the lower left corners.</figcaption></figure>

motions may also affect the line profiles. This can be confirmed by the spatial distributions of these lines. The depressions of the H<sup>13</sup>CO<sup>+</sup> emission toward VLA 1623B are found across a wide range of velocities, as seen in the channel maps. Therefore, it is likely that VLA 1623B is highly embedded within its envelope.

Alternatively, the depression toward the protostars of VLA 1623 A1 and A2 could be the result of absorption due to the background emission from the circumstellar disks around these protostars. In addition, it might also be possible that the depression is caused by lower abundance of the HCO<sup>+</sup> molecule in the dense regions of the circumstellar disks due to a decrease of the ionization degree.

## 4.2. The Rotation Axes of the Envelope and Circumbinary Disk

In this subsection, we investigate the rotation motions from the large scales of the envelope and outflow to the small scales of the circumbinary disk and the circumstellar disk around VLA 1623AB, and we derive their rotation axes. First, we show the rotation axes of the envelope (~1000 au) and circumbinary disk (~100 au) using the H<sup>13</sup>CO<sup>+</sup> emission. Then, we investigate the outflow direction and rotation by using both CS and CCH emission.

The intensity-weighted mean velocity (moment 1) maps of the H<sup>13</sup>CO<sup>+</sup> emission are shown in Figure 11. The mean velocity is derived by using a velocity range of 1.2–6.2 km s<sup>-1</sup> with an intensity of >0.015 Jy beam<sup>-1</sup> (corresponding to >6σ). The wide velocity range of 1.2–6.2 km s<sup>-1</sup> includes the high-velocity components of the circumbinary disk rotation in addition to the slow rotation of the envelope. The overall velocity structure is characterized by blueshifted emission on the northeast side of the circumbinary disk and redshifted emission on the southwest side with respect to the systemic velocity of the envelope (3.8 km s<sup>-1</sup>). This velocity gradient in

the direction of the large-scale elongation of the H<sup>13</sup>CO<sup>+</sup> emission is interpreted as rotation of the envelope.

Interestingly, Figure 11 seems to show different rotation axes on different scales of the envelope and circumbinary disk because the velocity pattern is twisted from the envelope to the circumbinary disk. The zoomed-in view of Figure 11 shows that the velocity gradient in the circumbinary disk is parallel to the major axis of the circumbinary disk, indicating that the rotation motion corresponds to that of the disk. The magenta dashed lines indicate a direction with a position angle of 115°, the minor axis of the circumbinary disk, which is derived from Gaussian fitting (Harris et al. 2018). The magenta dashed lines agree with the rotation direction of the circumbinary disk, which is tilted with respect to the envelope rotation axis.

To measure the rotation axis of the envelope, we use the mean velocity map made with a velocity range of 3.2–4.4 km s<sup>-1</sup>, which is narrower than that used in Figure 11 because the velocities around the systemic velocity better represent the envelope kinematics without contamination of the rotation of the circumbinary disk. The mean velocity map of the envelope is shown in Figure 12.

To measure the rotation axis, we perform a linear fit to the positions whose mean velocity is 3.8 ± 0.05 km s<sup>-1</sup>, assuming that the systemic velocity should be on the rotation axis. Furthermore, the fitting area is narrowed to the range from +2′′6 to −6′′0 in R.A. direction and from +4′′5 to −4′′0 in decl. direction from VLA 1623A1 to avoid the contamination from the infalling components. By linear fitting, the rotation axis of the envelope is derived to be a position angle of ~127° ± 1° as shown by the black dotted lines in Figure 12. Even if we select the mean velocity range of 3.8 ± 0.2 km s<sup>-1</sup>, we obtain a similar value of ~124° ± 1°. After these examinations, we conclude that the rotation axis of the envelope is ~127° ± 4°. Figure 12 also plots the rotation axis
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The same channel maps of the H<sup>13</sup>CO<sup>+</sup> (J = 3–2) emission as Figure 7, but zoomed-in views around the protostars. The black contours indicate the 0.87 mm dust continuum emission (Harris et al. 2018) to indicate the disks and stellar positions. The contour levels are the same as in Figure 1. The beam size is shown in the lower left corners.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Channel maps of the CS emission, overlaid with the 0.87 mm continuum emission (Harris et al. 2018) as black contours. The contour level is 0.93 mJy beam<sup>-1</sup>, indicating the disks of VLA 1623A, B, and W. The blue and red arrows indicate the directions of the CO outflows observed by Hara et al. (2021). The channel width and step are 0.6 km s<sup>-1</sup>, and the channel centroid velocities are labeled at the upper center of each panel. The 1σ noise level for each panel is 0.0012 Jy beam<sup>-1</sup> km s<sup>-1</sup>. The beam size is shown in the lower left corners.</figcaption></figure>

of the circumbinary disk in the magenta dashed lines. By comparing these two rotation axes, we confirm that the rotation axes of the circumbinary disk and envelope are misaligned by ~12° ± 6°.

Note that the derived rotation axis of the envelope changes slightly if we adopt a different mean velocity range and fitting area. However, the misalignment rotation axes between the envelope and circumbinary disk should be robust because the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Channel maps of the CCH (N = 3–2, J = 7/2–5/2, F = 4–3) and (N = 3–2, J = 7/2–5/2, F = 3–2) emission, overlaid on the 0.87 mm continuum emission (Harris et al. 2018) as black contours. The contour level is 0.93 mJy beam⁻¹, indicating the disks of VLA 1623A and VLA 1623B. The channel width and step are 0.6 km s⁻¹, and the channel centroid LSR velocities relative to the rest frequency of the N = 3–2, J = 7/2–5/2, F = 4–3 transition are labeled at the upper center of each panel. The 1σ noise level for each panel is 0.0014 Jy beam⁻¹ km s⁻¹. The beam size is shown in the lower left corner of each panel.</figcaption></figure>

mean velocity field shows the twisted pattern from the envelope to the circumbinary disk as shown in Figure 11. Such a twisted pattern of the velocity field can also be reproduced by infalling −rotating motion on ~100 au scales (e.g., Yen et al. 2013; Oya et al. 2016). However, the twisted pattern of this source is seen on a larger scale of ~1000 au. Therefore, we suggest that this velocity pattern indicates misalignment.

To investigate the kinematics of the envelope, we show the position–velocity (PV) diagrams of the H¹³CO⁺ emission, along with the major and minor axes of the envelope rotation, in Figure 13. The PV diagram along the major axis (Figure 13(b)) shows broad velocity components (V<sub>LSR</sub> ~2–6 km s⁻¹) within 1" from the protostars. These high-velocity components trace the rotation motion of the circumbinary disk. In addition, a velocity gradient on a large scale is found with an offset range from −20" to 20". This velocity gradient indicates the rotation motion of the envelope. The PV diagram along the minor axis (Figure 13(c)) also contains a high-velocity component in the blueshifted emission (V<sub>LSR</sub> ~2–3 km s⁻¹) near the protostars. This high-velocity component traces the rotation motion of the circumbinary disk because the minor axis is slightly tilted from the rotation axis of the circumbinary disk. Contrary to the PV diagram of the major axis, a velocity gradient along the minor axis is not as significant as in the

major axis. The emission is almost at the systemic velocity of 3.8 km s⁻¹.

As shown in Figure 13, the large velocity width is only seen near the protostars. Except for this part, the velocity width (Δv) on the envelope scale is only ~0.3–0.5 km s⁻¹. Furthermore, the velocity gradient is seen along the major axis. Therefore, we suggest that turbulence has only a minor effect on the velocity field of the envelope.

### 4.3. Is the Outflow Cavity Created by a Single Outflow or Multiple Outflows?

Before investigating the outflow direction and rotation, we discuss the launching source(s) of the CS and CCH outflow cavity. The integrated intensity maps of the CS and CCH emission show the simple “X-shaped” outflow cavity, but the several high-velocity outflows are detected by the CO observations by Santangelo et al. (2015) and Hara et al. (2021). These CO outflows show wider velocity ranges (V<sub>LSR</sub> = −34.5 to −0.5 km s⁻¹ and V<sub>LSR</sub> = +6.5 to 31.5 km s⁻¹) than the CS and CCH outflow cavity. Thus, the CO emission traces the high-velocity outflows, while the CS and CCH emission traces the low-velocity outflow cavity. Santangelo et al. (2015) reported that the high-velocity outflows are launched from VLA 1623A and VLA 1623B, and Hara et al. (2021) suggested that there are two outflows launched from A1
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Intensity-weighted mean velocity (moment 1) maps of the H$^{13}$CO$^{+}$ line toward the VLA 1623 region. The right panel is the zoomed-in view around the protostars. The mean velocity is derived using a velocity range of 1.2–6.2 km s$^{-1}$ with an intensity of >0.015 Jy beam$^{-1}$ (6σ). The black contour indicates the 0.87 mm dust continuum emission (Harris et al. 2018) with an intensity of 0.93 mJy beam$^{-1}$. Magenta dashed lines indicate the minor axis of the circumbinary disk. The protostellar positions are shown as the star markers. The beam size is shown in the lower left corners.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Same as Figure 11, but for a velocity range of 3.2–4.4 km s$^{-1}$ with an intensity of >0.015 Jy beam$^{-1}$. The right panel is the zoomed-in view around the protostars. Magenta dashed and black dotted lines indicate the rotation axes of the circumbinary disk and envelope, respectively. The beam size is shown in the lower left corners.</figcaption></figure>

and A2. Thus, the three protostars of VLA 1623A1, A2, and B are suggested to have the high-velocity outflows.

The CS and CCH channel maps (Figures 9 and 10) indicate that the emission has a velocity gradient from north to south centered around the systemic velocity of 3.8 km s$^{-1}$. Furthermore, the channel maps of the outflow cavity indicate a symmetric distribution around the systematic velocity of 3.8 km s$^{-1}$. In particular, the emission distributions around the systematic velocity show the clear “X-shaped” structure. This velocity gradient would be naturally explained by the rotation of a single outflow. If material from the circumbinary disk is only accreting onto the circumstellar disk of VLA 1623A1, the low-velocity outflow would only be launched from the outer edge of the VLA 1623A1 circumstellar disk. In the disk wind

model, the low- and high-velocity outflows are driven near the disk outer and inner edges, respectively (Machida 2014). In contrast, the high-velocity outflows can be launched from the inner edges of the disks owing to the accretion toward the protostars. The low-velocity outflow has a wider opening angle than the high-velocity outflows. Therefore, the CO high-velocity outflows are encompassed by the CS and CCH emission, which is consistent with the recent MHD simulation of a binary formation (Saiki & Machida 2020).

The other possibility is that the two independent outflows may produce such a velocity gradient of the cavity structure. Similar to the rotation scenario, Hara et al. (2021) showed that the blueshifted outflows are located in the north and the redshifted outflows are located in the south (Figure 9). These
<figure><img src="image_16.png" /><figcaption>**Figure 13.** (a) Moment 1 map from Figure 12, and PV diagrams of the H<sup>13</sup>CO<sup>+</sup> (J = 3–2) emission along (b) the major axis and (c) the minor axis. The shadows shown in panels (b) and (c) indicate the circumbinary disk region. Contours are 0.01, 0.05, and 0.1 Jy beam<sup>-1</sup>. Horizontal dashed lines show the systemic velocity (V<sub>sys</sub> = 3.8 km s<sup>-1</sup>).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Same channel maps as Figure 9, but zoomed-in views around the protostars. The contour levels are the same as in Figure 1. The channel width and step are 0.4 km s<sup>-1</sup>, and the channel centroid velocities are labeled in the upper right corner of each panel. The dashed lines indicate the outflow cavity structures seen at the channel maps of V<sub>LSR</sub> = 2.4 and 4.8 km s<sup>-1</sup>. The beam size is shown in the lower left corners.</figcaption></figure>

high-velocity outflows are encompassed by the CS and CCH outflow cavity. Therefore, it would be possible that the velocity gradient of the outflow cavity is produced by the entrainment of these high-velocity outflows.

Our current observations cannot conclude whether the velocity gradient of the outflow is caused by rotation or the entrainment of the high-velocity outflows. However, we discuss the outflow properties such as the outflow axis, angular momentum, and launching radius by assuming the outflow rotation scenario in Section 4.4.

Note that the observed velocity gradient would not be explained by the precession motion of the binary protostars because an “S-shaped” wiggling structure is not identified in the channel maps. If the outflow has a precession motion and produces the observed velocity gradient, a curved trajectory is

expected (e.g., Fendt & Zinnecker 1998). However, the channel maps indicate that the outflow cavity structure runs from north to south without any wiggling structures and with increasing velocity.

To find the launching position in detail, we reanalyzed the CS data by CASA task *clean* with a robust parameter of -2. The image is obtained with a slightly better spatial resolution of 0.′′474 × 0.′′435. Figure 14 shows the zoomed-in views of the CS channel maps around the protostars with a better spatial resolution. The velocity gradient can be seen in the vicinity of the protostars. Even though the spatial resolution of the CS emission is not enough to distinguish sources A1 and A2, the outflow seems to be launched in the vicinity of source A1. This can be seen in the channels at V<sub>LSR</sub> = 2.4, 4.4, and 4.8 km s<sup>-1</sup>. For example, the outflow cavity structures seen at V<sub>LSR</sub> = 2.4
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Integrated intensity maps of the (a) CS and (b) CCH emission overlaid with the 0.87 mm dust continuum emission (Harris et al. 2018) in blue thick contours. The velocity range for integration is 3.4–4.2 km s⁻¹. Thin black contours of the CS emission start from 5σ and have intervals of 15σ, where 1σ = 1.2 mJy beam⁻¹ km s⁻¹. Thin black contours of the CCH emission start from 5σ and have intervals of 7σ, where 1σ = 1.4 mJy beam⁻¹ km s⁻¹. Thick black contours of the CS and CCH emission encompass the outflow cavity structure and are 25σ and 7σ, respectively. The 0.87 mm dust continuum contours are at a level of 0.93 mJy beam⁻¹. The magenta dashed, black dotted, and gray solid lines indicate the minor axis of the circumbinary disk, rotation axis of the envelope, and outflow direction, respectively. The beam size is shown in the lower left corners.</figcaption></figure>

and 4.8 km s⁻¹ cross at the position of A1 (see the dashed lines in the channel maps of $V_{LSR} = 2.4$ and 4.8 km s⁻¹). However, our spatial resolution is not high enough to determine the launching source definitively. Further investigations with higher spatial resolution are needed to identify the launching position(s) of the CS and CCH outflow cavity and to distinguish the outflow rotation from the entrainment of the high-velocity outflows.

## 4.4. Outflow Axis, Angular Momentum, and Launching Radius

In the previous subsection, we discussed the outflow launching source. From the channel maps of the CS and CCH emission, we suggested a possibility that the velocity gradient is caused by the outflow rotation launched from VLA 1623A1. Based on this assumption, we investigate the outflow properties of the axis, angular momentum, and launching radius. We note that it may also be possible that the CS and CCH outflow cavity is caused by the entrainment of the two high-velocity outflows rather than rotation. In that case, the discussion of this subsection would need to be reinvestigated.

### 4.4.1. Outflow Axis

To measure the outflow direction, we use the integrated intensity maps of the CS and CCH emission shown in Figure 15. The velocity range for the integration is selected to be $V_{LSR} = 3.4–4.2$ km s⁻¹ around the systemic velocity of the outflow (~3.8 km s⁻¹).

The thick black contours in Figure 15 show the 25σ and 7σ noise levels of the CS and CCH integrated intensities, respectively. These contours encompass the outflow cavity structure. First, we derive the directions of the outflow cavity

edges (north and south) by performing the linear fitting of these contours. Then, we derive the outflow direction by averaging the position angles of these two edges.

We focus only on the southeast lobe, because the northwest outflow is affected by absorption against the VLA 1623B disk. In addition, the fitting area is narrowed to within 6.2′′ (corresponding to 850 au) from source A1 because the contours (25σ of CS and 7σ of CCH) do not follow the outflow structure outside of this area.

By fitting the contours, the position angles of the outflow cavity edges (north and south directions) are derived to be ~84° ± 0.4° and ~166° ± 1° for the CS outflow and ~91° ± 0.5° and ~167° ± 1° for the CCH outflow. Therefore, we derive the outflow directions of ~125° ± 1° for the CS outflow and ~129° ± 1° for the CCH outflow. These outflow directions are close to each other. Furthermore, the outflow direction of ~125°–129° is also consistent with that of the CO high-velocity outflow of 128° (outflow-2 launched from A1 or A2) measured by Hara et al. (2021). We note that the opening angle of the CS outflow cavity looks slightly wider than the CCH outflow cavity since the CS emission traces not only the outflow cavity structure but also the accretion flows. It also causes the slight difference of the outflow directions between the CS and CCH emission.

The outflow directions are shown with the thick gray lines in Figure 15. We also plot the rotation axis of the envelope and minor axis of the circumbinary disk in the magenta dashed and black dotted lines, respectively. Note that the minor axis of the circumbinary disk can be regarded as the rotation axis of the circumbinary disk. These directions indicate that the rotation axis of the envelope is consistent with the outflow direction, while the rotation axis of the circumbinary disk is slightly
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Mean velocity (moment 1) maps of the (a) CS and (b) CCH emission. The mean velocity of the CS emission is derived by using a velocity range of 0.6–5.4 km s⁻¹ with an intensity of >0.05 Jy beam⁻¹, while that of the CCH emission is derived by using a velocity range of 2.6–6.2 km s⁻¹ with an intensity of >0.02 Jy beam⁻¹ to avoid the contamination from the hyperfine emission. Magenta contours indicate 25σ and 7σ of the CS and CCH integrated intensity emission, respectively. Black lines indicate the outflow directions. The dashed lines perpendicular to the outflow axis indicate the cuts for the PV diagrams shown in Figures 17 and 18. The colors of the cuts across the eastern lobe are used in Figure 18. (b) Same as panel (a), but for the CCH emission. The beam size is shown in the lower left corners.</figcaption></figure>

different from the others. We discuss possible origins of these different rotation axes in Section 5.

### 4.4.2. Rotation and Launching Radius of the Outflow

In this subsection, we investigate outflow properties by analyzing the angular momentum and the launching radius of the outflow. The mean velocity maps of the CS and CCH emission are shown in Figure 16. We confirm that the velocity gradient is found to be perpendicular to the outflow direction. The north part of the outflow is blueshifted, while the south part is redshifted. We suggested that this velocity gradient indicates the outflow rotation in Section 4.3. However, we should note that there might be the possibility that the velocity gradient is produced by the two independent outflows.

The magenta contours indicate the 25σ and 7σ noise levels of the CS and CCH emission, respectively, encompassing the outflow cavity structure. The velocity field of the CCH emission shows a gradient within this contour. In contrast, the velocity field of the CS emission shows the gradient even outside the 25σ contour. This may indicate that the velocity gradient shown in the CS emission may be contaminated by the accretion flows or envelope materials because the velocity gradient in the CS is in the same direction as the rotation of the envelope. Therefore, we focus only on the CCH emission in the following analysis.

To investigate the rotation of the outflow, we analyze the CCH emission in a similar way to Zhang et al. (2018). We show the PV diagrams of the CCH emission in the southeast lobe and the base of the northwest lobe, along the lines perpendicular to the outflow axis (the lines are shown in Figure 16). We set the center at the outflow axis, and the

negative value is the northern part of the outflow cavity, while the positive value is the southern part of the outflow cavity. The outflow axis (z) passes through the VLA 1623A1 protostar and has a position angle of 129°. The lines of the PV diagrams are taken for every 0′′6 from $z = -1′′2$ to $+6′′0$ except for $z = 0$.

The PV diagrams in Figure 17 indicate that the higher-velocity components are found closer to the protostar and the emission has a wider spatial distribution at farther distances from the protostar. These features are consistent with the outflow from NGC 1333 IRAS 4C (Zhang et al. 2018). A red curve shown in each panel indicates a constant angular momentum of 100 au km s⁻¹ with respect to the outflow axis. The CCH emission seems to follow this curve, suggesting constant angular momentum of the outflow. Note that the PV diagrams with $z \ge 2′′4$ do not show any sign of higher-velocity components around the outflow axis position (offset ~0′′0). This is because the CCH emission traces the outflow cavity wall. Even if the CCH emission is observed apparently on the outflow axis, the emission region is the cavity wall at the near and far sides of the outflow axis. The slow rotation velocity of the cavity wall is due to a relatively large linear distance from the outflow axis. The median velocity of 3.8 km s⁻¹ (the horizontal dashed lines) does not change with z from the systemic velocity of the envelope, indicating a nearly edge-on configuration. We suggest that the CCH emission traces the walls of the outflow cavity instead of the inner part of the outflow. Thus, high-velocity components at large distances from the central source cannot be seen by the CCH emission.

To show the outflow rotation more clearly, the contours of the PV diagrams for the southeast lobe ($z \ge 0′′6$) are overlaid to form the composite PV diagram shown in Figure 18. The colors depict the position of the cut lines shown in Figure 16.
<figure><img src="image_20.png" /><figcaption>**Figure 17.** PV diagrams of the CCH ($N = 3-2, J = 5/2-3/2$) emission along the lines perpendicular to the outflow axis (shown in Figure 16). The distances of these lines to the central source are labeled in the upper left corner of each panel ($z > 0$ for the southeast lobe and $z < 0$ for the northwest lobe). The contours are 6, 20, 80, and 100 mJy beam$^{-1}$. The rest frequency of the $F = 4-3$ hyperfine line is used, so that the $F = 3-2$ hyperfine line appears at a blueshifted velocity. The red curves correspond to a constant angular momentum of 100 au km s$^{-1}$ (see Section 4.4.2). The black bar in the lower right corner of each panel indicates the beam size.</figcaption></figure>

The contours are set to be $0.6 \times I_{\text{peak}}$, where $I_{\text{peak}}$ is the maximum intensity in each panel. The contours follow the expected curve for constant angular momentum of 100 au km s$^{-1}$. Interestingly, the angular momentum of 100 au km s$^{-1}$ is quite similar to that for the NGC 1333 IRAS 4C outflow

(Zhang et al. 2018). Because both VLA 1623A1 and NGC 1333 IRAS 4C are low-mass (~0.2 $M_{\odot}$) protostars in the Class 0 stage, it is not surprising that these sources have similar outflow properties. However, the rotation feature on the PV diagrams in this outflow is not as clear as that of NGC 1333
<figure><img src="image_21.png" /><figcaption>Figure 18. PV diagram showing the emission peaks of different panels of the southeast lobe in Figure 17 in one panel. The color corresponds to the different cuts shown in Figure 16. The contours are at levels of 0.6 times the maximum intensity in each panel. The curves are the same as those in Figure 17.</figcaption></figure>

IRAS 4C. Further investigation with better spatial resolution observations will allow us to reveal the rotation motion more clearly.

Next, we estimate the launching radius of the outflow from the angular momentum. Following the method described by Anderson et al. (2003), the wind-launching radius can be constrained by the outflow rotation, the poloidal velocity of the outflow ($v_p$), and the protostellar mass. Because the outflow is aligned with the plane of sky, the poloidal velocity is difficult to measure. Therefore, we assume a typical value of $v_p = 3-10$ km s$^{-1}$, as is the case for NGC 1333 IRAS 4C. We also assume that the poloidal velocity of the outflow is much higher than the toroidal (rotation) velocity. The total mass of VLA 1623A1 and A2 is assumed to be 0.3–0.5 $M_\odot$, and VLA 1623A has been considered to be an equal-mass binary system (Harris et al. 2018; Kawabe et al. 2018). Therefore, we simply assume the half value of 0.2 $M_\odot$ for the VLA 1623A1 protostellar mass. With these assumptions, the launching radius is derived to be 5–16 au. Even if we use a higher mass of 0.4 $M_\odot$ assuming the launching point in the circumbinary disk, the launching radius is still 5–10 au. This value is consistent with the fact that the launching position is not resolved by our observations. Since the separation of the binary is 28 au, we conclude that the outflow is launched from the circumstellar disk of VLA 1623A1 (or maybe A2) rather than the circumbinary disk.

The derived launching radius (5–16 au) is comparable to other outflow launching radii such as TMC1-A (Bjerkeli et al. 2016), Orion Source I (Hirota et al. 2017), HH 212 (Tabone et al. 2017; Lee et al. 2018), HH 46/47 (Zhang et al. 2016), NGC 1333 IRAS 4C (Zhang et al. 2018), and Monoceros R2-IRS2 (Jiménez-Serra et al. 2020). The rotation in the inner jet, on the other hand, has been reported in some other sources with a launching radius of <0.1 au, such as B335 (Bjerkeli et al. 2019) and HH 212 (Lee et al. 2017). Theoretical models predict the launching of the fast jet and slow outflow at the inner disk edge and outer disk radius, respectively (Machida 2014). If this

is the case, our results are consistent with a wide-angle disk wind model launched from relatively large radii in the disk.

### 4.5. The Reversed Velocity Gradient of the Disk around VLA 1623B

Here we report the discovery of the velocity gradient of the VLA 1623B disk, which has the opposite sense to the others of the envelope, circumbinary disk, and the outflow.

As shown in the integrated intensity map of the CS emission (Figure 6), the peak emission is found toward the VLA 1623B disk owing to the contribution of the high-velocity component. Figure 19 shows channel maps of the blue- and redshifted high-velocity components of the CS emission. The velocity range for the integration and velocity step of the channels is 0.8 km s$^{-1}$. The channel maps of the blueshifted component show the velocity ranging from -6.4 to 1.6 km s$^{-1}$, while the channel maps of the redshifted component show the velocity ranging from 5.2 to 13.2 km s$^{-1}$. We confirm that the high-velocity emission is mainly detected in the VLA 1623B disk.

From the channel maps of Figure 19, we find a velocity gradient along the disk major axis. The blueshifted emission is detected in the southern part of the disk, while the redshifted emission is detected in the northern part, indicating the rotation motion. Interestingly, this velocity gradient is in the opposite sense to that of the envelope, outflow, and VLA 1623A circumbinary disk.

Figure 20 shows the mean velocity maps of the $H^{13}CO^+$ and CS emission overlaid on the 0.87 mm dust continuum image. The mean velocity map of the $H^{13}CO^+$ emission is the same as that of Figure 11, which is used for revealing the rotation motion of the circumbinary disk. We confirm that Figure 20 shows that the velocity gradient of the VLA 1623B is reverse to that of the circumbinary disk. The reversed velocity gradient is also confirmed by the PV diagram along the disk major axis shown in panel (c).

Here we roughly estimate the dynamical mass of VLA 1623B from the velocity gradient. The dynamical mass ($M_{dyn}$) is calculated as

$$M_{dyn} = \frac{r v_{rot}^2}{G} \sim 0.01 \left( \frac{r}{10 \text{ au}} \right) \left( \frac{v_{rot}}{1 \text{ km s}^{-1}} \right)^2 M_{\odot}, \quad (1)$$

where $r$ and $v_{rot}$ are the radius and rotation velocity of the disk, respectively. The disk is assumed to be edge-on. As shown in Figure 20, the velocity gradient of the rotation is roughly ~39 km s$^{-1}$ arcsec$^{-1}$ ($V_{LSR} = -5$ to $9$ km s$^{-1}$ within a $0''35$ distance). Although the disk structure is not resolved by our observation, we roughly assume a disk radius of $0''2$ (corresponding to 27 au) because the peak distance between the blue and red components is about ~$0''4$. Then, the dynamical mass ($M_{dyn}$) is derived to be ~1.7 $M_\odot$, which is much larger than the total mass of the VLA 1623A (~0.4 $M_\odot$). If the CS-emitting region is smaller than the above estimate, the derived mass may be smaller: it would be 0.2–0.4 $M_\odot$ if the radius of the emitting region were 13–17 au (corresponding to ~$0''1$). Higher spatial resolution observations are needed to derive a more precise value for the dynamical mass of VLA 1623B.

It is very surprising that the velocity gradient of the VLA 1623B disk is completely opposite to the other velocity gradients of the envelope, circumbinary disk, and outflow. The
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Channel maps of the CS emission from the high-velocity components, overlaid with the 0.87 mm dust continuum emission (Harris et al. 2018) in black contours. The channel width and step are 0.8 km s⁻¹, and the channel centroid velocities are labeled at the upper center of each panel. The beam size is shown in the lower right corners.</figcaption></figure>

position angle of the VLA 1623B disk is derived to be 41°.4 ± 0°.5 from Gaussian fitting of the dust continuum image (Harris et al. 2018). Thus, the rotation axis corresponds to 131°.4 ± 0°.5, which is similar to the outflow axis (~129°) and the envelope rotation axis (~127°). We discuss the origin of the reversed velocity gradient of the VLA 1623B disk in the next section.

Note that opposite rotations in close binary disks have recently been identified in a high-mass star-forming region (Tanaka et al. 2020), suggesting that these different rotations in close binary systems might be a common occurrence and are established in the early phases of star formation.

Although previously it has been proposed that VLA 1623B corresponds to a shocked cloudlet (Maury et al. 2012; Hara et al. 2021), the clear rotation pattern in these observations further confirms that the structure corresponds to a protostellar disk.

**5. The Origin of the Misalignment**

In the previous sections, we found that the rotation axes of the envelope and low-velocity outflow are consistent, while the rotation axis of the circumbinary disk (VLA 1623A) is slightly tilted by ~12° ± 6°. Even though the velocity gradient of the

low-velocity outflow will be explained by either rotation or entrainment of the two high-velocity outflows, the direction of the low-velocity outflow is approximately consistent with the envelope rotation axis as shown in Section 4.4.1. In addition, the velocity gradient of the VLA 1623B disk is opposite to the other velocity gradients. The several high-velocity outflows are reported by the CO observations, which would be launched from the circumstellar disks around the VLA 1623A1, A2, and B protostars (Santangelo et al. 2015; Hara et al. 2021). In contrast, we suggest that the low-velocity outflow may be launched from the outer edge of the VLA 1623A1 disk at a radius of $r \sim 5–16$ au. These high-velocity outflows and single low-velocity outflow are well aligned with the envelope rotation.

According to these results, Figure 21 illustrates a schematic view of the possible rotation motions on various spatial scales, and Table 2 lists the derived position angles of the rotations and disk inclinations.

To identify the origin of the misalignment of the circumbinary disk of VL1623A, the magnetic field direction is expected to be an important parameter. Recent ALMA dust continuum polarization observations detected the polarization in the circumbinary disk and the VLA 1623B disk (Harris et al. 2018; Sadavoy et al. 2018). Sadavoy et al. (2018) found that
<figure><img src="image_23.png" /><figcaption>**Figure 20.** (a) Intensity-weighted mean velocity (moment 1) map of the H$^{13}$CO$^{+}$ emission from Figure 11. (b) Mean velocity (moment 1) map of the CS emission toward the VLA 1623B disk. The mean velocity is derived using a velocity range of -8 to 15 km s$^{-1}$ with an intensity of >0.01 Jy beam$^{-1}$. White contours indicate the 0.87 mm dust continuum emission. (c) PV diagram of the CS emission toward the circumstellar disk of VLA 1623B. The direction is indicated as the black dashed line in panel (b). The contours are 0.01, 0.02, and 0.04 Jy beam$^{-1}$, respectively.</figcaption></figure>

the polarization morphology is well explained by a static, oblate spheroid model with a poloidal magnetic field, assuming that the dust polarization traces the magnetic fields in the disk, similar to other disks (Alves et al. 2018; Ohashi et al. 2018). Furthermore, the direction of the poloidal magnetic field is consistent with the circumbinary disk minor axis. Thus, the direction of the global magnetic field is misaligned with the rotation axis of the envelope and outflow by ~12° (Figure 21).

Based on our results, together with the previous observations, the three-dimensional structures of the disks are revealed. We note, however, that these disk inclinations are based on an interpretation of the velocity structures as purely due to outflow rotation. If the multiple protostars and multiple outflows create complex velocity structures, the three-dimensional structures shown in Figure 21 would be changed.

### 5.1. The Origin of the Misalignment of the Circumbinary Disk

The VLA 1623 region has complicated structures such as the multiple system. Our observations found that the circumbinary disk is misaligned with the envelope rotation. It may be possible that the binary motion or the third stellar component (VLA 1623B) perturbs the circumbinary disk. In this case, it is expected that the circumstellar disks and outflow direction are also perturbed by such orbital motions. The change of the outflow direction is recently reported by Okoda et al. (2021) in another young stellar object, IRAS 15398–3359. However, Hara et al. (2021) estimate a precession amplitude of 0.39 ± 0.02 for the outflow, which is much smaller than the misalignment of ~12°.

The misalignment of a disk from an envelope has been investigated with numerical simulations by taking into account the magnetic field, turbulence, and rotation of a dense core (e.g., Matsumoto & Tomisaka 2004; Hennebelle & Ciardi 2009; Li et al. 2013; Seifried et al. 2013). Recent hydrodynamical simulations suggest that a circumbinary disk and circumstellar disks around individual stars can be misaligned with each other in turbulent cores because the angular momentum of accreting material changes with evolution (Bate 2018). Thus, the orientations of the circumbinary disk and circumstellar disks could be more chaotic. Such misalignment between a circumbinary disk and circumstellar disks has recently been found by Alves et al. (2019) and Maureira et al. (2020).

However, this scenario seems to be inconsistent with the well-aligned structures of the envelope and outflow in this

source. The same rotation axes of the envelope and outflow suggest that the angular momentum of the accreting material does not change significantly with evolution. Furthermore, the turbulent motion does not appear to affect the velocity fields on scales of the envelope, as shown in Figure 13.

Nonideal MHD simulations have also investigated the outflow and disk structures in relation to the angular momentum and magnetic field directions (Matsumoto et al. 2017; Hirano & Machida 2019). Hirano & Machida (2019) investigated the relationship between outflows, rotation axes, and magnetic field directions in dense cores. They showed that low-density collapsing gas is tightly coupled to the magnetic field, resulting in the misalignment of the outer region of the disk with the rotation axis. In addition, they indicated that the inner part of the disk is well aligned with the core rotation axis because the rotation or centrifugal forces overcome the magnetic effects, due to magnetic dissipation in high-density regions. The degree of the misalignment is determined by the initial conditions. The misalignment of ~12° found by our study is comparable to these simulations. Thus, a misalignment or warped disk can be formed owing to the different directions of the rotation axis and magnetic fields.

This is indeed the case for VLA 1623 as shown in Figure 21. The circumstellar disks of the VLA 1623A1 and A2 protostars may be formed in the high-density inner region, where the magnetic field is dissipated, and these circumstellar disks are aligned with the rotation axis. Because the dissipation of the magnetic field becomes effective within $r \lesssim 10$ au (Machida et al. 2007), the launching radius of the outflow, 5–16 au, is consistent with the region where the rotation is dominant. In contrast, the circumbinary disk is formed in the lower-density outer region ($r \sim 100$ au), resulting in the circumbinary disk being aligned with the magnetic field rather than the rotation axis. This scenario can explain why the outflow launched from the circumstellar disk of VLA 1623A1 (or A2) is aligned with the rotation of the envelope rather than that of the circumbinary disk.

The misalignment between the angular momentum and magnetic field could be a reasonable scenario to explain the misalignment of the circumbinary disk. Further studies of angular momentum of the circumstellar disks of VLA 1623A1 and A2 will allow us to further constrain the detailed structures of the VLA 1623A binary system. The large-scale magnetic
<figure><img src="image_24.png" /><figcaption>**Figure 21.** A schematic view of the various rotation motions in the VLA 1623 region studied in this work.</figcaption></figure>

<table><caption>**Table 2**
List of Derived Position Angles at Various Scales</caption><thead><tr><th>Region</th><th>Probed by</th><th>Position Angle (deg)</th><th>Inclination (deg)</th><th>Scale (au)</th></tr></thead><tbody><tr><td>Envelope</td><td>H<sup>13</sup>CO<sup>+</sup></td><td>127 ± 4</td><td></td><td>~2000</td></tr><tr><td>Outflow</td><td>CS</td><td>125 ± 1</td><td></td><td>~2000</td></tr><tr><td>Outflow</td><td>CCH</td><td>129 ± 1</td><td></td><td>~2000</td></tr><tr><td>circumbinary disk<sup>a</sup><br>(VLA 1623A)</td><td>0.87 mm dust continuum</td><td>115 ± 4<sup>b</sup></td><td>~55<sup>b</sup></td><td>~300</td></tr><tr><td>circumstellar disk<sup>c</sup><br>(VLA 1623B)</td><td>0.87 mm dust continuum</td><td>131.4 ± 0.5<sup>b</sup></td><td>~74<sup>b</sup></td><td>~100</td></tr></tbody></table>

**Notes.**

field should also be investigated to reveal better the magnetic field direction in the envelope.

## 5.2. The Origin of the Opposite Velocity Gradient of the Circumstellar Disk around VLA 1623B

As discussed in the previous subsection, the VLA 1623B disk rotation has the opposite direction to the envelope, outflow, and circumbinary disk. The opposite rotations or misalignments of the angular momentum have been reported in several protostellar disks such as BHR71 (with a separation of ~3000 au; Tobin et al. 2019), IRAS 16293-2422 (with a separation of ~600 au; Zapata et al. 2013), and IRAS 04191+1523 (with a separation of ~860 au; Lee et al. 2017). These misaligned binary systems are suggested to be caused by turbulent fragmentation of their parent cores (e.g., Padoan & Nordlund 2002; Goodwin et al. 2004; Offner et al. 2010) rather than fragmentation with rotating cores (e.g., Price & Bate 2007; Boss & Keiser 2014).

One difficulty with the turbulent fragmentation for VLA 1623 is the separation between VLA 1623A and VLA 1623B.

In general, initial separations of the fragments are suggested to be a few thousand au (Offner et al. 2010; Tsukamoto & Machida 2013). In contrast, the projected distance of the VLA 1623A and VLA 1623B is ~150 au, which is much closer than the other misaligned binary systems noted. A separation of 150 au is similar to the case of the triplet system of L1448 IRS3B, which is thought to be formed by disk fragmentation via the gravitational instability (Tobin et al. 2016). However, it might be possible that the separation of VLA 1623A and VLA 1623B is much larger than 150 au. We found that VLA 1623B is located behind VLA 1623A, because the outflow emission is masked by the optically thick background emission of the VLA 1623B disk. The lack of any indication of gravitational interaction between the circumbinary disk and the VLA 1623B disk (e.g., through disturbances in the density or velocity structure at the edges of those disks) is also notable and suggests an association in projection only. Therefore, it should be noted that the separation of 150 au is a lower limit. Recent magnetohydrodynamical simulations show that a companion formed by the turbulent fragmentation can migrate

<footnote><sup>a</sup> The rotation is seen in H<sup>13</sup>CO<sup>+</sup>.</footnote>

<footnote><sup>b</sup> The position angle and inclination are derived by Harris et al. (2018).</footnote>

<footnote><sup>c</sup> The rotation is seen in CS.</footnote>
to a primary star with a separation of a ~100 au (Kuffmeier et al. 2019), which may be consistent with our results.

Another difficulty is that turbulent motions were not identified in the envelope. The typical line width of the H$^{13}$CO$^{+}$ is $\Delta v \sim 0.3-0.5$ km s$^{-1}$, and an ordered rotation pattern is clearly observed in the envelope. These observed features indicate that the turbulence is not significant in this region in comparison with other regions such as BHR71 (Tobin et al. 2019).

Nonideal MHD simulations show that the counterrotation can appear in the outflow on envelope scales owing to the Hall effect (e.g., Krasnopolsky et al. 2011; Tsukamoto et al. 2015; Wurster et al. 2016; Zhao et al. 2021). However, it is not clear whether such counterrotation can form a protostar and an accretion disk.

One possibility is that the mutual gravito-hydrodynamical interaction between the circumbinary disk and circumstellar disk disturbs the orbit of the tertiary disk, resulting in the reverse rotation (Takaishi et al. 2021). Takaishi et al. (2021) simulated the evolution of a triplet protostellar system by using three-dimensional smoothed particle hydrodynamics without magnetic fields. They found that the tertiary protostar is formed via the circumbinary disk fragmentation and the initial rotational directions of all three circumstellar disks are almost parallel to that of the orbital motion of the binary system. Then, the tertiary orbit becomes unstable owing to the three-body effect, and the circumstellar disk of the tertiary orbit receives the counterrotating gas accretion. As a result, the rotational direction reverses. This scenario seems to be consistent with our results because the close encounter or three-body interaction was suggested (Murillo et al. 2013; Harris et al. 2018).

An interaction with the outflow might also change the rotation direction of the VLA 1623B disk. The VLA 1623B disk has previously been interpreted as a knot created by outflow shock (Bontemps & André 1997; Maury et al. 2012; Hara et al. 2021) because VLA 1623B (and also source W) is in the outflow cavity. Furthermore, the CO brightness temperature is ~90 K at the position of VLA 1623B. The morphology and the local temperature enhancement indicate that the gas is heated, possibly due to a shock in the outflow. Therefore, the VLA 1623B protostar may be affected by the outflow launched from VLA 1623A1. The interaction with the outflow might have changed the disk direction even if the disk were originally aligned with the envelope and circumbinary disk. However, more detailed studies of the possible time evolution, through the proper motions and the temporal changes in the inclination of the disks, for example, are needed to confirm this scenario.

Another possible scenario might be a collision of the accretion flows from the envelope. Hsieh et al. (2020) identified several accretion flows connecting to VLA 1623B. They suggested that the blue- and redshifted accretion flows collide with each other at VLA 1623B, which may remove and/or change the angular momentum of the accretion. In this case, a disk with the opposite rotation might be formed. We need further molecular line observations with high sensitivity and high spatial resolution to investigate the accretion flows and their accretion directions in detail.

## 6. Summary

As part of the ALMA Large Program FAUST, we have investigated the kinematics of the triple protostellar system of

the VLA 1623–2427 low-mass star-forming region using H$^{13}$CO$^{+}$ ($J=3-2$), CS ($J=5-4$), and CCH ($J=7/2-5/2$) emission. We have identified the rotation motions on large scales of the envelope and outflows and on a smaller scale (50 au) of the circumbinary and circumstellar disks. The main results are listed below.

1. The H$^{13}$CO$^{+}$ emission traces the envelope and circumbinary disk structures, while the CS and CCH emission mainly traces a wide and low-velocity outflow cavity structure. In the mean velocity map of the H$^{13}$CO$^{+}$ emission we find a twisted pattern, which indicates misalignment by ~12° ± 6° between the rotation axis of the envelope and that of the circumbinary disk: the rotation axis of the envelope is ~127°, while that of the circumbinary disk is ~115°. The outflow direction is ~125° and ~129° (from CS and CCH), in agreement with the envelope and not with the circumbinary disk.

2. We find that the velocity gradient of the outflow cavity structure is traced by CS and CCH emission, which can be explained either by rotation of the low-velocity outflow or by entrainment of the two high-velocity outflows. If we assume the case of the outflow rotation, the launching radius of the low-velocity outflow is estimated to be $r \sim 5-16$ au from the protostar. This suggests that the low-velocity outflow is launched from the outer edge of the VLA 1623A1 disk by accreting materials to the disk from the circumbinary disk, while the high-velocity outflows are launched from the inner edges of the circumstellar disks. Future observations with higher resolution are needed to assess this scenario.

3. The magnetic field appears to be parallel to the rotation axis of the circumbinary disk, supporting the idea that the circumbinary disk was formed in a low-density outer region coupled with the magnetic field. In contrast, the circumstellar disk is formed in a high-density inner region. Therefore, the misalignment between the envelope rotation and magnetic field direction could cause the misalignment of the circumbinary disk. This picture is consistent with the nonideal MHD simulations studied by Hirano & Machida (2019).

4. We also find that the high-velocity CS emission traces the disk rotation around the VLA 1623B protostar. The velocity gradient of the VLA 1623B disk is opposite to those of the other components of the region. The dynamical mass of VLA 1623B is suggested to be higher than the total mass of the VLA 1623A protostars. Although detailed origins of such opposite velocity gradients are not clear and call for further investigations, the hydrodynamical interaction with the circumbinary disk and the interaction with the outflow are possible causes of the change in the disk direction.

We thank Masahiro N. Machida, Yu Saiki, and Hauyu Baobab Liu for fruitful discussions. This paper makes use of the following ALMA data set: ADS/JAO.ALMA# 2018.1.01205.L (PI: Satoshi Yamamoto). ALMA is a partnership of the ESO (representing its member states), the NSF (USA) and NINS (Japan), together with the NRC (Canada) and the NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, the AUI/NRAO, and the NAOJ. The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated Universities, Inc. The authors thank the ALMA staff for their excellent support.

This work is supported by the projects PRIN-NAF 2019 "Planetary systems at young ages (PLATEA)" and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program, for the Project "The Dawn of Organic Chemistry" (DOC), grant agreement No. 741002; the PRIN-NAF The Cradle of Life— GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA); and the European Union's Horizon 2020 research and innovation programs under projects "Astro-Chemistry Origins" (ACO), grant No. 811312.

This project is also supported by a Grant-in-Aid from Japan Society for the Promotion of Science (KAKENHI: Nos. 18H05222, 19H05069, 19K14753, 20K14533, 20H05845).

D.J. is supported by NRC Canada and by an NSERC Discovery Grant. I.J.-S. has received partial support from the Spanish State Research Agency (AEI; project No. PID2019-105552RB-C41).

Data analysis was in part carried out on the common-use data analysis computer system at the Astronomy Data Center (ADC) of the National Astronomical Observatory of Japan.

## ORCID iDs

Satoshi Ohashi https://orcid.org/0000-0002-9661-7958

Claudio Codella https://orcid.org/0000-0003-1514-3074

Nami Sakai https://orcid.org/0000-0002-3297-4497

Claire J. Chandler https://orcid.org/0000-0002-7570-5596

Cecilia Ceccarelli https://orcid.org/0000-0001-9664-6292

Felipe Alves https://orcid.org/0000-0002-7945-064X

Tomoyuki Hanawa https://orcid.org/0000-0002-7538-581X

Aurora Durán https://orcid.org/0000-0002-7622-5762

Cécile Favre https://orcid.org/0000-0002-5789-6931

Ana López-Sepulcre https://orcid.org/0000-0002-6729-3640

Laurent Loinard https://orcid.org/0000-0002-5635-3345

Seyma Mercimek https://orcid.org/0000-0002-0742-7934

Linda Podio https://orcid.org/0000-0003-2733-5372

Yichen Zhang https://orcid.org/0000-0001-7511-0034

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

Nadia Balucani https://orcid.org/0000-0001-5121-5683

Eleonora Bianchi https://orcid.org/0000-0001-9249-7082

Mathilde Bouvier https://orcid.org/0000-0003-0167-0746

Gemma Busquet https://orcid.org/0000-0002-2189-6278

Paola Caselli https://orcid.org/0000-0003-1481-7911

Emmanuel Caux https://orcid.org/0000-0002-4463-6663

Steven Charnley https://orcid.org/0000-0001-6752-5109

Spandan Choudhury https://orcid.org/0000-0002-7497-2713

Nicolas Cuello https://orcid.org/0000-0003-3713-8073

Marta De Simone https://orcid.org/0000-0001-5659-0140

François Dulieu https://orcid.org/0000-0001-6981-0421

Siyi Feng https://orcid.org/0000-0002-4707-8409

Logan Francis https://orcid.org/0000-0001-8822-6327

Tetsuya Hama https://orcid.org/0000-0002-4991-4044

Eric Herbst https://orcid.org/0000-0002-4649-2536

Shingo Hirano https://orcid.org/0000-0002-4317-767X

Tomoya Hirota https://orcid.org/0000-0003-1659-095X

Muneaki Imai https://orcid.org/0000-0002-5342-6262

Andrea Isella https://orcid.org/0000-0001-8061-2207

Izaskun Jímenez-Serra https://orcid.org/0000-0003-4493-8714

Doug Johnstone https://orcid.org/0000-0002-6773-459X

Claudine Kahane https://orcid.org/0000-0003-1691-4686

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Bertrand Lefloch https://orcid.org/0000-0002-9397-3826

Luke T. Maud https://orcid.org/0000-0002-7675-3565

Maria Jose Maureira https://orcid.org/0000-0002-7026-8163

Francois Menard https://orcid.org/0000-0002-1637-7393

Anna Miotello https://orcid.org/0000-0002-7997-2528

George Moellenbrock https://orcid.org/0000-0002-3296-8134

Shoji Mori https://orcid.org/0000-0002-7002-939X

Riouhei Nakatani https://orcid.org/0000-0002-1803-0203

Hideko Nomura https://orcid.org/0000-0002-7058-7682

Yasuhiro Oba https://orcid.org/0000-0002-6852-3604

Yuki Okoda https://orcid.org/0000-0003-3655-5270

Yoko Oya https://orcid.org/0000-0002-0197-8751

Jaime Pineda https://orcid.org/0000-0002-3972-1978

Albert Rimola https://orcid.org/0000-0002-9637-4554

Takeshi Sakai https://orcid.org/0000-0003-4521-7492

Dominique Segura-Cox https://orcid.org/0000-0003-3172-6763

Brian Svoboda https://orcid.org/0000-0002-8502-6431

Vianney Taquet https://orcid.org/0000-0003-0407-7489

Leonardo Testi https://orcid.org/0000-0003-1859-3070

Charlotte Vastel https://orcid.org/0000-0001-8211-6469

Serena Viti https://orcid.org/0000-0001-8504-8844

Naoki Watanabe https://orcid.org/0000-0001-8408-2872

Yoshimasa Watanabe https://orcid.org/0000-0002-9668-3592

Ci Xue https://orcid.org/0000-0003-2760-2119

Satoshi Yamamoto https://orcid.org/0000-0002-9865-0970

## References

Alves, F. O., Caselli, P., Girart, J. M., et al. 2019, *Sci*, 366, 90

Alves, F. O., Girart, J. M., Caselli, P., et al. 2017, *A&A*, 603, L3

Alves, F. O., Girart, J. M., Padovani, M., et al. 2018, *A&A*, 616, A56

Anderson, J. M., Li, Z.-Y., Krasnopolsky, R., et al. 2003, *ApJL*, 590, L107

André, P., Martin-Pintado, J., Despois, D., et al. 1990, *A&A*, 236, 180

André, P., Ward-Thompson, D., & Barsony, M. 1993, *ApJ*, 406, 122

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

Bjerkeli, P., Ramsey, J. P., Harsono, D., et al. 2019, *A&A*, 631, A64

Bjerkeli, P., van der Wiel, M. H. D., Harsono, D., et al. 2016, *Natur*, 540, 406

Bodenheimer, P. 1995, *ARA&A*, 33, 199

Bontemps, S., & André, P. 1997, in IAU Symp. 182, Herbig-Haro Flows and the Birth of Stars, ed. B. Reipburth & C. Bertout (Dordrecht: Kluwer), 63

Boss, A. P., & Keiser, S. A. 2014, *ApJ*, 794, 44

Caratti o Garatti, A., Giannini, T., Nisini, B., et al. 2006, *A&A*, 449, 1077

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

Codella, C., Cabrit, S., Gueth, F., et al. 2014, *A&A*, 568, L5

Dent, W. R. F., Matthews, H. E., & Walther, D. M. 1995, *MNRAS*, 277, 193

Endres, C. P., Schlemmer, S., Schilke, P., et al. 2016, *JMoSp*, 327, 95

Fendt, C., & Zinnecker, H. 1998, *A&A*, 334, 750

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

Furuya, R. S., Kitamura, Y., Wootten, A., et al. 2003, *ApJS*, 144, 71

Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, *ApJ*, 856, 23

Galametz, M., Maury, A., Girart, J. M., et al. 2020, *A&A*, 644, A47

Goodwin, S. P., Whitworth, A. P., & Ward-Thompson, D. 2004, *A&A*, 423, 169

Hara, C., Kawabe, R., Nakamura, F., et al. 2021, *ApJ*, 912, 34

Harris, R. J., Cox, E. G., Looney, L. W., et al. 2018, *ApJ*, 861, 91

Hennebelle, P., & Ciardi, A. 2009, *A&A*, 506, L29

Hirano, S., & Machida, M. N. 2019, *MNRAS*, 485, 4667

Hirota, T., Machida, M. N., Matsushita, Y., et al. 2017, *NatAs*, 1, 0146

Hsieh, C.-H., Lai, S.-P., Cheong, P.-I., et al. 2020, *ApJ*, 894, 23

Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, *ApJS*, 213, 13

Hull, C. L. H., & Zhang, Q. 2019, *FrASS*, 6, 3
Ikeda, N., Sunada, K., & Kitamura, Y. 2007, ApJ, 665, 1194

Jiménez-Serra, I., Báez-Rubio, A., Martín-Pintado, J., et al. 2020, ApJL, 897, L33

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Kawabe, R., Hara, C., Nakamura, F., et al. 2018, ApJ, 866, 141

Konigl, A., & Pudritz, R. E. 2000, in Protostars and Planets IV, ed. V. Mannings et al. (Tuscon, AZ: Univ. Arizona Press), 759

Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2011, ApJ, 733, 54

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

Lee, C.-F., Ho, P. T. P., Li, Z.-Y., et al. 2017, NatAs, 1, 0152

Lee, C.-F., Li, Z.-Y., Codella, C., et al. 2018, ApJ, 856, 14

Lee, J.-E., Lee, S., Dunham, M. M., et al. 2017, NatAs, 1, 0172

Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2013, ApJ, 774, 82

Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477

Machida, M. N. 2014, ApJL, 796, L17

Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2007, ApJ, 670, 1198

Machida, M. N., Matsumoto, T., Hanawa, T., et al. 2006, ApJ, 645, 1227

Matsumoto, T., Machida, M. N., & Inutsuka, S.-i. 2017, ApJ, 839, 69

Matsumoto, T., & Tomisaka, K. 2004, ApJ, 616, 266

Maureira, M. J., Pineda, J. E., Segura-Cox, D. M., et al. 2020, ApJ, 897, 59

Maury, A., Ohashi, N., & André, P. 2012, A&A, 539, A130

McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, in Astronomical Data Analysis Software and Systems XVI, ASP Conf. Ser. 376, ed. R. A. Shaw (San Francisco, CA: ASP), 127

Murillo, N. M., Bruderer, S., van Dishoeck, E. F., et al. 2015, A&A, 579, A114

Murillo, N. M., Harsono, D., McClure, M., et al. 2018, A&A, 615, L14

Murillo, N. M., & Lai, S.-P. 2013, ApJL, 764, L15

Murillo, N. M., Lai, S.-P., Bruderer, S., et al. 2013, A&A, 560, A103

Narayanan, G., & Logan, D. W. 2006, ApJ, 647, 1170

Offner, S. S. R., Kratter, K. M., Matzner, C. D., et al. 2010, ApJ, 725, 1485

Ohashi, S., Kataoka, A., Nagai, H., et al. 2018, ApJ, 864, 81

Ohashi, S., Sanhueza, P., Chen, H.-R. V., et al. 2016, ApJ, 833, 209

Okoda, Y., Oya, Y., Francis, L., et al. 2021, ApJ, 910, 11

Onishi, T., Mizuno, A., Kawamura, A., et al. 2002, ApJ, 575, 950

Oya, Y., Sakai, N., López-Sepulcre, A., et al. 2016, ApJ, 824, 88

Oya, Y., Sakai, N., Watanabe, Y., et al. 2017, ApJ, 837, 174

Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870

Price, D. J., & Bate, M. R. 2007, MNRAS, 377, 77

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 859, 165

Sai, J., Ohashi, N., Saigo, K., et al. 2020, ApJ, 893, 51

Saiki, Y., & Machida, M. N. 2020, ApJL, 897, L22

Sakai, N., Hanawa, T., Zhang, Y., et al. 2019, Natur, 565, 206

Santangelo, G., Murillo, N. M., Nisini, B., et al. 2015, A&A, 581, A91

Seifried, D., Banerjee, R., Pudritz, R. E., et al. 2013, MNRAS, 432, 3320

Shu, F. H. 1977, ApJ, 214, 488

Shu, F. H., Najita, J. R., Shang, H., et al. 2000, in Protostars and Planets IV, ed. V. Mannings (Tuscon, AZ: Univ. Arizona Press), 789

Snell, R. L., Loren, R. B., & Plambeck, R. L. 1980, ApJL, 239, L17

Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6

Takaishi, D., Tsukamoto, Y., & Suto, Y. 2021, PASJ, 73, L25

Tanaka, K. E. I., Zhang, Y., Hirota, T., et al. 2020, ApJL, 900, L2

Tobin, J. J., Bourke, T. L., Mader, S., et al. 2019, ApJ, 870, 81

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2011, ApJ, 740, 45

Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Natur, 538, 483

Tsukamoto, Y., Iwasaki, K., Okuzumi, S., et al. 2015, ApJL, 810, L26

Tsukamoto, Y., & Machida, M. N. 2013, MNRAS, 428, 1321

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037

Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2013, ApJ, 772, 22

Yu, T., & Chernin, L. M. 1997, ApJL, 479, L63

Zapata, L. A., Loinard, L., Rodríguez, L. F., et al. 2013, ApJL, 764, L14

Zhang, Y., Arce, H. G., Mardones, D., et al. 2016, ApJ, 832, 158

Zhang, Y., Higuchi, A. E., Sakai, N., et al. 2018, ApJ, 864, 76

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2018, MNRAS, 473, 4868

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2021, MNRAS, 505, 5142