© The Authors 2023

1 Introduction

Observational studies have shown that magnetic fields are present at all scales where star formation processes are at work, permeating molecular clouds, prestellar cores, and protostellar envelopes (e.g., Girart et al. 2006; Alves et al. 2014; Zhang et al. 2014; Soler 2019). While magnetized models of star formation were developed early on, it is only recently that complex physics has been included in numerical magnetohydrodynamic (MHD) models, and that the predictions of dedicated models are directly compared to observations of star-forming regions. One of the major achievements of magnetized models is the prediction of the small sizes of protostellar disks, regulated by magnetic braking (Dapp et al. 2012; Hennebelle & Inutsuka 2019; Zhao et al. 2020). Recent observations have indeed shown that disks are compact, and are almost an order of magnitude smaller than those produced in purely hydrodynamical models (Maury et al. 2019; Lebreuilly et al. 2021).

In magnetized collapse models, non-ideal MHD effects play a major role in the regulation of magnetic flux during the early stages of protostellar formation (Machida et al. 2011; Li et al. 2011; Marchand et al. 2020). Hence, these processes indirectly limit the angular momentum transported by the magnetic field from large envelope scales to the inner envelope, and set the resulting disk size. However, these resistive processes depend heavily on the local physical conditions, such as dust grain properties (electric charge and size), gas density and temperature, density of ions and electrons, and the cosmic-ray (CR) ionization rate, ζ_CR (namely the number of molecular hydrogen ionizations per unit time; see e.g., Zhao et al. 2020, for a review). Observationally, only very few measurements of these quantities have been obtained toward solar-type protostars: the effective coupling of magnetic fields to the infalling-rotating material, at scales where the circumstellar material feeds the growth in mass of the protostellar embryo and its disk, is therefore still largely unknown.

B335, located at a distance of 164.5 pc (Watson 2020), is an isolated Bok Globule which contains an embedded Class 0 protostar (Keene et al. 1983). The core 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; Stutz et al. 2008; Yen et al. 2010). Its isolation and relative proximity makes B335 an ideal object with which to test models of star formation. Asymmetric double-peaked line profiles observed in the molecular emission of the gas toward the envelope have been interpreted as optically thick lines tracing infalling motions, and have been extensively used to derive mass-infall rates from 10⁻⁷ to ~3 × 10⁻⁶ M_⊙ yr⁻¹ at radii of 100-2000 au (Zhou et al. 1993; Yen et al. 2010; Evans et al. 2015). However, new observations of optically thin emission from less abundant molecules reveal the presence of two velocity components tracing asymmetric motions at these envelope scales, which could contribute significantly to the double-peaked line profiles (Cabedo et al. 2021). These results suggest that simple spherically symmetric infall models might not be adequate to describe the collapse of the B335 envelope, and tentatively unveil the existence of preferential accretion streamlines along outflow cavity walls.

B335 is also an excellent prototype for the study of magnetized star-formation models, because it has been proposed as an example of a protostar where the disk size is set by a magnetically regulated collapse (Maury et al. 2018, hereafter Paper I). This hypothesis rests mainly on two observations. First, the rotation of the gas observed at large envelope radii is not found at small envelope radii (<1000 au) and no kinematic signature of a disk was reported down to ~ 10 au (Kurono et al. 2013; Yen et al. 2015). Second, observations of polarized dust emission have revealed an “hourglass” magnetic field morphology (Paper I) in the inner envelope: comparisons to MHD models of protostellar formation (see e.g., Masson et al. 2016; Hennebelle et al. 2020, for a whole description of the models) suggest that the B335 envelope is threaded by a rather strong magnetic field which is highly coupled to the infalling-rotating gas, preventing disk growth out to large radii. The purpose of the analysis presented here is to characterize the level of ionization and its origin in order to test the model scenario proposed in Paper I, and put observational constraints on the efficiency of the coupling of the magnetic field to the star-forming gas in the inner envelope of B335.

In Sect. 2, we present the ALMA observations used to constrain physical and chemical properties of the gas at envelope radii ≲1000 au. In Sect. 3, we derive the deuteration fraction, R_d, from DCO⁺ (J = 3–2) and H¹³CO⁺ (J = 3–2). In Sect. 4, we compute the ionization fraction, χ_e, and the ζ_CR. Finally, in Sect. 5, we discuss our results concerning deuteration processes, the ionization, the possible influence of a local source of radiation and its effect on the coupling between the gas and the magnetic field.

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. The centroid position of B335 is assumed to be α = 19 : 37 : 00.9 and δ = +07 : 34 : 09.6 (in J2000 coordinates) based on the dust continuum peak obtained by Maury et al. (2018).

We used DCO⁺ (J = 3–2) and H¹³CO⁺ (J = 3–2) to obtain R_D and χ_e. We used H¹³CO+ (J = 1–0) and C¹⁷O (J = 1–0) to derive the hydrogenation fraction, R_H. Additionally, we observed the dust continuum at 110 GHz to derive the estimated line opacities, CO depletion factor, f_D, and ζ_CR. We also targeted ¹²CO (J = 2–1) and N₂D⁺ (J = 3–2) for comparison purposes. All lines were targeted using a combination of ALMA configurations to recover the largest length-scale range possible. As our data for H¹³CO⁺ (J = 3–2) only include observations with the Atacama Compact Array (ACA), we used additional observations of this molecular line at smaller scales from the ALMA project 2015.1.01188.S. Technical details of the ALMA observations are shown in Appendix A.

Calibration of the raw data was done using the standard script for Cycle 4 ALMA data and the Common Astronomy Software Applications (CASA), version 5.6.1–8. The continuum emission was self-calibrated with CASA. Line emission from ¹²CO (J = 2–1) and N₂D⁺ (J = 3–2) was additionally calibrated using the self-calibrated continuum model at the appropriate frequency (231 GHz, not shown in this work).

Final images of the data were generated from the calibrated visibilities using the tCLEAN algorithm within CASA, using Briggs weighting with a robust parameter of 1 for all the tracers. As we want to compare our data to the C¹⁷O (J = 1–0) emission presented in Cabedo et al. (2021), we adjust our imaging parameters to obtain matching beam maps with similar angular and spectral resolution. For DCO⁺ (J = 3–2) and H¹³CO⁺ (J = 3–2), we restricted the baselines to a common “u, v” range between 9 and 140 kλ, and finally smoothed them to the same angular resolution. The preliminary analysis shows that both C¹⁷O (J = 1–0) and H¹³CO⁺ (J = 1–0) emission is barely detected in the most extended configurations. We proceed by applying a common “uv tapering” of 1.5″ to both lines. This procedure allows to decrease the weighting of the largest baselines in our analysis, giving more weight to the smaller baselines and allowing us to obtain a better signal-to-noise ratio (S/N). Furthermore, we smoothed the data to obtain exactly the same angular resolution. Additionally, we smooth C¹⁷O (J = 1–0) and the dust continuum at 110 GHz to obtain DCO⁺ matching beam maps to compute f_D. The imaging parameters used to obtain all the spectral cubes and their final characteristics are shown in Table 1. Even though the characteristics and properties of the C¹⁷O (J = 1–0) maps are slightly different from the ones presented in Cabedo et al. (2021) due to the differences in the imaging process, the line profiles show the same characteristics (i.e., line profiles are double-peaked and present the same velocity pattern at different offsets from the center of the source) confirming that the imaging process has no large effect on the shape of the line profile, and that the two velocity components can still be observed.

Data products

From the data cubes, we obtained spectral maps that present the spectrum (in velocity units) at each pixel of a determined region around the center of the source. These maps allow us to evaluate the spectra at different offsets from the center of the envelope where the dust continuum emission peak localizes the protostar, and to detect any distinct line profile or velocity pattern. The obtained spectral maps are discussed in Appendix B.1.

We derived the integrated intensity maps by integrating the molecular line emission over the velocity range in which it is emitting (the velocity range used for each molecule is shown in Table 1). Figure 1 shows the dust continuum emission map at 110 GHz and the integrated intensity contours of C¹⁷O (J = 1–0), DCO⁺ (J = 3–2), H¹³CO⁺ (J = 1–0), and H¹³CO⁺ (J = 3–2). The H¹³CO⁺ (J = 1–0) emission clearly appears more spatially extended than the other lines, roughly following most of the dust emission. The emission from this line extends further than the dust toward the southeast outflow cavity wall direction. The C¹⁷O (J = 1–0), DCO⁺ (J = 3–2), and H¹³CO⁺ (J = 3–2) emission present a similar, relatively compact morphology centered around the dust peak and elongated along the north-south equatorial plane. However, the spatial extent along the equatorial plane is slightly more compact for the H¹³CO⁺ emission (~2100 au) than the DCO⁺ emission (~2500 au), and both are more extended than C¹⁷O (~1700 au). The DCO⁺ peak intensity is clearly more displaced from the continuum peak than the H¹³CO⁺ (J = 3–2) peak; this could be a local effect around the protostar due to the large temperature and irradiation conditions, which destroy DCO⁺ (Wootten 1987; Butner et al. 1995). The C¹⁷O (J = 1–0) emission appears to be slightly extended along the outflow cavities.

As all the molecules present a double-peaked profile with similar velocity patterns to the ones observed in Cabedo et al. (2021), namely two different velocity components with varying intensity across the source spatial extent, we use two independent velocity components to fit the line profiles simultaneously. Following Cabedo et al. (2021), we modeled the spectrum at each pixel using the HfS fitting program (Estalella 2017). We obtained maps of peak velocity and velocity dispersion for each velocity component, deriving the corresponding uncertainties from the χ² goodness of fit (see Estalella 2017, for a thorough description of the derivation of uncertainties). A more detailed discussion of these maps is presented in Appendix B.2. In addition, an estimate of line opacities is given in Appendix C. In the following sections, the statistical uncertainties derived (shown in Figs. 2, 4, and 5) are computed using standard error propagation analysis from the uncertainties of each parameter obtained with the line modeling (peak velocity, velocity dispersion, and opacity). Other systematic uncertainties are discussed separately in the corresponding sections.

3 Deuteration fraction

The process of deuteration consists of an enrichment of the amount of deuterium with respect to hydrogen in molecular species. The deuteration fraction, R_D,= [D]/[H], in molecular ions, in particular HCO⁺, has been extensively used as an estimator of the degree of ionization in molecular gas (Caselli et al. 1998; Fontani et al. 2017). Here, we apply this method to obtain maps of R_D, computed as the column density ratio of DCO⁺ (J = 3–2) and its non-deuterated counterpart, H¹³CO⁺ (J = 3–2):

$$R_{\text{D}}=\frac{1}{f_{{}^{12/13}\text{C}}}\frac{N_{{\text{DCO}}^{+}}}{N_{{\text{H}}^{13}{\text{CO}}^{+}}},$$(1)

where N_i is the column density of each species and $f_{{}^{12/13}\text{C}}$ is the abundance ratio of ¹²C to ¹³C. The column density of both species is computed such as:

$$\begin{array}{lll}{N}_i\hfill & =\hfill & \frac{8\pi }{{\lambda }^{3}A}\frac{1}{J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)}\hfill \\ \hfill & \hfill & \times \frac{1}{1-\exp \left(-h\nu \text{/}{k}_{\text{B}}{T}_{\text{ex}}\right)}\frac{Q_{\text{rot}}}{g_u\exp \left(-E_l\text{/}{k}_{\text{B}}{T}_{\text{ex}}\right)}\hfill \\ \hfill & \hfill & \times {\displaystyle \int I_{0}dv,}\hfill \end{array}$$(2)

where λ is the wavelength of the transition, A is the Einstein coefficient, J_v(T) is the Planck function at the background (2.7 K) and at the excitation temperature of the line (assumed to be equal to the dust temperature at the observed scales, T_ex = T_d = 25 K, Galametz et al. 2019; Cabedo et al. 2021), k_B is the Boltzmann constant, g_u is the upper state degeneracy, Q_rot is the partition function at 25 K, E_l is the energy of the lower level, and ∫I₀dv is the integrated intensity. Values of these parameters for each molecule are listed in Table 2.

The top panels of Fig. 2 show the R_D maps of B335 for the blue- and redshifted velocity components (left and right column, respectively). The mean values of R_D are ~ 1–2%, being higher in the outer region of the envelope where the gas is expected to be colder, and decreasing toward the center where the protostar is located and the temperature is expected to rise. The bottom panels show that the statistical uncertainties are roughly one-tenth of the deuteration values. We note that the highest values of R_D are found for the blueshifted component, as high as 3.5%. However, these values have the largest associated errors due to the presence of a third velocity component in the line profiles of H¹³CO⁺ (J = 3–2) and should not be totally trusted (see Appendix B.1). Finally, we note that while the average R_D values for both velocity components are similar, the most widespread blueshifted gas component shows larger R_D dispersion, with values from 0.06% to 2%, while the more localized redshifted gas exhibits more uniform values, between 0.01% and 1%.

One of the main uncertainties of the deuteration analysis is the assumption that both DCO⁺ (J = 3–2) and H¹³CO⁺ (J = 3–2) have the same T_ex. This should be accurate enough, because both molecules present a similar geometry and dipolar moment and are assumed to be spatially coexistent. However, to check the effect that the variation of T_ex would produce in our values, we computed the column densities of both species at two additional temperatures, 20 and 30 K. While the DCO⁺ column density does not change significantly at these two temperatures with respect to the value adopted (25 K), the H¹³CO⁺ column density changes by about 15%. Nevertheless, the new values are still in agreement with the derived values of R_D from DCO⁺ reported in the literature (see Sect. 5.1).

Fig. 1 Map of the dust continuum emission at 110 GHz clipped for values smaller than 3σ (see Table 1) and integrated emission for the four molecular lines as indicated in the upper left corner (red contours). Contours are −2, 3, 5, 8, 11, 14, and 17σ (the individual σ values are indicated in Table 1). The velocity ranges used to obtain the integrated intensity maps are shown in Table 1. In the bottom-left of each plot, we show the beam size for each molecular line (in red) and for the dust continuum emission (in black). The physical scale is shown in the top-right corner of the figure.

Fig. 2 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the deuteration fraction maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The spatial scale is shown in the top-right corner of the upper panels.

4 Ionization fraction and cosmic-ray ionization rate

We followed the method presented in Caselli et al. (1998) to compute χ_e and ζ_CR from R_D. The two quantities are given by

$${\chi }_{\text{e}}=\frac{2.7\times {10}^{-8}}{R_{\text{D}}}-\frac{1.2\times {10}^{-6}}{f_{\text{D}}},$$(3)

and

$${\zeta }_{\text{CR}}=\left[7.5\times {10}^{-4}{\chi }_{\text{e}}+\frac{4.6\times {10}^{-10}}{f_D}\right]{\chi }_{\text{e}}{n}_{{\text{H}}_2}{R}_{\text{H}},$$(4)

where f_D is the depletion fraction of C and O, $n_{{\text{H}}_2}$ is the H₂ volume density, and R_H is the hydrogenation fraction, R_h = [HCO⁺]/[CO].

While this observational method to measure χ_e of the gas relies on a simple chemical network and a rather strong hypothesis, it is widely used in the literature in dense cores and also toward protostars. The validity and shortcomings of the chemical models used to constrain ionization fractions are concisely discussed in Sect. 5.3.

4.1 Estimate of the CO depletion in B335

As R_D, and therefore χ_e, depend on the level of C and O depletion from the gas phase, we hypothesize that the CO abundance is proportional to the CO column density and estimated f_D as the ratio between the “expected”, $N_{\text{CO}}^{\exp }$, and the “observed” CO column density, $N_{\text{CO}}^{\text{obs}}$:

$$f_D=\frac{N_{\text{CO}}^{\exp }}{N_{\text{CO}}^{\text{obs}}}.$$(5)

Here, $N_{\text{CO}}^{\exp }$ is computed as the product of the H₂ column density, $N_{{\text{H}}_2}$, and the expected CO to H₂ abundance ratio, X_CO = [CO]/[H₂] = 10⁻⁴ (Wilson & Rood 1994; Gerner et al. 2014), and NO is derived as the product of the C¹⁷O column density, $N_{{\text{C}}^{17}\text{O}}$, and the C¹⁷O to CO abundance ratio, $f_{{\text{C}}^{17}\text{O}}=\left[\text{CO}\right]/\left[{\text{C}}^{17}\text{O}\right]=2317$ (Wouterloot et al. 2008). Then, f_D is given by

$$f_D=\frac{N_{{\text{H}}_2}{X}_{\text{CO}}}{N_{{\text{C}}^{17}\text{O}}{f}_{{\text{C}}^{17}\text{O}}}.$$(6)

The H₂ column density has been estimated from the dust thermal emission at 110 GHz (see Fig. 1) corrected for the primary beam attenuation. Indeed, $N_{{\text{H}}_2}$ depends directly on the flux density measured on the continuum map, the beam solid angle Ω_beam, the Planck function B_ν(T_d) at the dust temperature, and the dust mass opacity κ_ν (in units of cm² g⁻¹) at the frequency at which the dust thermal emission was observed:

$$N_{{\text{H}}_2}=-\frac{1}{{\mu }_{{\text{H}}_2}{m}_{\text{H}}{\kappa }_v}\ln \left[1-\frac{S_v^{\text{beam}}}{{\text{Ω}}_{\text{beam}}{B}_v\left(T_d\right)}\right],$$(7)

where ${\mu }_{{\text{H}}_2}=2.8$, m_H is the atomic hydrogen mass, and for the dust opacity we used a power-law fit given by

$${\kappa }_v=\frac{{\kappa }_0}{{\chi }_{\text{d}}}{\left(\frac{v}{v_0}\right)}^{\beta }.$$(8)

We adopt a dust mass absorption coefficient κ₀ = 1.6 cm² g⁻¹ at λ₀ = 1.3 mm, following Ossenkopf & Henning (1994). We assume a standard gas-to-dust ratio χ_d = 100, a dust emissivity exponent β = 0.76 (Galametz et al. 2019), and a dust temperature T_d = 25 K. The resulting $N_{{\text{H}}_2}$ column densities probed in the ALMA dust continuum emission map range from ~10²⁰ cm⁻² at radii of ~1600 au up to a few 10²² cm⁻² at the protostar position. We note that the values of column density suffer from a systematic error due to assumptions on the parameters used to derive them (e.g., dust opacity, dust temperature). We estimate that the systematic error is about a factor of 3–5, being higher toward the peak position where standard dust properties and conditions of emission may not be met. We use Eq. (2) to estimate the observed CO column density $N_{\text{CO}}=N_{{\text{C}}^{17}\text{O}}{f}_{{\text{C}}^{17}\text{O}}$ from the C¹⁷O emission map shown in Fig. 1, assuming an excitation temperature of 25 K. The CO column density values are subject to systematic uncertainties due to the adopted gas temperature, as this latter is not well constrained and is expected to vary across the envelope. Altogether, the systematic uncertainty is about a factor of 2, and is higher toward the peak position where different gas temperatures are expected to be present along the line of sight.

Using Eq. (6), we obtain the CO depletion fraction map shown in Fig. 3. The depletion factor ranges from ~20 in the outer regions to ~70 in the center of the object. Propagating the errors from the column density maps they stem from, we expect these depletion values to have an uncertainty of a factor of 2–3, as the main source of errors – which are due to the assumed gas and dust temperatures – are lessened by the ratio. The CO depletion seems highly asymmetric around the protostar, being the lowest in the southeastern quadrant, while high values are associated to the protostar position and the northern region. We note that the regions with high f_D coincide with the regions of low deuteration shown in Fig. 2.

Fig. 3 Map of the CO depletion factor and contours of dust continuum emission at 110 GHz, for emission over 3σ. The spatial scale is shown in the top-right corner of the figure.

4.2 Ionization fraction of the gas

We obtained χ_e using our R_D maps, shown in Sect. 3 and the CO depletion map, presented in Fig. 3. The top panels of Fig. 4 show the derived χ_e and the bottom panels the associated uncertainties, for the blue- and red-shifted velocity components (left and right column, respectively). We obtained a mean value of χ_e = 2 × 10⁻⁶ for both components. Generally, these values are subject to statistical uncertainties of ~1.5 × 10⁻⁶. The range of values found for R_D and f_D indicates that the χ_e depends mostly on the level of deuteration, and therefore the systematic errors on the estimated χ_e from Eq. (3) are typically < 30% and statistical uncertainties dominate in this case. As R_D decreases toward the center of the source, we find that χ_e increases toward the center.

4.3 Derivation of the CR ionization rate

We produced ζ_CR maps using Eq. (4). We obtained $n_{{\text{H}}_2}$ as $N_{{\text{H}}_2}\text{/}L$, using the $N_{{\text{H}}_2}$ map obtained in Sect. 4.1 and considering a core diameter, L, of 1720 au derived from the C¹⁷O (J = 1–0) emission (Cabedo et al. 2021). With this approximation, we find values of the volume density ranging from ~10⁴ in the outer regions, to ~10⁶ cm⁻³ at the protostar position. We obtained R_h for both velocity components as the column density ratio of H¹³CO⁺ (J = 1–0) to C¹⁷O (J = 1–0), and accounting for $f_{{}^{12/13}\text{C}}$ and $f_{{\text{C}}^{17}\text{O}}$:

$$R_{\text{H}}=\frac{N_{{\text{H}}^{13}{\text{CO}}^{+}}{f}_{{}^{12/13}C}}{N_{{\text{C}}^{17}\text{O}}{f}_{{\text{C}}^{17}\text{O}}}.$$(9)

Values of R_H are relatively uniform across the source, with mean values between 2 and 3 × 10⁻⁷ for the two velocity components. The statistical uncertainties of these values are generally one order of magnitude smaller than R_H.

The two terms within the brackets in Eq. (4) have values of the order of ~10⁻⁹ and ~10⁻¹¹, respectively. Therefore, the derivation of ζ_CR is very sensitive to χ_e but not to f_D. The top panels of Fig. 5 show the derived ζ_CR and the bottom panels show the associated statistical uncertainties, for the blue- and redshifted velocity components (left and right column, respectively). As expected from the values of χ_e, ζ_CR increases toward the center, reaching values of 7 × 10⁻¹⁴ s−¹ at the peak of the dust continuum. The uncertainties on ζ_CR are large, due to all the errors coming from the modeling of the different molecules. Considering the systematic uncertainties on $n_{{\text{H}}_2}$, and R_H we estimate the global uncertainty on ζ_CR values, and find them to be of the same order of magnitude as the ζ_CR values.

To interpret the trend of ζ_CR as a function of the distance from the density peak in an unbiased way, we considered profiles of ζ_CR along directions from the density peak outward denoted by the position angle ϑ, with 0 ≤ ϑ ≤ π (see the lower right panel of Fig. 6). The envelope of the profiles is represented by the orange filled region in the upper left panel of Fig. 6. As the ζ_CR(r, ϑ) distributions between 40 and 700 au are skewed, for each radius $\overline{r}$ we considered the median value of ζ_CR(r, ϑ) and estimated the errors using the first and third quartiles. We find that the trend of ζ_CR can be parameterized by two independent, decreasing power-law profiles ζ_CR(r) ∝ r^s. At small radii, r < 270 au, the slope of the power law, s = −0.96 ± 0.04, is compatible with the diffusive regime of CR propagation. However, at r > 270 au, ζ_CR drops abruptly to s = −3.77 ± 0.30, which cannot be explained by any propagation regime. The possible explanations for this result are explained with more details in Sect. 5.3.

Fig. 4 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the ionization fraction maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The spatial scale is shown in the top-right corner of the upper panels.

5 Discussion

5.1 Local destruction of deuterated molecules

The deuteration fractions observed in Class 0 protostars range between 1% and 10%, with a large scatter from source to source and at the different protostellar scales (Caselli 2002; Roberts et al. 2002; Jørgensen et al. 2004). The [DCO⁺]/[H¹³CO⁺] values that we find for B335 are therefore in broad agreement with the typical range reported in the literature, even accounting for temperature effects. However, Butner et al. (1995) found a [DCO⁺]/[HCO⁺] of ≃3% in B335 at scales of ~10000 au. This value is larger than the range of values we find with ALMA at smaller scales, which indicates a decrease in R_D toward the center of the envelope. Indeed, R_D decreases down to < 1% toward the center and the northern region of the source. This decrease could be attributed to local radiation processes occurring during the protostellar phase, after the protostar has been formed. An increase in the local radiation could promote processes that would lead to the destruction of deuterated molecules due to the evaporation of neutrals to the gas phase (Caselli et al. 1998; Jørgensen et al. 2004), or the increased abundance of ionized molecules like HCO⁺ (Gaches et al. 2019).

The decreasing trend in R_D is further confirmed by our observations of N₂D⁺ (J = 3–2). Figure 7 shows the deuteration map overplotted with the N₂D⁺ (J = 3–2) integrated emission. A lack of N₂D⁺ is observed in regions where RD is lower, suggesting that the physical processes lowering the abundance of O-bearing deuterated molecules also lower the abundance of N-bearing deuterated species around the protostar.

Fig. 5 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the ionization rate maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The scale is shown in the top-right corner of the upper panels.

5.2 High values of C and O depletion in the source

At core scales, observations of the depletion in B335 have been rather inconclusive as some studies suggest very little CO depletion (JCMT 20-arcsec beam probing 3000 au, Murphy et al. 1998), while others suggest a CO depletion of one order of magnitude, f_D ~ 10, at similarly large scales (Walmsley & Menten 1987). Our spatially resolved observations reveal large values of CO depletion of the gas in the inner envelope of B335 at radii <1000 au. As shown in Fig. 3, we find f_D ranging from 20 to 40 at radii of 1000–100 au, except at the core center where it increases up to 80. Such CO depletion at protostellar radii 100–1000 au is somehow unexpected, if one assumes depletion is due to freeze-out onto dust grains. Indeed, gas and dust at such radii should have temperatures beyond ~20–30 K, where the evaporation of CO from the surface of the dust should replenish it in the gas phase (Anderl et al. 2016). The high depletion values found at the protostar position are highly uncertain due to unknown optical depth, and unconstrained dust and gas temperatures. We stress that the possible depletion increase toward the center in B335 should be confirmed with further observations of higher J transitions of CO.

Several observational studies of protostellar cores have found depletion factors f_D > 10 at scales of ~ 5000 au, with some spatially resolved studies pointing to values up to ~ 20 (Alonso-Albi et al. 2010; Christie et al. 2012; Yildiz et al. 2012; Fuente et al. 2012). Models of protostellar evolution with short warm-up timescales in the envelope (see e.g., Aikawa et al. 2012) show that f_D may remain high even inside the sublimation radius, especially at early stages, because of the conversion of CO into CH₃OH and CH₄ on the grain surfaces. Photodissociation due to the protostellar UV radiation could contribute to lowering the global CO abundance in the innermost layers of protostellar envelopes (Visser et al. 2009). Additionally, other sources of ionization, such as CRs, could decrease the abundance of CO in the gas phase by promoting its reaction for the formation of HCO⁺ (Gaches et al. 2019). This decrease in CO by transformation into HCO⁺ could also explain the observed coincidence between regions with large CO depletion and low R_D.

Another possible cause for the observed depletion could be the presence of relatively large dust grains, as observed in several protostars at similar scales by Galametz et al. (2019). Indeed, larger grains are less efficiently warmed up (Iqbal & Wakelam 2018), limiting the desorption of CO ices to the gas phase. We note that the presence of large dust grains in the inner envelope of B335 could also favor enhanced ionization of the gas (Walmsley et al. 2004), and thus also be consistent with our findings regarding ionization.

Fig. 6 Upper left panel: ζCR as a function of the radius, r. The two dotted black lines at 270 and 500 au identify the radius ranges used to compute the slopes (s1 and s2, respectively). The gray shaded region shows the expected X-ray ionization in the outflow cavity for a typical and a flaring T Tauri star (Rab et al. 2017). Lower left panel: hydrogen column density, $N_{{\text{H}}_2}$, from the source center outwards as a function of radius. Upper right panel: ζCR as a function of the H2 column density. The upper x-axis shows the minimum energy that CR protons must have in order not to be thermalized (Padovani et al. 2018). The orange, magenta, and green filled regions in the three upper panels represent the envelope of the ζCR and $N_{{\text{H}}_2}$ profiles (see Sect. 5.2), while the red, purple, and dark green solid lines show their median value. Lower right panel: ζCR map of the blueshifted component (zoom of the upper left panel of Fig. 5) in logarithmic scale (colored map). The purple hatched region shows a circle of radius 270 au, where CRs propagate according to the diffusive regime. For illustration purposes, the black arrow at the position angle # shows a direction used to extract the ζCR profile from the map. The black dashed lines identify the radius as well as the corresponding values of $N_{{\text{H}}_2}$ and ζCR, where the slope of ζCR changes from s1 to s2.

Fig. 7 Deuteration fraction maps for the blue- and redshifted velocity components (left and right columns, respectively). The N2D+ (J = 3–2) integrated emission is shown overlaid as black contours at 3, 5, 10, 15, and 20σ, where σ = 5.23 mJy beam−1. The spatial scale is shown in the top-right corners of each panel.

5.3 Origin of the ionization

From the single-dish observations of Butner et al. (1995), Caselli et al. (1998) find ionization fractions of protostellar gas at core scales of the order of 10⁻⁸–10⁻⁶. In the B335 outer envelope where n(H₂) ~ 10⁴ cm⁻³, these authors estimate χ_e of a few 10⁻⁶. With our interferometric data, we measure χ_e ~2 × 10⁻⁶ in the inner envelope, where n(H₂) ~ 10⁶ cm⁻³, with the largest values being χ_e ~ 9 × 10⁻⁶. We note that while all these aforementioned values are only providing line-of-sight averages for χ_e, the interferometric data only select a narrow range of spatial scales, making them less sensitive to line-of-sight averaging than single-dish data, and therefore they are more representative of local abundance variations.

We stress that, while the use of deuteration to measure χ_e is a standard practice in the literature, it has been noted that such an approach may not give accurate results in low-density medium experiencing small deuteration fractions (e.g., ~10³ cm⁻³, Vaupré et al. 2014). In this case, the abundances of DCO⁺ and HCO⁺ mostly depend on the ζ_CR/n_H ratio and not only on ζ_CR, and using R_D to get ζ_CR is valid only when deuteration is > 2%: the measured values in B335 are mostly all lower than this threshold at radii <500 au, which would suggest that our approach of using these observations to derive χ_e is deficient. However, recent models have shown that in warm dense gas (T ≳ 24 K and n_H ~ 10⁶ cm⁻³), which is typical of the conditions probed here, the relation between R_D and ζ_CR is more robust. Shingledecker et al. (2016) show that the scatter of ζ_CR predicted by the Caselli model at low ionization values seems to be related to variations in the ortho-to-para ratio (OPR) of H₂. While the OPR is an unobservable parameter for dense gas conditions, and therefore cannot be measured for our source, the work of these latter authors shows that the dependence of the relation between R_D and ζ_CR on OPR is weaker for evolved sources due to the stabilization of the OPR. This effect has been shown to happen sooner when n_H and the temperature are increased. This was also pointed out by Bron et al. (2021) who show that R_D is a robust estimator of ζ_CR when the gas is subject to high ionization fractions (χ_e ≳ 3 × 10⁻⁷). As we are probing a dense region of the core with temperatures larger than 20 K and ionization from an embedded source, this supports our use of the Caselli method to measure ζ_CR.

The large χ_e values could be produced by local UV and X-ray radiation produced in outflow shocks, altering the chemistry from the molecular content around the shocks (e.g., Viti & Williams 1999; Girart et al. 2002, 2005). However, ionization by far-UV (FUV) photons can be ruled out. Indeed, the maximum energy of FUV photons is around 13 eV, thus below the threshold for ionization of molecular hydrogen (15.44 eV). The second argument is related to the fact that the extinction cross-section at FUV wavelengths (≃0.1 μm) is σ_UV ≃ 2 × 10⁻²¹ cm² per hydrogen atom (Draine 2003). Thus, FUV photons are rapidly absorbed in a thin layer of column density equal to (2σ_UV)⁻¹ ≃ 3 × 10²⁰ cm⁻².

The absorption cross-section of X-ray photons, σ_X, in the range 1–10 keV is much smaller than at FUV frequencies (Bethell & Bergin 2011). For example, at 1 keV and 10 keV, σ_X ≃ 2 × 10⁻²² cm² and ≃8 × 10⁻²⁵ cm² per H atom, respectively.

While at 1 keV the corresponding absorption column density is still small, (2σ_X)⁻¹ ≃ 2 × 10²¹ cm⁻², at 10 keV it reaches about 6 × 10²³ cm⁻², much larger than the maximum H2 column density, $N_{{\text{H}}_2}$, that we found. The latter was computed by integrating the dust continuum at 110 GHz along directions identified by the position angle ϑ (see lower right panel of Fig. 6) and is representative of the column density in the outflow cavity. We therefore estimated the X-ray ionization rate, ζ_X, using the spectra described by Rab et al. (2017) for a typical and a flaring T Tauri star, assuming a stellar radius of 2R_⊙. Results for ζ_X are shown in the upper left panel of Fig. 6. It is evident that X-ray ionization cannot explain the values of the ionization rate estimated from observations. A fraction of these X-ray photons end up in the disk, which has much higher densities than the cavity, and so $N_{{\text{H}}_2}$ can easily reach values greater than 10²⁴ cm⁻², as found by Grosso et al. (2020) for a Class 0 protostar. At such high $N_{{\text{H}}_2}$, ζ_X decreases dramatically.

Once the high-energy radiation field is excluded as the main source of ionization rate estimated from observations, the only alternative left is CRs (see Sect. 4.3). However, the range of derived values is much higher than those expected from Galactic CRs (see Appendix C in Padovani et al. 2022, for an updated compilation of observational estimates of ζ_CR). Following McKee (1989), χ_e = 1.3 × 10⁻⁵n(H₂)^−1/2, and for the volume densities derived in Sect. 4.3, we should expect χ_e values of 10⁻⁷ to 10⁻⁸, which are between one and two orders of magnitude lower than the observed values in B335. As we also observed a significant increase in χ_e toward the center, and especially of ζ_CR, the origin of the high values should be of local origin. In particular, we are most likely witnessing the local acceleration of CRs in shocks located in B335 as predicted by theoretical models (Padovani et al. 2015, Padovani et al. 2016, 2020). The values of ζ_CR found in B335, between about 10⁻¹⁶ and 10⁻¹⁴ s⁻¹, are among the highest values reported in the literature toward star-forming regions (e.g., Maret & Bergin 2007; Ceccarelli et al. 2014b; Fuente et al. 2016; Fontani et al. 2017; Favre et al. 2017; Bialy et al. 2022). We stress that this is the first time that ζ_CR is measured at such small scales in a solar-type protostar and can be attributed to locally accelerated CRs, as well as the first time that a map of ζ_CR is obtained.

There are two possible origins for the production of local CRs close to a protostar. One is in strong magnetized shocks along the outflow (Padovani et al. 2015, 2016, 2021; Fitz Axen et al. 2021). Indeed, synchrotron radiation, which is the signature of the presence of relativistic electrons, has been detected in some outflows (Carrasco-González et al. 2010; Ainsworth et al. 2014; Rodríguez-Kamenetzky et al. 2016, 2019; Osorio et al. 2017). The second possibility is in accretion shocks near the stellar surface (Padovani et al. 2016; Gaches & Offner 2018). Both mechanisms are expected to generate low-energy CRs (≲ 100 MeV–30 GeV) through the first-order Fermi acceleration mechanism with a rate of up to ~ 10⁻¹³ s⁻¹.

The ionization rate slope found in B335 for r < 270 au is close to −1 and compatible with a diffusive regime, in agreement with predictions of theoretical models. Surprisingly, at radii larger than ~ 250 au, the slope decreases below −2, thus beyond the geometrical dilution regime. We propose two possible explanations. On the one hand, local CRs may have accumulated enough column density to start being thermalized. On the other hand, at radii above about 250 au, both the dust emission and ζ_CR maps gradually lose their central symmetry (see Figs. 1 and 5, respectively). For example, the continuum map shows two higher density arms toward the northeast and southeast. Therefore, the slope, calculated by averaging over the position angle distribution ϑ, could be less than −2. The loss of symmetry at larger radii in the ionization maps might also be additional evidence as to the importance of non-symmetrical motions during the protostellar collapse (Cabedo et al. 2021). However, we note that these values are also affected by our limited angular resolution of 1.5 arcsec (~250 au), implying that only three or four points in the plot are completely independent, and that observations at larger angular resolution are needed in order to confirm the observed trend.

Previous observations indicate that B335 hosts a powerful and variable jet, which might explain a large production of local CRs (Gâlfalk & Olofsson 2007; Yen et al. 2010). Nevertheless, the exact origin of the local CR source is difficult to pinpoint due to our limited angular resolution. The column density calculated from the dust continuum map at 110 GHz (see lower left panel of Fig. 6) is consistent with that expected in the outflow cavity. To calculate the minimum energy that protons must have in order to pass through a given column density without being thermalized, we can make use of the stopping range function (see Fig. 2 in Padovani et al. 2018). At $N_{{\text{H}}_2}\simeq {10}^{22}$ cm⁻², protons with energies of about 10 MeV are thermalized and lose their ionizing power (see also the upper x-axis of Fig. 6). As discussed above, if the shock is in the vicinity of the protostellar surface, a column density of the order of 10²⁴–10²⁵ cm⁻² can easily be accumulated (Grosso et al. 2020). In this case, the minimum proton energy required to prevent thermalization is 100 MeV and 400 MeV at 10²⁴ and 10²⁵ cm⁻², respectively. Models of local CR acceleration in low-mass protostars predict that accelerated protons can reach maximum energies of the order of 100 MeV and 10 GeV if the shock in which they are accelerated is located in the jet or close to the protostellar surface, respectively (Padovani et al. 2015, 2016). Thus, although the position of the shock that accelerates these local CRs in B335 is unknown, models suggest that the maximum energies of the protons are sufficiently high to explain the observed ionization.

Additionally, we note that B335 has been observed to exhibit an organized magnetic field at scales similar to the ones probed here (Maury et al. 2018). This could also create favorable conditions for enhanced CR production. Finally, our results are in agreement with CR ionization models on protostellar cores with an embedded radiation source (Gaches et al. 2019); these latter models also suggest that ionization models based solely on ${\text{H}}_3^{+}$ chemistry might underpredict χ_e by one order of magnitude, except in regions where CR ionization is dominant.

5.4 Magnetic field coupling and non-ideal MHD effects

The ionization fraction of the gas in a young accreting protostar is not only important to understand the early chemistry around solar-type stars, feeding the scales where disks and planets will form, but is also the critical parameter that determines (i) the coupling between the magnetic field and the infalling-rotating gas in the envelope, and (ii) the role of diffusive processes, such as ambipolar diffusion, which counteract the outward transport of angular momentum due to magnetic fields (a process known as magnetic braking). The higher the gas ionization, the more efficient the coupling, and the braking.

The large fraction of ionized gas unveiled by our observations in the inner envelope of B335 should lead to almost perfect coupling of the gas to the local magnetic field lines, generating a drastic braking of rotational motions. We note that these new results lend extra support to other observational evidence of a very efficient magnetic braking at work in B335, notable the highly pinched magnetic field lines observed at similar scales by Maury et al. (2018), and the failure to detect any disk larger than ~ 10 au in this object (Yen et al. 2015); although a new kinematic analysis may be motivated in light of the detection of several velocity components in the accreting gas at scales of ~500 au (Cabedo et al. 2021).

At the end of the pre-stellar stage, the fraction of ionized gas toward the central part of the core is expected to be very low, as Galactic CRs do not penetrate deeply and there is no local ionization source (Padovani & Galli 2013; Ceccarelli et al. 2014a; Silsbee & Ivlev 2019). Hence, the initial stages of protostellar evolution are proceeding under low ionization conditions at typical ζ_CR < 10⁻¹⁶ s⁻¹ (Padovani et al. 2018; Ivlev et al. 2019), which enhance the importance of non-ideal MHD processes (Padovani et al. 2014), with efficient diffusion of magnetic flux outwards during the very first phases of the collapse, and reduced magnetic braking. If the local ionization processes we observe in B335 are prototypical of solar-type Class 0 protostars, then CRs accelerated in the proximity of the protostar could be responsible for changing the ionization fraction of the gas at disk-forming scales once the protostar is formed, setting quasiideal MHD conditions in the inner envelope. The timescale and magnitude of this transition from non-ideal MHD conditions to quasi-ideal MHD conditions may depend on the protostellar properties: more detailed modeling and observations toward other protostars are required to address this question. Moreover, we note that the observed large χ_e is not in agreement with the values used to calculate the simplified chemical networks in non-ideal MHD models of protostellar formation and evolution (Marchand et al. 2016; Zhao et al. 2020), as gas ionization is usually predicted from the gas density following McKee (1989). Thus, our observations suggest that careful treatment of ionizing processes in magnetized models may be crucial to properly describe the gas-magnetic field coupling at different scales in embedded protostars.

In this scenario, the properties of protostellar disks could be tightly related to the local acceleration of low-energy CRs, and the development of a highly ionized region around the protostar. Therefore, we would not expect the ionization fraction of the gas present at large scale in the surrounding cloud to be a key factor in setting the disk properties, as proposed for example by Kuffmeier et al. (2020). The scenario we propose is also in agreement with recent ALMA observations of the Class II disks in Orion that do not find supporting evidence of local cloud properties affecting the disk properties (van Terwisga et al. 2022).

6 Conclusions and summary

This work provides new insights into the physico-chemical conditions of the gas in the young Class 0 protostar B335. For the first time, we derived a map of the gas χ_e and of ζ_CR at envelope scales < 1000 au in a Class 0 protostar. Our results highlight the interplay between physical processes responsible for gas ionization at disk-forming scales, and its consequences for magnetized models of solar-type star formation. Here, we summarize the main results of our analysis:

• We used ALMA data to create molecular line emission maps and used these to characterize the physico-chemical properties of the dense gas in B335. We used the chemical recipe from Caselli et al. (1998) to derive χ_e based on the deuteration fraction of HCO+ and of the f_D. This recipe was recently shown to be valid for dense molecular gas, ~10⁶ cm⁻³, with a high ionization fraction (≳3 × 10⁻⁷, Bron et al. 2021). It allows us to characterize the ionization of the gas at envelope radii ≲ 1000 au: we find large fractions of ionized gas, χ_e, between 1 and 8 × 10⁻⁶. These values are remarkably higher than the ones usually measured at core scales.

• We produced for the first time a map of ζ_CR, in the close environment of a young embedded protostar. Our analysis reveals very high values of ζ_CR between 10⁻¹⁶ and 10⁻¹⁴ s⁻¹, increasing toward the central protostellar embryo at small envelope radii. This suggests that local acceleration of CRs, and not the penetration of interstellar CRs, may be responsible for the gas ionization in the inner envelope, potentially down to disk forming scales.

• The large fraction of ionized gas we find suggests an efficient coupling between the magnetic field and the gas in the inner envelope of B335. This interpretation is also supported by the observations of highly organized magnetic field lines, and no detection of a large rotationally supported disk in B335.

If our findings are found to be prototypical of the low-mass star formation process, they might imply that the collapse at scales < 1000 au transitions from non-ideal to a quasi-ideal MHD once the central protostar starts ionizing its surrounding gas, and very efficient magnetic braking of the rotating-infalling protostellar gas might then take place. Protostellar disk properties may therefore be determined by local processes setting the magnetic field coupling, and not only by the amount of angular momentum available at large envelope scales and by the magnetic field strength in protostellar cores. We stress that the gas ionization we find in B335 significantly stands out from the typical values routinely used in state-of-the-art models of protostellar formation and evolution. Our observations suggest that a careful treatment of ionizing processes in these magnetized models may be crucial to properly describe the processes responsible for disk formation and early evolution. We also note that more observations of B335 at higher spatial resolution, and of other Class 0 protostars, are crucial to confirm our results.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme (MagneticYSOs project, grant agreement N° 679937). This work was also partially supported by the program Unidad de Excelencia Maria de Maeztu CEX2020-001058-M. J.M.G. also acknowledges support by the grant PID2020-117710GB-I00 (MCI-AEI-FEDER, UE). Additionally, this paper makes use of the following ALMA data: 2015.1.01188.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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

Appendix A Technical details of the ALMA observations

Table A.1 presents the technical details of the ALMA observations presented in this work, including the ALMA configuration, the date and time range, the central frequency, the obtained spectral resolution, and the different calibrators used.

Appendix B Line profile analysis

B.1 Spectral maps

The DCO⁺ (J=3-2) and H¹³CO⁺ (J=3-2) spectral maps are presented in Figs. B.1 and B.3. Both maps show emission in a region of 5.5″ × 5.5″ (900 au) around the center of the source, with each spectrum corresponding to a pixel of 0.5″ (82 au). The H¹³CO⁺ (J=1-0) shows only a region 4″ × 4″ (685 au) and each spectrum corresponds to 0.25″ (41 au).

The three tracers present the same pattern observed in C¹⁷O (J=1-0), presented in Cabedo et al. (2021). The line profiles are double-peaked with one of the two velocity components progressively disappearing at different offsets of the source, that is, the blueshifted component is dominant in the eastern region of the source while the redshifted component is present only in the western region. For the three transitions, the blueshifted component is always more intense than the redshifted one, which produces the blue asymmetric morphology observed in the integrated intensity images.

As none of the tracers present a morphology typical of the outflow cavity, we assume that they are not tracing outflowing gas and that the two velocity components are not caused by this effect. This is confirmed with the spectral maps where no signatures coming from outflowing gas —such as very broad lines or high velocity wings— can be observed in the line profiles. This is generally true, except for a small region on the southeast region of the H¹³CO⁺ (J=3-2) map, where an additional third velocity component, or “shoulder”, can be observed. This additional component appears to be very broad, and at velocities of ~7.6–7.7 km s⁻¹. We do attribute this signature to outflow contamination, which is observed neither in DCO⁺ (J=3-2) nor in H¹³CO⁺ (J=1-0). We consider that while the outflow might have some effect in the line profiles, this is negligible and is not the main cause of the double-peaked line profiles. The effect of large-scale filtering cannot be ruled out. This is particularly visible in the DCO⁺ (J=3-2) and H¹³CO⁺ (J=3-2) where negative emission is observed in the center and west regions of the source. While these two effects can be a source of uncertainty in our analysis, we conclude that they alone cannot produce the observed line profiles and the velocity pattern through the source, and proceed with the line modeling with the assumption of two independent velocity components.

B.2 Line profile modeling

In this case, as none of the transitions present a resolved hyperfine structure, the opacity of the line could not be obtained by modeling the line profile and only the peak velocity and velocity dispersion maps (as σ) are presented. The resulting maps for DCO⁺ (J=3-2), H¹³CO⁺ (J=1-0), and H¹³CO⁺ (J=3-2) are shown in Figs. B.4, B.5, and B.6. For the three tracers, the two separated velocity components occupy different regions of the source, indicating that they are probing gas with different motions. The morphology of the two components is similar for the three tracers, but it is slightly different from the fitting obtained with the C¹⁷O (Fig. 5 in Cabedo et al. 2021), where the redshifted component is almost always present in the C¹⁷O spectra.

Peak velocity value (top panels of Figs. B.4, B.5 and B.6) ranges are similar for all molecules (from ~8 to 9 km s⁻¹), although the blueshifted component goes to slightly lower velocities for H¹³CO⁺ (J=3-2), being the lowest velocities ~7.8 km s⁻¹. This decrease occurs toward the southeastern region of the core, at ~5″ from the center. We attribute this to the presence of the third velocity component produced by outflow contamination observed in Fig. B.3, which produces a bad model of the line profile. No gradient indicating rotation is observed in the velocity maps for any tracer.

Velocity dispersion (bottom panels of Figs. B.4, B.5 and B.6) ranges are also similar for all tracers (from 0.1 to 0.5 km s⁻¹). However, there are some main differences: velocity dispersion is generally larger for H¹³CO⁺ (J=3-2) with average values of ~0.4 km s⁻¹, whereas for DCO⁺ (J=3-2) and H¹³CO⁺ (J=1-0) mean values are around ~ 0.25 km s⁻¹. Very large values of the velocity dispersion (>0.5 km s⁻¹) are observed in the same southeastern region of the blueshifted component of H¹³CO⁺ (J=3-2) which is a consequence of the poor fitting due to the third velocity component. Values for the velocity dispersion are generally larger for the blueshifted than for the redshifted velocity component, indicating the two components trace gas with different kinematics, i.e., different infall velocities are affected differently by turbulence. For all tracers and both velocity components, the velocity dispersion increases toward the center, as the temperature rises. This increase is smaller for the H¹³CO⁺ (J=1-0) emission (up to ~0.35 km s⁻¹) than for the DCO+ (J=3-2), and H¹³CO⁺ (J=3-2) (up to ~0.5 km s⁻¹) indicating that H¹³CO⁺ (J=1-0) traces gas at colder temperatures and further to the center of the source.

In general, there is good agreement between the values obtained for DCO⁺ (J=3-2) and H¹³CO⁺ (J=3-2), both for peak velocity and velocity components. Considering that we are tracing similar scales, this confirms that both tracers are tracing similar gas. However, C¹⁷O presents some differences with respect to H¹³CO⁺ (J=1-0) which indicates that they are tracing slightly different gas, C¹⁷O tracing warmer and denser gas closer to the protostar.

Fig. B.1 Spectral map of the DCO+ (J = 3–2) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel corresponds to 0.5″(~ 82 au). The red spectrum refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1).

Fig. B.2 Spectral map of the H13CO+ (J=1–0) emission in the inner 658 au centered on the dust continuum emission peak. The whole map is 4″ × 4″ and each pixel corresponds to 0.25″ (~ 41 au). The red spectrum refers to the peak of the continuum emission, and the blue line indicates the systemic velocity (8.3 km s−1).

Fig. B.3 Spectral map of the H13CO+ (J = 3–2) emission in the inner 900 au centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel corresponds to 0.5″(~ 82 au). The red spectrum refers to the peak of the continuum emission and the blue line indicates the systemic velocity (8.3 km s−1)

Fig. B.4 DCO+ (J=3–2) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, 10, and 20 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.

Fig. B.5 H13CO+ (J=1–0) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, and 10 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.

Fig. B.6 H13CO+ (J=3–2) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, 10, and 20 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.

Fig. C.1 Expected optical depth maps for the four molecular tracers (see text in Section C for details). Top left: C17O (J=1–0). Top right: DCO+ (J=3–2). Bottom left: H13CO+ (J=1–0). Bottom right: H13CO+ (J=3–2).

Appendix C Estimation of the opacity effects

As the line profile modeling does not allow us to obtain the line opacities, we estimated the optical depth using the dust continuum emission and the molecular abundances in order to determine the impact on the parameters derived from the line modeling. We first compute the H₂ column density map from the dust continuum emission map at 110 GHz (shown in Fig. 1), assuming its emission is optically thin (e.g., tracing the full column density) and applying the following standard abundances: [C¹⁷O] = 5×10⁻⁸ (Thomas & Fuller 2008), [DCO⁺] = 10⁻¹¹ and [H¹³CO⁺] = 10⁻¹⁰ (Jørgensen et al. 2004). Figure C.1 shows the expected optical depth maps for the four different molecular lines we obtain with this method. Values of the optical depth are generally low (< 1), although they increase to larger values for all tracers in the central beam. We note that, as we are considering a constant temperature for the derivation of the H₂ column density, which is likely higher at the center, the H₂ column density values in the central compact source are probably smaller, and we are overestimating the optical depth in that region. It can be seen that both transitions of H¹³CO⁺ present larger values of optical depth than both C¹⁷O (J=1-0) and DCO⁺ (J=3-2). We acknowledge that these large values can potentially affect our deuteration measurements toward the central beam, but stress that the small optical depth values found in DCO⁺ confirm that the two velocity components are not mainly produced by self-absorption effects. Outside of that central beam, the expected line opacities are generally smaller than 1, so we found no need for applying an opacity correction to our data.

References

All Tables

All Figures

Fig. 1 Map of the dust continuum emission at 110 GHz clipped for values smaller than 3σ (see Table 1) and integrated emission for the four molecular lines as indicated in the upper left corner (red contours). Contours are −2, 3, 5, 8, 11, 14, and 17σ (the individual σ values are indicated in Table 1). The velocity ranges used to obtain the integrated intensity maps are shown in Table 1. In the bottom-left of each plot, we show the beam size for each molecular line (in red) and for the dust continuum emission (in black). The physical scale is shown in the top-right corner of the figure.

In the text

	Fig. 2 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the deuteration fraction maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The spatial scale is shown in the top-right corner of the upper panels.
In the text

	Fig. 3 Map of the CO depletion factor and contours of dust continuum emission at 110 GHz, for emission over 3σ. The spatial scale is shown in the top-right corner of the figure.
In the text

	Fig. 4 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the ionization fraction maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The spatial scale is shown in the top-right corner of the upper panels.
In the text

	Fig. 5 Dust continuum emission at 110 GHz for −2, 3, 5, 7, 10, 30, and 50σ (black contours) superimposed on the ionization rate maps (top panels) and on the corresponding statistical uncertainties (lower panels) for the blue- and redshifted velocity components (left and right columns, respectively). The scale is shown in the top-right corner of the upper panels.
In the text

Fig. 6 Upper left panel: ζCR as a function of the radius, r. The two dotted black lines at 270 and 500 au identify the radius ranges used to compute the slopes (s1 and s2, respectively). The gray shaded region shows the expected X-ray ionization in the outflow cavity for a typical and a flaring T Tauri star (Rab et al. 2017). Lower left panel: hydrogen column density, $N_{{\text{H}}_2}$, from the source center outwards as a function of radius. Upper right panel: ζCR as a function of the H2 column density. The upper x-axis shows the minimum energy that CR protons must have in order not to be thermalized (Padovani et al. 2018). The orange, magenta, and green filled regions in the three upper panels represent the envelope of the ζCR and $N_{{\text{H}}_2}$ profiles (see Sect. 5.2), while the red, purple, and dark green solid lines show their median value. Lower right panel: ζCR map of the blueshifted component (zoom of the upper left panel of Fig. 5) in logarithmic scale (colored map). The purple hatched region shows a circle of radius 270 au, where CRs propagate according to the diffusive regime. For illustration purposes, the black arrow at the position angle # shows a direction used to extract the ζCR profile from the map. The black dashed lines identify the radius as well as the corresponding values of $N_{{\text{H}}_2}$ and ζCR, where the slope of ζCR changes from s1 to s2.

In the text

	Fig. 7 Deuteration fraction maps for the blue- and redshifted velocity components (left and right columns, respectively). The N2D+ (J = 3–2) integrated emission is shown overlaid as black contours at 3, 5, 10, 15, and 20σ, where σ = 5.23 mJy beam−1. The spatial scale is shown in the top-right corners of each panel.
In the text

	Fig. B.1 Spectral map of the DCO+ (J = 3–2) emission in the inner 900 au, centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel corresponds to 0.5″(~ 82 au). The red 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. B.2 Spectral map of the H13CO+ (J=1–0) emission in the inner 658 au centered on the dust continuum emission peak. The whole map is 4″ × 4″ and each pixel corresponds to 0.25″ (~ 41 au). The red 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. B.3 Spectral map of the H13CO+ (J = 3–2) emission in the inner 900 au centered on the dust continuum emission peak. The whole map is 5.5″ × 5.5″ and each pixel corresponds to 0.5″(~ 82 au). The red 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. B.4 DCO+ (J=3–2) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, 10, and 20 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.
In the text

	Fig. B.5 H13CO+ (J=1–0) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, and 10 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.
In the text

	Fig. B.6 H13CO+ (J=3–2) maps obtained from modeling the line profiles with two velocity components. Overlaid contours show the integrated intensity at −3, 3, 5, 10, and 20 σ. The top map shows values for the peak velocity and bottom show values for the line width (given as σ). Left and right show the blue- and redshifted component, respectively.
In the text

	Fig. C.1 Expected optical depth maps for the four molecular tracers (see text in Section C for details). Top left: C17O (J=1–0). Top right: DCO+ (J=3–2). Bottom left: H13CO+ (J=1–0). Bottom right: H13CO+ (J=3–2).
In the text