Low-energy Cosmic Rays and Associated MeV Gamma-ray Emissions in the Protoplanetary System

1. Introduction

Cosmic rays (CRs) are high-energy particles from outer space, predominantly composed of protons and atomic nuclei [1]. These particles exhibit a wide energy range, from tens of MeV to over ${10}^{20}$ eV. Low-energy CRs (LECRs), defined as CRs with kinetic energies less than approximately 1 GeV per nucleon, have profound impacts on the ionization, heating, and chemical evolution of interstellar medium due to their high energy loss rates within the interstellar medium [2,3,4].

Protoplanetary disks are dense gas and dust disks that form around young stars during the collapse of molecular clouds [5]. These disks, with radii extending up to 1000 AU, are the birthplaces of planetary systems. LECRs significantly impact protoplanetary disks when interacting with the gases and dust in the disk. LECRs can effectively ionize the gas in the protoplanetary disk, altering its ionization level and electrical conductivity, which in turn affects the coupling efficiency of magnetic fields and the excitation of turbulence [6,7,8], changing the dynamic structure and evolutionary path of the disk. Additionally, LECRs provide an extra heating source for the disk, influencing its temperature distribution and geometric shape [9]. Changes in the disk temperature can affect the evolution of its chemical composition and dust properties. Furthermore, LECRs can induce various photochemical reactions in molecules and atoms, altering the chemical composition of the disk [10,11]. These chemical changes might impact the composition of planetary atmospheres and the synthesis of prebiotic molecules. LECRs interacting with surrounding material produce gamma-ray and X-ray radiation, modifying the radiation environment of the protoplanetary disk. This high-energy radiation can affect the ionization and chemical state of the disk’s material. Moreover, LECRs can interact with dust particles, changing their charge, size, and composition [12]. This affects the aggregation of dust and the formation process of planets.

Studying the propagation and spectral evolution of LECRs in protoplanetary disks is thus crucial for understanding the formation and early evolution of planetary systems. However, due to the slow propagation and energy loss of sub-relativistic particles, LECRs are concentrated near their acceleration sites, making them difficult to detect directly. They are often studied through indirect measurements of molecular ionization rate and non-thermal radiation, such as 6.4 KeV Fe K$\alpha$ line emission. Unfortunately, these methods cannot rule out the influence of CR electrons and UV photons, leading to considerable uncertainties in the results. In the protoplanetary disk environment, LECRs interacting with gas and dust may produce a significant amount of nuclear de-excitation line emission with energy mainly within the range of 0.1–10 MeV. The gamma-ray line emission emitted via nuclear de-excitation following the collisions of LECR nuclei with the ambient gases have strong penetration capability and are unaffected by magnetic fields and stellar winds, thus providing a very effective probe to investigate the acceleration and transport processes of LECRs and their composition [13,14,15] and offering important observational evidence for exploring the physical processes and chemical evolution within the protoplanetary disk.

The nuclear de-excitation gamma-ray lines generated during solar flares have been extensively investigated observationally [16,17,18] and theoretically [19,20]. The nuclear de-excitation gamma-ray line emission has also been proposed to be a promising probe to study the LECRs in the interstellar medium [21] and in the vicinity of their sources [22] such as in the supernova remnants [23] and in the Galactic center region [24]. In this study, we investigate the LECRs in protoplanetary disks by the simulation of LECR propagation in protoplanetary disks and the calculation of corresponding nuclear de-excitation line emission, which may shed light on the formation processes of protoplanetary systems from future observations.

The paper is structured as follows: Section 2 simulates the propagation and evolution of LECRs in protoplanetary disks. Section 3 calculates and examines the gamma-ray radiation characteristics resulting from LECR interactions. Section 4 discusses potential improvements to the model. Section 5 provides the summary and outlook.

2. Calculation of LECR Spectrum

This section introduces the numerical methods employed to solve the LECR spectrum, specifically the numerical algorithm and code implementation details for solving the Parker transport equation. The evolution of LECRs is crucial for understanding the chemical processes in protoplanetary disks and the ionization environment of early Earth. Therefore, accurate calculation of the LECR spectrum is essential.

2.1. Theoretical Basis

2.1.1. Parker Transport Equation

CRs from outside the astrosphere are modulated by the stellar wind during their transport in a solar-type stellar system. The transport of CR particles in the interplanetary medium can be described by the Parker equation [25]:

$$\frac{\partial f}{\partial t}+\mathbf{v}·\nabla f-\nabla ·(\kappa \nabla f)-\frac{1}{3}(\nabla ·\mathbf{v})p\frac{\partial f}{\partial p}=0,$$

(1)

where $f(\mathbf{r},p,t)$ is the phase space distribution function of CR particles, $\mathbf{v}$ is the stellar wind velocity field, $\kappa$ is the diffusion coefficient, and p is the momentum of CR particles. The last three terms of LHS in the equation represent spatial advection, spatial diffusion, and adiabatic energy loss, respectively.

While a complete description of the protoplanetary disk and wind system is beyond the scope of the present work, we assume that the stellar wind is spherically symmetric [26] as the lowest order approximation to the system, then the Parker equation reads:

$$\frac{\partial f}{\partial t}+v_r\frac{\partial f}{\partial r}-\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\kappa \frac{\partial f}{\partial r}\right)-\frac{1}{3}\frac{1}{r^2}\frac{\partial \left(r^2{v}_r\right)}{\partial r}p\frac{\partial f}{\partial p}=0,$$

(2)

where r is the radial distance from the central star, $v_r$ is the radial velocity of stellar wind, and p is the momentum of LECRs. We will employ numerical methods to solve this equation and obtain the temporal and spatial evolution as well as the spectral distribution of LECRs in the protoplanetary disk.

2.1.2. Stellar Wind Velocity

The radial velocity of the stellar wind, $v_r$, tends to increase with the stellar rotation rate $\Omega$[26]. At a distance of 1 AU from the star, $v_r$ is adapted in this study from the fitting formula in [27]:

$${log}_{10}{v}_r=a{\left(\frac{\Omega }{{\Omega }_{\odot }}\right)}^b+c{\left(\frac{\Omega }{{\Omega }_{\odot }}\right)}^d+e{\left(\frac{\Omega }{{\Omega }_{\odot }}\right)}^f+g,$$

(3)

where the coefficients $a\sim f$ are given in [27]. For $\Omega <1.4{\Omega }_{\odot }$, an additional term $g=-{log}_{10}1.12$ is included for a better fit to the data from [26]; for $\Omega \ge 1.4{\Omega }_{\odot }$, g is set to 0.

2.1.3. Diffusion Coefficient

The diffusion coefficient $\kappa$ is a critical parameter in the Parker equation, describing the scattering and diffusion process of CR particles in a turbulent magnetic field. For LECRs, due to their small gyroradius, they are more likely to be carried and scattered by stellar winds, resulting in a relatively low diffusion coefficient. This study uses a widely accepted semi-empirical model [26,28,29]:

$$\frac{\kappa (r,p,\Omega )}{\beta c}={\eta }_0{\left(\frac{p}{p_0}\right)}^{1-\gamma }{r}_{\mathrm{L}},$$

(4)

where $\beta =v/c$ is the relativistic factor, $\Omega$ is the stellar rotation angular velocity, ${\eta }_0={\left(\frac{B}{\delta B}\right)}^2$ represents the level of turbulence in the magnetic field, B is associated with the energy density of the large-scale magnetic field, and $\delta B$ is related to the total energy density in the small-scale magnetic turbulence mode. $p_0=3\mathrm{GeV}/c$, $\gamma$ is related to the turbulence power spectrum and is commonly set to $\gamma =1$ based on solar wind modulation observations [30], and $r_{\mathrm{L}}=p/\left[eB(r,\Omega )\right]$ is the Larmor radius of a proton. This model effectively captures the behavior of the diffusion coefficient in the inner interplanetary space.

2.1.4. Boundary Conditions

To solve the Parker equation, appropriate boundary conditions are required. In the momentum space p direction, due to the negligible contribution of extremely high energy CRs to the total spectrum, boundary condition $f(r,p=p_{\max },t)$ = 0 is adopted, where $p_{\max }$ is the maximal momentum in our numerical scheme.

For the spatial coordinate r, a reflective boundary condition is adopted at the inner boundary:

$$f(r=r_{\odot },p,t)=f(r=r_{\odot }+\Delta r,p,t).$$

(5)

This implies that CR particles near the stellar surface are reflected to outer space without being lost to the stellar surface.

At the outer boundary in the r direction, the local interstellar spectrum (LIS) obtained from observations is used as the boundary condition:

$$f(r=r_{\mathrm{LIS}},p,t)=f_{\mathrm{LIS}}\left(p\right).$$

(6)

The LIS refers to the CR spectrum outside the heliosphere, unaffected by stellar wind modulation. Typically, the outer boundary $r_{\mathrm{LIS}}$ is set at over 100 AU, far from the star to avoid stellar wind influence. In this study, the heliosphere radius is used as $r_{\mathrm{LIS}}$ [26]:

$$r_{\mathrm{LIS}}(\Omega )=R_{\mathrm{h}}(\Omega )=R_{\mathrm{h}}\left({\Omega }_{\odot }\right)\sqrt{\frac{\dot{M}(\Omega )v(\Omega )}{\dot{M}\left({\Omega }_{\odot }\right)v\left({\Omega }_{\odot }\right)}}$$

(7)

where $R_{\mathrm{h}}\left({\Omega }_{\odot }\right)=122$ AU is the current heliosphere radius; v is the radial velocity of the stellar wind at 1 AU; and $\dot{M}$ is the mass loss rate of the star, with the calculation formula provided in [27].

Measurements by Voyager 1 and 2 provide important observational data for the LIS, which can be fitted by the following differential intensity spectrum [31]:

$$j_{\mathrm{LIS}}\left(T\right)=2.70\frac{T^{1.12}}{{\beta }^2}{\left(\frac{T+0.67}{1.67}\right)}^{-3.93}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}{\mathrm{sr}}^{-1}{\mathrm{MeV}}^{-1},$$

(8)

where j is the differential intensity of CRs, representing the number of particles passing through a unit area, per unit time, per unit solid angle, and per unit energy interval; T is the kinetic energy of CRs in GeV; and $\beta =v/c$ is the ratio of particle speed to the speed of light. The differential intensity $j\left(T\right)=p^{2}f\left(p\right)$, where $f\left(p\right)$ is the phase space distribution function of CRs.

2.2. Numerical Methods

To numerically solve the Parker transport equation, it is necessary to discretize the three independent variables: time (t), space (r), and momentum (p).

Explicit difference schemes for diffusion have to take a sufficiently small timestep $\Delta t$ as required by the Courant–Friedrichs–Lewy condition, thus that it is time-consuming generally. Hence, the semi-implicit Crank-Nicolson scheme [32] for diffusion is adopted (see Appendix A for further details), while the explicit upwind schemes for advection and adiabatic loss are used. In addition, the operator splitting method is adopted [33].

To further improve computational efficiency, we utilize the Numba1[34] library to compile critical calculation steps in Python code into efficient machine code, significantly reducing runtime.

2.3. Results

By solving the Parker equation with the aforementioned algorithm and parameter settings, we obtain the temporal and spatial evolution of LECRs. Figure 1 presents the momentum spectra of LECRs at different distances for the stellar rotation rate $\Omega =1.0{\Omega }_{\odot }$ (${\Omega }_{\odot }$ is the current solar rotation rate) and $\Omega =2.1{\Omega }_{\odot }$ (it is estimated that the Sun’s age was $t=1\mathrm{Gyr}$ when life began on Earth, corresponding to $\Omega \sim 2.1{\Omega }_{\odot }$).

The figure shows a significant decrease in the flux of LECRs as the distance r decreases, especially in the low energy region where the flux is significantly attenuated in the inner interplanetary space. Additionally, the peak flux shifts to higher energies as r decreases, consistent with the expected effects of diffusion and stellar magnetic field attenuation. The spectrum in the inner interplanetary space shows a steep drop at low energies, while the spectrum in the outer interplanetary space remains relatively flat. This matches our estimates of the protoplanetary disk morphology, where the inner interplanetary space exhibits significant modulation and shielding effects, while the outer interplanetary space is less affected by stellar wind.

Furthermore, compared to the case with $\Omega =1.0{\Omega }_{\odot }$, the LECR flux at the same location with $\Omega =2.1{\Omega }_{\odot }$ is smaller. This is understandable—young stars have stronger stellar winds, which inhibit the entry of LECRs.

The shape of the calculated spectrum curves agrees with results from previous related studies [26], validating the reliability of the Parker transport model and the numerical algorithm used in this study. In conclusion, our method can accurately simulate the transport and diffusion behavior of LECRs in the environment of a protoplanetary disk.

3. Calculation of Nuclear De-Excitation Line Emission

LECR nuclei interacting with protoplanetary disk material can trigger a series of nuclear reactions, exciting atomic nuclei to high-energy states, which subsequently de-excite by emitting gamma rays in the 0.1–10 MeV range. These gamma rays mainly originate from two processes: collisions between protons and helium nuclei in LECRs with heavy elements in the environmental gas, and collisions between heavy nuclei in LECRs with hydrogen and helium atoms in the environment [14]. Narrow lines are produced when CR protons and $\alpha$ particles interact with ambient nuclei heavier than helium (such as C, N, O). Broad lines occur when heavy nuclei in CRs (such as C, N, O) are excited by the medium’s hydrogen and helium, subsequently de-exciting to produce gamma rays with broader line widths (fractional width of about 20%).

In this section, we compute the gamma-ray radiation in the protoplanetary disk based on the LECR energy spectrum distribution obtained in section 2, combined with models of stellar atmospheres and CR elemental abundances, laying the necessary foundation for subsequent discussions.

3.1. Calculation Method

Following the method developed by Ramaty [14], Murphy [20], Benhabiles-Mezhoud [21], we applied the same procedure as described in section 3.1 of Liu [22]. For the nuclear reaction cross sections, we use TALYS code (1.96 version) [38,39] for theoretical simulations as a complement of experimental data complied by Murphy [20]. Additionally, the relative abundances of different elements in the protoplanetary disk significantly affect the shape of the gamma-ray radiation spectrum. This study uses stellar atmospheric elemental abundances and CR elemental abundances models from [14], as shown in Table 1. By multiplying the contribution of each element by its corresponding abundance and summing the weighted contributions of all elements, the total gamma-ray emissivity spectra per hydrogen atom are obtained. Finally, by integrating (in practice, performing a discrete sum in the code) the gamma-ray emissivity throughout the entire protoplanetary disk (from $r=r_{\odot }\sim 0$ to $r=r_{\mathrm{LIS}}$), and assuming a hydrogen atom number density of 1 cm³ in the disk, we can calculate the observed gamma-ray flux assuming the star is 1 kpc away from the Earth.

3.2. Results

Figure 2 shows the calculated gamma-ray emissivity at different distances to the center star. Because LECRs are significantly attenuated by the stellar wind in the inner planetary space, the corresponding gamma-ray emissivity decreases as the distance r decreases.

The spectra display some distinct lines, such as the 4.44 MeV line emission mainly resulting from the de-excitation of ¹²C and the 6.13 MeV line emission mainly produced via the de-excitation of ¹⁶O. These spectral lines indicate significant nuclear reactions between LECRs and the protoplanetary disk material. The gamma-ray emissivity for the case with $\Omega =2.1{\Omega }_{\odot }$ is lower compared to $\Omega =1.0{\Omega }_{\odot }$ at the same location. This difference reflects the lower CR flux at that time.

Figure 3 presents the theoretical calculation of the gamma-ray flux spectra supposing the observer is located at a distance of 1 kpc from the star. It is evident that some of the original spectral lines are preserved, notably the 4.44 MeV de-excitation line of ¹²C and the 6.13 MeV de-excitation line of ¹⁶O.

These results highlight how the interaction of LECRs with different elements produces a characteristic gamma-ray emission spectrum, consisting of both narrow and broad lines, as well as a continuous radiation component. Therefore, the structure and intensity of the gamma-ray emission spectrum provide valuable information about the LECRs themselves. Through detailed analysis of the gamma-ray spectrum, the composition, energy, and spatial distribution of LECRs at that time can be coherently inferred, which is crucial for understanding the radiation and ionization environment of early planetary systems.

4. Discussion

In this article, we calculated the spectra of LECRs in the protoplanetary disk in Section 2 by using Parker’s Equation (2) under the assumption of a spherically symmetric wind, and the resulting nuclear de-excitation gamma-ray line emission spectra in Section 3 of LECRs interacting with the protoplanetary disk gases assuming a uniform distribution of material in the disk. Nonetheless, the protoplanetary disk is a much more complex system than that considered in the present study, and it can be very important to consider a more realistic model for making comparisons with future observations. Hence, we discuss several important points that can further improve our present model below.

It is important to consider the physical and chemical structure and material distribution of the protoplanetary disk. This study adopts a simplified spherical symmetric model, assuming a uniform distribution of material in the protoplanetary disk. The assumption of azimuthal symmetry means that any short-term modulation effects shorter than the stellar rotation period are ignored, but the actual protoplanetary disk may have a flattened disk structure with significantly uneven material distribution perpendicular to the disk plane. The geometric structure will affect the propagation of stellar wind and LECRs in the protoplanetary disk. For example, due to the higher material density in the disk plane, CRs are more likely to undergo scattering and energy loss in this region. Observations also reveal many structural details, such as planetary rings, gaps, and spiral arms [40]. Therefore, it is necessary to establish a three-dimensional protoplanetary disk model that is closer to the actual situation and study the effects of geometric and material distribution through numerical simulations, such that we can infer more information about the physical and chemical conditions of the protoplanetary disks from gamma-ray observations.

Magnetic fields play a significant role in the dynamics and evolution of protoplanetary disks and stellar winds [41]. Therefore, it is crucial to study the effects of magnetic fields on the transport of LECRs in the protoplanetary disks [6,42]. The current model assumes a uniform diffusion coefficient, ignoring the impact of magnetic field evolution. To address this, it is necessary to coherently simulate the dynamic evolution of the magnetic field in the protoplanetary disk environment and couple it with the calculation of LECR transport. This requires further development of relevant theoretical models and numerical algorithms.

The turbulence of protoplanetary disks [43,44] is highly important for the transport of LECRs. The parameter ${\eta }_0$ describes the turbulence level in the stellar wind. A higher value means that CRs propagate further before scattering. For example, young stars are believed to produce more coronal mass ejections because they are observed to have higher flare rates, which may lead to stronger turbulent components in the magnetic fields of planetary systems. This study uses ${\eta }_0=1$, representing the turbulence limit of CRs in the stellar wind, which may not accurately approximate the situation for stars. Assuming a constant ${\eta }_0$ means that the diffusion coefficient decreases with increasing magnetic field strength and therefore also decreases with increasing stellar rotation rate. The value of $\gamma$ is related to the power of turbulence. This study uses $\gamma =1$ because it matches observations of solar wind modulation; however, this approximation may not be accurate for other stars.

Currently, the LECR transport model and the gamma-ray radiation calculations are conducted separately. To obtain coherent results, the two need to be tightly integrated. This involves simultaneously calculating the LECR energy distribution at each location while solving the transport equation and then calculating the corresponding gamma-ray radiation spectrum. By iteratively solving these, we can ultimately obtain the coherent distribution of LECRs and gamma rays throughout space and time. This requires developing new numerical algorithms to improve computational efficiency.

5. Summary and Outlook

This study proposes a novel method for investigating LECRs in protoplanetary disks via gamma-ray spectra. By solving the Parker transport equation, we simulated the propagation of LECRs in the protoplanetary disk and calculated the energy spectra of LECRs at different distances. Additionally, we simulated the nuclear reactions between LECRs and the gas and dust in the protoplanetary disk, producing nuclear de-excitation lines, and calculated the gamma-ray radiation spectra. This established a connection between the properties of LECRs and observable gamma-ray spectra. This method enables us to infer the composition, energy distribution, and spatial distribution of LECRs during the early stages of planetary system formation, thereby reconstructing the high-energy particle radiation environment of that period.

LECRs play an important role in various stages of planetary system formation. They can ionize and heat the gas in the protoplanetary disk, affecting chemical reactions and dust evolution, and thus impacting the formation and evolution of planetary systems. Therefore, accurately detecting and understanding the properties of LECRs is vital for understanding the mechanisms of planetary system formation. Traditional direct measurements face challenges, such as modulation by stellar wind, making our proposed method a valuable alternative for in-depth LECR research.

The future of LECR research is promising with the emergence of new-generation MeV detector projects such as MeVGRO 2 [45], MASS [46], e-ASTROGAM [47], AMEGO 3 [48,49], COSI 4 [50], and MeVCube [51]. With the deployment of these new high-resolution MeV gamma-ray telescopes, we are entering a golden age for LECR research. We note that these planned instruments mentioned above can hardly archive an angular resolution better than $1^{\circ }$ in the MeV energy range, making the observation of protoplanetary systems difficult. Nevertheless, new techniques such as the Laue lens [52] can, in principle, provide an angular resolution of dozens of arcseconds, making the observation of the protoplanetary system in the MeV band feasible.

Through precise gamma-ray observations of protoplanetary disks, combined with advanced theoretical models, we will be able to reconstruct the real conditions of LECRs during the birth of planetary systems by obtaining information about the composition, energy spectrum, and spatial distribution of LECRs, thus understanding the high-energy particle environment of early planetary systems. We will also reveal the impact of CRs on planetary system formation by studying the effects of LECRs on chemical evolution, dust growth, and planet formation in protoplanetary disks, thereby understanding the origins of planetary system diversity. Furthermore, we will explore the CR environment of life’s origin [53,54,55] by assessing the impact of LECRs on the radiation environment and the origin of life on early Earth, and searching for habitable zones and clues to extraterrestrial life.

In conclusion, utilizing gamma-ray detection of LECRs not only provides a unique approach to investigating the cosmic ray environment during the early stages of planetary system formation but also holds significant implications for deepening our understanding of the formation processes of stellar and planetary systems. This research will advance our comprehension of the universe, life, and the origins of humanity to new heights.