© ESO, 2017

1. Introduction

A key issue in the early stages of star formation is the formation epoch of a rotationally supported disk, which mediates the accretion onto the forming star and evolves into a protoplanetary disk at later stages. Analytical and numerical studies have shown that the formation of a rotationally supported disk requires leaving enough angular momentum during the collapse. This can be a problem when a magnetic field is present because the field can efficiently transport angular momentum and prohibits the formation of a disk (e.g., Hennebelle & Fromang 2008). Earlier theoretical studies were mostly restricted to the simple case of a magnetic field parallel to the rotation axis, a configuration where magnetic braking can be very efficient and quench the formation of a disk. The B335 protostar seems to be the best case of such a geometry (Yen et al. 2015). The occurrence frequency of such systems is unknown, however. They may be the exception rather than the rule, since Lee et al. (2016) have shown that outflow axes of wide binary or multiple systems are not aligned in Perseus, suggesting that these systems are formed in more complex configurations. Different effects contribute to reducing the efficiency of magnetic braking and allowing early formation of disks: a rotation axis not aligned with the magnetic field (Ciardi & Hennebelle 2010), the level of turbulence (Joos et al. 2013), or deviations from ideal magnetohydrodynamics (MHD), especially the role of ambipolar diffusion (Tomida et al. 2015; Masson et al. 2016; Zhao et al. 2016; Hennebelle et al. 2016).

Disks have been commonly found around class I protostars and class II sources (Williams & Cieza 2011; Harsono et al. 2014), but their presence around younger objects, especially class 0 protostars, is still debated. Establishing the presence of rotationally supported disks in class 0 protostars (defined as protostars with a high ratio of the submillimeter to the bolometric luminosity; André et al. 1993) is especially important because objects at this stage are accreting material from their surrounding envelope at the highest rate. Therefore these objects are best suited for a comparison with theoretical collapse models, especially regarding the question of the complex interplay of velocity field and magnetic field on the final structure of the environment of the forming star. The number of such protostars that can be studied in detail is limited, and direct high angular resolution imaging in the near-infrared is often impossible because of the very high column densities of dust and gas in the surrounding envelopes. Only (sub)millimeter interferometers have the sensitivity and angular resolution to detect circumstellar material. Maury et al. (2010) performed a pilot study of multiplicity with the IRAM-PdBI at ~0.4′′ resolution that showed a low multiplicity fraction at small separations. At the achieved angular resolution, the class 0 protostars are detected, but the structure of the millimeter emission is not resolved. The young protostars B1b-S and B1b-N are promising objects for studying the early stages of low-mass star formation. These sources are located in the Barnard 1b core of the Barnard 1 dark cloud, a moderately active star-forming region in the Perseus molecular cloud at a distance of 230 pc (Hirota et al. 2008). The analysis of the spectral energy distribution (SED) based on Herschel and Spitzer data performed by Pezzuto et al. (2012) shows that both sources have low infrared luminosities (L_bol = 0.1 L_⊙ for B1b-N, L_bol = 0.3 L_⊙ for B1b-S), with the SED peaking longward of 100 μm. To further progress in the understanding of the nature of these sources and their relation to their parental core, Hirano & Liu (2014) obtained CO (J = 2 → 1), ¹³CO (J = 2 → 1) and H¹³CO⁺ (J = 1 → 0) observations with the SMA, which showed that both B1b-S and B1b-N are driving low-velocity molecular outflows. The continuum sources are largely unresolved down to an angular resolution of 0.5′′ at 7 mm wavelength (Hirano & Liu 2014). Huang & Hirano (2013) showed that N₂H⁺ (J = 3 → 2) and N₂D⁺ (J = 3 → 2) present a compact emission component associated with the continuum sources, together with extended emission in the surrounding envelope. The high intensity of these high-dipole moment species has been interpreted by Daniel et al. (2013) using a sophisticated radiative transfer code. As expected for a cold and dense region, the deuterium fractionation clearly increases with the density, although the absolute abundances slowly decrease. Using NOEMA with a 2.3′′ beam, Gerin et al. (2015) further confirmed the presence of slow molecular outflows driven by B1b-N and B1b-S, and determined the outflow dynamical times, masses, mechanical luminosities, and mass-loss rates. They concluded that the two sources are at most a few thousand years old, and that the properties of B1b-N are consistent with those of a first hydrostatic core (FHSC). The outflows are launched in different directions, which is probably a consequence of the high degree of turbulence in the Barnard 1b core, which is impacted by outflows from nearby young stellar objects (YSOs; Hirano & Liu 2014). So far, the inner structure of the compact components of B1b-N and B1b-S remained unresolved or only marginally resolved. Because these sources are young and the physical and chemical structure of the Barnard 1b core are well known (see the recent analysis of the depletion and ionization rate by Fuente et al. 2016), it is interesting to probe these sources in more detail, as a means to constrain current theoretical understanding of low-mass star formation.

2. ALMA observations

We present ALMA Cycle 3 observations, performed on 2015 November 22 during the Long Baseline Campaign, combined with supplementary Director’s Discretionay Time (DDT) observations obtained on 2016 June 16. For the first data set, a total of 46 antennas were used in the C36-7 configuration, with baselines from 82 m up to 15 km (projected on source), providing an angular resolution of ~0.06′′, corresponding to ~13 au at the distance of Perseus, 230 pc. Four continuum base-bands with 2 GHz bandwidth, 128 channels, and dual polarization received the digitized signals from the Band 7 receiver. Central sky frequencies in each band are 336.5 and 338.4 GHz in the lower sideband (LSB), and 348.5 and 350.5 GHz in the upper sideband (USB). The channel spacing is 15.6 MHz, resulting in a spectral resolution of ~27 km s^-1 (Hanning smoothed). Imaging with the first dataset revealed artifacts attributable to a lack of shorter spacings. For the second data set, a total of 38 antennas were used in the C40-4 configuration with baselines from 15 m to 704 m. The setups were otherwise identical.

Observations for both datasets were obtained in single sessions, providing 30 min on source combined. The two fields B1b-N and B1b-S were observed alternately for 55 s each, with a phase calibrator scan in between, observed for 20 s. The phase calibrator was J0336+3218, located at an angular separation of ~1.36° from the B1b cores. A second calibrator (J0359+3220, located at 5.7° distance), observed every 5 min, was also included in order to check the quality of the phase transfer (see Appendix A for more information). Standard calibration was performed in CASA 4.5.0, while for imaging and data analysis we used the GILDAS software and CASA 4.7.0. For the first dataset, two other quasars were used for flux (J0238+1636) and bandpass (J0237+2848) calibration. The adopted flux density of J0238+1636 is 1.06 Jy at 343.5 GHz. For the second dataset, Ceres was observed for flux calibration and J0237+2848 for bandpass calibration. J0238+1636 and J0237+2848 are both among the regularly monitored ALMA calibrators, and the obtained fluxes are consistent with the data in the ALMA Calibrator Database. For the DDT observations, the Ceres data could not be used for flux calibration and J0237+2848 (flux 1.10 Jy at 343.5 GHz) was used instead. Weather conditions were good and appropriate for Band 7 observations, with a precipitable water vapor of 0.8 (0.7) mm and a median system temperature of 170 (160) K for the long-baseline (DDT) sessions.

The flux calibration we achieved is very good, with the curves of flux versus uv distance overlapping for the two program sources, as illustrated in Fig. A.2 (although the gain calibrator’s flux varied somewhat). Extended flux is clearly present for both sources, as the recovered flux increases with decreasing baseline length. Spectral windows for each sideband were averaged for the combined data.

For the high-resolution data, the final clean beam is 0.064′′ × 0.048′′ at PA 168° with a noise level of 0.1 mJy/beam for each spectral window, corresponding to about 0.35 K in the LSB and 0.32 K in the USB. For the DDT data, the final clean beam is 0.39′′ × 0.27′′ at PA − 9.9°, with a noise level of 0.2 mJy/beam, corresponding to about 18 mK.

We have checked that all spectral windows are free of strong line emission. With the typical line widths in Barnard 1b of ~1 km s^-1, which can increase to ~5 km s^-1 when line wings are present, line emission will be diluted in the broad spectral channels used for continuum observations. Although we cannot exclude some line contamination, we are confident that the contribution of spectral lines to the detected continuum fluxes remains modest.

For the data reported here, we have concatenated the two datasets. The final clean beams are 0.086′′ × 0.069′′ for B1b-N, and 0.14′′ × 0.11′′ for B1b-S. We averaged the two spectral windows of each sideband to obtain the final images.

3. Results

3.1. Resolving the B1b-N and B1b-S continuum emission

The ALMA images are displayed in Fig.1, while Table 1 presents a summary of the observations. For each source we show as gray contours the 32.9 GHz emission mapped with a 0.2′′ beam with the Very Large Array (VLA) as a part of the VLA Nascent Disk and Multiplicity Survey of Perseus Protostars (VANDAM; Tobin et al. 2015, 2016). ALMA clearly resolves a compact slightly elongated structure toward both sources that is located within a more extended region of faint emission.

We list in Table 1 the values of flux densities at the peak and within a 1′′ radius (230 au). While both sources reach similar intensities at large scale, the compact component in the youngest source B1b-N is smaller and represents a smaller fraction of the total flux than for the slightly more evolved source B1b-S. The high flux densities of the compact components toward B1b-N and B1b-S, >10 mJy/beam, correspond to a brightness temperature of ~35 K, comparable to the expected gas and dust temperatures in these sources (Commerçon et al. 2012a). It is therefore likely that the continuum emission of the compact component is optically thick.

To obtain further insight into the spatial distribution of the emission, we fit the flux density isocontours with ellipses centered on the protostar positions, with the position angle (PA) and ellipticity as free parameters. These fits are illustrated in Fig. 2 with B1b-N data in blue and B1b-S data in red. We also used simple Gaussian profiles to determine the source sizes and shapes, and compared these estimates with the core sizes obtained by fitting the circularly averaged intensity profiles as described in Appendix B. Table 2 lists the sizes and PAs of the central sources. We obtain a FWHM of about 0.4′′, corresponding to ~90 au at the Perseus distance. The sizes of the compact component for both sources correspond to the size of the centimeter continuum emission, as illustrated in Fig. 1. The similarity between the submillimeter and centimeter images, and the low spectral index between 32 and 350 GHz of 2.4 and 2.8 for B1b-N and B1b-S, respectively, confirm the high opacity of the 350 GHz emission. It is remarkable that B1b-N is weaker than B1b-S at wavelengths shorter than about 3 mm, but becomes stronger for longer wavelengths (Gerin et al. 2015). This behavior is explained by a larger part of the total flux coming from the compact source at long wavelengths and the shallower spectral index of this compact component as compared with that of the extended envelope. Such a shallow spectral index down to at least 32 GHz is consistent with either very high opacities even at 32 GHz or the presence of large grains with a significantly higher emissivity at 32 GHz than standard interstellar grains. Studies of the submillimeter (Chen et al. 2016) and mid-infrared emission (Lefèvre et al. 2014) have shown that large dust grains are present in Barnard 1b, as testified by the high portion of “transitional β”¹ values, and the coreshine phenomenon. With a standard grain population, the dust opacity at 350 GHz is a factor 50–100 higher than at 32 GHz. High opacity at 32 GHz therefore requires extreme column densities, of N(H₂) ~ 10²⁷ cm^-2, and mean H₂ densities higher than ~10¹² cm^-3 at the 0.2′′ scale. We therefore propose that both phenomena contribute to the flattening of the spectral index, the high opacity of the continuum emission at 350 GHz and the presence of large dust grains.

3.2. Analysis of intensity, density, and temperature radial distributions

Figure 3 presents the observed distributions of the circularly averaged specific intensity as a function of the angular radius for B1b-N and B1b-S and their associated fits in the top panel, together with the variations of dust temperature, molecular hydrogen density, and the cumulative gas masses in the other panels. While the dust continuum emission from the compact component (disk) is probably optically thick, the sharp decrease in flux densities outside the central component indicates that the opacity of the dust thermal emission drops. The assumption of optically thin emission is therefore likely valid for the envelope region. Using the procedure outlined in Appendix B, we derived the radial distribution of the gas density. This assumes spherical symmetry, that the gas is molecular, and that the dust has constant optical properties, cm² g^-1, a value commonly used for cold and dense cores, from Ossenkopf & Henning (1994), which is close to the composite aggregate disk opacities computed by Semenov et al. (2003). A further assumption is the good thermal coupling of the dust and gas, which leads to equal gas and dust temperatures, supported by the high densities encountered in the cores (n(H₂) > 10⁵ cm^-3, Daniel et al. 2013). Finally, we used a simple equation of state to relate the dust temperature and the gas density: i) the gas is isothermal and the kinetic temperature T₀ is 12 K (derived by Lis et al. 2010, from NH₃) for all densities lower than a density threshold n₀(H₂) = 10¹⁰ cm^-3; ii) the gas follows a polytropic equation of state for densities higher than n₀(H₂), with an exponent γ, implying . We have chosen γ = 1.6, adapted for moderately warm molecular gas. For the “disk region” within a radius of 0.3′′, we derived the dust temperature directly from the brightness temperatures. As the 350 GHz emission is optically thick, we used the 32 GHz flux densities to determine the gas masses, assuming an emissivity of κ = 0.006 cm² g^-1 at 1 cm, adapted to a population of large grains with sizes reaching 3 mm (Tazzari et al. 2016), and using the temperature determined from the 350 GHz data. The total masses were derived by combining the determinations at 32 GHz and at 350 GHz.

The parameters used for the fit and the resulting physical conditions are given in Table 2. For both sources, we used a two-component fit with a central component with a steep flux density profile (r^-4 or r^-6) and a more extended envelope with flux decreasing as r^-2. Because of the lack of short spacings, some flux may be missing at scales larger than ~2′′. Therefore the parameters of the second component are not well constrained. The exponent and size scale should not be considered as having a physical meaning, but only as a convenient mathematical description of the observed data. Assuming optically thin 350 GHz emission, the maximum H₂ density reached in the central region is about 10¹⁰ cm^-3 for both sources. The true densities are expected to be significantly higher. The peak dust temperature is 40 and 60 K, somewhat higher than the mean temperature used for the SED fit by Hirano & Liu (2014, 16 and 18 K for B1b-N and B1b-S), but still consistent with faint emission at wavelengths shorter than 70 μm. For both sources, the warm region where the dust temperature is higher than 12 K is fairly small, less than 0.3′′ in radius. The cumulative mass distributions show that the sources may still be accreting, since a large part of the core mass is located at radii larger than 0.3′′. B1b-S is more condensed and presents a steeper increase in mass with radius than B1b-N, supporting its more evolved classification (established from the larger extent of its molecular outflow and warmer SED; Gerin et al. 2015). The mass of the compact component is at least 0.1 M_⊙ for both B1b-N and B1b-S, given our hypothesis on the dust grain emissivity at 32 GHz.

To check whether the compact component detected with ALMA could simply represent the inner region of a protostellar envelope, we compared the cumulative flux distribution of a simple isothermal spherical envelope model with that of B1b-N and B1b-S. We used a density profile , where r is the radius in arcseconds, a uniform temperature of 12 K, and the dust emissivity law described above. These parameters lead to a predicted flux in a 20′′ beam of about 2 Jy, comparable to the measured flux density of B1b-N and B1b-S with SCUBA, about 3 Jy (Pezzuto et al. 2012; Hirano & Liu 2014; I-Hsiu Li et al. 2017), and a high opacity of the central region within a radius of 0.1 arcsec, as illustrated in Fig. 4.

However, this model presents a cumulative flux distribution very different from that of B1b-N and B1b-S, which have a more spatially extended region of high opacity and a higher temperature, supporting the association of the compact component detected with ALMA with a disk. To further understand the nature of this compact component, we have used numerical simulations of collapsing cores.

4. 3D collapse model and comparison with the ALMA data

4.1. Description of the model

In this section, we compare our ALMA data with simulated observations obtained by post-processing a 3D numerical simulation from Hennebelle et al. (2016). We have used the RAMSES code (Teyssier 2002), which integrates the MHD equations using a constrained transport scheme (Fromang et al. 2006; Teyssier et al. 2006) and accounts for the ambipolar diffusion (Masson et al. 2012). The model is almost identical to the model used in Masson et al. (2016), except that we account for initial turbulent velocity fluctuations instead of a coherent global solid body rotation (e.g., Joos et al. 2013). The model consists of a collapsing 1 M_⊙ magnetized and turbulent dense core of radius R₀ = 2500 au. The initial dense core density and temperature are uniform, ρ₀ = 9.4 × 10^-18 g cm^-3 corresponding to n(H₂) = 2.4 × 10⁶ cm^-3 and T = 10 K. The corresponding ratio between the thermal energy and the gravitational energy is α_th = 0.25. The initial mass-to-flux ratio parameter is μ = 2 (strongly magnetized model). The initial turbulent velocity fluctuations are introduced following Joos et al. (2012, see their Sect. 2.2). The velocity fluctuations are scaled to match an initial Mach number ℳ = 1.2. The thermal behavior of the collapsing material is reproduced using a barotropic equation of state. The initial resolution is 32³ and the mesh is refined to always describe the Jeans length with at least eight points, down to a resolution of ~0.6 au.

Figure 5 shows the protostellar disk structure in the models 6.3 kyr after the FHSC formation. The rotationally supported disk is identified following the criteria derived in Joos et al. (2012). The disk extends up to ~75 au, and 98% of its mass is contained within a radius <60 au. The FHSC mass is 0.09 M_⊙ and the disk mass 0.07 M_⊙. Figure 6 shows an azimuthal average of the density and velocity fields around the rotation axis. The outflow extends up to ~1100 au.

4.2. Synthetic observations

We post-processed the simulation using the RADMC-3D radiative transfer code (Dullemond 2012) to compute the expected 350 GHz emission. Following Commerçon et al. (2012a,b), we assumed that the gas and dust temperature are perfectly coupled and that the dust-to-gas ratio is uniform and equal to 1%. We used the dust opacity from Semenov et al. (2003) for the homogeneous spheres model (with Fe/Fe + Mg = 0.3 “normal” silicate composition). We used the Barnard 1 distance and performed the computation for a range of inclinations i between face-on to edge-on. The emission structure is not very sensitive to the inclination between about 20° and 60°.

4.3. Comparison of models and observations

Figure 7 presents the comparison of the radially averaged flux toward B1b-N and B1b-S with the predictions obtained from the model at two epochs. We used the views with i = 20° for two different epochs: an early stage of the collapse, 400 yr before the formation of the FHSC, when the gas structure remains close to that of an isothermal sphere and no disk is formed yet, and a late stage (6.3 × 10³ yr after the formation of the FHSC) when a rotationally supported disk has already formed, as illustrated in Figs. 5 and 6. The emission structure of the late-stage model is similar to that of the B1b protostars: both the disk radius and the predicted flux density radial distribution agree within better than a factor of two with the observations. The early-stage model presents a very different structure, with a more localized emission peak of about 0.04′′, and a rapidly decreasing flux significantly lower than that detected toward B1b-N and B1b-S, as expected from the similarity between this early stage and the simple model shown in Fig. 4. The difference of flux radial distributions between the two epochs stems from the different temperature and density distributions. In particular, the high surface brightness region is more localized when no disk is formed (0.04′′ vs. 0.2′′) because the very high density region (n(H₂) > 10⁹ cm^-3) is reduced. The excellent agreement of the observed 350 GHz emission profile with the model profile confirms the high opacity of the dust emission and indicates that the dust temperature remains moderate.

To further explore the similarity of the model with the B1b-N data, we now consider the outflow spatial extend and maximum velocity, which are illustrated in Fig. 6. The values predicted by the model are fairly similar to those of B1b-N (Gerin et al. 2015). The good correspondence of the central sources’ size and shape, combined with their detection at 32 GHz, the extremely high molecular hydrogen densities and their estimated masses, support the identification of these sources with rotationally supported disks. This identification remains tentative in the absence of kinematic data, however. The sizes of the tentative disks are within a factor of a few from the predictions using the analytical expression derived by Hennebelle et al. (2016). Although fairly uncertain because we lack information about the grain emissivity law, especially when large grains are present (Dunham et al. 2014), the disk masses and temperatures are consistent with theoretical predictions of low-mass star collapse models (Masunaga & Inutsuka 2000; Machida et al. 2010; Dunham et al. 2014; Masson et al. 2016). Zhao et al. (2016) have shown that the removal of small grains can favor the formation of disks by reducing the coupling of the matter with the magnetic field, hence the efficiency of magnetic braking. The presence of large dust grains in Barnard 1b argues for a lower fraction of small dust grains in this core, which may have favored the formation of disks. We conclude that although no kinematic information is available so far, the sizes (~50 au) and brightnesses are consistent with those of rotationally supported disks formed during the collapse of dense cores in current MHD simulations. B1b-N and B1b-S thus appear to be excellent objects for testing our understanding of disk formation during low-mass protostar formation, which deserve further studies of their kinematics and spatial structure.

---

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00025.S and ADS/JAO.ALMA#2015.1.00025.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), 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. We acknowledge funding support from PCMI-CNRS/INSU, PNPS-CNRS/INSU, AS ALMA, Observatoire de Paris, ANR project ANR-13-BS05-0008 (IMOLABS) and the European Research Council (ERC Grant 610256: NaNOCOSMOS). We thank IRAM ARC for help with the data processing, and J. Tobin for access to the VANDAM images. We thank the referee for a careful reading of the manuscript which led to a significant improvement of our analysis. This work was performed using HPC resources from GENCI-CINES (Grant 2016-047247).

References

Appendix A: Calibration

In order to check the quality of the phase-calibration transfer to the source targets for the long-baseline dataset, we first imaged the control source (J0359+3220) using the phase solution obtained on the phase calibrator (J0336+3218). However, the solution was only moderately satisfying because J0359+3220 is relatively far from the phase calibrator (5.7°), and the quasar appears slightly extended. In such conditions, the imaging results do not reflect the real performances on the targets. Therefore we performed a second test on the phase calibrator J0336+3218. In these tests, we used half of the scans observed on the phase calibrator to obtain the calibration tables, which then were used to transfer the phase calibration onto the other half of scans. For the second half, we therefore performed phase transfer instead of self-calibration. The images are acceptable for the phase-transferred scans, which shows that we can rely on the phase calibration of the targets. We also performed the cleaning in GILDAS in the same way as for the B1b sources. The resulting images are displayed in Fig. A.1.

Figure A.2 presents the amplitude versus UV distance for the combined data sets. The good overlap of the curves of flux versus UV distance in the domain from 20 m up to 600 m shows that the flux calibration is very good. B1b-N and B1b-S are heavily resolved, with most of their fluxes on short baselines.

Appendix B: Inversion of radial profile

As shown by Krčo & Goldsmith (2016), it is possible to use the radial profile of the projected 2D image of a spherically symmetrical source to retrieve the information on the initial 3D radial distribution. This appendix describes the method we applied to the submillimeter images. Assuming a spherical source of radius R_max, at a distance D, the flux density can be expressed as the integrated emission along the line of sight: (B.1)where n is the dust density, κ_ν is the dust emissivity, and Ω the beam size. The total flux density is then (B.2)Since the density and temperature only depend on the radius, it is possible to inverse the radial distribution of the flux density to derive the product of the density and Planck function n(r)B_ν(T(r)). Inversion methods are rather sensitive to the noise level. We therefore chose another approach and fit the radial distributions of the flux density with an analytical radial profiles. We used the analytical solution of the inversion to determine n(r)B_ν(T(r)). The two variables can be separated in a following step, using a simple assumption for the temperature profile as described in Sect. 3.2. We used a simple attenuated power-law profile for the radial distribution of the flux density (see, e.g., Krčo & Goldsmith 2016), (B.3)which can be easily inverted to yield (B.4)where Γ(x) = (x − 1) ! with x an integer number ≥2.

Developing B_ν(T(r)) yields in the Rayleigh-Jeans limit (B.5)We note that at 350 GHz, the Rayleigh-Jeans limit is only valid for T ≥ 35 K, hence it is better to use the actual Planck function. The fitting parameters for both sources are displayed in Table 2.

All Tables

All Figures