Two-zone Diffusion of Electrons and Positrons from Geminga Explains the Positron Anomaly

The cosmic-ray positron excess has been known for nearly a decade (Adriani et al. 2009; Ackermann et al. 2012; Aguilar et al. 2013), but its origin is still a mystery. These extra positrons may originate either from astrophysical sources like pulsars or the dark matter annihilation/decay. Among various kinds of astrophysical sources, the Geminga pulsar (PSR J0633+1746) has been widely believed to be a very promising candidate to produce the positron excess (Hooper et al. 2009, 2017; Yüksel et al. 2009; Yin et al. 2013). Geminga is one of the nearest pulsars, with a distance of ${250}_{-62}^{+120}$ pc (Faherty et al. 2007). Its age is estimated to be about 3.42 × 10⁵ yr, and the derived spin-down energy is 1.23 × 10⁴⁹ erg (Manchester et al. 2005; Fang et al. 2018). All these parameters suggest that Geminga can probably dominate the high-energy positron flux observed on the Earth. The extended TeV γ-ray halo around the Geminga pulsar observed by Milagro (Abdo et al. 2007) and HAWC (Abeysekara et al. 2017b) gives straightforward evidence supporting that Geminga can indeed generate very-high-energy (VHE) electrons and positrons (e^±).

However, the detailed morphological study of the VHE γ-ray emission from Geminga and PSR B0656+14 (Monogem) by HAWC suggests that the e^± produced by these pulsars diffuses out significantly slower than those in the average interstellar medium (ISM) as inferred from the boron-to-carbon ratio (B/C) measurements (Abeysekara et al. 2017a). In such a case, the e^± produced by Geminga or Monogem can hardly reach the Earth, thus these two pulsars may be unlikely to account for the positron excess. On the other hand, Hooper & Linden (2017) pointed out that this slow-diffusion scenario should not be representative even in the local environment. Since H.E.S.S. has detected high-energy e^± up to ∼20 TeV,⁵ these e^± can only travel for 10 ∼ 20 pc within the cooling time given such a slow-diffusion condition. We can hardly find any high-energy e^± sources within such a small distance around the solar system.

Since the HAWC data can only probe a region of ∼30 pc around Geminga and Monogem, it is very likely that the slow-diffusion region is actually limited to a small region around the sources, beyond which particles diffuse faster as typical Galactic cosmic rays. This scenario would be consistent with the B/C data and the H.E.S.S e^± spectrum (Hooper et al. 2017; Hooper & Linden 2017). In this work, we investigate whether this two-zone diffusion model can explain both the HAWC γ-ray data of Geminga and the positron excess. The two-zone diffusion model is solved with a numerical method. Because the diffusion coefficient has a jump at the boundary between the two zones, the differencing scheme should be carefully dealt with.

In Section 2, we introduce the two-zone diffusion model and the numerical treatment to the propagation equation. In Section 3, we calculate the positron spectrum of Geminga in the two-zone scenario, and compare the result with the AMS-02 data (Aguilar et al. 2014). Then we conclude our work with some discussion in the last section.

The propagation process of e^± can be described by the diffusion-cooling equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial N}{\partial t}-{\rm{\nabla }}(D{\rm{\nabla }}N)-\displaystyle \frac{\partial }{\partial E}({bN})=Q,\end{eqnarray} \tag{ 1 }$$

where N is the differential number density of e^±, D denotes the diffusion coefficient, b is the energy-loss rate, and Q is the source term. In the present work, we are interested in the energy range higher than 10 GeV, so the convection and reacceleration terms that affect the low-energy spectrum are neglected (Delahaye et al. 2009). The energy-loss rate has the form of b(E) = b₀(E)E², which describes the synchrotron and inverse Compton radiation cooling of e^±. The interstellar magnetic field in the Galaxy is set to be 3 μG to get the synchrotron term (Minter & Spangler 1996). For the inverse Compton scattering term, it is necessary to consider a relativistic correction to the scattering cross section. We follow the calculation of Schlickeiser & Ruppel (2010), where b₀ is energy-dependent.

The diffusion coefficient is usually assumed to be D(E) = β^ηD₀(R/1 GV)^δ, where D₀ and δ are both constants, β is the velocity of particles in units of light speed, η is a low-energy correction parameter of the velocity dependence, and R is the rigidity. Assuming D is spatially constant, Yuan et al. (2017) constrained the propagation parameters with the B/C data of AMS-02 (Aguilar et al. 2016), and found that the best-fit model is the diffusion model with reacceleration. The obtained propagation parameters are D₀ = (2.08 ± 0.28) × 10²⁸ cm² s⁻¹ and δ = 0.500 ± 0.012 (hereafter Y17 model). The thickness of the propagation halo is z_h = 5.02 ± 0.86 kpc, and the low-energy correction parameter η is not relevant for this study.

For the e^± produced by Geminga, the propagation distance can be estimated as $2\sqrt{D(E){t}_{E}}$, where ${t}_{E}=\min \{{t}_{{\rm{g}}},1/({b}_{0}E)\}$, i.e., the smaller one of the age of Geminga t_g and the cooling time (Abeysekara et al. 2017a). If we adopt the propagation parameters of Y17, e^± of ≃1 TeV energies can diffuse to a distance of ≃1.7 kpc. This scale is much smaller than the z_h obtained in Y17. If the particles diffuse slower in the region around Geminga, the total diffusion distance is even less. Therefore, it should be fine to assume a spherically symmetrical geometry of the propagation of e^± from Geminga. The diffusion coefficient for the two-zone model is then

$$\begin{eqnarray}D(E,r)=\left\{\begin{array}{ll}{D}_{1}(E), & r\lt {r}_{\star }\\ {D}_{2}(E), & r\geqslant {r}_{\star }\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where r is the distance from Geminga, r_⋆ is the discontinuity shell of the diffusion coefficient, D₁ is the diffusion coefficient around Geminga, which is inferred by the HAWC γ-ray data (Abeysekara et al. 2017a), and D₂ is the average diffusion coefficient in the Milky Way given by Y17. In a more realistic picture the diffusion coefficient may change gradually. However, no observations can constrain the transitional zone. Furthermore, since r_⋆ is a free parameter we can always get a equivalent result for the local e^± flux by adjusting r_⋆ in this simplified picture.

The spatially dependent diffusion equation is difficult to solve analytically, so we adopt a numerical method instead. We assume that Geminga is a burst-like point source, then Equation (1) can be rewritten as follows, along with the initial condition and the boundary conditions:

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{\partial N}{\partial t}={ \mathcal L }N,\\ N(0,E,r)=Q(E)\delta (r),\\ N(t,{E}_{\max },r)=0,\\ {\left.\tfrac{\partial N}{\partial r}\right|}_{r=0}=0,\\ N(t,E,{r}_{\max })=0,\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where E_max = 500 TeV, and r_max = 4 kpc. The operator ${ \mathcal L }$ is the sum of the diffusion operator ${{ \mathcal L }}_{r}$ and the energy-loss operator ${{ \mathcal L }}_{E}$, which are

$$\begin{eqnarray}{{ \mathcal L }}_{r} & = & \displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}D(r)\displaystyle \frac{\partial }{\partial r}\right],\\ {{ \mathcal L }}_{E} & = & b\displaystyle \frac{\partial }{\partial E}+\displaystyle \frac{\partial b}{\partial E}.\end{eqnarray} \tag{ 4 }$$

We apply the operator splitting method to deal with the two operators separately. For the discretization of the energy-loss scheme, we follow the method given by Kissmann (2014), while for the diffusion operator, we use two different discretizations depending on if r_i lies on the discontinuity surface (i is the spatial step index). If not, the well-known Crank–Nicolson scheme for constant D is a good choice. For the case of r_i = r_⋆, the difference equation should be rewritten to ensure the conservation of the flux at both sides of the discontinuity surface. We adopt the finite volume method to derive the differencing scheme for this case. All the details of the discretization are presented in the Appendix. Note that the differencing scheme derived by the finite volume method (Equation (7)) is different from those used in GALPROP (Strong & Moskalenko 1998) or DRAGON (Evoli et al. 2008); the numerical schemes of GALPROP and DRAGON are basically applicable for continuously and slowly changing diffusion coefficients.

Then we describe the parameter settings. In the difference equation, the time-step Δt is set to be 1000 years, which is much smaller than the age of Geminga. The spatial step Δr is 1 pc, then the initial condition given in Equation (17) describes a spherical source with a radius of 1 pc. For an old source like Geminga, the spectrum at the Earth can be very close to that of a point source, especially when E ≳ 100 GeV (Thoudam & Hörandel 2012). We use a logarithmic scale with a ratio of 1.2 for the energy grids, that is, E_l+1/E_l = 1.2, where l is the energy step index. The e^± injection spectrum takes the form of Q(E) = Q₀E^−γ, where γ ≃ 2.2, as indicated by the HAWC observation (Abeysekara et al. 2017a).

If D is spatially uniform, Equation (1) can be solved analytically using the Green's function method (Ginzburg & Syrovatskii 1964). We verify that our numerical solution for the one-zone case matches well with the analytical solution. For the two-zone diffusion case, there is no analytical way to test the results. To check that there is no "swallowing" or "spitting" of particles in the discontinuity surface, we integrate the number of particles, N(r_i), for the two-zone diffusion cases with different r_⋆. We find that the number of particles is always the same as that of the one-zone cases. The present radial distributions of 1 TeV e^± are shown in Figure 1 for the two-zone diffusion scenarios and also the one-zone cases. Here, we assume that all the spin-down energy of the Geminga pulsar is converted to the energy of injected e^± (conversion efficiency), to determine the normalization.

Zoom In Zoom Out Reset image size

The positron (or electron) spectra at the Earth generated by Geminga are shown in the left panel of Figure 2 for different r_⋆. The conversion efficiency of the Geminga pulsar is also set to be 100%. The spectrum becomes harder for a larger r_⋆, since fewer low-energy e^± can reach the Earth when the slow-diffusion zone is larger. The extreme cases of one-zone diffusion with Y17 or HAWC diffusion coefficients are also shown for comparison. We find that at high energies, the fluxes of the two-zone diffusion models are considerably higher than that of the one-zone diffusion case with Y17 parameters. For a fast-diffusion model like Y17, the propagation scale of 1 TeV e^± of Geminga is ≃1.7 kpc. This scale is about 7 times larger than the distance between Geminga and the Earth, which implies that the particles have diffused to a considerably large region. For the two-zone diffusion scenario, the propagation of particles is hindered within the slow-diffusion zone. After being confined for some time in the inner region, these particles enter the fast propagation region and diffuse to the Earth. Since they spend less time in the fast-diffusion region, they can be effectively accumulated in a relatively smaller region and result in a higher flux at the Earth's location.

Zoom In Zoom Out Reset image size

We find that the two-zone model with r_⋆ = 50 pc can reproduce the AMS-02 data of the positron flux. The results for r_⋆ > 50 pc are too hard to explain the data, while the spectrum is too soft for the case with a smaller r_⋆ or the one-zone, fast-diffusion model. Note that the final spectra depend on the assumed injection spectrum of e^±, which is γ ≃ 2.2 in this work. The comparison between the model and the data is shown in the right panel of Figure 2. Here, the conversion efficiency to e^± is assumed to be 75% in order to match the data. For the secondary contribution from pp collisions, one can refer to Delahaye et al. (2009) and Fang et al. (2018).

The e^± spectrum of Monogem under the two-zone diffusion can be obtained in the same way. The distance, age, and spin-down energy of Monogem are 290 pc, 1.1 × 10⁵ yr, and 1.58 × 10⁴⁸ erg, respectively (Manchester et al. 2005). The observation of HAWC indicates a spectral index of ∼2.0 for the injection spectrum of Monogem (Abeysekara et al. 2017a). As can be seen from Figure 3, the positron spectra of Monogem is too hard to fit the AMS-02 data in the two-zone diffusion cases, due to its farther distance, younger age, and harder injection spectrum compared with Geminga. Thus Monogem cannot make a major contribution to the local positron fluxes even in the two-zone diffusion model.

Zoom In Zoom Out Reset image size

Finally, we would like to give a simple self-consistency check if the two-zone region can reproduce the γ-ray profile observed by HAWC. The HAWC γ-ray data around 20 TeV are produced by ∼100 TeV electrons, whose cooling time is about 7000 yr. This indicates that the parent e^± of the γ-ray observed by HAWC should be younger than 7000 yr. We show in Figure 4 the density profiles of 100 TeV e^± from Geminga with a diffusion of 7000 yr. Within 30 pc from Geminga, we find very good consistency between the one-zone slow-diffusion scenario of Abeysekara et al. (2017a) and the two-zone diffusion scenario with r_⋆ = 50 pc, which means the latter, proposed in this work, can also accommodate the HAWC γ-ray data. The reason for this consistency is that 7000 yr is too short for e^± to escape from the slow-diffusion region of the two-zone model.

Zoom In Zoom Out Reset image size

In this work we propose that a two-zone propagation scenario of cosmic-ray e^± to account for the HAWC observations of the extended halo around the Geminga pulsar and the locally observed positron flux by AMS-02. The diffusion of e^± is assumed to be significantly slow within a distance of r_⋆ from the pulsar, as inferred by the spatial brightness profile of γ-ray emission observed by HAWC (Abeysekara et al. 2017a). The particles diffuse beyond r_⋆ with a diffusion coefficient inferred from the B/C data. A numerical method is adopted to solve the propagation equation, and the differencing scheme of the diffusion operator is derived by the finite volume method, which is different from those used in GALPROP or DRAGON.

The exact value of r_⋆ to fit the AMS-02 positron flux degenerates with the injection spectrum of e^±, and hence can neither be too big nor too small. For larger (smaller) values of r_⋆, low-energy particles would be confined in the slow-diffusion region for longer (shorter), and the final positron spectrum is harder (softer). For an injection spectral index of ∼2.2 and r_⋆ ∼ 50 pc, we find that the e^± from Geminga can reasonably account for the positron excess, with about 75% of the spin-down energy converted into e^±. Actually, such a result is even more natural than that in the one-zone diffusion scenario, which requires an efficiency even slightly larger than 100% (Yin et al. 2013).

The origin of the slow-diffusion zone around Geminga is still unclear. It is possible that the zone pre-existed. For example, the shock of the parent supernova remnant (SNR) of Geminga, which is non-observable today, may have swept the ISM and made it more turbulent. Considering an age of 3 × 10⁵ yr, the shocked region can reach a size of ∼100 pc (Yamazaki et al. 2006). We note that Geminga has a transverse velocity of ${205}_{-47}^{+90}$ km s⁻¹ (Faherty et al. 2007), which means a 70 pc offset from its birthplace. So Geminga may still be within the shocked region of its SNR, where the diffusion coefficient is smaller. Meanwhile, if most of e^± are injected in the early age of Geminga, when the offset of Geminga is small and the scale of its SNR is not large as it is today, our assumption for the symmetrical slow-diffusion zone may not be impacted.

Alternatively, the slow-diffusion zone can be generated intrinsically by the e^± injected from Geminga. Near the cosmic-ray sources, the spatial gradient of particle density is significantly larger than that of the average ISM, which can lead to a growth of the streaming instability (Ptuskin et al. 2008; Malkov et al. 2013; D'Angelo et al. 2016). In this case, the magnetohydrodynamic turbulence may be considerably stronger than that in the ISM. Therefore, it is plausible that the e^± leaving Geminga are confined in the nearby zone for a longer time by the waves induced by themselves.

The slow-diffusion regions around pulsars may be common, which may result in extended γ-ray halos of pulsars and can probably explain the diffuse γ-ray excess observed in the Galactic plane (Abdo et al. 2008; Guo et al. 2016; Linden & Buckman 2018).

Finally, a similar study (Profumo et al. 2018) appeared on arXiv after this paper, which develops a Monte Carlo method to solve the propagation equation, and discusses several kinds of spatial distributions of the propagation coefficient. Their conclusion is consistent with the present work.

This work is supported by the National Key Program for Research and Development (No. 2016YFA0400200) and by the National Natural Science Foundation of China under grants No. U1738209, 11475189, 11475191, 11722328. Q. Yuan acknowledges the 100 Talents Program of the Chinese Academy of Science.

We adopt the operator splitting method to divide the propagation equation into two one-dimensional sub-problems, and choose the Strang splitting which is second-order accurate (Strang 1968). We depict the flow chart as

$$\begin{eqnarray*}&&{N}^{n}\underset{({t}_{n},\,{t}_{n+1/2})}{\overset{{{ \mathcal L }}_{E}}{\longrightarrow }}{\tilde{N}}^{n+1/2}\underset{[({t}_{n},{t}_{n+1})]}{\overset{{{ \mathcal L }}_{r}}{\longrightarrow }}{N}^{n+1/2}\underset{[({t}_{n+1/2},{t}_{n+1})]}{\overset{{{ \mathcal L }}_{E}}{\longrightarrow }}{N}^{n+1},\end{eqnarray*}$$

where n is the time-step index. Kissmann (2014) used trapezoidal integration for the discretization of the energy-loss scheme, which leads to a second-order accuracy. The algebraic equation is written as

$$\begin{eqnarray}{N}_{l,i}^{n+1} & = & \left(\displaystyle \frac{{b}_{l+1}{\rm{\Delta }}t-{\rm{\Delta }}E}{{b}_{l}{\rm{\Delta }}t+{\rm{\Delta }}E}\right){N}_{l+1,i}^{n+1}-\left(\displaystyle \frac{{b}_{l}{\rm{\Delta }}t-{\rm{\Delta }}E}{{b}_{l}{\rm{\Delta }}t+{\rm{\Delta }}E}\right){N}_{l,i}^{n}\\ & & +\,\left(\displaystyle \frac{{b}_{l+1}{\rm{\Delta }}t+{\rm{\Delta }}E}{{b}_{l}{\rm{\Delta }}t+{\rm{\Delta }}E}\right){N}_{l+1,i}^{n},\end{eqnarray} \tag{ 5 }$$

where l is the energy step index, i is the spatial step index, Δt denotes the time-step, and ΔE = E_l+1 − E_l.

For the diffusion operator, we use two different discretizations, depending on if r_i lies on the discontinuity surface. If not, we apply the Crank–Nicolson scheme for constant D, which is second-order accurate and unconditionally stable. The corresponding algebraic equation is

$$\begin{eqnarray} & & \,-\displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i-1}^{n+1}+\left(1+\displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{{\rm{\Delta }}{r}^{2}}\right){N}_{l,i}^{n+1}\\ & & \,-\displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i+1}^{n+1}\\ & = & \displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i-1}^{n}+\left(1-\displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{{\rm{\Delta }}{r}^{2}}\right){N}_{l,i}^{n}\\ & & +\,\displaystyle \frac{{D}_{l,i}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i+1}^{n},\end{eqnarray} \tag{ 6 }$$

where Δr is the radial step size, and r_i = iΔr. Equation (6) is a tridiagonal system, which can be solved with LU decomposition. For the case of r_i = r_⋆, the difference equation should be rewritten to ensure the conservation of the flux at both sides of the discontinuity surface. We adopt the finite volume method to derive the differencing scheme for this case, which reads

$$\begin{eqnarray} & & -\,\displaystyle \frac{{D}_{l,i-1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i-1}^{n+1}\\ & & +\,\left[1+\displaystyle \frac{{D}_{l,i-1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right)+\displaystyle \frac{{D}_{l,i+1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right)\right]{N}_{l,i}^{n+1}\\ & & -\,\displaystyle \frac{{D}_{l,i+1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i+1}^{n+1}\\ & = & \displaystyle \frac{{D}_{l,i-1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i-1}^{n}\\ & & \,+\left[1-\displaystyle \frac{{D}_{l,i-1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1-\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right)-\displaystyle \frac{{D}_{l,i+1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right)\right]{N}_{l,i}^{n}\\ & & \,+\displaystyle \frac{{D}_{l,i+1}{\rm{\Delta }}t}{2{\rm{\Delta }}{r}^{2}}\left(1+\displaystyle \frac{{\rm{\Delta }}r}{{r}_{i}}\right){N}_{l,i+1}^{n}.\end{eqnarray} \tag{ 7 }$$

The derivation of this scheme is presented as follows.

Consider a diffusion equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial N}{\partial t}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}D(r)\displaystyle \frac{\partial N}{\partial r}\right]=0,\end{eqnarray} \tag{ 8 }$$

where D(r) is defined in Equation (2). We set F = r²D(r)∂N/∂r, which is the particle flux in r. Then Equation (8) can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{Nr}}^{2}}{\partial t}=\displaystyle \frac{\partial F}{\partial r}.\end{eqnarray} \tag{ 9 }$$

Assuming r_i = r_⋆, the continuity of particle density and flux should be satisfied:

$$\begin{eqnarray}&&\left\{\begin{array}{l}N(t,{r}_{i}-0)=N(t,{r}_{i}+0),\\ F(t,{r}_{i}-0)=F(t,{r}_{i}+0).\end{array}\right.\end{eqnarray} \tag{ 10 }$$

We integrate Equation (9) in the time range of [t_n, t_n+1] and the radial range of [r_{i − 1/2}, r_i+1/2], where r_i±1/2 = r_i ± Δr/2, and we obtain

$$\begin{eqnarray} & & \,{\int }_{{r}_{i-1/2}}^{{r}_{i+1/2}}[N({t}_{n+1},r)-N({t}_{n},r)]{r}^{2}{dr}\\ & = & {\int }_{{t}_{n}}^{{t}_{n+1}}\left[{\int }_{{r}_{i-1/2}}^{{r}_{i+1/2}}\displaystyle \frac{\partial }{\partial r}F(t,r){dr}\right]{dt}\\ & = & {\int }_{{t}_{n}}^{{t}_{n+1}}\left[{\int }_{{r}_{i-1/2}}^{{r}_{i}}\displaystyle \frac{\partial }{\partial r}F(t,r){dr}+{\int }_{{r}_{i}}^{{r}_{i+1/2}}\displaystyle \frac{\partial }{\partial r}F(t,r){dr}\right]{dt}\\ & = & {\int }_{{t}_{n}}^{{t}_{n+1}}\{[F(t,{r}_{i}-0)-F(t,{r}_{i-1/2})]\\ & & +\,[F(t,{r}_{i+1/2})-F(t,{r}_{i}+0)]\}{dt}\\ & = & {\int }_{{t}_{n}}^{{t}_{n+1}}[F(t,{r}_{i+1/2})-F(t,{r}_{i-1/2})]{dt}.\end{eqnarray} \tag{ 11 }$$

On the other hand, we have ∂N/∂r = F/[r²D(r)] as defined above. Integrating this relation in [r_i, r_i+1], we get

$$\begin{eqnarray}N(t,{r}_{i+1})-N(t,{r}_{i}) & = & {\int }_{{r}_{i}}^{{r}_{i+1}}\displaystyle \frac{F(t,r)}{{r}^{2}D(r)}{dr}\\ & \approx & F(t,{r}_{i+1/2}){\int }_{{r}_{i}}^{{r}_{i+1}}\displaystyle \frac{1}{{r}^{2}D(r)}{dr},\end{eqnarray} \tag{ 12 }$$

which means

$$\begin{eqnarray}&&F(t,{r}_{i+1/2})\approx \displaystyle \frac{N(t,{r}_{i+1})-N(t,{r}_{i})}{{\int }_{{r}_{i}}^{{r}_{i+1}}\tfrac{1}{{r}^{2}D(r)}{dr}}.\end{eqnarray} \tag{ 13 }$$

The expression of F(t,r_{i − 1/2}) can be found similarly. The integral terms in Equation (11) can be approximated by

$$\begin{eqnarray}&&{\int }_{{r}_{i-1/2}}^{{r}_{i+1/2}}N({t}_{n},r){r}^{2}{dr}\approx N({t}_{n},{r}_{i}){r}_{i}^{2}{\rm{\Delta }}r\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{\int }_{{t}_{n}}^{{t}_{n+1}}F(t,{r}_{i+1/2}){dt}\approx [F({t}_{n+1},{r}_{i+1/2})\\ &&\quad +F({t}_{n},{r}_{i+1/2})]{\rm{\Delta }}t/2,\end{eqnarray} \tag{ 15 }$$

where the latter corresponds to the Crank–Nicolson scheme. Combining Equations (11) and (13)–(15), we finally obtain

$$\begin{eqnarray} & & \displaystyle \frac{N({t}_{n+1},{r}_{i})-N({t}_{n},{r}_{i})}{{\rm{\Delta }}t}\\ & = & \displaystyle \frac{{D}_{i+1}{r}_{i+1}}{2{\rm{\Delta }}{r}^{2}{r}_{i}}[N({t}_{n+1},{r}_{i+1})-N({t}_{n+1},{r}_{i})\\ & & \,+N({t}_{n},{r}_{i+1})-N({t}_{n},{r}_{i})]\\ & & -\,\displaystyle \frac{{D}_{i-1}{r}_{i-1}}{2{\rm{\Delta }}{r}^{2}{r}_{i}}[N({t}_{n+1},{r}_{i})-N({t}_{n+1},{r}_{i-1})\\ & & +\,N({t}_{n},{r}_{i})-N({t}_{n},{r}_{i-1})],\end{eqnarray} \tag{ 16 }$$

which can be rewritten as Equation (7).

Finally, the initial condition is approximated by a step function

$$\begin{eqnarray}N(0,E,{r}_{i})=\left\{\begin{array}{ll}Q^{\prime} (E), & i\leqslant 1\\ 0, & i\gt 1\end{array}\right.,\end{eqnarray} \tag{ 17 }$$

where $Q^{\prime} (E)=Q(E)/(4\pi {r}_{1}^{3})$. For the inner spatial boundary condition, we handle it with N(r₁) = N(r₋₁). Plugging this relation into Equation (6), we can then obtain N(r₀).