Understanding the multiwavelength observation of Geminga's TeV halo: the role of anisotropic diffusion of particles

TL;DR

This paper explains Geminga's TeV halo and X-ray observations through anisotropic diffusion of electrons/positrons, emphasizing the role of local magnetic field alignment and turbulence without extreme environmental assumptions.

Contribution

It introduces a model of anisotropic diffusion aligned with the magnetic field to interpret Geminga's multiwavelength observations, avoiding extreme environmental conditions.

Findings

Anisotropic diffusion explains the observed TeV and X-ray emissions.

Sub-Alfvénic turbulence with aligned magnetic fields suffices for the model.

No need for extremely low diffusion coefficients or weak magnetic fields.

Abstract

In this letter we propose that the X-ray and the TeV observations in the vicinity of Geminga can be understood in the framework of anisotropic diffusion of injected electrons/positrons. This interpretation only requires the turbulence in the vicinity of Geminga to be sub-Alfv\'enic with the local mean magnetic field direction approximately aligned with our line of sight towards Geminga, without invoking extreme conditions for the environment, such as an extremely small diffusion coefficient and a weak magnetic field of strength $<1\mu$G as suggested in previous literature.

Figures

Equations

$$\frac{\partial N}{\partial t}=\nabla \cdot (\mathcal{D}\cdot \nabla N)-\frac{\partial }{\partial E_e}\left({\dot{E}}_{e}N\right)+Q$$

$$\begin{array}{rl}\frac{\partial N}{\partial t}=& \frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_{rr}\frac{\partial N}{\partial r}\right)+D_{zz}\frac{{\partial }^{2}N}{\partial z^2}\\ & -\frac{\partial }{\partial E_e}\left({\dot{E}}_{e}N\right)+Q(E_e)S(t)\delta (r)\delta (z).\end{array}$$

$$D_{zz}=D_{\parallel }=D_0(E_e/1\mathrm{GeV}{)}^q$$

$$D_{rr}=D_{\perp }=D_{zz}{M}_A^4$$

$$x_E=(d_{\mathrm{gem}}-l\cos \theta )\sin \varphi -l\sin \theta \cos \zeta \cos \varphi$$

$$y_E=l\sin \theta \sin \zeta$$

$$z_E=(d_{\mathrm{gem}}-l\cos \theta )\cos \varphi +l\sin \theta \cos \zeta \sin \varphi$$

$$dN(E,z,r)=N(E,z,r)(l\sin \theta d\zeta \cdot ld\theta \cdot dl)$$

$$\begin{array}{rl}d{F}_{\mathrm{syn}}(ϵ)=& \int {\mathcal{F}}_{\mathrm{syn}}\left\{dN(E,z,r),B\sin \alpha \right\}\\ & \times f({\cos }^2\delta )d{\cos }^2\delta /4\pi l^2\end{array}$$

$$d{F}_{\mathrm{IC}}(ϵ)={\mathcal{F}}_{\mathrm{IC}}\left\{dN(E,z,r),n_{\mathrm{ph}}\right\}/4\pi l^2$$

$$\begin{array}{rl}{I}_{\mathrm{syn}}(ϵ,\theta ,\zeta )=& \frac{1}{4\pi }{\int }_{{\cos }^2\delta }{\int }_{l_{\min }}^{l_{\max }}{\mathcal{F}}_{\mathrm{syn}}\left\{N(E,z,r),B_0\sin \alpha \right\}\\ & \times f({\cos }^2\delta )d{\cos }^2\delta dl\end{array}$$

$$I_{\mathrm{IC}}(ϵ,\theta ,\zeta )=\frac{1}{4\pi }{\int }_{l_{\min }}^{l_{\max }}{\mathcal{F}}_{\mathrm{IC}}\left\{N(E,z,r),n_{\mathrm{ph}}\right\}dl,$$

$$\frac{\partial N}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_{rr}\frac{\partial N}{\partial r}\right)+D_{zz}\frac{{\partial }^{2}N}{\partial z^2}-\frac{\partial }{\partial E_e}\left({\dot{E}}_{e}N\right)+Q(E_e)S(t)\delta (r)\delta (z),$$

$$N^l\underset{(t_l,t_{l+1/2})}{\overset{{\mathcal{L}}_E}{\to }}{\tilde{N}}^{l+1/2}\underset{(t_l,t_{l+1})}{\overset{{\mathcal{L}}_r}{\to }}{N}^{l+1/2}\underset{(t_{l+1/2},t_{l+1})}{\overset{{\mathcal{L}}_E}{\to }}{N}^{l+1},$$

$$\begin{split}\frac{N_{i,j,k}^{l+1}-N_{i,j,k}^{l}}{\Delta t}&=\frac{1}{2}\Bigg{[}\frac{-b_{i,j,k+2}N_{i,j,k+2}^{l+1}+4b_{i,j,k+1}N_{i,j,k+1}^{l+1}-3b_{i,j,k}N_{i,j,k}^{l+1}}{2\Delta E}\\ &+\frac{-b_{i,j,k+2}N_{i,j,k+2}^{l}+4b_{i,j,k+1}N_{i,j,k+1}^{l}-3b_{i,j,k}N_{i,j,k}^{l}}{2\Delta E}\Bigg{]}.\end{split}$$

$$\begin{split}N_{i,j,k}^{l+1}&=\left[\frac{-b_{i,j,k+2}N_{i,j,k+2}^{l+1}+4b_{i,j,k+1}N_{i,j,k+1}^{l+1}-3b_{i,j,k}N_{i,j,k}^{l+1}-b_{i,j,k+2}N_{i,j,k+2}^{l}+4b_{i,j,k+1}N_{i,j,k+1}^{l}}{4\Delta E}\Delta t+N_{i,j,k}^{l}\right]\\ &\Big{/}\left[1+\frac{3b_{i,j,k}\Delta t}{4\Delta E}\right].\end{split}$$

$$\begin{array}{rl}& \frac{N_{i,j,k}^{l+1}-N_{i,j,k}^l}{\Delta t}=\frac{D_{rr,k}}{2r_i}\left(\frac{N_{i+1,j,k}^l-N_{i-1,j,k}^l}{2\Delta r}+\frac{N_{i+1,j,k}^{l+1}-N_{i-1,j,k}^{l+1}}{2\Delta r}\right)\\ & +\frac{D_{rr,k}}{2}\left(\frac{N_{i+1,j,k}^l-2N_{i,j,k}^l+N_{i-1,j,k}^l}{\Delta r^2}+\frac{N_{i+1,j,k}^{l+1}-2N_{i,j,k}^{l+1}+N_{i-1,j,k}^{l+1}}{\Delta r^2}\right)\\ & +\frac{D_{zz,k}}{2}\left(\frac{N_{i,j+1,k}^l-2N_{i,j,k}^l+N_{i,j-1,k}^l}{\Delta z^2}+\frac{N_{i,j+1,k}^{l+1}-2N_{i,j,k}^{l+1}+N_{i,j-1,k}^{l+1}}{\Delta z^2}\right)\end{array}$$

$$\begin{array}{rl}\frac{N_{i,j,k}^{l+1/2}-N_{i,j,k}^l}{\Delta t}& =\frac{D_{rr,k}}{2r_i}\left(\frac{N_{i+1,j,k}^{l+1/2}-N_{i-1,j,k}^{l+1/2}}{2\Delta r}\right)+\frac{D_{rr,k}}{2}\left(\frac{N_{i+1,j,k}^{l+1/2}-2N_{i,j,k}^{l+1/2}+N_{i-1,j,k}^{l+1/2}}{\Delta r^2}\right)\\ & +\frac{D_{zz,k}}{2}\left(\frac{N_{i,j+1,k}^l-2N_{i,j,k}^l+N_{i,j-1,k}^l}{\Delta z^2}\right)\end{array}$$

$$\begin{array}{rl}\frac{N_{i,j,k}^{l+1}-N_{i,j,k}^{l+1/2}}{\Delta t}& =\frac{D_{rr,k}}{2r_i}\left(\frac{N_{i+1,j,k}^{l+1/2}-N_{i-1,j,k}^{l+1/2}}{2\Delta r}\right)+\frac{D_{rr,k}}{2}\left(\frac{N_{i+1,j,k}^{l+1/2}-2N_{i,j,k}^{l+1/2}+N_{i-1,j,k}^{l+1/2}}{\Delta r^2}\right)\\ & +\frac{D_{zz,k}}{2}\left(\frac{N_{i,j+1,k}^{l+1}-2N_{i,j,k}^{l+1}+N_{i,j-1,k}^{l+1}}{\Delta z^2}\right)\end{array}$$

Full text

Introduction- Recent observation of the High-Altitude Water Cherenkov Observatory (HAWC) has revealed a TeV gamma-ray halo around the Geminga pulsar, with a spatial extension of $\gtrsim 30$ pc (Abeysekara et al., 2017). The TeV emission is believed to arise from cosmic-ray electrons/positrons (hereafter we do not distinguish positrons from electrons unless specified) injected from the pulsar wind nebula (PWN), via inverse-Compton (IC) scattering off cosmic microwave background (CMB) photons. The detection of such a diffuse TeV emission has been interpreted as the presence of a slow diffusion zone around the pulsar (Abeysekara et al., 2017; Fang et al., 2018; Profumo et al., 2018; Xi et al., 2018) in the framework of 1D isotropic diffusion. Under the same framework, Ref.(Liu et al., 2019) studied the X-ray observation by XMM-Newton and Chandra in the vicinity of Geminga and an upper limit of $5\times 10^{-15}\rm erg~{}cm^{-2}s^{-1}$ in $0.7-1.3\,$keV has been obtained for a region within $600^{\prime\prime}$ around the pulsar. This translates to an upper limit for the magnetic field strength in the TeV halo, i.e., $\leq 0.8\mu$G, which is significantly weaker than the typical interstellar medium (ISM) magnetic field. Furthermore, the combination of a small diffusion coefficient and a weak magnetic field would imply the saturation of the turbulence ($\delta B_{g}/B\simeq 1$ where $\delta B_{g}$ is the fluctuation amplitude of magnetic field at the gyro-scale of particles and $B$ is the mean magnetic field) and the Bohm limit of diffusion, which is, however, very difficult to achieve. For interstellar turbulence, the energy injection scale is $\sim 100$pc and much larger than the gyroscale of the TeV-emitting electrons, where the resonant scattering happens. It is unlikely that $\delta B_{g}/B$ approaches to unity. A plausible scenario is the small-scale waves generated by instabilities. The electron flux at $~{}100$TeV is, nonetheless, too small to generate strong enough streaming instability to overcome Landau damping Kulsrud and Pearce (1969) as well as damping by the background turbulence Farmer and Goldreich (2004); Yan and Lazarian (2004); Lazarian (2016); Fang et al. (2019).

The magnetic field in ISM generally has a mean direction within one coherent length, which is typically $\sim 50-100\,$pc Cho and Ryu (2009); Chepurnov and Lazarian (2010); Beck et al. (2016) and comparable to the size of the TeV halo. 1D particle diffusion actually cannot hold in this scenario, since particles diffuse faster along the mean magnetic field than they diffuse perpendicular to the mean magnetic field in the case of sub-Alfvénic turbulence (Yan and Lazarian, 2008; Giacinti and Sigl, 2012; Nava and Gabici, 2013; López-Coto and Giacinti, 2018). Due to the anisotropy of turbulence in this case, the perpendicular diffusion coefficient is given by $D_{\perp}=D_{\|}M_{A}^{4}$ (Yan and Lazarian, 2008; Xu and Yan, 2013), where $D_{\|}$ is the diffusion coefficient parallel to the magnetic field, $M_{A}\equiv\delta B_{\rm inj}/B$ is the Alfvénic Mach number, which is not far from unity for ISM (i.e., $M_{A}>0.1$), and $\delta B_{\rm inj}$ is the magnetic perturbation at the injection scale of magnetohydrodynamic (MHD) turbulence or coherence length of magnetic field. Also, the synchrotron radiation intensity becomes anisotropic. Electrons that move along the magnetic field will radiate much less efficiently than those move perpendicular to the magnetic field. Therefore, if the mean magnetic field in the vicinity of Geminga has small inclination toward our line of sight (LOS), the observed synchrotron radiation flux would be much reduced compared to that with the assumption of an isotropic magnetic field, while the diffusion perpendicular to the LOS is slow as suggested by the TeV observation. Besides, the small inclination is also beneficial to reproduce the isotropic morphology of Geminga’s TeV halo (López-Coto and Giacinti, 2018).

In this letter, we show that both X-ray and TeV observations can be explained with typical conditions for ISM, such as the magnetic field, the diffusion coefficient and the field perturbation level, by considering anisotropic particle diffusion which is a natural outcome in the presence of sub-Alfvénic turbulence. We will see that the viewing angle plays an important role in determining the observation signals.

Method- The temporal evolution of particle number density in space and energy space $N$ is governed by the transport equation

[TABLE]

where $\dot{E}_{e}$ is the cooling rate of electrons due to synchrotron radiation in ISM magnetic field which is assumed to be $B=3\mu$G in this work, and IC radiation **in the interstellar radiation field with considering the Klein-Nishina effect **. Following (Abeysekara et al., 2017), in addition to CMB, we also consider an infrared photon field (with temperature 20 K and energy density 0.3 $\rm eVcm^{-3}$), and an optical photon field (with temperature 5000 K and energy density 0.3 $\rm eVcm^{-3}$). $Q$ is the source term depicting the electron injection from the pulsar. $\bm{\mathcal{D}}$ is the diffusion tensor. For simplicity, we solve the equation in the cylinder coordinate, defining the $z$-axis to be the direction of the mean magnetic field and the pulsar location to be the origin. By further assuming the system to be symmetric with respect to the $z$-axis (i.e., $\partial/\partial\theta=0$), we can write the transport equation into

[TABLE]

where the diffusion coefficient parallel to the mean magnetic field and perpendicular to it can be set, respectively, by

[TABLE]

Here we neglect the drift effect which would cause asymmetric diffusion, and take $D_{\|}$ to be the typical ISM diffusion coefficient throughout the work which is $D_{0}=3.8\times 10^{28}\rm cm^{2}s^{-1}$ and $q=1/3$ Trotta et al. (2011). Based on our motivation in this study, we will only look into the case with $M_{A}=0.1,0.2,0.3$ respectively, since a larger $M_{A}$ would result in a less anisotropic magnetic field topology. The rightmost term in Eq. 13 consists of three parts: the Dirac functions $\delta(r)$ and $\delta(z)$ specify the injection location, $Q(E_{e},t)$ represents the injection spectrum of electron, and $S(t)$ shows the temporal behavior of the injection rate. More specifically, we assume the injection spectrum of electron to follow a power-law distribution $Q(E_{e})=N_{0}E_{e}^{-p}e^{E_{e}/E_{\rm max}}$, starting from 1 GeV. Here, $N_{0}$ is the normalization constant, $p$ is the spectral index and $E_{\rm max}$ is the high-energy cutoff energy in the spectrum. (Xi et al., 2018) reported a null detection of the diffuse multi-GeV emission from the vicinity of Geminga by Fermi-LAT, suggesting a hard injection spectrum of electron from the PWN. We fix the value of $p$ to be 1.6 in this work, noting that the value of $p$ is actually not important to the predicted X-ray flux as long as the TeV observation is reproduced. $E_{\rm max}$ is assumed to be $200\,$TeV to produce a proper spectral shape measured by HAWC under the hard injection spectrum. Assuming the pulsar to be a pure dipole radiator with a braking index of 3, we have $S(t)=(1+t/\tau)^{-2}$ with $\tau=12\,$kyr being the spin-down timescale of the pulsar. The value of $N_{0}$ is then determined by $\int\int S(t)Q(E_{e})dE_{e}dt=W_{e}$, i.e., the total injected energy in CRe. The evolution of the differential electron density $N(E_{e},r,z,t)$ is solved by a finite difference method (see Supplement for details).

Next, we calculate the emissivity of electrons in the Cartesian coordinate system. Again, we put the pulsar at the origin ($x_{P}=0,y_{P}=0,z_{P}=0$) and define the direction of the mean magnetic field as the $z$-axis. We define the $x-$axis so that the line connecting the pulsar and the observer, i.e., $\overline{PO}$, is in the $xz$ plane (see Fig. 1 for a sketch). Then we envisage a random point $E$ in the space and denote the distance between the point to the observer (i.e., the length of $\overline{EO}$) by $l$, and denote the angle between the line $\overline{EO}$ and the line $\overline{PO}$ by $\theta$. Now let us further consider a circle perpendicular to both the $xz$ plane and the line $\overline{PO}$, with its center, denoted by $C$, attaching to the line $\overline{PO}$ and with the point $E$ on the ring. The two intersection points of the circle and the $xz$ plane are called point $A$ and $B$ respectively, and we call the angle between the line $\overline{AC}$ and the line $\overline{CE}$ angle $\zeta$. The coordinates of point $E$ can then be given by

[TABLE]

The electron density at an arbitrary point $E$ can be found by interpolation given $z_{E}$ and $r_{E}=\sqrt{x_{E}^{2}+y_{E}^{2}}$, based on the obtained electron density distribution $N(E,z,r)$. The number of electron in the element volume in the neighbourhood of point $E$ can then be given by

[TABLE]

where the quantity in the bracket represents the element volume around point $E$. Note that, although the angular distribution of electron is still isotropic (i.e., $dN/d\Omega=N(E,z,r)/4\pi$ since the mean scattering time of an electron $\sim D/c^{2}$ is much shorter than the cooling timescale in the energy range of interest in this work), the synchrotron radiation is anisotropic given a mean orientation of magnetic field considered in this work. Due to the relativistic beaming effect, we can only receive the radiation of electrons moving towards us and the radiation power highly depends on the pitch angle $\alpha$ with respect to the local magnetic field. The latter one is determined by the viewing angle $\phi$ and the position of the point $E$. If there is no magnetic field perturbation, the pitch angle can be given by $\cos\alpha_{0}\equiv\overrightarrow{EO}\cdot\vec{u}_{z}/\overline{EO}=(z_{O}-z_{E})/l=\cos\theta\cos\phi-\sin\theta\sin\phi\cos\zeta$, where $\vec{u}_{z}$ is the unit vector along $z$-axis. In the presence of perturbation, the local magnetic field direction will deviate from the $z$-axis by an angle $\delta$. The average cosine of the pitch angle then becomes $\cos\alpha=\cos\alpha_{0}\cos\delta$. The distribution of $\cos^{2}\delta$, i.e., $f(\cos^{2}\delta)$, where $\int f(\cos^{2}\delta)d\cos^{2}\delta=1$, is obtained from MHD simulations for different $M_{A}$ 111see Supplemental Material [url] for brief description, which includes Ref.[21-26].. The average value of $\cos^{2}\delta$ is close to unity. The flux of synchrotron radiation by electrons in the element volume in the neighbourhood of point $E$ can then be given by

[TABLE]

The IC radiation is isotropic and flux can be given by

[TABLE]

where $n_{\rm ph}$ is the differential density of the background photon field.

For observers at Earth, radiation of any electrons in LOS adds up and is projected onto the celestial sphere. The intensity at any given direction depicted by $\theta$ and $\zeta$ can then be found by $I(\epsilon,\theta,\zeta)=\int dF/\sin\theta d\theta d\zeta$. More specifically, the intensity of synchrotron radiation and the intensity of IC radiation can be given by

[TABLE]

and

[TABLE]

respectively, where $n_{\rm ph}$ is the photon number density of the background radiation. The total flux within certain angle $\theta_{0}$ from the pulsar can be obtained by $F(\epsilon,\theta<\theta_{0})=\int_{0}^{\theta_{0}}\int_{0}^{2\pi}I(\epsilon,\theta,\zeta)\sin\theta d\theta d\zeta$, where $I=I_{\rm syn}+I_{\rm IC}$.

Result- We firstly show the predicted 8-40 TeV gamma-ray morphology for different Alfvénic Mach number $M_{A}$ and different viewing angle $\phi$ in Fig. 2. The Geminga pulsar is located at the center of each panel or the coordinate (0,0). The horizontal axis is parallel to the line $\overline{AB}$ while the vertical axis is parallel to $y-$axis in Fig. 1. The projected distance is calculated based on a nominal distance of 250 pc for Geminga. We can see that the morphology is too compact in the case of $\phi=0^{\circ},M_{A}=0.1$. This is because the perpendicular diffusion coefficient is only $D_{\perp}=3.8\times 10^{24}(E_{e}/1\rm GeV)^{1/3}\rm cm^{2}s^{-1}$ for $M_{A}=0.1$, and the perpendicular diffusion distance is correspondingly only $\sim 5\,$pc within the TeV-emitting electron’s cooling timescale which is $t_{c}\lesssim 10^{12}$s. For a viewing angle of $\phi=0^{\circ}$, such a perpendicular diffusion length is translated to only $\sim 1^{\circ}$ extension in the celestial sphere. The morphology is highly anisotropic in the case of $\phi=5^{\circ},M_{A}=0.1$ which is obviously inconsistent with the observation. This is because the LOS towards the left side of the pulsar (e.g. direction of $\overline{OA}$) passes through more electrons than the LOS towards the right side of the pulsar (e.g., direction of $\overline{OB}$). On the other hand, in both two cases with $M_{A}=0.3$, the morphology does not show a sufficient gradient as that observed by HAWC. We therefore focus on the cases of $\phi=0^{\circ}$ and $\phi=5^{\circ}$ with $M_{A}=0.2$ below. A larger $\phi$ would result in a more anisotropic morphology and a higher X-ray flux so we do not consider it here.

We integrate the TeV emission and X-ray emission over a circular region with a radius of $10^{\circ}$ and $600^{\prime\prime}$ centred at Geminga, respectively, to compare with predicted fluxers with HAWC’s observation and XMM-Newton upper limit. In the calculation, we adjust the value of parameter $W_{e}$ to normalize the predicted TeV flux to the measured one. As we can see from Fig. 3, the predicted X-ray fluxes are lower than the upper limit of XMM-Newton in both two cases. The predicted surface brightness profile (SBP) in $8-40\,$TeV is in good agreement with the observation for $\phi=0^{\circ}$. For $\phi=5^{\circ}$, the predicted SBP is a little flatter than the observation. The reduced Chi-square test returns $\chi^{2}/dof=1.73$ with $dof=14$ being the degrees of freedom in the fitting. It corresponds to a $p$-value of 0.044, implying the fitting is marginally consistent with the data at $2\sigma$ level. We therefore conclude that there should be a magnetic field alignment $\lesssim 5^{\circ}$ with our LOS to explain the observation for $M_{A}\simeq 0.2$. The small inclination between the mean magnetic field and LOS is consistent with the synchrotron polarization measurement by Gao et al. (2010) on the region of $\sim 1^{\circ}$ around Geminga showing very small plane-of-sky magnetic field component. Our model can be tested in the future after many such TeV halos being detected with well determined morphologies, since it would expect (in a flux-limited sample) to observe a large number of elongated systems ($B$ not aligned along LOS) and a minority of roughly spherical ones ($B$ aligned along LOS) provided that $M_{A}$ being significantly smaller than unity. For a larger $M_{A}$, the particle diffusion becomes more isotropic and consequently the morphology of the TeV halo would be less dependent on $\phi$.

Lastly, we note that given the distance from Geminga to Earth being 250 pc, the intervening ISM is expected to contain two or three coherent magnetic fields. The mean magnetic field is unlikely always aligned with our LOS between Geminga and Earth. The pitch angle between the electron moving towards us and the mean field is supposed to be larger outside the TeV halo. The potential increase of synchrotron radiation flux is, however, limited 222Even in the extreme case of a mean magnetic field perpendicular to LOS (or $90^{\circ}$ pitch angle) outside the TeV halo, the synchrotron flux can still be consistent with the upper limit by employing a slightly weaker magnetic field of $2\mu$G for the intervening medium.. Another related issue is the contribution of Geminga to the positron excess measured in many experiments above 10 GeV (Adriani et al., 2009; Ackermann et al., 2012; Aguilar et al., 2014). There is a debate on the contribution of Geminga to the positron excess in the framework of inefficient (isotropic) diffusion of particles(Hooper et al., 2017; Fang et al., 2018; Profumo et al., 2018; Xi et al., 2018; Tang and Piran, 2019; Di Mauro et al., 2019). In the global frame from Geminga to Earth, the injected positrons from Geminga are likely to diffuse fast with the typical ISM diffusion coefficient so that the resulting positron flux at Earth can be enhanced compared to the isotropic diffusion scenario, whereas the exact flux will depend on properties of turbulence between Earth and Geminga.

Solving the transport equation

We employ the operator splitting technique to solve the transport equation

[TABLE]

such that the problem can be simplified into solving a convection equation in energy space and a diffusion equation in real space separately. The Strang splitting scheme is used to decouple the energy operator and the spatial operator, with the flow chart as

[TABLE]

with $l$ being the index of the time step. The spatial operator $\mathcal{L}_{r}$ contain both the $r$-derivative terms and $z$-derivative terms. We further employ the alternating-direction implicit method to divide each time step into two steps of size $\Delta t/2$, and each step can be solved using tridiagonal matrix algorithm.

The implicit second-order upwind scheme is used to discretize the convection equation in energy space as

[TABLE]

where $i$ and $j$ are the indexes of the spatial step in $r$ direction and $z$ direction respectively, while $k$ is the index of the energy step. The above equation can be reduced into

[TABLE]

Given the boundary condition $N_{i,j,k_{\rm max}}=0$ for any $l$, we can solve $N_{i,j,k}^{l+1}$ from $k=k_{\rm max}-1$ to $k=0$. Note that when adopting the Strang splitting scheme, one should replace $\Delta t$ by $\Delta t/2$.

To solve the diffusion equation, we first discretize the equation implementing the Crank-Nicolson scheme, and obtain

[TABLE]

By applying the ADI method, the above equation can be divided into

[TABLE]

and

[TABLE]

For the inner boundary condition where $i=0$ and/or $j=0$, we utilize the spatial symmetry, i.e., $N_{-1,j,k}=N_{1,j,k}$ and $N_{i,-1,k}=N_{i,1,k}$. We impose the CR density to be 0 at the outer boundary which is set to be $r=150\,$pc and/or $z=2.5\,$kpc.

Simulation of the turbulent magnetic field topology in ISM

The synchrotron radiation power highly dependent on the pitch angle of electrons that moving towards the observer with respect to the magnetic field direction. This pitch angle, as derived in the main text, is mainly determined by the observer’s viewing angle, i.e., the angle between the observer’s line of sight and mean magnetic field direction of the interstellar medium (ISM) around the TeV halo, but the turbulence in the magnetic field would cause a small deviation of the local magnetic field direction from the mean magnetic field direction. This would increase the resulting synchrotron flux especially in the case of very small viewing angle. The angle of the deviation (denoted by $\delta$) were obtained from MHD turbulence simulations with different Alfvénic Mach number $M_{A}$ [21-25]. Two state-of-the-art MHD codes (Pluto, http://plutocode.ph.unito.it and Pencil, http://pencil-code.nordita.org) for cross check purpose. The generated magnetic field topologies by these two codes are statistically identical in terms of the distribution of $\cos^{2}\delta$, which is the relevant quantity in calculating the synchrotron radiation. We show the distribution of $\cos^{2}\delta$ for $M_{A}=0.1,0.2,0.3$ in Fig.4. Another cross check for the generated $\theta$ comes from an analytic estimation of the average $\cos^{2}\theta$ [26], which reads $\overline{\cos^{2}\delta}\simeq 1/(1+M_{A}^{2})$. The average value of $\cos^{2}\delta$ from simulation and analytical estimation are also labelled in Fig. 1. We can see that the differences between the simulation and the analytical estimations for all three $M_{A}$ are less than $1\%$.