# The ALMA-PILS survey: 3D modeling of the envelope, disks and dust filament of IRAS 16293-2422

S. Jacobsen, J. Jørgensen, H. van Der Wiel, H. Calcutt, L. Bourke, C. Brinch, Audrey Coutens, M. Drozdovskaya, L. Kristensen, S. P. Müller, et al.

► To cite this version:

S. Jacobsen, J. Jørgensen, H. van Der Wiel, H. Calcutt, L. Bourke, et al. The ALMA-PILS survey: 3D modeling of the envelope, disks and dust filament of IRAS 16293-2422. Astronomy and Astrophysics - A&A, 2018, 612, pp.id.A72. 10.1051/0004-6361/201731668. hal-01797148

HAL Id: hal-01797148

https://hal.science/hal-01797148

Submitted on 1 Dec 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The ALMA-PILS survey: 3D modeling of the envelope, disks and dust filament of IRAS 16293-2422

S. K. Jacobsen<fnref n="1" />, J. K. Jørgensen<fnref n="1" />, M. H. D. van der Wiel<fnref n="1" />,<fnref n="2" />, H. Calcutt<fnref n="1" />, T. L. Bourke<fnref n="3" />, C. Brinch<fnref n="4" />, A. Coutens<fnref n="5" />,<fnref n="6" />, M. N. Drozdovskaya<fnref n="7" />, L. E. Kristensen<fnref n="1" />, H. S. P. Müller<fnref n="8" />, and S. F. Wampfler<fnref n="7" />

Received 28 July 2017 / Accepted 2 December 2017

## ABSTRACT

**Context.** The Class 0 protostellar binary IRAS 16293-2422 is an interesting target for (sub)millimeter observations due to, both, the rich chemistry toward the two main components of the binary and its complex morphology. Its proximity to Earth allows the study of its physical and chemical structure on solar system scales using high angular resolution observations. Such data reveal a complex morphology that cannot be accounted for in traditional, spherical 1D models of the envelope.

**Aims.** The purpose of this paper is to study the environment of the two components of the binary through 3D radiative transfer modeling and to compare with data from the Atacama Large Millimeter/submillimeter Array. Such comparisons can be used to constrain the protoplanetary disk structures, the luminosities of the two components of the binary and the chemistry of simple species.

**Methods.** We present $^{13}$CO, $^{17}$O and $^{18}$O J = 3-2 observations from the ALMA Protostellar Interferometric Line Survey (PILS), together with a qualitative study of the dust and gas density distribution of IRAS 16293-2422. A 3D dust and gas model including disks and a dust filament between the two protostars is constructed which qualitatively reproduces the dust continuum and gas line emission.

**Results.** Radiative transfer modeling in our sampled parameter space suggests that, while the disk around source A could not be constrained, the disk around source B has to be vertically extended. This puffed-up structure can be obtained with both a protoplanetary disk model with an unexpectedly high scale-height and with the density solution from an infalling, rotating collapse. Combined constraints on our 3D model, from observed dust continuum and CO isotopologue emission between the sources, corroborate that source A should be at least six times more luminous than source B. We also demonstrate that the volume of high-temperature regions where complex organic molecules arise is sensitive to whether or not the total luminosity is in a single radiation source or distributed into two sources, affecting the interpretation of earlier chemical modeling efforts of the IRAS 16293-2422 hot corino which used a single-source approximation.

**Conclusions.** Radiative transfer modeling of source A and B, with the density solution of an infalling, rotating collapse or a protoplanetary disk model, can match the constraints for the disk-like emission around source A and B from the observed dust continuum and CO isotopologue gas emission. If a protoplanetary disk model is used around source B, it has to have an unusually high scale-height in order to reach the dust continuum peak emission value, while fulfilling the other observational constraints. Our 3D model requires source A to be much more luminous than source B; $L_A \sim 18 L_\odot$ and $L_B \sim 3 L_\odot$.

**Key words.** astrochemistry – stars: formation – stars: protostars – ISM: individual objects: IRAS 16293-2422 (except planetary nebulae) – ISM: molecules – submillimeter: ISM

## 1. Introduction

Stars form inside the dense cores of interstellar cold clouds of gas and dust. During the gravitational collapse of the dense core material, the slow rotation of the core itself will eventually cause a protoplanetary disk to appear, due to the conservation of angular momentum (Cassen & Moosman 1981). When exactly such disks are formed and how they evolve through their early stages remain poorly understood (see, e.g., Li et al. 2014, for a review). Well established disks around T Tauri stars or Class II young

stellar objects have been widely reported and studied extensively (e.g., see Andrews et al. 2009; Ricci et al. 2010). Due to the low spatial resolution of the observations used in earlier studies, the emission from the disks in Class 0 and Class I sources have been mixed with the emission from the massive amounts of dust in the ambient protostellar envelope, making detection and characterization of disks in Class 0 sources very difficult (Jørgensen et al. 2005a, 2009; Chiang et al. 2008), with only a small number of Class 0 disks with indicated Keplerian rotation being detected so far (Tobin et al. 2012; Murillo et al. 2013; Ohashi et al. 2014;

<footnote n="1">¹ Centre for Star and Planet Formation, Niels Bohr Institute and Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen K, Denmark
e-mail: steffen.jacobsen@ni ku.dk</footnote>

<footnote n="2">² ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA Dwingeloo, The Netherlands</footnote>

<footnote n="3">³ SKA Organization, Jodrell Bank Observatory, Lower Withington, Macclesfield, Cheshire SK11 9DL, UK</footnote>

<footnote n="4">⁴ Niels Bohr International Academy, Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen Ø, Denmark</footnote>

<footnote n="5">⁵ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="6">⁶ Department of Physics and Astronomy, University College London, Gower St., London, WC1E 6BT, UK</footnote>

<footnote n="7">⁷ Center for Space and Habitability (CSH), Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland</footnote>

<footnote n="8">⁸ I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany</footnote>
Lindberg et al. 2014; Codella et al. 2014). Jørgensen et al. (2009) detect compact mm continuum emission around several Class 0 objects, but it is unknown whether these dust components represent rotationally supported Keplerian disks. The exact transition from early density structures to a rotationally supported Keplerian disk is still uncertain, though some tentative observations of the so-called centrifugal barrier, where the infalling, rotating envelope transitions into a rotationally supported disk, have been made (Sakai et al. 2014; Oya et al. 2016). These inner regions around embedded protostars are also particularly interesting from a chemical point of view. Observational studies of Class 0 objects have revealed complex organic molecules (e.g., Bottinelli et al. 2004b; Kuan et al. 2004; Bottinelli et al. 2004a; Jørgensen et al. 2005a; Bisschop et al. 2008; Jørgensen et al. 2011; Maury et al. 2014), as well as active chemistry in the so-called hot corino region, where water and organic molecules are sublimated into the gas-phase from the dust grain surface (Schöier et al. 2002).

A particularly interesting source in this context is IRAS 16293-2422 (hereafter IRAS 16293), a nearby binary protostar in the ρ Ophiuchus cloud complex, composed of at least two components, IRAS 16293A and IRAS 16293B (Wootten 1989; Mundy et al. 1992; Looney et al. 2000). IRAS 16293 is one of the most well-studied low-mass protostars, particularly in terms of astrochemistry; IRAS 16293 was the first low-mass star for which complex organics (van Dishoeck et al. 1995) and prebiotic molecules (Jørgensen et al. 2012) were detected. The source IRAS 16293A appears to be highly active with several detected outflows, while IRAS 16293B is more quiescent, with no known outflows (e.g., Chandler et al. 2005; Kristensen et al. 2013). Alves et al. (2012) find that the magnetic field strengths in IRAS 16293 are comparable to the outflow ram pressure, and thus that the magnetic field is dynamically important for this system. Both components are embedded within a larger, circumbinary envelope of dust and gas, which has been well studied in previous works using single-dish data and 1D radiative transfer modeling (Schöier et al. 2002; Jørgensen et al. 2005b; Crimier et al. 2010). These models all use a single central radiation source, where disks or complex structures, such as a dust and gas filament, have never been included in the radiative transfer calculations. The single-star approximation is not a problem at large distances from the interbinary region as the two sources can then be approximated as a single source, especially if one source is more luminous than the other. In contrast, in the interbinary region, two sources are needed in order to estimate the temperature structure. The total luminosity is $L = 21 \pm 5 L_\odot$ (Jørgensen et al. 2016), while the luminosity ratio of the two protostars is unknown. IRAS 16293 has been, and continues to be, under intense observational scrutiny due to its rich chemistry and any chemical modeling must assume the parameters of the physical environment in terms of dust and $H_2$ density, temperature structure and the local radiation field. An assessment of the individual luminosities of the sources in IRAS 16293 is crucial for chemical modeling, especially for warm gas-phase chemistry in the hot corinos.

In this paper, we construct a 3D model, including disks, a filament and two radiation sources using high angular resolution data from the ALMA Protostellar Interferometric Line Survey (PILS) as observational constraints. The impact of the binary nature of IRAS 16293 on the temperature distribution is evaluated as opposed to that of a single star. As stellar binaries are abundant, especially in the Class 0 stage (Tobin et al. 2016), this temperature evaluation of IRAS 16293 is useful as a case study

on the temperature distribution in the environment of protostellar binaries in wide orbits, found in large fractions of star-forming regions (e.g., Duchêne et al. 2004), as well as a retrospective analysis on the validity of previous works approximating a binary as a single protostar in their radiative transfer modeling, which will be useful for data-mining. This work is part of a new generation of complex 2D and 3D radiative transfer models describing Class 0 YSOs (see e.g., Yang et al. 2017), which were previously modeled with 1D envelope models.

After the observations and the general radiative transfer setup are described in Sects. 2 and 3, the paper is divided into three parts: Sect. 3.5 supplies a general investigation of temperature structures in Class 0 objects and the difference between having one or two central radiation sources in a 1D model of a collapsing envelope. Section 4 shows a 3D model of the dust density structure of IRAS 16293, where we used a dust continuum image at 868 µm and the spectral energy distribution (SED) as constraints. Finally, new ALMA-PILS observations of the CO isotopologues $^{13}$CO, $C^{17}$O and $C^{18}$O are used in Sect. 5 to further constrain the physical conditions in the innermost region of IRAS 16293, ending with an evaluation of the nature of the disks, their vertical distribution and the luminosity ratio between the sources, offering a significant increase in physical complexity from earlier models as well as the first luminosity evaluation of the individual sources.

## 2. Observations

The PILS program (PI: Jes K. Jørgensen, project-id: 2013.1.00278.S) is an unbiased spectral survey of IRAS 16293, covering a significant part of ALMA's Band 7, corresponding to a wavelength of roughly 0.8 mm, as well as selected windows of ALMA's Bands 3 and 6, corresponding to ~3 mm and ~1.3 mm, respectively. The Band 7 data cover the frequency range from 329.15 to 362.90 GHz with ~0.2 km s⁻¹ spectral resolution and ~0.5” angular resolution (60 AU beam diameter, assuming $d = 120$ pc). The pointing center was set between source A and B at $\alpha_{J2000} = 16^h32^m22^s72; \delta_{J2000} = -24^\circ28'34''.3$. Further details about the survey, data reduction and an overview of the observations can be found in Jørgensen et al. (2016).

In this paper, we focused on the 868 µm continuum and CO isotopologue data, which include data from both the array of 12 m telescopes and the Atacama Compact Array (ACA), allowing both compact and extended structures to be detected. The continuum image was taken as an average of the Band 7 continuum in each channel (Jørgensen et al. 2016). The 868 µm dust continuum shows broad extended emission, enveloping both sources. Around each source, emission structures consistent with disks are observed. Source A appears to host a disk close to edge-on, while the emission around source B appears more circular, that is, suggestive of a face-on disk. The continuum is much brighter toward source B, peaking at 2.0 Jy beam⁻¹, while source A has a more extended and weaker continuum, peaking at 1.0 Jy beam⁻¹ (Fig. 1). The dust continuum images from Bands 3 and 6 were not used in this work due to lower spatial resolutions than the Band 7 data, which confuses the disk structures with the filament.

The gas line emission from the $J = 3-2$ transitions of three CO isotopologues, $^{13}$CO (330.588 GHz, Klapper et al. 2000), $C^{18}$O (329.331 GHz, Klapper et al. 2001) and $C^{17}$O (337.061 GHz, Klapper et al. 2003) were extracted from the dataset, using the entire pixel map, after the lines were identified. The source A LSR velocity of 3.1 km s⁻¹
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Dust continuum at 868 µm, in log-scale. Contour levels are logarithmically divided between 0.02 and 2 Jy beam⁻<fnref n="1" />. The RA and Dec offsets are relative to the phase center of the observations. Beamsize is shown in bottom right corner.</figcaption></figure>

(Jørgensen et al. 2011) was used for all lines. Channels within ±5 km s⁻<fnref n="1" /> of each line were combined into an integrated intensity (moment zero) map for each isotopologue, using the image analysis software MIRIAD. Using this approach, the vast majority of the observed line emission was extracted, while avoiding major line-blending from other species.

Figure 2 shows the integrated intensity maps of the three isotopologues. C¹⁷O appears to trace the dust distribution fairly well, as does the C¹⁸O emission. ¹³CO, however, shows a markedly different emission structure, which is likely attributable to optical depth effects, where only the outer layer of the ¹³CO structure is probed. Absorption toward source B is evident in all lines but to a smaller extent in C¹⁷O, revealing that C¹⁷O is the line with the lowest optical depth, as expected since it has the lowest abundance (Wilson & Rood 1994). An interesting feature is the broad emission arc between the sources, which is also present in the dust continuum emission, suggesting a common origin of both the dust and gas emission in this region. We focused on the overall structure of the emission and left a more detailed analysis of the gas kinematics to a future paper (Van der Wiel et al., in prep.).

### 3. Dust radiative transfer modeling

This section describes the setup of the radiative transfer code RADMC-3D, including a description of the dust opacity, model grid and the envelope dust density model used. At the end of the Section, we present our investigation of the temperature structure difference between a single- and a binary protostar, where the envelope dust density model was used. We constructed 3D dust and gas line radiative transfer models consisting of a large envelope with a radius of 8 × 10³ AU based on Schöier et al. (2002), two disks around the protostars, as well as a dust filament between the two protostars. The disk's dust densities were modeled with both the model protoplanetary disk (PP-disk) structure known from Class II sources (e.g., Andrews et al. 2012) and the density solution from the angular momentum conservation of an infalling, rotating collapse (Terebey et al. 1984). This infalling, rotating collapse density structure is hereafter

referred to as a rotating toroid model, based on the toroidal shape appearing at the centrifugal radius (Terebey et al. 1984), which is the radius where the infalling material with the largest angular momentum impacts the midplane. In terms of luminosity, all models, whether they used one or two radiation sources, were limited to a total luminosity of 21 L⊙. This luminosity is based on integration of the SED using data from Spitzer, ISO-LWS, Herschel/SPIRE and submillimeter data (Chandler et al. 2005; Correia et al. 2004; Makiwa 2014; Schöier et al. 2002), where the luminosity is estimated to be 21 ± 5 L⊙ (Jørgensen et al. 2016). This luminosity assumes a distance of 120 pc, based on recent distance estimates of the bulk of the ρ Oph cloud from extinction measurements (Knude & Hog 1998; Lombardi et al. 2008) and VLBI parallax (Loinard et al. 2008). The interbinary distance was fixed to 636 AU for all models, based on the angular distance of ~5.3" between the peaks of the 868 µm ALMA image, with the caveat that this distance, in reality, is only the projected, that is, the minimum distance between the sources.

#### 3.1. RADMC-3D

To determine the dust temperature and to produce continuum images, the radiative transfer code RADMC-3D (Dullemond et al. 2012, v. 0.40)<fnref n="1" /> was used. RADMC-3D uses the Bjorkman & Wood (2001) algorithm to perform a Monte Carlo simulation of photon propagation from one or more radiation sources, through a predefined grid of dust density cells, in order to estimate the dust equilibrium temperature. This equilibrium is a balance between the incoming radiative energy and the blackbody emission from the dust cell itself. The code setup in this work assumed that all the dust energy originated from the protostars through irradiation, and not from any internal energy source such as viscous accretion, or any external radiation source. The precision of the resulting dust equilibrium temperature depends, among other parameters, on the number of photons used in the radiation simulation (this differed between model setups, see Table 1). We ignored scattering, since the absorption opacity dominates over the scattering opacity in the submillimeter regime (Dunham et al. 2010). RADMC-3D was also used to calculate the SED (see Appendix A) and to produce synthetic continuum images (see Sect. 4.3).

#### 3.2. Grid

We constructed a spherical grid of cells for use in the RADMC-3D radiative transfer calculations, resolved linearly in the polar and azimuthal angle ranges, while the radial dimension was divided into an inner and outer region. The inner region with constant dust density (i.e., $r \le r_{\text{plat}}$, see Eq. (1)) was resolved linearly, while $r > r_{\text{plat}}$ (containing the bulk mass of the cloud) was resolved logarithmically. In the disk components, of the global model, octree refinement was performed (splitting a grid cell into eight new grid cells), with two refinement levels, to ensure that the density gradients of the disk component models were resolved, see Table 1. If a PP-disk model was used, then the inner 10 AU of the disk was refined two times further, for a total of four octree refinements, to resolve the inner gap and nearby structure properly, as geometric effects such as self-shadowing onto the outer disk became important.

<footnote n="1">¹ http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/</footnote>
<table><caption>**Table 1.** RADMC-3D model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td colspan="3"><strong>Grid</strong></td></tr><tr><td>r<sub>out</sub></td><td>Envelope radius</td><td>8 × 10<sup>3</sup> AU</td></tr><tr><td>r<sub>in</sub></td><td>Grid start radius</td><td>5 AU</td></tr><tr><td>r<sub>plat</sub></td><td>Radius of the density plateau</td><td>600 AU</td></tr><tr><td>n<sub>θ</sub></td><td>Number of grid cells in the polar range</td><td>131</td></tr><tr><td>n<sub>φ</sub></td><td>Number of grid cells in the azimuthal range</td><td>131</td></tr><tr><td colspan="3"><strong>Temperature analysis</strong></td></tr><tr><td>n<sub>r, out</sub></td><td>Number of grid cells in the outer radial region (r > 600 AU)</td><td>88</td></tr><tr><td>n<sub>r, in</sub></td><td>Number of grid cells in the inner radial region (r ≤ 600 AU)</td><td>110</td></tr><tr><td>n<sub>photons</sub></td><td>Number of photons used in thermal Monte Carlo process</td><td>5 × 10<sup>7</sup></td></tr><tr><td>L<sub>A</sub></td><td>Luminosity of source A</td><td>3–20 L<sub>☉</sub></td></tr><tr><td colspan="3"><strong>One radiation source</strong></td></tr><tr><td>r<sub>oct</sub></td><td>Octree refinement radius around source</td><td>400 AU</td></tr><tr><td>n<sub>octree levels</sub></td><td>Octree refinement number</td><td>1</td></tr><tr><td colspan="3"><strong>Two radiation sources</strong></td></tr><tr><td>r<sub>oct</sub></td><td>Octree refinement radius around each source</td><td>250 AU</td></tr><tr><td>n<sub>octree levels</sub></td><td>Octree refinement number</td><td>2</td></tr><tr><td colspan="3"><strong>IRAS 16293 model</strong></td></tr><tr><td>n<sub>r, out</sub></td><td>Grid cells in the outer radial region (r > 600 AU)</td><td>22</td></tr><tr><td>n<sub>r, in</sub></td><td>Grid cells in the inner radial region (r ≤ 600 AU)</td><td>110</td></tr><tr><td>ρ<sub>0, env</sub></td><td>Envelope reference density at r<sub>0, env</sub></td><td>2.5 × 10<sup>-14</sup> g cm<sup>-3</sup></td></tr><tr><td>r<sub>0, env</sub></td><td>Envelope reference radius</td><td>1 AU</td></tr><tr><td>p<sub>env</sub></td><td>Envelope density power-law exponent</td><td>1.7</td></tr><tr><td>M<sub>A</sub></td><td>Source A stellar mass</td><td>1 M<sub>☉</sub></td></tr><tr><td>M<sub>B</sub></td><td>Source B stellar mass</td><td>0.1 M<sub>☉</sub></td></tr><tr><td>†T<sub>eff</sub></td><td>Star surface temperature</td><td>5000 K</td></tr><tr><td>L<sub>A</sub></td><td>Source A luminosity</td><td>3–20 L<sub>☉</sub></td></tr><tr><td>r<sub>0, disk</sub></td><td>PP-disk reference radius</td><td>10 AU</td></tr><tr><td>r<sub>inner disk</sub></td><td>Disk inner radius</td><td>1 AU</td></tr><tr><td>r<sub>disk, A</sub></td><td>Disk A radius</td><td>150 AU</td></tr><tr><td>r<sub>disk, B</sub></td><td>Disk B radius</td><td>50 AU</td></tr><tr><td colspan="3"><strong>Disk constraints</strong></td></tr><tr><td>F<sub>peak, A</sub></td><td>Peak flux density toward source A</td><td>1.0 Jy beam<sup>-1</sup></td></tr><tr><td>F<sub>peak, B</sub></td><td>Peak flux density toward source B</td><td>2.0 Jy beam<sup>-1</sup></td></tr><tr><td>Disk A aspect ratio</td><td>...</td><td>1.64</td></tr><tr><td>N<sub>H<sub>2</sub>, B</sub></td><td>Disk B column density</td><td>≥1.2 × 10<sup>25</sup> cm<sup>-2</sup></td></tr><tr><td colspan="3"><strong>PP-disk models</strong></td></tr><tr><td>Σ<sub>A</sub></td><td>Parameter space for reference surface density of PP-disk A at 10 AU</td><td>0.1–5.0 g cm<sup>-2</sup></td></tr><tr><td>Σ<sub>B</sub></td><td>Parameter space for reference surface density of PP-disk B at 10 AU</td><td>1.0–12.0 g cm<sup>-2</sup></td></tr><tr><td>Ψ</td><td>Parameter space of the PP-disks flaring constants</td><td>0–0.25</td></tr><tr><td>p</td><td>Parameter space of the PP-disks dust density radial power-law exponents</td><td>0.5–1.5</td></tr><tr><td>H<sub>0, A</sub></td><td>Parameter space of PP-disk A scale-height at 10 AU</td><td>0.85–2.0 AU</td></tr><tr><td>H<sub>0, B</sub></td><td>Parameter space of PP-disk B scale-height at 10 AU</td><td>0.85–3.0 AU</td></tr><tr><td>i<sub>A</sub></td><td>Parameter space of PP-disk A inclination, 90° is edge-on</td><td>30–90°</td></tr><tr><td>n<sub>lev, disk</sub></td><td>Octree refinement level of entire disk</td><td>2</td></tr><tr><td>n<sub>lev, in disk</sub></td><td>Octree refinement level in innermost region</td><td>4</td></tr><tr><td>n<sub>photons</sub></td><td>Number of photons used in thermal Monte Carlo process</td><td>10<sup>7</sup></td></tr><tr><td colspan="3"><strong>Rotating toroid models</strong></td></tr><tr><td>r<sub>c, A</sub></td><td>Parameter space of the centrifugal radius of rotating toroid A</td><td>50–100 AU</td></tr><tr><td>r<sub>c, B</sub></td><td>Parameter space of the centrifugal radius of rotating toroid B</td><td>5–40 AU</td></tr><tr><td>M&#775;<sub>A</sub></td><td>Parameter space of the mass accretion rate into source A</td><td>1.0 × 10<sup>-6</sup>–5.5 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>M&#775;<sub>B</sub></td><td>Parameter space of the mass accretion rate into source B</td><td>1.0 × 10<sup>-6</sup>–5.5 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>n<sub>lev, disk</sub></td><td>Octree refinement level of entire disk</td><td>2</td></tr><tr><td>n<sub>lev, in disk</sub></td><td>Octree refinement level in innermost region</td><td>2</td></tr><tr><td>n<sub>photons</sub></td><td>Number of photons used in thermal Monte Carlo process</td><td>10<sup>6</sup></td></tr><tr><td colspan="3"><strong>Dust arc</strong></td></tr><tr><td>ρ<sub>0</sub></td><td>Reference density</td><td>2.5 × 10<sup>-17</sup> g cm<sup>-3</sup></td></tr><tr><td>r<sub>out</sub></td><td>Filament radius from filament center</td><td>130 AU</td></tr><tr><td>Filament semi-minor axis</td><td>...</td><td>170 AU</td></tr><tr><td>Filament semi-major axis</td><td>...</td><td>300 AU</td></tr></tbody></table>

Notes. Model parameters for each setup. L<sub>B</sub> is always defined as L<sub>B</sub> = 21.0 L<sub>☉</sub>–L<sub>A</sub>. (†) The radiation spectra of both source A and B are treated as blackbodies with a temperature of 5000 K, following Schöier et al. (2002).
<figure><img src="image_3.png" /><figcaption>Fig. 2. Zeroth moment maps of the ¹³CO, C¹⁸O, and C¹⁷O isotopologues from the PILS observations. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam⁻¹ km s⁻¹. The RA and Dec offsets are relative to the phase center of the observations.</figcaption></figure>

### 3.3. Dust opacity and temperature

We used dust opacities from Ossenkopf & Henning (1994), which are good approximations for dense cores (van der Tak et al. 1999; Shirley et al. 2002, 2011). Opacity tables of dust with thin ice-mantles and another with bare dust grains were used, assuming coagulated dust grains with an ambient gas number density of 10⁶ cm⁻³. These dust opacity tables are generally used for models of dust around deeply embedded protostars, including IRAS 16293 (e.g., Schöier et al. 2002; Jørgensen et al. 2005b; Crimier et al. 2010). The dust opacity was allocated in the cells in a self-consistent manner, with the first RADMC-3D thermal Monte Carlo photon propagation done in a model of purely icy-dust opacities. Grid cells with temperatures above the assumed water sublimation temperature were then allocated bare-grain dust opacities for the next photon propagation. We did not change the sublimation temperature to account for the local environment in terms of pressure but instead used a single sublimation temperature of 90 K (Sandford & Allamandola 1993). As the cell opacities were redefined, the cell dust equilibrium temperatures after the next photon propagation will be different, with some dust cells shifting from below to above 90 K and vice versa, that is, an opacity allocation error. For both disk component model types, this error was found to be negligible after two recursions (5 and 0.1% relative to the total bare-grain dust mass, at 1σ, for the PP-disk model and the rotating toroid model, respectively). The higher error in the PP-disk model, even with ten times more photon packages than the rotating toroid model (Table 1), came from the very high dust densities necessary to match the observed peak flux density, especially in disk B. In contrast, the relatively more uniform dust density distribution in a rotating toroid model ensured good photon statistics in RADMC-3D, with a much lower photon number.

Therefore, two recursions of the bare-grain opacity allocation, that is, three Monte Carlo dust temperature estimations in total were performed in RADMC-3D, for a single parameter set. The dust temperature estimated in the final photon propagation was then used for all later computations with that parameter set.

### 3.4. Dust envelope model

The dust envelope was modeled with a 1D radial density power-law density distribution, as the envelope has been well-fitted previously by such models (e.g., Schöier et al. 2002). More complicated circumbinary envelope structures, such as an infalling,

rotating circumbinary envelope were not considered, as we focused on the innermost region of IRAS 16293 in this work. Recent observations of GG Tau (Dutrey et al. 2014), as well as simulations of the same system (Nelson & Marzari 2016) suggest that binary star formation results in partial clearing of the material in the innermost region where the protostellar binary resides. Based on the 868 µm continuum observation, we do not expect the interior to be completely evacuated of dust between the disks, but depleted instead, due to accretion into the two protostars and dynamical effects, as described by Nelson & Marzari (2016). A dust density plateau was used to imitate the mass depletion inside a given radius $r_{plat}$, meaning that we used a constant, rather than a radial power-law, density distribution in this inner region, while keeping the dust density distribution in the outer envelope as a power-law:

$$ \rho(r) = \begin{cases} \rho_0 \left(\frac{r}{r_0}\right)^{-p}, & \text{if } r \ge r_{plat} \\ \rho_0 \left(\frac{r_{plat}}{r_0}\right)^{-p}, & \text{if } r < r_{plat}. \end{cases} \qquad (1) $$

Where $\rho$ is dust density, $r$ is the distance to the model center, $\rho_0$ is the density at distance $r_0$ and $r_{plat}$ is the radius of the density plateau. $r_{plat}$ was fixed to 600 AU in all models, a distance we chose to ensure that the entire dust filament structure resided in this density plateau. We fixed $p = 1.7$, taken from Schöier et al. (2002).

### 3.5. Envelope temperature structure

In order to compare the temperature structure between a single star and a binary system, we used the dust envelope model from Eq. (1) with one or two radiation sources in RADMC-3D, without including the filament and disk component models. We were interested in the volume of regions with dust above 30, 50 and 90 K, as these temperatures roughly correspond to the sublimation temperatures of CO, H₂CO and H₂O, respectively (Jørgensen et al. 2015; Ceccarelli et al. 2001; Fraser et al. 2001; Sandford & Allamandola 1993).

These gaseous species are often used as probes of the gas and dust in observations of star-forming regions. We followed Jørgensen et al. (2015) and used the sublimation temperature of CO at 30 K, assuming that the CO is mixed with water-ice on the grain surface, which is more realistic in the protostellar environment than pure CO ice, which has a lower binding energy and thus a lower sublimation temperature of 20 K (Jørgensen et al. 2015).
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Relative differences in volumes of dust above different temperatures. The volume difference between the single star temperature volume, $V_1$, and binary star temperature volume, $V_2$, is shown, normalized to $V_1$, i.e. -0.2 means that the total temperature volume is 20% lower if there are two stars compared to one. The envelope mass range is 0.7–24 $M_\odot$, consistent with the dense core mass range observed in Alves et al. (2007). The relevant mass range for IRAS 16293 is 2.4–4.8 $M_\odot$, (our IRAS 16293 envelope model has a mass of ~4 $M_\odot$). There is a clear tendency for a single star to result in more dust heated above ~50 K.

A walk through the parameter space of the dust reference density $\rho_0$, with $p = 1.7$, was performed with $r_{\text{plat}} = 600$ AU, $L_* = 21 L_\odot$ (if one source) and $L_1 = 10.5, 14.0$ or $18.0 L_\odot$, which means that, accordingly $L_2$ would be $10.5, 7.0$ or $3.0 L_\odot$, in the case of two sources. The most important value is the central density plateau, $\rho_{\text{plat}}$, as the 90 K region, almost all of the 50 K region and the bulk of the 30 K region resides here (Fig. 3). A given volume difference in two different $\rho_{\text{plat}}$ scenarios corresponds to a higher absolute difference in the mass of dust above the given temperature threshold (Fig. 3).

Across the sampled parameter space, as seen in Fig. 3, the single star scenario has larger 90 K volumes than the binary and larger 50 K volumes for regions of higher density. The differences in 50 K go up to ~25%, normalized to the volume of dust above 50 K for a single star. This is more pronounced for the 90 K volumes, where the relative difference is almost up to 40%. In the 30 K region, the binary system generally produced a slightly larger volume than the single star, with modest differences of 2–5% in the total volume. In the case of IRAS 16293, with a total model mass of ~4 $M_\odot$ in this work, the differences are on the order of 5–20% in the 50–90 K volumes.

Therefore, for observations of gaseous species such as CO, H₂CO and H₂O, where the abundances in the gas-phase are expected to increase with sublimation from the dust ice-mantles, two radiation sources should be included in the modeling of IRAS 16293, and of binaries with similar separations, as the total amount of sublimated gas will otherwise be overestimated. If one star of the binary dominates the luminosity, the single-source approximation becomes more valid (Fig. 3) and the overestimation of the sublimated gas becomes less important. The implications of these results for previous 1D chemical models of the hot corinos in IRAS 16293 are discussed further in Sect. 6.

## 4. Dust density structures of IRAS 16293

In this section we present our analysis of the ALMA dust continuum image of IRAS 16293 at 868 µm, where we investigated the presence of disks, or disk-like structures, and a dust filament, using a 3D dust density model. The disk component models were fixed to 150 and 50 AU around source A and B, respectively, based on visual inspection of the angular sizes and emission morphology of the dust continuum. The protostellar masses of source A and B are unknown. Bottinelli et al. (2004b) and Caux et al. (2011) estimate $M_A \sim 1 M_\odot$, assuming the observed gas lines originate from infalling material, while the narrow linewidths observed by Caux et al. (2011) toward source B provide an upper limit of 0.1 $M_\odot$, under the same assumption. Pineda et al. (2012) resolve the region around source A with a 2.2" × 1.0" beam and interpret a Position-Velocity (PV) diagram of methyl formate (CH₃OCHO) and ketene (H₂CCO) as Keplerian rotation around a central object of 0.53 $M_\odot$ and a disk inclination close to edge-on, based on the linewidths (though the inclination is not constrained). Oya et al. (2016) use gas line observations toward source A of higher spatial resolution (0.6" × 0.5" beam) to conclude that a large part of the material around source A is undergoing a rotating collapse, down to 40–60 AU, within which a Keplerian disk is speculated to lie. The inclination of the rotating collapse is not constrained, with best fits to their ballistic model between 30° and 70° (where 90° is edge-on). The derived mass range by Oya et al. (2016) of 0.5–1.0 $M_\odot$, is dependent on the inclination and centrifugal barrier radius. The mass of 0.53 $M_\odot$ reported by Pineda et al. (2012) assumes exactly edge-on configuration and has lower spatial resolution than Oya et al. (2016), thus arguably interpreting material undergoing a rotating collapse as part of a Keplerian disk. Also, the analysis of the velocity gradients on 50–400 AU
scales around source A by Favre et al. (2014), shows a rotating structure which cannot be explained by simple Keplerian rotation around a point mass but needs to take the enclosed mass into account. We used $M_A = 1.0 M_\odot$ in this work. For source B, Pineda et al. (2012) find that the narrow linewidths from the gas in this region are consistent with a face-on disk. The aspect ratio upper limit of 1.1 from the disk emission around source B (Jørgensen et al. 2016) is also consistent with a nearly face-on disk structure. Thus, the analysis by Caux et al. (2011) with a derived mass of $\le 0.1 M_\odot$ seems valid and was used in this work. There have been reports of source A itself being a binary, based on cm emission (Wootten 1989; Pech et al. 2010) and submm emission (Chandler et al. 2005) with a projected distance of roughly 0.4". The PILS program, with a comparable beam (Jørgensen et al. 2016), does not find evidence of the split into separate submm sources reported by Chandler et al. (2005), but this could be due to optically thick dust emission, as the PILS continuum image (868 $\mu$m) is at a shorter wavelength than the Chandler et al. (2005) image (~1 mm). While it is unclear whether source A is a binary or not, we modeled it as a single radiation source.

### 4.1. Disk component model constraints

The PILS interferometric beam of $0.5'' \times 0.5''$ is comparatively large compared to the continuum emission at 868 $\mu$m around the two radiation sources, with the strongest continuum emission extending roughly 0.5" around source B and with a semi-major axis of 1" for the elliptical emission around source A (Fig. 1). Since disk A appears close to edge-on, it should be possible to investigate the vertical distribution of dust with the continuum image. The model continuum emission around source A was fitted with a 2D Gaussian in a $2'' \times 2''$ box centered on the peak flux density pixel on source A, from which the aspect ratio of the semi-major axis to the semi-minor axis was derived from the fitted FWHM in the two dimensions. The observed 868 $\mu$m continuum emission in a $2'' \times 2''$ box around source A was found to have an aspect ratio of 1.64, which we required the model to match within 10% along with the peak central flux density of source A within 10%. The box of $2'' \times 2''$ was chosen to avoid the observed warped structure around source A at larger scales (likely belonging to the dust filament, see Fig. 1), a disk morphology we did not attempt to recreate in this work. We only required our continuum modeling to match the peak flux density of source B within 10%, as more detailed structures are hard to distinguish with our resolution. The column density limit toward the optically thick source B peak continuum emission, which Jørgensen et al. (2016) derive to be $N_{H_2} > 1.2 \times 10^{25} \text{ cm}^{-2}$, was used as well. Any model which satisfied these constraints also needed to reproduce the general features in both the dust continuum emission arc between the sources and the three CO isotopologue maps. Including these, six constraints were available per disk component model (the four emission maps, plus the two individual disk constraints).

### 4.2. Dust filament

The arc of continuum emission (Fig. 1) is assumed to arise from a local dust density enhancement. We modeled this dust filament as an elliptical tube, extending between and enveloping the sources. Around each source, the dust filament was extended spherically up to 130 AU. The density was modeled as a radially dependent power-law distribution of $\rho = \rho_0(r/AU)^{-1/4}$ as

measured from the tube center, where $\rho_0$ is the reference density at 1 AU. As such the filament dust mass density was made to be independent of the distance to either source. The disk components were carved out inside the dust filament model, in a sphere of 150 and 50 AU around source A and B, respectively. No attempts were made to satisfy boundary conditions, resulting in density jumps in the model filament-disk boundary regions. For simplicity, the dust filament structure was kept in the plane of the sky, as were the two protostars. We fixed the filament model with a $\rho_0$ value, which matched the observed central arc emission with several different luminosity ratios (Fig. 4), see Table 1. We performed a sanity check using the maximum observed emission of $0.056 \text{ Jy Beam}^{-1}$ in the filament to calculate a column density based on Eq. (5) in Jørgensen et al. (2007), assuming a dust temperature of 30 K. The dust column density of the filament was derived to be $9.9 \text{ g cm}^{-2}$, while the model column density was $9.7 \text{ g cm}^{-2}$, assuming a gas-to-dust mass ratio of 100. This match is not surprising as we used the same opacities as Jørgensen et al. (2007) and the optical depth through the model filament is low (~0.2).

### 4.3. Disk component models

The disk components were modeled with both a PP-disk and a rotating toroid model. In the PP-disk model, the dust distribution is based on hydrostatic equilibrium with the gravity of the host star (e.g., Andrews et al. 2012).

The following equations all use cylindrical coordinates. The disk scale-height, $H_d$, is

$$H_d(r) = H_0(r/r_0)^{1+\Psi}, \qquad (2)$$

where $r$ is the radius, $H_0$ is the scale-height at radius $r_0$ and $\Psi$ is the flaring constant, controlling the vertical distribution of the dust above the disk midplane. The disk surface density, $\Sigma$, is

$$\Sigma(r) = \Sigma_0(r/r_0)^{-p}, \qquad (3)$$

where $p$ is a gradient parameter. Finally the dust density is

$$\rho_d(r, z) = \frac{\Sigma(r)}{H_d \sqrt{2\pi}} e^{-z^2/(2H_d^2)}, \qquad (4)$$

where $z$ is the height above the disk midplane.

The disk inner radius was fixed to 1 AU for both disks, with the dust density within 1 AU set to $10^{-25} \text{ g cm}^{-3}$ ($n_{H_2} \sim 3 \text{ cm}^{-3}$), effectively assuming a dust and gas evacuated inner gap. After the dust temperature of a given model was estimated in RADMC-3D, subsequent ray tracing produced a 868 $\mu$m continuum image, which, after convolution with a $0.5'' \times 0.5''$ beam using MIRIAD, was compared with the 868 $\mu$m ALMA continuum image in Fig. 1. The disk orientation for source B was fixed to face-on, while the inclination of PP-disk A was a free parameter (Table 1).

Due to the high number of grid cells required to resolve the relevant regions and high photon numbers required for reliable dust temperature results, but most importantly the high densities in the inner disk regions of this work, the thermal Monte Carlo calculation execution time, on the available computing facilities, for a single model was very long, sometimes on the order of days. This problem was compounded by the iterative allocation of dust opacities, which requires three separate dust temperature estimations in RADMC-3D. Thus, standard model optimization approaches such as a Markov chain Monte Carlo (MCMC) or
<table><caption>**Table 2.** Example dust model parameters.</caption><thead><tr><th>PP-disk model</th><th>L<sub>A</sub> [L<sub>☉</sub>]</th><th>L<sub>B</sub> [L<sub>☉</sub>]</th><th>Σ<sub>0,A</sub> [g cm<sup>-2</sup>]</th><th>Σ<sub>0,B</sub> [g cm<sup>-2</sup>]</th><th>H<sub>0,A</sub> [AU]</th><th>H<sub>0,B</sub> [AU]</th><th>i<sub>A</sub> [°]</th></tr></thead><tbody><tr><td>DM 1</td><td>18.0</td><td>3.0</td><td>1.0</td><td>3.0</td><td>1.25</td><td>1.5</td><td>55</td></tr><tr><td>DM 2</td><td>14.0</td><td>7.0</td><td>1.0</td><td>0.05</td><td>0.85</td><td>2.0</td><td>50</td></tr><tr><td>DM 3</td><td>10.5</td><td>10.5</td><td>1.0</td><td>0.05</td><td>1.25</td><td>2.0</td><td>60</td></tr></tbody></table><table><thead><tr><th>Rotating toroid model</th><th>L<sub>A</sub> [L<sub>☉</sub>]</th><th>L<sub>B</sub> [L<sub>☉</sub>]</th><th>M&#x0307;<sub>A</sub> [M<sub>☉</sub> yr<sup>-1</sup>]</th><th>M&#x0307;<sub>B</sub> [M<sub>☉</sub> yr<sup>-1</sup>]</th><th>r<sub>c,A</sub> [AU]</th><th>r<sub>c,B</sub> [AU]</th></tr></thead><tbody><tr><td>RTM 1</td><td>18.0</td><td>3.0</td><td>3.75 &times; 10<sup>-6</sup></td><td>2.75 &times; 10<sup>-6</sup></td><td>75</td><td>20</td></tr><tr><td>RTM 2</td><td>10.5</td><td>10.5</td><td>2.1 &times; 10<sup>-6</sup></td><td>2.1 &times; 10<sup>-6</sup></td><td>50</td><td>20</td></tr><tr><td>RTM 3</td><td>20.0</td><td>1.0</td><td>3.8 &times; 10<sup>-6</sup></td><td>4.6 &times; 10<sup>-6</sup></td><td>75</td><td>20</td></tr><tr><td>RTM 4</td><td>7.0</td><td>14.0</td><td>3.0 &times; 10<sup>-6</sup></td><td>2.0 &times; 10<sup>-6</sup></td><td>75</td><td>20</td></tr><tr><td>RTM 5</td><td>3.0</td><td>18.0</td><td>3.0 &times; 10<sup>-6</sup></td><td>2.0 &times; 10<sup>-6</sup></td><td>75</td><td>20</td></tr></tbody></table>

a $\chi^2$ analysis of a large, uniform grid of model parameters were not feasible in this work and no parameter confidence intervals were derived. Instead, we identified parameter spaces not compatible with the current constraints and derived conclusions from these parameter space exclusions.

### 4.4. Model analysis

Starting with the PP-disk model, large steps were taken in the parameter space of $\Sigma_0$ and $H_0$ for the disk around each source, and inclination for disk A (Table 1), while $p$ and $\Psi$ were fixed to 0.5 and 0.25, respectively (see Appendix D for more details). The final number of free parameters in the PP-disk model was four for disk A ($\Sigma$, $H_0$, $i$ and source A luminosity, $L_A$) and two for disk B ($\Sigma$ and $H_0$), since the source B luminosity, $L_B$, was determined by $L_A$.

When only using the continuum image, the PP-disk model did not offer any meaningful constraints on disk A and B. Due to the degeneracy between $H_0$ and $i$, only models with inclinations below 55° could be excluded. The peak flux values were degenerate with different $\Sigma_0$, $H_0$, $L_*$ and $i$ values and the column density toward source B did not offer any constraints either, as almost all of our investigated parameter space satisfied the $N_{H_2}$ column density condition.

Due to previous observations of infall toward both sources (Caux et al. 2011; Chandler et al. 2005; Pineda et al. 2012; Zapata et al. 2013; Oya et al. 2016), we considered the rotating toroid model. The dust density distribution of this model can be found semi-analytically by defining discrete streamlines from the rotating envelope, which spiral inwards and impact opposing streamlines in the midplane, creating a disk structure. From Hartmann (2009),

$$ \rho(r, \theta) = \frac{\dot{M} \left(1 + \frac{\cos \theta}{\cos \theta_0}\right)^{-1/2}}{4\pi \sqrt{GM_* r^3}} \left( \frac{\cos \theta}{\cos \theta_0} + \frac{2\cos^2 \theta_0}{r/r_c} \right)^{-1}, \quad (5) $$

where $\dot{M}$ is the mass accretion rate onto the central object, $\theta_0$ is the starting polar angle of the streamline in question and $r_c$ is the centrifugal radius. The density in each grid cell in our model was found by predefining $n_{stream}$ streamlines with discrete $\theta_0$ values and evaluating them at $n_r$ radial points, along the streamline. A nearest-neighbor search was performed, with the dust density of the RADMC-3D grid cell defined to be that of the nearest defined streamline point. The rotating toroid model is effectively a 2D density model and a 3D velocity model (Hartmann 2009).

The rotating toroid model has a radially extended dust distribution with $\rho \propto r^{-0.5}$ and a much more vertically extended dust

distribution than the exponentially decreasing vertical distribution of the PP-disk model (Eq. (4)). Consequently, the rotating toroid model has a much lower optical depth through the mid-plane, hence higher dust temperatures, than a PP-disk model with, for example, $H_0 = 0.85$ AU (this example value was taken from the range in Andrews et al. 2010).

The centrifugal radius of the collapse around source B ($r_{c,B}$) was a free parameter between 5–25 AU, to allow a compact toroid within the defined model disk B regime of 50 AU, fixed to face-on inclination. We fixed the inclination of the rotating toroid model around source A to edge-on and used $r_{c,A} = 75$ AU, as it matched the aspect ratio at this inclination. The free parameters of the source A rotating toroid model were $\dot{M}_A$ and $L_A$, as $r_{c,A}$ and $M_A$ were fixed (Table 1), while the source B rotating toroid model had $\dot{M}_B$ and $r_{c,B}$ as free parameters (since $L_B$ was determined by $L_A$). We found that the observed peak flux density toward source A and B could not constrain the luminosity or the mass accretion rate of either source, as these parameters are degenerate. The minimum column density toward source B did not offer any constraint. In conclusion, the rotating toroid model can satisfy the disk constraints, but cannot constrain $L_A$ or $L_B$, using continuum data alone (Table 1).

## 5. Gas line emission modeling

In order to constrain the nature of the disk-like emission further, and to obtain more constraints on the source luminosities, the $J = 3-2$ line emission from $^{13}$CO, C$^{17}$O and C$^{18}$O of select dust models were modeled in LIME (Brinch & Hogerheijde 2010), a line radiative transfer code. The dust temperature output from RADMC-3D was imported to LIME (Appendix C) and by using the standard gas-to-dust mass ratio of 100 (Bohlin et al. 1978) together with the given dust density model, the molecular hydrogen number density $n_{H_2}$ was found. The gas line emission depends on the number density of the molecule undergoing a rotational transition, the CO isotopologue in our case, the number density of the other molecules in the gas which excite the molecule collisionally and the kinetic temperature which changes the collisional rates (Brinch & Hogerheijde 2010) and thus the excitation of the CO isotopologue. In LIME, the number density of the colliders is approximated as $n_{H_2}$, as $H_2$ is by far the most abundant species in the dense molecular cores where star formation takes place (e.g., Genzel & Stutzki 1989). LIME performs non-LTE (Local Thermodynamic Equilibrium) line radiative transfer modeling, where the rotational level populations are found iteratively. The CO isotopologue data files containing collisional rate coefficients (Yang et al. 2010)
<figure><img src="image_5.png" /><figcaption>Fig. 4. Observed and synthetic dust continuum emission at 868 µm for the different physical models summarized in Table 2. Image is in log-scale and the contours are logarithmically divided between 0.02 and 2 Jy beam⁻¹. DM denotes the PP-disk model, while RTM denotes the rotating toroid model.</figcaption></figure>

were downloaded from the LAMDA<fnref n="2" /> database (Schöier et al. 2005). The collisional rate coefficients of the CO isotopologues with H₂ in these files include collisions with both ortho- and para-H₂.

In this work, we assumed that the dust and gas temperatures are coupled, which is usually a good approximation in the dense inner environments of the envelope (Schöier et al. 2002), that were modeled here. The radius of the LIME model was set to 2 × 10³ AU, as more distant, cold regions in the envelope will not contribute much emission due to freeze-out of CO.

Optimally, the absorption features around each source and the major features of the interbinary CO isotopologue emission seen in Fig. 2 were recreated in the synthetic moment zero maps, which would imply that our dust density model is consistent with the observed gas line emission. We did not attempt to match the exact absorption level as it depends not only on the

molecule abundance but also on the gas-to-dust mass ratio and dust opacity, neither of which were constrained in this work.

In the case of a PP-disk model, the velocity in the disk around each source was modeled as Keplerian rotation. Here the central stellar mass becomes important if there is any disk inclination toward the observer, as radial velocities cause spectral broadening, reducing the line optical depth. This, in turn, will result in higher integrated emission if there is any line self-absorption, which we expect around each source, based on the observations seen in Fig. 2. In the case of a rotating toroid model, the 3D velocity structure from Terebey et al. (1984) was used. For all regions outside the disk components, the velocity was modeled as simple free-fall,

$$v(r) = \sqrt{GM_{\text{tot}}/r}, \qquad (6)$$

where r is the distance to the model center, G is the gravitational constant and Mₜot is the combined mass of source A and B. This approximation is valid at kAU scales but breaks down in the innermost region. The dust in the interbinary region should

<footnote n="2">² http://home.strw.leidenuniv.nl/~moldata/</footnote>
<table><caption>**Table 3.** LIME model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>Grid outer radius</td><td>...</td><td>2 &times; 10<sup>3</sup> AU</td></tr><tr><td>Grid start radius</td><td>...</td><td>30 AU</td></tr><tr><td>n<sub>nodes</sub></td><td>Number of grid nodes in the LIME calculations</td><td>7 &times; 10<sup>4</sup></td></tr><tr><td>CO/H<sub>2</sub></td><td>Standard ISM abundance of CO</td><td>10<sup>-4</sup></td></tr><tr><td><sup>13</sup>CO/H<sub>2</sub></td><td>Standard ISM abundance of <sup>13</sup>CO</td><td>1.30 &times; 10<sup>-6</sup></td></tr><tr><td>C<sup>18</sup>O/H<sub>2</sub></td><td>Standard ISM abundance of C<sup>18</sup>O</td><td>1.79 &times; 10<sup>-7</sup></td></tr><tr><td>C<sup>17</sup>O/H<sub>2</sub></td><td>Standard ISM abundance of C<sup>17</sup>O</td><td>5.58 &times; 10<sup>-8</sup></td></tr><tr><td>Freeze out factor</td><td>...</td><td>0.01</td></tr><tr><td>Gas-to-dust-mass ratio</td><td>Used for all conversions from dust to gas mass</td><td>100</td></tr><tr><td>CO isotopologues abundances of all models</td><td></td><td></td></tr><tr><td>X<sub>13CO</sub></td><td><sup>13</sup>CO abundance</td><td>2.6 &times; 10<sup>-7</sup></td></tr><tr><td>X<sub>C<sup>18</sup>O</sub></td><td>C<sup>18</sup>O abundance</td><td>3.58 &times; 10<sup>-8</sup></td></tr><tr><td>X<sub>C<sup>17</sup>O</sub></td><td>C<sup>17</sup>O abundance</td><td>1.12 &times; 10<sup>-8</sup></td></tr></tbody></table>

**Notes.** Given abundances assume that C and O atoms maximally combine into CO.

arguably be drifting toward either source, in accordance with the local gravitational field. Since only one dust opacity model can be set in LIME (v. 1.5), bare-grain dust opacities were used in the entire model, as the dust grains in the disk regions are mostly above 90 K in the rotating toroid models, and in the PP-disk model atmosphere. If we had used icy-grain opacities in the LIME models instead, the reduced absorption from the dust would result in increased flux density globally in the integrated gas emission maps, with the largest difference in the disk regions.

## 5.1. Abundance modeling

It is assumed that the bulk of the CO isotopologue emission originates from molecules sublimated from the dust grains, represented with an abundance jump model in this work, similar to some of the previous chemical modeling of IRAS 16293 (Schöier et al. 2002; Schöier et al. 2004; Coutens et al. 2012). The individual CO isotopologue abundances relative to H₂, the main collision partner in our modeling, were found by using the C and O isotopologue abundances given by Wilson & Rood (1994; Table 3). Standard abundance values from Wilson & Rood (1994) were applied to regions with T<sub>dust</sub> ≥ 30 K where the CO is sublimated from the dust grains, and a drop factor of 10<sup>-2</sup> was used in regions with T < 30 K, to imitate freeze-out. Examining the temperature complexity of a typical rotating toroid model (Fig. F.4), when disks and the dust filament are added to the model, we can see the expected emission shape from the CO isotopologues in Fig. 5, based on the T = 30 K contour.

## 5.2. Modeled CO isotopologue maps

To constrain the luminosity ratio of the two protostars, we modeled the rotating toroid model at different luminosity values (Table 5), with parameter sets that satisfied the dust continuum constraints. The resulting synthetic CO isotopologue maps (Fig. 6) indicate that L<sub>A</sub> ≥ 18 L<sub>☉</sub> in order to reproduce the observed tendency of concentrated emission toward source A. We can deduce the luminosity dominance of source A from the observations alone, as the truncation in the observed C<sup>17</sup>O and C<sup>18</sup>O maps occurs closer to source B than source A, which

requires source A to be more luminous, if we assume this truncation is due to the dust grain temperature dropping below the sublimation temperature of CO. L<sub>A</sub> ≥ 18 L<sub>☉</sub> is also consistent with the dust continuum model results (Fig. 4). Alternatively, the increase in continuum and gas line emission toward source A could be due to geometric effects, such as an inclination of the stellar plane so source A is closer to us than source B, or simply that more dust material exists in the filament toward source A. However, an increase in dust density would require L<sub>A</sub> to be lowered, in order to match the observed continuum emission as well as increasing the optical depth. Both instances lower the gas temperature in the region near source A, which creates a truncated region wherein CO is sublimated closer to source A than source B. In this sense, introducing more material in our filament model toward source A should require L<sub>A</sub> to be even higher than we currently find in order to match the observed CO isotopologue emission. The inferred low luminosity of source B implies that disk B is denser than disk A in order to dominate the dust submillimeter emission, which fits well with the observed absorption in the CO isotopologue emission around source B, compared to the vicinity of source A.

A representative rotating toroid model which satisfied the dust continuum constraints with L<sub>A</sub> = 18 L<sub>☉</sub> was modeled in LIME, together with a representative of the best matching PP-disk model at L<sub>A</sub> = 14 L<sub>☉</sub> (since L<sub>B</sub> < 7 L<sub>☉</sub> could not reproduce the observed peak flux density in our H<sub>0</sub> parameter space, though this does not rule out L<sub>B</sub> < 7 L<sub>☉</sub> by itself, see Fig. F.2) can be seen in Fig. 5. The first attempt with the standard abundances described in Table 3, resulted in significant absorption in a region toward source B as expected, but also with some unexpected absorption toward source A. For the rotating toroid model, the C<sup>18</sup>O model emission in the disk A component has evident self-absorption throughout the midplane. After lowering the abundances by a factor of five, the model C<sup>18</sup>O emission structure was more reminiscent of the observations, while C<sup>17</sup>O did not reproduce the observed peaked emission in the disk A component at any of the attempted abundances (Fig. 5).

The model emission qualitatively matches the observed <sup>13</sup>CO bimodal emission structure above and below the disk midplane, as well as the two observed C<sup>18</sup>O emission peaks near the corners of the disk A component, but with less absorption in C<sup>18</sup>O toward source A than observed. The rotating toroid model cannot, however, reproduce the C<sup>17</sup>O centrally concentrated
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Zeroth moment maps of the observed and synthetic CO isotopologue gas line emission. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam⁻¹ km s⁻¹. The RA and Dec offsets are relative to the phase center of the observations. See Table 4 for abundances. DM denotes the PP-disk model and RTM denotes the rotating toroid model.</figcaption></figure>

emission structure, but instead results in a bimodal distribu-
tion as seen in C¹⁸O, near r<sub>c,A</sub> (r<sub>disk</sub> = 75 AU). The PP-disk
model fails to produce the emission peaks around source A, in
all isotopologues, as the PP-disk A region heavily absorbs the
line emission, even when the standard abundances are lowered
by a factor of five, despite the fact that we used the PP-disk
model parameters with the lowest optical depth while matching
continuum constraints.

When the PP-disk model standard abundances were lowered
by a factor of ten or 15 instead (Fig. F.1) the C¹⁷O emission from
the filament became much weaker than observed, while the disk
A region remained in absorption. The model absorption is then
likely not due to line self-absorption, but rather originates from
the optically thick dust. This suggests that the disk-like emission
around source A is composed of a relatively diffuse, vertically
extended structure.

However, if larger amounts of CO are in front of disk A,
on similar spatial scales, then it could explain the extra line

emission in the observations compared to the PP-disk model
emission. Since the continuum emission is consistent with a PP-
disk model, we cannot exclude the possibility that the disk-like
emission structure around source A could be from a flat disk,
based on our model and observations alone. It is also possible
that the presence of material in front of disk A could allow the
rotating toroid models to use the standard ISM CO abundance
and still qualitatively match the observations.

6. Discussion

The combined constraints from the dust continuum and CO
isotopologue emission, necessitate very high PP-disk model
scale-heights around source B, H₀ > 3.0 AU at 10 AU, which
is suggestive of a nonsettled disk with a large vertical distribu-
tion, meaning that it is indicative of a structure similar to that
of the diffuse, rotating toroid model. These scale-heights are
unprecedented for Class II disks, but similar to the one reported
<table><caption>**Table 4.** LIME model parameters of Fig. 5.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td colspan="3"><em>Rotating toroid model</em></td></tr><tr><td>L<sub>A</sub></td><td>Luminosity of source A</td><td>18 L<sub>☉</sub></td></tr><tr><td>L<sub>B</sub></td><td>Luminosity of source B</td><td>3 L<sub>☉</sub></td></tr><tr><td>M&#x0307;<sub>A</sub></td><td>Mass accretion rate</td><td>3.75 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>M&#x0307;<sub>B</sub></td><td>Mass accretion rate</td><td>2.75 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>r<sub>c,A</sub></td><td>Centrifugal radius</td><td>75 AU</td></tr><tr><td>r<sub>c,B</sub></td><td>Centrifugal radius</td><td>20 AU</td></tr><tr><td colspan="3"><em>PP-disk model</em></td></tr><tr><td>L<sub>A</sub></td><td>Luminosity of source A</td><td>14 L<sub>☉</sub></td></tr><tr><td>L<sub>B</sub></td><td>Luminosity of source B</td><td>7 L<sub>☉</sub></td></tr><tr><td>Σ<sub>0,A</sub></td><td>Dust surface density at 10 AU</td><td>0.6 g cm<sup>-2</sup></td></tr><tr><td>Σ<sub>0,B</sub></td><td>Dust surface density at 10 AU</td><td>5.0 g cm<sup>-2</sup></td></tr><tr><td>H<sub>0,A</sub></td><td>Scale-height at 10 AU</td><td>1.5 AU</td></tr><tr><td>H<sub>0,B</sub></td><td>Scale-height at 10 AU</td><td>2.0 AU</td></tr><tr><td>i<sub>A</sub></td><td>Disk A inclination</td><td>60°</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Observed and modeled C<sup>17</sup>O zeroth moment maps, contour levels are divided logarithmically from 0.5 to 7.3 Jy beam<sup>-1</sup> km s<sup>-1</sup>. All models are rotating toroid models matching the dust continuum peak flux densities for disk A and B.</figcaption></figure>

by Tobin et al. (2013), who find a best-fit PP-disk model to the Class 0 disk around L1527 with a scale-height of 48 AU at 100 AU. In comparison, when converting $H_0$ to the scale-height at 100 AU, the PP-disk model around source B needs a scale-height of more than 53 AU at 100 AU, when $L_B = 3 L_☉$, in order to reach the peak flux density value (Fig. F.2). This value is significantly higher than the scale-heights at 100 AU of Class I/II sources, which typically are 3–20 AU (Wolf et al. 2008; Andrews et al. 2009, 2010). Tobin et al. (2013) propose

that the large scale-height of their best fitting PP-disk model of L1527 is approximating material from the envelope falling onto the disk and that this unusually vertically extended structure could reflect a stage before effective dust settling. The disk-like structure around source B could be similar to the L1527 disk in this regard, as infall toward source B has previously been detected.

The CO isotopologue moment maps are not entirely qualitatively matched by a rotating toroid, but the fact that the standard
<table><caption>**Table 5.** LIME model parameters of Fig. 6.</caption><thead><tr><th>Model</th><th>M&#x0307;<sub>A</sub> [M<sub>&#x2213;</sub> yr<sup>-1</sup>]</th><th>M&#x0307;<sub>B</sub> [M<sub>&#x2213;</sub> yr<sup>-1</sup>]</th></tr></thead><tbody><tr><td>L<sub>A</sub> = 3 L<sub>&#x2213;</sub></td><td>9.0 &times; 10<sup>-6</sup></td><td>1.3 &times; 10<sup>-6</sup></td></tr><tr><td>L<sub>A</sub> = 7 L<sub>&#x2213;</sub></td><td>5.5 &times; 10<sup>-6</sup></td><td>1.4 &times; 10<sup>-6</sup></td></tr><tr><td>L<sub>A</sub> = 10.5 L<sub>&#x2213;</sub></td><td>4.6 &times; 10<sup>-6</sup></td><td>1.6 &times; 10<sup>-6</sup></td></tr><tr><td>L<sub>A</sub> = 14 L<sub>&#x2213;</sub></td><td>4.2 &times; 10<sup>-6</sup></td><td>1.85 &times; 10<sup>-6</sup></td></tr><tr><td>L<sub>A</sub> = 18 L<sub>&#x2213;</sub></td><td>3.75 &times; 10<sup>-6</sup></td><td>2.75 &times; 10<sup>-6</sup></td></tr><tr><td>L<sub>A</sub> = 20 L<sub>&#x2213;</sub></td><td>3.8 &times; 10<sup>-6</sup></td><td>4.6 &times; 10<sup>-6</sup></td></tr></tbody></table>

**Notes.** All models are rotating toroids, using $r_{c,A} = 75$ AU and $r_{c,B} = 20$ AU.

CO isotopologue abundances needed to be lowered by a factor of five to avoid self-absorption in C¹⁸O is interesting. Anderl et al. (2016) use abundances an order of magnitude lower than the standard CO isotopologue abundances in order to match their observations of regions interior to the CO snowline in a survey of nearby Class 0 protostars. In the vicinity of protostars and accretion shocks, the radiation field and chemistry is markedly different from the ISM, so deviations from the standard ISM abundances in the inner regions of a disk are not surprising (e.g., see Visser et al. 2009). However, because we base this necessity for a lower abundance on a single species, and since the self-absorption might simply reflect that the model is a poor approximation, a more thorough, quantitative investigation would be needed before making any firm conclusions on the CO abundance in the inner regions of IRAS 16293.

It has been suggested that a fraction of the CO is desorbed together with H₂O into the gas-phase at 90 K (Anderl et al. 2016), due to trapping in the water-ice, which would allow higher abundances in the gas model for regions with $T_{\text{dust}} \ge 90$ K than regions with 30 K $\le T_{\text{dust}} < 90$ K. This should lead to an emission structure closer to that of the observed, with a more abrupt peak toward source A and should be investigated in future modeling efforts.

Considering the unusual scale-heights of the PP-disk model, necessary to produce the peak flux density around source B, a rotating toroid model is perhaps more representative of the actual dust distribution around source B than a PP-disk model, as it better approximates the infalling material, unambiguously detected toward source B (Chandler et al. 2005; Pineda et al. 2012; Zapata et al. 2013) as well as shocks (Jørgensen et al. 2011) which would be expected with a rotating collapse. It is also possible that both disks A and B are transition structures, between an early rotating toroid structure and a later Keplerian disk, perhaps similar to the structure suggested by Oya et al. (2016).

In terms of the impact of the binarity of IRAS 16293 on earlier works using a single radiation source, we can consider Schöier et al. (2002), who assume a single radiation source with a spherical density model down to scales where $T_{\text{dust}} = 300$ K. The inferred luminosity dominance of source A from our work makes the single source approximation more valid than a luminosity ratio closer to unity.

When looking at Fig. F.4, it appears that the approximation of a single central source is not dramatically erroneous, as source A dominates the radiation, but rather that a spherically symmetric density model breaks down in the inner regions where the temperatures exceed 90 K. At the scales of 10''–15'' used by Schöier et al. (2002), spherically symmetrical density models are

good approximations of the outer envelope of IRAS 16293 and the binary protostar is well approximated to a single source at the distances of several kAU, as well as being unresolved. However, for radiative transfer modeling of the innermost region including the two sources, their disks as well as the conjoining dust filament, subarcsec resolution data is needed together with full 3D radiative transfer modeling to capture the complexity of the temperature structure. While previous radiative transfer models of IRAS 16293 do not describe the full picture, they can arguably be useful approximations to the hot corino material around source A. If source A is indeed dominant in terms of luminosity, then previous chemical modeling using a single star focusing on hot corino chemistry (Schöier et al. 2002; Bisschop et al. 2008; Coutens et al. 2012) will effectively have described the chemistry around source A, with the hot corino chemistry around source B acting as a “contaminant” to the signal from source A. However, for an updated model of the hot corino chemistry involving COMs in IRAS 16293, a new chemical model is needed which includes the disks, filamentary structure and radiation fields with $L_A \sim 18 L_\odot$ and $L_B \sim 3 L_\odot$.

## 7. Conclusions

We have constructed a 3D model of the dust and gas environment of IRAS 16293 and obtained a qualitative match to the 868 µm dust continuum and $J = 3-2$ line emission of C¹³CO, C¹⁸O and C¹⁷O from the PILS survey. This is the first model of IRAS 16293 to include multiple radiation sources as well as complex dust structures, such as disks and a dust filament. We included two different disk-like density structures; a disk structure (a true Keplerian disk) and the density solution to an infalling, rotating collapse, a disk-like structure resembling a rotating toroid. We determine that the dust density structure around source B must be vertically extended and puffed up, as a PP-disk model around source B with low to normal scale-heights cannot fulfill all the observational constraints. The disk structure around source A could not be constrained. The results from this 3D dust and gas model suggest that the source A luminosity is far higher than that of source B, based on the emission morphology of the dust and CO isotopologue line emission between the sources. This model is an important increase in complexity compared to previous works and should be followed up with higher angular resolution observations of the disks and the dust filament, for further investigation and confirmation of our results. We also performed a more generic investigation of the temperature structure differences between a single star and a binary in a simple dense envelope with a density plateau within 600 AU. The main conclusions are the following:

- The temperature distribution depends on the density structure and source luminosities. A single star will have higher volumes and masses of dust above ~50 K compared to a binary with the same combined luminosity, while the differences in dust with 30–50 K are negligible in a spherical 1D density model.

- Previous radiative transfer models of IRAS 16293 using a spherically symmetric density distribution and a single, central radiation source are valid at the ~kAU scales, but not in the interbinary region, where complex 3D dust structures and two radiation sources need to be included.

- A 3D model with disks or disk-like structures around source A and B in IRAS 16293 and a simple dust filament
between the sources, can qualitatively match the spatial dis-
tribution of the 868 µm dust continuum and CO isotopologue
line emission.

- The density structure around source A in our model is not constrained, as both a flat, inclined disk and an edge-on puffed-up disk structure (a rotating toroid model) can match the continuum peak emission and observed aspect ratio, while neither could match the observed C¹⁷O isotopologue emission.

- The dust density structure around source B in our model has to be vertically extended. Such a structure can be obtained with both a PP-disk model with an unusually high scale-height and with a rotating toroid model, which naturally has a puffed-up disk-like structure.

- The rotating toroid model abundances of the CO isotopologues need to be a factor of five lower than typical ISM values, in order to match the C¹⁸O observations, similar to the low inferred CO abundance for other embedded protostars (e.g., Anderl et al. 2016). This result needs to be investigated further, as the low abundances could be model dependent.

- Source A is likely much more luminous than source B. Our 3D density model (with coarse parameter space steps) requires the source A luminosity to be a factor of six larger; $L_A \sim 18 L_\odot$ and $L_B \sim 3 L_\odot$.

This work illustrates the importance of using complex,
full 3D radiative transfer modeling together with subarcsec-
ond observations of young stellar objects such as IRAS 16293
in order to constrain the physical parameters of such Class 0
disks and their environment. It also illustrates the need for
even higher resolution observations of the gas and dust in
the disks around source A and B in order to do a full quan-
titative investigation with parameter confidence intervals, to
determine the exact CO isotopologue abundances and the param-
eters of the dust filament enveloping both sources. We also
suggest that gas dynamics, including outflows and magnetic
fields, are investigated in future modeling efforts. Having a
detailed and robust physical model is critical for the inter-
pretation of ALMA spectral data of the various hot corino
species.

Acknowledgements. The authors are grateful to Søren Frimann, Jonathan Ramsey and Attila Juhasz for useful discussions about 3D radiative transfer modeling. We would also like to thank the anonymous referee for helpful suggestions, which improved the presentation of this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00278.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The group of J.K.J. acknowledges support from a Lundbeck Foundation Group Leader Fellowship as well as the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 646908) through ERC Consolidator Grant "S4F". Research at the Centre for Star and Planet Formation is funded by the Danish National Research Foundation. M. N. Drozdovskaya acknowledges the financial support of the Center for Space and Habitability (CSH) Fellowship and the IAU Gruber Foundation Fellowship. A. Coutens postdoctoral grant is funded by the ERC Starting Grant 3DICE (grant agreement 336474). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013). This research has made use of NASA's Astrophysics Data System.

References

Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17

Alves, F. O., Vlemmings, W. H. T., Girart, J. M., & Torrelles, J. M. 2012, A&A, 542, A14

Anderl, S., Maret, S., Cabrit, S., et al. 2016, A&A, 591, A3

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241

Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162

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

Bisschop, S. E., Jørgensen, J. K., Bourke, T. L., Bottinelli, S., & van Dishoeck, E. F. 2008, A&A, 488, 959

Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bottinelli, S., Ceccarelli, C., Lefloch, B., et al. 2004a, ApJ, 615, 354

Bottinelli, S., Ceccarelli, C., Neri, R., et al. 2004b, ApJ, 617, L69

Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25

Cassen, P., & Moosman, A. 1981, Icarus, 48, 353

Caux, E., Kahane, C., Castets, A., et al. 2011, A&A, 532, A23

Ceccarelli, C., Loinard, L., Castets, A., et al. 2001, A&A, 372, 998

Chandler, C. J., Brogan, C. L., Shirley, Y. L., & Loinard, L. 2005, ApJ, 632, 371

Chiang, H.-F., Looney, L. W., Tassis, K., Mundy, L. G., & Mouschovias, T. C. 2008, ApJ, 680, 474

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

Correia, J. C., Griffin, M., & Saraceno, P. 2004, A&A, 418, 607

Coutens, A., Vastel, C., Caux, E., et al. 2012, A&A, 539, A132

Crimier, N., Ceccarelli, C., Maret, S., et al. 2010, A&A, 519, A65

Duchêne, G., Bouvier, J., Bontemps, S., André, P., & Motte, F. 2004, A&A, 427, 651

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015]

Dunham, M. M., Evans, II, N. J., Terebey, S., Dullemond, C. P., & Young, C. H. 2010, ApJ, 710, 470

Dutrey, A., di Folco, E., Guilloteau, S., et al. 2014, Nature, 514, 600

Favre, C., Jørgensen, J. K., Field, D., et al. 2014, ApJ, 790, 55

Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Williams, D. A. 2001, MNRAS, 327, 1165

Genzel, R., & Stutzki, J. 1989, ARA&A, 27, 41

Hartmann, L. 2009, Accretion Processes in Star Formation, 2nd edn. (Cambridge: Cambridge University Press)

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2005a, ApJ, 632, 973

Jørgensen, J. K., Lahuis, F., Schöier, F. L., et al. 2005b, ApJ, 631, L77

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007, ApJ, 659, 479

Jørgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2009, A&A, 507, 861

Jørgensen, J. K., Bourke, T. L., Nguyen Luong, Q., & Takakuwa, S. 2011, A&A, 534, A100

Jørgensen, J. K., Favre, C., Bisschop, S. E., et al. 2012, ApJ, 757, L4

Jørgensen, J. K., Visser, R., Williams, J. P., & Bergin, E. A. 2015, A&A, 579, A23

Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, A117

Klapper, G., Lewen, F., Gendriesch, R., Belov, S., & Winnewisser, G. 2000, J. Mol. Spectrosc., 201, 124

Klapper, G., Lewen, F., Gendriesch, R., Belov, S. P., & Winnewisser, G. 2001, Z. Naturforsch. A, 56

Klapper, G., Surin, L., Lewen, F., et al. 2003, ApJ, 582, 262

Knude, J., & Hog, E. 1998, A&A, 338, 897

Kristensen, L. E., Klaassen, P. D., Mottram, J. C., Schmalzl, M., & Hogerheijde, M. R. 2013, A&A, 549, L6

Kuan, Y.-J., Huang, H.-C., Charnley, S. B., et al. 2004, ApJ, 616, L27

Li, Z.-Y., Banerjee, R., Pudritz, R. E., et al. 2014, in Protostars and Planets VI (Tucson, AZ: Univ. Arizona Press), 173

Lindberg, J. E., Jørgensen, J. K., Brinch, C., et al. 2014, A&A, 566, A74

Loinard, L., Torres, R. M., Mioduszewski, A. J., & Rodríguez, L. F. 2008, ApJ, 675, L29

Lombardi, M., Lada, C. J., & Alves, J. 2008, A&A, 480, 785

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

Makiwa, G. 2014, Ph.D. thesis, University of Lethbridge (Canada)

Maury, A. J., Belloche, A., André, P., et al. 2014, A&A, 563, L2

Mundy, L. G., Wootten, A., Wilking, B. A., Blake, G. A., & Sargent, A. I. 1992, ApJ, 385, 306

Murillo, N. M., Lai, S.-P., Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2013, A&A, 560, A103

Nelson, A. F., & Marzari, F. 2016, ApJ, 827, 93

Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ApJ, 796, 131

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

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

Pech, G., Loinard, L., Chandler, C. J., et al. 2010, ApJ, 712, 1403

Pineda, J. E., Maury, A. J., Fuller, G. A., et al. 2012, A&A, 544, L7
Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

Saftly, W., Baes, M., & Camps, P. 2014, A&A, 561, A77

Sakai, N., Oya, Y., Sakai, T., et al. 2014, ApJ, 791, L38

Sandford, S. A., & Allamandola, L. J. 1993, ApJ, 417, 815

Schöier, F. L., Jørgensen, J. K., van Dishoeck E. F., & Blake, G. A. 2002, A&A,
390, 1001

Schöier, F. L., Jørgensen, J. K., van Dishoeck, E. F., & Blake, G. A. 2004, A&A,
418, 185

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005,
A&A, 432, 369

Shirley, Y. L., Evans, II, N. J., & Rawlings, J. M. C. 2002, ApJ, 575, 337

Shirley, Y. L., Huard, T. L., Pontoppidan, K. M., et al. 2011, ApJ, 728, 143

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2012, Nature, 492, 83

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2013, ApJ, 771, 48

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016, ApJ, 818, 73

van der Tak, F. F. S., van Dishoeck, E. F., Evans, II, N. J., Bakker, E. J., & Blake,
G. A. 1999, ApJ, 522, 991

van Dishoeck, E. F., Blake, G. A., Jansen, D. J., & Groesbeck, T. D. 1995, ApJ,
447, 760

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Wolf, S., Schegerer, A., Beuther, H., Padgett, D. L., & Stapelfeldt, K. R. 2008,
ApJ, 674, L101

Wootten, A. 1989, ApJ, 337, 858

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718,
1062

Yang, Y.-L., Evans, II, N. J., Green, J. D., Dunham, M. M., & Jørgensen, J. K.
2017, ApJ, 835, 259

Zapata, L. A., Loinard, L., Rodríguez, L. F., et al. 2013, ApJ, 764, L14
# Appendix A: SED

After the dust density components were inserted into the model, the SED was used as a sanity check that our added dust components do not alter the SED too much compared to the Schöier et al. (2002) model. The envelope at $r > 10^3$ AU has temperatures of 10–30 K, which dominates the submm and mm range, while deviations in the warmer dust mass will alter the near-infrared SED in the 10–100 µm range. In Fig. A.1 the final SED after including all structures, can be seen. While the SED appears reasonable, it should be noted that SEDs generally have very degenerate solutions, and one cannot exclude other structures or multidimensional solutions, based on the SED alone. The SEDs of the best rotating toroid model and PP-disk model were almost indistinguishable, as the SED is dominated by the cold dust in the envelope.

# Appendix B: CO isotopologue spectra

While we did not attempt to fit the CO isotopologues spectral lines with our models, we provide the plots in Fig. A.2, for completeness. A box was used around both the observations and model emission to extract the same region around the filament and two protostars, in order to avoid some of the observed extended emission south of the source A, that we did not attempt to recreate in the models.

# Appendix C: Interfacing between RADMC-3D and LIME

A fundamental difference between RADMC-3D and LIME is the grid used. RADMC-3D uses a structured grid (curvilinear is chosen for this work), while LIME uses Delaunay triangulation to produce an unstructured grid. Custom source code changes were made to LIME (v. 1.5) to allow for two local logarithmically distributed grids around the center of each of the stars in order to resolve the disks, while also changing the default density weighting value of the node positions<fnref n="3" />, to allow the global uniform grid to resolve the dust filament, with fewer nodes, as the default setting resulted in many nodes being spread in the less interesting envelope around the disks and the filament. The RADMC-3D grid density, temperature and velocity values (in the case of a rotating toroid) are written as arrays to a C header file, and a nearest-neighbor search is done for each LIME node, which is then assigned the physical values of the nearest cell. Due to the size of the octree refined RADMC-3D grid (upward of 5 million grid cells), a brute-force search is ill-advised, as this quickly becomes a speed bottleneck. Instead, a k-d tree search algorithm<fnref n="4" /> (e.g., Saftly et al. 2014) was inserted into our custom LIME code. This k-d tree is used for quick lookup of the nearest neighbor, allowing us to have very large RADMC-3D grids without execution time issues. The exponential mode (fastest mode)<fnref n="5" /> as well as parallel computing was engaged in LIME to reduce execution time. However, due to the large number of nodes necessary to resolve the relevant disk and filament regions, the execution time is typically around two hours per calculation, that is, six hours to produce all spectral cubes for a single model, which compounded the issue of proper parameter sampling.

# Appendix D: PP-disk model *p* and *Ψ* parameters

After running the PP-disk models over a broad parameter space of *p*, *Ψ*, *Σ*₀ and *H*₀ for disk A and B (Table 1), it became clear that low values of *p* and high values of *Ψ* result in higher peak flux densities, both in combination and by themselves, compared to the opposite part of the parameter space. A low *p*-value ensures larger amounts of dust at larger midplane radii, while a higher flaring coefficient *Ψ* raises the dust vertically. This combination results in larger amounts of dust with relatively high temperatures. Thus, *p* and *Ψ* were fixed to 0.5 and 0.25, respectively, for both disks, but neither has been constrained as optimal parameters for the two disks. However, with these fixed values, we can reproduce the maximal peak flux densities of the two disks. Another reason for these fixed values is the added computational time in RADMC-3D for higher *p*-values, due to higher densities in the inner region, causing significant slowdown of the photon propagation calculation, as the photon can become trapped in high-density regions, taking hundreds or thousands of absorption and reemission events to escape the cell. Since the goal is to see if a PP-disk model type can satisfy the constraints, but not necessarily constrain the dust distribution, this approach is acceptable.

# Appendix E: Dust opacity in disk B

A dust opacity model could possibly accommodate the peak flux density of source B, with lower scale-heights, making the PP-disk model more plausible for source B. But introducing a full investigation into the dust opacity in the disks of IRAS 16293 A and B, in terms of grain size, compositions and morphology, was beyond the scope of this work and not feasible considering our computational limitations. We note that our bare grain opacity at 850 µm is consistent with those used by Tobin et al. (2013), who uses 3.5 cm² g⁻¹ at 850 µm, in order to fit the L1527 SED slope from 0.85 mm to 7 mm. In comparison, the bare-grain opacities that we used, which dominate the dust emission in the beam toward A and B, have a value of 3.53 cm² g⁻¹ at 850 µm. As such, while unconstrained, the dust opacities used in this work are consistent with the independent modeling of another Class 0 disk.

# Appendix F: Auxiliary CO isotopologue maps and dust temperature contours

Due to the absorption toward source A in the modeled CO isotopologue maps, when using both the ISM standard abundances (Table 4) and when dividing by a factor of five (Fig. 5), the standard ISM abundances were divided by 10 and 15 in two new model CO isotopologue maps (Fig. F.1). We notice the presence of absorption in source A and B in Fig. F.1, regardless of the CO isotopologue number density in the gas-phase, suggesting that optically thick dust is responsible for the absorption in the disk components.

<footnote n="3">³ https://lime.readthedocs.io</footnote>

<footnote n="4">⁴ The k-d tree was implemented using code modified from https://rosettacode.org/wiki/K-d_tree#C implemented in C.</footnote>

<footnote n="5">⁵ https://lime.readthedocs.io/en/latest/usermanual.html#command-line-options</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. A.1. SED of the best rotating toroid model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. A.2. Modeled CO isotopologue spectrum and observed CO isotopologue spectrum. Vertical dashed lines show the limits of the zeroth moment maps. The model CO isotopologues are from a rotating toroid model satisfying the continuum constraints with $L_A = 14 L_\odot$ and $L_B = 7 L_\odot$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. F.1.** Zeroth moment maps of CO isotopologue observations and synthetic gas line emission from PP-disk models. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam⁻¹ km s⁻¹. The RA and Dec offsets are relative to the phase center of the observations. Standard ISM abundances were divided by ten and 15, in the *middle and lower panel rows*, respectively.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

**Fig. F.2.** Peak flux density of disk B at 868 µm with a PP-disk model, with $L_B = 3 L_\odot$ and $L_B = 7 L_\odot$. The black dashed, horizontal line marks the peak continuum position of source B, while the blue area represents the uncertainty. We note that the disk mass reaches 0.35 $M_\odot$ and 0.5 $M_\odot$ at 7 g cm$^{-2}$ and 10 g cm$^{-2}$, respectively. The highest Class 0 disk mass found in a mm survey by Jørgensen et al. (2009) is ~0.46 $M_\odot$, with typical Class 0 disk masses of 0.05 $M_\odot$.

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

**Fig. F.3.** Temperature contours of a single star in the left panel with 21 $L_\odot$ in an envelope with a radial power-law density profile, $\rho_{0,\text{env}} = 3.1 \times 10^{-14} \text{ g cm}^{-3}$, $p_{\text{env}} = 1.7$, a central density plateau within 600 AU, and two protostars in the right panel, each with $10.5 L_\odot$. See Table 1 for more details.
<figure><img src="image_13.png" /><figcaption>**Fig. F.4.** Temperature contours for a rotating toroid around each source and a dust filament between the sources in a slice through the model at $z$ (height) = 0. Dust density contours are shown in a gray to green colormap. The radial power-law density distribution is visible along with the dust filament and rotating toroids. Here $L_A = 18 L_\odot$ and $L_B = 3 L_\odot$. The source A rotating toroid is edge-on while the source B rotating toroid is face-on. The 90 K contour around source B coincides with the transition from disk to dust filament. Due to the high number of grid cells, some temperature noise is visible in the plot center (due to very small cells with bad photon statistics in RADMC-3D). Interpolation has been performed along the temperature contours to make them appear more smooth for visual purposes, as the contours were otherwise slightly noisy in certain regions.</figcaption></figure>