# The Physical Properties of the SVS 13 Protobinary System: Two Circumstellar
Disks and a Spiraling Circumbinary Disk in the Making

Ana K. Diaz-Rodriguez Instituto de Astrofísica de Andalucía, CSIC, Glorieta de
la Astronomía s/n, 18008 Granada, Spain  UK ALMA Regional Centre Node, Jodrell
Bank Centre for Astrophysics, Department of Physics and Astronomy,
The University of Manchester, Oxford Road, Manchester M13 9PL, UK  Guillem
Anglada Instituto de Astrofísica de Andalucía, CSIC, Glorieta de la Astronomía
s/n, 18008 Granada, Spain  Guillermo Blázquez-Calero Instituto de Astrofísica
de Andalucía, CSIC, Glorieta de la Astronomía s/n, 18008 Granada, Spain  Mayra
Osorio Instituto de Astrofísica de Andalucía, CSIC, Glorieta de la Astronomía
s/n, 18008 Granada, Spain  José F. Gómez Instituto de Astrofísica de
Andalucía, CSIC, Glorieta de la Astronomía s/n, 18008 Granada, Spain  Gary A.
Fuller UK ALMA Regional Centre Node, Jodrell Bank Centre for Astrophysics,
Department of Physics and Astronomy,
The University of Manchester, Oxford Road, Manchester M13 9PL, UK  Instituto
de Astrofísica de Andalucía, CSIC, Glorieta de la Astronomía s/n, 18008
Granada, Spain  I. Physikalisches Institut, University of Cologne, Zülpicher
Str. 77, 50937 Köln, Germany  Robert Estalella Departament de Física Quàntica
i Astrofísica, Institut de Ciències del Cosmos, Universitat de Barcelona,
IEEC-UB, Martí i Franquès, 1, 08028 Barcelona, Spain  José M. Torrelles
Institut de Ciències de l’Espai (ICE, CSIC), Can Magrans s/n, E-08193
Cerdanyola del Vallès, Catalonia, Spain  Institut d’Estudis Espacials de
Catalunya (IEEC), E-08034 Barcelona, Catalonia, Spain  Sylvie Cabrit Univ.
Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France  Observatoire de Paris, PSL
University, Sorbonne Université, CNRS, LERMA, F-75014, Paris, France  Luis F.
Rodríguez Instituto de Radioastronomía y Astrofísica, Universidad Nacional
Autónoma de México, P.O. Box 3-72, 58090, Morelia, Michoacán, México  Charlène
Lefèvre Institut de RadioAstronomie Millimétrique (IRAM), 300 rue de la
Piscine, 38406 Saint Martin d’Hères, France  Enrique Macías Joint ALMA
Observatory, Alonso de Córdova 3107, Vitacura, Santiago 763-0355, Chile
European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura,
Santiago 763-0355, Chile  Carlos Carrasco-González Instituto de
Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, P.O.
Box 3-72, 58090, Morelia, Michoacán, México  Luis A. Zapata Instituto de
Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, P.O.
Box 3-72, 58090, Morelia, Michoacán, México  Itziar de Gregorio-Monsalvo
European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura,
Santiago 763-0355, Chile  Paul T. P. Ho Academia Sinica Institute of Astronomy
and Astrophysics, P.O. Box 23-141, Taipei 106, Taiwan  East Asian Observatory,
Hilo 96720, HI, USA

###### Abstract

We present VLA and ALMA observations of the close ($0\farcs 3$ = 90 au
separation) protobinary system SVS 13. We detect two small circumstellar disks
(radii $\sim$12 and $\sim$9 au in dust, and $\sim$30 au in gas) with masses of
$\sim$0.004-0.009 $M_{\sun}$ for VLA 4A (the western component) and
$\sim$0.009-0.030 $M_{\sun}$ for VLA 4B (the eastern component). A
circumbinary disk with prominent spiral arms extending $\sim$500 au and a mass
of $\sim$0.052 $M_{\sun}$ appears to be in the earliest stages of formation.
The dust emission is more compact and with a very high optical depth toward
VLA 4B, while toward VLA 4A the dust column density is lower, allowing the
detection of stronger molecular transitions. We infer rotational temperatures
of $\sim$140 K, on scales of $\sim$30 au, across the whole source, and a rich
chemistry. Molecular transitions typical of hot corinos are detected toward
both protostars, being stronger toward VLA 4A, with several ethylene glycol
transitions detected only toward this source. There are clear velocity
gradients, that we interpret in terms of infall plus rotation of the
circumbinary disk, and purely rotation of the circumstellar disk of VLA 4A. We
measured orbital proper motions and determined a total stellar mass of 1
$M_{\sun}$. From the molecular kinematics we infer the geometry and
orientation of the system, and stellar masses of $\sim$0.26 $M_{\sun}$ for VLA
4A and $\sim$0.60 $M_{\sun}$ for VLA 4B.

Circumstellar disks (235) — Close binary stars (254) — Exoplanet formation
(492) — Molecular spectroscopy (2095) — Multiple stars (1081) — Planetary
system formation (1257) — Protoplanetary disks (1300) — Protostars (1302) —
Star formation (1569)

††facilities: ALMA, VLA††software: Astropy (http://www.astropy.org; Astropy
Collaboration et al. 2013, 2018), CASA (http://casa.nrao.edu; McMullin et al.
2007), CASSIS (http://cassis.irap.omp.eu; Vastel et al. 2015)

##  1 Introduction

Although most stars belong to binary systems (Leinert et al., 1993; Raghavan
et al., 2010; Duchêne & Kraus, 2013), the way in which such systems form still
remains poorly known, especially in the first stages of the process. Two main
different mechanisms have been suggested for binary formation (e.g., Bate
2015), namely disk fragmentation (e.g., Bonnell 1994; Bonnell & Bate 1994a, b;
Kratter et al. 2010) and the large-scale fragmentation of cores and filaments,
either without turbulence (Matsumoto & Hanawa, 2003) or through turbulent
fragmentation (e.g., Padoan & Nordlund 1999; Offner et al. 2010). Although
several observational studies have investigated which mechanism is dominant
(e.g., Tobin et al. 2018), the debate is not yet settled. On the other hand,
dissipative dynamical interactions within a cluster are also a major factor to
be considered, especially for understanding the formation of close binary
systems (Bate, 2012, 2015).

At the early stages, a binary protostellar system will usually be embedded in
a more massive infalling envelope from which it accretes to its final mass.
Accretion of this material changes the masses of the components of the binary,
the binary’s orbit, and it is also largely responsible for forming the final
configuration of disks in the system (Bate, 2015). Thus, besides the
investigation of the fragmentation mechanism, a fundamental issue is to study
the early evolution of the binary during the formation process, as this
process determines not only the final stellar masses, but also the formation
of disks and, consequently, the formation of planets (e.g., Bate 2015, 2018).

In a binary system, three disks are possible: a circumprimary disk, a
circumsecondary disk, and a circumbinary disk. According to numerical
simulations, the specific angular momentum of the accreted gas relative to
that of the binary is the most important parameter to determine the evolution
of the binary properties and, in particular, the formation of a circumbinary
disk (e.g., Bate 2015).

Since circumbinary disks are not geometrically thin, high resolution 3D
simulations are required to study these structures in young binary systems.
The recent 3D simulations of Matsumoto et al. (2019) provide fine details of
the formation and evolution of circumstellar disks, spiral arms, and
circumbinary structures in protobinary systems. These simulations show the
formation of spirals, a bridge structure between the circumstellar disks, and
the development of a slowly rotating asymmetric pattern in the circumbinary
disk. Some of these features have been observationally detected (e.g., in the
systems L1551 NE; Takakuwa et al. 2014, 2017, and L1551 IRS 5; Cruz-Sáenz de
Miera et al. 2019; Takakuwa et al. 2020), but this kind of work is still
scarce and more observational studies providing good estimates of physical
parameters suitable for making detailed comparisons between observations and
3D simulations are needed.

Therefore, it is imperative to conduct in-depth observational studies of the
earliest stages of the binary star formation process, i.e., to study
protobinary systems. Binaries with separations $\sim$1000 au are called wide
binaries and are attributed to core fragmentation, while those with
separations of $\sim$100 au are in the boundary between wide and close
binaries. In the latter, the interaction between the two stars significantly
affects their evolution. This is particularly relevant when the two stars have
different initial properties. Thus, protobinaries with separations of the
order of 100 au or less, whose observation is currently achievable with the
resolution of the large radio interferometers, are excellent targets to study
their properties and to compare them with those predicted by numerical
simulations (e.g., Bate 2018; Matsumoto et al. 2019).

So far, only for a limited number of these systems there is available (but
partial) observational information about their physical properties, such as
temperatures, kinematics (from radial velocities and proper motions), presence
of circumbinary disks, stellar and disk masses (from dust and molecular
emission and from dynamics), chemical properties, interaction with the
surrounding medium (e.g., outflows and interaction with the ambient molecular
cloud, flyby encounters and tidal interaction with nearby young stars). Some
examples where these studies have been performed are the binary systems L1551
IRS 5 (Rodríguez et al., 1998; Cruz-Sáenz de Miera et al., 2019; Takakuwa et
al., 2020), RW Aur (Cabrit et al., 2006; Rodriguez et al., 2018), L1551 NE
(Takakuwa et al., 2014, 2017), XZ Tau (Carrasco-González et al., 2009; Zapata
et al., 2015; Osorio et al., 2016; Ichikawa et al., 2021), GG Tau (Phuong et
al., 2020), BHB2007 (Alves et al., 2019), UX Tau (Zapata et al., 2020), and
IRAS 16293$-$2422 (Maureira et al., 2020).

A particularly remarkable case is SVS 13, the proposed exciting source of the
HH 7-11 system (Strom et al., 1976). SVS 13 is a nearby ($d$ = 300 pc; Ortiz-
León et al. 2018; Gaia Collaboration et al. 2021, see also Hirota et al. 2008)
optically visible star (e.g., Eisloeffel et al. 1991) but also a strong mm
source (e.g., Looney et al. 2000). VLA 3.6 cm observations (Anglada et al.,
2000) revealed that SVS 13 is a close binary system, with the two components
(VLA 4A and VLA 4B) separated in projection by $0\farcs 3$ (90 au). Further
observations (Anglada et al., 2004) revealed that only one of the components,
VLA 4B, presents a strong 7 mm emission, that was attributed to dust. These
results were interpreted as implying that the development of a protoplanetary
disk occurs preferentially in one of the components, with the disk being much
weaker in the other component. Strong maser emission has been detected in the
proximity of SVS 13 and toward other nearby sources in the field (Rodríguez et
al., 2002). More recently, the close environment of the SVS 13 system has been
explored in the near-IR using laser adaptive-optics at the Keck telescope
(Hodapp & Chini, 2014), revealing a $0\farcs 2$ long microjet traced by [Fe
II] emission and a sequence of H2 bubbles extending a few arcsec away from the
position of VLA 4B. The system has also been observed at mm wavelengths with
the 30 m and PdBI of IRAM (e.g., Lefèvre et al. 2017; Bianchi et al. 2019;
Belloche et al. 2020b) tracing the dust and several molecular transitions with
angular resolutions ranging from $0\farcs 3$ to $\sim$30$\arcsec$. The source
has also been observed as part of the VANDAM survey at cm wavelengths with the
VLA (Tychoniec et al., 2018) and at 1.3 mm continuum and several molecular
line transitions with a resolution of $\sim 0\farcs 2$-$0\farcs 3$ with ALMA
(Tobin et al., 2018).

In this paper, we present new VLA results in the continuum from 3 cm to 7 mm,
as well as ALMA 0.9 mm continuum and line results with improved angular
resolution and sensitivity. These results, combined with recent Gaia
astrometry, shed new light on this protobinary system probing scales as small
as $\sim 0\farcs 05$ ($\sim$15 au).

Our paper is focused on characterizing as well as possible the physical and
chemical properties of the SVS 13 protobinary system. We defer to a subsequent
paper a detailed comparison with both model simulations and observations from
the literature of other sources reaching a similar degree of detail. The rest
of this paper is organized as follows. In Section 2 we describe the
observations and data reduction procedures. In Section 3 we report the sources
detected and discuss in detail our main results on the SVS 13 protobinary: In
Section 3.1 we describe our main VLA results on the ionized emission and
proper motions of the protostars; in Section 3.2 we describe the circumstellar
and circumbinary dust observed with ALMA; in Section 3.3 we analyze the
overall spectral energy distribution; in Section 3.4 we present our main ALMA
results regarding the distribution and kinematics of the molecular gas, as
well as the molecular transitions identified; in Section 3.5 we estimate the
physical parameters of the system (temperature, disk masses, and stellar
masses), and we discuss about their implications; in Section 3.6 we discuss
the chemical properties and in particular the hot corino nature of the
protostellar sources. Finally, Section 4 summarizes our main results and
conclusions.

##  2 Observations and Data Reduction

We present continuum observations in the 0.7-3 cm (10-43 GHz) range carried
out with the Karl G. Jansky Very Large Array (VLA) of the National Radio
Astronomy Observatory (NRAO)111The National Radio Astronomy Observatory is a
facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc., as well as continuum and line
observations at $\sim$0.9 mm ($\sim$333 GHz) taken with the Atacama Large
Millimeter/submillimeter Array (ALMA). Except when noted otherwise, the images
presented in this work have been corrected by the primary beam response.

###  2.1 VLA Continuum and H2O Maser Observations

The VLA observations were carried out as part of our Project 15A-260 (PI: Ana
K. Diaz-Rodriguez) in the B configuration at Ka (9 mm) and Q (7 mm) bands on
2015 February 14 and March 7, respectively, and in the A configuration at X (3
cm), K (1.3 cm), Ka, and Q bands on 2015 June 28, June 27, July 7, and June
20, respectively. Data at X band had poor quality, and we repeated the
observations at this band in the A configuration on 2016 November 22 as part
of DDT 16B-417 project (PI: Ana K. Diaz-Rodriguez). The results presented in
this paper at X band are obtained using only the 2016 data.

The bandwidth was 4 GHz at X band, and 8 GHz at K, Ka, and Q bands. We used
wideband mode setup for X, Ka, and Q bands, and spectral line mode setup for
the K band to observe the H2O maser emission at 22.23508 GHz. For this, we
included a narrow subband of 4 MHz covering the maser line, with 128 channels
of 31.25 kHz (0.42 km s-1) each; the other 63 subbands had 128 MHz bandwidth
each, and were used for the continuum imaging. The phase center was at
coordinates $\alpha$(J2000) = 03h29m3$\fs$75,
$\delta$(J2000)=+31°16′3$\farcs$9 in all VLA observations. Absolute flux
calibrators were J0137+3309 (3C48) at X, K, and Ka bands, and J0542+4951
(3C147) at Q band. Bandpass calibration was performed using J0137+3309 (3C48)
at X band, and J0319+4130 (3C84) for the rest of the bands. J0336+3218 was
used for complex gain calibration at all bands. The time on source in A
configuration was 16.5 min at Q band, 23.6 min at Ka band, 39.7 min at K band,
and 39.2 min at X band, while in B configuration was 20.4 min at Q band and
21.7 min at Ka band.

The data were first pipeline calibrated and then manually edited and
calibrated again using the Common Astronomy Software Applications package
(CASA) version 4.3.1. Self-calibration was attempted but no improvement of the
signal-to-noise ratio was obtained. Images were made with the task CLEAN of
CASA version 4.6.0. The synthesized beams of the continuum images range from
$\sim 0\farcs 20$ (X band) to $\sim 0\farcs 06$ (Q band) and the rms noise
from 6 $\mu$Jy beam-1 to 23 $\mu$Jy beam-1, respectively (see Fig. 1). In the
K-band cube we detect strong H2O maser emission at a few Jy level in the
proximity of SVS 13, indicating that the previously reported masers (Rodríguez
et al., 2002) are still active; these results will be presented in a
forthcoming paper.

Figure 1: VLA images of SVS 13 at X (3 cm), K (1.3 cm), Ka (9 mm), and Q (7
mm) bands obtained with the A configuration. The plus sign marks the optical
position of the star SVS 13 from Gaia EDR3 (Gaia Collaboration et al., 2021)
extrapolated to the epoch of the VLA observations (a correction $<0\farcs
014$; see text); after this correction the Gaia position is accurate to within
a fraction of mas, well below the size of the plus sign. Contour levels are
$-$3, 3, 5, 7, 10, 14, 19, 25, 35, and 50 times the rms of each image. The
synthesized beams and the rms are shown in the bottom left and right corners,
respectively, of the images. We used natural weighting for images at K, Ka,
and Q bands, and Briggs weighting with the robust parameter set to 0 for the
image at X band. The synthesized beams are $0\farcs 22\times 0\farcs 16$, PA =
$-80.3\arcdeg$ at X band; $0\farcs 11\times 0\farcs 10$, PA = $-80.7\arcdeg$
at K band; $0\farcs 10\times 0\farcs 07$, PA = $-72.9\arcdeg$ at Ka band; and
$0\farcs 06\times 0\farcs 05$, PA = $-83.4\arcdeg$ at Q band.

###  2.2 ALMA Continuum and Line Observations

The ALMA observations were carried out in Band 7 (0.9 mm) as part of our
projects 2015.1.01229.S (Cycle 3; PI: Guillem Anglada) and 2016.1.01305.S
(Cycle 4; PI: Guillem Anglada), with angular resolutions of $\sim 0\farcs 1$
and $\sim 0\farcs 3$, respectively, and a maximum recoverable scale of
$2\farcs 6$ in both projects. Data from Cycle 3 were obtained during two
identical runs on 2016 September 9 and 10, with 36 antennas operating with
baselines ranging from 15 m to 3247 m. The total time spent on source was 1.6
h. The precipitable water vapor was 0.44 mm and 0.41 mm, during the first and
second day of observations, respectively. Data from Cycle 4 were obtained on
2016 November 24, with 40 antennas operating with baselines ranging from 15 m
to 704 m. The time on source was 4.15 min and the precipitable water vapor was
0.63 mm. In all cases, the calibrators J0237+2848, J0238+1636, and J0336+3218
were used for bandpass, absolute flux, and complex gain calibration,
respectively.

In the Cycle 3 observations, the correlator was configured with one baseband
set to low spectral resolution continuum mode, and three basebands set to high
spectral resolution line mode. Each line baseband was subdivided into two
equal spectral windows of 0.234 GHz bandwidth. The continuum baseband was
centered at 333 GHz and had 1.875 GHz of bandwidth. The line spectral windows
were centered on the sky frequencies of the 13CO (3-2), SO2 (11-12), CS (7-6),
HC15N (4-3), H13CN (4-3), and 12CO (3-2), whose rest frequencies are
330.587965 GHz, 331.580244 GHz, 342.882857 GHz, 344.200109 GHz, 345.339756
GHz, and 345.795990 GHz, respectively. The channel width was 0.1058 km s-1
(0.122 MHz) for the CO spectral windows, 0.212 km s-1 (0.244 MHz) for H13CN,
CS, and SO2, and 0.426 km s-1 (0.488 MHz) for the HC15N spectral window. The
line free-channels of the line spectral windows (totaling $\sim$ 0.3 GHz) and
of the continuum baseband were used to form the continuum image, covering a
total bandwidth of $\sim$2.175 GHz.

In the Cycle 4 observations, the correlator was configured with three
basebands set to low spectral resolution continuum mode and one set to high
spectral resolution line mode. Continuum basebands had 1.875 GHz of bandwidth
each and were centered at 331.5, 333.5, and 344.0 GHz. The high spectral
resolution baseband was centered on the sky frequency of 12CO (3-2) (rest
frequency 345.795990 GHz), with a channel width of 0.426 km s-1 (0.488 MHz)
and 0.9375 GHz bandwidth.

After the standard pipeline calibration, the continuum data from Cycle 3 were
phase-only self-calibrated to increase the signal to noise ratio. We used two
rounds in which the solutions were found in intervals that spanned the entire
scan length of 6.07 minutes (first round), and 6.05 s, which is the length of
a single integration (second round). The self-calibration solutions were
applied to the spectral line data and the continuum emission was subtracted.
Briggs-weighted (robust = 0.5) images obtained in this way have a synthesized
beam of $0\farcs 170\times 0\farcs 084$, PA = $-2.3\arcdeg$, and an rms of
0.075 mJy beam-1 for the continuum, and 2-3.5 mJy beam-1 per channel
(depending on the spectral window) for the line images (cubes). Self-
calibrated Cycle 3 data were concatenated with standard-calibrated Cycle 4
data, and we followed the same steps to self-calibrate the whole measurement
continuum data set. The Briggs-weighted (robust = 0.5) image obtained with the
concatenated data has a synthesized beam of $0\farcs 176\times 0\farcs 088$,
PA = $-2.0\arcdeg$ and an rms of 0.10 mJy beam-1. In our discussion, we use
the image obtained with only the Cycle 3 data which has slightly better
angular resolution and rms.

The ALMA observations were pipeline calibrated with CASA versions 4.7.01
(Cycle 3) and 4.7.2 (Cycle 4). For both self-calibration and imaging, CASA
version 4.5.2 was used. The images were made using the task CLEAN of CASA.

##  3 Results and Discussion

###  3.1 VLA Results

####  3.1.1 Radio Continuum Sources

Figure 1 shows VLA images of the continuum emission toward SVS 13 at 3, 1.3,
0.9, and 0.7 cm. The two compact sources previously identified as VLA 4A (the
western one) and VLA 4B (the eastern one) by Anglada et al. (2000, 2004) and
by Tobin et al. (2016b) are detected at all the observed wavelengths. Table
3.1.1 lists the positions of the two sources and Table 3.1.1 their flux
densities obtained from our images and from the literature.

Additional sources detected in the observed field are VLA 2, VLA 3, VLA 17,
VLA 20, and VLA 22. All these sources were previously reported at cm
wavelengths by Rodríguez et al. (1999). Our VLA images are shown in Figures 17
to 21, and their positions and flux densities are given in Table 3. Source VLA
22 splits into two components, VLA 22N and VLA 22S, separated 40 au in
projection (Fig. 21), and both present a high degree of circular polarization.
The high degree of circular polarization of VLA 22 was already noted by
Rodríguez et al. (1999), but the binary was not resolved in those
observations. Interestingly, Rodríguez et al. (1999) reported a left-handed
circularly polarized source, while in our observations right-handed
polarization is dominant. In the following sections we discuss in more detail
the properties of the sources VLA 4A and VLA 4B, which are associated with SVS
13. Further study of the remaining sources detected in the field is deferred
to a forthcoming paper (A. K. Diaz-Rodriguez, in prep.).

Table 1: Positions of the VLA 4A and VLA 4B Radio SourcesaaPositions derived
from elliptical Gaussian fits. Uncertainties correspond to those of absolute
positions, and are calculated by adding in quadrature a systematic error of
0$\farcs$01 (for VLA data, Dzib et al. 2017) or 0$\farcs$005 (for ALMA data,
Remijan et al. 2020) to the formal error of the fit, $0.5\,\theta/{\rm SNR}$
(Reid et al., 1988). The systematic error accounts for uncertainties
introduced by the phase calibration process, and the error of the fit for
those due to the convolved source size ($\theta$) and signal-to-noise ratio
(SNR).

Epoch | Wavelength | VLA 4A | VLA 4B |
---|---|---|---|---
(yr) | (mm) | RA (J2000) | DEC (J2000) | RA (J2000) | DEC (J2000) | Ref.
1989.04 | 36 | 03 29 03.730$\pm$0.002 | 31 16 04.02$\pm$0.02 | 03 29 03.751$\pm$0.002 | 31 16 04.08$\pm$0.02 | 1
1994.31 | 36 | 03 29 03.730$\pm$0.002 | 31 16 03.95$\pm$0.02 | 03 29 03.758$\pm$0.004 | 31 16 04.01$\pm$0.02 | 1
1996.98 | 36 | 03 29 03.731$\pm$0.002 | 31 16 03.94$\pm$0.02 | 03 29 03.757$\pm$0.002 | 31 16 04.01$\pm$0.02 | 1
1998.32bbFrom an image obtained concatenating data from epochs J1998.24 and J1998.40 reported in Anglada et al. (2000) and Carrasco-González et al. (2008). | 36 | 03 29 03.7376$\pm$0.0017 | 31 16 03.957$\pm$0.026 | 03 29 03.7575$\pm$0.0015 | 31 16 03.968$\pm$0.023 | 1, 6
1999.50 | 36 | 03 29 03.734$\pm$0.002 | 31 16 03.96$\pm$0.02 | 03 29 03.757$\pm$0.001 | 31 16 03.91$\pm$0.02 | 1
2001.34 | 7 | 03 29 03.7405$\pm$0.0020 | 31 16 03.946$\pm$0.030 | 03 29 03.7594$\pm$0.0008 | 31 16 03.937$\pm$0.012 | 2ccAfter recalibration of the originally published data.
2002.29 | 13 | 03 29 03.7357$\pm$0.0015 | 31 16 03.935$\pm$0.022 | 03 29 03.7596$\pm$0.0011 | 31 16 03.914$\pm$0.016 | 3
2003.40 | 13 | 03 29 03.7386$\pm$0.0011 | 31 16 03.876$\pm$0.017 | 03 29 03.7609$\pm$0.0007 | 31 16 03.908$\pm$0.011 | 4
2014.56 | 9 | 03 29 03.7430$\pm$0.0007 | 31 16 03.790$\pm$0.010 | 03 29 03.7660$\pm$0.0007 | 31 16 03.810$\pm$0.010 | 5
2015.47 | 7 | 03 29 03.7448$\pm$0.0007 | 31 16 03.780$\pm$0.010 | 03 29 03.7676$\pm$0.0007 | 31 16 03.788$\pm$0.010 | 6
2015.49 | 13 | 03 29 03.7442$\pm$0.0007 | 31 16 03.785$\pm$0.010 | 03 29 03.7674$\pm$0.0007 | 31 16 03.788$\pm$0.010 | 6
2015.52 | 9 | 03 29 03.7455$\pm$0.0008 | 31 16 03.794$\pm$0.012 | 03 29 03.7676$\pm$0.0007 | 31 16 03.789$\pm$0.011 | 6
2016.69 | 0.9 | 03 29 03.7460$\pm$0.0003 | 31 16 03.769$\pm$0.005 | 03 29 03.7691$\pm$0.0003 | 31 16 03.761$\pm$0.005 | 6ddFrom the ALMA image obtained using baselines $>$750 k$\lambda$, where the extended emission has been filtered out (Figure 8).
2016.89 | 30 | 03 29 03.7457$\pm$0.0010 | 31 16 03.771$\pm$0.015 | 03 29 03.7689$\pm$0.0008 | 31 16 03.775$\pm$0.012 | 6

References. — (1) Carrasco-González et al. (2008); (2) Anglada et al. (2004);
(3) G. Anglada, private communication; (4) C. Carrasco-González, private
communication; (5) Tobin et al. (2016b); (6) This work.

Table 2: Flux Densities of the VLA 4A and VLA 4B Radio Sources

Frequency | Wavelength | Flux density (mJy)aaUncertainties account for the rms noise of the images and 10% of absolute flux calibration uncertainty, except for the data reported by (1), the 231 GHz data point reported by (7), and the data reported by (8), whose calibration uncertainties are 5%, 15%, and 20%, respectively. | EpochbbWhen a range of epochs is indicated, flux densities have been measured in an image made by concatenating data obtained at different epochs. | Resolution |
---|---|---|---|---|---
(GHz) | (mm) | VLA 4A | VLA 4B | (yr) | (″) | Ref.
4.68 | 64 | 0.065$\pm$0.008 | 0.059$\pm$0.006 | 2014.2-2014.3 |  $\sim$0.4 | 1
7.31 | 41 | 0.109$\pm$0.008 | 0.093$\pm$0.007 | 2014.2-2014.3 |  $\sim$0.4 | 1
8.33 | 36 | 0.17$\pm$0.03 | 0.12$\pm$0.03 | 1994.3-1999.5 |  $\sim$0.2 | 2
10.0 | 30 | 0.045$\pm$0.009 | 0.084$\pm$0.013 | 2016.9 |  $\sim$0.2 | 3
22.1 | 13 | 0.24$\pm$0.03 | 0.61$\pm$0.06 | 2015.5 |  $\sim$0.1 | 3
30.0 | 10 | 0.38$\pm$0.05 | 1.12$\pm$0.11 | 2013.7-2015.9 |  $\sim$0.1 | 4
33.0 | 9 | 0.46$\pm$0.05 | 1.36$\pm$0.14 | 2015.1-2015.5 |  $\sim$0.2 | 3
33.0 | 9 | 0.41$\pm$0.06 | 1.22$\pm$0.13 | 2015.5 |  $\sim$0.1 | 3
37.5 | 8 | 0.51$\pm$0.06 | 1.76$\pm$0.18 | 2013.7-2015.9 |  $\sim$0.1 | 4
42.8 | 7 | 0.34$\pm$0.12 | 1.29$\pm$0.19 | 2001.3 |  $\sim$0.1 | 5ccAfter recalibration of the originally published data.
44.0 | 7 | 0.61$\pm$0.07 | 1.96$\pm$0.20 | 2015.2-2015.5 |  $\sim$0.1 | 3
44.0 | 7 | 0.37$\pm$0.06 | 2.08$\pm$0.22 | 2015.5 |  $\sim$0.05 | 3
86.8 | 3.5 | 26$\pm$3ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 1998.2 |  $\sim$3.5 | 6
94 | 3.2 | 21$\pm$2ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 2011.9 |  $\sim$1.2 | 7
107.7 | 2.7 | 38$\pm$10ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 1996.3-1998.2 |  $\sim$0.6 | 8
230.5 | 1.3 | 320$\pm$34ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 1998.2 |  $\sim$1.5 | 6
231 | 1.3 | 120$\pm$14ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 2011.1 |  $\sim$0.4 | 7
232.5 | 1.3 | 400$\pm$44ddTotal observed flux density (the binary is unresolved). Assumed to be dominated by the emission of the circumbinary disk. | 2015.7 |  $\sim$0.2 | 9
338.2 | 0.9 | 83$\pm$9eeMeasured in the image made only with the Cycle 3 ALMA data and baselines $>$750 k$\lambda$, where the extended emission has been filtered out (Fig. 8). Assumed to be dominated by the emission of the circumstellar disk. | 110$\pm$11eeMeasured in the image made only with the Cycle 3 ALMA data and baselines $>$750 k$\lambda$, where the extended emission has been filtered out (Fig. 8). Assumed to be dominated by the emission of the circumstellar disk. | 2016.7 |  $\sim$0.08 | 3
338.2 | 0.9 | 1102$\pm$110ffMeasured in an image made using both the Cycle 3 and 4 ALMA data, to better recover the extended emission. Assumed to be dominated by the emission of the circumbinary disk. | 2016.7-2016.9 |  $\sim$0.1 | 3

References. — (1) Tychoniec et al. (2018); (2) Anglada et al. (2000), (3) This
work, (4) Tobin et al. (2016b), (5) Anglada et al. (2004), (6) Bachiller et
al. (1998), (7) Maury et al. (2019), (8) Looney et al. (2000), (9) Tobin et
al. (2018)

Table 3: Positions and Flux Densities of Other Radio Sources in the SVS 13 FieldaaPositions derived from elliptical Gaussian fits to the ALMA image at 0.9 mm, except for sources VLA 2 and VLA 22, for which the VLA image at K band (1.3 cm) has been used. Positional uncertainties are calculated as in Table 3.1.1. VLA flux densities are measured in A-configuration (3 and 1.3 cm) and combined A+B-configuration (9 and 7 mm) total intensity images, except for VLA 22 where polarized A-configuration (3 cm and 7 mm) and B-configuration (9 mm) images are used. Flux density uncertainties include a 10% absolute flux calibration uncertainty. | Position | Flux Density (mJy)
---|---|---
Source | RA (J2000) | DEC (J2000) | 3 cm | 1.3 cm | 9 mm | 7 mm | 1.3 mmbbFrom the ALMA observations by Tobin et al. (2018). | 0.9 mm
VLA 2 | 03 29 01.9719$\pm$0.0007 | 31 15 38.044$\pm$0.010 | 0.64$\pm$0.07 | 1.59$\pm$0.17 | 1.76$\pm$0.18 | 2.11$\pm$0.25 |  |  $\cdots$ ccSource outside the imaged field of view.
VLA 17 | 03 29 03.0813$\pm$0.0007 | 31 15 51.701$\pm$0.010 | 0.023$\pm$0.006 | 0.31$\pm$0.04 | 1.18$\pm$0.12 | 2.12$\pm$0.22 | 222.2$\pm$22.3 | 400$\pm$40
VLA 3 | 03 29 03.3888$\pm$0.0007 | 31 16 01.571$\pm$0.010 | 0.067$\pm$0.013 | 0.183$\pm$0.024 | 0.34$\pm$0.04 | 0.28$\pm$0.04 | 13.9$\pm$1.7 | 27$\pm$3
VLA 20 | 03 29 04.2545$\pm$0.0007 | 31 16 09.067$\pm$0.010 | 0.087$\pm$0.016 | 0.127$\pm$0.018 | 0.16$\pm$0.03 | 0.13$\pm$0.04 | 4.4$\pm$1.1 | 8.4$\pm$0.9
VLA 22N | 03 29 05.7621$\pm$0.0008 | 31 16 39.617$\pm$0.012 | 0.009$\pm$0.007ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. | 0.090$\pm$0.017 | 0.25$\pm$0.05ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. | 0.39$\pm$0.15ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. |  |  $\cdots$ ccSource outside the imaged field of view.
VLA 22S | 03 29 05.7649$\pm$0.0008 | 31 16 39.503$\pm$0.012 | 0.029$\pm$0.007ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. | 0.099$\pm$0.018 | 0.05$\pm$0.05ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. | 0.23$\pm$0.16ddThe two sources appear blended in the total intensity ($I$) image at this wavelength. Flux densities are obtained from two-Gaussian fits to the LCP and RCP images shown in Figure 21, assuming unresolved sources at the positions (taken as fixed parameters) of VLA 22N and VLA 22S. The total intensity is calculated as $I=(LCP+RCP)/2$. |  |  $\cdots$ ccSource outside the imaged field of view.

####  3.1.2 Binary Protostars Traced by the Ionized Gas

As shown in Figure 1, the two sources detected toward SVS 13 are separated by
$\sim 0\farcs 30$, corresponding to 90 au in projection at a distance of 300
pc. The sources appear compact, although there is a hint of weak emission,
which is more evident at X and Ka bands, connecting the two sources. As noted
previously by Anglada et al. (2004), source VLA 4B is significantly stronger
than VLA 4A at mm wavelengths (intensity ratios of 3.8 at 7 mm and 3.7 at 9
mm), while at longer wavelengths the difference in the observed intensities of
the two sources is less significant (intensity ratios of 2.4 at 1.3 cm and 1.6
at 3 cm, see further discussion in Section 3.3).

In the figures, we plot the optical position of the star SVS 13 as obtained
from Gaia EDR3 Catalog (Gaia Collaboration et al. 2021; given at epoch
J2016.0), corrected by the reported Gaia proper motions (see Table 4) to the
epoch of the VLA observations (a small correction of not more than $0\farcs
014$) and converted to equinox J2000 coordinates with the task ICRSToJ2000 of
the CASA Extensions AnalysisUtils Python package. As can be seen in the
figure, the optical position falls closer to VLA 4B than to VLA 4A. This
association is in disagreement with the initial identification of the optical
star with VLA 4A made by Anglada et al. (2000) from the limited astrometric
precision ($\sim$0$\farcs$3) optical data available at that time. The
identification of VLA 4B with the optically visible star was already proposed
by Hodapp & Chini (2014), who noted the positional agreement of the near-IR
position measured by 2MASS (with an astrometric accuracy of $0\farcs 1$,
better than that of the previous optical data) and the radio source VLA 4B,
assuming that the 2MASS source and the optical star were the same object. The
identification was corroborated by Lefèvre et al. (2017) by comparing more
precise absolute astrometry ($\sim$0$\farcs$02) PdBI mm data of the base of
the SiO/CO jet, which was found to be associated with the near-IR jet imaged
by Hodapp & Chini (2014).

Source VLA 4A is closer to the densest part of the ambient cloud, as traced by
the moderately high angular resolution ($\sim 4^{\prime\prime}$) VLA ammonia
images (Rudolph et al., 2001; Diaz-Rodriguez, 2021). These images show that
the ammonia emission extends preferentially toward the west and south-west,
suggesting that the VLA 4A component of the binary could be more extincted by
foreground material than VLA 4B, consistent with the association of VLA 4B
with the optical star.

Nevertheless, we note that the coincidence between the optical and the radio
position is not perfect (see Fig. 1 and Table 3.1.1). The discrepancy is $\sim
0\farcs 07$, larger than the estimated astrometric uncertainty, which is $\sim
0\farcs 01$ (see discussion in Section 3.1.3).

The VLA Ka-band image in Figure 1 also shows a hint of a SE elongation of the
source VLA 4B, which is difficult to confirm since the beam is elongated in
the same direction. In Figure 2 we show the Ka-band image obtained from the
concatenated A+B configuration data. Although with a lower angular resolution,
this image shows quite clearly an extension of the emission from the position
of VLA 4B toward the SE. Interestingly, this SE protuberance coincides very
well with the axis, at PA = 145∘, of the $0\farcs 2$ long [Fe II] jet imaged
in the near-IR by Hodapp & Chini (2014); this microjet pierces inside the
first of a series of bubbles and fragments of bubbles traced in the H2 1–0
S(1) line to the SE of SVS 13. A false-color image of the microjet and its
associated bubble is overlaid on the contour image of the 9 mm emission shown
in Figure 2. The IR and radio images were aligned by matching the IR emission
of the SVS 13 star (in light blue in the image) with the Gaia EDR3 optical
position of SVS 13 corrected to the epoch of our VLA observations (plus sign).
We note that Hodapp & Chini (2014) found evidence of expanding motion of the
bubble, with its front end moving at a projected velocity of 20 mas yr-1 (28.5
km s-1); since their image was obtained $\sim$ 1.5 yr earlier than our VLA
image, the actual size of the bubble at the epoch of our observations could be
$\sim 0\farcs 05$ larger than what is shown in the figure, which does not
affect the accuracy of our comparison.

Figure 2: VLA image of SVS 13 at Ka band (9 mm) obtained from the concatenated
data of the A and B configurations using natural weighting (beam = 0$\farcs$28
$\times$ 0$\farcs$18, PA = $-64.3\arcdeg$), with an overlay of the bottom
panel of Figure 3 of Hodapp & Chini 2014, showing the IR continuum emission of
the SVS 13 star (light blue), the bubble traced by H2 S(1) (red) and the jet
traced by [Fe II] (green). The plus sign marks the optical position of the
star SVS 13 from Gaia EDR3 (Gaia Collaboration et al., 2021) extrapolated to
the epoch of the Ka-band observations (see text). The crosses mark the
positions of the two radio sources. Note that the extension to the SE of the
(free-free) radio emission traces the direction of the [Fe II] jet, and
possibly part of the NE wall of the H2 S(1) bubble. Contour levels are $-$3,
3, 5, 7, 10, 14, 19, 25, 35, 50, 70, and 95 times the rms of the image. The
synthesized beam and the rms are shown in the bottom right and left corners,
respectively, of the image.

Therefore, it seems that the 9 mm continuum emission of our observations
traces the free-free emission of the ionized radio jet associated with the
high excitation shock-excited [Fe II] jet. Additionally, in Figure 2, there is
a hint that the 9 mm emission extending to the SE of VLA 4B, further bends
toward the east, partially tracing the NE edge of the bubble observed by
Hodapp & Chini (2014) in the lower excitation emission of shock-excited H2.
Unfortunately, our A-configuration data alone are not sensitive enough to
trace further details of this jet-bubble system, but additional, more
sensitive VLA A-configuration observations could reveal with a high degree of
detail the emission of the ionized material associated with this shock
interaction, and trace the jet back to its origin to clarify whether the
optical/IR star SVS 13 is the same object as the radio source VLA 4B.

Figure 3: Overlay of the enhanced angular resolution images of SVS 13 at Q
band (7 mm; red contours) and Ka band (9 mm; black contours) obtained from the
A configuration VLA data using robust = 0 weighting. Only baselines $>$1000
k$\lambda$ (beam = 0$\farcs$036 $\times$ 0$\farcs$032, PA = $89.2\arcdeg$; rms
= 34 $\mu$Jy beam-1) and $>$600 k$\lambda$ (beam = 0$\farcs$066 $\times$
0$\farcs$044, PA = $-76.3\arcdeg$; rms = 26 $\mu$Jy beam-1) have been used for
the Q-band and Ka-band images, respectively. The beam-deconvolved sizes and
position angles obtained from Gaussian fits to the source VLA 4B in these
images are 60.3$\pm$3.5 mas $\times$ 49.0$\pm$3.1 mas, PA = 51${}^{\circ}\pm
25^{\circ}$ (7 mm) and 71.3$\pm$8.8 mas $\times$ 53.4$\pm$7.5 mas, PA =
129${}^{\circ}\pm 14^{\circ}$ (9 mm). Contour levels are $-$3, 3, 5, 7, 10,
13, and 17 times the rms of the images. The synthesized beams are shown in the
bottom left corner of the image.

Figure 3 shows an overlay of the Q (7 mm) and Ka (9 mm) band images obtained
with robust=0 weighting. This improved angular resolution image illustrates
the change in the orientation of the major axis of the emission associated
with VLA 4B, from PA = 51${}^{\circ}\pm 25^{\circ}$ (7 mm) to PA =
129${}^{\circ}\pm 14^{\circ}$ (9 mm). Thus, the morphology at $<0\farcs 1$
scales further suggests that the 9 mm image of VLA 4B is dominated primarily
by free-free emission arising from the jet, observed in [Fe II] at PA = 145∘
(Hodapp & Chini, 2014), while the 7 mm image is dominated by the dust emission
of its associated, roughly perpendicular, circumstellar disk. VLA 4A remains
unresolved at all bands in Figures 1 and 3, and a similar morphological
analysis cannot be performed.

####  3.1.3 Proper motions

Proper motions of SVS 13 and nearby radio sources have been reported by
Carrasco-González et al. (2008) and by Lefèvre et al. (2017). VLA data are
currently available for SVS 13 spanning nearly 30 years, and with a high
enough angular resolution to separate the two components, VLA 4A and VLA 4B,
of the binary. In Table 3.1.1 we list the positions of the two radio sources
obtained from the literature and from the data presented in this paper. From
these data we can obtain more accurate proper motions and try to measure
differences in the relative positions of the two sources to trace, or at least
to constrain, their orbital motion. In Figure 4 we plot the right ascension
and declination coordinates of VLA 4A and VLA 4B as a function of the epoch of
observation, from 1989 to 2017. In addition to the cm VLA data we use the ALMA
positions at 0.9 mm obtained from an image made using baselines $>$750
k$\lambda$ in order to minimize the contribution of the extended emission (see
Fig. 8 in Section 3.2.2). Table 4 gives the measured proper motions as derived
from the linear least-square fit to the data shown in Figure 4. The absolute
proper motions are 12.33$\pm$0.52 mas yr-1 at PA = $138.2\pm 2.4^{\circ}$ for
source VLA 4A, and 13.77$\pm$0.38 mas yr-1 at PA = $141.6\pm 1.5^{\circ}$ for
source VLA 4B. These values are in agreement with previous estimates
(Carrasco-González et al., 2008; Lefèvre et al., 2017) within the
uncertainties, but our values are more accurate and allow us to search for
evidence of possible orbital motions.

Figure 4: J2000 right ascension and declination of the two radio components of SVS 13, VLA 4A (blue) and VLA 4B (red), as a function of the epoch of observation. The positions used in this plot are listed in Table 3.1.1. The proper motions are obtained through linear least-squares fits (dashed lines) to the radio positions. The near-IR 2MASS position of SVS 13, $\alpha$(J2000) = 03h29m03.759${}^{s}\pm 0.005^{s}$, $\delta$(J2000) = 31∘16′03$\farcs 99\pm 0\farcs 07$, at epoch J1999.90 (Skrutskie et al., 2006) and the optical Gaia EDR3 position, $\alpha$(J2000) = 03h29m03.77088${}^{s}\pm 0.00002^{s}$, $\delta$(J2000) = 31∘16′03$\farcs 7237\pm 0\farcs 0002$, at epoch J2016.0 (Gaia Collaboration et al., 2021) are also plotted (green) as a reference. The Gaia EDR3 proper motions (Table 4) have been plotted as green dashed lines. For consistency, the Gaia EDR3 coordinates, originally referred to ICRS have been transformed to J2000 using the task ICRSToJ2000 of the AnalysisUtils Python package.  Table 4: Proper Motions and 3D Velocities  | Proper MotionsaaProper motions in the plane of the sky, derived from the linear least-square fit to the data shown in Table 3.1.1 and Figure 4. | 3D VelocitybbFull 3D velocity. $V_{\alpha}$ and $V_{\delta}$ are the velocity components in the direction of the right ascension and declination axes, respectively, obtained from the proper motions in the plane of the sky, assuming a distance of 300 pc; $V_{\rm LSR}$ is the line-of-sight velocity component, obtained from our ALMA molecular line emission data toward each component (Section 3.4.2).
---|---|---
| $\mu_{\alpha}\cos\delta$ | $\mu_{\delta}$ | $V_{\alpha}$ | $V_{\delta}$ | $V_{\rm LSR}$
| (mas yr-1) | (mas yr-1) | (km s-1) | (km s-1) | (km s-1)
VLA 4A | 8.22$\pm$0.63 |  $-$9.19$\pm$0.41 | 11.69$\pm$0.90 |  $-$13.07$\pm$0.58 | 7.36$\pm$0.21
VLA 4B | 8.56$\pm$0.34 |  $-$10.79$\pm$0.44 | 12.17$\pm$0.48 |  $-$15.34$\pm$0.58 | 9.33$\pm$0.18
Gaia EDR3 | 9.98$\pm$0.26 |  $-$13.02$\pm$0.20 | 14.19$\pm$0.37 |  $-$18.52$\pm$0.28 |  $\cdots$
VLA (4B$-$4A)ccMotion of VLA 4B relative to VLA 4A. | 0.34$\pm$0.72 |  $-$1.60$\pm$0.58 | 0.48$\pm$1.02 |  $-$2.28$\pm$0.82 | 1.97$\pm$0.28

Figure 4 and Table 4 show indications that the declination of VLA 4B relative
to VLA 4A decreases with time (VLA 4B moving to the south relative to VLA 4A)
at a rate of 1.60$\pm$0.58 mas yr-1 (2.28$\pm$0.82 km s-1) with a total
displacement of 48.0$\pm$17.4 mas (14.40$\pm$5.22 au) in 30 yr, while the rate
and the total displacement in right ascension is considerably smaller,
0.34$\pm$0.72 mas yr-1 (0.48$\pm$1.02 km s-1) with a total displacement of
10.2$\pm$21.5 mas (3.06$\pm$6.44 au). Relative positions should not be
affected by errors in absolute astrometry from one epoch to another, and a
high precision can be reached, in principle. However, small differences in the
protostellar positions are difficult to measure from their cm radio emission
because these sources usually trace free-free jets that can suffer apparent
changes in their position because of morphological changes due to the ejection
of new ionized knots (see, e.g., the multiepoch study of the radio sources in
NGC 2071 by Carrasco-González et al. 2012). In this respect, positions
obtained at Q band and shorter wavelengths (tracing the dust contribution of
the disk) are more stable and provide more reliable proper motion
measurements. This has been noticed, for example, in the well-studied A1-A2
Class 0 close binary of the IRAS 16293$-$2422 triple system, whose proper
motions are affected by the positional changes due to the ejection of ionized
knots at cm wavelengths, while positions from data at shorter wavelengths
appear more reliable (see discussion in Loinard et al. 2007; Hernández-Gómez
et al. 2019; Maureira et al. 2020). In the case of SVS 13, the consistent
departure of the position of VLA 4A in declination from the overall motion
(Fig. 4) suggests that such an ejection event may have taken place around
epoch J2000, resulting in a reversal of the relative positions of the two
sources in the declination axis. For a better visualization, we plot in Figure
5 the positions of the two sources after subtraction of the proper motion
derived from the linear least-squares fit to the VLA 4A positions (blue dashed
lines in Fig. 4), as a function of time. As this figure better illustrates,
despite the anomaly around epoch J2000, the currently available data show a
systematic trend in the relative proper motions, with VLA 4B moving to the
south relative to VLA 4A. This result deserves further investigation and can
be easily confirmed by future Q-band VLA and mm/sub-mm ALMA observations at
high angular resolution.

Figure 5:  Residual offsets in right ascension and declination as a function
of time for each radio component of SVS 13, VLA 4A (blue) and VLA 4B (red),
relative to the fitted linear proper motion of VLA 4A (blue dashed lines in
Fig. 4). An offset of $0\farcs 2$ has been added to the difference in
declination of VLA 4B for improved visualization. The dashed lines correspond
to the fits obtained in Fig. 4. The figure is intended to better illustrate
the relative motion between the two sources using a co-moving reference system
that is reasonably independent of sudden changes in the position of one of
them (e.g., produced by ejection events, as it possibly occurred in source VLA
4A around epoch J2000). The center of mass would be a natural reference for
relative motions, but is not known here. The use of VLA 4A instead of VLA 4B
as a reference is arbitrary.

The observed trend in the proper motions, with VLA 4B moving to the south
relative to VLA 4A, corresponds to a counter-clockwise orbital motion as seen
(projected on the plane of the sky) from the Earth. Since the observed blue-
shifted jet and molecular outflow lobe point southeast, and assuming that its
associated disk in VLA 4B is roughly co-planar with the binary orbit, this
implies that the angular momentum vector points away from the molecular cloud,
toward the observer. That is, we are observing the orbital motion from the
head of the angular momentum vector. These results are in agreement with what
is inferred from the ALMA data (see Sects. 3.2.3, 3.4.1, and 3.4.2).

Absolute proper motions can, in principle, inform us about the relative masses
of the two components in the binary, since the center of mass is closer to the
more massive component, and therefore this component is expected to be the one
that shows smaller absolute proper motions. However, when the system has a
significant overall proper motion, which is not independently known, as in our
case, disentangling both types of motions becomes extremely challenging unless
the data cover a significant fraction of the orbital period (e.g., in
Carrasco-González et al. 2009; Ichikawa et al. 2021). We note that measuring
the difference in the absolute proper motions of the two components of the
A1-A2 IRAS 16293$-$2422 binary is insufficient to derive the relative masses
and that resolved gas kinematics is required for this purpose (see Maureira et
al. 2020 for details).

In Figure 4 we also plot, as a reference, the proper motions of the optical
position of the SVS 13 star reported by Gaia (green dashed line), together
with the Gaia position (epoch J2016.0) and the 2MASS IR position (epoch
J1999.90) plotted as green points; the 2MASS position is consistent with the
Gaia proper motions but the 2MASS positional errors are large and the
agreement is not significant enough. We note that there is a significant
difference between the radio proper motions of VLA 4B and the Gaia EDR3 proper
motions of SVS 13 (Table 4). This difference is 1.42$\pm$0.43 mas yr-1 in
right ascension and $-$2.23$\pm$0.46 mas yr-1 in declination, which is larger
than the dispersion in proper motions in the whole NGC 1333 region,
0.92$\pm$0.05 mas yr-1 in right ascension and 0.60$\pm$0.03 mas yr-1 in
declination, as estimated by Ortiz-León et al. (2018). Thus, the observed
differences in position and proper motions between VLA 4B and the optical star
SVS 13 (Tables 3.1.1 and 4, and Fig. 4) would be compatible with the two being
independent objects belonging to the same molecular cloud. However, the
uncertainties of the Gaia EDR3 catalog may possibly be underestimated in the
case of SVS 13 (see below). Indeed, the position, proper motions and parallax
given in Gaia EDR3 for SVS 13 fall outside the uncertainty range of the
previous values in the DR2 (Gaia Collaboration et al., 2018).

A reasonable explanation for this discrepancy is that what is seen at optical
wavelengths is not the direct light of the SVS 13 star but scattered light
from its circumstellar environment, namely, from the walls of a cavity created
by the jet imaged by Hodapp & Chini (2014). It has been known for a long time
that jets are frequently driven by embedded stars, whose actual position is
better traced by the cm emission peak, while the observed optical intensity
peak appears displaced by a fraction of arcsec (e.g., HL Tau and HH 30:
Weintraub et al. 1995; Ray et al. 1996). The observed difference of $0\farcs
07$ between the optical and radio positions of SVS 13 can be easily accounted
by this fact. Indeed, HST images of SVS 13 suggest that what is seen at
optical wavelengths is not a point-like star but an elongated source
consistent with scattered light from a cavity carved by the jet (G. Blázquez-
Calero et al., in prep.). Although Gaia provides extremely accurate positions
for point-like stellar sources, it does not work well when substructure is
present at scales below $0\farcs 7$, and works very poorly below $0\farcs 4$
(Fabricius et al., 2021). Therefore, Gaia cannot identify the presence of
extended emission below these scales, which could evidence the non-stellar
nature of some sources. For extended sources this would also result in
increased errors in the Gaia-determined proper motions and parallax, as may be
occurring in SVS 13.

We also note that the relative proper motion of the Gaia star with respect to
the radio source VLA 4B is oriented along PA = 147∘, which is very close to
the [Fe II] jet PA of 145∘. This PA agreement makes it even less likely that
the Gaia position traces a different object, since optical proper motion along
jet axis could be easily caused by a tiny change in jet/cavity brightness
distribution (e.g., because of the ejection of new knot, as it is also
frequently observed in free-free radio jets; see discussion above).

In summary, the Gaia results show that the optical position of SVS 13 is very
close to VLA 4B and different from that of VLA 4A, definitely favoring VLA 4B,
and not VLA 4A, as the radio counterpart of the optically identified star SVS
13. However, the optical and radio positions differ by an amount larger than
the formal uncertainties, and the possibility that they trace different
objects (and, therefore, a triple system is present) cannot be completely
discarded. We favor the identification of the optical star and the radio
source VLA 4B as a single object, and we attribute the small difference in the
observed optical/radio positions to the fact that SVS 13 is an embedded object
(likely a Class I source) whose direct light is not optically visible but is
scattered in its circumstellar environment. As a consequence, its true
position is better traced by the radio emission, as it has already been
observed to occur in other cases, such as the examples mentioned in a previous
paragraph.

###  3.2 ALMA Continuum Results

####  3.2.1 Millimeter Continuum Sources

Figure 6: ALMA image of the continuum emission at 0.9 mm of a region
$\sim$18′′$\times$20′′ in size in the vicinity of SVS 13. The region shown
corresponds essentially to the ALMA field of view (the FWHM of the primary
beam is $18^{\prime\prime}$, centered on SVS 13), and includes the continuum
dust emission from the nearby YSOs VLA 20, VLA 3, and VLA 17 (MMS 2 = SVS
13B), whose positions are given in Table 3. A close-up of these nearby sources
is shown as an inset in the bottom-left corner. For better visualization, the
wide-field image has not been corrected for the primary beam response, but the
values given below for the contour levels have been corrected. Contours are
$-$12, 12, 24, 40, 56, 72, 88, 120, 200, 400, 800, and 1200 times 0.075 mJy
beam-1, the rms at the center of the image (SVS 13); $-$6, 6, 12, 20, 30, 50,
70, 100, 150, and 200 times 0.094 mJy beam-1 (VLA 3); $-$3, 3, 6, 12, 20, 30,
and 40 times 0.135 mJy beam-1 (VLA 20); and $-$6, 6, 12, 20, 30, 50, and 70
times 0.75 mJy beam-1 (VLA 17). The synthesized beam ($0\farcs 170\times
0\farcs 084$, PA = $-2.3\arcdeg$) is shown in the bottom right corner of the
images.

Figure 6 illustrates the 0.9 mm dust continuum emission detected in the ALMA
field of view (the FWHM of the ALMA primary beam is $\sim$18′′, centered on
SVS 13). The image is quite similar to the 1.3 mm ALMA image obtained by Tobin
et al. (2018) with a larger synthesized beam. Our field of view encloses the
SVS 13 close binary, which stands out prominently as the strongest and more
extended continuum source near the center of the field, but also dust emission
from the nearby YSOs VLA 20, VLA 3, and VLA 17 (MMS 2 = SVS 13B) is detected.
All these sources are detected in our VLA observations (Sect. 3.1.1). Sources
VLA 3, VLA 17, and VLA 20 were also detected in the 1.3 mm ALMA images of
Tobin et al. (2018).

We note that, although the Class 0 object MMS 2 (= VLA 17) has been considered
in the past a (wide) binary companion of SVS 13 (e.g., Bachiller et al. 2000;
hence its alternative name, SVS 13B, and the use by some authors of SVS 13A to
refer to the star SVS 13), at first glance, there is no direct connection
between these two objects in the continuum emission images, nor with the rest
of the sources which, indeed, are closer in projection to SVS 13 than VLA 17
(MMS 2 = SVS 13B). For the sake of clarity in the identification of the true
(close) binary SVS 13, we will avoid, as far as possible, this terminology in
terms of SVS 13A/B.

The positions and flux densities of the ALMA sources detected in the field are
given in Tables 3.1.1, 3.1.1, and 3. As noted in Section 3.1.1, in this paper
we focus our discussion on SVS 13. A study of the proper motions and SEDs of
the remaining sources will be presented in a forthcoming paper (A. K. Diaz-
Rodriguez, in prep.).

####  3.2.2 Circumstellar and Circumbinary Dust

Figure 7: ALMA image at 0.9 mm toward SVS 13. The two emission peaks are
associated with the previously identified cm sources, VLA 4A (the western one)
and VLA 4B (the eastern one). The image has been made from the highest angular
resolution data (Cycle 3) using all the baselines and Briggs weighting with
robust parameter of 0.5. Contours are $-6$, 6, 12, 24, 40, 56, 72, 88, 120,
200, 400, 800, and 1200 times 0.075 mJy beam-1, the rms of the image. Contour
levels in units of beam-averaged brightness temperature are $-$0.3, 0.3, 0.6,
1.3, 2.2, 3.0, 3.9, 4.8, 6.5, 10.8, 21.6, 43.2, and 64.9 K. The synthesized
beam ($0\farcs 170\times 0\farcs 084$, PA = $-2.3\arcdeg$) is shown in the
bottom left corner of the image. The image has been corrected by the primary
beam response.

Figure 7 shows a close-up toward SVS 13 of the 0.9 mm continuum emission
imaged with ALMA. The image shows two strong dust peaks, coincident with the
sources VLA 4A and VLA 4B, surrounded by extended emission in the form of
spirals. A similar morphology can be seen in the ALMA image at 1.3 mm obtained
by Tobin et al. (2018) (see their Fig. 1). We interpret the two compact
emission peaks as tracing circumstellar disks of dust associated with the two
protostars, VLA 4A and VLA 4B, at scales $\lesssim 0\farcs 2$. A large
fraction of the 0.9 mm continuum emission is distributed at larger scales
(within $\lesssim 1\farcs 5$ from the protostars), probably tracing dust from
a circumbinary disk that is still in the process of forming. The dust emission
shows an inner and an outer spiral arms, which are better differentiated to
the east of the binary. Most of this dust emission seems to converge toward
the position of source VLA 4A through the outer arm. For the inner arm it is
unclear whether it is connected to VLA 4A or to an intermediate position in
between VLA 4A and VLA B. Overlapped with the emission from these two spiral
features, a slight emission enhancement is observed extending to the southeast
of VLA 4A, which appears to follow a straight path along a PA=160∘ angle,
similar to the PA of the CO/H2 jet that appears to originate in VLA 4B (Hodapp
& Chini 2014; Lefèvre et al. 2017; G. Blázquez-Calero et al. in prep.). With
the present data, it is unclear whether this is just the result of the
superposition of the emission from the spiraling features observed south of
the protostars, or whether it could potentially be related to the southern
lobe of the outflow. Further analysis is needed to clarify the nature of this
feature. There is also some weaker fluffy emission extending further out to
the southwest, which seems to be associated with the ambient high density
molecular core that continues to the west and southwest of SVS 13, as is seen
in VLA ammonia images (Rudolph et al., 2001; Diaz-Rodriguez, 2021).

Figure 8: Close-up of the 0.9 mm continuum emission near the two components of
the SVS 13 protobinary. The image has been obtained from the same data and
with the same weighting as in Fig. 7, but using only baselines $>$750
k$\lambda$ to filter out the extended emission. The black plus sign marks the
J2000 position of the Gaia EDR3 (Gaia Collaboration et al., 2021) optical
source at the epoch of our observations (see Section 3.1.3, Fig. 4). Contours
are $-$6, 6, 14, 30, 60, 100, and 140 times 0.40 mJy beam-1, the rms of the
image. The synthesized beam ($0\farcs 120\times 0\farcs 053$, PA =
$-2.4\arcdeg$) is shown in the bottom left corner of the image.

If the extended emission is suppressed by restricting the uv-range of the data
($>$750 k$\lambda$), we obtain the image shown in Figure 8. This image is
expected to preferentially trace the compact circumstellar dust associated
with each protostar. From this image, accurate positions for the two sources
are obtained, coincident within $\sim$10 mas with the positions inferred from
the cm data (see Table 3.1.1). It is interesting to note that, although in the
image using all the baselines (Fig. 7) the 0.9 mm emission concentrates toward
the western side, apparently mainly associated with VLA 4A, when excluding the
extended emission the peak toward VLA 4B is brighter than the peak toward VLA
4A. This suggests that there is a larger amount of emitting dust surrounding
source VLA 4A, extending over a larger region than in VLA 4B. The emitting
dust associated with VLA 4B appears to be distributed on a more compact scale
and, at this scale, the column density is higher and/or the dust is hotter
than near VLA 4A.

Specifically, in the 0.9 mm image obtained using only baselines $>$750
k$\lambda$ (Fig. 8), a Gaussian fit gives a flux density toward VLA 4A of
82.9$\pm$3.0 mJy (Table 5), arising from a compact region with a deconvolved
size (FWHM of the major axis) of 80.6$\pm$3.8 mas (or 24$\pm$1 au), at
PA=83$\arcdeg$$\pm$8$\arcdeg$. Given the small size obtained, we interpret
this emission as dust emission from an inclined circumstellar disk. The total
flux density that is actually associated with VLA 4A is difficult to estimate,
as there is not a clearly defined boundary in the image obtained using all the
baselines (Fig. 7), but we estimate it as $\sim$220$\pm$25 mJy (within a
region of $\sim$100 au in size), from this image.

Table 5: Parameters of the 0.9 mm Circumstellar Dust EmissionaaSizes, position angles and flux densities are obtained from elliptical Gaussian fits. Uncertainties in these parameters only account for the formal error of the fit. Sizes correspond to the deconvolved FWHM of the major and minor axes. Since the sources are elongated roughly along the direction of the minor axis of the beam, where the angular resolution is higher, we consider that the values of the size of the major axis are reliable but the uncertainty of the minor axis can be quite underestimated. |  $>$750 k$\lambda$ |  $>$500 k$\lambda$
---|---|---
Parameter | VLA 4A | VLA 4B | Beam | VLA 4A | VLA 4B | Beam
Size (mas $\times$ mas) | 80.6$\pm$3.8$\times$34.2$\pm$20.2 | 61.3$\pm$2.6$\times$45.7$\pm$8.5 | 120$\times$53 | 99.7$\pm$4.0$\times$53.5$\pm$13.3 | 65.8$\pm$3.5$\times$45.4$\pm$15.1 | 150$\times$60
Size (au $\times$ au) | 24.2$\pm$1.1$\times$10.3$\pm$6.1 | 18.4$\pm$0.8$\times$13.7$\pm$2.5 | 36$\times$16 | 29.9$\pm$1.2$\times$16.1$\pm$4.0 | 19.7$\pm$1.1$\times$13.6$\pm$4.5 | 45$\times$18
PA ($\arcdeg$) | 83$\pm$8 | 93$\pm$17 |  $-$2 | 80$\pm$7 | 93$\pm$20 |  $-$2
$i$ ($\arcdeg$)bbInclination angle of the disk plane with respect to the plane of the sky, inferred from the aspect ratio. Uncertainties are derived following the procedures explained in Appendix C. Only the formal errors of the Gaussian fits have been taken into account in the error calculation. Therefore, the uncertainties of the inclination angles can be notably underestimated. | 64$\pm$16 | 40$\pm$12 |  | 57$\pm$9 | 43$\pm$18 |
$S_{\nu}$ (mJy) | 82.9$\pm$3.0 | 109.7$\pm$2.8 |  | 121.6$\pm$3.4 | 122.0$\pm$3.5 |
$T_{\rm dust}$ (K)ccOptically thick lower limit to the dust temperature obtained from $S_{\nu}/\Omega_{S}=B_{\nu}(T_{\rm dust})$, where $B_{\nu}$ is the Planck function and $S_{\nu}$ the flux density of the source, assumed to be a Gaussian of solid angle $\Omega_{S}=\pi\,a\,b\,/(4\ln 2)$, where $a$ and $b$ are its FWHM dimensions. Uncertainties in the $S_{\nu}/\Omega_{S}$ ratio are estimated as explained in Appendix C. | 430$\pm$190 | 430$\pm$80 |  | 260$\pm$60 | 480$\pm$145 |

On the other hand, in the image obtained using only baselines $>$750
k$\lambda$, the flux density toward VLA 4B is 109.7$\pm$2.8 mJy, coming from a
region with a deconvolved size (FWHM of the major axis) of 61.3$\pm$2.6 mas
(or 18$\pm$1 au), at PA=93$\arcdeg$$\pm$17$\arcdeg$ (Table 5). We interpret
this 0.9 mm emission as dust emission from a circumstellar disk associated
with VLA 4B. From the image obtained using all the baselines, we estimate that
the total flux density of the dust associated with VLA 4B (within a region of
$\sim$80 au in size) is $\sim$160$\pm$16 mJy.

In summary, the flux density of the whole dust emission associated with VLA 4B
appears to be $\sim$ 30% lower than in VLA 4A, but the emission at very small
scales (within a radius of $\sim$9 au) toward the position of VLA 4B is
$\sim$30-50% higher than toward VLA 4A. This suggests that, while VLA 4A
appears to be associated with a (total) larger mass of dust in its
surroundings, the dust column density (and therefore the dust opacity) at very
small scales from the protostar is significantly higher toward VLA 4B.

We note that the synthesized beam of the $>$750 k$\lambda$ image is elongated
in the north-south direction, with a minor axis of 50 mas in the east-west
direction, very close to the direction of the major axes of the two disks
(Table 5). Therefore, we consider that, although the sizes of the minor axes
are rather uncertain, the deconvolved sizes of the major axes of the disks are
quite well determined and reliable. Another indication of the significance of
the disk sizes is the fact that we obtain comparable values of the sizes and
flux densities by fitting an image obtained with a less restrictive uv-range,
using baselines $>$500 k$\lambda$, and a somewhat larger beam (Table 5). In
this last case, we obtain a major axis size of 99.7$\pm$4.0 mas (30$\pm$1 au),
at PA=80$\arcdeg$$\pm$7$\arcdeg$, for VLA 4A; and 65.8$\pm$3.5 mas (20$\pm$1
au), at PA=93$\arcdeg$$\pm$20$\arcdeg$, for VLA 4B (Table 5). This suggests
that there is a fairly well-defined separation between the boundary of the
circumstellar disks and the surrounding envelope, and that the angular sizes
of $\sim$80 and $\sim$60 mas (radii of 12 au and 9 au, respectively) for the
circumstellar disks of dust obtained from the $>$750 k$\lambda$ images are
therefore significant. The position angles of the major axes are also
consistently determined to be $\sim$90∘, further supporting the significance
of the fits.

However, the sizes of the disk minor axes are poorly constrained, as the beam
is highly elongated in the north-south direction. From the source aspect
ratios, and assuming they trace circular disks, the inclination angle of the
disk plane with respect to the plane of the sky can be calculated (see results
in Table 5, and Appendix C for the procedures followed in the error
calculations). But we want to caution that the values of the disk inclination
we have obtained are very uncertain, probably considerably more than indicated
by the uncertainties given in the table, as these uncertainties only take into
account the formal errors of the Gaussian fits.

Thus, from the results described above, we conclude that the radii of the
circumstellar disks of dust associated with VLA 4A and VLA 4B are likely of
the order of $\sim$12 au and $\sim$9 au, respectively. These values of the
radii are smaller than those typically found for disks around single stars. In
binary systems, tidal truncation to disk radii of $\sim$1/3 of the binary
separation is expected (Papaloizou & Pringle, 1977; Artymowicz & Lubow, 1994),
resulting in smaller disks. Since the projected distance between the two
components of the binary is 90 au, disks with radii of $\sim$30 au are
expected to result from tidal truncation assuming circular orbits. Thus, the
circumstellar disks in SVS 13 seem to be somewhat smaller than expected, if
only tidal truncation is at work and unless the orbit of one (or both) of the
binary components has a high eccentricity. In any case, we note that we are
comparing the values of the radii at half the maximum intensity, but the dust
emission certainly extends beyond this radius; it should also be noted that,
as has been observed in other disks, the gas is usually distributed up to
larger radii than the dust. Thus, in the case of the circumstellar disks of
SVS 13, the deviation from expectations is not severe and is roughly
consistent with a tidal truncation origin. An extreme case of a still smaller
disk was reported by Osorio et al. (2016) in the XZ Tau A/B young system
(separation between components of 42 au), where a disk of dust with a radius
of only 3 au was imaged.

The small sizes and large flux densities measured for the 0.9 mm emission
toward VLA 4A and VLA 4B imply relatively high brightness temperatures, and
therefore warm dust. From the flux densities and sizes measured in the highest
resolution images described above (Table 5), and assuming optically thick
emission (which gives a lower limit for the dust temperature), we obtain dust
temperatures of the order of $\sim$300 K toward VLA 4A and $\sim$450 K toward
VLA 4B. Since implausibly higher temperatures would be obtained if the
emission were optically thin, and because there is independent evidence that
the dust emission is optically thick toward the protostars (Sect. 3.4.1), we
conclude that the dust emission is very optically thick in the SVS 13
circumstellar disks, and that the dust reaches kinetic temperatures of the
order of $\sim$300-450 K at radii of $\sim$10 au from the protostars.

Under certain conditions, and at these small radii, these relatively high dust
temperature values are not unrealistic. Indeed, similar values have been
inferred from the modeling of the SEDs and mm images of other low-mass
protobinaries such as L1551 IRS 5 (Osorio et al., 2003). In general, flared
disks are warmer at the disk surface and cooler near the equatorial plane. But
in sources with relatively high mass accretion rates, such as the
circumstellar disks of L1551 IRS 5, models predict temperatures on the order
of 300-400 K near the equatorial plane, at radii of 5-10 au (Osorio et al.,
2003). The model predictions on the configuration and physical properties of
the L1551 IRS 5 system by Osorio et al. (2003) have been nicely confirmed
through recent high-resolution ($\sim$$0\farcs 1$) ALMA observations at 0.9 mm
by Takakuwa et al. (2020). In particular, these authors measured brightness
temperatures reaching $\sim$260 K, implying high dust temperatures in the
L1551 IRS 5 circumstellar disks, similar to the values we infer in the
circumstellar disks of SVS 13. However, the current observational
uncertainties in SVS 13 remain large (e.g., the disk minor axes and, thus, the
solid angle of the sources are poorly constrained with the current angular
resolution) and, also, more specific model predictions, tailored to the
properties of the SVS 13 disks, are needed to test whether their properties
are compatible with a high accretion regime and which scales above the disk
midplane the optically thick 0.9 mm emission traces.

As one moves away from the protostars, the brightness temperature decreases
sharply, down to values of $\sim$80 K at scales of the order of the beam size
of the robust = 0.5 image shown in Figure 7. Since the temperatures inferred
from the molecular emission at these scales are $\sim$140 K (see Sect. 3.5.1),
we attribute the decrease in brightness temperature to the combined effect of
a true decrease of the kinetic temperature together with a decrease in the
optical depth and/or the beam filling factor.

Toward other positions, in the circumbinary material and in the spirals of SVS
13, we obtain dust brightness temperature values in the range 5-12 K. Since
the temperatures for the circumbinary disk inferred from the molecular
emission (see Sect. 3.5.1) are of the order of 140 K, we conclude that the 0.9
mm dust emission is optically thin (optical depth $\sim$0.05-0.1) across the
circumbinary structure and the spiral arms. This is in contrast with the
warmer temperatures and high optical depths inferred for the close environment
of the protostars.

####  3.2.3 Spirals of Dust

An outstanding feature of the ALMA continuum emission (Fig. 7) are the
prominent spiral arms associated with the SVS 13 binary. Trailing spiral arms
should point away from the direction of rotation, therefore implying an SVS 13
disk counter-clockwise rotation, as seen from the Earth. This is in agreement
with our conclusion about the sense of the orbital motion inferred from the
observed proper motions (Section 3.1.3) and with the analysis of radial
velocities from the molecular line emission (Section 3.4.1 and Fig. 12).

Figure 9: Same as in Fig. 7 but with the intensity scaled by a factor $\propto
r^{1.5}$, where $r$ is the projected distance to the midpoint between VLA 4A
and VLA 4B in arbitrary units, to compensate for the decreasing intensity and
to enhance the visibility of the spiral features at large distances. The
resulting intensity has been normalized to the maximum. No correction has been
applied for projection effects. Three spirals (S1, S2, and S3) are tentatively
identified in the image and marked with dashed lines. The plus signs mark the
positions of VLA 4A (western) and VLA 4B (eastern) sources.

To better emphasize the spiral structures, in Figure 9 we present the same 0.9
mm ALMA image as in Figure 7 but normalized by a factor $\propto r^{1.5}$,
where $r$ is the distance to the midpoint between VLA 4A and VLA 4B, to
compensate for the natural decrease of intensity with distance to the source.
We identify three features, S1, S2, and S3, in order of apparent intensity
(see Fig. 9). S1 seems to originate near the location of VLA 4A. For S2, it is
unclear whether it originates near VLA 4A or at an intermediate position
between VLA 4A and VLA 4B. S3, the weakest feature, seems to originate near
the location of VLA 4B. S1 and S2 are quite clearly defined spirals, but for
S3 it is unclear whether it is a spiral arm or a feature of the northern part
of the circumbinary disk.

Spirals have been observed in a number of protoplanetary disks and various
mechanisms, which can possibly act together, have been suggested to explain
their origin, including planet-disk interactions (e.g., Ogilvie & Lubow 2002;
Muto et al. 2012; Dong et al. 2015; Phuong et al. 2020), gravitational
instabilities (e.g., Lodato & Rice 2005) or shadows (Montesinos et al., 2016).
In the hydrodynamical simulations of the formation of binary systems, spirals
also appear as a natural outcome (e.g., Bate & Bonnell 1997; Bate 2018;
Matsumoto et al. 2019).

In particular, the overall evolution of the fiducial model in the Matsumoto et
al. (2019) simulations shows two spiral arms on the sides of the primary and
secondary stars, where the secondary arm is connected directly to the
circumsecondary disk, while the primary spiral arm runs around the
circumprimary disk and is connected to the bridge between the two
circumstellar disks. This result from the numerical simulations appears to be
broadly consistent with our observations of SVS 13 if the secondary were VLA
4A (the component with the more prominent spiral arm, S1) and the primary were
VLA 4B (with the S2 spiral surrounding this protostar and being connected to
some undetermined point between VLA 4A and VLA 4B). This scenario is indeed
consistent with the estimates of the individual protostellar masses obtained
in Sect. 3.5.3 and Appendices D and E since we obtain a VLA 4B mass higher
than that of VLA 4A. Unfortunately, our present results only provide rough
mass estimates. Future observations of accurate proper motions at mm
wavelengths, such as those outlined in Section 3.1.3, combined with small-
scale kinematics obtained from new molecular data with even improved angular
resolution, should be able to shed more light on the masses of the SVS 13
binary and to reliably identify the primary and secondary components. This
will make it possible to test the models of the evolution of a binary system
which is undergoing the development of a circumbinary disk and spiral arms, as
seems to be the case for SVS 13.

###  3.3 The Spectral Energy Distribution

Our VLA and ALMA observations angularly resolve the SVS 13 system over a range
of wavelengths from 6 cm to 0.9 mm. There are not many observational data on
close ($<$100 au) protobinary systems that angularly separate the emission of
the two components. Our data, combined with previous data from the literature
allow us to characterize in depth the nature of the emission and the physical
properties of this system. In Table 3.1.1 we list the measured flux densities
and in Figure 10 we plot the resulting SED.

Figure 10: Spectral energy distribution (SED) of SVS 13 from cm to mm
wavelengths. Data used are reported in Table 3.1.1; blue and red data points
represent the flux density toward the individual components (VLA 4A and VLA
4B), while green data points represent extended emission toward the whole
system (identified with footnotes $d$ and $f$ in Table 2, and labeled as “CB”
in this figure), tracing mainly the circumbinary dust. The spectral indices
($\alpha$, where $S_{\nu}\propto\nu^{\alpha}$) are obtained from least-squares
fits (in log-log scale) that are shown as solid lines. For the individual
components, VLA 4A and VLA 4B, the free-free spectral indices are obtained
using only data with a frequency below 40 GHz ($\lambda>7.5$ mm). The dashed
lines are extrapolations from our 7 mm measurements, assuming optically thin
dust emission with the same spectral index obtained in the circumbinary dust.
Note that, in the case of VLA 4B, the observed 0.9 mm flux density is well
below the optically thin extrapolation, indicating a high dust optical depth
at this wavelength (see Sect. 3.2.2).

Regarding the VLA data, we note that in 2016.9 there was an important decrease
of the X-band (3 cm) flux density of VLA 4A, by a factor close to 4 with
respect to previous dates ($\sim$1995-2014), which is atypical for thermal
free-free emission, where the flux density variability usually does not exceed
$\sim$20% (see Anglada et al. 2018). This variability is usually attributed to
sudden variations in the mass accretion rate, which trigger subsequent
variations in the outflow mass loss rate. But, in general, they imply sudden
increases in flux density, sometimes associated with the ejection of new
ionized knots, rather than sudden decreases. Further monitoring of the sources
with high sensitivity and angular resolution are necessary to improve the
previous data and to clarify the origin of the X-band variability. At shorter
wavelengths, $\sim$7 mm, a smaller increase in flux density, by a factor of
$\sim$1.5, took place in 2015 in both sources, but only in VLA 4B persisted in
the subsequent observations.

Anglada et al. (2004) noticed that the flux density of VLA 4B was much
stronger than that of VLA 4A at 7 mm, while the fluxes of the two sources were
comparable at longer wavelengths. We confirm this observational result, which
can also be seen in the data of Tychoniec et al. (2018) and Tobin et al.
(2016b). Since the observed 7 mm flux of VLA 4B coincided with the
extrapolation of the mm data available at that time (which did not separate
the binary), Anglada et al. (2004) attributed the excess of 7 mm emission to
dust and proposed that, while VLA 4B was associated with a bright
circumstellar disk of dust, VLA 4A was not. This observational result was
interpreted by Anglada et al. (2004), in terms of the Bate & Bonnell (1997)
theoretical simulations, as suggestive that the formation of the SVS 13 binary
system occurred via the accretion of material with low specific angular
momentum. In this case, the primary may have a large circumstellar disk, while
the secondary is essentially naked; thus, VLA 4B was proposed as the primary,
more massive, component and VLA 4A as the secondary. Also, according to Bate &
Bonnell (1997) the systems formed by this mechanism should not develop a
significant circumbinary disk.

Our present results clarify the previous conclusions by Anglada et al. (2004)
and shed new light on the interpretation of the cm and mm continuum data.
Despite the difficulties in distinguishing between circumstellar and
circumbinary contributions at such early stages of the formation process of a
binary system, as is the case of SVS 13, we identify three main dust
contributions. We see that both VLA 4A and VLA 4B protostars are associated
with circumstellar dust, although with different properties regarding its
distribution, as discussed in Section 3.2.2. We also clearly see that there is
dust (and gas; see Sect. 3.4) in between and surrounding the two protostars,
and with prominent spiral arms. We interpret this dust as material being
assembled into a proto-circumbinary disk, which we will refer to simply as
“the circumbinary disk”.

In Figure 10 we plot the SED in the cm/mm regimes, where we identify several
contributions, namely: (i) The free-free emission associated with VLA 4A (with
a spectral index of 1.0 at wavelengths longer than $\sim$7 mm), likely tracing
a thermal radio jet; (ii) The free-free emission associated with VLA 4B (that
dominates the emission of this object at long wavelengths), also consistent
with the presence of a thermal radio jet; (iii) The dust emission of the
circumstellar disk associated with VLA 4A (traced by the 0.9 mm blue data
point) that implies a change to a steeper spectral index, $\sim$2.5, at short
wavelengths; (iv) The circumstellar disk around VLA 4B (traced by the 0.9 mm
red data point), with a flux density consistent with the extrapolation of the
long-wavelength data with a single slope of 1.7. Since a spectral index $<$2
cannot be produced by dust emission alone, this value suggests a free-free
contribution at long wavelenghts and optically thick dust emission at short
wavelengths; (v) The dust emission of the circumbinary disk and its associated
spiral arms (green diamonds with a spectral index of 2.5) which dominates the
total emission observed at short wavelengths.

Based on the previous discussion on the morphology of the sources at Q and Ka
bands (Sect. 3.1.2) and on the dust brightness temperatures (Sect. 3.2.2), it
is reasonable to assume that the 9 mm emission is dominated by free-free
emission, while the 7 mm and 0.9 mm emissions are dominated by dust. The 0.9
mm emission in the circumbinary disk is optically thin dust, while the dust
emission is optically thick in the VLA 4A circumstellar disk, and very
optically thick in the disk of VLA 4B. The 7 mm emission of both circumstellar
disks is expected to be optically thin, but it can contain some residual free-
free contamination.

###  3.4 ALMA Line Emission Results

####  3.4.1 Distribution and Overall Kinematics of the Molecular Gas

In this Section, we are focusing on the molecular line emission that appears
to be more directly associated with the SVS 13 binary and associated disks,
both spatially (within a few arcsec) and in velocity (within a few km s-1 of
the adopted systemic velocity of the cloud, +8.5 km s-1; Rudolph et al. 2001;
Diaz-Rodriguez 2021), and that is most useful to infer the physical properties
of the system and its close environment. The ALMA observations also reveal
complex molecular structures at large scales, $\gtrsim$$10^{\prime\prime}$
($\gtrsim$3000 au), and high-velocity features likely associated with the blue
outflow lobe that will be discussed in forthcoming papers (A. K. Diaz-
Rodriguez et al., in prep., and G. Blázquez-Calero et al. in prep.).

Figure 11: Continuum-subtracted velocity-integrated emission (zeroth-order
moment; contours), intensity-weighted mean velocity (first-order moment; color
scale in top panels) and velocity dispersion (second-order moment; color scale
in bottom panels) of the CS (7-6), H13CN (4-3), and HC15N (4-3) molecular
lines observed with ALMA toward SVS 13. Contours are $-$3, 3, 5, 10, 15, 20,
30, 40, 50, 60, and 70 times the rms of each map (25 mJy beam-1 km s-1 for CS,
and 10 mJy beam-1 km s-1 for H13CN and HC15N). The synthesized beam (indicated
in the bottom left corner of the maps) is $0\farcs 15\times 0\farcs 08$, PA =
$-$2∘. The positions of VLA 4A (western) and VLA 4B (eastern) are marked with
plus signs.

Figure 11 shows the continuum-subtracted velocity-integrated emission (zeroth-
order moment), the intensity-weighted mean velocity (first-order moment) and
the velocity dispersion (second-order moment) of the CS (7-6), H13CN (4-3),
and HC15N (4-3) molecular lines. These lines are the centering lines of three
of the six spectral windows in our observations and have been selected as
representative of the emission most directly associated with the binary
system. As can be seen in the figures, the emission appears distributed
surrounding the binary, and we interpret it as mostly tracing circumbinary gas
from which the stars are being assembled. The integrated emission clearly
peaks toward source VLA 4A, while there is an apparent decrease of emission
toward source VLA 4B that we think is probably not real but the result of an
imperfect continuum subtraction, because of the non-additivity of optically
thick intensities. As discussed in Section 3.2.2, the dust continuum emission
toward VLA 4B is very compact and intense, so it probably has a very high
optical depth and is saturated, making it difficult for line emission to
emerge above the continuum.

The integrated intensity maps presented in Figure 11 (in contours) show a
roughly elliptical molecular emission ($\sim$$1\farcs 0\pm 0\farcs 1\times
0\farcs 8\pm 0\farcs 1$, or $\sim$300$\pm$30 au $\times$ 240$\pm$30 au in size
for CS, and slightly smaller, $\sim$$0\farcs 8\pm 0\farcs 1\times 0\farcs 6\pm
0\farcs 1$, or $\sim$240$\pm$30 au $\times$ 180$\pm$30 au for the HCN
isotopologues), with the major axis roughly E-W. If the real structure is
assumed to be roughly circular, the observed ellipticity implies an
inclination angle of the disk (angle of the disk rotation axis and the line-
of-sight, where 0∘ corresponds to face-on) of $i\simeq 38^{\circ}\pm
14^{\circ}$ (see Appendix C for uncertainty calculations).

Figure 12: Sketch of the proposed geometrical and kinematical configuration of
the SVS 13 protobinary system and its circumbinary material. The red and blue
colors represent the appropriate velocity shift of the gas with respect to the
systemic velocity with the inferred angular momentum axis also shown.

The first-order moment maps (in color in the top panels of Fig. 11) show a
clear velocity gradient roughly in the SW-NE direction, with blue-shifted
velocities toward the SW side. The ALMA maps obtained by Tobin et al. (2018)
for other molecules and lower frequency transitions show differences from one
molecule to another, but are roughly consistent with ours in the shape and in
the direction of the overall velocity gradient. However, our higher angular
resolution and better sensitivity allow us to perform a more detailed
analysis. Actually, our maps show a kind of “yin-yang”-shaped velocity
distribution that we attribute to the simultaneous presence of both infall and
rotation with comparable velocities (see Fig. 12). Given that the southern
lobe of the outflow associated with SVS 13 is blue-shifted (i.e., pointing
toward the observer; Bachiller et al. 2000; Raga et al. 2013; Hodapp & Chini
2014), and assuming that the circumbinary emission is mainly distributed in a
disk-like structure close to the equatorial plane of the circumstellar disk
associated with (and thus perpendicular to) the outflow, this implies that we
are observing these disks from “below”. In other words, we are observing the
circumbinary disk in such a way that its northern half is the one closest to
the observer. Thus, the observed velocity gradient, with the eastern side red-
shifted, implies a counter-clockwise rotation as seen from the observer. This
direction of the circumbinary disk rotation is in agreement with that inferred
from the observed orbital proper motions (Sect. 3.1.3) and the dust spirals in
the 0.9 mm continuum emission (see Sect. 3.2.3 and Figs. 7, 9). Given the
orientation of the disk-like structure, the presence of infall motions would
produce additional red-shifted emission in its northern part and blue-shifted
emission in its southern part, resulting in the observed “yin-yang”-shaped
velocity distribution (see Fig. 12). However, given that the northern lobe of
the bipolar outflow is red-shifted and the southern one is blue-shifted (e.g.,
Bachiller et al. 2000; Raga et al. 2013; Hodapp & Chini 2014), a contribution
to the observed velocity field from the outflow cannot be discarded.

The second-order moment maps (in color in Fig. 11 bottom panels) show a large
increase in velocity dispersion toward VLA 4A. The second-order moment
increases from a typical value of 1.2 km s-1 in the circumbinary structure up
to $\sim$3.2 km s-1 toward source VLA 4A. This local increase in the velocity
dispersion could be due to infalling motions, increased turbulence associated
with accretion shocks caused by the impact of infalling material onto the
disk, or unresolved rotation in the disk. The origin of the main spiral arm
(S1), near the position of VLA 4A, also strongly suggests that this object is
a major accretion focus (see Sect. 3.2.3). In the proximity of VLA 4B there is
only a modest increase in the velocity dispersion, with a local decrease
toward the position of the star. However, it is unclear whether this apparent
decrease in the second-order moment toward VLA 4B is real, with the turbulence
and/or accretion being smaller toward this object, or is an artifact due to
the suppression of line emission in the presence of very optically thick
continuum emission.

####  3.4.2 Detection and Identification of Molecular Transitions

Figure 22 shows the full 2.8 GHz spectral band, obtained over a region with
the size of the beam ($\sim$$0\farcs 16\times 0\farcs 08$, PA = 0$\arcdeg$)
toward source VLA 4A (blue) and toward source VLA 4B (red), as well as toward
the circumbinary material (CB) in a position centered between the two
protostars (black). It is evident from the figure that there is a shift in the
velocity of the centroids of the molecular lines, with the lines toward VLA 4A
being blue-shifted and the lines toward VLA 4B being red-shifted with respect
to the LSR velocity of the lines observed toward the central position CB. This
velocity shift is fully consistent with the overall velocity gradient observed
in the first-order moment shown in the top panels of Figure 11. By measuring
the LSR velocity centroid of 15 lines toward VLA 4A and VLA 4B, and 20 lines
toward CB, we estimate mean values of $V_{\rm LSR}$(VLA 4A) = 7.36$\pm$0.21 km
s-1, $V_{\rm LSR}$(VLA 4B) = 9.33$\pm$0.18 km s-1, and $V_{\rm LSR}$(CB) =
8.77$\pm$0.16 km s-1, where the quoted uncertainties correspond to the
standard deviation of these centroids among the measured lines (see also Table
4).

In general, the intensity of a given line is stronger toward the central (CB)
position than in the rest of the molecular structure, and lines observed
toward VLA 4A are stronger than toward VLA 4B. On the other hand, lines are
broader toward VLA 4A than toward other positions of the molecular structure,
as was already inferred from the second-order moment images shown in Figure 11
(bottom panels). As can also be seen in Figure 11, the lines are broader and
more intense toward the western part of the circumbinary structure than toward
the eastern one, and are narrower toward the southern part than the northern
one. Additionally, the intensity of the lines in general decreases as we move
farther from the center of the circumbinary disk. The strongest lines often
have a double peaked profile, characteristic of absorption, which is
particularly evident toward VLA 4A (see Fig. 22). This suggests that these
lines are optically thick, but it could also be indicative of two velocity
components. Lines are weaker toward VLA 4B, but the emission at the stellar
position could be saturated by the intense, and likely very optically thick
continuum emission, making the lines hardly detectable above the continuum.

To identify the detected transitions we compared the observed spectra with
synthetic ones modeled in CASSIS (Vastel et al., 2015) applying the LTE
radiative transfer equations to a homogeneous column of gas along the line-of-
sight, where the distribution of the population of all the rotational energy
levels for a given molecule is described by a single excitation temperature,
also known as the rotational temperature. To obtain an initial guess of the
expected intensities of the different transitions, in our model we considered
Tex = 100 K based on previous works (Belloche et al., 2020b), a source size of
$0\farcs 16$ (the region where the spectrum was taken), and an initial FWHM of
2 km s-1, and changed manually the species column density and the FWHM to
match the observed spectrum. We note that our goal during this procedure was
the identification of the lines, not to derive parameters from the model. In a
first round, we used the species in the “Lowmass hot corino” and “ISM” filters
of CASSIS. We only included in our model transitions with upper level energies
($E_{\rm up}$) below 500 K. To identify the remaining unknown lines after this
first stage, we included in the model any species with transitions close to
the unidentified lines and with upper level energies up to 700 K.

Our proposed identifications are listed in Table 8. We report 180 transitions
of 29 species, counting the different isotopologues, as well as 28
unidentified lines. Among the identified species, 13 are complex organic
molecules (COMs), seven of which are reported here for the first time in the
SVS 13 system (see further details in Section 3.6.1).

We note that in some cases (marked as “nLTE” in Table 8) the intensity of the
lines is not well reproduced by the model (differences $\gtrsim 50\%$), which
could be due to departures from LTE or other effects such as optical depth or
gradients.

####  3.4.3 Distribution and Kinematics of the Circumstellar Molecular Gas

We interpreted the presence of compact 0.9 mm continuum emission directly
associated with VLA 4A and VLA 4B as tracing circumstellar disks of dust with
sizes of $\sim$ 20 au associated with these objects (see Fig. 8 and discussion
in Section 3.2.2). On the other hand, the maps presented in Figure 11 and
discussed in Section 3.4.1 show that both the integrated intensity (zero-order
moment) and the velocity dispersion (second-order moment) of the observed
molecular emission clearly peak toward VLA 4A and (less clearly, likely
because of the high optical depth of the dust) tend to increase in the
proximity of VLA 4B. This suggests that the two protostars are also associated
with circumstellar molecular emission.

Figure 13: Continuum-subtracted velocity-integrated emission (zeroth-order
moment; contours) and intensity-weighted mean velocity (first-order moment;
color scale in top panels) of the ethylene glycol (aGg’-(CH2OH)2 conformer)
lines directly associated with component VLA 4A of SVS 13 that have been
detected with ALMA. Contours are $-$3, 3, 5, 9, 13, 18, 25, 30, and 35 times
the rms of each map (4.5 mJy beam-1 km s-1 for the left and middle top panels,
and 15 mJy beam-1 km s-1 for the right top panel). Synthesized beam (indicated
in the bottom left corner of the maps) is $0\farcs 15\times 0\farcs 08$, PA =
$-$2∘. The positions of VLA 4A (western) and VLA 4B (eastern) are marked with
plus signs. The E-W velocity gradients are clearly observed. We note that
there is a 3-$\sigma$ negative contour toward VLA 4B in the left and central
top panels (see also Fig. 14); it is unclear whether it is related to a weak
absorption signature or just the result of a somewhat imperfect continuum
subtraction. The lower panels show the corresponding position-velocity
diagrams along the E-W direction centered on the position of VLA 4A. Offsets
are in right ascension coordinates relative to VLA 4A. Transitions are
indicated in the upper panels. Note that the rightmost panels include the
emission of two blended transitions.  Figure 14: Continuum-subtracted emission
of the ethylene glycol aGg’-(CH2OH)2 conformer (34 4 31 1 - 33 4 30 0)
transition (left panels in Fig. 13), associated with the component VLA 4A of
SVS 13. The emission has been integrated over two velocity ranges, from 4.46
to 7.23 km s-1 (blue contours) and from 7.23 to 9.37 km s-1 (red contours), to
illustrate the two velocity components observed toward VLA 4A. Contours are
$-3$, 3, 4, 5, 7, 10, 15, 20, 25, and 30 times 2.5 mJy beam-1 km s-1. The
positions of VLA 4A (western) and VLA 4B (eastern) are marked with plus signs.
The synthesized beam is indicated in the bottom left corner.

This association of molecular gas at circumstellar scales is nicely
illustrated by some transitions of the ethylene glycol molecule which are only
detected toward VLA 4A, as Figure 13 shows. The top panels in Figure 13 show
the integrated emission of different transitions of ethylene glycol which is
clearly associated with VLA 4A at circumstellar scales of $\sim$30 au. These
panels also show the presence of a clear velocity gradient in the E-W
direction which is consistently detected in all these transitions (note that
the rightmost panels include two blended transitions). The bottom panels in
Figure 13 show position-velocity cuts along a PA = $90^{\circ}$ to better
unveil the structure of this gradient, which is suggestive of rotation motions
associated with VLA 4A. Figure 14 further illustrates the distribution of the
blue- and red-shifted emission, integrated over two separate velocity ranges,
of the transition displayed in the left panel of Figure 13.

###  3.5 Physical Parameters of the SVS 13 System

####  3.5.1 Temperature

We can use the detected molecular transitions to estimate the kinetic
temperature by inferring the excitation temperature of the rotational levels,
also known as the rotational temperature, $T_{\rm rot}$. For this, we used
CASSIS to build rotational diagrams to derive the temperature and the column
density, assuming optically thin lines and LTE. In this case, the upper-level
column density

$N_{u}=\frac{8\pi k\nu^{2}}{hc^{3}A_{ul}}\int{T_{\rm mb}\,dV},$ (1)

and the rotational temperature $T_{\rm rot}$ are related as

$\ln{\frac{N_{u}}{g_{u}}}=\ln{N_{\rm mol}}-\ln{Q\left(T_{\rm
rot}\right)}-\frac{E_{u}}{kT_{\rm rot}},$ (2)

where $k$ is Boltzmann’s constant, $h$ is Planck’s constant, $c$ is the speed
of light, $\nu$ is the frequency of the transition and $A_{ul}$ its Einstein
coefficient for spontaneous emission, $g_{u}$ is the statistical weight of the
upper level and $E_{u}$ its energy, $N_{\rm mol}$ is the total column density
of the molecule, $Q(T_{\rm rot})$ is the partition function for $T_{\rm rot}$,
and $\int{T_{\rm mb}\,dV}$ is the integrated main beam brightness temperature
of the line, which is derived from a Gaussian fit to the observed molecular
line.

We use only molecular lines that have been positively identified and well
reproduced by the LTE model, with energies that span over a relatively large
range, and that are non-blended, or, in case they are blended, that can be
effectively isolated by a multi-component Gaussian fitting taking into account
all the blended lines involved. These selection criteria allow us to populate
the rotational diagrams only for a few molecules (CH3OCHO, CH2DOH, and
aGg’-(CH2OH)2); in the case of aGg’-(CH2OH)2 this is possible only toward VLA
4A.

Figure 15: Rotational diagrams for methyl formate (CH3OCHO), deuterated
methanol (CH2DOH) and ethylene glycol (aGg’-(CH2OH)2), for VLA 4A (blue), VLA
4B (red), and the circumbinary disk (black), using spectra taken in beam-size
regions ($0\farcs 16\times 0\farcs 08$, PA = 0$\arcdeg$) toward each component
and toward the geometric center of the system (shown in Fig. 22). We only
include in the plot lines that have been positively identified and well
reproduced by the LTE model, that are non-blended, or, in case they are
blended, that we were able to isolate by a multi-component Gaussian fit. In
the case of aGg’-(CH2OH)2 this is possible only toward VLA 4A.  Table 6:
Rotational Temperatures, Molecular Column Densities and Abundances

| VLA 4A | VLA 4B | CBaaCentral position, between VLA 4A and VLA 4B, taken as representative of the circumbinary material.
---|---|---|---
|  $T_{\rm rot}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  $N_{\rm mol}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  |  $T_{\rm rot}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  $N_{\rm mol}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  |  $T_{\rm rot}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  $N_{\rm mol}$bbRotational temperature and molecular column density obtained from the rotational diagram fits shown in Fig. 15, except when indicated. |  |  $X_{\rm H_{2}}$ddMolecular abundance relative to H2. The adopted H2 column density is 1.8$\times$1024 cm-2, calculated from the dust emission, assuming optically thin emission, a dust temperature of 140 K, and a gas-to-dust ratio of 100. See text for further details.
Molecule | (K) | ($10^{17}$ cm-2) |  $X_{\rm CH_{3}OH}$ccMolecular abundance relative to CH3OH obtained from the CH2DOH column density, assuming a deuteration fraction D/H = 0.05, estimated as the ratio of the CH2DOH and CH3OH column densities toward CB (see notes e, g). | (K) | ($10^{17}$ cm-2) |  $X_{\rm CH_{3}OH}$ccMolecular abundance relative to CH3OH obtained from the CH2DOH column density, assuming a deuteration fraction D/H = 0.05, estimated as the ratio of the CH2DOH and CH3OH column densities toward CB (see notes e, g). | (K) | ($10^{17}$ cm-2) |  $X_{\rm CH_{3}OH}$ccMolecular abundance relative to CH3OH obtained from the CH2DOH column density, assuming a deuteration fraction D/H = 0.05, estimated as the ratio of the CH2DOH and CH3OH column densities toward CB (see notes e, g). | ($\times$$10^{-7}$)
CH3OCHO | 127$\pm$15 | 4.2$\pm$1.8 | 2.6$\times$$10^{-2}$ | 90$\pm$6 | 2.5$\pm$0.8 | 6.9$\times$$10^{-2}$ | 139$\pm$6 | 2.3$\pm$0.3 | 2.3$\times$$10^{-2}$ | 1.3
aGg’-(CH2OH)2 | 145$\pm$39 | 1.1$\pm$0.8 | 7.1$\times$$10^{-3}$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |
CH2DOH | 141$\pm$15 | 7.8$\pm$2.8 |  $\cdots$ | 160$\pm$39 | 1.8$\pm$1.2 |  $\cdots$ | 138$\pm$17 | 5.0$\pm$2.0eeFrom the CH2DOH and CH3OH column density ratio a deuteration fraction D/H = 0.05 is estimated, which is used for the calculation of the molecular abundances relative to methanol. | 5$\times$$10^{-2}$ | 2.7
13CH3OH |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | 0.8ffObtained from the observed spectra assuming $T_{\rm rot}$ = 140 K. | 8$\times$$10^{-3}$ | 0.44
CH318OH |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | 0.3ffObtained from the observed spectra assuming $T_{\rm rot}$ = 140 K. | 3$\times$$10^{-3}$ | 0.17
CH3OH |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | 100e,ge,gfootnotemark:  | 1 | 56

ggfootnotetext: The CH3OH lines are optically thick and saturated, and
therefore unsuitable for deriving the column density. The value listed in the
table is obtained as the geometrical mean of the CH3OH column densities of
0.6$\times$1019 cm-2 and 1.7$\times$1019 cm-2, derived from the 13CH3OH and
CH318OH isotopologues, assuming isotopic ratios of 12C/13C = 77 and 16O/18O =
560 (Wilson & Rood, 1994), respectively.

In Figure 15 we show the rotational diagrams toward the positions of VLA 4A,
VLA 4B, and toward the geometric center of the system (CB, that we adopt as
representative of the circumbinary disk). In Table 6 we report the rotational
temperatures and column densities derived from these diagrams. From CH3OCHO
and CH2DOH we obtain similar temperatures of $\sim$140 K toward both VLA 4A
and the central position CB. From ethylene glycol, aGg’-(CH2OH)2, a molecule
emitting in a small region of $\sim$$0\farcs 2$ ($\sim$60 au) in size around
VLA 4A, we find a slightly higher value, $T_{\rm rot}\simeq 145$ K, toward
this position. Therefore, we conclude that there are no substantial
differences in the rotational temperatures toward the CB and the VLA 4A
positions, measured at scales of the order of the beam size ($\sim$48 au
$\times$ 24 au), which is the size of the region where the spectra used in the
$T_{\rm rot}$ analysis have been obtained.

On the other hand, our two temperature measurements toward VLA 4B give values
of 90$\pm$6 K and 160$\pm$39 K. Since the observed molecular emission toward
this position could be more affected by dust opacity (see Sects. 3.4.1 and
3.4.2), we consider these values of $T_{\rm rot}$ toward VLA 4B as less
reliable. As noted in Section 3.2.2, the brightness temperature of VLA 4B
inferred from our highest angular resolution continuum image (Fig. 8 and Table
5) is $\sim$450 K, suggesting a possible temperature increase at distances
much closer to the star ($r$ $\lesssim$ 10 au), compared to the scales tested
in the $T_{\rm rot}$ analysis. In fact, we expect molecular emission to be
unobservable above the continuum toward the compact region of high optical
depth in VLA 4B, identified in the highest angular resolution continuum images
(baselines $>$750 k$\lambda$; Fig. 8), so we expect $T_{\rm rot}$ to be
insensitive to the region with the brightest and hottest circumstellar dust.
We will adopt 140 K as a representative value for the kinetic temperature, at
scales of the order of the beam size ($\sim$25-50 au), throughout the source
but with the caveat that in the circumstellar disks, at very small radii, dust
can reach higher temperatures.

Belloche et al. (2020b), from a large number of molecular transitions, but
with data at lower frequencies and poorer angular resolution ($\sim 0\farcs
33-1\farcs 75$) that does not separate the binary well, obtained a wide range
of rotational temperatures, from $\sim$30 K to $\sim$350 K. However, when
considering only the values with signal-to-noise ratio higher than 3.5$\sigma$
their temperatures range from $\sim$100 to $\sim$200 K, which are in better
agreement with ours. For the two molecules we have in common, CH3OCHO and
aGg’-(CH2OH)2, the temperatures obtained by Belloche et al. (2020b) are
$\sim$100 K, which are somewhat lower than our values.

The molecular column density of CH3OCHO obtained by these authors is
$2.1\times 10^{17}$ cm-2, which is a factor $\sim$2 lower than the value we
find toward VLA 4A, but similar to the column density we find in the rest of
the source (Table 6). For aGg’-(CH2OH)2, which we only detect toward the
position of VLA 4A, Belloche et al. (2020b) found a column density of
$1.3\times 10^{16}$ cm-2, which is about one order of magnitude lower than the
value we obtain. We attribute these discrepancies to the non-uniform
distribution of this molecular emission, as revealed by our higher angular
resolution data, which is smeared within the larger beam of the Belloche et
al. (2020b) observations, resulting in lower averaged values. We point out
that the rotational diagrams presented in Belloche et al. (2020b) have an
error in the labeling of the y-axis units (see Belloche et al. 2020a). Our
diagrams are consistent with those after the proper correction is performed.
The rotational temperatures obtained from these diagrams, and the column
densities reported by these authors, which were obtained from radiative
transfer models, are not affected by the error.

####  3.5.2 Disk Masses from Dust

Once we have an estimate of the temperature, we can derive the masses of the
circumstellar and circumbinary disks of SVS 13 from the observed 0.9 mm
continuum emission, which traces thermal dust emission (since the expected
free-free contribution at this short wavelength is negligible). Given the
uncertainties involved and the fact that we do not observe large variations of
rotational temperature from one position to another across the circumbinary
structure (Sect. 3.5.1), we will assume a single value of 140 K for the dust
temperature in the circumbinary disk.

For the circumstellar disks, the highest angular resolution 0.9 mm continuum
results discussed in Section 3.2.2 suggest dust temperatures of $\sim$300 K
for VLA 4A and of $\sim$450 K for VLA 4B, at scales $\lesssim$10 au (Table 5).
The molecular data are consistent with rotational temperatures of $\sim$140 K
toward the positions of both VLA 4A and VLA 4B, at beam scales of 25-50 au.
So, in the calculation of dust masses we will consider a range of dust
temperatures from 140 K to 300 K (VLA 4A) or from 140 K to 450 K (VLA 4B).

In addition, as noted in Section 3.2.2, the observed 0.9 mm dust brightness
temperatures indicate that this dust emission is optically thin in the
circumbinary disk, but not in the circumstellar disks, particularly in VLA 4B
where the dust emission can be very optically thick.

For optically thin emission from an isothermal source, the observed flux
density is directly related to the total mass of emitting material as:

$M=\frac{S_{\nu}D^{2}}{\kappa_{\nu}B_{\nu}(T_{\rm dust})},$ (3)

where $D$ is the distance to the source, $S_{\nu}$ is the flux density,
$\kappa_{\nu}$ is the dust opacity at the observed frequency $\nu$,
$B_{\nu}(T_{\rm dust})$ is the Planck function, and $T_{\rm dust}$ is the dust
temperature. Therefore, we can use eq. (3) to estimate the mass of
circumbinary material from the observed flux density at 0.9 mm.

If the dust emission is not optically thin, as we suspect is the case for the
circumstellar disks, then eq. (3) gives a lower limit of the mass. A suitable
upper limit to the circumstellar masses can be obtained from the 7 mm
emission, which is expected to be optically thin and dominated by the dust
emission (see Sect. 3.1.2). We take this as an upper limit because of a
possible residual free-free contamination. Therefore, using the 0.9 mm and 7
mm data we can constrain the masses of the circumstellar disks.

We assume a gas-to-dust ratio of 100 (Bohlin et al., 1978) and adopt a dust
opacity $\kappa_{\nu}$(0.9 mm) = 1.65 cm2 g-1 (per unit mass of dust). This
value of $\kappa_{\nu}$(0.9 mm) was obtained by interpolating the values at
0.7 mm and 1 mm from Ossenkopf & Henning (1994) and should be adequate for
protostellar objects where the dust grain growing is still incipient. These
dust opacities have been obtained by including thin ice mantles in an MRN
model (Mathis et al., 1977), consisting of a mixture of spherical silicate and
graphite grains with separated power-law particle size distributions.

Since the longest wavelength in the Ossenkopf & Henning (1994) opacities is
1.3 mm, we need to extrapolate the opacities up to 7 mm. For this, we assume a
power-law dependence of the dust opacity with frequency
($\kappa_{\nu}\propto\nu^{\beta}$, with 0$<$$\beta$$<$2), where $\beta$ is
related to the dust grain-size distribution and composition. In the case of
optically thin dust emission, $\beta$ can be inferred from the observed SED as
$\beta=\alpha-2$, where $\alpha$ is the observed spectral index
($S_{\nu}\propto\nu^{\alpha}$). Assuming that the overall emission of the
circumbinary dust of SVS 13 (green diamonds in Fig. 10) is optically thin, we
can use the spectral index obtained from this SED ($\alpha\simeq 2.5$) to
infer a value of $\beta$ = 0.5 for the circumbinary dust. By further assuming
that $\beta$ has the same value in the circumstellar and circumbinary
material, we can use $\beta$ = 0.5 to extrapolate the Ossenkopf & Henning
(1994) opacity $\kappa_{\nu}$(1.3 mm) = 0.899 cm2 g-1 to $\kappa_{\nu}$(7 mm)
= 0.387 cm2 g-1 (per unit mass of dust).

Using eq. (3), a dust temperature of $\sim$140 K (from the rotational diagrams
presented in Section 3.5.1), and a flux density at 0.9 mm of 909 mJy, obtained
by subtracting the estimated flux densities of VLA 4A and VLA 4B from the
total flux density given in Table 3.1.1, we obtain a total (dust+gas) mass for
the circumbinary disk $M_{\rm disk}$(CB) $\simeq$ 0.052 $M_{\odot}$.

For the circumstellar disks, assuming a range of dust temperatures 140-300 K
(VLA 4A) and 140-450 K (VLA 4B), we obtain strict lower limits of the disk
masses of $M_{\rm disk}$(VLA 4A) $\simeq$ 0.002-0.005 $M_{\odot}$ and $M_{\rm
disk}$(VLA 4B) $\simeq$ 0.002-0.006 $M_{\odot}$, from the likely optically
thick 0.9 mm emission. From the optically thin 7 mm emission, we obtain strict
upper limits of $M_{\rm disk}$(VLA 4A) $\simeq$ 0.004-0.009 $M_{\odot}$ and
$M_{\rm disk}$(VLA 4B) $\simeq$ 0.009-0.030 $M_{\odot}$. The masses obtained
from the 7 mm emission are higher by a factor of $\sim$2 for VLA 4A (partially
overlapping with the 0.9 mm mass range) and $\sim$5 for VLA 4B, strongly
suggesting very optically thick emission in this last case. Given that the
very high optical depth of the 0.9 mm emission would likely affect the mass
calculations more than the possible free-free contamination at 7 mm (which is
expected to be small, as discussed in Sect. 3.1.2), in the following we will
adopt the range of circumstellar disk masses obtained from the 7 mm data,
$\sim$0.004-0.009 $M_{\odot}$ for VLA 4A and $\sim$0.009-0.030 $M_{\odot}$ for
VLA 4B, as representative of the true values.

####  3.5.3 Dynamical Masses of the Stars

A powerful tool to obtain the stellar masses is from the measured 3D
velocities. For a gravitationally bound system of two objects in orbital
motion, the orbital period is given by:

$\left(\frac{P}{\rm yr}\right)=\left(\frac{a}{\rm
au}\right)^{3/2}\left(\frac{M_{\rm tot}}{M_{\odot}}\right)^{-1/2},$ (4)

where $M_{\rm tot}$ is the total mass of the system and $a$ is the separation
between the two components.

Equivalently,

$\left(\frac{M_{\rm tot}}{M_{\odot}}\right)=1.12\times
10^{-3}\left(\frac{a}{\rm au}\right)\left(\frac{V_{\rm rel}}{\rm
km\>s^{-1}}\right)^{2},$ (5)

where $V_{\rm rel}$ is the relative velocity between the two objects.

From the observed proper motions in right ascension and declination (Sect.
3.1.3) and the line-of-sight velocities inferred from the molecular lines
(Sect. 3.4.2), we can obtain the 3D relative velocity between the VLA 4A and
VLA 4B protostars, which is given in Table 4. This velocity is $\vec{V}_{\rm
rel}$ = (+0.48$\pm$1.02 km s-1, $-$2.28$\pm$0.82 km s-1, +1.97$\pm$0.28 km
s-1).

By adding quadratically the three velocity components we obtain $V_{\rm rel}$
= 3.06$\pm$0.65 km s-1. Then, taking into account that the separation is $a$ =
90 au (actually, we measure the projected separation but, since the two stars
are aligned in a direction close to the observed major axis of the
circumbinary structure, we do not expect an important difference between the
real and the projected separation), from eq. (5) we obtain a total mass of the
system of 1.0$\pm$0.4 $M_{\odot}$, and, from eq. (4), an orbital period of 850
yr. Also, the inclination angle of the binary orbital plane with respect to
the plane of the sky can be calculated as the arctan of the ratio between the
line-of-sight (LSR) and the plane-of-the-sky (proper motion) velocity
components, resulting $i$ = $41.7^{\circ}\pm 9.5^{\circ}$ (see Appendix C for
the uncertainty calculation).

We note that the total mass of the stars is much higher than those of the
circumstellar and circumbinary disks estimated from the dust emission in the
previous Section 3.5.2 (the total mass of the disks is $\sim$6-9% of the total
stellar mass). Therefore, it is reasonable to neglect these disks in the
calculation of the masses of the stars from the orbital parameters.

In principle, one could obtain a dynamical estimate of the total mass
(stars+disks) from the observed gradients in the first-order moment maps (Fig.
11) of the circumbinary molecular emission. However, the likely high optical
depth of the brighter molecular transitions, the complex overall structure and
kinematics of the source, with two central objects, spiral arms, and other
asymmetries prevent us from attempting this calculation. Nevertheless, we
still can obtain interesting results for the individual objects, if we
restrict our calculation to some specific molecular lines and the regions that
are less perturbed by asymmetries.

First, we discuss the implications of the kinematics of the four ethylene
glycol transitions that are only detected toward VLA 4A (Figs. 13 and 14). The
first-order moments and position-velocity diagrams in Figure 13 show the
presence of a clear velocity gradient in the E-W direction and centered on VLA
4A, which is consistently detected in all these transitions (note that the
rightmost panels include two blended transitions). This gradient is suggestive
of rotation motions associated with VLA 4A, with a rotation axis projected on
the plane of the sky at PA $\simeq$ 0∘. This is consistent with all rotation
axes having a similar position angle, PA $\simeq$ 0∘, as suggested by the fact
that the major axes of dust and molecular disks (observed in projection on the
plane of the sky) all have a PA $\simeq$ 90∘, as inferred in Sections 3.2.2,
3.4.1, and 3.4.3. As can be seen in Figures 13 and 14, toward VLA 4A there are
two sharply differentiated velocity components that we interpret may be
tracing the receding (eastern) and approaching (western) sides of a
circumstellar rotating disk around VLA 4A.

Thus, assuming a rotating disk around the VLA 4A protostar, from Figures 13
and 14 we infer a line-of-sight projected rotation velocity of $\sim$1.5 km
s-1 at a radius of $\sim$15 au, resulting in a mass $M_{*}$ = 0.04
$\sin^{-2}i~{}M_{\odot}$ for VLA 4A. Unfortunately, the value of this mass
depends strongly on the inclination angle, which is not known. Assuming
coplanarity between the circumstellar disk and the binary orbit ($i$ $\simeq$
40∘) we obtain $M_{*}$ = 0.09 $M_{\odot}$, but this value of the mass
increases significantly for smaller inclination angles (e.g., $M_{*}$ = 0.2
$M_{\odot}$ for $i$ = 25∘).

Alternatively, since the characteristic shape of a Keplerian velocity profile
is not clearly seen with the angular resolution of our observations, we
considered the possibility that the two velocity components of the ethylene
glycol transitions observed toward VLA 4A (Figs. 13 and 14) are tracing two
different objects orbiting one around each other. In this case, VLA 4A itself
would be a binary with a difference of $\sim$3 km s-1 in the line-of-sight
velocity of the two subcomponents and a projected separation $a\simeq 30$ au
(Fig. 13). Using eq. (5) and assuming an orbit coplanar with that of the main
SVS 13 binary ($i=40^{\circ}$), we obtain a total mass of $\sim$0.73
$M_{\odot}$ for VLA 4A. Furthermore, the presence of dust toward VLA 4A and
the relatively wide velocity range ($\sim$2.5 km s-1) of each subcomponent in
the position-velocity diagrams suggest that each of them could be associated,
in its turn, with an angularly unresolved ($\lesssim 0\farcs 1$ in size) disk,
implying an individual mass of each subcomponent of $\lesssim$0.064
$M_{\odot}$, if all the orbits were coplanar. Thus, the total mass of VLA 4A
obtained from the dynamics of the disk of each subcomponent separately
($\lesssim$0.13 $M_{\odot}$) would be significantly smaller than the mass
inferred from their overall orbital motion ($\sim$0.73 $M_{\odot}$). A
misalignment of $\sim$25∘ between this orbital plane and the circumstellar
disks of the potential VLA 4A binary would be required to reach an agreement
in the total mass, making this binary scenario for VLA 4A somewhat
speculative.

In order to obtain more accurate values of the mass and an estimate of
additional parameters such as the inclination angle, in Appendix D we perform
a fit of the observed channel maps of the ethylene glycol transitions to a
Keplerian thin disk model (see Zapata et al. 2019 for further details). The
best-fitted physical and geometrical parameters of the Keplerian disk
surrounding VLA 4A are presented in Tables 11 and 12. The best fit is obtained
for a disk model with an inner radius of $\sim$2 au, an outer radius of
33$\pm$2 au, an inclination of 22∘$\pm$7∘, and a Keplerian rotation velocity
of 2.7 km s-1 at a reference radius of 30 au, implying a central mass of
0.25$\pm$0.15 $M_{\odot}$. It is remarkable that the three fits to different
ethylene glycol transitions (one of them consisting of two blended lines)
consistently give very similar results (see Tables 11 and 12).

We note, however, that the disk inclination inferred from the gas rotation in
VLA 4A is significantly smaller than the inclination inferred from the aspect
ratio of dust continuum sizes obtained from a Gaussian fit (Sect. 3.2.2). We
consider the former to be more reliable, as it is obtained from a model
fitting over several spectral channels, but it is still very uncertain. Since
the beams are elongated in the direction of the minor axes of the disks, the
actual uncertainties on these axes, and on the inclinations derived from them,
are larger than the formal fitting errors for both the dust and the line
images.

Figure 16: Observed position-velocity diagrams of the continuum-subtracted
CS(7-6) line emission. (a) Along the E-W direction (PA = 90∘) toward the
central position (CB) midway between VLA 4A and VLA 4B. Offsets are in
increasing right ascension relative to this central position. (b) Along a N-S
direction (PA = 0∘), toward VLA 4A, with offsets in increasing declination
relative to this position. (c) Same as (b) but for VLA 4B. The dashed line
indicates the systemic velocity of the cloud (8.5 km s-1). The solid curves,
which are plotted as a reference, correspond to the line-of-sight velocity of
a particle in free-fall toward a point-like central mass of 0.1, 0.3, and 0.5
$M_{\odot}$, as a function of the projected separation.

Unfortunately a similar analysis in terms of disk rotation cannot be performed
in VLA 4B since there are no molecular transitions that trace its potential
circumstellar disk, likely because of the high dust opacity at these very
small scales ($\sim$10-20 au). However, there are signs of gravitational
acceleration at scales somewhat larger ($\sim$100-200 au), traced by the
molecular emission of the infalling material onto the disk. The signature of
this gravitational acceleration, with the highest observed velocity along a
given line-of-sight increasing as the (projected) distance to the central
protostar decreases, can be identified in the CS(7-6) position-velocity plots
shown in Figure 16 (the free-fall velocity of a particle for different values
of the central mass is plotted for reference). This behavior can be more
clearly identified on the western side of VLA 4A and on the eastern side of
VLA 4B, which appear less perturbed by the presence of the respective
companion. In addition, the presence of a jet originating in VLA 4B and
pointing to the southeast (see Fig. 2) has apparently swept away the material
in this direction, which appears devoid of molecular emission; therefore, for
the study of VLA 4B we have considered only the northeastern part of the PV
diagram.

In order to obtain an estimate of the masses of the protostars VLA 4A and VLA
4B, we calculated the CS (7-6) line intensity resulting from a simple model
consisting of a spherical envelope in gravitational collapse onto a central
compact object with velocity and temperature power-law dependences with the
radius with indices of $-0.5$. We are aware that such a simple model will
provide only a crude estimate of the physical parameters, as the lack of
spherical symmetry is evident. However since we restrict the fitting to the
less perturbed parts of the source and to a relatively narrow range of radii,
we consider that this local approximation can give satisfactory results as a
local measurement of the velocity and gravitational acceleration is enough to
measure the enclosed mass. The parameter space was explored by comparing the
results of the model, after convolution with the observing beam, with the
observations. A range of values of the free parameters (central mass, inner
and outer radii of the envelope) for which the model agrees with the observed
position-velocity diagrams were obtained for VLA 4A and VLA 4B (see details in
Appendix E).

The parameters of the best-fit models for VLA 4A and VLA 4B are given in
Tables 13 and 14, respectively. In Figures 26, 27, and 28 we show a comparison
of the model and observed PV plots. From the best-fit models we obtain for VLA
4A an infalling envelope with an inner radius of $\sim$15 au, an outer radius
of $\sim$105 au, and a central mass of 0.27$\pm$0.10 $M_{\odot}$, where the
quoted uncertainty is only the formal fitting error and does not include the
effect of possible non-sphericity and optical depth. For VLA 4B we obtain an
envelope inner radius of $\sim$12 au, an outer radius of $\sim$150 au, and a
central mass of 0.6$\pm$0.2 $M_{\odot}$. Interestingly, these results are
fully consistent with the dynamical masses obtained above from the ethylene
glycol transitions and from the orbital motion of the SVS 13 system. Also, we
note that the inner radii of the envelopes nicely fit the outer radii of the
circumstellar disks, so we consider these results as satisfactory.

Finally, we can assess the importance of self-gravity in the SVS 13 disks
based on the analysis of the Toomre Q parameter. Theoretically, $Q>1$
indicates that the self-gravity of the disk is negligible compared with the
gravity of the central object. Such a disk would be stable against local
gravitational instabilities and would not be expected to undergo
fragmentation. For Keplerian disks $Q\simeq 2M_{*}H/(M_{\rm disk}R)$ (Kratter
& Lodato, 2016; Tobin et al., 2016a, 2020), where $M_{*}$ and $M_{\rm disk}$
are the masses of the central protostar/s and the disk, respectively, and $H$
is the scale height of the disk at radius $R$. For a Keplerian disk, $H$ =
$c_{s}/\Omega_{K}$, where $c_{s}$ is the sound speed in the gas (in this case
H2) and $\Omega_{K}$ is the Keplerian angular velocity at radius $R$. Using
the values derived from our observations (see also Sections 3.2.2 and 3.5.2),
we obtain for the circumbinary disk, $c_{s}$ = 0.9 km s-1 (for $T$ = 140 K),
$H$ = 40 au (at the disk outer radius of 120 au), and $Q$ = 13. For the VLA 4A