Constraining the magnetohydrodynamic turbulence around Geminga by observing the γ-ray halo beyond 100 TeV

TL;DR

Scientists observed high-energy gamma rays from a pulsar wind nebula to study cosmic ray acceleration and magnetic turbulence in space.

Contribution

The study presents the first high-precision observation of a gamma-ray halo above 100 TeV and constrains subparsec-scale MHD turbulence.

Findings

Electron/positron flux suppresses exponentially beyond 100 TeV.

The energy dependence of the diffusion coefficient aligns with Kolmogorov turbulence.

Magnetic turbulence at scales smaller than 1 parsec was measured.

Abstract

The extended γ-ray halo surrounding an old pulsar wind nebula (PWN) is an ideal place to investigate both the acceleration and diffusion processes of cosmic rays. In this work, we report the γ-ray halo observed with high precision by the Tibet ASγ experiment at the world’s highest-energy above 100 tera–electron volt (TeV). We determined the acceleration limit at about 100 TeV, beyond which the electron/positron flux suppresses exponentially, while we measured the morphology of the γ-ray halo at three different energies. We also found that the energy dependence of the diffusion coefficient is consistent with the Kolmogorov magnetohydrodynamic (MHD) turbulence and determined the turbulence property at scales smaller than 1 parsec. The measurements indicate that the energy of electrons/positrons is insufficient to amplify the strong MHD turbulence around itself, thereby suppressing the diffusion coefficient.

The Tibet ASγ experiment measured the subparsec-scale MHD turbulence in the Geminga γ-ray halo above 100 TeV.

Figures

Energy spectrum of the Geminga halo γ rays measured by the Tibet ASγ experiment.The blue curve indicates the best-fit with the ECPL model, whereas the light blue band indicates the 1σ confidence interval of the fit. The spectrum is derived by fitting to the data below 158 TeV, whereas the red arrow at 251 TeV represents the flux upper limit at 90% C.L. The green points are the result by HAWC (15), and the gray point is the result by Milagro (23). It can be seen that the HAWC measurements (15) agree with this work below 50 TeV, whereas the Tibet ASγ experiment demonstrates superior sensitivity around 100 TeV. This is primarily due to its enhanced ability to reject high-energy CR backgrounds, enabled by a 2.5-m-deep underground MD array. Alternatively, there is a possibility of some systematic error between the Tibet ASγ and HAWC experiments.

Two-dimensional significance maps of γ-ray emission from Geminga in three energy bins.(A) Each map is smoothed by a circular window with a radius of 3°. The green triangle in each map indicates the location of the Geminga pulsar, whereas the green X mark indicates the location of an another pulsar PSR B0656+14. (B) Data points display the observed radial profile of the surface brightness of halo γ rays, whereas blue curves show the best-fit results with the ECPL injection model in three energy bins. The diffusion coefficient in each energy bin is obtained individually. deg, degrees.

Energy dependence of the diffusion coefficient measured by the Tibet ASγ experiment.Red points display the diffusion coefficient obtained by fitting to the observed radial profile of halo γ rays in each of three energy bins shown in Fig. 2B. The electron energy of each point is calculated as the weighted average energy of electrons producing γ rays in each of three energy bins, i.e., 15.8 TeV < E < 39.8 TeV, 39.8 TeV < E < 100 TeV, and 100 TeV < E < 251 TeV, respectively. The green shaded area represents the 1σ confidence interval, evaluated by taking account of uncertainties in D100 and δ. The light blue dashed line is the best-fit curve obtained by assuming δ = 1/3. The diffusion coefficient measured by HAWC (15) (blue point) is obviously lower than the result in this work, which may be due to differences in CR background estimation. The extrapolation from the average diffusion coefficient at the scale of whole Galaxy is also shown by the gray dashed line, which is derived from the boron/carbon ratio (36).

The power spectrum of the turbulent magnetic field in the interstellar medium.The power spectrum of magnetic turbulence measured in this work in comparison with those measured by Minter and Spangler (25) (blue solid line) and Han et al. (24) (red solid line). The yellow point is derived from the diffusion coefficient measured by HAWC (15). The turbulence is three-dimensional at a scale below ~4 pc, whereas it gradually turns to be two-dimensional above ~4 pc. For three-dimensional turbulence, the spectral index of −5/3 is expected from the Kolmogorov theory.

Tables

Table 1.. Fitting results to the γ-ray spectrum and the SBPs of the Geminga halo.The PL model is disfavored as the energy conversion efficiency required by this model is higher than 100% at 3σ C.L. The cutoff energy representing the acceleration limit of the Geminga PWN is determined to be about 100 TeV. d.o.f., degrees of freedom.

	PL	ECPL
α	2.75 ± 0.16	–*
Ec †	–	98.0−19.5+24.3
logη	2.48 ± 0.77	−0.94 ± 0.13
D100 ‡	1.37−0.47+0.71×1028	1.26−0.42+0.64×1028
χ2/d.o.f.	0.73	0.84

Full text

INTRODUCTION

Cosmic rays (CRs) are high-energy particles from space, mainly composed of protons and atomic nuclei with a small fraction of electrons, γ rays, and neutrinos, as well as fewer antiparticles (1). The antiparticles, such as antiprotons or positrons, are secondary particles produced by the interaction of CR protons/nuclei with the interstellar matter (ISM) during the CR propagation through the Galaxy (2). The charged protons and nuclei propagate diffusively in the Galaxy due to multiple interactions with the disordered magnetic field embedded in the interstellar medium (3, 4). The diffusion picture describes propagation of CRs successfully. It explains most of the observed CR data on Earth, including the precise measurements of secondary particle spectra (5, 6), except for a puzzling positron excess recently reported (7, 8). It is found that the observed positron flux is much higher than that expected from the propagation model above about 10 GeV. This positron excess was measured by PAMELA (7) and confirmed by AMS-02 (8) with higher precision and up to higher energies.

It is generally believed that nearby pulsar wind nebulas (PWNs) are responsible for the extra high-energy positrons (9–11). A PWN contains a central rapidly spinning neutron star surrounded by extremely strong magnetic field, which blows out electron- and positron-winds interacting with the ISM and forming shock waves at the front (12). The injected electrons and positrons, which are further accelerated by the shock wave to very high energies, propagate diffusively to Earth and contribute to the local positron flux. However, a recent observation of TeV γ-ray halos around two nearby PWNs by the High-Altitude Water Cherenkov Observatory (HAWC) experiment made things more complicated. By analyzing the morphology of γ-ray halos, the HAWC experiment found that the electron/positron diffusion rate is too small to interpret the positron excess observed at Earth in terms of the contribution from nearby PWNs (13). Later, the High Energy Stereoscopic System (HESS) also measured the γ-ray spectrum from the halo center (14). Recently, HAWC updated the observation of the γ-ray halos with much higher precision (15).

How charged CRs propagate in the ISM is an essential aspect of astrophysics. A dominant interaction that CRs experience during their propagation is the pitch-angle scattering by the random magnetic field in the Galaxy (4). The interaction between CRs and magnetohydrodynamic (MHD) turbulence leads to diffusive propagation of CRs, which is characterized by the diffusion coefficient D. The widely accepted diffusion coefficient in the Galaxy is, however, constrained only globally based on the secondary nuclei abundance in CRs (16). The morphology of γ-ray halos around a point source, like a PWN, provides a unique measure of the local diffusion coefficient and, thereafter, the property of local MHD turbulence. The energy dependence of diffusion coefficient is especially useful for investigating the power spectrum of MHD turbulence, such as one predicted by the Kolmogorov theory (17). The Tibet ASγ experiment recently observed the γ-ray halo surrounding Geminga, a nearby pulsar at about 250 pc away from Earth (18). The morphology and energy spectrum of the Geminga γ-ray halo are measured from about 10 TeV up to beyond 100 TeV. Here, we derive the acceleration limit at the Geminga PWN. We also present the direct measurements of the CR diffusion coefficient as a function of energy and the property of MHD turbulence at tiny scales in the Galactic disk.

RESULTS AND DISCUSSION

The electron/positron injection spectrum dominantly determines the energy spectrum of halo γ rays and is not very sensitive to the assumption of the energy dependence of diffusion coefficient D(E). Therefore, the electron injection spectrum is determined at first by assuming the diffusion coefficient expected from the Kolmogorov MHD turbulence, i.e., $[eqn]$ . We then get the electron spectrum by best-fitting this model to the observed radial distribution and energy spectrum of γ rays simultaneously.

The energy spectrum of halo γ rays measured by the Tibet ASγ experiment is very soft (see Fig. 1). It could be fitted by either a soft or hard power-law (PL) spectrum with a high-energy cutoff. We therefore perform best-fit calculations with two different forms of the injection spectrum: a PL and a PL with an exponential cutoff (ECPL), $[eqn]$ , where E is the energy of an electron or positron, α is the PL index, $[eqn]$ is the cutoff energy, and $[eqn]$ is the number of electrons/positrons at 1 TeV, far below $[eqn]$ . We use the best-fit $[eqn]$ to estimate η, the conversion efficiency of pulsar’s spin-down energy dissipation to electron/positron energy. In the best-fitting with ECPL, the soft spectrum observed in a limited energy range results in low $[eqn]$ , and the PL index α cannot be well constrained. In this case, therefore, we set α = 1 following the indication of x-ray measurements of the Geminga PWN (19–21).

Energy spectrum of the Geminga halo γ rays measured by the Tibet ASγ experiment.The blue curve indicates the best-fit with the ECPL model, whereas the light blue band indicates the 1σ confidence interval of the fit. The spectrum is derived by fitting to the data below 158 TeV, whereas the red arrow at 251 TeV represents the flux upper limit at 90% C.L. The green points are the result by HAWC (15), and the gray point is the result by Milagro (23). It can be seen that the HAWC measurements (15) agree with this work below 50 TeV, whereas the Tibet ASγ experiment demonstrates superior sensitivity around 100 TeV. This is primarily due to its enhanced ability to reject high-energy CR backgrounds, enabled by a 2.5-m-deep underground MD array. Alternatively, there is a possibility of some systematic error between the Tibet ASγ and HAWC experiments.

We present the fitting results in Table 1. The PL model requires a soft spectrum with α ≈ 2.75, which results in a very large conversion efficiency η ≫ 1 and η < 1 is excluded at 3σ confidence level (C.L.). The PL model is strongly disfavored because η cannot exceed 1.0 unless there is an unknown energy source other than the Geminga pulsar itself, and thus the PL model is strongly disfavored. As for the ECPL model, on the other hand, the conversion efficiency is quite reasonable at η ∼ 0.1 (10%). The hard electron/positron spectrum with =1.0 below $[eqn]$ in the ECPL model is also supported by Fermi-LAT data (22), which set an upper limit on GeV γ-ray flux from the Geminga halo indicating a hard spectrum with α < 1.7. The cutoff energy in the ECPL model is $[eqn]$ TeV, which gives the acceleration limit of the Geminga PWN. The Geminga PWN acceleration limit is determined from the γ-ray observation, thanks to the very high sensitivity of γ rays beyond 100 TeV of the Tibet ASγ experiment. As seen in Fig. 1, the γ-ray spectrum of the Geminga halo observed by the Tibet ASγ experiment is well fitted by the ECPL model. The HAWC measurement of the γ-ray spectrum (15) and the Milagro result (23) are also shown in the plot. It can be seen that the HAWC measurements (15) agree with this work below 50 TeV, whereas the Tibet ASγ experiment demonstrates superior sensitivity around 100 TeV. This enhanced performance stems from its superior background rejection capability for CRs at high energies, afforded by a 2.5-m-deep underground muon detector (MD) array. Alternatively, there is a possibility of some systematic error between the Tibet ASγ and HAWC experiments.

After getting electron/positron spectrum, we derive the diffusion coefficient D(E), from best-fit calculations to the radial profiles of halo γ rays in three energy bins, as shown in Fig. 2B. The best-fit D(E) at the three energy bins are shown in Fig. 3. The energy dependence of diffusion coefficient is determined by fitting a PL form of $[eqn]$ with two free parameters, $[eqn]$ and δ. In Fig. 3, a 1σ confidence interval of the fit is indicated by the light green shaded area. The corresponding 1σ confidence interval of δ is 1.15 ± 0.55, whereas the possibility distribution of $[eqn]$ is integrated. The green dashed line in Fig. 3 is the best fit with δ = 1/3, which is expected from the Kolmogorov turbulence.

Two-dimensional significance maps of γ-ray emission from Geminga in three energy bins.(A) Each map is smoothed by a circular window with a radius of 3°. The green triangle in each map indicates the location of the Geminga pulsar, whereas the green X mark indicates the location of an another pulsar PSR B0656+14. (B) Data points display the observed radial profile of the surface brightness of halo γ rays, whereas blue curves show the best-fit results with the ECPL injection model in three energy bins. The diffusion coefficient in each energy bin is obtained individually. deg, degrees.

Energy dependence of the diffusion coefficient measured by the Tibet ASγ experiment.Red points display the diffusion coefficient obtained by fitting to the observed radial profile of halo γ rays in each of three energy bins shown in Fig. 2B. The electron energy of each point is calculated as the weighted average energy of electrons producing γ rays in each of three energy bins, i.e., 15.8 TeV < E < 39.8 TeV, 39.8 TeV < E < 100 TeV, and 100 TeV < E < 251 TeV, respectively. The green shaded area represents the 1σ confidence interval, evaluated by taking account of uncertainties in D100 and δ. The light blue dashed line is the best-fit curve obtained by assuming δ = 1/3. The diffusion coefficient measured by HAWC (15) (blue point) is obviously lower than the result in this work, which may be due to differences in CR background estimation. The extrapolation from the average diffusion coefficient at the scale of whole Galaxy is also shown by the gray dashed line, which is derived from the boron/carbon ratio (36).

As seen in Fig. 3, the energy dependence of the diffusion coefficient obtained in this work is consistent with the expectation of Kolmogorov’s theory after considering the errors. This is the direct measurement of the energy dependence of the diffusion coefficient using the Tibet ASγ experimental data in an astrophysical environment. As actual tracers of high-energy electrons/positrons, the observed γ rays tell us how electrons/positrons diffuse away from the central PWN. It should be noted that the shaded region in Fig. 3 reflects the 68% confidence interval of the joint distribution of the two parameters, δ and $[eqn]$ . It is usually wider than the marginal distribution for a single parameter δ.

Furthermore, we derive from the obtained D(E) the power spectrum of MHD turbulence, which is closely related to D(E) as expected from the resonant pitch angle scattering of electrons/protons by the turbulence. In Fig. 4, we summarize the magnetic turbulence power spectra derived by different methods. Han et al. (24) obtained an average power spectrum over one-third of the Galactic disk by using the rotation measure and dispersion measure of Galactic pulsars. Minter and Spangler (25), on the other hand, derived the magnetic turbulence over a galactic anticenter sky including the Galactic halo. In Fig. 4, the theoretical extrapolation of the result by Han et al. (24) to a smaller scale is compared with the result derived in the present work. The power spectrum obtained in the present work below the ~4-pc scale is extrapolated with the spectral index of −5/3 expected from Kolmogorov turbulence (black dashed line). The spectral index above the ~4-pc scale is uncertain as the turbulence gradually turns to be two-dimensional with the spectral index lying between −5/3 and −2/3 (black dotted line) (25). As seen in Fig. 4, the spectrum derived in this work is well consistent with the theoretical extrapolation from the spectrum by Han et al. (24). It is evident in this figure, on the other hand, that the magnetic turbulence power measured by Minter and Spangler (25) in the Galactic halo is much weaker than the average values by Han et al. (24) at the Galactic disk scale and that by the present work at a much smaller scale. The amplitude of magnetic turbulence in the Galactic disk may be much larger than in the Galactic halo. Figure 4 also implies that the MHD turbulence environment around Geminga may be common in the Galactic disk. We reveal the power spectrum of MHD turbulence at such a small scale in the Galactic disk.

The power spectrum of the turbulent magnetic field in the interstellar medium.The power spectrum of magnetic turbulence measured in this work in comparison with those measured by Minter and Spangler (25) (blue solid line) and Han et al. (24) (red solid line). The yellow point is derived from the diffusion coefficient measured by HAWC (15). The turbulence is three-dimensional at a scale below ~4 pc, whereas it gradually turns to be two-dimensional above ~4 pc. For three-dimensional turbulence, the spectral index of −5/3 is expected from the Kolmogorov theory.

We lastly add that the present result has another important implication about the origin of strong MHD turbulence and small diffusion coefficient. A straightforward explanation of such a strong turbulence might be that electrons and positrons accelerated and released from the PWN stimulate MHD turbulence in the region around the PWN (26, 27). Large flux of electron-positron pairs in the vicinity of the PWN may amplify the resonant magnetic waves, and the amplified waves may, in turn, confine CR particles around the PWN. As seen in Fig. 1, however, the electron/positron flux above 100 TeV is obviously suppressed whereas the diffusion coefficient at E > 100 TeV is consistent with the behavior of Kolmogorov-type turbulence. We calculated the energy density of magnetic turbulence at the highest energy bin and found that the electron/positron energy density in this energy bin is about two orders of magnitude smaller than the magnetic turbulence energy. Therefore, external energy supplies are needed to produce enough turbulence to suppress diffusion. A possible energy supply is expected from the shock wave in its progenitor supernova remnant, as suggested by Fang et al. (28).

On the basis of the above results, the Tibet ASγ experiment has made a beautiful measurement of a γ-ray halo surrounding the Geminga PWN in an energy region up to the highest energy above 100 TeV. The energy spectrum of electrons/positrons injected from the PWN and the diffusion coefficient at different energies are determined simultaneously. The acceleration limit for the Geminga PWN is determined in this work. This measurement also indicates that MHD turbulence in the region around the Geminga PWN is consistent with that of Kolmogorov type. This work presents the direct measurement of the CR diffusion coefficient as a function of energy and the property of MHD turbulence at tiny scales in the Galactic disk. The turbulence property obtained in this work is consistent with the extrapolation from measurements of the Galactic disk turbulence at larger scales. The Galactic disk seems to have a stronger magnetic turbulence than that in the Galactic halo. The present work also shows that the origin of strong turbulence around the Geminga PWN might be an environmental effect.

MATERIALS AND METHODS

Experiment

The Tibet air shower (AS) array has operated since 1990 in Yangbajing (90.522°E, 30.102°N; 4300 m above sea level) in Tibet, China and observing CRs and γ rays above TeV energies. The Tibet ASγ experiment consists of a surface AS array with an underground MD array (29). We equip the AS array with 597 plastic scintillation detectors, each covering 0.5 m^2^, achieving a total detection area of 65,700 m^2^. We design the MD array to span 3400 m^2^, using water-Cherenkov detectors to capture Cherenkov light from shower muons passing through the water layer. We construct each MD detector with water pools partitioned by concrete and bury it 2.4 m underground. The combined soil overburden and concrete shielding equivalent to ~19 radiation lengths, which restricts access to the water layer to only those muons with energies exceeding 1 GeV (for vertical incidence). By using the MD array, we enhance the experiment’s sensitivity to 100-TeV gamma rays by over an order of magnitude, by rejecting more than 99.9% of hadronic CRs while keeping 90% of gamma rays around the 100-TeV energy range (29, 30). The angular resolutions of arrival direction (50% containment) is estimated to be ~0.5° and 0.2° for 10- and 100-TeV primary γ rays, respectively (29). We reconstruct the primary particle’s energy using measured particle densities from AS detectors. For energies above 10 TeV, we fit the lateral density distribution with the Nishimura-Kamata-Greisen (NKG) function and derive the particle density at 50 m from the shower axis (S50). We then determine the reconstructed gamma-ray energy (Erec) as a function of S50 and the zenith angle. We estimate the energy resolution to be about 40% at 10 TeV, improving to 20% at 100 TeV for primary gamma rays, respectively (31).

In this work, we use the data collected by the Tibet AS array combined with the underground MD array during 719 live days between February 2014 and May 2017. Our data selection criteria are the same as those described in our previous paper (29). For more details including our muon-cut threshold, please see our previous paper (29).

Data analysis

To estimate the background contribution from CRs when analyzing γ-ray emission from Geminga pulsar ( $[eqn]$ , $[eqn]$ ), we use the “equi-declination method” in which an event observed at the same declination as Geminga but at a different right ascension is regarded as a background event. In a declination band including Geminga, we first set adjacent 11 circular windows, each with an angular radius of 10°, and define the middle window centered at Geminga as the “ON-source” window and the remaining 10 windows on both sides of Geminga as the “OFF-source” windows. We then divide each window into the inner window with 3° radius and the outer ring (the outer window) outside the inner window. In this method, we assume that the background contribution in the inner “ON-source” window is the same as that in the outer “ON-source” window. We estimate the anisotropy contribution to the outer “ON-source” window by $[eqn]$ , where $[eqn]$ is the event number in the outer “ON-source” window and $[eqn]$ is the average background event number in outer “OFF-source” windows given by $[eqn]$ with $[eqn]$ denoting the event number in the ith outer “OFF-source” window. By using this k, we calculate the signal event number $[eqn]$ in the inner “ON-source” window as $[eqn]$ where $[eqn]$ is the event number in the inner “ON-source” window and $[eqn]$ is the average background event number in inner “OFF-source” windows given by $[eqn]$ with $[eqn]$ denoting the event number in the ith inner “OFF-source” window.

We divide datasets into the following three energy bins above 10 TeV. Figure 2A shows the significance map of the observed γ-ray emission around the Geminga pulsar in three energy bins, each as a function of the declination and right ascension. Each map is smoothed with a circular window with an angular radius of 3°. In 15.8 to 39.8 TeV, γ-ray emission is seen around the Geminga pulsar with 4.7σ Li-Ma significance. The total $[eqn]$ and $[eqn]$ in this energy bin are 12,721 and 12,177, respectively. We see a separate ~3σ excess around the PSR B0656+14 (Monogem) (32–34). As energy increases, in 39.8 to 100 TeV, the average significance of γ rays from Geminga increases to 5.3σ with $[eqn]$ = 1594 and $[eqn]$ = 1380. As energy increases further, in 100 to 251 TeV, the average significance decreases to 2.3σ with $[eqn]$ = 168 and $[eqn]$ = 138 due to a rapid decrease in event number. The combined significance of the Geminga halo reaches 6.3σ.

Model construction

Electrons and positrons are accelerated by the PWN and injected into the ISM surrounding the PWN. These particles propagate in the ISM diffusively and produce the γ-ray halo through inverse Compton scattering (ICS) with the background radiation. By solving the diffusive propagation equation (eq. S1 in the Supplementary Materials), we get the current electron/positron density distribution $[eqn]$ in the halo, where E and r are the electron/positron energy and position, respectively. Then, we obtain the electron/positron surface density $[eqn]$ at an angular distance θ from Geminga by integrating $[eqn]$ along the line of sight as

$$[eqn]$$

where $[eqn]$ , and d is the distance between Geminga and Earth. Subsequently, the energy-differential surface brightness profile (SBP) around Geminga can be computed by

$$[eqn]$$

where $[eqn]$ is energy of emitted γ rays, ν is the background photon energy, and $[eqn]$ is the emitted γ-ray energy distribution for incident electrons/positrons with energy E and target photons with energy ν (35).

The energy-integrated SBP in a given energy range of $[eqn]$ , as shown in Fig. 2, is calculated by

$$[eqn]$$

The γ-ray SBP reflects the spatial distribution of electrons/positrons and is sensitive to the diffusion coefficient, indicating that the diffusion coefficient can be determined by fitting the SBP to the data.

The total γ-ray spectrum of the halo, as shown in Fig. 1, is calculated by

$$[eqn]$$

The γ-ray spectrum is predominantly determined by the injection energy spectrum of electrons/positrons. This implies that the spectral data can be used to derive the parameters of the injection spectrum.