Effective Shielding of $\lesssim$ 10 GeV Cosmic Rays from Dense Molecular Clumps

We chose Taurus and Perseus clouds to search for the effect of possible shielding of CRs in dense clumps. Located in different layers of the Taurus-Auriga-Perseus complex, the largest local associations of dark nebula [ungerechts87], Taurus and Perseus are among the nearest GMCs to the Earth, whose distances are estimated as $100-200$ pc and $\sim~{}300$ pc, respectively [Dame01]. In the vicinity of these clouds, both low mass and high mass star formations have been observed [walter91]. They are amongst the GMCs with the largest $M/d^{2}$ values, where $M$ is the mass and $d$ is the distance to the solar system. This parameter is crucial for the predicted gamma-ray flux of the GMCs [yang14], which is proportional to $M/d^{2}$. High gamma-ray detection significance allows more accurate spatial and spectral analysis, which is required in this work.

We use Planck dust opacity maps to trace the gas distributions. A detailed description can be found in the Methods. The derived gas column density map is shown is Fig. LABEL:fig:gasdis, in which we identified six dense clumps with average cubic density higher than $5000\,\rm cm^{-3}$. Of those clumps, five were selected for our investigation (dubbed C1-5), whereas the sixth one is disregarded, as it coincides with the star formation region IC 348. The gamma-ray emission in the direction of IC 348 shows different properties and shall be considered independently (Methods).

We analysed about 12 years of Fermi-LAT gamma-ray data (above $100~{}\rm MeV$ ) towards the Taurus-Perseus region (see Methods for details on the analysis). The observed gamma-ray emission in the field of view results from the CR sea interacting with the molecular content (diffuse background) and individual sources. The former includes the information we want to derive, where we use a template obtained from the Planck dust opacity [planck13_11]. The latter is modelled using the fourth source catalog of Fermi-LAT (4FGL) catalogue [4fgl]. A detailed description of how the diffuse background model was obtained can be found in the Methods. The diffuse gamma-ray model should be much more complex than the gas distributions traced by the Planck dust opacity, but in the high latitude region as the one in here, and given the proximity of the GMCs, the total gas in the line of sight is dominated by the GMCs. In such a scenario, the assumption that the gas distribution has a similar spatial distribution to the gamma-ray emission is a good approximation. The usage of dust opacity as the background template also allows us to separate the dust emission from a different line of sight and derive the CR distribution in different regions. To validate our results, we compared our results with the ones obtained using the galactic diffusion emission model provided by the Fermi collaboration and used in the standard analysis. The consistent existence of the negative residuals showed in Supplementary Fig. LABEL:fig:resmap_fermi and LABEL:fig:resmap_dust demonstrates that the adopted model based on the dust distribution is correct for this region.

We separated the Planck dust opacity template into diffuse envelope and dense cores, of which the gas column density is below and above $1.8\times 10^{22}~{}\rm cm^{-2}$, respectively. We found such separation can improve the likelihood fitting significantly when comparing with the fitting without such separation, which implies that the gamma-ray emission from the diffuse envelope and dense cores may have different spectra. This separation also allows us to derive the diffuse spectra of gamma-rays in regions of lower density (envelope) and dense ones (cores). Then we normalized the gamma-ray emission to the emissivities per H atom, which should be proportional to the CR density. The results are shown in the left panel of Fig. LABEL:fig:sed. We found that the gamma-ray emissivities derived in the diffuse envelope are consistent with the predicted gamma-ray emissivities, assuming the CR spectra are the same as the local interstellar spectra (LIS). For dense cores, however, the measured gamma-ray emissivities exhibit significant localized deficits at lower energies ($E<2~{}\rm GeV$).

One possible explanation of such a deficit of low-energy gamma-rays is that the gas column density of the dense core is overestimated, which lowers the gamma-ray flux consequentially. This would imply a significantly lower gas-to-dust ratio in dense cores. Indeed, such a variation of the gas-to-dust ratio has already been observed inside the dense core $\rm\rho~{}Oph~{}A$ [liseau15]. However, the very low gas-to-dust ratio is only observed in the very central region of the core ($<0.05~{}\rm pc$), which is far below the resolution of gamma-ray instruments. On the scale of the core of several parsecs, the average gas-to-dust ratio estimated in ref [liseau15] is consistent with the average value in our Galaxy. Stronger evidence against the interpretation is that in almost all our cases, the deficit appears to be stronger for CRs of lower energies, and an overestimation of the gas column density would affect the normalization of the spectrum but not the shape. The different spectral shapes (as shown in Fig. LABEL:fig:sed) point rather to different CR spectra at different regions.

We further derived the CR spectra from the observed gamma-ray emissivities. Here, we assumed a smoothed broken power-law spectrum for parent CR protons, that is:

where $N_{p}$ is the proton density, $E$ is the kinetic energy of CR protons, $A_{0}$ is the normalization factor, $a_{1}$ and $a_{2}$ are the spectral indices below and above the cutoff energy $E_{c}$, respectively. We note that the low-energy data points are not very constraining if $a_{1}$ is free, thus, we fixed $a_{1}$ to be 1 in the fitting. The derived $a_{2}$ and $E_{c}$ for diffuse envelope and dense cores are $a_{2}=3.01\pm 0.06,E_{c}=2.71\pm 0.65$ and $a_{2}=2.87\pm 0.10,E_{c}=3.66\pm 1.20$, respectively. The full covariance matrix is used when deriving the CR proton flux, which is shown in the right panel of Fig. LABEL:fig:sed. The indications of the CR spectra are consistent with those of the gamma-ray emissivities, that is, the CR in the diffuse envelope has nearly identical spectrum as the LIS, while in the dense cores, the CR density is significantly lower below a few GeV, and became consistent with LIS at higher energies.

The CR deficit at low energies has three possible origins. Firstly, similar to the ’solar modulation’ effects in the solar system, this may be due to the young stars formed already inside these dense cores producing strong stellar winds, which prevent the CRs in the ISM to enter the dense cores. Such an effect has also recently been adopted to explain the low CR density in the central molecular zone and other regions[yang14, huang21]. The effectiveness of this mechanism depends on the overall strength of the stellar winds driven by young stars. The star formation in Perseus and Taurus differs vastly from one region to other[Mercimek2017], which would result in major differences in the CR deficit from cloud to cloud. With the current data, we did not find significant difference in the gamma-ray spectra for different clumps (Methods), although given the current statistics, we cannot rule out such a scenario yet. Additionally, protostellar outflows are collimated objects, with opening angles ranging from a few degrees for young objects to tens of degrees for evolved objects [2007prpl.conf..245A]. Although these regions might contain a few of these outflows, the outflow cavities are not volume-filling, and their role in CR shielding should be limited. In view of the above arguments, we would not discuss such a scenario in detail.

The magnetic mirroring effects can also prevent CRs entering the dense clumps. As estimated in ref [owen21], the CR flux can be decreased by about $50\%$ in the clumps compared with the cores. However, such an effect is also energy independent and thus contradicts to the different spectral shapes observed by gamma-rays (also see the Methods for details). Thus, the magnetic mirroring cannot account for the entire CR deficit. We note that the effect of the magnetic mirroring is estimated based on the method in [owen21], in which the time evolution of magnetic fields is neglected. Detailed calculations based on more realistic assumptions is a subject of future studies.

A third explanation that seems to be favoured by our observations is a slower diffusive transport inside the dense cores, which has already been studied in ref [gabici07]. Under this scenario, the CR deficit is caused by a combined effect of increased proton-proton absorption and slower diffusion caused by an increased magnetic field [2016MNRAS.459.2432V, 2017MNRAS.464.4096L, 2017MNRAS.464.4096L, 2018MNRAS.474.2167L]. Since these factors are independent on star-forming process, it can explain why CR deficit occurs to all of our clumps. Additionally CRs of higher energies have larger gyroradii, and should diffuse faster [Berezinskii1990book], therefore the deficit would naturally be stronger for protons of lower energies.

To apply this hypothesis to our observations, we describe the CR transport adopting the equation in ref [gabici07], that is, CR transport is dominated by the diffusion equation with energy loss:

where $N$ is the space- and energy-dependent particle distribution function of CRs, $R$ is the distance to the center of the clump, $D$ is the diffusion coefficient, and $\dot{E}$ is the energy loss rate of CR protons including both ionization loss and nuclear interaction loss. We assume a spherical clump with density profile of $n(R)=\frac{8000~{}\rm cm^{-3}}{0.04+(R/1~{}\rm pc)^{2}}$, and impose the boundary condition $N(R_{\rm mc})=N_{\rm lis}$, where $R_{\rm mc}$ presents the clump boundary, and $F_{\rm lis}$ is the LIS CR density and a reflective (symmetric) boundary condition at $R=0$. The normalizations are adjusted such that a total of 1800 $M_{\odot}$ is contained in a region of $\sim 1\;\rm pc$, which is consistent with the clump parameters in Table 1. Applying the above conditions, and for a stationary solution ($\frac{dN}{dt}=0$), we can constrain the diffusion coefficient by fitting the derived CR spectrum. In this work we used a generic form $D(E)=D_{0}(\frac{E}{1~{}\rm GeV})^{\Gamma}$,where $D_{0}$ is the diffusion coefficient 1 GeV at 1 GeV and $\Gamma$ is the index reflecting the energy dependence.. Since $\Gamma$ is highly unknown, we varied it within the reasonable limits to constrain the value of $D_{0}$. Namely, changing $\Gamma$ within the physically realistic limits from 0 (energy-independent diffusion) to 1 (Bohm-type diffusion) we found that $D_{0}$ cannot be outside the interval $2\times 10^{26}-6\times 10^{26}~{}\rm cm^{2}/s$. Thus, we arrive a robust conclusion that the diffusion coefficient around 1 GeV should be smaller than the one in the ISM by two orders of magnitude. Three schematic fittings with different $\Gamma$ are shown in the right panel of Fig. LABEL:fig:sed.

We conclude that our observations required a slower diffusion characterized by a smaller diffusion coefficient. This stands in contrast to what was expected in these regions. According to the prevailing view[cesarsky78], the magnetic turbulence should be damped out in the dense neutral gas environment, such that the propagation inside these dense cores should be dominated by free streaming or advection. The faster transport leads to an effective diffusion coefficient that is larger than those in the ISM [cesarsky78, morlino15], which contracts the observational results. However, recent studies reveal that the magnetic field strength increases with the gas density [Crutcher2012], with, for example, $B\propto\rho^{1/2}$. This can be explained if the energy density of the magnetic field is a fixed fraction ($f_{B}$) of that of the kinematic motion, which is further tied to the potential energy if the system is virialized [2016MNRAS.459.2432V, 2017MNRAS.464.4096L, 2018MNRAS.474.2167L]. The increased magnetic field in denser regions leads to slower CR diffusion. This scenario can potentially explain the decrease of CRs in dense regions. Another possibility is the magnetic turbulence induced by CR streaming instability, which will also result in a slower diffusion inside dense clouds[dogiel18]. Possible tests to distinguish these two scenarios are the GeV gamma-ray observations on clouds near CR accelerators, such as the molecular clouds interacting with supernova remnants[Jiang2010], where the CR density are expected to be higher then the nearby GMCs. If the magnetic turbulence is increased mainly due to CR streaming instability, high CR density would imply more severe shielding. However, such supernova remnant-molecular cloud systems are all quite distant and the angular resolutions of current GeV gamma-ray instruments make such study rather difficult.

The Galactic molecular ISM contains density fluctuations over- all scales. Having established the mechanism of magnetic shielding, we establish a picture of the effectiveness of magnetic shielding under different conditions. We assume that the clumps are virialised and that the magnetic energy density is comparable to the kinetic energy density. The condition for effective shielding can be derived through $\tau_{\rm diffusion}=\tau_{\rm pp}$, as other effects such as advective transport can be neglected [gabici07]. Here $\tau_{\rm diffusion}$ and $\tau_{\rm pp}$ are the diffusion and cooling time scale, respectively. $\tau_{\rm diffusion}\approx r_{\rm clump}^{2}/D$, where $D$ is the diffusion coefficient and $r_{\rm clump}$ is the physical size of the the clump. Since the inelastic scattering cross section of proton-proton collision is only weakly energy-dependent, $\tau_{\rm pp}\propto n_{\rm clump}^{-1}$, where $n_{\rm clump}$ is the average gas density in the clump. We assume $D\propto B^{-\gamma}E^{\gamma}$, where $\gamma<1$ is needed to account for a possible suppression of the diffusion coefficient inside the turbulent cloud (in the following we adopt $\gamma=1/2$) [Berezinskii1990book, Casse2001, Fatuzzo2010].

Using Zeeman measurement, observations [Crutcher1999] found out a positive correlation between the gas density and magnetic field strength, a trend that was later refrained by further observations [Crutcher2012]. This positive correlation can be understood in framework where the magnetic energy density follows the kinetic energy density [2018MNRAS.474.2167L] $B^{2}/8\pi\approx 1/2\,n_{\rm clump}\sigma_{\rm v}^{2}$, where $\sigma_{\rm v}^{2}$ is the velocity dispersion in the clump. If clumps are self-gravitating and virialized, for example, $\sigma_{\rm v}=(Gm_{\rm clump}/r_{\rm clump})^{1/2}$, this scaling relation allows us to estimate the magnetic field strength for clumps of different masses and radii. Thus from $\tau_{\rm diffusion}=\tau_{\rm pp}$ we derive a relation between shielding energy $E_{\rm s}$ and the clump parameters, that is, $E_{\rm s}\propto\frac{m_{\rm clump}^{3}}{r_{\rm clump}^{4}}$ (where $m_{\rm clump}$ is the clump mass), below which the shielding effect is substantial. The details can be found in the section ‘A general model for CR diffusion’ in the Methods. Shielding becomes effective above the threshold of clump mass

In Fig. LABEL:fig:general:diffusion, we plot the thresholds for effective shielding on the mass-size plane, where locations of molecular clouds [2020MNRAS.493..351C] and clumps [2014MNRAS.443.1555U] are also overlaid. The shielding effect is not important on the cloud scales, but becomes notable for the clumps with significantly increased densities.

The CRs inside GMCs determine the ionization rate and temperature of the molecular gas and affect the star formation process. The CR density is a crucial parameter in ISM modeling, yet its value is poorly constrained due to our limited understanding of CR propagation. In this Article, for the first time, we find evidence that GeV and sub-GeV CRs, a dominant source of molecular ionization and a major driver of astrochemistry, do not penetrate deep into the dense core of GMCs, leading to localized negative regions (holes) in the gamma-ray residual maps. The effect is more pronounced at lower CR energies. A direct implication of our results here is the lower ionization rate in these dense clumps. Indeed, the anti-correlation between ionization rates and gas column densities has already been revealed from previous astrochemistry observations[albertsson18], although not for the exact clumps we investigated in this work.

We propose that the CR deficit is caused by slower diffusion in the cloud. This slower diffusion of CRs contradicts the prevailing view where magnetic turbulence is suppressed by neutral-ion damping [cesarsky78], which enables CRs to penetrate freely into the dense cores. On the contrary, the diffusion should be slower in the dense, star-cluster-forming clumps, where the magnetic field is stronger because the energy density of magnetic field follows the energy density of, for example, turbulent motions [2018MNRAS.474.2167L]. This scenario also explains the energy-dependent CR deficit as indicated by our observations.

Another possibility is that there is a small region separating the dense clumps from the more diffuse cloud, possibly due to CR self-generated waves [morlino15, ivlev18, dogiel18], in which the CR diffusion is suppressed. In this case the CR propagation inside clumps can still be dominated by free streaming. The CRs deficit is caused by a barrier at the clump boundary rather than the slow diffusion inside dense clump.

Molecular clouds are objects with significant density fluctuations caused by turbulence and gravity, forming various barriers for CR transport. To describe this process, we propose a self-consistent model to describe the transport of CRs in different regions and find effective shielding occurs at the upper left part of the mass-size plane, above the threshold estimated in Eq.3. The shielding effect is important in dense, gravitationally bound regions which are believed to be star-cluster progenitors. In cloud evolution, gravitational collapse creates high-density regions, where the CR density is reduced and the initial condition of star formation is modified.

The Fermi-LAT data used in this work is provided online by the NASA-GSFC Fermi Science Support Center, and can be downloaded from the data serve https://fermi.gsfc.nasa.gov/ssc/data/access/lat/.The Planck dust opacity map is publicly available in Planck Legacy Archive (http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP_ID=COM_CompMap_ThermalDust-commander_2048_R2.00.fits) .

Fermi-LAT data used in our study were reduced and analysed using the standard FERMITOOLS V1.0.1 software package available from https://github.com/fermi-lat/Fermitools-conda/wiki.

Rui-zhi Yang is supported by the NSFC under grants 12041305, 11421303 and the national youth thousand talents program in China. Guang-Xing Li acknowledges supports from NSFC grant 1227030463, W820301904 and 12033005. Bing Liu acknowledges the support from the NSFC under grant 12103049. Yudong Cui is supported by Ministry of Education of China, and the China Manned Space Project (No. CMS-CSST- 2021-B09).

Rui-zhi Yang, Guang-Xing Li and Bing Liu contributed to the paper in equivalent fractions. Rui-zhi Yang, Bing Liu and Emma de Ona Wilhelmi performed the data analysis, Guang-Xing Li, Rui-zhi Yang, Yu-dong Cui and Felix Aharonian were responsible for the interpretation part, all authors contributed to the manuscript writing.

The authors declare no competing interests.