© V. Cabedo et al. 2021

1 Introduction

Low-mass stars are known to form in dense molecular gas clouds. Class 0 objects represent the first stage of the star formation process, when most of the mass is still contained in the envelope surrounding the protostar (André et al. 1993; André 1995). Models of protostellar collapse (Shu et al. 1987) suggest that it is during this phase that the circumstellar gas is transported to the central object thanks to accretion processes. During this stage, angular momentum needs to be removed from the envelope and stored in the central object or dissipated through viscous processes to allow for the formation of the star. Moreover, the accretion mode, the rate at which it happens, and the duration of possible accretion episodes during this phase will determine the final stellar mass (André 1995; Basu & Jones 2004; Bate & Bonnell 2005; Myers 2012). Therefore, studying this phase is crucial, as it allows us to understand what the kinematics and dynamics of the gas are at the onset of collapse and to determine how those affect the outcome of the star formation process.

The gas making up the majority of protostellar envelopes is typically probed using molecular gas line profiles, which trace the gas kinematics in the dense envelope and measure gas motions such as rotation or infall. Observations of the molecular line emission from embedded protostars have been reported to trace widespread infall signatures in the inner ~2000 au of some protostellar envelopes (Zhou et al. 1993; Rawlings 1996; Di Francesco et al. 2001; Mottram et al. 2013). Most of those studies rely on the detection and interpretation of the infall spectral signature known as blue asymmetry or inverse P Cygni profile. In the center of a dense protostellar core, where gas densities become larger as the collapse proceeds (>10⁵ cm⁻³, André 1995), the large optical depth of some molecular line emission produces a self-absorbed line profile, with the dip centered on the emission at systemic velocity, where most of the circumstellar gas emits. Because the core is collapsing under the effect of gravity, line emission from dense gas tracers will present the typical blue-asymmetry line profile. Most works that are focused on modeling these line profiles to put constraints on protostellar infall models rely on the assumption of a symmetrical cloud collapsing, which produces such a double-peaked profile on optically thick emission. From the models, key parameters can be extracted, such as central protostellar mass, infall velocities, and mass accretion rates (Zhou et al. 1993; Di Francesco et al. 2001; Evans et al. 2005, 2015). However, blue-asymmetries in line profiles are not unique to infall motions. Complex gas kinematics, such as asymmetric collapse (Tokuda et al. 2014; Maureira et al. 2017), accretion streamers (Pineda et al. 2020), or outflow-entrained gas, can produce separated velocity components along the same line-of-sight (LOS) that are observed as double-peaked line profiles that are not caused by optical thickness.

The isolated Bok globule B335, which contains an embedded Class 0 protostar (Keene et al. 1983), is located at a distance of 164.5 pc, (Watson 2020) and has been the prototypical object for testing symmetrical collapse infall models, since blue-asymmetries were first detected in the molecular emission of the source at core scales (Zhou et al. 1993; Choi et al. 1995; Evans et al. 2005). Double-peaked line profiles have also been observed with interferometric observations of the molecular emission from the inner envelope (Chandler & Sargent 1993; Saito et al. 1999; Yen et al. 2010; Kurono et al. 2013; Evans et al. 2015). Attempts have been made for infall models of an optically thick line emission to compute infall mass rates (Yen et al. 2010, 2015; Evans et al. 2015), obtaining values affected by large uncertainties ranging from 10⁻⁷ M_⊙ yr⁻¹ to ~3 × 10⁻⁶ M_⊙ yr⁻¹ at radii of 100–2000 au, as well as infall velocities from 1.5 km s⁻¹ to ≈0.8 km s⁻¹ at radii of ~100 au. New models based on continuum emission, and using the revised distance of 164.5 pc, have determined an infall rate from the envelope to the disk of 6.2 × 10⁻⁶ M_⊙ yr⁻¹ (Evans et al., in prep.). For the estimated age of 4 × 10⁴ yr, this implies a total mass at the center (star + disk) of 0.26 M_⊙. However, the physical cause of these double-peaked line profiles has been debated. For example, Kurono et al. (2013) pointed out that despite expecting the H¹³CO⁺ emission to be optically thin, the inverse P Cygni profile and the position-velocity diagram they observe can be reproduced with models of moderately optically thick infalling gas. B335 is associated with an east-west outflow that is prominently detected in ¹²CO with an inclination of 10° on the plane of the sky and an opening angle of 45° (Hirano et al. 1988; Hirano et al. 1992; Yen et al. 2010). The eastern lobe is slightly oriented on the near side (Stutz et al. 2008), producing blueshifted emission on the eastern side and redshifted emission on the western side. While the core has been found to be slowly rotating at large scales (>2500 au) (Frerking et al. 1987; Saito et al. 1999; Yen et al. 2010, 2011), no clear rotation was found at smaller radii (<1000 au) and no kinematic signature of a disk was reported down to ~10 au (Yen et al. 2015, 2018). Recent observations of the hour-glass shaped magnetic field at small scales have suggested that B335 is an excellent candidate for magnetically regulated collapse (Maury et al. 2018).

In this work, observations of the molecular lines C¹⁷O (1–0) and C¹⁸O (1–0), which trace the dense circumstellar gas of the inner envelope of B335, are presented along with the ¹²CO (2–1) line emission tracing the outflow cavity. The molecular line profiles are analyzed and interpreted, providing new constraints on the gas kinematics close to the protostar.

2 Observations and data reduction

Observations of the Class 0 protostellar object B335 were carried out with the ALMA interferometer during the cycle 4 observation period from October 2016 to September 2017, as part of the 2016.1.01552.S project (PI: A. Maury). In the entirety of the work, it is assumed that the centroid position of B335 is at α = 19:37:00.9 and δ = +07:34:09.6 in J2000 coordinates, corresponding to the peak of dust continuum obtained from high-resolution maps (Maury et al. 2018). All lines were targeted using a combination of ALMA configurations: C¹⁷O (1–0) and C¹⁸O (1–0) were targeted using two configurations, C40-2 and C40-5, and ¹²CO (2–1) was targeted using C40-1 and C40-4. Technical details of the observations are shown in Table A.1.

A preliminary analysis of the data was done with the product images delivered by ALMA to check if any emission was detected and to check the shape of the line profiles. The C¹⁷O emission was only detected in the most compact configuration (C40-2), therefore, it is only this configuration that was used to produce C¹⁷O and C¹⁸O maps, while ¹²CO was detected in both configurations; thus, a combination of the two data sets was used to produce the maps. Calibration of the raw data was done using the standard script for cycle-4 ALMA data using the Common Astronomy Software Applications (CASA) version 5.6.1–8. The continuum emission was self-calibrated with CASA. Line emission was calibrated using the self-calibrated model derived from the continuum data when it was possible. Final images of the data were generated from the calibrated visibilities using the tCLEAN algorithm within CASA, using Briggs weighting with robust parameter set to 2 for C¹⁷O and C¹⁸O and 1 for ¹²CO. After imaging, the ¹²CO maps were smoothed to reach the same angular resolution as C¹⁷O and C¹⁸O. The resulting map characteristics are shown in Table 1.

3 Results and analysis

Figure 1 shows the moment 0 of the C¹⁷O (1–0) emission in red contours. The mean radius of the 3σ emission is 860 au, indicating that it traces the dense envelope. Top image shows in intensity the dust continuum map observed at 110 GHz, showing emission over 3σ. The bottom image shows in intensity themoment 0 of ¹²CO (2–1) emission. The ¹²CO emission probes the outflowing gas, therefore, it can be confirmed that the C¹⁷O emission is tracing the envelope and it is not affected by the outflow.

In order to understand the dynamics of the gas that is being probed, the C¹⁷O spectra were taken at every 0.5′′ pixel of the emission cube, producing a spectral map, as shown in Fig. 2. The line profile patterns show two distinct velocity components with a dip centered around the systemic velocity (8.3 km s⁻¹). Their respective intensities vary depending on the direction of the offset from the continuum peak, with the blue component shown to be more intense in the eastern part of the core, while the red component is dominant in the western part. This behavior is true for all the hyperfine components detected. A clear broadening of the line can be observed near the dust continuum emission peak, and in the north-east region part of the core. The former can be due to natural thermal broadening of the two velocity components as the temperature rises in the center of the object, while the latter might be the consequence of the overlapping of the two components due to other dynamical processes. The possibility of the dip being caused by interferometric filtering is discarded since the recovered emission size is on the order of the largest recoverable scale (LRS, see Table 1). Moreover, because C¹⁷O is a rare isotopologue, it is not expected to be abundant at largest scales and, therefore, no emission can be filtered at the systemic velocity. Figure 3 shows the velocity channel maps of the C¹⁷O (1–0) emission, covering a velocity range from 7.7 to 9.1 km s⁻¹ and only showing the main hyperfine component (Fig. B.1 shows a range covering from 5.0 to 9.6 km s⁻¹, showing the two observed hyperfine components). It can be seen that the blue- and red-shifted emission are confined to the east and west regions, respectively, suggesting that the two components are probing gas with different dynamics.

Three independent methods have been used to investigate the origins of the C¹⁷O spectral profiles and to examine whether the two distinct peaks could be produced by an optically thick line emission. In the following sections, the analysis of these methods and the modeling of the velocity field in the B335 inner envelope are presented.

Fig. 1 Red contours showing moment 0 of C17O emission, integrated over the velocity range 4.8–6.2 km s−1 and 7.6–9.4 km s−1. Contours show emission at −3, 3, 5, 10 and 30σ, where σ is 10.0 mJy beam−1. Top: intensity shows dust continuum emission map at 110 GHz for emission over 3σ, where σ is 8.56 × 10−2 mJy beam−1. Bottom: intensity shows 12CO (2–1) moment 0 emission integrated over the velocity range 1.4–16.2 km s−1, for emission over 3σ, where σ is 2016.9 mJy beam−1. The two arrows show the direction of the E-W outflow.

3.1 Line opacity estimation

The maximum opacity at the center of the source has been estimated from the H₂ column density, and assuming a standard C¹⁷O abundance using Eq. (1) (Jansen 1995): (1)

where = 6.695 × 10⁻⁸ s⁻¹ is the Einstein coefficient for the C¹⁷ O J = 1–0 transition (Gordon et al. 2017; Müller et al. 2005), c is the speed of light, ν is the frequency of this transition, = 3.1 × 10²² cm⁻² is the peak column density of H₂ in B335 at a radius of 3600 au (Launhardt et al. 2013), [C¹⁷O] = 5 × 10⁻⁸ is the C¹⁷O abundance relative to H₂ abundance (Thomas & Fuller 2008), and ΔV ≈ 1 km s⁻¹ is the average observed linewidth for the two components together. Using these values, the obtained opacity is τ₀ = 0.77, which corresponds to an opacity typical from an optically thin line.

3.2 Intensity ratio

Because of their similar mass and molecular structure, C¹⁷O and C¹⁸O are expected to probe gas under similar physical conditions. The [C¹⁷O]/[C¹⁸O] isotope ratio does not appear to be affected by fractionation and, in addition, if the emission from both C¹⁷O and C¹⁸O is associated with dense gas shielded from external ultraviolet radiation, selective photo-dissociation is not likely to affect the relative abundances (van Dishoeck & Black 1988). Thus, the only difference in the emission from these two isotopes would result from opacity effects because C¹⁸O is a factor of 3.6–3.9 more abundant than C¹⁷O (Penzias 1981; Jørgensen et al. 2002). The ratio of integrated intensities is much less sensitive to linewidth effects than the ratio of peak intensity; therefore, we use the latter to rule out any abundances effect on the C¹⁷O emission.

We produced beam-matching maps for C¹⁷O and C¹⁸O to compare the gas at similar scales in both isotopes. The obtained synthesized beams are given in Table 1. Intensities are integrated over the two velocity ranges of 4.8–6.2 km s⁻¹ and 7.4–9.4 km s⁻¹ for C¹⁷O, taking into account the two observed hyperfine components and 7.4–9.4 km s⁻¹ for C¹⁸O. The fact that both lines emit in the same range of velocities and are spatially coincident indicates that they are probing the same reservoir of circumstellar gas.

The integrated intensity ratio was computed as W/W, where W_i is the integrated intensity of each C¹⁷O and C¹⁸O individual spectra. Figure 4 shows the obtained integrated intensity ratio map, with values ranging from 0.15 to 0.50 and a mean value centered at 0.28. We note that the intensity ratio is quite homogeneous in most of the extension of the emission, but gets larger at the north and south-west regions. This is attributed to a line broadening of the C¹⁷O in those regions when compared tothe C¹⁸O emission (see C¹⁸O spectral map, Fig. C.1). The origin of this broadening is unknown, but we attribute it to complex dynamical processes that might be taking place in these regions.

The expected ratio if both transitions are optically thin is computed as [C¹⁷O]/[C¹⁸O] ≈ 0.25, where [C¹⁷O] and [C¹⁸O] are respectively the C¹⁷O and C¹⁸O abundance with respect to H₂ ([C¹⁷O] = 5 × 10⁻⁸ and [C¹⁸O] = 2 × 10⁻⁷, Thomas & Fuller 2008). Our observations are therefore in general agreement with the expected ratio if both transitions are optically thin. No increase of the intensity ratio is seen on the observations towards the center of the source, where opacity is expected to be higher, further confirming the later hypothesis. We note that this is also in agreement with the conclusions reached by the analysis of single-dish observations of the C¹⁸O (2−1) and C¹⁷O (2−1) (, Evans et al. 2005).

Fig. 2 Spectral map of the C17O (1–0) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel correspond to 0.5″ (~82 au). For clarity, the spectral range only shows the main hyperfine component. The green spectrum refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1).

3.3 Modeling of the molecular line profiles

The spectrum at each pixel has been modeled using a program that allows to fit the hyperfine structure of spectral lines with multiple velocity components (HfS, Estalella 2017). It also allows to compute the opacity of the line from the relative intensity of the different hyperfine components of the given transition. For every velocity component, the general HfS procedure fits four independent parameters simultaneously: the linewidth assumed to be the same for each hyperfine component; the main line central velocity; the main line peak intensity; and the optical depth. The fitting procedure samples the space parameters to find the minimum value of the fit residual χ². Because the C¹⁷O (1–0) emission is expected to be optically thin, we attempted to model the double-peaked profiles with two velocity components: one blue- and another red-shifted. Initial expected values were introduced and the program was allowed to proceed to fit emission with a minimum signal-to-noise ratio (S/N) of 3σ. The peak velocity, the linewidth, and the opacity maps obtained from the fitting are shown in Fig. 5.

The velocity maps (top images of Fig. 5) show that the two different components, blue and red-shifted, occupy two separated regions, at east and west offsets from the center of the source respectively, with some regions overlapping, where the double peakcan be observed in the spectra. The mean average, velocity for the two components are 8.1 and 8.6 km s⁻¹ respectively. The velocitydispersion maps (middle images of Fig. 5) show a mean velocity dispersion for the two components of 0.7 and 0.5 km s⁻¹, which get broader, up to 0.9 ~ 1, closer to thecenter of the object where the two components overlap (see central spectra in Fig. 2). The opacity maps (bottom images of Fig. 5) show that the opacity is generally less than 1, and only goes up to 3 in some specific pixels. These larger values are also associated with a very large error (on the order of the value itself) so, they are not significant and are not shown in the plots. Figure 6 shows the histograms for the opacity values for both fitted components. The average opacity has been estimated from the HfS modeling and is found to be τ_bs = 0.464 and τ_rs = 0.474, for the blue- and red-shifted components, respectively, which is in concordance with the upper limit estimated before.

Fig. 3 Contour channel maps of the emission of the F = 5/2–(−5/2) hyperfine component of the C17O J = 1–0 transition,overlapped with the gray scale image of the 1.3 mm dust emission. Contours are − 4, − 2, and from 2 to 18 by steps of 2σ, where σ is 10 mJy beam−1. The synthesized beam is shown in the bottom left panel as a filled ellipse. The vLSR channel velocities are shown in the top left part of the panels.

Fig. 4 C17O to C18O integrated intensity ratio. Red contours show integrated intensity for C17O at −3, 3, 5, 10 and 30σ, where σ is 10.0 mJy beam−1. Black contours show integrated intensity for C18O at −3, 3, 5, 10 and 30σ, where σ is 10.7 mJy beam−1. The yellow cross indicates the centroid position of B335.

4 Discussion

4.1 Linewidths and kinetic temperature

The main radius probed with the C¹⁷O emission was computed from the 3σ contour and found to be 860 au. A region enclosing radii from 100 (about half of the FWHM beam) to 860 au was chosen for the gas kinematics analysis. The kinetic temperature of the gas has been estimated from the formula for dust temperature in an optically thin regime assuming only central heating by the B335 protostar derived by Shirley et al. (2011). The underlying assumption is that dust and gas are expected to be in thermal equilibrium, being coupled via collisions at the densities probed here (> 10⁵ cm⁻³). We use Eq. (2) (Evans et al. 2015) which we adapted to the new distance of 164.5 pc: (2)

The gas kinetic temperature was computed for the two radii probed, with values in the range of T_k (100 au) = 46 K and T_k(860 au) = 20 K. The observed linewidths obtained in the previous section have been compared with the expected thermal linewidth, given by: (3)

where T_k is the gas kinetic temperature, m is the molecular mass (29.01 amu for C¹⁷O) and k_b is the Boltzmann constant. The expected thermal linewidth for C¹⁷O has been computed for the temperatures at the two different radii: Δv_th(100 au) = 0.27 km s⁻¹ and Δv_th (860 au) = 0.17 km s⁻¹. The observed linewidths are larger than the thermal ones for both velocity components. This indicates that the observed linewidth is the result of the thermal component plus a non-thermal contribution (e.g., turbulence and large-scale motions such as infall and outflow, = + ). The non-thermalcontribution of the line has been computed for both velocity components and the results are shown in Table 2. The non-thermal component at the inner and outer radius are indistinguishable because of the limited spectral resolution (0.2 km s⁻¹).

The sound speed, c_s, is in the 0.2–0.3 km s⁻¹ range for temperatures between 20 and 46 K and γ ~ 7/5. This means that the non-thermal contribution to the linewidth is supersonic.

Simulations have shown that in star forming cores systematic large-scale motions (such as infall) can contribute significantly (~50%) to the non-thermal component of the linewidth (Guerrero-Gamboa & Vázquez-Semadeni 2020). We can do a rough estimation of the contribution from infall, by measuring the infall velocity difference from two different radii, Δσ = v_ff(r1) − v_ff(r2), where , G is the gravitational constant and M_B335(r) is the mass enclosed at the considered radius r. To compute the mass, we consider that the total mass at different radii is the sum of the mass in the envelope plus the mass of the central object. The mass of the gas in the envelope contained up to a certain radius is computed by integrating the 110 GHz dust continuum emission and using standard assumptions (see Eq. (4) in Jørgensen et al. 2007). The maximum velocity difference along the LOS would occur at the core’s center. For the two adopted radius, 860 and 100 au, the envelope mass, M_env, is ~0.16 M_⊙ and ~0.012 M_⊙, respectively. We also adopt a mass for the central object between 0.05 and 0.26 M_⊙, which are the predicted values from infall models (see Yen et al. 2015; Evans et al. 2015, and in prep.). Therefore, the total mass enclosed within a 100 and 860 au radii are in the 0.06–0.16 M_⊙ and 0.27–0.41 M_⊙ ranges, respectively. The estimated range of Δν_ff is 0.46–1.01 km s⁻¹. Given the errors coming from the computation of all the previous parameters, the non-thermal contribution of the linewidth is of the same order as the broadening due to free-fall motions along the line of sight: it is, hence, possible that the observed velocity pattern is due to infall.

Fig. 5 C17O (1–0) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at 5, 10, 20 and 30σ, where σ is 10.0 mJy beam−1. Top: peak velocity. Middle: linewidt. Bottom: opacity. Each column shows one of the two velocity components. Red crosses indicate the centroid position of B335.

Fig. 6 Opacity histograms from the fitting of the two individual velocity components.

4.2 Possible origins of the observed gas motions

Our ALMA observations of the C¹⁷O (1–0) emission suggest an optically thin emission at all scales probed by the observations (100–860 au). Overall, the spatial extent of the C¹⁷O emission is similar to the one of the dust continuum emission, but C¹⁷O is less peaked and decreases more smoothly with decreasing density outwards; this suggests that the gas traced with C¹⁷O is not entirely related to the outflow cavity. The C¹⁷O emission maximum isnot coincident with the dust continuum peak position, which might suggest slight abundance variations of the C¹⁷O at high densities close to the protostar. Nevertheless, the prominence of the double-peaked velocity pattern is not correlated with the intensity of the dust continuum emission, demonstrating that those profiles are not due to red-shifted absorption against a strong continuum or a source of C¹⁷O. Thus, these two velocity components trace distinguished gas motions. A simple isotropic inside-out envelope collapse can not easily reproduce the gas motions we observe in the B335 envelope. In this section, we discuss various hypothesis for the physical origin of thegas motions observed.

Despite being isolated, B335 is embedded in an extended molecular gas cloud of density ~ 10³ cm⁻³ (Frerking et al. 1987). However, C¹⁷O is a rare isotopologue which is mostly confined to a high-density central region and its low abundance at large-scales would prevent observing such a tenuous layer, suggesting there is no missing flux coming from large-scale C¹⁷O emission. To confirm that this is the case, we estimated the missing flux from the C¹⁸O (1–0) ALMA observations by comparing it with the 45 m Nobeyama data of the same transition presented in Saito et al. (1999). Our C¹⁸O map was smoothed to match the beam size of the Nobeyama telescope, which at the frequency of this transition (109.782 GHz) is 16″. We obtained the spectra on a region of one beam size around the center of B335 and transformed the flux density to brightness temperature using the Rayleigh-Jeans law. We obtained a peak temperature of T_MB = 0.79 ± 0.08 K and an integrated temperature of ∫ T_MBdv = 0.53 K km s⁻¹. This means that our ALMA observations are recovering around 14% of the total flux detected with single-dish data. However, because C¹⁷O is expected to be much more compact than C¹⁸O we expect the missing flux to be much less for the former. Frerking et al. (1987) presented single-dish data of the C¹⁷O (1–0) transition and concluded that all their emission is coming from the center of the source in a region smaller than the beam of their telescope (1′.6 for the C¹⁷O (1–0) transition). This extension is much smaller than the one observed for C¹⁸O (1–0) detected in both works, which is about 4′. This is consistent with the fact that C¹⁷O is much less abundant than C¹⁸O, especially at large scales, and that it is mainly tracing the core and not the envelope. Therefore, we expect the missing flux of C¹⁷O in our ALMA data to be much less and, in addition, to be recovering at least twice the recovered flux of C¹⁸O, namely, around 30%. We also note that while our observation might be missing flux, this should not be enough to produce the huge dip in our data and it cannot explain the structured velocity pattern we observe in the spectral maps, since the missing flux will be at the systemic velocity and would not be able to completely absorb only one of the two components at different offsets.

A possible cause for the blueshifted and redshifted gas motions observed in protostellar envelopes could be organized core rotation. Our observations do not support this hypothesis as they do not show a clear velocity gradient in the equatorial plane, where rotation motions would mostly contribute to the observed velocity field. Instead, both redshifted and blueshifted velocities are observed in both the northern and southern regions (see Fig. 5). While rotation motions have only been detected at larger envelope radii in B335 (>2500 au, Saito et al. 1999; Yen et al. 2011), we stress that the conclusion regarding the absence of small-scale rotation (e.g., Yen et al. 2010)should be further investigated using the new insights on gas motions in the envelope that our observations have uncovered.

B335 has a well-studied outflow, with its axis close to the plane of the sky and with a well-defined X-shaped biconical shape in ¹²CO (Bjerkeli et al. 2019). Although some contamination by the gas from the outflow cannot be completely ruled out, we present here arguments supporting the hypothesis that our C¹⁷O maps can be used to study the envelope gaskinematics. C¹⁷O is a rare isotopologue which is known to trace dense envelope gas and is not expected to be detected in more tenuous outflow cavities. The morphology of C¹⁷O emission is very different from the one observed in typical outflow cavities tracers, such as C₂H (Imai et al. 2016) or ¹²CO (see bottom image in Fig. 1 and Bjerkeli et al. 2019). Moreover, no spectral signature of outflow is observed, such as large wings observed in ¹²CO (Bjerkeli et al. 2019), and the maximum velocity shift from the rest velocity remains quite small (±1 km s⁻¹). Therefore, the kinematic pattern observed in our C¹⁷O maps cannot be produced by outflow alone, and it does provide a strong evidence of distinguished velocity contributions from the gas in the inner region of the B335 protostellar envelope.

The C¹⁷O velocity maps in Fig. 5 show that the largest gas velocities are found ~ 1′′ from the center along the two northern outflow cavity walls, tracing gas at reverse velocities with respect to the outflow velocities. Considering the 10° inclination of the system, the spatial distribution of C¹⁷O emission following closely that of the dust and other typical dense gas tracers; and the fact that the linewidths of the two velocity components are in general agreement with the expected linewidths from infall motions, the most likely hypothesis is that these high-velocity (±1 km s⁻¹) features trace, accreting gas flowing along the outflow cavity walls and onto the central protostar. The peak velocities tentatively increase towards the central protostellar objects, for the features along the eastern outflow cavity walls, but no clear velocity gradient could be resolved in the current observations: additional observations with better spatial resolution may allow to test further this hypothesis. Finally, we note that the strongly redshifted emission at the north-east was already detected in ALMA C¹⁸O observations reported by Yen et al. 2015 (see Fig. 2 in their work).

Dust continuum emission observed with ALMA at various millimeter and sub-millimeter wavelengths all show a striking excess of dust emission associated to the outflow cavity walls. While this could be a temperature effect due to increased heating from the central protostar of these walls, it could also be a true density increase in compact features easily picked up by interferometric observations. Magnetized models of protostellar formation (for a review see Zhao et al. 2020) suggest cavity walls could be preferential sites to develop accretion streamers, as observed in the non-ideal magneto-hydrodynamic (MHD) models of protostellar accretion and outflow launching (Machida 2014; Figs 8 and 9 in Tomida et al. 2012). Indeed, these are locations where the poloidal magnetic field is mostly parallel to the inflow direction and, therefore, would exert less magnetic braking for material infalling along the walls. This hypothesis is also in agreement with the dust polarization observations of magnetic field lines in B335 (the redshifted gas feature we observe along the north-eastern cavity wall is associated to highly organized B-field lines aligned with the tentative gas flow) and the scenario of magnetically-regulated infall proposed in Maury et al. (2018).

We note that the observed non-thermal components of the linewidths are found to be supersonic. If the observed gas motions we detect indeed trace localized accretion motions, these could be supersonic. While the development of supersonic filamentary accretion features were reported in numerical models of protostellar formation (Padoan et al. 2005; Banerjee et al. 2006; Kuffmeier et al. 2019) and the observations suggested supersonic infall is occurring in a few protostellar envelopes at larger scales (> 1000 au, Tobin et al. 2010; Mottram et al. 2013), it is the first time such anisotropic supersonic infall motions have been tentatively reported in the B335 inner envelope.

4.3 Impact on the characterization of protostellar mass accretion rates

In the following, we briefly discuss the implications of our work and whether the localized accretion features detected in B335 are common, allthe while remaining mostly unresolved in many observations of accreting protostars.

Self-similar solutions for analytical models of the collapse of an isothermal sphere, including only thermal pressure and gravity, predict typical mass accretion rates on the order of ~ 10⁻⁴ M_⊙ yr⁻¹ (Larson 1969; Penston 1969; Shu 1977). Turbulent models and MHD numerical models have produced slightly lower mass accretion rates of ~ 10⁻⁶–10⁻⁵ M_⊙ yr⁻¹. Episodic accretion with highly variable rates (from a Ṁ ~ 10⁻⁵ M_⊙ yr⁻¹ down to Ṁ < 10⁻⁶ M_⊙ yr⁻¹) is often observed in both hydro and MHD numerical models of protostellar formation, in the accretion of envelope material onto the disk and the protostar itself (Lee et al. 2021), and of disk material to the central growing protostar (Dunham & Vorobyov 2012; Vorobyov & Basu 2015). Robust observational estimates of protostellar accretion rates are crucial to distinguish between models, but also to shed light on several open questions on star formation, since they are key quantities for our interpretation of the protostellar luminosities and of the typical duration of the main protostellar accretion phase (Evans et al. 2009; Maury et al. 2011). Indeed, observations may have revealed a discrepancy between the observed protostellar bolometric luminosities, and the protostellar accretions rates: this is the so-called “luminosity problem” (Kenyon et al. 1990). Protostellar accretion rates derived from molecular line emission – and more particularly from the modeling of inverse P Cygni profiles with analytical infall models – should produce luminosities 10–100 times larger than the typically observed bolometric luminosities (for a review, see Dunham et al. 2014). Observations of the molecular line profiles in B335 (Evans et al. 2005, 2015) have been used to fit models of protostellar infall suggesting Ṁ ~ 6.2 × 10⁻⁶ M_⊙ yr⁻¹ (assuming an effective sound speed of 0.3 km s⁻¹, Evans et al., in prep.), although arguably these estimates are associated to large error bars. The observed bolometric luminosity of B335 lies an order of magnitude below the accretion luminositym L_acc, that such accretion rates should produce (Evans et al. 2015); despite being a prototype for protostellar infall models, B335 also suffers from the luminosity problem.

If the “actual” protostellar mass accretion rates stem from localized collapsing gas at small scales, potentially affected by unresolved multiple velocity components, the true linewidths associated with infalling gas feeding the growth of the protostar would be quadratically smaller than the ones measured at larger scales where these components would be entangled (or if the individual velocity components are interpreted as being part of a single velocity component with a central dip due to optically thickness). Smaller intrinsic linewidths of the infalling gas, at small radii, may result in smaller effective sound speed and, hence, lower mass accretion rate derived from the analytical infall models, since . It is, therefore, possible that the observed bolometric luminosity of B335 is, ultimately, compatible with its accretion luminosity L_acc. In the revised B335 scenario we propose here, the mass accretion rate on the protostar could be dictated by localized supersonic infall rather than by the large scale infall rate of the envelope: this may open a window to partially solving the “luminosity problem”, although episodic vigorous accretion would probably continue to be necessary to explain the relatively short statistical lifetimes for the Class 0 phase. Future observations should be used to carry out a more detailed characterization of whether a significant fraction of the final stellar mass could be fed to the central object through highly localized anisotropic infall. Recently, Pineda et al. (2020) reported the detection of an “accretion streamer” connecting the dense core to disk scales and found a streamer infall rate ~ 10⁻⁶ M_⊙ yr⁻¹, of the same order of magnitude as the global mass accretion rate inferred from molecular line observations in B335. Hence, it is possible that many previous studies have failed to grasp the full complexity of the gas motions making up the accretion onto the central protostars. It is the large angular and spectral resolution, along with the great sensitivity of the presented ALMA observations which allowed to detect the two distinct velocity components on the line emission profiles in B335. More observations at the same small-scales and spectral resolution of optically thin emission in different protostars are needed in order to determine if these localized gas motions are common. Moreover, more refined protostellar infall models will have to be carried out in the future to take this new small scales into account and include more complex geometries with, for instance, asymmetric structures and preferential accretion along outflow cavities.

5 Summary

In this work, we describe ALMA observations of the C¹⁷O emission tracing gas kinematics in the B335 envelope and we show that the line emission exhibits widespread double-peaked profiles. Based on the analysis, we present the following conclusions:

• Derivations of the line opacity have shown that the emission of the line is optically thin and, therefore, the observed double-peaked profiles cannot be produced by self-absorption. Thus, inverse P Cygni profiles coming from symmetrical inside-out collapse cannot fully explain the observed complex velocity field.

• After discarding filtering of large-scale emission or other types of motions, such as rotation or outflow, we determined that only distinct gas motions contributing to the same line-of-sight (LOS) could explain the observed line profile pattern.

• The linewidth analysis determined that the two velocity components are compatible with infall motions and could be due to localized infall in preferential directions. The main hypothesis presented here is that the collapse of the envelopeonto the protostar is occurring along the equatorial plane, but also along the outflow cavity walls, where the magnetic field topology is more favorable.

More observations at similar scales and spectral resolution are needed to determine whether these double-peaked profiles are common in protostellar objects at similar evolutionary states. Moreover, further modeling of the B335 envelope with more complex collapse models, such as anisotropic collapse, are needed to determine the exact physical origin of the observed velocity field.

Acknowledgments

The authors acknowledge the very useful discussions with N. Evans and Y. Yang. This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (MagneticYSOs project, grant agreement N° 679937). J.M.G. is supported by the grant AYA2017-84390-C2-R (AEI/FEDER, UE). This publication is based on data of ALMA data from the project 2016.1.01552.S (PI: A. Maury).

Appendix A Technical details of the ALMA observations

Appendix B C¹⁷O full channel maps

Fig. B.1 Contour channel maps of the C17O (1-0) emission showing the two main hyperfine components, overlapped with the gray scale image of the 1.3 mm dust emission. Contours are − 4, − 2, and from 2 to 18 by steps of 2σ, where σ is 10 mJy Beam−1. The synthesized beam is shown in the bottom left panel as a filled ellipse. The vLSR channel velocities are shown in the top-left part of the panels.

Appendix C C¹⁸O spectral map

Fig. C.1 Spectral map of the C18O (1-0) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5′′ × 5.5′′ and each pixel correspond to 0.5′′ (~ 82 au). The green refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1).

References

All Tables

All Figures

Fig. 1 Red contours showing moment 0 of C17O emission, integrated over the velocity range 4.8–6.2 km s−1 and 7.6–9.4 km s−1. Contours show emission at −3, 3, 5, 10 and 30σ, where σ is 10.0 mJy beam−1. Top: intensity shows dust continuum emission map at 110 GHz for emission over 3σ, where σ is 8.56 × 10−2 mJy beam−1. Bottom: intensity shows 12CO (2–1) moment 0 emission integrated over the velocity range 1.4–16.2 km s−1, for emission over 3σ, where σ is 2016.9 mJy beam−1. The two arrows show the direction of the E-W outflow.

In the text

	Fig. 2 Spectral map of the C17O (1–0) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel correspond to 0.5″ (~82 au). For clarity, the spectral range only shows the main hyperfine component. The green spectrum refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1).
In the text

Fig. 3 Contour channel maps of the emission of the F = 5/2–(−5/2) hyperfine component of the C17O J = 1–0 transition,overlapped with the gray scale image of the 1.3 mm dust emission. Contours are − 4, − 2, and from 2 to 18 by steps of 2σ, where σ is 10 mJy beam−1. The synthesized beam is shown in the bottom left panel as a filled ellipse. The vLSR channel velocities are shown in the top left part of the panels.

In the text

	Fig. 4 C17O to C18O integrated intensity ratio. Red contours show integrated intensity for C17O at −3, 3, 5, 10 and 30σ, where σ is 10.0 mJy beam−1. Black contours show integrated intensity for C18O at −3, 3, 5, 10 and 30σ, where σ is 10.7 mJy beam−1. The yellow cross indicates the centroid position of B335.
In the text

	Fig. 5 C17O (1–0) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at 5, 10, 20 and 30σ, where σ is 10.0 mJy beam−1. Top: peak velocity. Middle: linewidt. Bottom: opacity. Each column shows one of the two velocity components. Red crosses indicate the centroid position of B335.
In the text

	Fig. 6 Opacity histograms from the fitting of the two individual velocity components.
In the text

	Fig. B.1 Contour channel maps of the C17O (1-0) emission showing the two main hyperfine components, overlapped with the gray scale image of the 1.3 mm dust emission. Contours are − 4, − 2, and from 2 to 18 by steps of 2σ, where σ is 10 mJy Beam−1. The synthesized beam is shown in the bottom left panel as a filled ellipse. The vLSR channel velocities are shown in the top-left part of the panels.
In the text

	Fig. C.1 Spectral map of the C18O (1-0) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5′′ × 5.5′′ and each pixel correspond to 0.5′′ (~ 82 au). The green refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1).
In the text