Structured velocity field in the inner envelope of B335: ALMA observations of rare CO isotopologues

Studying Class 0 objects is very important, as it allows to characterize dynamical processes at the onset of the star formation process, and to determine the physical mechanisms responsible for the outcome of the collapse. Observations of dense gas tracers allow the characterization of key kinematics of the gas directly involved in the star-formation process, such as infall, outflow or rotation. This work aims at investigating the molecular line velocity profiles of the Class 0 protostellar object B335 and attempts to put constraints on the infall motions happening in the circumstellar gas of the object.} Observations of C$^{17}$O (1-0), C$^{18}$O (1-0) and $^{12}CO$ (2-1) transitions are presented and the spectral profiles are analyzed at envelope radii between 100 and 860 au. C$^{17}$O emission presents a double peaked line profile distributed in a complex velocity field. Both peaks present an offset of 0.2 to 1 km s$^{-1}$ from the systemic velocity of the source in the probed area. The optical depth of the C$^{17}$O emission has been estimated and found to be less than 1, suggesting that the two velocity peaks trace two distinct velocity components of the gas in the inner envelope. After discarding possible motions that could produce the complex velocity pattern, such as rotation and outflow, it is concluded that infall is producing the velocity field. Because inside-out symmetric collapse cannot explain those observed profiles, it is suggested that those are produced by non-isotropic accretion from the envelope into the central source along the outflow cavity walls.


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 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 to understand which are the kinematics and dynamics of the gas at the onset of collapse, and to determine how those affect the outcome of the star formation process.
The gas making most of the protostellar envelopes is typically probed using molecular gas line profiles, which trace the gas kinematics in the dense envelope and are used to measure gas motions such as rotation or infall. Observations of the molecular line emission from embedded protostars have been suggested 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 5 cm −3 , André 1995), the large optical depth of some molecular line emission produce 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 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 double-peaked profile on optically thick emission, to extract key information such as central protostellar mass, infall velocities and mass accretion rates (Zhou et al. 1993, Di Francesco et al. 2001, Evans et al. 2005, Evans et al. 2015. However, blue-asymmetries in line profiles are not unique to infall motions. Complex gas kinematics (Maureira et al. 2017), such as asymmetric collapse (Tokuda et al. 2014), accretion streamers , Segura-Cox et al. 2020 or outflow-entrained gas, can produce separated velocity components on the same line-of-sight which are observed as doublepeaked line profiles 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 to test symmetrical collapse infall models, since blueasymmetries were first detected in 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). Infall models in an optically thick line emission have attempted to compute infall mass rates (Yen et al. 2010, Evans et al. 2015, Yen et al. 2015, obtaining values affected by large uncertainties that range from 10 −7 M yr −1 to ∼3×10 −6 M yr −1 at radii of 100-2000 au, and infall velocities from 1.5 km s −1 to ≈0.8 km s −1 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 −6 M yr −1 (Evans et al., in prep.). For the estimated age of 4×10 4 yr, this implies a total mass at the center (star + disk) of 0.26 M . The physical cause of these double peaked line profiles has been questioned, however. For example Kurono et al. (2013) pointed out that despite expecting the H 13 CO + emission should 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, prominently detected in 12 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, Yen et al. 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, Yen et al. 2018. Recent observations of the hourglass 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 17 O (1-0) and C 18 O (1-0), which trace the dense circumstellar gas of the inner envelope of B335, are presented along with the 12 CO (2-1) line emission, tracing the outflow cavity. The molecular line profiles are analyzed and interpreted, giving new constraints on the gas kinematics close to the protostar.

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 project 2016.1.01552.S (PI. A. Maury). In all 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 17 O (1-0) and C 18 O (1-0) were targeted using two configurations, C40-2 and C40-5, and 12 CO (2-1) was targeted using C40-1 and C40-4. Technical details of the observations are shown in Table A.1. Preliminary analysis of the data was done with the product images delivered by ALMA to check if emission was detected and to check the shape of the line profiles. The C 17 O emission was only detected in the most compact configuration (C40-2), therefore only this configuration has been used to produce C 17 O and C 18 O maps, while 12 CO was detected in both configurations so a combination of the two data sets has been 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. 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 17 O and C 18 O and 1 for 12 CO. After imaging, 12 CO maps were smoothed to reach the same angular resolution as C 17 O and C 18 O. The resulting map characteristics are shown in Table 1.   Figure 1 shows the moment 0 of C 17 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 the moment 0 of 12 CO (2-1) emission. The 12 CO emission probes the outflowing gas, therefore it is confirmed that the C 17 O emission is tracing the envelope and it is not affected by the outflow.

Results and analysis
In order to understand the dynamics of the gas that is being probed, the C 17 O spectra were taken at every 0.5 pixel of the emission cube, producing a spectral map shown on Fig. 2. The line profile patterns show two distinct velocity components with a dip centered around the systemic velocity (8.3 km s −1 ). Their respective intensities vary depending on the direction of the offset from the continuum peak, being the blue component 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 detected hyperfine components. 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 Contours show emission at −3, 3, 5, 10 and 30 σ, where σ is 10.0 mJy/beam. Top: Intensity shows dust continuum emission map at 110 GHz for emission over 3σ, where σ is 8.56×10 −2 mJy/beam. Bottom: Intensity shows 12 CO (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. The two arrows show the direction of the E-W outflow.
components due to other dynamical processes. The possibility of the dip being caused by interferometric filtering is discarded since the recovered emission size is of the order of the Largest Recoverable Scale (see Table 1). Moreover, because C 17 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 17 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 redshifted 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 17 O spectral profiles and if 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.

Line opacity estimation
The maximum opacity at the center of the source has been estimated from the H 2 column density and assuming a standard C 17 O abundance using Eq. 1 (Jansen 1995). (Gordon et al. 2017, Müller et al. 2005, c is the speed of light, ν is the frequency of this transition, N H 2 0 = 3.1 × 10 22 cm −2 is the peak column density of H 2 in B335 at a radius of 3600 au (Launhardt et al. 2013), [C 17 O] = 5 × 10 −8 is the C 17 O abundance relative to H 2 abundance (Thomas & Fuller 2008) and ∆V ≈ 1 kms −1 is the average observed linewidth for the two components together. Using these values the obtained opacity is τ 0 = 0.77, which corresponds to an opacity typical from an optically thin line.

Intensity ratio
Because of their similar mass and molecular structure, C 17 O and C 18 O should probe gas under similar physical conditions. The [C 17 O]/[C 18 O] isotope ratio does not appear to be affected by fractionation and, if emission from both C 17 O and C 18 O is associated to dense gas shielded from external ultraviolet radiation, selective photo-dissociation probably does not affect the relative abundances (van Dishoeck & Black 1988). Thus, the only difference in the emission from these two isotopes should result from opacity effects, because C 18 O is a factor of 3.6-3.9 more abundant than C 17 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 17 O emission.
We produced beam-matching maps for C 17 O and C 18 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 −1 and 7.4-9.4 km s −1 for C 17 O hence taking into account the two observed hyperfine components, and 7.4-9.4 km s −1 for C 18 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 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 17 O in those regions when compared to the C 18 O emission (see C 18 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  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. Note that this is also in agreement with the conclusions reached by the analysis of single-dish observations of the C 18 O (2 − 1) and C 17 O (2 − 1) (τ C 18 O (2−1) ∼ 0.8, Evans et al. 2005).

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 simultaneously four independent parameters: 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 χ 2 . Because the C 17 O (1-0) emission is expected to be optically thin, it was at-tempted to model the double-peaked profiles with two velocity components, one blue-and another red-shifted. Initial expected values were introduced and the program is allowed to proceed to fit emission with a minimum signal-to-noise ratio of 3σ. The peak velocity, the linewidth and the opacity maps obtained from the fitting are shown in Fig. 5.
Velocity maps (top images of Fig. 5) show that the two different components, blue and red-shifted, occupy generally two separated regions, at East and West offsets from the center of the source respectively, with some regions overlapping, where the double peak can be observed in the spectra. The mean average, velocity for the two components are 8.1 and 8.6 km s −1 respectively. The velocity dispersion maps (middle images of Fig. 5) show a mean velocity dispersion for the two components of 0.7 and 0.5 km s −1 , which get broader, up to 0.9∼1, closer to the center 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 also have a very large error associated, of 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.

Linewidths and kinetic temperature
The main radius probed with the C 17 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 has been chosen to analyze the gas kinematics. 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 5 cm −3 ). We use Eq. 2 (Evans et al. 2015) which we adapted to the new distance of 164.5 pc. The gas kinetic temperature was computed for the two radii probed, with values in 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: where T k is the gas kinetic temperature, m is the molecular mass (29.01 amu for C 17 O) and k b is the Boltzmann constant. The expected thermal linewidth for C 17 O has been computed for the temperatures at the two different radii: ∆v th (100 au) = 0.27 km s −1 and ∆v th (860 au) = 0.17 km s −1 . 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 like infall and outflow, v 2 obs = v 2 th + v 2 non−th ). The non-thermal contribution 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 −1 ).
The sound speed, c s , is in the 0.2-0.3 km s −1 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, ∆σ = , 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 line-of-sight 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. 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 ∆ν f f is 0.46 -1.01 km s −1 . 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.

Possible origins of the observed gas motions
Our ALMA observations of the C 17 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 17 O emission is similar to the one of the dust continuum emission, but C 17 O is less peaked and decreases more smoothly with decreasing density outwards: this suggests that the gas traced with C 17 O is not mostly related to the outflow cavity. The C 17 O emission maximum is not coincident with the dust continuum peak position, which might suggest slight abundance variations of the C 17 O at high densities close to the protostar. Nevertheless, the prominence of the double-peaked velocity pattern does not correlate with the intensity of the dust continuum emission, proving that those profiles are not due to red-shifted absorption against a strong continuum and/or C 17 O source. These two velocity components thus 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 the gas motions observed. Despite being isolated, B335 is embedded in an extended molecular gas cloud of density ∼ 10 3 cm −3 (Frerking et al. 1987). However, C 17 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 17 O emission. To confirm that this is the case, we estimated the missing flux from the C 18 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 18 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 MB dv = 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 17 O is expected to be much more compact than C 18 O we expect the missing flux to be much less for the former. Frerking et al. (1987) presented single-dish data of the C 17 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 17 O (1-0) transition). This extension is much smaller than the one observed for C 18 O (1-0) detected in both works, which is about 4'. This is consistent with the fact that C 17 O is much less abundant than C 18 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 17 O in our ALMA data to be much less, and to be recovering at least twice the recovered flux of C 18 O, i.e. around 30 %. We also note that while our observation might be missing flux, this should not be enough as to produce the huge dip in our data, and it can not 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 observing blueshifted and redshifted gas motions 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. 4). 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 12 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 17 O maps can be used to study the envelope gas kinematics. C 17 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 17 O emission is very different from the one observed in typical outflow cavities tracers, such as C 2 H (Murillo et al. 2018) or 12 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 12 CO (Bjerkeli et al. 2019), and the maximum velocity shift from the rest velocity remains quite small (± 1 km s −1 ). Therefore, the kinematic pattern observed in our C 17 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 17 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 17 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 −1 ) features trace accreting gas flowing along the outflow cavity walls, 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 18 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 or Figure 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 Bfield 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 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 are tentatively reported in the B335 inner envelope.

Impact on the characterization of protostellar mass accretion rates
In the following, we briefly discuss the implications of our work, if the localized accretion features detected in B335 are common 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 of the order ∼ 10 −4 M yr −1 (Larson 1969;Penston 1969;Shu 1977). Turbulent models and MHD numerical models have produced slightly lower mass accretion rates ∼ 10 −6 -10 −5 M yr −1 . Episodic accretion with highly variable rates (from aṀ ∼ 10 −5 M yr −1 down toṀ < 10 −6 M yr −1 ) 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(Evans et al. , 2015 have been used to fit models of protostellar infall suggestingṀ ∼ 6.2 × 10 −6 M yr −1 (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 luminosity L acc 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 'true' protostellar mass accretion rates stems from localized collapsing gas at small scales, potentially affected by unresolved multiple velocity components, the true linewidths associated to 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Ṁ ∝ c 3 s /G. 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 solve the 'luminosity problem', although episodic vigorous accretion would probably still remain 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 −6 M yr-1, 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, e.g., asymmetric structures and preferential accretion along outflow cavities.

Summary
ALMA observations of the C 17 O emission tracing gas kinematics in the B335 envelope have been presented in this work. It is shown that the line emission exhibits widespread double-peaked profiles. From the analysis, the following conclusions have been obtained: -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 selfabsorption. Therefore, 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, it is determined that only distinct gas motions contributing to the same lineof-sight could explain the observed line profiles pattern. -Linewidth analysis have 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 is that the collapse of the envelope onto 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 if 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 which is the exact physical origin of the observed velocity field. (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 spectra identifies the spectrum at the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s −1 ).