# Probing the physics of star formation (ProPStar). I. First resolved maps of the electron fraction and cosmic-ray ionization rate in NGC 1333

Jaime E. Pineda, Olli Sipilä, Dominique M. Segura-Cox, Maria Teresa Valdivia-Mena, Roberto Neri, Michael Kuffmeier, Alexei V. Ivlev, Stella S. R. Offner, Maria Jose Maureira, Paola Caselli, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Jaime E. Pineda, Olli Sipilä, Dominique M. Segura-Cox, Maria Teresa Valdivia-Mena, Roberto Neri, et al. Probing the physics of star formation (ProPStar). I. First resolved maps of the electron fraction and cosmic-ray ionization rate in NGC 1333. Astronomy & Astrophysics - A&A, 2024, 686, 10.1051/0004-6361/202347997. insu-04836863

HAL Id: insu-04836863

https://insu.hal.science/insu-04836863v1

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing the physics of star formation (ProPStar)

## I. First resolved maps of the electron fraction and cosmic-ray ionization rate in NGC 1333 ⋆, ⋆⋆

Jaime E. Pineda<fnref n="1" />, Olli Sipilä<fnref n="1" />, Dominique M. Segura-Cox<fnref n="2" />,<fnref n="1" />****, Maria Teresa Valdivia-Mena<fnref n="1" />, Roberto Neri<fnref n="3" />, Michael Kuffmeier<fnref n="4" />,<fnref n="1" />, Alexei V. Ivlev<fnref n="1" />, Stella S. R. Offner<fnref n="2" />, Maria Jose Maureira<fnref n="1" />, Paola Caselli<fnref n="1" />, Silvia Spezzano<fnref n="1" />, Nichol Cunningham<fnref n="5" />, Anika Schmiedeke<fnref n="1" />, and Mike Chen<fnref n="6" />

Received 17 September 2023 / Accepted 15 February 2024

## ABSTRACT

*Context.* Electron fraction and cosmic-ray ionization rates in star-forming regions are important quantities in astrochemical modeling and are critical to the degree of coupling between neutrals, ions, and electrons, which regulates the dynamics of the magnetic field. However, these are difficult quantities to estimate.

*Aims.* We aim to derive the electron fraction and cosmic-ray ionization rate maps of an active star-forming region.

*Methods.* We combined observations of the nearby NGC 1333 star-forming region carried out with the NOEMA interferometer and IRAM 30 m single dish to generate high spatial dynamic range maps of different molecular transitions. We used the DCO⁺ and H¹³CO⁺ ratio (in addition to complementary data) to estimate the electron fraction and produce cosmic-ray ionization rate maps.

*Results.* We derived the first large-area electron fraction and cosmic-ray ionization rate resolved maps in a star-forming region, with typical values of 10⁻<fnref n="6" />.<fnref n="5" /> and 10⁻¹⁶.<fnref n="5" /> s⁻<fnref n="1" />, respectively. The maps present clear evidence of enhanced values around embedded young stellar objects (YSOs). This provides strong evidence for locally accelerated cosmic rays. We also found a strong enhancement toward the northwest region in the map that might be related either to an interaction with a bubble or to locally generated cosmic rays by YSOs. We used the typical electron fraction and derived a magnetohydrodynamic (MHD) turbulence dissipation scale of 0.054 pc, which could be tested with future observations.

*Conclusions.* We found a higher cosmic-ray ionization rate compared to the canonical value for N(H₂) = 10²<fnref n="1" />–10²<fnref n="3" /> cm⁻<fnref n="2" /> of 10⁻¹⁷ s⁻<fnref n="1" /> in the region, and it is likely generated by the accreting YSOs. The high value of the electron fraction suggests that new disks will form from gas in the ideal-MHD limit. This indicates that local enhancements of ζ(H₂), due to YSOs, should be taken into account in the analysis of clustered star formation.

**Key words.** astrochemistry – techniques: interferometric – stars: formation – ISM: abundances – cosmic rays – ISM: molecules

## 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, $\zeta(H_2)$. 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) \propto 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) \approx 10^{-8}–10^{-6}$ (Guelin et al. 1977; Caselli et al. 1998;

<footnote n="1">¹ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany
e-mail: jpineda@mpe.mpg.de</footnote>

<footnote n="2">² Department of Astronomy, The University of Texas at Austin, 2500 Speedway, Austin, TX 78712, USA</footnote>

<footnote n="3">³ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 Saint-Martin d'Hères, France</footnote>

<footnote n="4">⁴ Department of Astronomy, University of Virginia, Charlottesville, VA, 22904, USA</footnote>

<footnote n="5">⁵ IPAG, Université Grenoble Alpes, CNRS, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Department for Physics, Engineering Physics and Astrophysics, Queen's University, Kingston, ON, K7L 3N6, Canada</footnote>

<footnote>* The reduced DCO⁺ and H¹³CO⁺ maps and the maps of derived physical parameters are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/686/A162</footnote>

<footnote>** Based on observations carried out under project number S21AD with the IRAM NOEMA Interferometer and 090-21 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).</footnote>

<footnote>*** NSF Astronomy and Astrophysics Postdoctoral Fellow.</footnote>
Bergin et al. 1999; Maret & Bergin 2007, while in the transition from the diffuse to the dense medium in TMC1 ($n(H_2) \sim 10^3 \text{ cm}^{-3}$), an electron fraction of $X(e) \approx 9.8 \times 10^{-8} - 3.6 \times 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) \approx 1-8 \times 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 $H_3^+$ absorption feature toward background sources, which provide a value of $\zeta(H_2) \approx 2 \times 10^{-16} \text{ 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 $\zeta(H_2) \approx 10^{-17} \text{ s}^{-1}$ (Cummings et al. 2016) at column densities between $10^{21}-10^{23} \text{ 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, $\zeta(H_2) \approx 10^{-17} \text{ 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 $NH_3$ or $N_2H^+$, 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 $H^{13}CO^+$ and DCO$^+$ lines. The whole A162, page 2 of 13

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Color image of the NGC 1333 region using Spitzer 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.</figcaption></figure>

<table><caption>**Table 1.** Lines observed in the high-spectral resolution band.</caption><thead><tr><th>Transition</th><th>Rest Freq.<br>(MHz)</th><th>Beam<sup>(a)</sup><br>("x")</th><th>rms<br>(mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td>DCO<sup>+</sup> (1–0)</td><td>72 039.3122</td><td>6.2 × 6.0(−37°)</td><td>23</td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup> (1–0)</td><td>86 754.2884</td><td>5.0 × 4.9(−43°)</td><td>15</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Beam size before smoothing and matching to common beam and grid.

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.

### 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 $H^{13}CO^+$ (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 package<fnref n="1" />. The beam efficiency, $B_{\text{eff}}$, was obtained using the Ruze formula (available in CLASS),

<footnote n="1">¹ http://www.iram.fr/IRAMFR/GILDAS</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Maps of the DCO$^+$ and H$^{13}$CO$^+$ (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 H$_2$ column map; the first level is at $log_{10}(N(H_2)) = 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 *Spitzer* observations.</figcaption></figure>

and it was used to convert the observations into main beam temperatures, $T_{\text{mb}}$. The noise for the H$^{13}$CO$^+$ (1–0) and DCO$^+$ (1–0) cubes is 48 and 79 mK, in $T_{\text{mb}}$ scale, respectively.

## 2.2. NOEMA interferometer

The observations carried out with the IRAM NOthermal 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 $\alpha_{J2000} = 03^h29^m10.2^s$, $\delta_{J2000} = 31^\circ13'49.4''$. We used the PolyFix correlator with a LO 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 H$^{13}$CO$^+$ (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. H$_2$ column density

We used the total H$_2$ column density, $N(H_2)$, derived using the *Herschel* observations (Pezzuto et al. 2021), which are available in the *Herschel* Gould Belt Survey repository$^2$. The effective angular resolution of the $N(H_2)$ 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 C$^{18}$O observations

We used the C$^{18}$O (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 C$^{18}$O data to match the angular resolution of the *Herschel*-based $N(H_2)$. 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,

<footnote>$^2$ http://www.herschel.fr/cea/gouldbelt/en/Phocea/Vie_des_labos/Ast/ast_visu.php?id ast=66</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Map of the Herschel-based H₂ column map, N(H₂). The contours are drawn starting at log₁₀(N(H₂)) = 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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Integrated intensity map of the C¹⁸O (3–2) transition line. The contours correspond to the Herschel-based H₂ 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.</figcaption></figure>

X(e), and the cosmic-ray ionization rate, ζ(H₂), 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) = \frac{2.7 \times 10^{-8}}{R_D} - \frac{1.2 \times 10^{-6}}{f_D}, \qquad (1)$$

$$\zeta(\text{H}_2) = \left[ 7.5 \times 10^{-4} X(e) + \frac{4.6 \times 10^{-10}}{f_D} \right] X(e) n(\text{H}_2) R_H, \quad (2)$$

where $R_D \equiv [\text{DCO}^+]/[\text{HCO}^+]$; $f_D \equiv [\text{[}^{12}\text{CO}]/[\text{H}_2] / ([\text{[}^{12}\text{CO}]/[\text{H}_2])_{\text{fiducial}}$; $([\text{[}^{12}\text{CO}]/[\text{H}_2])_{\text{fiducial}}$ is the expected ¹²CO abundance; $R_H \equiv [\text{HCO}^+]/[\text{[}^{12}\text{CO}]$; and $n(\text{H}_2)$ is the average H₂ number density.

However, since HCO⁺ is usually optically thick and it also traces outflow emission (not traced by DCO⁺), we used H¹³CO⁺ and the canonical isotropic ratio of [¹²C]/[¹³C] = 68 (Milam et al. 2005) to derive $R_D$. Similarly, we used C¹⁸O and the canonical isotropic ratio of [¹²CO]/[C¹⁸O] = 560 (Wilson & Rood 1994) to estimate the ¹²CO column density.

## 3.2. Column densities

### 3.2.1. DCO⁺ and H¹³CO⁺

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:

$$N_{\text{tot}} = \frac{8\pi v^3 Q(T_{\text{ex}})}{c^3 A_{\text{ul}} g_u e^{-E_{\text{up}}/hT_{\text{ex}}} (e^{hv/k_B T_{\text{ex}}} - 1) J(T_{\text{ex}}) - J(T_{\text{bg}})}, \quad (3)$$

where v is the frequency of the transition observed, $A_{ul}$ is the Einstein coefficient for the transition from level u to l, $g_u$ is the degeneracy of level u, $E_{up}$ is the energy of level u, $T_{ex}$ is the excitation temperature, $Q(T)$ is the partition function at temperature T, and $J(T) \equiv \hbar v/k_B (\exp(\hbar v/k_B T) - 1)$ is the Rayleigh-Jeans equivalent temperature. We used the $A_{ul}$, $g_u$, and $E_{up}$ listed in the LAMDA database, and implemented in the molecular_columns package<fnref n="3" />.

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 ($T_{\text{peak}} > 7.5 \times \text{rms}$) in order to obtain robust column densities.

### 3.2.2. C¹⁸O

The C¹⁸O column density was derived using the optically thin approximation,

$$N(\text{C}^{18}\text{O}) = 5 \times 10^{12} \frac{T_{\text{ex}}}{\text{K}} \exp(31.6 \text{ K}/T_{\text{ex}}) \frac{\int T_{\text{mb}} dv}{\text{km s}^{-1}} \text{cm}^{-2}, \quad (4)$$

and an excitation temperature of $T_{ex} = 12$ K, which were previously used by Curtis et al. (2010). These excitation temperature values are similar to those derived using ¹²CO (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 $T_{ex}$ 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(C¹⁸O). Appendix A shows details regarding the optical depth calculation and a figure showing the map.

<footnote n="3">³ https://github.com/jpinedaf/molecular_columns</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Map of the CO depletion fraction, $f_D$. The beam size and scale bar are shown in the top-left and bottom-right corners, respectively. The contours correspond to the Herschel-based H$_2$ column map; the contours are the same as in Fig. 2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Volume density map estimated from Herschel H$_2$ 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 H$_2$ column map; the contours are the same as in Fig. 2.</figcaption></figure>

### 3.3. CO depletion factor

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

$$N(\mathrm{H}_2) = N(\mathrm{C}^{18}\mathrm{O}) \frac{[\mathrm{H}_2]}{[\mathrm{C}^{18}\mathrm{O}]} + \delta_{18}. \quad (5)$$

This offset represents the H$_2$ column in the molecular cloud that is not traced by the C$^{18}$O (3–2) emission, and therefore it needed to be removed from the abundance calculations. We estimated the offset, $\delta_{18} = 2.7 \times 10^{21} \mathrm{cm}^{-2}$, as the median value of the difference observed assuming an undepleted abundance in the column density range of $3 \times 10^{21} \mathrm{cm}^{-2} < N(\mathrm{H}_2) < 5 \times 10^{21} \mathrm{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(\mathrm{H}_2) = \frac{N(\mathrm{H}_2) - \delta_{18}}{\delta L}, \quad (6)$$

where the N(H$_2$) is the H$_2$ column density, $\delta_{18}$ is the derived offset with C$^{18}$O, and $\delta L$ is the depth of the structure. We assumed

$\delta L = 0.4 \mathrm{pc} (\sim 60 \times 10^3 \mathrm{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 $^{12}$CO and $^{13}$CO (1–0) observations (Pineda et al. 2008). The typical density in places with a detection of DCO$^+$ and H$^{13}$CO$^+$ is $\approx 4 \times 10^3 \mathrm{cm}^{-3}$, which suggests that this density is likely a lower limit since this region is bright in NH$_3$ (Friesen et al. 2017; Dhabal et al. 2019) and N$_2$H$^+$ (Hacar et al. 2017), which are reliable high-density tracers.

### 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 $\log_{10}$) of $\langle X(e) \rangle = 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 $R_D$ 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 $\log_{10}$) of $\langle \zeta(\mathrm{H}_2) \rangle = 10^{-16.3} \mathrm{s}^{-1}$ from these measurements. We estimate that the error associated
<figure><img src="image_9.png" /><figcaption>**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 H$_2$ 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 log$_{10}$) of $10^{-6.5}$ was derived from these measurements and is marked with a red dashed vertical line.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Derived cosmic-ray ionization rate in the region. Left: map of the cosmic-ray ionization rate, $\zeta(H_2)$. 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 H$_2$ column map; the contours are the same as in Fig. 2. Right: Kernel density estimate of the $\zeta(H_2)$ distribution across the mapped region. The median value of the log$_{10} \zeta(H_2)/s^{-1} = -16.5$ and the canonical value of $10^{-17} s^{-1}$ are marked with the blue and black dashed vertical lines, respectively.</figcaption></figure>

with using the wrong excitation temperature (with the temperature within the range of 5 and 30 K) $R_H$ could be underestimated by up to a factor of two, in addition to the uncertainty of less than 10% in $R_D$ and $X(e)$ mentioned in Sect. 3.5. This indicates that the $\zeta(H_2)$ could be underestimated by a factor of two toward hot regions ($T_{ex} > 20$K), but elsewhere, these uncertainties are of the order 10–30%.

# 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 $\zeta(H_2)$, 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(H_2)$ and $f_D$ were derived using *Herschel* data, we assumed that they smoothly vary between 16" and $\approx 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(H_2)$, 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(H_2) \gtrsim 10^{3.7} \text{ cm}^{-3}$, the relation steepens. This steeper relation is mostly driven by the north-west 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 ($\approx 10^{3.7} \text{ cm}^{-3}$) even though there is a clear measurement of the depletion of metals ($f_D$).

The derived densities are comparable to previous estimates. The mean density in the region probed with $^{12}$CO and $^{13}$CO (1–0) transitions are $\approx 10^3 \text{ 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} \text{ 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.,

<figure><img src="image_11.png" /></figure>

**Fig. 9.** Electron fraction, $X(e)$, as a function of $H_2$ density. The KDE of the underlying distribution was obtained from the data. The contour levels are drawn starting at $0.5\sigma$ and progressing outward in steps of $0.5\sigma$, where the $\sigma$ 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.

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 $\zeta(H_2)$ occurring in the northwest (upper) section of the map. Furthermore, there is a noticeable enhancement of $\zeta(H_2)$ 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 $\le 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 $\zeta(H_2)$ 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 cosmic-ray 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 $\sim 10^{23} \text{ cm}^{-2}$ (Igea & Glassgold 1999), which is about the measured maximum value of $N(H_2)$ in the upper section (where most of
<figure><img src="image_12.png" /><figcaption>Fig. 10. Outflows and magnetic field orientation over $\zeta(H_2)$. The background image is the derived $\zeta(H_2)$ 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.</figcaption></figure>

the sources are located). Therefore, X-rays would result in a relatively isotropically enhanced $\zeta(H_2)$ around the sources, which is clearly incompatible with the distribution of $\zeta(H_2)$ 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 $\zeta(H_2)$ around most of the YSOs. Given that we strongly underestimated the gas density in the proximity to these objects, the actual values of $\zeta(H_2)$ are expected to be much higher than those shown in Fig. 8. Second, the elongated region with the highest local ionization rate of $\ge 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 distances<fnref n="4" />,

$$ \frac{1}{d_{\text{YSO}}^2} = \frac{1}{N} \sum_{i=1}^{N} \frac{1}{d_i^2}, $$

where $d_i$ is the distance to the different YSOs marked in Fig. 10. Figure 11 shows the KDE for the map, and it shows that within

<figure><img src="image_13.png" /><figcaption>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, $d_{YSO}$, 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.</figcaption></figure>

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 Spitzer 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 SVS13 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 $\zeta(H_2) \sim 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 $\zeta(H_2) \sim 2-5 \times 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

<footnote n="4">⁴ This is akin as to deriving the typical flux at a given distance.</footnote>
measurements of H$_{3}^{+}$ 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\sigma$ upper limit of $9 \times 10^{-17}$ s$^{-1}$ for the ionization rate. Given that the obtained spectra of H$_{3}^{+}$ lines show no visible traces of absorption (see their Fig. 10), we conclude that the actual value of $\zeta(H_2)$ 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.

### 4.4. Evidence of shocks and interaction

The spatial distribution for both of the main parameters estimated in this work, $X(e)$ and $\zeta(H_2)$, 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 $\zeta(H_2)$ 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 $\approx 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 $\langle\zeta(H_2)\rangle > 10^{-14}$ s$^{-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 $\langle\zeta(H_2)\rangle \approx 10^{-17}$ s$^{-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

$$ \lambda_{\min} = \sqrt{\frac{\pi}{4\mu m_{\text{H}} n(H_2)} \frac{B_0}{n(H_2) X(e) \langle\sigma v\rangle}}, \quad (7) $$

where $B_0$ is the unperturbed magnetic field and $\langle\sigma v\rangle$ is the rate coefficient for elastic collisions. The rate coefficient term is approximated by the Langevin term,

$$ \langle\sigma v\rangle = 1.69 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}, \quad (8) $$

for HCO$^{+}$–H$_{2}$ collisions (McDaniel & Mason 1973). We estimated the magnetic field strength, $B_0$, using the relation of Crutcher et al. (2010), a typical volume density of $10^{3.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 $\approx 0.054$ pc (or $\approx 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., $N_2H^+$ and $NH_3$) and angular resolutions showed that the southern end filament presents subsonic levels of turbulence at scales of $\approx 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 $\approx 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

$$ \lambda_{\text{mag}} \equiv \frac{M_c}{\pi R_c^2 B_0} 2\pi \sqrt{G}, \quad (9) $$

where $B_0$ is the magnetic field strength, and $M_c$ and $R_c$ are the mass and radius of the core.

For the typical values of the ionization and cosmic rays found in this work, $\langle X(e) \rangle = 10^{-6.5}$ and $\langle\zeta(H_2)\rangle = 10^{-16.5}$ s$^{-1} = 3 \times 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 $\lambda_{\text{mag}} = 5$ and an initial ratio of rotational to gravitational energy, $\beta$, of the
spherical core between $\beta = 0.025$ and $\beta = 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 $\langle\zeta(H_2)\rangle = 10^{-17} \text{ s}^{-1}$, of a collapsing clump of $1000 \text{ M}_\odot$. They found an average disk radius of $\approx 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 $\zeta(H_2)$

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 $\zeta(H_2)$. However, the $H^{13}CO^+$ and $DCO^+$ (1–0) lines have similar parameters ($E_{up}$ and $A_{ul}$), 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 $H^{13}CO^+$ might trace different volumes or be affected by high optical depths, $H_2D^+$ is a better estimator for the cosmic-ray ionization rate. Unfortunately, this approach is severely limited in low-mass star-forming regions since the $H_2D^+$ 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, $\zeta(H_2)$, “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 $\zeta(H_2)$ 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 $H^{13}CO^+$ (1–0) emission at $\approx 6''$ resolution ($\approx 1$ 800 au).

2. We used the analytic relation from Caselli et al. (1998) relating the measured $DCO^+$ and $H^{13}CO^+$ column density ratio in addition to archival $C^{18}O$ and $H_2$ in order to derive

the first large-area map for the electron fraction, $X(e)$, and cosmic-ray ionization rate, $\zeta(H_2)$, 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} \text{ 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 ($\approx 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 $\zeta(H_2)$ 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.

## References

Beltrán, M. T., Cesaroni, R., Moscadelli, L., et al. 2016, A&A, 593, A49

Bergin, E. A., & Langer, W. D. 1997, ApJ, 486, 316

Bergin, E. A., Plume, R., Williams, J. P., & Myers, P. C. 1999, ApJ, 512, 724

Bron, E., Roueff, E., Gerin, M., et al. 2021, A&A, 645, A28

Cabedo, V., Maury, A., Girart, J. M., et al. 2023, A&A, 669, A90

Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234

Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002a, ApJ, 572, 238

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 565, 344

Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37

Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703

Ceccarelli, C., Dominik, C., López-Sepulcre, A., et al. 2014, ApJ, 790, L1

Ceccarelli, C., Codella, C., Balucani, N., et al. 2023, ASP Conf. Ser., 534, 379

Chen, M. C.-Y., Di Francesco, J., Rosolowsky, E., et al. 2020, ApJ, 891, 84

Chuang, C.-Y., Aso, Y., Hirano, N., Hirano, S., & Machida, M. N. 2021, ApJ, 916, 82

Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466
Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18

Curtis, E. I., Richer, J. S., & Buckle, J. V. 2010, MNRAS, 401, 455

Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47

Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35

De Simone, M., Codella, C., Ceccarelli, C., et al. 2022, MNRAS, 512, 5214

Dhabal, A., Mundy, L. G., Chen, C.-y., Teuben, P., & Storm, S. 2019, ApJ, 876, 108

Diaz-Rodriguez, A. K., Anglada, G., Blázquez-Calero, G., et al. 2022, ApJ, 930, 91

Dionatos, O., Kristensen, L. E., Tafalla, M., Güdel, M., & Persson, M. 2020, A&A, 641, A36

Doi, Y., Hasegawa, T., Furuya, R. S., et al. 2020, ApJ, 899, 28

Doi, Y., Tomisaka, K., Hasegawa, T., et al. 2021, ApJ, 923, L9

Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11

Favre, C., Ceccarelli, C., López-Sepulcre, A., et al. 2018, ApJ, 859, 136

Fitz Axen, M., Offner, S. S. S., Gaches, B. A. L., et al. 2021, ApJ, 915, 43

Fontani, F., Ceccarelli, C., Favre, C., et al. 2017, A&A, 605, A57

Friesen, R. K., Di Francesco, J., Myers, P. C., et al. 2010, ApJ, 718, 666

Friesen, R. K., Di Francesco, J., Bourke, T. L., et al. 2014, ApJ, 797, 27

Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63

Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105

Gaches, B. A. L. & Offner, S. S. R. 2018, ApJ, 861, 87

Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73

Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223

Guelin, M., Langer, W. D., Snell, R. L., & Wootten, H. A. 1977, ApJ, 217, L165

Guelin, M., Langer, W. D., & Wilson, R. W. 1982, A&A, 107, 107

Gutermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ, 674, 336

Hacar, A., Tafalla, M., & Alves, J. 2017, A&A, 606, A123

Harju, J., Haikala, L. K., Lehtinen, K., et al. 2006, A&A, 454, L55

Hodapp, K. W., & Chini, R. 2014, ApJ, 794, 169

Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13

Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848

Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91

Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009

Koumpia, E., Evans, L., Di Francesco, J., van der Tak, F. F. S., & Oudmaijer, R. D. 2020, A&A, 643, A61

Krumholz, M. R., Crocker, R. M., & Offner, S. S. R. 2023, MNRAS, 520, 5126

Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7

Kuffmeier, M., Zhao, B., & Caselli, P. 2020, A&A, 639, A86

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2016, ApJ, 820, L2

Lefèvre, C., Cabrit, S., Maury, A. J., et al. 2017, A&A, 604, L1

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956

Maret, S., Bergin, E. A., Neufeld, D. A., et al. 2009, ApJ, 698, 1244

McDaniel, E. W., & Mason, E. A. 1973, The Mobility and Diffusion of Ions in Gases, Wiley Series in Plasma Physics (Hoboken: John Wiley & Sons)

McKee, C. F. 1989, ApJ, 345, 782

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126

Mouschovias, T. C., Ciolek, G. E., & Morton, S. A. 2011, MNRAS, 415, 1751

Nakano, T., & Tademaru, E. 1972, ApJ, 173, 87

Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163

Offner, S. S. R., Taylor, J., Markey, C., et al. 2022, MNRAS, 517, 885

Osorio, M., Díaz-Rodríguez, A. K., Anglada, G., et al. 2017, ApJ, 840, 36

Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619

Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8

Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111

Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, Space Sci. Rev., 216, 29

Padovani, M., Marcowith, A., Galli, D., Hunt, L. K., & Fontani, F. 2021, A&A, 649, A149

Padovani, M., Bialy, S., Galli, D., et al. 2022, A&A, 658, A180

Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, ASP Conf. Ser., 534, 193

Pezzuto, S., Benedettini, M., Di Francesco, J., et al. 2021, A&A, 645, A55

Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481

Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, ApJ, 712, L116

Pineda, J. E., Schmiedeke, A., Caselli, P., et al. 2021, ApJ, 912, 7

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, ASP Conf. Ser., 534, 233

Plunkett, A. L., Arce, H. G., Corder, S. A., et al. 2013, ApJ, 774, 22

Punanova, A., Vasyunin, A., Caselli, P., et al. 2022, ApJ, 927, 213

Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109

Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P., et al. 2021, A&A, 648, A120

Rodríguez-Kamenetzky, A., Carrasco-González, C., Araudo, A., et al. 2017, ApJ, 851, 16

Sabatini, G., Bovino, S., & Redaelli, E. 2023, ApJ, 947, L18

Sanna, A., Moscadelli, L., Goddi, C., et al. 2019, A&A, 623, L3

Segura-Cox, D. M., Harris, R. J., Tobin, J. J., et al. 2016, ApJ, 817, L14

Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161

Sokolov, V., Pineda, J. E., Buchner, J., & Caselli, P. 2020, ApJ, 892, L32

Soler, J. D., Zucker, C., Peek, J. E. G., et al. 2023, A&A, 675, A206

Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971

Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Hoboken: Wiley-VCH)

Stephens, I. W., Dunham, M. M., Myers, P. C., et al. 2017, ApJ, 846, 16

Storm, S., Mundy, L. G., Lee, K. I., et al. 2016, ApJ, 830, 127

Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, ASP Conf. Ser., 534, 317

Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617

Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103

van der Tak, F. F. S., Caselli, P., & Ceccarelli, C. 2005, A&A, 439, 195

Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198

Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Winston, E., Megeath, S. T., Wolk, S. J., et al. 2010, AJ, 140, 266

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037

Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 476, 2063

Zhang, Y., Higuchi, A. E., Sakai, N., et al. 2018, ApJ, 864, 76

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050

Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43

Zucker, C., Schlafly, E. F., Speagle, J. S., et al. 2018, ApJ, 869, 83
# Appendix A: Optical depth of C$^{18}$O

We estimated the C$^{18}$O (3–2) line optical depth as

$$ \tau(\mathrm{C}^{18}\mathrm{O}) = -\ln\left(1.0 - \frac{T_p(\mathrm{C}^{18}\mathrm{O})}{J_v(T_{ex}) - J_v(T_{cmb})}\right), \quad (\mathrm{A.1}) $$

where $T_p(\mathrm{C}^{18}\mathrm{O})$ is the peak brightness of the line, and $J_v(T)$ is the brightness temperature of a black body with temperature $T$ at frequency $v$. In the case of C$^{18}$O (3–2), we assumed $T_{ex} = 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).

<figure><img src="image_14.png" /><figcaption>Fig. A.1. Optical depth of C$^{18}$O (3–2). The red line shows the $\tau(\mathrm{C}^{18}\mathrm{O}) = 1$ contour. The contours correspond to the Herschel-based H$_2$ 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.</figcaption></figure>

# 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 H$^{13}$CO$^+$ and DCO$^+$, which we used to determine the $R_D$ and $R_H$ quantities. We explored the effect of both species having the same excitation temperature but that temperature being different from 10 K ($5 K < T_{ex} < 12$ K), and we found that the variation in $R_D$ is less than 5%, while for $R_H$, 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: $\Delta \log_{10} \zeta$

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), $\Delta \log_{10} \zeta = \log_{10} \zeta_{obs} - \log_{10} \zeta_{model}$. Figure C.1 shows the map of $\Delta \log_{10} \zeta$, 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.

<figure><img src="image_15.png" /><figcaption>Fig. C.1. Difference between observed CRIR and the model-L to better illustrate the local effects, $\Delta \log_{10} \zeta = \log_{10} \zeta_{obs} - \log_{10} \zeta_{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 H$_2$ 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.</figcaption></figure>

# 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><caption>**Table D.1.** Properties of Class 0 and I protostars within the field of view.</caption><thead><tr><th>Source</th><th>Class</th><th>R<sub>disk</sub><sup>(a)</sup><br>(au)</th><th>PA<sub>outflow</sub><sup>(b)</sup><br>(deg)</th><th>Other Name</th><th>Ref.</th></tr></thead><tbody><tr><td>NGC 1333 IRAS4A1<sup>(c)</sup></td><td>0</td><td>49.7</td><td>19.2</td><td>Per-emb-12</td><td>(1, 5)</td></tr><tr><td>NGC 1333 IRAS4A2<sup>(c)</sup></td><td>0</td><td>–</td><td>20</td><td>Per-emb-12</td><td>(1, 7)</td></tr><tr><td>NGC 1333 IRAS4B<sup>(c)</sup></td><td>0</td><td>–</td><td>176</td><td>Per-emb-13</td><td>(1, 2, 5)</td></tr><tr><td>NGC 1333 IRAS4B2<sup>(c)</sup></td><td>0</td><td>–</td><td>76</td><td>Per-emb-13</td><td>(1, 8)</td></tr><tr><td>NGC 1333 IRAS2B</td><td>I</td><td>–</td><td>24</td><td>Per-emb-36</td><td>(1, 9)</td></tr><tr><td>NGC 1333 IRAS4C</td><td>0</td><td>37.2</td><td>96</td><td>Per-emb-14</td><td>(2, 6)</td></tr><tr><td>SVS13 VLA4A<sup>(d)</sup></td><td>0/I</td><td>12</td><td>8</td><td>SVS13A, Per-emb-44</td><td>(3, 4)</td></tr><tr><td>SVS13 VLA4B<sup>(d)</sup></td><td>0/I</td><td>9</td><td>-50</td><td>SVS13A, Per-emb-44</td><td>(3, 4)</td></tr><tr><td>SVS13 VLA17</td><td>0</td><td>31.7</td><td>170</td><td>SVS13B</td><td>(2, 5)</td></tr><tr><td>SVS13C</td><td>0</td><td>42.9</td><td>8</td><td>–</td><td>(1, 4)</td></tr><tr><td>Per-emb-15</td><td>0</td><td>–</td><td>-35</td><td>SK14</td><td>(4)</td></tr><tr><td>Per-emb-3</td><td>0</td><td>–</td><td>97</td><td>–</td><td>(1, 9)</td></tr></tbody></table>

**References.** (1) Segura-Cox et al. (2018), (2) Segura-Cox et al. (2016), (3) Diaz-Rodriguez et al. (2022), (4) Plunkett et al. (2013), (5) Lee et al. (2016), (6) Zhang et al. (2018), (7) Chuang et al. (2021), (8) Hull et al. (2014), (9) Stephens et al. (2017)

**Notes.** <sup>(a)</sup>The disk radius measurements are dust based, and they were recomputed to the same 300 pc distance, if needed. <sup>(b)</sup>Outflow PA was measured using the red lobe from due east from north. <sup>(c)</sup>These objects form a binary system but do not share circumbinary material. <sup>(d)</sup>These objects belong to a close binary system surrounded by circumbinary material or disk.