Free Access
Issue
A&A
Volume 612, April 2018
Article Number A72
Number of page(s) 20
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201731668
Published online 27 April 2018

© ESO 2018

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; 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 ± 5 L (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 H2 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 13CO, C17 O and C18 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−1 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 αJ2000 = 16h32m22.s72; δJ 2000 = 24°2834.′′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−1, while source A has a more extended and weaker continuum, peaking at 1.0 Jy beam−1 (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 = 32 transitions of three CO isotopologues, 13CO (330.588 GHz, Klapper et al. 2000), C18O (329.331 GHz, Klapper et al. 2001) and C17O (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−1 (Jørgensen et al. 2011) was used for all lines. Channels within ±5 km s−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. C17 O appears to trace the dust distribution fairly well, as does the C18O emission. 13CO, however, shows a markedly different emission structure, which is likely attributable to optical depth effects, where only the outer layer of the 13CO structure is probed. Absorption toward source B is evident in all lines but to a smaller extent in C17O, revealing that C17O 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.).

thumbnail Fig. 1

Dust continuum at 868 μm, in log-scale. Contour levels are logarithmically divided between 0.02 and 2 Jy beam−1. The RA and Dec offsets are relative to the phase center of the observations. Beamsize is shown in bottom right corner.

Open with DEXTER

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 × 103 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)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).

Table 1

RADMC-3D model parameters.

thumbnail Fig. 2

Zeroth moment maps of the 13CO, C18 O, and C17 O isotopologues from the PILS observations. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. The RA and Dec offsets are relative to the phase center of the observations.

Open with DEXTER

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., rrplat, see Eq. (1)) was resolvedlinearly, while r > rplat (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.

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 106 cm−3. 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 temperaturewere 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 rplat, 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: (1)

Where ρ is dust density, r is the distance to the model center, ρ0 is the density at distance r0 and rplat is the radius of the density plateau. rplat 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, H2 CO and H2 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).

A walk through the parameter space of the dust reference density ρ0, with p = 1.7, was performed with rplat = 600 AU, L* = 21 L (if one source) and L1 = 10.5, 14.0 or 18.0 L, which means that, accordingly L2 would be 10.5, 7.0 or 3.0 L, in the case of two sources. The most important value is the central density plateau, ρplat, as the 90 K region, almost all of the 50 K region and the bulk of the 30 K region resides here (Fig. F.3). A given volume difference in two different ρ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 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, H2 CO and H2 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.

thumbnail 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, 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, (our IRAS 16293 envelope model has a mass of ~4 M). There is a clear tendency for a single star to result in more dust heated above ~ 50 K.

Open with DEXTER

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 MA   ~ 1 M, 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, 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 (CH3OCHO) and ketene (H2 CCO) as Keplerian rotation around a central object of 0.53 M 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, is dependent on the inclination and centrifugal barrier radius. The mass of 0.53 M 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 MA = 1.0 M 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 ≤ 0.1 M 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 μ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′′× 0.5′′ is comparatively large compared to the continuum emission at 868 μ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′′× 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 μm continuum emission in a 2′′× 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′′ × 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 NH2 > 1.2  × 1025 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 as measured from the tube center, where ρ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 ρ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 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 g cm−2, while the model column density was 9.7 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, Hd, is (2)

where r is the radius, H0 is the scale-height at radius r0 and Ψ is the flaring constant, controlling the vertical distribution of the dust above the disk midplane. The disk surface density, Σ, is (3)

where p is a gradient parameter. Finally the dust density is (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 g cm−3 ( ~ 3 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 μm continuum image, which, after convolution with a 0.5′′× 0.5′′ beam using MIRIAD, was compared with the 868 μ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 a χ2 analysis of alarge, 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.

Table 2

Example dust model parameters.

4.4 Model analysis

Starting with the PP-disk model, large steps were taken in the parameter space of Σ 0 and H0 for the disk around each source, and inclination for disk A (Table 1), while p and Ψ 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 (Σ, H0, i and source A luminosity, LA) and two for disk B (Σ and H0), since the source B luminosity, LB, was determined by LA.

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 H0 and i, only models with inclinations below 55° could be excluded. The peak flux values were degenerate with different Σ 0, H0, 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 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 andimpact opposing streamlines in the midplane, creating a disk structure. From Hartmann (2009), (5)

where is the mass accretion rate onto the central object, θ0 is the starting polar angle of the streamline in question and rc is the centrifugal radius. The density in each grid cell in our model was found by predefining nstream streamlines with discrete θ0 values and evaluating them at nr 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 ρ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 midplane, hence higher dust temperatures, than a PP-disk model with, for example, H0 = 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 (rc,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 rc,A = 75 AU, as it matched the aspect ratio at this inclination. The free parameters of the source A rotating toroid model were A and LA, as rc,A and MA were fixed (Table 1), while the source B rotating toroid model had B and rc,B as free parameters (since LB was determined by LA). 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 LA or LB, 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 =32 line emission from 13CO, C17O and C18 O of select dust models were modeled in LIME (Brinch & Hogerheijde 2010), a line radiative transfer code. The dust temperature output fromRADMC-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 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 , as H2 is by far the mostabundant 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) were downloaded from the LAMDA2 database (Schöier et al. 2005). The collisional rate coefficients of the CO isotopologues with H2 in these files include collisions with both ortho- and para-H2.

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 × 103 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 consistentwith 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, (6)

where r is the distance to the model center, G is the gravitational constant and Mtot 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 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.

thumbnail 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−1. DM denotes the PP-disk model, while RTM denotes the rotating toroid model.

Open with DEXTER

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 H2, 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 Tdust ≥ 30 K where the CO is sublimated from the dust grains, and a drop factor of 10−2 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.

Table 3

LIME model parameters.

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 LA   ≥ 18 L 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 C17O and C18 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. LA   ≥ 18 L 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 LA 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 LA 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 LA = 18 L was modeled in LIME, together with a representative of the best matching PP-disk model at LA = 14 L (since LB < 7 L could not reproduce the observed peak flux density in our H0 parameter space, though this does not rule out LB < 7 L 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 C18 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 C18 O emission structure was more reminiscent of the observations, while C17O 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 13CO bimodal emission structure above and below the disk midplane, as well as the two observed C18 O emission peaks near the corners of the disk A component, but with less absorption in C18 O toward source A than observed. The rotating toroid model cannot, however, reproduce the C17O centrally concentrated emission structure, but instead results in a bimodal distribution as seen in C18 O, near rc,A (rdisk = 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 C17O 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.

thumbnail 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−1 km s−1. 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.

Open with DEXTER
Table 4

LIME model parameters of Fig. 5.

thumbnail Fig. 6

Observed and modeled C17O zeroth moment maps, contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. All models are rotating toroid models matching the dust continuum peak flux densities for disk A and B.

Open with DEXTER
Table 5

LIME model parameters of Fig. 6.

6 Discussion

The combined constraints from the dust continuum and CO isotopologue emission, necessitate very high PP-disk model scale-heights around source B, H0 > 3.0 AU at 10 AU, which is suggestive of a nonsettled disk with a large vertical distribution, 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 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 H0 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 LB = 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 CO isotopologue abundances needed to be lowered by a factor of five to avoid self-absorption in C18 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 H2 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 Tdust ≥ 90 K than regions with 30 K ≤ Tdust < 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 Tdust = 300 K. The inferredluminosity 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. Ifsource 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 LA ~ 18 L and LB ~ 3 L.

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 = 32 line emission of 13CO, C18 O and C17O 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 3Dmodel 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 distribution 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 C17O 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 C18O 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; LA ~ 18 L and LB ~ 3 L.

This work illustrates the importance of using complex, full 3D radiative transfer modeling together with subarcsecond observationsof 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 quantitative investigation with parameter confidence intervals, to determine the exact CO isotopologue abundances and the parameters 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 interpretation 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.

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

thumbnail Fig. A.1

SED of the best rotating toroid model.

Open with DEXTER
thumbnail 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 LA = 14 L and LB = 7 L.

Open with DEXTER

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 positions3, 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 algorithm4 (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)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, Ψ, Σ 0 and H0 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 cm2 g−1 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 cm2 g−1 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.

thumbnail Fig. F.1

Zeroth moment maps of CO isotopologue observations and synthetic gas l ne emission from PP-disk models. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. 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.

Open with DEXTER
thumbnail Fig. F.2

Peak flux density of disk B at 868 μm with a PP-disk model, with LB = 3 L and LB = 7 L. 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 and 0.5 M 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, with typical Class 0 disk masses of 0.05 M.

Open with DEXTER
thumbnail Fig. F.3

Temperature contours of a single star in the left panel with 21 L in an envelope with a radial power-law density profile, ρ0,env = 3.1 × 10−14 g cm−3, penv = 1.7, a central density plateau within 600 AU, and two protostars in the right panel, each with 10.5 L. See Table 1 for more details.

Open with DEXTER
thumbnail 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 LA = 18 L and LB = 3 L. 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.

Open with DEXTER

References


4

The k-d tree was implemented using code modified from https://rosettacode.org/wiki/K-d_tree#C implemented in C.

All Tables

Table 1

RADMC-3D model parameters.

Table 2

Example dust model parameters.

Table 3

LIME model parameters.

Table 4

LIME model parameters of Fig. 5.

Table 5

LIME model parameters of Fig. 6.

All Figures

thumbnail Fig. 1

Dust continuum at 868 μm, in log-scale. Contour levels are logarithmically divided between 0.02 and 2 Jy beam−1. The RA and Dec offsets are relative to the phase center of the observations. Beamsize is shown in bottom right corner.

Open with DEXTER
In the text
thumbnail Fig. 2

Zeroth moment maps of the 13CO, C18 O, and C17 O isotopologues from the PILS observations. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. The RA and Dec offsets are relative to the phase center of the observations.

Open with DEXTER
In the text
thumbnail 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, 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, (our IRAS 16293 envelope model has a mass of ~4 M). There is a clear tendency for a single star to result in more dust heated above ~ 50 K.

Open with DEXTER
In the text
thumbnail 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−1. DM denotes the PP-disk model, while RTM denotes the rotating toroid model.

Open with DEXTER
In the text
thumbnail 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−1 km s−1. 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.

Open with DEXTER
In the text
thumbnail Fig. 6

Observed and modeled C17O zeroth moment maps, contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. All models are rotating toroid models matching the dust continuum peak flux densities for disk A and B.

Open with DEXTER
In the text
thumbnail Fig. A.1

SED of the best rotating toroid model.

Open with DEXTER
In the text
thumbnail 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 LA = 14 L and LB = 7 L.

Open with DEXTER
In the text
thumbnail Fig. F.1

Zeroth moment maps of CO isotopologue observations and synthetic gas l ne emission from PP-disk models. Contour levels are divided logarithmically from 0.5 to 7.3 Jy beam−1 km s−1. 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.

Open with DEXTER
In the text
thumbnail Fig. F.2

Peak flux density of disk B at 868 μm with a PP-disk model, with LB = 3 L and LB = 7 L. 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 and 0.5 M 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, with typical Class 0 disk masses of 0.05 M.

Open with DEXTER
In the text
thumbnail Fig. F.3

Temperature contours of a single star in the left panel with 21 L in an envelope with a radial power-law density profile, ρ0,env = 3.1 × 10−14 g cm−3, penv = 1.7, a central density plateau within 600 AU, and two protostars in the right panel, each with 10.5 L. See Table 1 for more details.

Open with DEXTER
In the text
thumbnail 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 LA = 18 L and LB = 3 L. 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.