Exploration of Cosmic-ray Acceleration in Protostellar Accretion Shocks and a Model for Ionization Rates in Embedded Protoclusters

The physics of CR acceleration has been relatively well understood for decades (i.e., Umebayashi & Nakano 1981; Drury 1983). First-order Fermi acceleration, also known as diffusive shock acceleration (DSA), can work in jet shocks and protostellar accretion shocks to produce high-energy CRs (Padovani et al. 2016). Under this mechanism, CR protons gain energy every time they pass across the shock. If the flow is turbulent, magnetic field fluctuations scatter these protons back and forth across the shock many times, allowing them to continuously gain energy. However, several important conditions must be met for DSA to occur. The flow must be supersonic and super-Alfvénic for there to be sufficient magnetic fluctuations. The acceleration rate must be greater than the collisional loss rate, the wave-dampening rate, and the rate of diffusion in the transverse direction of the shock. Finally, the acceleration time must be shorter than the timescale of the shock. Each of these timescale conditions limits the energy at which the CRs can be accelerated (as discussed below and in Appendix A). We verify that all of these conditions are met throughout our parameter space (see also Padovani et al. 2016).

We describe in detail the relevant physics for accretion shocks in Appendix A following Padovani et al. (2016). Here we give a brief summary of the model. First-order Fermi acceleration leads to a power-law momentum distribution, f(p), where

$$\begin{eqnarray}&&f(p)\propto {p}^{-q}.\end{eqnarray} \tag{ 10 }$$

The physical quantity of most interest in this work is the CR flux spectrum, j(E). The flux spectrum is related to the accelerated number density spectrum, ${ \mathcal N }(E)$, by

$$\begin{eqnarray}&&j(E)=\displaystyle \frac{v(E){ \mathcal N }(E)}{4\pi }\,\,(\mathrm{particles}\,{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}),\end{eqnarray} \tag{ 11 }$$

where v(E) is the velocity as a function of energy in the relativistic limit. The number density spectrum is related to the more fundamental momentum distribution,

$$\begin{eqnarray}&&{ \mathcal N }(E)=4\pi {p}^{2}f(p)\displaystyle \frac{{dp}}{{dE}}\,\,(\mathrm{particles}\,{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-3}),\end{eqnarray} \tag{ 12 }$$

where the momentum distribution power-law index, q, depends on the underlying shock properties. The flux spectrum is defined in an energy range, ${E}_{\mathrm{inj}}\lt E\lt {E}_{\max }$, where E_inj is the injection energy scale of thermal CR particles and E_max is the maximum energy possible for acceleration to be efficient. Here E_inj depends on the strength of the shock and the hydrodynamic properties, such that stronger shocks and stiffer equations of state lead to an injection energy increase. Here E_max is determined by a combination of the magnetic field, ionization fraction, and shock acceleration efficiency. It is important to note that as long as E_max > 1 GeV, any additional increases only weakly affect our results below. We are mainly interested in the effects of ionization produced by CRs, which is dominated by CRs with energies between 100 MeV and 1 GeV (Indriolo et al. 2010). Higher-energy CRs are important to understand and characterize gamma rays or similar high-energy phenomena.

Neutral gas is not only ionized by the primary CRs. Electrons produced by CR ionization can have sufficient energy to cause additional ionizations. We account for secondary electron ionizations following Ivlev et al. (2015), which we discuss in detail in Appendix B. We ignore the effects of primary electron acceleration. Electrons couple more strongly to the magnetic field and have a significantly lower mass than protons, resulting in a maximum energy orders of magnitude below that of protons (Padovani et al. 2016).

The most uncertain parameters are the magnetic field, B, and the shock efficiency parameter, η. The latter represents the fraction of particles accelerated from the thermal population. Before selecting the fiducial values, we explore the impact of each on the results. The magnetic field at the surface of a protostar has not been accurately measured. Theoretical studies of Class II/T Tauri stars predict a surface magnetic field of a few G to 1 kG (Johns-Krull 2007). We vary η between 10⁻⁵ and 10⁻³ and find that the CRIR scales linearly with η and changes E_max by factors of a few. However, a value of η = 10⁻³ is an extreme case. We fix B = 10 G and η = 10⁻⁵ following Padovani et al. (2016). Table 1 summarizes our fiducial physical parameters. We discuss the effects of varying these parameters in Section 4.

Table 1. Model Parameters

Parameter	Fiducial Value	Range
${\mu }_{{\rm{I}}}$	0.6 (ionized)	⋯
${\mu }_{M}$	2.8 (neutral molecular)	⋯
ρ (Section 2.2)	103 μmH cm−3	⋯
Σcl	1.0 g cm−2	⋯
f	0.1	0.1–0.9
B	10 G	10 G−1 kG
η	10−5	10−5−10−3

Note. Values for the parameters we assume in the model and the ranges we discuss in Section 4.

Download table as:  ASCIITypeset image

The CRIR from protons and secondary electrons as a function of gas column density is

$$\begin{eqnarray}&&\zeta (N)=2\pi \int [j(E,N){\sigma }_{p}^{\mathrm{ion}}(E)+{j}_{e}^{\sec }(E,N){\sigma }_{e}^{\mathrm{ion}}(E)]{dE},\end{eqnarray} \tag{ 13 }$$

where $j(E,N)$ is the CR flux at energy E after traveling through the column density, N, and ${\sigma }_{k}^{\mathrm{ion}}$ is the ionization cross section (Padovani et al. 2009). Krause et al. (2015) proposed relativistic corrections to the $p \mbox{-} {{\rm{H}}}_{2}$ ionization cross sections, applicable when E > 1 GeV. We confirmed that the correction factor to the cross section has no impact on our results due to the small population of GeV CRs. Therefore, we use the nonrelativistic cross section for simplicity. At the shock, the CRIR is expected to be considerable. However, as CRs propagate away from the protostars, they undergo two different processes: energy losses due to collisions with matter and geometric dilution. The former directly modifies the spectrum's shape, since the energy loss is not a gray process (with respect to energy). The latter reduces the overall flux. We use the formalism of Padovani et al. (2009) to account for the energy losses from interactions with matter. The loss function is defined by

$$\begin{eqnarray}&&L(E)=-\displaystyle \frac{{dE}}{{dN}({{\rm{H}}}_{2})}.\end{eqnarray} \tag{ 14 }$$

We can calculate the new energy, E_k, after losses as a function of the initial energy, E₀, for a specific column density, N(H₂),

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})=n({{\rm{H}}}_{2})[R({E}_{0})-R({E}_{k})],\end{eqnarray} \tag{ 15 }$$

with the range, R(E), defined as

$$\begin{eqnarray}&&R(E)=\displaystyle \frac{1}{n({{\rm{H}}}_{2})}{\int }_{0}^{E}\displaystyle \frac{{dE}}{L(E)}.\end{eqnarray} \tag{ 16 }$$

The attenuated spectrum is calculated assuming the number of particles is conserved:

$$\begin{eqnarray}&&j^{\prime} ({E}_{k},N({{\rm{H}}}_{2}))=j({E}_{k},N=0)\displaystyle \frac{L({E}_{0})}{L({E}_{k})}.\end{eqnarray} \tag{ 17 }$$

Equation (17) only takes into account interactions with matter. However, the CRs are generated by a point source, so we must also take into account the spatial dilution of the flux. This is in contrast to Padovani et al. (2009), who considered a plane-parallel slab geometry. We account for the spatial dilution by modifying the attenuated flux as

$$\begin{eqnarray}&&j({E}_{k},N({{\rm{H}}}_{2}))=j^{\prime} ({E}_{k},N({{\rm{H}}}_{2})){\left(\displaystyle \frac{{R}_{* }}{({R}_{* }+r(N))}\right)}^{a},\end{eqnarray} \tag{ 18 }$$

where $r(N({{\rm{H}}}_{2})$ is the radius at which the gas has column density N(H₂), and a is the power-law index for how fast the flux is diluted. A full solution of the CR transport equation is needed to properly determine a. However, we take a = 2, corresponding to free streaming, as a lower limit for the CRIR, which is a common assumption (e.g., Turner & Drake 2009; Rab et al. 2017). Observations of ions in protostellar envelopes may be able to constrain the transport further, primarily whether they undergo free streaming (a = 2) or diffusive (a = 1) transport. We discuss the implications of different transport regimes in Section 4.1.4.

The H₂ column density, N(H₂), from the embedded protostar to the edge of the core is the final piece needed to relate the protostellar feedback to the natal cloud environment. We use the McKee & Tan (2003) model describing protostellar accretion from a turbulent core to calculate appropriate column densities. For a dense core embedded in a turbulent star-forming clump, the envelope column density and core radius are given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{core}}=1.22{{\rm{\Sigma }}}_{\mathrm{cl}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&N{({{\rm{H}}}_{2})}_{\mathrm{core}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{core}}}{{\mu }_{M}{m}_{{\rm{H}}}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{core}}=0.057{{\rm{\Sigma }}}_{\mathrm{cl}}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{{m}_{f}}{30{M}_{\odot }}\right)}^{\tfrac{1}{2}}\,\mathrm{pc},\end{eqnarray} \tag{ 21 }$$

where Σ_cl is the surface density of the embedding clump, which is the normalization factor in ${\dot{m}}_{\mathrm{TTC}}$, and N(H₂) is the H₂ column density.

Our cluster results do not depend on an assumed density profile. However, we adopt a density profile to calculate how the CRIR changes within a protostellar envelope (see Section 3.2.1). We calculate the radius for a given column density by assuming a polynomial density distribution, $n(r)={n}_{s}{\left(\tfrac{{R}_{\mathrm{core}}}{r}\right)}^{-{\kappa }_{\rho }}$, where n_s is the number density at the surface of the core and ${\kappa }_{\rho }=\tfrac{3}{2}$ is motivated by McKee & Tan (2003). The column density as measured from the protostar follows from ${N}_{{{\rm{H}}}_{2}}(r)=N{({{\rm{H}}}_{2})}_{\mathrm{core}}-{\int }_{r}^{{R}_{\mathrm{core}}}n(r){dr}$. Inversion results in $r(N({{\rm{H}}}_{2}))={\left(\tfrac{2{n}_{s}{R}_{\mathrm{core}}^{\tfrac{3}{2}}}{N({{{\rm{H}}}_{2})}_{\mathrm{core}}-2{n}_{s}{R}_{\mathrm{core}}-N({{\rm{H}}}_{2})}\right)}^{2}$.

The total CRIR produced by a forming star cluster is calculated by

$$\begin{eqnarray}&&\zeta ({N}_{* })=\displaystyle \sum _{i}^{{N}_{* }}\int [{j}_{i}(E,{N}_{i}){\sigma }_{p}^{\mathrm{ion}}(E)+{j}_{i,e}^{\sec }(E,{N}_{i}){\sigma }_{e}^{\mathrm{ion}}(E)]\,{dE},\end{eqnarray} \tag{ 22 }$$

where N_i is the H₂ column density from the protostar to the surface of the core.

The CR spectrum of protostars is an observational unknown. The unattenuated spectrum is impossible to constrain because CRs quickly interact with matter, both neutral (in the form of excitations and ionizations) and ionized (through Coulomb interactions). Since protostars are embedded within their natal envelope, their radiation is heavily reprocessed by the surrounding dust. Current observations cannot differentiate between the CRs accelerated by Galactic sources and the Sun versus protostellar sources. Previous studies of young stellar objects have therefore used scaled versions of the local solar spectrum (Rab et al. 2017; Rodgers-Lee et al. 2017). In this section, we present predictions for the CR flux spectrum both at the protostellar surface and at the edge of its core.

Figure 1 shows the CR spectrum generated by the accretion shock for a protostar with an instantaneous mass m = 0.5 M_⊙ as a function of its final mass, taking into account both protons and secondary electrons assuming Σ_cl = 1.0 g cm⁻². The unattenuated spectrum shows a clear power-law behavior with an index of −1.9. In the strong shock regime, the index in Equation (10) asymptotically approaches q = 4. The energy spectrum scales as ${p}^{2}f(p)$, such that $j(E)\propto {p}^{-2}$ (Amato 2014), which is consistent with our result. The unattenuated spectrum is well described as a product of an efficient strong shock at the protostellar surface. We note that the final mass dependence acts largely to scale the spectrum through the accretion rate.

Zoom In Zoom Out Reset image size

The injection energy of 1 keV corresponds to an ionized plasma with a temperature of roughly $1.5\times {10}^{6}$ K. The proton spectrum shows that the maximum energy weakly scales with the final mass of the protostar (or the accretion rate). The spectrum tapers to q = 3.0 at high energies due to the acceleration inefficiency at such relativistic speeds. The energy corresponding to the turnover for all spectra, E ≈ 1 GeV, is the transition where $E={m}_{p}{c}^{2}$. The secondary electron spectrum likewise shows qualitatively similar behavior.

The attenuated spectrum in Figure 1 shows very different behavior. Interactions with the dense core greatly alter the shape of the spectrum, and the radial distance traveled significantly reduces the flux. Shortward of 1 GeV, ionizations and excitations effectively flatten the spectrum and shift higher-energy CRs to lower energies. Losses due to pions are minimal due to the lack of CRs above 10 GeV. From 100 MeV to 1 GeV, the proton flux spectrum exhibits a power-law index of q = 2.5.

The secondary electron spectrum shows similar features from collisional losses. However, the interactions of higher-energy CRs enhance the secondary electron spectrum such that there are significantly more lower-energy electrons. At E = 1 keV, for every CR proton, there are 10⁴ secondary electrons due to the interactions of higher-energy CRs, which are less affected by collisional losses.

Figure 2 shows the maximum energy of CR protons as a function of protostellar and final mass and the dominant constraint on acceleration. We find for protostars with $m\,\gt 1$ M_⊙, CR protons have maximum energies greater than 10 GeV. Only protostars with M < 0.1 M_⊙ have sub-GeV maximum energies. The maximum energy for solar- and supersolar-mass protostars is a constant E_max = 17 GeV. This behavior changes at the transition between different acceleration constraints. The CR acceleration in subsolar-mass protostars is constrained by upstream diffusion. In this process, CRs are lost by diffusion upstream at the shock, thus inhibiting the maximum possible energy. At greater masses, the constraint is set by interactions with neutral gas near the shock. The collisional timescale becomes less than the time to accelerate CRs with E > 17 GeV.

Zoom In Zoom Out Reset image size

The attenuated spectrum is only weakly affected by the cluster mass surface density, Σ_cl (Equations (5) and (19)). While a drop by a factor of 10 in Σ_cl produces a reduction by a similar factor in the unattenuated spectrum, the lower column results in a decline of a factor of only a few in the attenuated spectrum. It is also important to note that the unattenuated spectrum here is from the protostellar surface, while previous theoretical models have calibrated their CR spectrum from terrestrial or space-based measurements (Rab et al. 2017; Rodgers-Lee et al. 2017) and correct for geometric attenuation. However, it is difficult to correct for the effects of matter interactions.

In order to properly model protostellar cores and describe their dynamical state, various pressures must be taken into account. We calculate the CR pressure, P_CR, from the energy flux spectrum,

$$\begin{eqnarray}&&{P}_{\mathrm{CR}}=\displaystyle \frac{4\pi }{3}\int p(E)j(E){dE},\end{eqnarray} \tag{ 23 }$$

where p(E) is the relativistic momentum. Figure 3 shows the CR pressure across the parameter space of instantaneous mass, m, and final mass, ${m}_{f}$, assuming Σ_cl = 1.0 g cm⁻². The unattenuated CR pressure is of order 1 dyne cm ⁻² for most of the parameter space. There is a discrete change in the pressure at 3 M_⊙ caused by the similarly discrete radius change (itself brought on by a change in the internal structure of the protostar; Offner & McKee 2011). The pressure, in general, increases with mass. The maximum occurs when $m\approx 10$ and ${m}_{f}=100$ M_⊙. The attenuated spectrum shows a significant decrease in the pressure: by 13 orders of magnitude. The gradient inverts, and the highest CR pressures occur toward the $m={m}_{f}$ boundary and highest $(m,{m}_{f})$. This is due to the change in the radius of the core. For a given final instantaneous mass, the core is physically smallest when the final mass is smallest. The discrete change at 3 M_⊙ is still apparent, although it is less significant.

Zoom In Zoom Out Reset image size

The ratio of the CR pressure to the kinetic pressure is an important test of the model. Figure 4 shows the ratio ${P}_{\mathrm{CR}}/{P}_{\mathrm{kin}}$ across the $(m,{m}_{f})$ parameter space. We compare the unattenuated CR pressure to the ram pressure of the accreting matter, ${P}_{\mathrm{kin}}={\rho }_{s}{v}_{s}^{2}$. Across the parameter space, P_CR is approximately a millionth of the kinetic pressure. Therefore, it is negligible compared to the accretion, as expected. The important kinetic pressure for the attenuated CR pressure is the surrounding molecular cloud turbulent pressure, ${P}_{\mathrm{kin}}\,={{\rm{\Phi }}}_{\mathrm{core}}{{\rm{\Phi }}}_{s}\,G{{\rm{\Sigma }}}_{\mathrm{cl}}^{2}$, with ${{\rm{\Phi }}}_{\mathrm{core}}=2$ and Φ_s = 0.8 following McKee & Tan (2003). The maximum value of the ratio throughout the parameter space is only ${P}_{\mathrm{CR}}/{P}_{\mathrm{kin}}\approx {10}^{-6}$. The CR pressure of CRs leaving the core is negligible to the dynamics of the surrounding molecular cloud, as expected.

Zoom In Zoom Out Reset image size

The CRIR is one of the key parameters of any astrochemical model, controlling the ionization fraction of gas with ${A}_{V}\gt$ a few, where the external FUV cannot penetrate. Figure 5 shows the CRIR, ζ, as a function of $(m,{m}_{f})$ for a single protostar. The same discrete jump at 3 M_⊙ appears due to the radius discontinuity discussed in Section 3.1.2.

Zoom In Zoom Out Reset image size

The unattenuated CRIR, near the protostellar surface, is incredibly high. Most of the parameter space exhibits ζ = 0.1–1 s⁻¹. This value serves as an initial condition to scale the CRIR throughout the protostellar core. The attenuated CRIR, on the right side of Figure 5, shows much more modest values. The attenuated CRIR at the surface of the core varies between 10⁻¹⁷ and 10⁻¹⁹ s⁻¹. The reduction is due in part to the radial dilution (decreasing the overall flux) and the collisional losses (moving 100 MeV–1 GeV protons to lower energies ionize less efficiently). At the surface of the core, the CRIR produced by an individual protostar becomes comparable to the attenuated external CRIR (Padovani et al. 2009). Therefore, it is likely that in star-forming clouds, gas near embedded protostars may be equally affected by external and internal CR sources.

There is a large difference between the unattenuated and attenuated CRIR: 17 orders of magnitude. Figure 6 shows the CRIR as a function of column density for a protostar with m = 0.5 and ${m}_{f}=1.0$ M_⊙ as the solid black line. There is a near power-law behavior showing a 6 dex decrease in ζ with a 5 dex decrease in N(H₂). The column density is a proxy for the distance from the central protostar. As such, different molecules used to constrain ζ may suggest very different values of ζ depending on what radial surface they trace.

Zoom In Zoom Out Reset image size

Figures 1–6 assume a fixed value of Σ_cl = 1 g cm⁻³. Here, Σ_cl has a linear relationship with the unattenuated CRIR. A decline of a factor of 10 in Σ_cl incurs a similar factor of 10 decrease in the unattenuated ζ and P_CR due to the decrease in the accretion rate through ${\dot{m}}_{\mathrm{TC}}$. However, there is a much weaker dependence of Σ_cl on the attenuated CRIR and CR pressure. The core radius depends on ${{\rm{\Sigma }}}_{\mathrm{cl}}^{\tfrac{1}{2}}$, and the core molecular column density N(H₂) scales linearly with Σ_cl. Together, these factors result in only a factor of a few decrease in ζ and P_CR with an order of magnitude decrease in Σ_cl. We present ζ(N) for the whole $(m,{m}_{f}$) space and for different values of Σ_cl in an interactive online tool.³

We have so far presented results for individual protostars within their natal core. However, molecular clouds form many stars simultaneously. We have shown in Section 3.2.1 that at the protostellar core surface, the CRIR can be on par with the attenuated external CRIR. This suggests that it is important to consider both CRIR components in order to understand cloud chemistry in forming clusters.

Figure 7 shows the attenuated CRIR due to all of the embedded protostars in a cluster. The size of the points indicates the number of protostars, and the color of the points indicates the assumed star formation efficiency, which impacts the result through Σ_cl. Figure 7 represents 400 mock clusters covering a large range of $({N}_{* },{\epsilon }_{g}$). The error bars indicate the 1σ spread due to sampling the bivariate PMF. For ${N}_{* }\gt 500$, the error bars are smaller than the data points due to more complete sampling of the bivariate PMF.

Zoom In Zoom Out Reset image size

The two parameters, _g and N_*, produce opposite trends in the CRIR. A reduction in _g leads to a higher Σ, thus causing a greater CRIR. However, this effect is sublinear—a dex change in _g leads to less than a dex change in the cluster CRIR. The CRIR depends more strongly on N_* than _g. Increasing N_* leads to a slightly superlinear increase in ζ due to the inclusion of more high-mass protostars.

We fit ζ(N_*, _g) with a two-dimensional linear function in log space, which represents the model results well:

$$\begin{eqnarray}&&\mathrm{log}\zeta =-0.24\mathrm{log}{\epsilon }_{g}+1.24\mathrm{log}{N}_{* }-19.56.\end{eqnarray} \tag{ 24 }$$

We plot this function in Figure 7 as the dashed lines, and we add lines of constant N_* for reference. For clusters with ${N}_{* }\gt$ few hundred protostars, the CRIR due to embedded protostars is greater than the typically assumed fiducial rate of ζ₀ = 3 × 10⁻¹⁷ s⁻¹, which is shown as the gray horizontal line.

The magnetic field strength at the protostellar surface is not well constrained. Typically, it is thought to range between a few G and 1 kG. In our fiducial model, we assume B = 10 G. However, this is on the smaller end of the possible range. The magnetic field plays a dominant role in setting the maximum energy due to wave dampening. Wave dampening is very sensitive to the magnetic field. A small increase in the magnetic field leads to significantly more dampening of self-produced Alfvén waves by the CRs, described in detail in Appendix A. The dampening criterion in Equation (34) depends on ${B}^{-4};$ a factor of 100 difference in magnetic field strength yields a substantial change in this criterion.

We recalculated the CR spectrum and CRIR for the $(m,{m}_{f})$ parameter space with B = 1 kG. Figure 8 shows the maximum energy and acceleration constraints as a function of $(m,{m}_{f})$. We find that there are swaths of the parameter space where the shock density, temperature, and velocity are such that wave dampening becomes the dominant constraint in acceleration. Figure 8 shows that in these regions, E_max is reduced to 50–100 MeV, significantly below the GeV energy scale in our fiducial model. This has relatively little effect on the unattenuated CRIR. However, for high column densities, the collisional losses are sufficient to significantly reduce the high CR flux. Therefore, there are regions within the $(m,{m}_{f})$ parameter space where the attenuated CRIR will be negligible due to wave dampening. These regions only account for a few percent of the parameter space, i.e., mainly low-mass protostars, so that our cluster results are largely independent of the assumed magnetic field strength.

Zoom In Zoom Out Reset image size

The accretion flow filling fraction, f, directly influences the shock density as $n\propto {f}^{-1}$. Our fiducial model assumes f = 0.1. However, Class 0 sources, which likely have higher accretion rates, may undergo more spherical accretion. We investigate the effect of increasing the shock filling fraction to f = 0.9. Figure 9 shows the maximum energy and acceleration constraint mechanisms for f = 0.9. The behavior is very similar to the fiducial values in Figure 2, although the region where the acceleration is constrained by matter interactions is smaller and pushed toward higher masses. We find that a factor of 9 increase in f leads to a factor of 3–4 decrease in E_max for protostars with masses below 3 M_⊙. The maximum energy for the rest of the parameter space remains above 10 GeV. For this higher filling fraction, we find a factor of 9 decrease in the unattenuated and attenuated CRIR due to a change in the normalization of f(p) (see Equation (10)).

Zoom In Zoom Out Reset image size

While variation in the filling fraction leads to a respectively linear change in the CRIR, in a cluster environment, higher f values for young sources may cancel with lower values exhibited by older sources, whose accretion has declined.

The shock efficiency parameter, η, which describes the fraction of thermalized protons that are accelerated to relativistic speeds, is also poorly constrained. We assume η = 10⁻⁵ for the fiducial model following Padovani et al. (2016). The normalized CR pressure, ${\tilde{P}}_{\mathrm{CR}}=\tfrac{{P}_{\mathrm{CR}}}{{\rho }_{s}{v}_{s}^{2}}$, depends linearly on η (see Appendix A for details). Stronger CR pressure, in relation to the shock ram pressure, decreases the effectiveness of wave dampening. As such, the maximum energy is constrained mainly by interactions with neutral gas. Figure 10 shows the maximum energy and acceleration constraints for η = 10⁻³. We find that the maximum energy for the whole parameter space is approximately 17 GeV. A 2 dex increase in η leads to a 2 dex increase in the unattenuated and attenuated CRIR; i.e., ζ depends linearly on η. Consequently, uncertainties in η lead to large uncertainties in the expected CRIR. However, evidence for CRs with energies greater than 10 GeV would be an indication of a higher η.

Zoom In Zoom Out Reset image size

How CRs are transported through a protostellar core from its central protostar has not been modeled in detail. Transport of CRs is determined by factors relating to the gas density, magnetic field configuration, diffusion coefficients, and CR energy (Padovani et al. 2013; Rodgers-Lee et al. 2017). Figure 11 shows a comparison between the two limiting cases of transport through the core for a subsolar Class 0 protostar. The CRIR is five orders of magnitude higher in the diffusive regime than the free streaming. At the edge of the core, the CRIR is ζ = 10⁻¹¹ s⁻¹. Balancing CR heating, ${{\rm{\Gamma }}}_{\mathrm{cr}}$, with atomic and molecular line cooling, given by Goldsmith (2001), predicts temperatures of $T\gt {10}^{3}$ K for densities of $n\,={10}^{3}$ cm⁻³ and a CRIR ζ = 10⁻¹¹ s⁻¹. Such temperatures at the core edge are inconsistent with observations. However, ζ = 10⁻¹⁷ s⁻¹, the case of free streaming, produces temperatures of $T\approx 10$ K. Observations of molecular ions can measure the CRIR in the outer regions of cores, constraining the transport mechanism. We discuss this in Section 4.2.

Zoom In Zoom Out Reset image size

The case of transport in protostellar disks is much more complicated. The transport through the disk strongly depends on assumptions about the magnetic field morphology (Padovani et al. 2018). In contrast, protostellar core magnetic fields are thought to exhibit an hourglass morphology (Machida et al. 2007; Crutcher 2012). Such a morphology will allow free streaming, although not fully isotropically.