Inefficient Cosmic-Ray Diffusion around Vela X: Constraints from H.E.S.S. Observations of Very High-energy Electrons

Recent measurements of an increasing cosmic-ray (CR) positron (e⁺) fraction above 10 GeV by the Payload for Antimatter Exploration and Light-nuclei Astrophysics (PAMELA) and Alpha Magnetic Spectrometer (AMS-02) identify an excess of high-energy positrons relative to the standard predictions for secondary electron–positron production in the interstellar medium (ISM; Adriani et al. 2010; Aguilar et al. 2013), suggesting that these positrons are produced as primary particles. The CR e⁺+e⁻ spectrum has been measured up to a few TeV by Fermi-Large Area Telescope (-LAT), the Very Large Energetic Radiation Imaging Array System (VERTIAS), the Dark Matter Particle Explorer (DAMPE) (Abdollahi et al. 2017; Staszak & VERITAS Collaboration 2015; DAMPE Collaboration et al. 2017), and to ∼20 TeV by HESS recently (H.E.S.S. Collaboration et al. 2017). As high-energy electrons (hereafter, we do not distinguish positrons from electrons) efficiently cool in the ISM through synchrotron and inverse-Compton radiation, they cannot travel beyond a kiloparsec distance before depleting their energies. Therefore, these high-energy electrons must be produced by nearby sources, such as pulsars (or PWNe), annihilating dark matter particles, and supernova remnants (SNR). While the real sources are still under debate, nearby pulsars such as Geminga, PSR B0656+14, and Vela, which are at distances of only ∼200 pc, are the most attractive candidates (Shen 1970; Atoyan et al. 1995; Hooper et al. 2017).

Vela X, powered by a ∼10⁴ year old pulsar (i.e., the Vela pulsar or PSR B0833-45), is one of the most well-studied PWN. It consists of an extended radio halo of the size of 2° × 3° and a small collimated structure, e.g., the "cocoon." High-energy gamma-ray emission has been detected from both the halo and the cocoon by Fermi-LAT and HESS. de Jager et al. (2008) propose that there are two distinct populations of electrons: one responsible for the radio and GeV gamma-ray emission and the other for the X-ray and TeV emission. To explain the steep GeV spectrum measured by Fermi-LAT and the dimness of the TeV nebula relative to the spin-down power of the Vela pulsar, Hinton et al. (2011) argue that a significant diffusive escape of electrons from the halo must have occurred. The same conclusion is reached more recently by Tibaldo et al. (2018) with an analysis of ∼9.5 years of data from Fermi-LAT observations. While the confinement of particles in PWNe is thought to be effective in the early stage, the interaction with the SNR reverse shock, which seems to have appeared in Vela X several thousand years ago (Blondin et al. 2001; Gelfand et al. 2009), may have brought an end to the confinement. The asymmetric structure of the PWN with respect to the pulsar supports such an interpretation. The interaction is expected to disrupt the PWN sufficiently enough that diffusion of particles out of the PWN becomes possible.

The escaped electrons will contribute to the CRe spectrum measured at the Earth. Using a diffusion coefficient of 1.07 × 10²⁷ (E/1 GeV)^0.6 cm² s⁻¹ and a total electron energy of 6.8 × 10⁴⁸ erg, Hinton et al. (2011) predict a distinct bump in the CR electron spectrum at several TeV. This diffusion coefficient is actually smaller than the standard one for the ISM, which is ${D}_{\mathrm{ISM}}\simeq 3.86\times {10}^{28}{({E}_{e}/\mathrm{GeV})}^{0.33}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, as inferred from the boron-to-carbon ratio and other CR secondary-to-primary ratios for electrons with E ≲ 10 TeV (see http://galprop.stanford.edu/). If one uses the standard diffusion coefficient, the bump would be even higher. However, the latest results from HESS do not show such a bump (H.E.S.S. Collaboration et al. 2017), which invalidates the simple diffusion model for the escape of particles from Vela X.

Interestingly, recent TeV observations of Geminga PWN with the High-Altitude Water Cheronkov Observatory (HAWC) telescope have been interpreted as evidence that diffusion of high-energy electrons within PWN is highly inefficient compared to the standard value for the ISM (Abeysekara et al. 2017; Hooper et al. 2017). Motivated by this discovery, we here study whether a spatially dependent diffusion around Vela X can lead to a CR electron flux consistent with that measured by HESS. We make a two-zone approximation for the diffusion coefficient surrounding Vela X, i.e., an inner inefficient diffusion zone and an outer standard diffusion zone. We present an analytic approach to calculate the CR electron flux at the Earth. Although this analytic approach is a crude approximation, it provides a convenient way to estimate the CR flux for the two-zone model. In Section 2, we first present the difficulties for the simple one-zone model with a standard diffusion coefficient. Then, we study the two-zone model and obtain the constraint on the diffusion coefficient of the inner zone around Vela X in Section 3. Finally, we give the discussions and conclusions in Section 4.

We first study whether a simple one-zone diffusion model, which is usually used for CR propagation in the ISM, can produce a flux consistent with the HESS measurement. This simple model assumes that the diffusion is homogeneous along the path from the PWN to the Earth. The diffusion coefficient is given by the standard one, ${D}_{\mathrm{ISM}}\simeq 3.86\times {10}^{28}{({E}_{e}/\mathrm{GeV})}^{0.33}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$.

The transport of CR electrons can be described by the equation:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}{n}_{e}({E}_{e},\overrightarrow{x},t) & = & \overrightarrow{{\rm{\nabla }}}\cdot [D({E}_{e},\overrightarrow{x})\overrightarrow{{\rm{\nabla }}}{n}_{e}({E}_{e},\overrightarrow{x},t)]\\ & & +\ \displaystyle \frac{\partial }{\partial {E}_{e}}\left[\displaystyle \frac{{{dE}}_{e}}{{dt}}{n}_{e}({E}_{e},\overrightarrow{x},t)\right]\\ & & +\ Q({E}_{e},\overrightarrow{x},t),\end{array}\end{eqnarray} \tag{ 1 }$$

where n_e(E) is the differential number density of electrons, D(E_e) is the diffusion coefficient, and Q(E_e, $\overrightarrow{x}$, t) is the source term. Energy losses caused by inverse-Compton (IC) and synchrotron processes are described as

$$\begin{eqnarray}\begin{array}{rcl}-\displaystyle \frac{{{dE}}_{e}}{{dt}}(r) & = & \displaystyle \sum \displaystyle \frac{4}{3}{\sigma }_{{\rm{T}}}{\rho }_{i}(r){S}_{i}({E}_{e}){\left(\displaystyle \frac{{E}_{e}}{{m}_{e}}\right)}^{2}\\ & & +\ \displaystyle \frac{4}{3}{\sigma }_{{\rm{T}}}{\rho }_{\mathrm{mag}}(r){\left(\displaystyle \frac{{E}_{e}}{{m}_{e}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

where σ_T is the Thomson cross section, ρ_i denotes the radiation energy density of background photons, and ρ_mag denotes the energy density of the magnetic field. Various components of the radiation backgrounds are taken into consideration, including the cosmic microwave background (CMB), starlight (star), ultraviolet emission (UV), and infrared emission (IR). Parameters are adopted as follows for the area surrounding Vela: ${\rho }_{\mathrm{CMB}}=0.260\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$, ${\rho }_{\mathrm{star}}=0.44\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$, ρ_UV = 0.10 eV cm⁻³, ${\rho }_{\mathrm{IR}}=0.44\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$, ${\rho }_{\mathrm{mag}}=0.622\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$ (corresponding to B = 5 μG) T_CMB = 2.7 K, T_star = 7500 K, T_UV = 20000 K, and T_IR = 25 K (de Jager et al. 2008; Grondin et al. 2013). We also check the results adopting other values of these parameters (Hooper et al. 2017; Fang et al. 2018b), but negligible differences are found. When ${E}_{e}\gtrsim {m}_{e}^{2}/2T$, the suppression of the inverse-Compton scattering by the Klein–Nishina effect cannot be ignored, which can be parameterized by (Longair 2011)

$$\begin{eqnarray}&&{S}_{i}({E}_{e})\approx \displaystyle \frac{45{m}_{e}^{2}/64{\pi }^{2}{T}_{i}^{2}}{(45{m}_{e}^{2}/64{\pi }^{2}{T}_{i}^{2})+({E}_{e}^{2}/{m}_{e}^{2})}.\end{eqnarray} \tag{ 3 }$$

For a burst-like injection of $Q({E}_{e},t)=\delta (t){Q}_{0}{E}_{e}^{-\alpha }\exp (-{E}_{e}/{E}_{c})$, the solution to Equation (1) is given by

$$\begin{eqnarray}&&{n}_{e}({E}_{e},r,t)=\displaystyle \frac{{Q}_{0}{E}_{0}^{2-\alpha }{e}^{-{E}_{0}/{E}_{c}}}{8{\pi }^{3/2}{E}_{e}^{2}{L}_{\mathrm{dif}}^{3}({E}_{e},t)}\exp \left[-\displaystyle \frac{{r}^{2}}{4{L}_{\mathrm{dif}}^{2}({E}_{e},t)}\right],\end{eqnarray} \tag{ 4 }$$

where E₀ is the initial energy of the electron of energy E_e at t and L_dif is the diffusion length scale given by

$$\begin{eqnarray}&&{L}_{\mathrm{dif}}\equiv {\left[{\int }_{{E}_{0}}^{{E}_{e}}\displaystyle \frac{D({E}^{{\prime} })}{-{{dE}}_{e}/{dt}({E}^{{\prime} })}{{dE}}^{{\prime} }\right]}^{1/2}.\end{eqnarray} \tag{ 5 }$$

We note that Equation (4) permits a fraction of particles to propagate faster than the speed of light, which is the so-called "superluminal diffusion problem" (Dunkel et al. 2007; Aloisio et al. 2009; Liu et al. 2016). From the mathematical point of view, the superluminal propagation always exists in the solutions of the nonrelativistic diffusion equations, but very frequently the contribution of unphysical regions to the solution is negligibly small. In these cases, one can regard that the diffusion equation gives a correct description of the considered physical phenomenon. However, when t ∼ r/c, the problem of superluminal diffusion in Equation (4) becomes severe. Aloisio et al. (2009) find a solution to this problem: if the cooling of particles is unimportant (which is applicable to our case), the probability to find one electron at a distance, r, from the source at a time, t, after its injection can be described by

$$\begin{eqnarray}&&P({E}_{e},t,r)=\displaystyle \frac{\theta ({ct}-r)}{{({ct})}^{3}Z(\tfrac{{c}^{2}t}{2D}){\left[1-{\left(\tfrac{r}{{ct}}\right)}^{2}\right]}^{2}}\exp \left[-\displaystyle \frac{\tfrac{{c}^{2}t}{2D}}{\sqrt{1-{\left(\tfrac{r}{{ct}}\right)}^{2}}}\right],\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&Z(y)=4\pi {K}_{1}(y)/y,\end{eqnarray} \tag{ 7 }$$

with θ(r) as the Heaviside function and K₁(y) as the modified Bessel function. Note that the above formula only works when t is much smaller than the cooling timescale of electrons, which is true for our following calculation. The electron density at time t after the injection then can be obtained by

$$\begin{eqnarray}&&{n}_{e}({E}_{e},r,t)={\int }_{-\infty }^{t}P({E}_{e},t-t^{\prime} ,r)\delta (t^{\prime} ){Q}_{0}{E}_{e}^{-\alpha }{e}^{-{E}_{e}/{E}_{c}}{dt}^{\prime} .\end{eqnarray} \tag{ 8 }$$

Following Hinton et al. (2011), we take α = 1.8 and E_c = 6 TeV for the electron spectrum. The total energy injected into the initially confined PWN depends on the birth-period, P₀, of the pulsar. The energy in relativistic electrons is about $\sim 2\times {10}^{49}\epsilon {({P}_{0}/30\mathrm{ms})}^{-2}$, where is the fraction of spin-down power converted into relativistic electrons. All previous estimates of the total energy give a value of several 10⁴⁸ erg (de Jager et al. 2008; Abdo et al. 2010; Hinton et al. 2011; Grondin et al. 2013). This is consistent with the estimate of the birth-period of P₀ = 40 ms that is invoked to account for the ratio between the PWN radius and the SNR radius (van der Swaluw & Wu 2001). In the following calculation, we use a total energy of 6.8 × 10⁴⁸ erg, as obtained from the spectral energy distribution (SED) fit of Vela X by Hinton et al. (2011). We use r = 287 pc as the distance between the Earth and Vela. We assume that the electrons are released instantaneously at disruption time t of the initial PWN. The exact disruption time of Vela X is unknown, but the reasonable values should be at least several kyr (e.g., Blondin et al. 2001; Gelfand et al. 2009). Note that the time span between the disruption time of Vela X and the current time is exactly the propagation time of injected electrons t.

The results of the CR electron spectrum for various t are shown in Figure 1, where Equation (8) is used to avoid the possible superluminal diffusion problem.⁶ Obviously, the CR flux produced by Vela X exceed the measured flux by several order of magnitudes for any reasonable value of t. Only for a very small t (i.e., t < 1100 years), the CR flux can be consistent with the HESS measurement (see the blue line in Figure 1). Such a small t seems unreasonable since the timescale when the PWN collides with the SNR reverse shock explosion is usually several thousand years (e.g., Blondin et al. 2001; Gelfand et al. 2009). This suggests that the simple one-zone diffusion model does not work for Vela X.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recent HAWC observations of the PWN regions around Geminga and Monogem provide detailed information about the spatially extended emission of TeV gamma-rays. The spectrum and morphology of the TeV emission can be used to infer the features of underlying high-energy electrons responsible for the IC photons. The angular profiles of the TeV emission observed from Geminga and Monogem indicate that the diffusion is highly inefficient in the regions surrounding these sources (Abeysekara et al. 2017). This is also the first empirical determination of a diffusion coefficient in the region of tens of pc around pulsars in the local Galaxy. The inferred diffusion coefficient is more than two orders of magnitude smaller than the standard diffusion coefficient in the ISM. Furthermore, assuming a spatially dependent diffusion, Fang et al. (2018a) and Profumo et al. (2018) show that nearby pulsars, such as Geminga, could contribute significantly to the CR electron spectrum in 0.1–1 TeV. Motivated by these results, we study whether a spatially dependent diffusion model with inefficient diffusion in the inner region surrounding Vela X could resolve the inconsistency between the predicted CR electron flux and the observed one by HESS. We approximate the diffusion coefficient of the entire space to be a step function of the distance to Vela X, where the diffusion is suppressed with a coefficient of D₁ within a few tens of parsecs from Vela X and a standard diffusion coefficient D₂ for the outside region, i.e.,

$$\begin{eqnarray}D(r)=\left\{\begin{array}{ll}{D}_{1}, & r\leqslant {r}_{0}\\ {D}_{2}\,= & {D}_{\mathrm{ISM}},\,r\gt {r}_{0}\end{array}\right..\end{eqnarray} \tag{ 9 }$$

We develop an analytic approach to solve Equation (1) in the above two-zone diffusion model, as shown in Figure 2. We first use the solution of Equation (4) with a diffusion coefficient, D₁, to estimate the number density of electrons at the interface of the two zones (i.e., at the spherical surface with a radius, r₀). We note that Equation (4) is strictly correct only when D₁ = D₂. So the above obtained number density is a crude estimate when ${D}_{1}\ne {D}_{2}$. The flux density at the sphere is ${F}_{D}=-D\tfrac{\partial {n}_{e}}{\partial r}$, and the number of particles in unit time passing an surface element on the sphere outwardly is ${\rm{\Delta }}\dot{N}={F}_{D}{\rm{\Delta }}S$, with ΔS as the area of the element surface. We then regard the surface element as a point source with an injection rate, ${\rm{\Delta }}\dot{N}$. Since the flux at the sphere is a function of time, ${\rm{\Delta }}\dot{N}$ is also a function of time. We then regard the sphere consist of many surface elements, each of which is treated as a point source. A convenient choice is to divide the sphere into annular rings perpendicular to the line connecting the Earth and Vela X, so that the distance between the Earth and each surface element on the ring is the same. The radius of the ring is r₀ · sinθ, and the distance to the Earth is $d={({r}^{2}-2r\cdot {r}_{0}\cos \theta +{r}_{0}^{2})}^{1/2}$. The area of this ring is ${\rm{\Delta }}{S}_{\mathrm{ring}}=2\pi {r}_{0}^{2}\sin \theta d\theta$. The integral can be done over the annular rings on the spherical surface. We further regard the continuous injection from the sphere as the sum of a series of discrete injection so that the electron number density at the radius of the Earth can be calculated with Equation (4). The precision of this analytic method is tested in the Appendix. We find that the difference of the CR electron flux at the Earth between our analytic results and the numerical approach used by Fang et al. (2018b) is, at most, ∼30%. Since the uncertainty of the CR flux measured by HESS at ∼10 TeV is about a factor of two, we consider that the precision of our analytic approach is sufficient for this study.

We note that the burst-like injection is an assumption based on the reverse shock interaction scenario. To be more complete, we also consider the case of continuous injection. Since the injection rate is most likely to be monotonically decreasing, a flat injection profile (i.e., a constant injection rate of electrons) would provide the most conservative estimate of the electron flux arriving at the Earth, given that other parameters are the same. That is, the burst-like profile and the flat-time profile can be regarded as the two extremes for any injection patterns. To deal with the flat injection, we decompose the injection function into many small time bins, treat each time bin as a burst-like injection, and last sum up the contribution from each of them. We compare the results of the burst-like injection and the flat injection in Figure 3 for a different time, t, at which electron injection started. We find that, in the case of the flat-profile injection, a larger diffusion coefficient is obtained. This is due to that a significant part of electrons are injected at a later time in this case, and these electrons have not arrived at the Earth. The peak of the spectrum moves toward lower energy for earlier injection (i.e., larger t). This is because electrons with lower energy can reach the Earth while electrons with higher energy will be cooled down. We find that, for all the t considered here, the difference in the obtained D₁ for the two different injection profiles is within a factor of 2. We thus conclude that our result is not sensitive to the injection profile, and in what follows, we only consider the burst-like injection profile for simplicity.

Zoom In Zoom Out Reset image size

The disruption time t of the PWN and the radius of the inner zone are two unknown parameters. We use r₀ = 30 pc and r₀ = 50 pc as two reference values for the radius of the inner zone. The CR fluxes for the two-zone model with r₀ = 30 pc and r₀ = 50 pc are, respectively, shown in the left and right panels of Figure 4. For a reasonable range of t (i.e., t > 2.5 kyr), we find that only when ${D}_{1}(10\,\mathrm{TeV})\lesssim {10}^{28}\,\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ can the CR electron flux be consistent with the HESS measurement. This value is about two orders of magnitude smaller than the standard diffusion coefficient at ∼10 TeV, which is ${D}_{\mathrm{ISM}}(10\,\mathrm{TeV})\simeq 0.8\,\times {10}^{30}\,\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$. The constraint on the diffusion coefficient becomes more stringent for a larger t. This is because a longer diffusion time requires a smaller diffusion coefficient in the inner zone.

Zoom In Zoom Out Reset image size

In the above calculation, we have assumed D(E) ∝ E^0.33. As the physics of the diffusion is not well-known, we also consider another case of the energy dependence, i.e., D(E) ∝ E^0.548 for both the ISM and the region surrounding the PWN (Fang et al. 2017). The results are shown in the left and right panels of Figure 5 for r₀ = 30 pc and r₀ = 50 pc, respectively. We find that, for this kind of diffusion, we also have ${D}_{1}(10\,\mathrm{TeV})\lesssim {10}^{28}\,\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, which is almost the same as the D(E) ∝ E^0.33 case.

Zoom In Zoom Out Reset image size

In the above, we have taken a total injection energy of 6 × 10⁴⁸, as suggested by Hinton et al. (2011). To be more conservative and motivated by HAWC observations of other TeV halos, we also take a total injection energy corresponding to ∼10% of the total energy budget, i.e., 1.9 × 10⁴⁸ erg. The results are shown in Figure 6. As shown, the required D₁ does not change significantly, as the electron/positron flux is more sensitive to D₁ than the normalization of the injection spectrum.

Zoom In Zoom Out Reset image size

The cutoff energy of the electron spectrum has been taken to be E_c = 6 TeV following Hinton et al. (2011). We note that the spectra observed from Geminga and B0656+14 by HAWC favor much larger values of E_c, i.e., several tens of TeV. Results with different E_c are shown in Figure 7. For a larger E_c, more electrons with higher energy are injected. As these energetic electrons can more easily escape from the slow-diffusion region and reach the Earth, the received flux of the higher-energy electrons are significantly larger. Therefore, the constraint on the diffusion coefficient becomes even more stringent, and a smaller D₁ is required to reconcile with the HESS data.

Zoom In Zoom Out Reset image size

In all of these cases, ${D}_{1}(10\,\mathrm{TeV})\lesssim {10}^{28}\,\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ can be obtained. Thus, the limit on the diffusion coefficient is robust and independent of uncertainties in the input parameters. We conclude that the diffusion in the immediate vicinity of Vela X must be highly inefficient.

There have been suggestions that electrons must have undergone diffusive escape from the the Vela X PWN, indicated by the evidence for a roll-over of the electron spectrum at energies of a few tens of GeV. These escaped electrons may contribute to the CR electron flux at the Earth. In this paper, we have shown that recent HESS data of the CR electron flux at ∼10 TeV place interesting constraints on the diffusion coefficient around Vela X. We find that a highly inefficient diffusion region in the immediate vicinity of Vela X must be present, with D(10 TeV) ≲ 10²⁸ cm² s⁻¹. The result is consistent with the recent finding that there are inefficient diffusion regions surrounding the Geminga and PSR B0656+14 PWNe, suggesting that such inefficient diffusion regions may be common around PWNe of various ages.

Previous theoretical studies have suggested that CRs can be scattered by the self-generated Alfvén waves induced by streaming instability (Yan & Lazarian 2004; Yan et al. 2012; Malkov et al. 2013; Quenby 2018), whereby the particles are self-confined. If the CR flux is sufficiently high, the growth rate of the streaming instability can dominate the nonlinear damping of background turbulence. We speculate that such a process enhances the CR scattering rate inside r₀ and hence reduces the diffusion coefficient to the required level. Another possible mechanism is the influences of the Vela SNR surrounding Vela X. The turbulence in the areas swept up by the shock of the SNR can be strong, which may induce a smaller diffusion coefficient (Bell 1978; Li & Chen 2010, 2012).

We thank Dan Hooper and Huirong Yan for helpful discussions. X.Y.W. is supported by the National Key R & D program of China under the grant 2018YFA0404200, 973 program under grant 2014CB845800, and the NSFC under grants 11625312 and 11851304.

We first compare the CR flux obtained with our analytic approach for the two-zone model with D₁ = D₂ and the simple one-zone model with the same diffusion coefficient (i.e., D = D₁). The result is shown in Figure 8. One can see that the difference between the two models is negligibly small. This demonstrates that the treatment described by Figure 2 is quite accurate.

Zoom In Zoom Out Reset image size

However, this test does not account for the influence caused by the difference in the diffusion coefficients of the two zones (i.e., ${D}_{1}\ne {D}_{2}$). There are two factors that may cause differences between the results in our analytical model and that of the realistic case. The density gradient of electrons at r₀ calculated by Equation (4) becomes smaller when D₂ > D₁, since a larger diffusion coefficient in the outer zone leads to a faster diffusion outward. Thus, our model underestimates the flux escaping from the spherical surface at r₀ in this regard. On the other hand, part of the electrons located at the spherical surface at r₀ will go back to the inefficient diffusion zone. These electrons will take a longer time to arrive at the Earth. Our analytic estimate does not take into account this effect, so our result overestimate the CR flux at the Earth in this regard.

To test the precision for the case of ${D}_{1}\ne {D}_{2}$, we compare our analytic result with the numerical result obtained by Fang et al. (2018a), which solves Equation (1) for the case of ${D}_{1}\ne {D}_{2}$ with a numerical method. The results are shown in Figure 9. For the cases of t = 2500 years, the difference between our analytic result and that obtained with the numerical method is, at most, 30%. For the cases of t = 5000 years, the difference is even smaller.

Zoom In Zoom Out Reset image size