EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 686 (June 2024) A&A, 686 (2024) A162 Full HTML
Open Access
Issue
A&A
Volume 686, June 2024
Article Number A162
Number of page(s) 13
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202347997
Published online 07 June 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open Access funding provided by Max Planck Society.

1 Introduction

Stars form in cold, dense cores within molecular clouds (Pineda et al. 2023). Probing these regions enables observations of the initial conditions for star and disk formation. Some of the key aspects still poorly constrained in such processes are the electron fraction, X(e), and cosmic-ray ionization rate, ζ(H2). Both of these quantities play an important role in astrochemical models (Ceccarelli et al. 2023) and in the coupling between magnetic fields and dense gas (Pineda et al. 2021; Pattle et al. 2023; Tsukamoto et al. 2023).

In the case of a constant cosmic-ray ionization rate, the electron fraction should follow a relation with a density of X(e) ∝ n−0.5 (McKee 1989; Caselli et al. 1998). The normalization of the relation depends substantially on whether the gas shows no depletion of metals (McKee 1989), appropriate for a low-density environment, or if the depletion of metals is taken into account in the modeling (Bergin & Langer 1997; Caselli et al. 1998), appropriate for denser regions.

Observationally, it is not possible to directly measure the electron fraction. Instead, it must be inferred from the combined analysis of several molecules. Toward dense cores, different works have derived typical values in the dense regions of X(e) ≈ 10−8−0−6 (Guelin et al. 1977; Caselli et al. 1998; Bergin et al. 1999; Maret & Bergin 2007, while in the transition from the diffuse to the dense medium in TMC1 (n(H2) ~ 103 cm−3), an electron fraction of X(e) ≈ 9.8 × 10−8−3.6 × 10−7 has been found (Fuente et al. 2019; Rodríguez-Baras et al. 2021). Recently, high angular resolution ALMA observations of the protostellar core B335 have provided an estimate of the electron fraction, X(e) ≈ 1−8 × 10−6, within 1000 au (Cabedo et al. 2023). This last work showed a radial variation of the electron fraction, with increased values toward the protostar.

Similarly, the cosmic-ray ionization rate is estimated toward low column density line of sights by observing the H3+ absorption feature toward background sources, which provide a value of ζ(H2) ≈ 2 × 10−16 s−1 in diffuse molecular gas (Indriolo & McCall 2012; Neufeld & Wolfire 2017; Padovani et al. 2020). In comparison, the ionization rate calculated using the particle spectra measured by the Voyager spacecrafts is 10 times lower ζ(H2) ≈ 10−17 s−1 (Cummings et al. 2016) at column densities between 1021 −1023 cm−2 (see Padovani et al. 2022), suggesting some degree of variation throughout the Galaxy. Toward higher column densities, other methods are required, and therefore, different molecular ratios and/or detailed modeling are used to derive these values. Theoretical models propose that the local value of the cosmic-ray ionization rate can be increased by the accretion process onto the protostar (Padovani et al. 2016; Gaches & Offner 2018) over the “standard” cosmic-ray ionization rate often adopted for molecular clouds, ζ(H2) ≈ 10−17 s−1 (Padovani et al. 2009). Observations of the OMC-2 FIR4 region have revealed enhanced values of the cosmic-ray ionization rate (Ceccarelli et al. 2014; Fontani et al. 2017; Favre et al. 2018). Similarly, ALMA observations of the protostellar core B335 have also suggested an enhanced value of the cosmic-ray ionization rate (Cabedo et al. 2023). These works show that a better understanding of the role of protostars and their surroundings is important to properly understand the star- and disk-formation processes.

Isolated cores, such as L1544 and B68, provide a great opportunity to study in detail the physical conditions of star and disk formation (Maret & Bergin 2007; Redaelli et al. 2021). However, more clustered environments provide a chance to study the mode of star formation where most stars form, although at the price of a more difficult environment to model. The NGC 1333 region in the Perseus cloud, at a distance of ≈300 pc (Zucker et al. 2018), offers a good opportunity to study an active star-forming region in a clustered environment (Kirk et al. 2006; Winston et al. 2010; Gutermuth et al. 2008). Notably, studies on the dense gas in the region, as traced by NH3 or N2H+, have revealed narrow filamentary structures with a coherent velocity (Friesen et al. 2017; Hacar et al. 2017; Dhabal et al. 2019; Chen et al. 2020).

In this work, we present the first results of the Probing the physics of star formation (ProPStar) survey, which attempts to study the physical conditions of star-forming regions with interferometric observations in order to connect the molecular cloud and disk scales. We present new observations of the NGC 1333 region, which allowed us to present the first resolved large-area maps of the electron fraction and cosmic-ray ionization rate of a low-mass star-forming region. We compare the electron fraction and cosmic-ray ionization rate maps with theoretical expectations and explore possible explanations for the observed spatial variations. Finally, we discuss the implications of the results.

2 Data

In the following, we summarize the single dish and interferometric observations for the H13CO+ and DCO+ lines. The whole NGC 1333 star-forming region is shown in Fig. 1, and the coverage area of the combined observations is marked by a dashed box, which includes the SVS 13 and NGC 1333 IRAS 4 systems.

thumbnail Fig. 1

Color image of the NGC 1333 region using Spitɀer IRAC 3.6, 4.5, and 8.0 µm. The image highlights the presence of outflows in a green color. The area covered by the combined NOEMA and 30 m observations is shown with the dashed box. The scale bar is shown in the bottom-right corner.
Table 1

Lines observed in the high-spectral resolution band.

2.1 IRAM 30 m telescope

The observations were carried out with the IRAM 30 m telescope at Pico Veleta (Spain) on 2021 November 9, 10, and 11 and on 2022 February 19 and 20 under project 091-21. The EMIR E090 receiver and the FTS50 backend were employed. We used two spectral setups to cover the H13CO+ (1−0) and DCO+ (1−0) lines at 72.0 and 86.7 GHz (see Table 1). We mapped a region of ≈150″ × 150″ with the On-the-Fly mapping technique and using position switching. The data reduction was performed using the CLASS program of the GILDAS package1. The beam efficiency, Beff, was obtained using the Ruze formula (available in CLASS), and it was used to convert the observations into main beam temperatures, Tmb. The noise for the H13CO+ (1−0) and DCO+ (1−0) cubes is 48 and 79 mK, in Tmb scale, respectively.

thumbnail Fig. 2

Maps of the DCO+ and H13CO+ (1−0) emission obtained with NOEMA and 30 m telescopes. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the first level is at log10(N(H2)) = 22.2, and the step between levels is 0.5 dex. The stars mark the positions of the YSOs identified by Dunham et al. (2015) using Spitɀer observations.

2.2 NOEMA interferometer

The observations carried out with the IRAM NOrthern Extended Millimeter Array (NOEMA) interferometer within the S21AD program using the Band 1 receiver were obtained on 2021 July 18, 19, and 21; August 10, 14, 15, 19, 22, and 29; and September 1 in the D configuration. We observed a total of 96 pointings, which were separated into four different scheduling blocks. The mosaic’s center is located at αJ2000 = 03h29m10.2s, δj2000 = 31º13′49.4″. We used the PolyFix correlator with aLO frequency of 82.505 GHz and an instantaneous bandwidth of 31 GHz spread over two sidebands (upper and lower) and two polarizations. The centers of the two 7.744 GHz-wide sidebands were separated by 15.488 GHz. Each sideband is composed of two adjacent basebands of ~3.9 GHz width (inner and outer basebands). In total, there are thus eight basebands that were fed into the correlator. The spectral resolution is 2 MHz throughout the 15.488 GHz effective bandwidth per polarization. Additionally, a total of 112 high-resolution chunks were placed, each with a width of 64 MHz and a fixed spectral resolution of 62.5 kHz. Both polarizations (H and P) are covered with the same spectral setup, and therefore the high-resolution chunks provide 66 dual polarization spectral windows. The high spectral resolution windows used in this work are listed in Table 1.

2.3 Image combination

We resampled the original 30 m data to match the spectral setup of the NOEMA observations. We used the task uvshort to generate the pseudo-visibilities from the 30 m data for each NOEMA pointing. The imaging was done with natural weighting, a support mask, and using the SDI deconvolution algorithm.

The H13CO+ (1−0) was then convolved to match the DCO+ (1−0) beam size using the spectral-cube and RadioBeam Python packages. The noise level of the combined images is reported in Table 1. Both cubes were converted to units of K, and the integrated intensity maps were calculated between 5.4 and 10 km s−1, which covers all the emission seen in both molecules. The integrated intensity maps are shown in Fig. 2.

2.4 H2 column density

We used the total H2 column density, N(H2), derived using the Herschel observations (Pezzuto et al. 2021), which are available in the Herschel Gould Belt Survey repository2. The effective angular resolution of the N(H2) maps is 18.2″. The map was regridded to match the DCO+ integrated intensity map. The resulting map is shown in Fig. 3.

2.5 JCMT C18O observations

We used the C18O (3−2) observations of NGC 1333 taken with HARPS at JCMT (Curtis et al. 2010). The angular resolution of these observations is 17.7−, and the main beam efficiency is 0.66 (as used in Curtis et al. 2010). We smoothed the C18O data to match the angular resolution of the Herschel-based N(H2). The integrated intensity map we calculated is between 6 and 10 km s−1 and shown in Fig. 4. This range covers all the emission seen in the cube.

3 Analysis

3.1 Estimate of ionization fraction and cosmic-ray ionization rate

One of the most commonly used ionization fraction tracers is [DCO+]/[HCO+] (Guelin et al. 1977, 1982; Dalgarno & Lepp 1984; Caselli et al. 1998). We estimated the ionization fraction, X(e), and the cosmic-ray ionization rate, ζ(H2), following Caselli et al. (1998). The method uses the main reaction paths for the formation and destruction for HCO+ and DCO+. Although other line ratios and techniques have been explored (Bron et al. 2021), we used a method that combines the bright lines available in the observations and the CO depletion. The analysis relations between the observables and the desired physical parameters are X(e)=2.7×10−8RD−1.2×10−6fD,${X(e) = {{2.7 \times {{10}^{ - 8}}} \over {{R_{\rm{D}}}}} - {{1.2 \times {{10}^{ - 6}}} \over {{f_D}}},}$(1) ζ(H2)=[ 7.5×10−4X(e)+4.6×10−10fD ]X(e)n(H2)RH,${\zeta \left( {{{\rm{H}}_2}} \right) = \left[ {7.5 \times {{10}^{ - 4}}X(e) + {{4.6 \times {{10}^{ - 10}}} \over {{f_D}}}} \right]X(e)n\left( {{{\rm{H}}_2}} \right){R_{\rm{H}}},}$(2)

where RD ≡ [DCO+]/[HCO+]; fD ≡ [12CO]/[H2]/ ([12CO]/[H2])flduclal; ([12CO]/[H2])flduclal is the expected 12CO abundance; Rh ≡ [HCO+]/[12CO]; and n(H2) is the average H2 number density.

However, since HCO+ is usually optically thick and it also traces outflow emission (not traced by DCO+), we used H13CO+ and the canonical isotropic ratio of [12C]/[13C] = 68 (Milam et al. 2005) to derive RD. Similarly, we used C18O and the canonical isotropic ratio of [12CO]/[C18O] = 560 (Wilson & Rood 1994) to estimate the 12CO column density.

thumbnail Fig. 3

Map of the Herschel-based H2 column map, N(H2). The contours are drawn starting at log10(N(H2)) = 22.2 and with a step between levels of 0.5 dex. These contours are the same as in Fig. 2. The beam size and scale bar are shown in the top left and bottom right, respectively.
thumbnail Fig. 4

Integrated intensity map of the C18O (3−2) transition line. The contours correspond to the Herschel-based H2 column map, with the same contours as in Fig. 2. The beam size and scale bar are shown in top left and bottom right, respectively.

3.2 Column densities

3.2.1 DCO+ and H13CO+

We used the optically thin approximation to derive the column densities (Mangum & Shirley 2015). In this case, we used the following expression to calculate the total column densities: Ntot=8π  ν3c3Q(Tex)Aulgue−Eup/hTex1(ehv/kBTex−1)∫​TmbdυJ(Tex)−J(Tbg),${N_{{\rm{tot}}}} = {{8\pi \,\,{\nu ^3}} \over {{c^3}}}{{Q\left( {{T_{{\rm{ex}}}}} \right)} \over {{A_{{\rm{ul}}}}{g_{\rm{u}}}{e^{ - {E_{{\rm{up}}}}/h{T_{{\rm{ex}}}}}}}}{1 \over {\left( {{e^{hv/{k_{\rm{B}}}{T_{{\rm{ex}}}}}} - 1} \right)}}{{\mathop \smallint \nolimits^ {T_{{\rm{mb}}}}{\rm{d}}\upsilon } \over {J\left( {{T_{{\rm{ex}}}}} \right) - J\left( {{T_{{\rm{bg}}}}} \right)}},$(3)

where v is the frequency of the transition observed, Aul is the Einstein coefficient for the transition from level u to l, 𝑔u is the degeneracy of level u, Eup is the energy of level u, Tex is the excitation temperature, Q(T) is the partition function at temperature T, and J(T) ≡ hv/kB (exp(hv/kB T) − 1) is the Rayleigh-Jeans equivalent temperature. We used the Aul, 𝑔u, and Eup listed in the LAMDA database, and implemented in the molecular_columns package3.

In this work, we used a constant excitation temperature of 10 K for both transitions. Moreover, we only used estimates of the column density toward pixels with a signal-to-noise ratio of the line above 7.5 (Tpeak > 7.5 × rms) in order to obtain robust column densities.

3.2.2 C18O

The C18O column density was derived using the optically thin approximation, N(C18O)=5×1012TexK exp (31.6 K/Tex)∫​Tmbdυkm s−1cm−2,$N\left( {{{\rm{C}}^{18}}{\rm{O}}} \right) = 5 \times {10^{12}}{{{T_{{\rm{ex}}}}} \over {\rm{K}}}\exp \left( {31.6{\rm{K}}/{T_{{\rm{ex}}}}} \right){{\mathop \smallint \nolimits^ {T_{{\rm{mb}}}}{\rm{d}}\upsilon } \over {{\rm{km}}\,\,{{\rm{s}}^{ - 1}}}}{\rm{c}}{{\rm{m}}^{ - 2}},$(4)

and an excitation temperature of Tex =12 K, which were previously used by Curtis et al. (2010). These excitation temperature values are similar to those derived using 12CO (1−0) but with lower angular resolution in this region (Pineda et al. 2008). We also estimated the optical depth in the map for the assumed Tex and obtained a typical value of ~0.4, which suggests that the optical depth is not an issue across the map. However, it might produce lower limits to the derived N(C18O). Appendix A shows details regarding the optical depth calculation and a figure showing the map.

thumbnail Fig. 5

Map of the CO depletion fraction, fD. The beam size and scale bar are shown in the top-left and bottom-right corners, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2.

3.3 CO depletion factor

We estimated the CO depletion factor, fD, using the C18O column density, N(C18O), the isotopologue abundance ratio of 12CO/C18O = 560 (Wilson & Rood 1994), the H2 column density derived from Herschel, and the expected [12CO]/[H2] = 2.7 × 10−4 value, as determined in Taurus (Punanova et al. 2022) and comparable to the value of 2.5 × 10−4 obtained in NGC1333 at a lower angular resolution (Pineda et al. 2008). When comparing the N(H2) and N(C18O) column densities maps, it was clear that the H2 map traces more extended emission, and therefore, an offset is present, N(H2)=N(C18O)[ H2 ][ C18O ]+δ18.$N\left( {{{\rm{H}}_2}} \right) = N\left( {{{\rm{C}}^{18}}{\rm{O}}} \right){{\left[ {{{\rm{H}}_2}} \right]} \over {\left[ {{{\rm{C}}^{18}}{\rm{O}}} \right]}} + {\delta _{18}}.$(5)

This offset represents the H2 column in the molecular cloud that is not traced by the C18O (3−2) emission, and therefore it needed to be removed from the abundance calculations. We estimated the offset, δ18 = 2.7 × 1021 cm−2, as the median value of the difference observed assuming an undepleted abundance in the column density range of 3 × 1021 cm−2 < N(H2) < 5 × 1021 cm−2. This is equivalent to determining the offset required to achieve undepleted abundances at the lower column densities in the region. The resulting map is shown in Fig. 5, with depletion fractions increasing at higher column densities.

3.4 Volume density

We estimated the volume density as n(H2)=N(H2)−δ18δL,$n\left( {{{\rm{H}}_2}} \right) = {{N\left( {{{\rm{H}}_2}} \right) - {\delta _{18}}} \over {\delta L}},$(6)

where the N(H2) is the H2 column density, δ1g is the derived offset with C18O , and δL is the depth of the structure. We assumed δL = 0.4 pc (~60 × 103 au), which is the depth for HCO+ derived for the nearby L1451 (Storm et al. 2016). The corresponding volume density map is shown in Fig. 6, and it displays densities similar to those constrained from 12CO and 13CO (1−0) observations (Pineda et al. 2008). The typical density in places with a detection of DCO+ and H13CO+ is ≈4 × 103 cm−3, which suggests that this density is likely a lower limit since this region is bright in NH3 (Friesen et al. 2017; Dhabal et al. 2019) and N2H+ (Hacar et al. 2017), which are reliable high-density tracers.

thumbnail Fig. 6

Volume density map estimated from Herschel H2 column density. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2.

3.5 Electron fraction

We used the relation of Eq. (1) together with the previously presented measurements to estimate the electron fraction across the active star-forming region, and it is presented in the left panel of Fig. 7. The underlying distribution of the electron fraction in the region was estimated using the kernel density estimate (KDE), and it is shown in the right panel of Fig. 7. A typical value (median of the log10) of 〈X(e)〉 = 10−6.5 was derived from these measurements. We estimate that the error associated with using the wrong excitation temperature (in the range of 5 to 30 K) in the column density determination yields an uncertainty of less than 10% in RD and X(e).

3.6 Cosmic-ray ionization rate

We used the relation of Eq. (2) together with the previously presented measurements to estimate the cosmic-ray ionization rate across the active star-forming region, and the estimate is presented in the left panel of Fig. 8. The underlying distribution of the cosmic-ray ionization rate in the region was estimated using the KDE, and it is shown in the right panel of Fig. 8. We derived a typical value (median of the log10) of 〈ζ(H2)〉 = 10−16.3 s−1 from these measurements. We estimate that the error associated with using the wrong excitation temperature (with the temperature within the range of 5 and 30 K) RH could be underestimated by up to a factor of two, in addition to the uncertainty of less than 10% in RD and X(e) mentioned in Sect. 3.5. This indicates that the ζ(H2) could be underestimated by a factor of two toward hot regions (Tex > 20K), but elsewhere, these uncertainties are of the order 10–30%.

thumbnail Fig. 7

Derived electron fraction in the region. Left: map of the electron fraction, X(e). The map shows a smooth distribution but with a clear enhancement in the northwest section. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. Right: Kernel density estimate of the X(e) distribution across the mapped region. A typical value (median of the log10) of 10−6.5 was derived from these measurements and is marked with a red dashed vertical line.
thumbnail Fig. 8

Derived cosmic-ray ionization rate in the region. Left: map of the cosmic-ray ionization rate, ζ(H2). The map shows a smooth distribution but with the value clearly increasing in the northwest section. The beam size and scale bar are shown in top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. Right: Kernel density estimate of the ζ(H2) distribution across the mapped region. The median value of the log10 ζ(H2)/s−1 = −16.5 and the canonical value of 10−17 s−1 are marked with the blue and black dashed vertical lines, respectively.

4 Discussion

4.1 Electron fraction and cosmic-ray ionization rate distributions

The spatial distribution for both of the main parameters estimated in this work, X(e) and ζ(H2), are shown in the left panel of Figs. 7 and 8. These maps show the first resolved maps of these quantities, which present little substructure, except close to the cloud or structure’s edge and close to the young protostars. These maps were derived at the angular resolution of the DCO+ map, but we note that since n(H2) and ƒD. were derived using Herschel data, we assumed that they smoothly vary between 16″ and ≈6″, which might be less accurate toward the young protostars.

4.2 Electron fraction as a function of density

Several attempts have been carried out to parametrize the variation of the electron fraction as a function of density. McKee (1989) presented an analytic relation for the electron fraction as a function of density assuming no depletion of metals and a constant cosmic-ray ionization rate. This relation is expected to work even at high densities (McKee & Ostriker 2007). After the depletion of metals was identified as an important process in the denser regions within star-forming regions, an updated electron fraction relation with the density relation was presented by Caselli et al. (2002b), where a constant cosmic-ray ionization rate was assumed. Further detailed modeling of the prestellar core L1544 using different molecular transitions constrained the radial profile electron fraction (Redaelli et al. 2021). This derived relation is closer to that derived by Caselli et al. (1998), while (McKee 1989) provides a clear overestimation.

We compared the KDE distribution of the values determined from the data and the analytic relations in Fig. 9. In addition, we added a comparison with the relation obtained from modeling the L1544 prestellar core (Redaelli et al. 2021). A decrease in the electron fraction, X(e), at higher densities, n(H2), is a common trend from the analytical (and modeling) results. Notably, Fig. 9 shows a consistent increase in the electron fraction toward higher density areas, while beyond n(H2) ≳ 103.7 cm−3, the relation steepens. This steeper relation is mostly driven by the northwest area in the map, which is coincident with a previously suggested bubble interacting with the cloud (e.g., Dhabal et al. 2019; De Simone et al. 2022).

A key difference between the analytical relations for electron fraction as a function of density and data is the assumption of no local sources of cosmic rays. Recently, different works have suggested that enhanced levels of cosmic rays (thanks to their local generation by YSOs, e.g., Ceccarelli et al. 2014; Padovani et al. 2016) could provide a higher value of electrons than previously expected. This important difference might be the main physical reason as to why the analytical relations have such a poor match to the data, even in the case of McKee (1989), who reported a similar predicted value at the typical density (≈ 103.7 cm−3) even though there is a clear measurement of the depletion of metals (ƒD).

The derived densities are comparable to previous estimates. The mean density in the region probed with 12CO and 13CO (1−0) transitions are ≈103 cm−3 (Pineda et al. 2008). Therefore, although our density estimate is rather simple and likely a lower limit, it should not be the main reason for the disagreement between the measured electron fractions and the theoretical models.

4.3 Local generation of cosmic rays

Recently, several indirect observations have suggested an extremely high ionization rate over 10−15 s−1 in protostellar environments (e.g., Ceccarelli et al. 2014; Fontani et al. 2017; Favre et al. 2018; Cabedo et al. 2023). Also, synchrotron emission, the fingerprint of the presence of relativistic electrons, has been detected in the shocks that develop along protostellar jets (e.g., Beltrán et al. 2016; Rodríguez-Kamenetzky et al. 2017; Osorio et al. 2017; Sanna et al. 2019). Neither signature can be explained by interstellar cosmic rays since they are strongly attenuated at the high column densities typical of protostellar environments. As argued by Padovani et al. (2016) and Padovani et al. (2021), shocks forming along the jets and on the protostellar surface could be the regions where the locally accelerated cosmic rays are produced.

The left panel of Fig. 8 shows a dramatic local increase of ζ(H2) occurring in the northwest (upper) section of the map. Furthermore, there is a noticeable enhancement of ζ(H2) when in proximity to the embedded YSOs, which are also in the central area of the studied region. In fact, the volume density in dense cores around YSOs (typically within <0.1 pc, Kirk et al. 2006) is much higher than the “large-scale” estimate based on Eq. (6). Therefore, according to Eq. (2), the actual value of ζ(H2) is expected to be further drastically enhanced around these YSOs. In order to gain insights into the underlying mechanisms of the observed phenomenon, in Fig. 10, we have overlaid the cosmicray ionization rate map with the outflow orientation of the different embedded sources (see Table D.1), the dust polarization vectors obtained from the JCMT SCUBA2-POL Bistro Survey (Ward-Thompson et al. 2017) on this region (Doi et al. 2020, 2021), and the X-ray sources detected with Chandra (Winston et al. 2010).

First, we ruled out the X-ray sources as a possible explanation of the observed ionization enhancement. A significant attenuation of X-ray photons produced by source emitting in the keV energy range starts at column densities of ~1023 cm−2 (Igea & Glassgold 1999), which is about the measured maximum value of N(H2) in the upper section (where most of the sources are located). Therefore, X-rays would result in a relatively isotropically enhanced ζ(H2) around the sources, which is clearly incompatible with the distribution of ζ(H2) seen in this region, revealing no correlation with the source locations.

However, we have strong indicators that the sources of the ionization enhancement are the young embedded objects. The first obvious fact is the noticeable increase of ζ(H2) around most of the YSOs. Given that we strongly underestimated the gas density in the proximity to these objects, the actual values of ζ(H2) are expected to be much higher than those shown in Fig. 8. Second, the elongated region with the highest local ionization rate of ≳10−15 s−1 in the upper section of the map is approximately parallel to the magnetic field lines, practically connecting the SVS13 system and Per-emb-15. If these objects are able to generate energetic charged particles, their further propagation and the resulting ionization must occur along the local field lines (Fitz Axen et al. 2021). We note that in this work, we also include a figure of the difference between the observed CRIR and the expected values (Padovani et al. 2018) in Appendix C, which shows the same trends.

We compared the CRIR as a function of the typical YSO distance, calculated as the harmonic mean of the individual distances4, 1dYSO2=1N∑i=1N1di2,${1 \over {d_{{\rm{YSO}}}^2}} = {1 \over N}\mathop {\mathop \sum \nolimits^ }\limits_N^{i = 1} {1 \over {d_i^2}},$

where di is the distance to the different YSOs marked in Fig. 10. Figure 11 shows the KDE for the map, and it shows that within 15 000 au, the CRIR is relatively constant, while there is a reduction in the CRIR at larger distances. This is evidence that the YSOs are related to the increase of the CRIR. Although, we cannot confirm a correlation between YSO luminosity and CRIR due to the small sample size and the issues determining accurate CRIR toward the YSOs.

It is worth noting that the three-color Spitɀer image in Fig. 1 shows a strong outflow emission (green) from the SVS13 region, which overlaps with a fraction of the high enhancement in the cosmic-ray ionization rate. The SVS 13 region hosts several strong and well-studied outflows (Maret et al. 2009; Hodapp & Chini 2014; Lefèvre et al. 2017; Dionatos et al. 2020). It is possible that part of the local enhancement in CRIR is related to the outflow cavity, especially if we are probing the emission closer to the outflow cavity itself. Models following CR propagation indicate that the CR distribution near accelerating protostellar sources is likely inhomogeneous with the CR enhancement occurring along the direction of the outflow (Fitz Axen et al. 2021).

We would like to emphasize that most of the observed region (in particular, the lower section of the map in Fig. 8) is dominated by ζ(H2) ~ 10−17 s−1. This level of ionization, likely generated by interstellar cosmic rays, is comparable to the “standard” value of the ionization rate in diffuse molecular clouds Glassgold & Langer (1974). This value is also in good agreement with constraints on the CRIR derived from energy arguments. For the Milky Way, Krumholz et al. (2023) inferred a mean primary ionization budget due to star-formation activity of ζ(H2) ~ 2–5 × 10−17 s−1, significantly below the value measured in diffuse molecular gas. These values can be reconciled if CR acceleration sources are distributed inhomogeneously.

Remarkably, we can independently estimate the expected upper bound of the ionization rate near NGC 1333 produced by interstellar cosmic rays in the surrounding diffuse gas. One of the measurements of H3+ ions reported by Indriolo & McCall (2012) was performed in the direction of the target star HD 21483. Along this line of sight, the measurements probe the cosmic-ray ionization at the same radial distance of 300 pc, which is only 4–5 pc away from NGC 1333 in the plane of sky. The resulting non-detection yields a 1σ upper limit of 9 × 10−17 s−1 for the ionization rate. Given that the obtained spectra of H3+ lines show no visible traces of absorption (see their Fig. 10), we conclude that the actual value of ζ(H2) in this region is substantially lower and may well be close to the expected “standard” value.

In a separate publication, we will present a detailed analysis of the ionization rate map. We will also perform a quantitative comparison of the results with available theories of propagation of interstellar and locally accelerated cosmic rays.

thumbnail Fig. 9

Electron fraction, X(e), as a function of H2 density. The KDE of the underlying distribution was obtained from the data. The contour levels are drawn starting at 0.5σ and progressing outward in steps of 0.5σ where the σ levels are equivalent to that of a bivariate normal distribution. The analytic relations from McKee (1989) and Caselli et al. (2002b) as well as the result from chemical modeling of L1544 (Redaelli et al. 2021) are shown with orange, blue, and green curves, respectively. The analytic relations shown here do not take into account the local generation of cosmic rays.
thumbnail Fig. 10

Outflows and magnetic field orientation over ζ(H2). The background image is the derived ζ(H2) map, as shown in Fig. 8. The region identified with a local cosmic-ray ionization rate enhancement is marked with an orange contour. The outflow lobes (as reported in previous works) are shown by red and blue arrows. The magnetic field orientation, measured from continuum polarization observations with SCUBA2-POL, are shown by the segments of white lines. The X-ray sources are marked with red crosses. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.
thumbnail Fig. 11

Cosmic-ray ionization rate as a function of distance to YSOs. The KDE of the CRIR and the typical distance to the YSOs, dYSO, was calculated over the map. The dotted and dashed lines show the inner constant section, and the two different power-law decays help guide the eye.

4.4 Evidence of shocks and interaction

The spatial distribution for both of the main parameters estimated in this work, X(e) and ζ(H2), are shown in the left panel of Figs. 7 and 8. Both quantities present a clear one-order-of-magnitude enhancement toward the northwest in the mapped region. This could be consistent with the suggestion that the region is affected by interaction with an expanding shell (De Simone et al. 2022) coming from the northwest section of the map (the SVS13 region). A similar interaction, although at a much larger scale, is suggested in the nearby Taurus cloud with the Local Bubble and the Per-Tau shell (Soler et al. 2023), which could also be occurring in the Perseus Cloud.

A second region along the southwest of the map has previously been suggested as showing evidence of interaction between the cloud and a turbulent cell (Dhabal et al. 2019) or of a second bubble impacting the cloud (De Simone et al. 2022). However, we do not see clear evidence in X(e) or ζ(H2) of any of these features, although we cannot rule them out.

4.5 The relevance of nonideal-MHD effects

On scales of a couple of beams away from the embedded sources, the ionization fraction is ≈10−6.5 (see Fig. 7). It is usually assumed that at a relatively high electron fraction, the nonideal-MHD terms are negligible. Studying the initial core collapse phase Wurster et al. (2018) showed that above cosmic-ray ionization rates of 〈ζ(H2)〉 > 10−14s−1, the results are practically indistinguishable from results obtained by assuming zero resistivity (i.e., ideal MHD). The values measured in NGC 1333 are slightly higher with respect to the often-assumed canonical value for prestellar cores of 〈ζ(H2)〉 ≈ 10−17s−1 (Spitzer & Tomasko 1968; Umebayashi & Nakano 1990), and therefore, the nonideal-MHD effects are expected to not play an important role in current and future disk formation.

4.6 MHD wave cutoff

In a molecular cloud where the ionization level is low, the Alfvé waves cannot propagate when the collision frequency between ions and neutrals is comparable to or smaller than the MHD wave frequency. Therefore, a critical length for wave propagation beyond which waves do not transmit can be defined. The critical length for wave propagation (Stahler & Palla 2005; Mouschovias et al. 2011) is obtained from λmin=π4μ  mH  n(H2)Bon(H2) X(e) 〈 σ  υ 〉,${\lambda _{\min }} = \sqrt {{\pi \over {4\mu \,\,{m_{\rm{H}}}\,\,n\left( {{{\rm{H}}_2}} \right)}}} {{{B_o}} \over {n\left( {{{\rm{H}}_2}} \right)X(e)\left\langle {\sigma \,\,\upsilon } \right\rangle }},$(7)

where Bo is the unperturbed magnetic field and 〈σ υ〉 is the rate coefficient for elastic collisions. The rate coefficient term is approximated by the Langevin term, 〈 σ  υ 〉=1.69×10−9cm3 s−1,$\left\langle {\sigma \,\,\upsilon } \right\rangle = 1.69 \times {10^{ - 9}}{\rm{c}}{{\rm{m}}^3}{{\rm{s}}^{ - 1}},$(8)

for HCO+−H2 collisions (McDaniel & Mason 1973). We estimated the magnetic field strength, Bo, using the relation of Crutcher et al. (2010), a typical volume density of 103.6 cm−3, and the median value of X(e) = 10−6.5. With these values, we then obtained an MHD wave cutoff scale of 0.054 pc (11 100 au). Therefore, we predicted that if the turbulence is MHD in nature, then a dissipation scale should be observed at ≈0.054 pc (or ≈37″ at the distance of Perseus).

Previously, this scale was proposed as a possible origin for the transition to coherence between the supersonic cloud and subsonic cores (~0.1 pc; Goodman et al. 1998; Caselli et al. 2002a; Pineda et al. 2010). Recent observations of dense gas with different tracers (e.g., N2H+ and NH3) and angular resolutions showed that the southern end filament presents subsonic levels of turbulence at scales of ≈40″ (Friesen et al. 2017; Hacar et al. 2017; Dhabal et al. 2019; Sokolov et al. 2020), which corresponds to 0.054 pc, or 11 200 au, at the distance of Perseus. This is consistent with the estimated dissipation scale; however, studies of turbulence at higher angular resolution than 30″ should reveal this scale.

4.7 Disk sizes

Of the 12 Class 0 and I protostars in this region, dust emission observations of disks have resolved only six sources with either VLA or ALMA continuum observations (Segura-Cox et al. 2016, 2018; Diaz-Rodriguez et al. 2022). The mean dust disk radii measured for these six resolved sources is ≈30 au (see Table D.1). These values are lower than the typical values obtained for Class 0 and I objects (median values of 61 au and 81 au for Class 0 and I, respectively; see Tsukamoto et al. 2023, for a recent review).

Theoretical models suggest disk properties are sensitive to the degree of coupling between ions and neutrals, with perfect coupling corresponding to the “ideal” MHD limit in which angular momentum is efficiently transported away from the circumstellar region by magnetic braking (Zhao et al. 2020). Higher levels of CR ionization effectively enhance magnetic braking, thereby leading to smaller disks (Dapp et al. 2012; Zhao et al. 2016; Wurster et al. 2016). For example, Kuffmeier et al. (2020) carried out a suite of simulations that followed the collapse of spherical cores and showed that the disk size decreases with an increasing cosmic-ray ionization rate for identical initial conditions because of more efficient magnetic braking. Although prestellar cores appear as more dynamic entities in larger scale simulations (Kuffmeier et al. 2017; Lebreuilly et al. 2021; Offner et al. 2022), the results from Kuffmeier et al. (2020) provide a well-controlled experiment with different mass-to-flux-ratios λmag≡McπRc2B02πG,${\lambda _{{\rm{mag}}}} \equiv {{{M_{\rm{c}}}} \over {\pi R_{\rm{c}}^2{B_0}}}2\pi \sqrt G ,$(9)

where B0 is the magnetic field strength, and Mc and Rc are the mass and radius of the core.

For the typical values of the ionization and cosmic rays found in this work, 〈X(e)〉 = 10−6.5 and 〈ζ(H2)〉 = 10−16.5 s−1 = 3 × 10−17 s−1, the typical gas disk sizes found in this region are well matched by the cases in Kuffmeier et al. (2020) with λmag = 5 and an initial ratio of rotational to gravitational energy, β, of the spherical core between β = 0.025 and β = 0.1 (see also Wurster et al. 2018). Similarly, Lebreuilly et al. (2021) investigated the disk size distribution in nonideal-MHD models, assuming an ionization rate of 〈ζ(H2)〉 = 10−17 s−1, of a collapsing clump of 1000 M⊙. They found an average disk radius of ≈50 au, which is about a factor of two larger than what is found in our observations. This might indicate that the relatively higher cosmic-ray ionization rate in NGC 1333 leads to smaller disk sizes, although a lower angular momentum or higher magnetization cannot be ruled out.

We highlight that the cosmic-ray ionization rate is not a constant value throughout the entire region, and it may fluctuate, especially toward higher densities closer to the actively accreting protostars. At these high densities, this rate could decrease due to the shielding of cosmic rays (Nakano & Tademaru 1972; Umebayashi & Nakano 1981; Padovani et al. 2009), or it could increase due to local enhancement of cosmic-ray ionization by the individual protostars themselves (Padovani et al. 2016; Gaches & Offner 2018).

4.8 Improvements on X(e) and ζ(H2)

Although the method for determining the electron fraction and cosmic-ray ionization rate is widely used, there are clear paths for improvements. The assumption of a single excitation temperature could impact the total column densities derived and the derived X(e) and ζ(H2). However, the H13CO+ and DCO+ (1−0) lines have similar parameters (Eup and Aul), and therefore, this assumption is not expected to cause major issues. Also, other transitions of these species, combined with an excitation analysis, could better constrain the optical depth and excitation temperatures of these species and therefore derive more reliable molecular column densities.

Recently, Sabatini et al. (2023) has shown that in high-mass star-forming regions, where DCO+ and H13CO+ might trace different volumes or be affected by high optical depths, H2D+ is a better estimator for the cosmic-ray ionization rate. Unfortunately, this approach is severely limited in low-mass star-forming regions since the H2D+ emission is compact and only traces the highest density regions in dense cores (e.g., Caselli et al. 2003, 2008; van der Tak et al. 2005; Harju et al. 2006; Vastel et al. 2006; Friesen et al. 2010, 2014; Koumpia et al. 2020).

A complementary approach might involve analyzing multiple molecular species and/or transitions and comparing the line emission with a suite of chemical models to directly constrain the typical physical parameters of the region (e.g., density, radiation field, ζ(H2), “chemical” age). This could be carried out with single-zone models (0D) or more complex ones using some physical model for the region. Regardless of the approach, a new and exciting possibility of deriving robust and high-quality resolved maps of X(e) and ζ(H2) is upon us.

5 Conclusions

We have presented the first results of the ProPStar survey. We studied the southeast section of the NGC 1333 region using combined observations of the 30 m and NOEMA telescopes to recover all scales. The main results are listed below:

  1. We mapped the DCO+ and H13CO+ (1−0) emission at ≈6″ resolution (≈1 800 au).

  2. We used the analytic relation from Caselli et al. (1998) relating the measured DCO+ and H13CO+ column density ratio in addition to archival C18O and H2 in order to derive the first large-area map for the electron fraction, X(e), and cosmic-ray ionization rate, ζ(H2), in a low-mass star-forming region.

  3. The electron fraction and cosmic-ray ionization rate distributions are peaked at about 10−6.5 and 10−16.5 s−1, respectively, showing local enhancements around the embedded YSOs. These values are higher than expected from the “standard” cosmic-ray ionization rate in molecular clouds.

  4. The northwest region in the map shows highly elevated values for both the electron fraction and cosmic-ray ionization rate, which might be a sign of the interaction between the molecular cloud and a local bubble, or it could be due to cosmic rays being locally generated by YSOs.

  5. We used the typical value of X(e) in the region to estimate the critical length for wave propagation of MHD turbulence of 0.054 pc (≈37″). If the turbulence in the molecular clouds is of MHD nature, then this dissipation scale could be measured with observations of enough sensitivity and resolution.

  6. The increased values of X(e) and ζ(H2) are sufficiently high such that when compared to numerical simulations, the newly formed disks should be consistent with ideal-MHD simulations (e.g., small). This also suggests that disk formation in clustered regions could strongly depend on the feedback from the first stars that are formed.

In summary, these observations show the great potential offered by combined NOEMA and IRAM 30 m observations to determine the electron fraction and cosmic-ray ionization rate in the NGC 1333 region. Both of these quantities are important to understanding the physical conditions and processes involved in star and disk formation. More maps of the cosmic-ray ionization rate and electron fraction, combined with complementary disk size measurements, should provide a much better understanding of the role of stellar feedback (e.g., locally generated cosmic rays) and the importance of nonideal MHD in the whole star- and disk-formation processes.

Acknowledgements

The authors kindly thank the anonymous referee for the comments that helped improve the manuscript. Part of this work was supported by the Max-Planck Society. D.M.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102405. We gratefully acknowledge the helpful discussions with Elena Redaelli, Cecilia Ceccarelli, Giovanni Sabatini, and Sylvie Cabrit, which improved the article. We thank Elena Redaelli for providing the electronic version of the electron fraction and density data points from modeling L1544. We thank Yasuo Doi for providing the electronic version of the polarization vectors. This research has made use of NASA’s Astrophysics Data System.

Appendix A Optical depth of C18O

We estimated the C18O (3−2) line optical depth as τ(C18O)=−ln(1.0−Tp(C18O)Jv(Tex)−Jv(Tcmb)),$\tau \left( {{{\rm{C}}^{18}}{\rm{O}}} \right) = - \ln \left( {1.0 - {{{T_p}\left( {{{\rm{C}}^{18}}{\rm{O}}} \right)} \over {{J_v}\left( {{T_{ex}}} \right) - {J_v}\left( {{T_{cmb}}} \right)}}} \right),$(A.1)

where Tp(C18O) is the peak brightness of the line, and Jv(T) is the brightness temperature of a black body with temperature T at frequency v. In the case of C18O (3−2), we assumed Tex = 12 K.

Figure A.1 shows the map of the optical depth. The figure clearly shows that the emission is mostly optically thin, which is also consistent with the findings from Curtis et al. (2010).

thumbnail Fig. A.1

Optical depth of C18O (3−2). The red line shows the τ(C18O) = 1 contour. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.

Appendix B The effect of excitation temperature

One of the most important simplifications in this work was the use of a single excitation temperature for the column density calculation of H13CO+ and DCO+, which we used to determine the RD and RH quantities. We explored the effect of both species having the same excitation temperature but that temperature being different from 10 K (5 K < Tex < 12 K), and we found that the variation in RD is less than 5%, while for RH, the variation is less than 10%. Therefore, an improved excitation analysis of both species would provide a substantial check to the calculation used in this work.

Appendix C ∆ log10 ζ

We calculated the difference between the derived CRIR and the expected values from the L model of CRIR as a function of column density from Padovani et al. (2018), ∆log10 ζ = log10 ζobs − log10 ζmodel. Figure C.1 shows the map of ∆log10 ζ, which is similar to Fig. 8, with lower values than expected in the southeast section, local increases toward YSOs, and a general excess in the rest of the map.

thumbnail Fig. C.1

Difference between observed CRIR and the model-ℒ to better illustrate the local effects, ∆ log10 ζ = log10 ζobs − log10 ζmodel The red and blue colors show the excess and deficiency of CRIR with respect to the expected value following Padovani et al. (2018). The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.

Appendix D Disk sizes and outflow directions of covered YSOs

We list the source name, disk radii, and outflow directions found in the literature in Table D.1. Except for SVS13 VLA4A and SVS13 VLA4B, the disk sizes are the result of modeling VLA 8.1 mm continuum observations, which have a FWHM beam size of ~16 au (Segura-Cox et al. 2016, 2018). For sources without a reported disk radii, the disks were detected but unresolved in Segura-Cox et al. (2018), and the disk radii at VLA wavelengths could only be estimated as being less than or equal to 8 au. For SVS13 VLA4A and SVS13 VLA4B, which form a close binary that shares dusty circumbinary material, the dust disk radii of the individual circumstellar disks are reported from Gaussian fits to ALMA 0.9 mm continuum observations, with a FWHM beam size of ~26 au (Diaz-Rodriguez et al. 2022).

Table D.1

Properties of Class 0 and I protostars within the field of view.

References

  1. Beltrán, M. T., Cesaroni, R., Moscadelli, L., et al. 2016, A&A, 593, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Bergin, E. A., & Langer, W. D. 1997, ApJ, 486, 316 [Google Scholar]
  3. Bergin, E. A., Plume, R., Williams, J. P., & Myers, P. C. 1999, ApJ, 512, 724 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bron, E., Roueff, E., Gerin, M., et al. 2021, A&A, 645, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Cabedo, V., Maury, A., Girart, J. M., et al. 2023, A&A, 669, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [Google Scholar]
  7. Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002a, ApJ, 572, 238 [Google Scholar]
  8. Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 565, 344 [Google Scholar]
  9. Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37 [CrossRef] [EDP Sciences] [Google Scholar]
  10. Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Ceccarelli, C., Dominik, C., López-Sepulcre, A., et al. 2014, ApJ, 790, L1 [Google Scholar]
  12. Ceccarelli, C., Codella, C., Balucani, N., et al. 2023, ASP Conf. Ser., 534, 379 [NASA ADS] [Google Scholar]
  13. Chen, M. C.-Y., Di Francesco, J., Rosolowsky, E., et al. 2020, ApJ, 891, 84 [NASA ADS] [CrossRef] [Google Scholar]
  14. Chuang, C.-Y., Aso, Y., Hirano, N., Hirano, S., & Machida, M. N. 2021, ApJ, 916, 82 [NASA ADS] [CrossRef] [Google Scholar]
  15. Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466 [NASA ADS] [CrossRef] [Google Scholar]
  16. Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18 [CrossRef] [Google Scholar]
  17. Curtis, E. I., Richer, J. S., & Buckle, J. V. 2010, MNRAS, 401, 455 [NASA ADS] [CrossRef] [Google Scholar]
  18. Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47 [NASA ADS] [CrossRef] [Google Scholar]
  19. Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. De Simone, M., Codella, C., Ceccarelli, C., et al. 2022, MNRAS, 512, 5214 [CrossRef] [Google Scholar]
  21. Dhabal, A., Mundy, L. G., Chen, C.-y., Teuben, P., & Storm, S. 2019, ApJ, 876, 108 [CrossRef] [Google Scholar]
  22. Diaz-Rodriguez, A. K., Anglada, G., Blázquez-Calero, G., et al. 2022, ApJ, 930, 91 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dionatos, O., Kristensen, L. E., Tafalla, M., Güdel, M., & Persson, M. 2020, A&A, 641, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Doi, Y., Hasegawa, T., Furuya, R. S., et al. 2020, ApJ, 899, 28 [Google Scholar]
  25. Doi, Y., Tomisaka, K., Hasegawa, T., et al. 2021, ApJ, 923, L9 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11 [Google Scholar]
  27. Favre, C., Ceccarelli, C., López-Sepulcre, A., et al. 2018, ApJ, 859, 136 [Google Scholar]
  28. Fitz Axen, M., Offner, S. S. S., Gaches, B. A. L., et al. 2021, ApJ, 915, 43 [NASA ADS] [CrossRef] [Google Scholar]
  29. Fontani, F., Ceccarelli, C., Favre, C., et al. 2017, A&A, 605, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Friesen, R. K., Di Francesco, J., Myers, P. C., et al. 2010, ApJ, 718, 666 [NASA ADS] [CrossRef] [Google Scholar]
  31. Friesen, R. K., Di Francesco, J., Bourke, T. L., et al. 2014, ApJ, 797, 27 [NASA ADS] [CrossRef] [Google Scholar]
  32. Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63 [NASA ADS] [CrossRef] [Google Scholar]
  33. Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Gaches, B. A. L. & Offner, S. S. R. 2018, ApJ, 861, 87 [NASA ADS] [CrossRef] [Google Scholar]
  35. Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73 [NASA ADS] [CrossRef] [Google Scholar]
  36. Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223 [Google Scholar]
  37. Guelin, M., Langer, W. D., Snell, R. L., & Wootten, H. A. 1977, ApJ, 217, L165 [NASA ADS] [CrossRef] [Google Scholar]
  38. Guelin, M., Langer, W. D., & Wilson, R. W. 1982, A&A, 107, 107 [NASA ADS] [Google Scholar]
  39. Gutermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ, 674, 336 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hacar, A., Tafalla, M., & Alves, J. 2017, A&A, 606, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Harju, J., Haikala, L. K., Lehtinen, K., et al. 2006, A&A, 454, L55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Hodapp, K. W., & Chini, R. 2014, ApJ, 794, 169 [Google Scholar]
  43. Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13 [NASA ADS] [CrossRef] [Google Scholar]
  44. Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 [Google Scholar]
  45. Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
  46. Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009 [Google Scholar]
  47. Koumpia, E., Evans, L., Di Francesco, J., van der Tak, F. F. S., & Oudmaijer, R. D. 2020, A&A, 643, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Krumholz, M. R., Crocker, R. M., & Offner, S. S. R. 2023, MNRAS, 520, 5126 [NASA ADS] [CrossRef] [Google Scholar]
  49. Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kuffmeier, M., Zhao, B., & Caselli, P. 2020, A&A, 639, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2016, ApJ, 820, L2 [NASA ADS] [CrossRef] [Google Scholar]
  53. Lefèvre, C., Cabrit, S., Maury, A. J., et al. 2017, A&A, 604, L1 [CrossRef] [EDP Sciences] [Google Scholar]
  54. Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266 [Google Scholar]
  55. Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956 [Google Scholar]
  56. Maret, S., Bergin, E. A., Neufeld, D. A., et al. 2009, ApJ, 698, 1244 [NASA ADS] [CrossRef] [Google Scholar]
  57. McDaniel, E. W., & Mason, E. A. 1973, The Mobility and Diffusion of Ions in Gases, Wiley Series in Plasma Physics (Hoboken: John Wiley & Sons) [Google Scholar]
  58. McKee, C. F. 1989, ApJ, 345, 782 [NASA ADS] [CrossRef] [Google Scholar]
  59. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [Google Scholar]
  60. Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126 [Google Scholar]
  61. Mouschovias, T. C., Ciolek, G. E., & Morton, S. A. 2011, MNRAS, 415, 1751 [NASA ADS] [CrossRef] [Google Scholar]
  62. Nakano, T., & Tademaru, E. 1972, ApJ, 173, 87 [NASA ADS] [CrossRef] [Google Scholar]
  63. Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163 [Google Scholar]
  64. Offner, S. S. R., Taylor, J., Markey, C., et al. 2022, MNRAS, 517, 885 [NASA ADS] [CrossRef] [Google Scholar]
  65. Osorio, M., Díaz-Rodríguez, A. K., Anglada, G., et al. 2017, ApJ, 840, 36 [Google Scholar]
  66. Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, Space Sci. Rev., 216, 29 [CrossRef] [Google Scholar]
  70. Padovani, M., Marcowith, A., Galli, D., Hunt, L. K., & Fontani, F. 2021, A&A, 649, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Padovani, M., Bialy, S., Galli, D., et al. 2022, A&A, 658, A189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, ASP Conf. Ser., 534, 193 [NASA ADS] [Google Scholar]
  73. Pezzuto, S., Benedettini, M., Di Francesco, J., et al. 2021, A&A, 645, A55 [EDP Sciences] [Google Scholar]
  74. Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481 [Google Scholar]
  75. Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, ApJ, 712, L116 [Google Scholar]
  76. Pineda, J. E., Schmiedeke, A., Caselli, P., et al. 2021, ApJ, 912, 7 [NASA ADS] [CrossRef] [Google Scholar]
  77. Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, ASP Conf. Ser., 534, 233 [NASA ADS] [Google Scholar]
  78. Plunkett, A. L., Arce, H. G., Corder, S. A., et al. 2013, ApJ, 774, 22 [NASA ADS] [CrossRef] [Google Scholar]
  79. Punanova, A., Vasyunin, A., Caselli, P., et al. 2022, ApJ, 927, 213 [CrossRef] [Google Scholar]
  80. Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P., et al. 2021, A&A, 648, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Rodríguez-Kamenetzky, A., Carrasco-González, C., Araudo, A., et al. 2017, ApJ, 851, 16 [Google Scholar]
  83. Sabatini, G., Bovino, S., & Redaelli, E. 2023, ApJ, 947, L18 [NASA ADS] [CrossRef] [Google Scholar]
  84. Sanna, A., Moscadelli, L., Goddi, C., et al. 2019, A&A, 623, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Segura-Cox, D. M., Harris, R. J., Tobin, J. J., et al. 2016, ApJ, 817, L14 [NASA ADS] [CrossRef] [Google Scholar]
  86. Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161 [Google Scholar]
  87. Sokolov, V., Pineda, J. E., Buchner, J., & Caselli, P. 2020, ApJ, 892, L32 [NASA ADS] [CrossRef] [Google Scholar]
  88. Soler, J. D., Zucker, C., Peek, J. E. G., et al. 2023, A&A, 675, A206 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971 [Google Scholar]
  90. Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Hoboken: Wiley-VCH) [Google Scholar]
  91. Stephens, I. W., Dunham, M. M., Myers, P. C., et al. 2017, ApJ, 846, 16 [Google Scholar]
  92. Storm, S., Mundy, L. G., Lee, K. I., et al. 2016, ApJ, 830, 127 [CrossRef] [Google Scholar]
  93. Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, ASP Conf. Ser., 534, 317 [NASA ADS] [Google Scholar]
  94. Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617 [NASA ADS] [Google Scholar]
  95. Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103 [NASA ADS] [CrossRef] [Google Scholar]
  96. van der Tak, F. F. S., Caselli, P., & Ceccarelli, C. 2005, A&A, 439, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198 [NASA ADS] [CrossRef] [Google Scholar]
  98. Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66 [Google Scholar]
  99. Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [Google Scholar]
  100. Winston, E., Megeath, S. T., Wolk, S. J., et al. 2010, AJ, 140, 266 [NASA ADS] [CrossRef] [Google Scholar]
  101. Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [Google Scholar]
  102. Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 476, 2063 [CrossRef] [Google Scholar]
  103. Zhang, Y., Higuchi, A. E., Sakai, N., et al. 2018, ApJ, 864, 76 [NASA ADS] [CrossRef] [Google Scholar]
  104. Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050 [Google Scholar]
  105. Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43 [CrossRef] [Google Scholar]
  106. Zucker, C., Schlafly, E. F., Speagle, J. S., et al. 2018, ApJ, 869, 83 [NASA ADS] [CrossRef] [Google Scholar]
4

This is akin as to deriving the typical flux at a given distance.

All Tables

Table 1

Lines observed in the high-spectral resolution band.

In the text
Table D.1

Properties of Class 0 and I protostars within the field of view.

In the text

All Figures

thumbnail Fig. 1

Color image of the NGC 1333 region using Spitɀer IRAC 3.6, 4.5, and 8.0 µm. The image highlights the presence of outflows in a green color. The area covered by the combined NOEMA and 30 m observations is shown with the dashed box. The scale bar is shown in the bottom-right corner.
In the text
thumbnail Fig. 2

Maps of the DCO+ and H13CO+ (1−0) emission obtained with NOEMA and 30 m telescopes. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the first level is at log10(N(H2)) = 22.2, and the step between levels is 0.5 dex. The stars mark the positions of the YSOs identified by Dunham et al. (2015) using Spitɀer observations.
In the text
thumbnail Fig. 3

Map of the Herschel-based H2 column map, N(H2). The contours are drawn starting at log10(N(H2)) = 22.2 and with a step between levels of 0.5 dex. These contours are the same as in Fig. 2. The beam size and scale bar are shown in the top left and bottom right, respectively.
In the text
thumbnail Fig. 4

Integrated intensity map of the C18O (3−2) transition line. The contours correspond to the Herschel-based H2 column map, with the same contours as in Fig. 2. The beam size and scale bar are shown in top left and bottom right, respectively.
In the text
thumbnail Fig. 5

Map of the CO depletion fraction, fD. The beam size and scale bar are shown in the top-left and bottom-right corners, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2.
In the text
thumbnail Fig. 6

Volume density map estimated from Herschel H2 column density. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2.
In the text
thumbnail Fig. 7

Derived electron fraction in the region. Left: map of the electron fraction, X(e). The map shows a smooth distribution but with a clear enhancement in the northwest section. The beam size and scale bar are shown in the top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. Right: Kernel density estimate of the X(e) distribution across the mapped region. A typical value (median of the log10) of 10−6.5 was derived from these measurements and is marked with a red dashed vertical line.
In the text
thumbnail Fig. 8

Derived cosmic-ray ionization rate in the region. Left: map of the cosmic-ray ionization rate, ζ(H2). The map shows a smooth distribution but with the value clearly increasing in the northwest section. The beam size and scale bar are shown in top left and bottom right, respectively. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. Right: Kernel density estimate of the ζ(H2) distribution across the mapped region. The median value of the log10 ζ(H2)/s−1 = −16.5 and the canonical value of 10−17 s−1 are marked with the blue and black dashed vertical lines, respectively.
In the text
thumbnail Fig. 9

Electron fraction, X(e), as a function of H2 density. The KDE of the underlying distribution was obtained from the data. The contour levels are drawn starting at 0.5σ and progressing outward in steps of 0.5σ where the σ levels are equivalent to that of a bivariate normal distribution. The analytic relations from McKee (1989) and Caselli et al. (2002b) as well as the result from chemical modeling of L1544 (Redaelli et al. 2021) are shown with orange, blue, and green curves, respectively. The analytic relations shown here do not take into account the local generation of cosmic rays.
In the text
thumbnail Fig. 10

Outflows and magnetic field orientation over ζ(H2). The background image is the derived ζ(H2) map, as shown in Fig. 8. The region identified with a local cosmic-ray ionization rate enhancement is marked with an orange contour. The outflow lobes (as reported in previous works) are shown by red and blue arrows. The magnetic field orientation, measured from continuum polarization observations with SCUBA2-POL, are shown by the segments of white lines. The X-ray sources are marked with red crosses. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.
In the text
thumbnail Fig. 11

Cosmic-ray ionization rate as a function of distance to YSOs. The KDE of the CRIR and the typical distance to the YSOs, dYSO, was calculated over the map. The dotted and dashed lines show the inner constant section, and the two different power-law decays help guide the eye.
In the text
thumbnail Fig. A.1

Optical depth of C18O (3−2). The red line shows the τ(C18O) = 1 contour. The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.
In the text
thumbnail Fig. C.1

Difference between observed CRIR and the model-ℒ to better illustrate the local effects, ∆ log10 ζ = log10 ζobs − log10 ζmodel The red and blue colors show the excess and deficiency of CRIR with respect to the expected value following Padovani et al. (2018). The contours correspond to the Herschel-based H2 column map; the contours are the same as in Fig. 2. The beam size and scale bar are shown in the top-let and bottom-right corners, respectively.
In the text

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

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

Initial download of the metrics may take a while.