Anisotropic diffusion of high-energy cosmic rays in magnetohydrodynamic turbulence

Cosmic rays (CRs) are charged particles moving with relativistic speeds, which can be directly detected in and near the solar system and indirectly observed from the emissions they produce far from the solar system (such as in the Milky Way and another galactic disk; Armillotta et al. 2022). Exploring the propagation processes of the Galactic CRs, including advection, diffusion, scattering, acceleration/reacceleration, and energy loss processes, will help to understand their origins (Castellina & Donato 2005; Strong et al. 2007; Blasi 2013; Grenier et al. 2015). Since propagating CRs will inevitably interact with turbulent magnetic fields, it also provides a useful tool for probing the properties of the interstellar medium (ISM).

Magnetohydrodynamic (MHD) turbulence is ubiquitous in astrophysical plasmas ranging from interplanetary space to interstellar and intergalactic media. It has been pointed out that MHD turbulence directly or indirectly affects the propagation of CRs. Studying the propagation of CRs in turbulent magnetic fields can help us to comprehend many important issues both in space and astrophysics, e.g., solar modulation of CRs (Jokipii & Parker 1969), driving Galactic winds (Wiener et al. 2017; Krumholz et al. 2020), CR anisotropy (Qiao et al. 2023; Li et al. 2024), X-rays binaries (Zhang et al. 2018), diffuse $\gamma-$ray emission (Yan et al. 2012), as well as the confinement and reacceleration of CRs (Chandran 2000; Yan & Lazarian 2002; Zhang & Xiang 2021; Gao & Zhang 2024).

The High-Altitude Water Cherenkov Observatory (HAWC) collaboration first reported the TeV $\gamma-$ray halos around the pulsars Geminga and Monogem (PSR B0656+14), with a spatial extension of $\sim$ 10-50 pc (Abeysekara et al. 2017). They explained these observations by the inverse Compton (IC) scattering process, that is, the high-energy CR electrons and positrons are injected from the pulsar wind nebula (PWN) and diffuse into the ISM to scatter off cosmic microwave background photons. Abeysekara et al. (2017) modeled the surface brightness profile of the $\gamma-$ray emission of these two TeV halos to constrain a low diffusion coefficient of $D=(4.5\pm 1.2)\times 10^{27}\ {\rm cm}^{2}\ \rm s^{-1}$ at 100 TeV. This diffusion coefficient is hundreds of times lower than the average value in the Milky Way as inferred from the secondary-to-primary ratio measurements in the CR spectrum (i.e., the boron-to-carbon ratio (B/C); Aguilar et al. 2016). This suggests that the diffusion coefficient may be highly inhomogeneous on a small scale. It could be meaningful to investigate the origin of this slow-diffusion region around the pulsars, which can help us to understand many important issues, such as the CR position excess (Hooper et al. 2017; Fang et al. 2018, 2019b; Tang & Piran 2019; Xi et al. 2019; Manconi et al. 2020), diffusion TeV $\gamma-$ray excess (Linden & Buckman 2018; Yan et al. 2024), and particle propagation near CR source (Recchia et al. 2021; López-Coto et al. 2022).

At present, there are mainly three plausible explanations for this slow diffusion as follows. Firstly, this slow diffusion is explained by the projection effects under the framework of anisotropic diffusion, which requires the nearby turbulent environment to be sub-Alfvénic and a small viewing angle between the local mean magnetic field and the line of sight (LOS) (Liu et al. 2019; De La Torre Luque et al. 2022; Fang et al. 2023). Secondly, it may be caused by the turbulent environment generated by the parent supernova remnants of the pulsars (Fang et al. 2019a) or the escaped electrons themselves (Evoli et al. 2018; Mukhopadhyay & Linden 2022). Alternatively, the mirror diffusion¹¹1Lazarian & Xu (2021) found that in MHD turbulence CRs would not be trapped, due to the perpendicular super-diffusion of turbulent magnetic fields (Xu & Yan 2013; Lazarian & Yan 2014; Hu et al. 2022). On the contrary, they bounce between different magnetic mirrors and move along the local magnetic field, which leads to a new diffusion mechanism called mirror diffusion. With the stochastic change of pitch angles due to gyroresonant scattering, CRs undergo slow mirror diffusion at relatively large pitch angles and fast scattering diffusion at smaller pitch angles, resulting in a Levy-flight-like propagation that is slower than that induced by scattering alone (Xu 2021; Zhang & Xu 2023; Barreto-Mota et al. 2024; Xiao et al. 2024). As estimated in Barreto-Mota et al. (2024), the physical value of the diffusion coefficient is $D_{\parallel}\leq 10^{27}\ {\rm cm^{2}\ s^{-1}}$ for 100 TeV CRs within an extended zone with tens of pc, which is consistent with the expected value to explain the observations of extended sources such as Geminga and Monogem (Abeysekara et al. 2017). may be a channel for understanding this slow diffusion process because the combination of resonance scattering with mirror diffusion can cause a slower diffusion than that induced by scattering alone (Lazarian & Xu 2021; Xu 2021; Zhang & Xu 2023; Barreto-Mota et al. 2024).

Due to the presence of ubiquitous magnetic turbulence, the resulting CR anisotropic diffusion may be reasonable in the key astrophysical environments (Holod et al. 2005; Parrish & Stone 2005; Parrish et al. 2008) and has obtained theoretical support (Evoli et al. 2012; Cerri et al. 2021). Recently, to explain the slow diffusion coefficient, the size, and morphology of observed TeV halos, Liu et al. (2019) and De La Torre Luque et al. (2022) limited the Alfvén Mach number to a small value of $M_{\rm A}\lesssim 0.3$, and the viewing angle to $\psi\lesssim 5^{\circ}$. To some extent, these two small phase spaces may impede the application of the anisotropic diffusion model.

The main motivations for the current work are the need for slow diffusion to explain the observations mentioned above and for the narrow parameter space constraints in the anisotropic diffusion model. Therefore, in the framework of CR anisotropic diffusion, we explore the interactions of high-energy electrons with the ambient turbulent magnetic fields, considering PWN, i.e., the Crab-like nebula, as a general research case. We want to know whether the slow-diffusion phenomenon could happen in a reasonable parameter space. The relevant CR diffusion law related to MHD turbulence theory to what extent is applicable for understanding the transport of CR particles.

This paper is organized as follows. In Sect. 2, we briefly introduce the theoretical description of diffusion properties of CRs in MHD turbulence, including the energy-dependent parallel diffusion and anisotropic diffusion model. We perform a parameterized simulation and describe our initial simulation setup in Sect. 3. We present the numerical results in Sect. 4. Discussion and summary are provided in Sects. 5 and 6, respectively.

The propagation process of CRs can be described with the Fokker-Planck equation (Gaisser 1990; Berezinskii et al. 1990; Schlickeiser 2002), which is defined by

where $f(\bm{r},E,t)$ is the CR distribution function in the position space $\bm{r}=(x,\ y,\ z)$ and the kinetic energy space $E$. $r_{\parallel}=x$ is the coordinate along the mean magnetic field, and $r_{\perp}=\sqrt{y^{2}+z^{2}}$ is the coordinate perpendicular to the mean magnetic field. $D_{\parallel}$ and $D_{\perp}$ are the spatial diffusion coefficients parallel to the mean magnetic field and perpendicular to it, respectively. $D_{\rm EE}=E^{2}v_{\rm A}^{2}/9D_{\parallel}$ is the energy diffusion coefficient (Michałek & Ostrowsky 1996), and $\dot{E}$ is the rate of energy loss. In the case of power-law injection, the source function is defined by $S=\frac{d\dot{n}}{dE}\delta(\bm{r})$, where $d\dot{n}/dE\propto E^{-q}$ with a power-law index $q$ over the range from $E_{0}$ to $E_{1}$. While for mono-energy injection, the source function is defined by $S=\frac{\mathcal{L}}{E}\delta(\bm{r})\delta(E-E_{\rm in})$, where $\mathcal{L}=10^{38}\ {\rm erg\ s^{-1}}$, $E_{\rm in}$, and $m_{\rm e}$ are the source luminosity, initial energy, and mass of the injected electrons, respectively.

To explore the interactions of high-energy particles with ambient magnetic turbulence, we numerically solve Eq. (1) using the CR transport code CRIPTIC (Krumholz et al. 2022). With a PWN environment, we consider a uniform region containing gas with a density of $\rho\simeq 10^{-24}\ {\rm g\ cm^{-3}}$ and velocity of $v\simeq$ 300 $\rm km\ s^{-1}$, threaded by a uniform magnetic field $B_{0}$ in the $x$-axis direction. For the general environment of PWN, we have the typical temperature of thermal gas of $T\simeq 10^{5}$ K (Zyuzin et al. 2021) and the sonic speed of $c_{\rm s}\simeq$ 200 $\rm km\ s^{-1}$, resulting in a slightly supersonic Mach number of $M_{\rm s}=v/c_{\rm s}\simeq 1.5$. When involving radiative loss processes, we consider synchrotron radiative losses of relativistic electrons due to the presence of the ISM magnetic field and IC scattering losses due to both interstellar radiation field with a temperature of 20 K and cosmic microwave one in the Klein-Nishina regime.

With the purpose of exploring how the magnetization parameter $M_{\rm A}$ and turbulence spectral index $\gamma$ affect particles’ diffusion behavior²²2In practice, we control the value of $M_{\rm A}$ by regulating the turbulent magnetic field strength $B_{0}$. For example, setting $B_{0}=210\ \rm\mu G$, one has $M_{\rm A}=0.51$. , together with various radiative processes, we divide our simulations into three groups:

• •

  ${\rm Group}\ \mathcal{A}$: the influence of radiative losses on CR diffusion. In this case, with the fixed $\gamma=$ 5/3 and $M_{\rm A}=0.51$, we change the energy range and distribution of CR electrons. We consider two distribution types of CR electrons: a power-law energy distribution of $d\dot{n}/dE\propto E^{-q}$ with $q=$ 2.2 over the range from $E_{0}$ to $E_{1}$, and another mono-energy distribution.

• •

  ${\rm Group}\ \mathcal{B}$: the influence of turbulence spectral index on CR diffusion. With the fixed $M_{\rm A}=0.51$ and a series of mono-energy distributions, we change the spectral index of magnetic turbulence, corresponding to the different turbulence types.

• •

  ${\rm Group}\ \mathcal{C}$: the influence of the magnetization parameter on CR diffusion. After fixing $\gamma=$ 5/3, we mainly consider the changes of the magnetization parameter $M_{\rm A}=v/v_{\rm A}=v/(B_{0}/\sqrt{4\pi\rho})$, which will match with different magnetic field strength.

In the framework of the three groups described above, we run 90 sets of simulations up to $t=3\ {\rm kyr}$ with an output step of about 10 yr, considering the limited saving abilities of the data. Given a constant injection with the particles injection rate $\Gamma=1.6\times 10^{2}\ \rm yr^{-1}$, the simulations include amounts of $4.8\times 10^{5}$ particles. For more details, the interested reader can refer to Appendix A for numerical procedures.

Through numerically solving the CR transport equation in a Crab-like nebula environment, we explored the interaction of high-energy CR electrons with magnetic turbulence in the framework of an anisotropic diffusion model. Specifically, regarding the magnetization parameter and turbulence spectral index as the adjustable parameters, we explore their influence on CR spectral, spatial diffusion, and anisotropy morphology.

The propagation process of CRs can be described with the time-dependent Fokker-Planck equation (Gaisser 1990; Berezinskii et al. 1990; Schlickeiser 2002). Provided that one considers different simplifications, the transport equation can be solved analytically (Webber et al. 1992; Maurin et al. 2002; Shibata et al. 2004). However, this simplification process may lead to the loss of CR propagation information. Therefore, the numerical/simulation approach is the main channel for revealing the CR propagation. At present, there are a number of numerical codes developed, e.g., GALPROP (Strong & Moskalenko 1998; Moskalenko & Strong 1998), PICARD (Kissmann 2014), DRAGON (Evoli et al. 2008), CRPROPA (Merten et al. 2017), and CRIPTIC (Krumholz et al. 2022), to understand CR observations. In the current work, we explored the diffusion properties of CRs using the CRIPTIC code, which is a new code for calculating both CR transport and observable emission.

In the framework of anisotropic diffusion, we investigated the influence of the magnetization parameter $M_{\rm A}$ and turbulence index $\gamma$ on the CR spatial distributions. In the sub-Alfvénic regime, we found that the resulting morphology and size of CR spatial distribution are consistent with that of emission images (Liu et al. 2019; De La Torre Luque et al. 2022) within small viewing angles, while in the case of super-Alfvénic turbulence, similar spatial distributions are also formed in the direction perpendicular to the mean magnetic field. Note that the stronger the anisotropy, the larger $M_{\rm A}$ is, and the smaller the spatial extent, the larger $\gamma$. Therefore, when changing the viewing angle with both the large $M_{\rm A}$ and $\gamma$, the expected spatial structures can be reproduced. For instance, with $D_{\perp}/D_{\parallel}\sim M_{\rm A}^{3}$ in the super-Alfvénic turbulence⁶⁶6This work mainly focuses on the strong turbulence region, that is, the spatial scale is smaller than the transition scale $l_{\rm A}$. For instance, considering the magnetic field strength of 50 $\mu$G (corresponding to the super-Alfvénic regime) and the maximal energy of CRs of 10 PeV, we can simply estimate the maximal gyroradius of CRs with a velocity of $v_{\rm e}$ as (Barreto-Mota et al. 2024) $R_{\rm g,max}\approx\frac{(E/1\ {\rm PeV})}{(B/1\ \mu{G})}\frac{v_{\rm e}}{c}\approx\frac{(10\ {\rm PeV}/1\ {\rm PeV})}{(50\ \mu{\rm G}/1\ \mu{\rm G})}\approx 0.2\ {\rm pc}$, which is much smaller than the parallel diffusion scales $l_{\perp}\approx$ hundreds of pc, as seen in Fig. 4(c) for one super-Alfvénic case. Given such a small gyroradius, CR diffusion is strongly subjected to the anisotropy of magnetic fields., when $M_{\rm A}\gtrsim 4.7$, the parallel diffusion coefficient is about hundreds of times lower than the perpendicular one. This provides an alternative way for explaining the slow diffusion and morphology of observed TeV halos, with the LOS approximately aligned with the mean magnetic field.

For the other plausible explanation of the slow diffusion, i.e., the mirror diffusion, Lazarian & Xu (2021) analytically predicted the mirror diffusion is generally slower than the diffusion of non-bouncing CRs with small pitch angles between the local mean magnetic field and velocity of particles that undergo gyroresonant scattering. Subsequently, many numerical studies focused on exploring the diffusion and acceleration properties of CRs in the frame of mirror diffusion (e.g., Xu 2021; Xu & Li 2023; Zhang & Xu 2023, 2024; Barreto-Mota et al. 2024). They found that the mirror diffusion of CRs with sufficiently large pitch angles can explain the slow diffusion around the CR source. CRs may stochastically undergo slow mirror and fast scattering diffusion in the Galactic diffuse medium away from CR sources. We also expect that the mirror diffusion may be a plausible explanation for the slow diffusion near the source, which needs further research.

In the current work, in addition to the continuous injection used throughout this paper, we also tested the case of instantaneous, finding that CR’s overall evolution is similar to the continuous one. The only difference is the lack of the power-law relationship between the energy density and the effective radius. We also tested the streaming interaction for a wide distribution of CR energy (up to 10 PeV) and found that the streaming interaction almost does not change our results. We note that the streaming affects only the diffusion of the GeV CRs (see also the review Cesarsky (1980) and Farmer & Goldreich (2004) for the theoretical expectation). Therefore, the results we presented in the current work do not consider the streaming interaction. Since our study focuses on the CR electrons, it is enough to take the cooling effects of synchrotron radiation and IC scattering into account. When studying diffusion processes of the CR protons, one has to consider more cooling effects such as the coulomb, ionization, bremsstrahlung, and positron annihilation. Considering the PWN environment as a general research object, this work mainly focused on understanding how CR electrons interact with turbulent magnetic fields. Studies explaining real observations will be discussed elsewhere.

Considering PWN as a general research case in this paper, we numerically study the interactions of CRs with the ambient turbulent magnetic field. Focusing on the CR spectral distribution, anisotropy morphology, and spatial diffusion coefficients, our work is devoted to understanding CR propagation processes in the presence of magnetic turbulence. The main findings are briefly summarized as follows:

1. 1.

   We find that the energy density $U_{\rm CR}$ without radiative losses displays the power-law relationships of $U_{\rm CR}\propto r^{-6/5}$ and $U_{\rm CR}\propto 1/E$, which are modified as $U_{\rm CR}\propto r^{-11/5}$ and $U_{\rm CR}\propto E^{-8/5}$ in the presence of radiative losses, respectively.

2. 2.

   The morphology of the CR spatial distribution strongly depends on the properties of turbulence. Namely, the magnetization parameter $M_{\rm A}$ affects the spatial anisotropic distribution, while the turbulence spectral index $\gamma$ its spatial extent. Note that the viewing angle also affects the spatial anisotropic distribution, consistent with the previous studies (e.g., Liu et al. 2019; De La Torre Luque et al. 2022).

3. 3.

   The CR electrons suffer from a slow diffusion near the source and a fast/normal diffusion away from the source. The effect of radiative losses, which happen far from the source, can suppress CR diffusion processes, resulting in a sub-diffusion behavior.

4. 4.

   The averaged energy density $U_{\rm CR}$ distributes a power-law relationship of $U_{\rm CR}\propto M_{\rm A}^{-6.0\pm 0.5}$, independent of both CR kinetic energy and radiative losses.

5. 5.

   Our simulation demonstrates that the relations of $D_{\parallel}\sim E^{2-\gamma}$ and $D_{\perp}/D_{\parallel}\sim M_{\rm A}^{\alpha}$ maintain their features throughout a long evolution period. Since these two relations establish a robust connection between turbulent magnetic fields and CR diffusion coefficients, one can apply them to understand complex astrophysical processes related to turbulence cascades.