Multi-messenger modeling of the Monogem pulsar halo

Abstract

The High-Altitude Water Cherenkov Telescope (HAWC) has detected TeV halos associated with two nearby pulsars/pulsar wind nebulae (PWN) — Geminga and B0656+14. These TeV halos extend up to tens of pc from the central accelerators, indicating that the diffusion of ultrarelativistic electrons and positrons in the interstellar medium has been suppressed by two orders of magnitude. Although Geminga and B0656+14 are at similar distances and in the same field of view, they have distinct histories. Notably, B0656+14 probably still resides within its parent supernova remnant, the Monogem Ring, which can be observed in X-rays. In this work, we perform high-resolution simulations of the propagation and emission of relativistic lepton pairs around B0656+14 using a two-zone diffusion model using the GALPROP numerical code. We compared the predicted inverse-Compton spectrum to the observations made by HAWC and Fermi-LAT and found physically plausible model parameters that resulted in a good fit to the data. Additionally, we estimated the contribution of this TeV-halo to the positron flux observed on Earth and found it to be smaller than 10% of the measured flux. We conclude that future observations of the TeV halo and its synchrotron emission counterpart in radio and X-ray frequencies will be crucial to distinguish between various possible models.

1 Introduction

B0656+14 is a middle-aged pulsar with a characteristic age of approximately 110,000 years. It is located about 288 pc away from Earth, as determined through parallax measurements (Brisken et al., 2003b). The estimated pulsar age and distance is consistent with those of the “Monogem Ring” supernova remnant, a ring-shaped structure visible in X-rays with an apparent extension of 25^∘, implying that they likely share a common origin (Thorsett et al., 2003). We refer to B0656+14 as the “Monogem pulsar,” hereafter. Monogem, first discovered as a radio pulsar (Manchester et al., 1978), has been extensively studied at all wavelengths from radio to gamma rays (e.g., Zharikov et al. 2021; Durant et al. 2011; Schwope, Axel et al. 2022; Abdo et al. 2010). Monogem has a rotation period of $P\approx 385~{}\mathrm{ms}$, and a spin-down power $(\dot{E})$ of $3.8\times 10^{34}$ erg $\mathrm{s}^{-1}$. It has a surface magnetic field $(B_{\text{surf }})$ of $4.65\times 10^{12}\;\mathrm{G}$. Pulsars like Monogem emit highly energetic winds of electrons and positrons escaping along open magnetic-field lines from the magnetosphere. The wind ends in a termination shock at the radius where the ram pressure of the wind equals the pressure of the medium surrounding the pulsar (cf. Fig. 1). These shocks create pulsar wind nebulae (PWNe), which can be observed through the non-thermal radiation emitted by the injected particles. The brightness of PWNe depends on the pulsar’s spin-down power which declines with pulsar age. So older pulsars are typically more challenging to detect. Interestingly, the PWN associated with Monogem is one of the few PWNe that have been identified. It was detected in between $\sim 0.005\text{--}0.2$ pc radius around Monogem using Chandra X-ray data, showing a roughly round morphology (Bîrzan et al., 2016). The X-rays are thought to be caused by synchrotron radiation of $e^{\pm}$ in the downstream region of the termination shock (Bîrzan et al., 2016), which are expected to produce very-high-energy (VHE) counterparts at TeV energies through inverse-Compton (IC) scattering of the ambient photon field (Kargaltsev et al., 2013).

The first evidence of TeV emission around the Monogem pulsar was discovered through Milagro observations (Abdo et al., 2009). Subsequently, more details about the TeV emission were revealed by HAWC with the detection of TeV halos around the PWNe of both the Monogem pulsar and another middle-aged pulsar known as Geminga (Abeysekara et al., 2017). Pulsar halos (Linden et al., 2017) are emission envelopes stretching up to tens of parsecs, which are purported to have been produced by a population of relativistic $e^{\pm}$ that have diffused out of the PWNe into the ISM. So far, these haloes have only been identified in gamma-ray emission above a few TeV. Recent studies (Abeysekara et al., 2017; Profumo et al., 2018; Hooper et al., 2017; Di Mauro et al., 2019; Manconi et al., 2020) suggest that to explain the vast extension of these TeV halos, it is necessary to have a reduction in the diffusion coefficient by approximately two orders of magnitude over distances ranging from tens to hundreds of parsecs around the parent pulsar wind nebulae. An alternative model suggests that the observations can be explained by considering the transition between the quasi-ballistic regime and the diffusion regime in the propagation of electrons and positrons after they leave the acceleration site (Recchia et al., 2021). However, this model has been critiqued for requiring an acceleration efficiency exceeding 100% of the pulsar’s spin-down power (Bao et al., 2022). Therefore, the slower diffusion process of cosmic rays (CRs) in the nearby ISM of the parent PWN remains the most plausible explanation for the TeV halos.

The biggest challenge in understanding TeV halos is to explain why the diffusion of particles with energies in the multi-TeV range is suppressed around the source for hundreds to thousands of years (Recchia et al., 2021; Martin et al., 2022; López-Coto & Giacinti, 2018). Some studies have put forward a well-motivated theory that suggests the origin of magnetic turbulence is environmental, caused by the host supernova remnant or a pre-existing feature in the ISM in which the TeV-halo develops (Fang et al., 2019). On the other hand, independent studies have proposed that the self-induced resonant and non-resonant streaming instability caused by the propagation of the electrons and positrons may lead to strong self-confinement in the vicinity of the sources reaching 50–100 pc from the PWN (Malkov et al., 2013; Nava et al., 2016; Evoli et al., 2018; Nava et al., 2019; Schroer et al., 2021).

These hypotheses may leave different imprints on the halo emission morphology. In the case of electron/positron self-confinement, the diffusion suppression region is expected to evolve dynamically with the pulsar. In the presence of environmental turbulence, the confinement region is expected to share the center of the parent SNR and evolve in conjugation with the expansion of the SNR’s shock wave and its interaction with the pulsar wind. These effects have been demonstrated phonologically by Jóhannesson et al. (2019) and Di Mauro et al. (2019) for the case of the Geminga pulsar. The slow diffusion phenomena around Geminga pulsar and Monogem have been explored analytically by Martin et al. (2022) and Schroer et al. (2023), but both have assumed a static system, leaving the impact of the proper motion of the pulsar and origin of the slow diffusion out of scope.

In this paper, we investigate the propagation of CRs in the Monogem PWN’s surrounding interstellar medium (ISM). We employ a numerical approach on a two-zone diffusion model to achieve this. The propagation of CRs and non-thermal emissions are calculated using the GALPROP numerical code (v57; Porter et al., 2022a). This package enables accurate forecasts of the gamma-ray emission morphology and radiation spectra by using advanced modeling of the interstellar radiation field (Porter et al., 2017), spatial dependent magnetic field, and three-dimensional interstellar gas (Jóhannesson et al., 2018). It also allows a non-equidistant spatial grid, which provides exceptionally high spatial resolution around the astrophysical accelerator. We compare our simulated non-thermal maps with HAWC observations and Fermi-LAT flux upper limits in the 8–40 TeV and 10–1000 GeV energy range, respectively. We explore how the pulsar’s age, injection spectra, slow diffusion zone size, magnetic field, and proper motion of the Monogem pulsar impact the IC and synchrotron emission properties. In the near future, advanced gamma-ray detectors such as the Cherenkov Telescope Array (CTA) could verify our predictions of the spatial shape and spectrum of Monogem, as analyzing Monogem using the conventional ON/OFF method may prove difficult due to the large extension of its emission and the energy-dependent nature of its gamma-ray spatial morphology (Aharonian et al., 2023). In addition, we offer hints about possible future searches of the synchrotron halo emission across radio to X-ray frequencies. Finally, utilizing models that align with the electromagnetic observations across multiple wavelengths, we put forth projections for the Monogem pulsar’s contribution to the positron flux that AMS-2 has measured.

This paper is structured as follows: In section 2, we introduce the properties of Monogem pulsar and the modeling of the injection and transportation of leptonic particles from the pulsar. In section 3, we compare the diffusive synchrotron and IC emission around Monogem to the HAWC and Fermi-LAT data. We also compare the expected positron flux from Monogem to the positron flux measured by AMS-2. The injection and propagation properties of CRs around Monogem pulsar are discussed.

2 Modeling of Non-thermal Emission from the Monogem Pulsar Halo

Figure 1: Illustration of PWN two-zone diffusion model. TS (termination shock) marks the acceleration site of the electron positrons escaped from the PWN. $r_{1}$ is the transition radius, and $r_{2}$ is the radius of the SDZ.

Source	2HWC J0700+143
Associated pulsar	B0656+14 (Monogem pulsar)
Present period $(P)$	384.94 ms
Period derivative ($\dot{P}$)	$5.5\times 10^{-14}\mathrm{~{}s}\mathrm{~{}s}^{-1}$
Initial spin-down power ($\dot{E_{0}}$)	$1.84\times 10^{35}\;\mathrm{erg}\;\mathrm{s}^{-1}$
Characteristic age $\left(\tau_{c}\right)$	$110\;\mathrm{kyr}$
Current age assumed $(T^{\star})$	$99\;\mathrm{kyr}$, $60.5\;\mathrm{kyr}$, $11\;\mathrm{kyr}$
Distance $(d)$	$288\;\mathrm{pc}$
Galactic coordinates ($l,b$)	$(201.1^{\circ},8.3^{\circ})$
Galactocentric coordinates $(X)$	$-8.77\;\mathrm{kpc}$
Galactocentric coordinates $(Y)$	$-0.10\;\mathrm{kpc}$
Galactocentric coordinates $(Z)$	$0.04\;\mathrm{kpc}$
Diffusion coefficient at 4 GV in the ISM $\left(D_{0}\right)$	$4.5\times 10^{28}$ cm2 s-1
Diffusion coefficient at 4 GV inside the bubble $\left(D_{\rm{SDZ}}\right)$	$1.3\times 10^{26}$ cm2 s-1
Diffusion bubble radius $\left(r_{1}\right)$	$30\;\mathrm{pc}$, $50\;\mathrm{pc}$
Diffusion bubble core $\left(r_{2}\right)$	$50\;\mathrm{pc}$, $70\;\mathrm{pc}$
Energy break $\left(E_{\rm b}\right)$	10 GeV
Smoothness $\left(s\right)$	0.5
Low-energy spectral index $\left(\gamma_{0}\right)$	1.0
High-energy spectral index $\left(\gamma_{1}\right)$	1.8–2.2
Spin-down power efficiency $\left(\eta\right)$	$3$–$16\%$

Table 1: Parameters used to model the pulsar and the slow diffusion zone.

2.1 Particle Injection from Monogem

To better understand the extended non-thermal emissions surrounding the Monogem pulsar, we must first determine the injection properties of the pairs responsible for these emissions. According to recent HAWC results, the TeV halo around Monogem is caused by the IC emission of high-energy $e^{\pm}$ particles that have escaped from the PWN and made their way into the ISM. These $e^{\pm}$ particles are accelerated at the PWN, which is powered by the pulsar’s spin-down. Currently, the spin-down power of the pulsar is estimated to be $\dot{E}=3.8\times 10^{34}$ erg s^-1, as reported by Manchester et al. (2005). Pulsar wind material, which consists mainly of $e^{\pm}$, enters the termination shock of the pulsar wind and the particles will obtain a nonthermal energy distribution as the result of diffusive shock acceleration at the shock, or due to acceleration by magnetic-field reconnection processes in the wind itself. Multi-wavelength studies suggest that the $e^{\pm}$ injection spectrum of PWNe follows a broken power-law distribution, with a low-energy index in the range $[1.2,1.7]$, a high-energy index in the range $[2.1,2.8]$, and a spectra break occurring at $\sim 10\text{--}500$ GeV (Bucciantini et al., 2010).

In our simulation, we do not consider the spatial extent of the PWN, because it only covers a radius of approximately 0.2 pc. This is much smaller than the extent of the halo emission we are interested in. Additionally, we assume that electrons and positrons undergo the same acceleration and propagation process. The pairs are injected into the ISM from the pulsar’s location in an isotropic manner, following a smooth, broken power-law spectrum. This same spectrum was used in the case of Geminga discussed by Jóhannesson et al. (2019):

$$\frac{dn_{e^{\pm}}}{dE}\propto E^{-\gamma_{0}}\left[1+\left(\frac{E}{E_{\rm{b}}}\right)^{\frac{\gamma_{1}-\gamma_{0}}{s}}\right]^{-s},$$ (1)

where $n_{e^{\pm}}$ is the number density of the $e^{\pm}$, and $E$ is the particle kinetic energy. The power-law indices $\gamma_{0}$ and $\gamma_{1}$ are the spectral index at low and high energies, respectively. We stress that these are not the same parameters as the broken power-law indices mentioned in the previous paragraph as we have employed a different formalism for the injection spectrum. We set $\gamma_{0}=1.0$, and explored scenarios when $\gamma_{1}=$ 1.8, 2.0, 2.2. We inject electrons from 100 MeV to an energy cutoff of 1 PeV, with a break energy at $E_{\rm{b}}=10$ GeV.

The normalization of the injected spectrum is a free parameter to ensure that the total luminosity of the injected pairs $L_{e^{\pm}}$ at time $t$ remains a fraction of the spin-down power $\dot{E}(t)$ of the pulsar (Pacini & Salvati, 1973):

$$L_{e^{\pm}}(t)=\eta\dot{E_{0}}\left(1+\frac{t}{\tau_{0}}\right)^{-\frac{n+1}{n-1}},$$ (2)

where $\dot{E_{0}}$ is the initial spin-down power of the pulsar, and $\eta$ is the fraction of the pulsar spin-down power that is converted into the luminosity of injected pairs. The $n$ in the index is the braking index of a pulsar. Until now, only a few pulsars have their braking index accurately measured. Previous studies show $n=15$ for Monogem pulsar, but the error bar spans $\sim 120$ in value, leaving the actual value ambiguous (Johnston & Galloway, 1999). We adopt the pulsar braking index corresponding to the magnetic-dipole model, thus $n=3$, which is also shown to yield a good agreement between the pulsar characteristic age and the kinetic age for middle-aged pulsars (Bailes, 1989; Lorimer et al., 1997). The initial spin-down time-scale of the pulsar $\tau_{0}$ is connected with the current characteristic age $\tau_{c}$ of the pulsar by the following expression:

$$\tau_{0}\equiv\frac{P_{0}}{(n-1)\dot{P_{0}}}=\frac{2\tau_{c}}{(n-1)}-t_{\rm age}.$$ (3)

For a braking index of $n=3$, equation (3) becomes $\tau_{0}=\tau_{c}-t_{\rm age}$. It is worth noting that the true age of a pulsar often has a considerable discrepancy with its characteristic age (Brisken et al., 2003a).

One of the key arguments that the Monogem Ring and the Monogem pulsar have the same origin is that the pulsar characteristic age agrees with the estimated SNR age based on Sedov modeling (Thorsett et al., 2003). This only holds if: 1) the characteristic age of Monogem pulsar is a representation of the true age of the pulsar, 2) if the Monogem Ring is indeed in the Sedov phase. Searching for evidence of the pulsar’s true age in non-thermal diffused emissions can provide insight into the origin of the Monogem pulsar. In this study, we investigate how the initial timescale affects the diffuse emission by testing three different scenarios where $\tau_{0}$ is set to 0.1$\tau_{c}$, 0.45$\tau_{c}$ and 0.9$\tau_{c}$. In the first scenario, it is assumed that the actual age of the pulsar is close to its characteristic timescale (which suggests that the pulsar had much higher spin-down power in the past).

We utilized GALPROP v57 to perform our calculations. Our simulations were conducted on a non-equidistant grid with high spatial resolution. The grid size varies based on a non-linear (tangent) transformation, as described in detail in Porter et al. (2022b). This grid setup allows for maximum resolution at a specific point in 3D space, decreasing at greater distances and optimizing simulation memory usage and speed. The grid size at the current location of the Monogem pulsar is 2 pc, while at a reference distance of 700 pc, it increases to 100 pc. The resolution is roughly uniform in the $\sim 60$ pc region from the source, encompassing the entire evolution of the Monogem PWN. To solve the transport equation in the time domain, we used a time step of 50 years, which is sufficient to capture the pairs’ energy loss. We ensured that the number of time steps was adjusted to cover the entire system’s lifetime based on the different true pulsar ages we tested. Our simulations were executed on the Dutch National supercomputer, Snellius, using a single node with 224 GB of memory and 128 CPUs. Each simulation took around $\sim 5$ hours with the computing node operating at total capacity.

The GALPROP coordinate system is right-handed, centering at the Galactic center, with the Sun locating at $(X_{s},\,Y_{s},\,Z_{s})=(8.50,\,0,\,0)$ kpc. The Monogem pulsar is positioned at the Galactic coordinates $(l,\,b)=(201.1^{\circ},\,8.3^{\circ})$, which corresponds to $(X,\,Y,\,Z)=(8.766,\,0.103,\,0.042)$ kpc in GALPROP coordinates. Its proper motion is measured to be $\mu_{\alpha}\cos(\delta)=44.16$ mas/yr, and $\mu_{\delta}=-2.43$ mas/yr (Hobbs et al., 2005). Assuming that the pulsar velocity is constant, the line of sight velocity is zero, and the pulsar age is 99 kyr (corresponds to $\tau_{0}=0.1\tau_{c}$), the Monogem pulsar moved roughly from $(X_{0},\,Y_{0},\,Z_{0})=(8.768,\,0.100,\,0.036)$ kpc with a velocity of $(V,\,U,\,W)=(-17.53,\,24.34,\,52.09)$ km s^-1 over its lifetime, where $V$, $U$, and $W$ are velocities along the $X$, $Y$ and $Z$ axis respectively.

We used the “Sun_ASS_RING_2” Galactic magnetic-field model and the R12 interstellar radiation field model to calculate the synchrotron and inverse Compton emissions. The magnetic-field model is based on Sun, X. H. et al. (2008) which consists of regular and random fields. At the location of the Sun, the regular field is $2\;\mu$G, and the random field is $3\;\mu$G. The magnetic-field strength also has an axisymmetric spiral component with reversals in rings. However, it is not flexible enough to consider the local variance of the magnetic field around the PWN due to magnetic turbulence. To test how the magnetic-field strength affects the non-thermal emission, we compare cases where the random field is $3\;\mu$G, $5\;\mu$G, and $10\;\mu$G at the location of the Sun. It is important to note that even though this modification impacts the magnetic field on a Galactic scale, it is still valuable for estimating the impact on the non-thermal emission around the source. In section 3.2, we will show that the additional synchrotron loss during propagation, caused by an enhanced magnetic field, leads to some suppression in our positron flux estimation at Earth.

2.2 Slow Diffusion Zone

We assume the diffusion coefficient within a confined region around Monogem is lower than in the general ISM. This assumption is based on recent studies (e.g., Profumo et al., 2018; Jóhannesson et al., 2019) of similar systems and the expected increased magnetic-field turbulence within the region surrounding the pulsar. As illustrated in Figure 1, the confined zone around PWN is called the “slower diffusion zone” (SDZ). Additionally, we suppose that the increased turbulence in this zone does not alter the spectra of the particle population. The spatial dependence of the diffusion coefficient is then determined by

$$D=\left(\frac{E}{E_{0}}\right)^{\delta}\times\begin{cases}D_{\rm SDZ},&r<r_{1},\\ D_{\rm SDZ}\left(\frac{D_{\rm ISM}}{D_{\rm SDZ}}\right)^{(r-r_{1})/(r_{2}-r_{1})},&r_{1}\leq r\leq r_{2},\\ D_{\rm{ISM}},&r>r_{2}.\end{cases}$$ (4)

where the diffusion coefficient has a power law energy dependence with slope $\delta$. We used $\delta=0.35$, which leads to a Kolmogorov diffusion. As for the reference energy $E_{0}$, we use 4 GeV following Jóhannesson et al. (2019). The diffusion coefficient is shown in Table 1. The transport of particles within a bubble (of radius of $r_{1}$) of slow diffusion surrounding the pulsar is determined by the coefficient $D_{\rm SDZ}$. The diffusion coefficient gradually changes from $D_{\rm SDZ}$ to $D_{\rm ISM}$ between the transitional radius $r_{1}$ and the outer radius $r_{2}$, where $D_{\rm ISM}$ is the average diffusion coefficient in the general interstellar medium. On top of the standard diffusion of the pairs in propagation, we include the re-acceleration due to energy transfer by scattering on the Alfvén wave of the ISM plasma. We use a velocity of the Alfvén wave $V_{A}=17$ km/s, as inferred from the CR isotope measurement of the primary to secondary species (Boschini et al., 2017). However, we expect the effect of the re-acceleration to be negligible.

Recent research on the self-confinement mechanisms that result from the propagation of charged particles around PWNs suggests that the suppression of diffusion extends from 50 to 100 pc, according to theoretical studies by Malkov et al. (2013); Nava et al. (2016); Evoli et al. (2018); Schroer et al. (2022). There may be a link between the Monogem pulsar and the magnetic disturbance it produces in the surrounding interstellar plasma. Due to this, the size and location of the slow-diffusion bubble may change over time. To account for the expansion of the bubble over the lifetime of the PWN, we estimate that the radius of the bubble grows as the diffusion of the charged particles. Therefore, $r_{1,2}=\mu\sqrt{t}$, where $\mu$ is a constant. To align with the observations made by HAWC, the current values of $r_{1}$ and $r_{2}$ must be within a few tens of parsecs.

Motivated by scenarios involving confinement due to the turbulence of a supernova remnant, we explore the possibility that the Monogem Ring defines the SDZ. In this case, the center of the SDZ stays static at the birth location of the Monogem pulsar. We assume the same expansion history of the slow-diffusion bubble. The value of $\mu$ is chosen to assure the current-day SDZ size is close to approximately 70 pc in radius, consistent with current size estimates of the expanding supernova remnant’s shell (Thorsett et al., 2003).

In order to cover different scenarios, we have set up the following simulation configurations: A) In the first configuration, both the pulsar and the SDZ remain stationary and are centered at the current location of B0656+14. The size of the SDZ remains unchanged over time. B) In the second configuration, both the pulsar and SDZ move with the current proper motion of the pulsar and the size of the SDZ evolves over time. C) In the third configuration, the PWN moves at the pulsar proper motion. However, the SDZ center remains static at the birthplace of the pulsar. The size of the SDZ expands over time and has a radius of 70 pc at the current time. It is clear that the stationary scenario is not a realistic representation of the actual system. However, with such setup, it is possible to for us to identify the effect of the size of the SDZ and the relative motion of the PWN and SDZ on the non-thermal emission.

3 Results

In this section, we explore model parameters by comparing the IC emission with HAWC and Fermi-LAT observations from the Monogem pulsar region. We also discuss how these parameters influence the IC and Synchrotron emission morphology. Our base propagation parameter setup and slow-diffusion-zone model parameters are presented in Table 1. To estimate the pulsar’s injection efficiency ($\eta$) for each model, we fit the surface brightness profile in the 8-40 TeV range to the HAWC observation using a traditional $\chi^{2}$ analysis.

3.1 Non-thermal Radiation

3.1.1 Senario A: Stationary PWN and Slow-diffusion Bubble

Figure 2: Surface brightness of IC emission in 8-40 TeV range for stationary pulsar and SDZ. The initial time-scale of pulsar $\tau_{0}=0.1\tau_{c}$. The magnetic-field strength $B=3\mu$G. The escape efficiency $\eta$ is tuned to fit blue data points for HAWC observation.

We assume that the Monogem pulsar and its slow-diffusion bubble are stationary and located at the Mongem pulsar’s present-day position. To construct non-thermal emission maps, we ran various GALPROP simulations using two different sizes of the SDZ: $(r_{1},r_{2})=(30,50)$ pc and $(r_{1},r_{2})=(50,70)$ pc, as well as three different high-energy injection indices: $\gamma_{1}=1.8,2.0,2.2$. We created radial profiles from the resulting gamma-ray maps and used them to fit the HAWC data points in the energy range of 8–40 TeV. Our analysis, presented in Figure 2, reveals no noticeable variation in the surface brightness with the selected sizes of the SDZ and injection indices.

The time scale ($\tau_{\rm cool}$) for a electron/positron cooling due to inverse Compton emission and synchrotron emission is:

$$\tau_{\rm cool}=\frac{4m_{\rm{e}}c^{2}}{3c\sigma_{T}\gamma}\left(U_{B}+\begin{cases}U_{\rm ph},&\gamma<\gamma_{\rm KN},\\ \frac{U_{\rm ph}}{(1+4\gamma\varepsilon_{0})^{3/2}},&\gamma_{\rm KN}\leq\gamma,\end{cases}\right)^{-1},$$ (5)

where, $\sigma_{T}$ is the Thompson scattering cross section, $U_{B}$ and $U_{\rm{ph}}$ are the energy density of the magnetic field and the photon field respectively, and $\varepsilon_{0}$ is the normalized CMB photon energy. When electrons/positrons have a Lorentz factor $\gamma$ greater than $\gamma_{KN}\equiv\frac{m_{e}c^{2}}{4h\nu_{0}}$, the Klein-Nishina (KN) effect becomes significant and suppresses their cooling rate. Assuming that the average background photon energy $h\nu_{0}$ is the same as that of CMB photons, the KN effect becomes significant when the energy of electrons/positrons is $\gtrsim 140$ TeV. The diffusion length under the slow diffusion regime is determined by the expression $R_{\rm{diff}}=2\sqrt{D_{100~{}{\rm TeV}}\textrm{min}\{\tau_{\rm{cool}},\tau_{\rm{inj}}\}}$, which amounts to approximately 28pc. Since the size of the SDZ is more prominent than this value, electrons/positrons cool down efficiently through IC and synchrotron radiation before diffusing away from the SDZ into the general Interstellar Medium.

As depicted in Figure 2, the required escape efficiency to fit the TeV-halo observation does not differ significantly for SDZ sizes of 30–50 pc and 50–70 pc. If the SDZ size is smaller than the e^± diffusion length, we expect a steeper surface brightness profile than those obtained in Figure 2, since the e^± diffuse away before they can be efficiently cooled. We have also found that for values of $\gamma_{1}$ between 1.8 to 2.2, a fraction of approximately 5% to 15% of the current pulsar spin-down energy must be carried out by the electrons/positrons that escape from the PWN into the ISM to explain the observations made by HAWC. A higher escape efficiency is required for a softer injection spectrum.

Figure 3 shows the energy spectra of IC emissions in the energy range of 10 GeV to 25 TeV from the $10^{\circ}$ region around the pulsar. We tested different parameter combinations and found that SDZ sizes of (30,50) pc and (50,70) pc are consistent with the HAWC observations, but the latter leads to slightly higher flux. For softer injection spectra, the flux falls off steeper as function of distance from the pulsar at TeV energies. When comparing the GeV emission with the Fermi-LAT upper-limit derived by Di Mauro et al. (2019), we found that all the parameter combinations we tested show GeV spectra well below the Fermi-LAT range, except for when the SDZ size is (50,70) pc and the injection index $\gamma_{1}$ is 2.2. In this case, the emission almost touches the Fermi-LAT upper-limit at $\sim 40\text{--}100$ GeV, and any higher $\gamma_{1}$ values are disfavored.

Figure 3: IC spectra of $\gamma_{1}=$1.8, 2.0, 2.2 of stationary pulsar and SDZ, compared with Fermi-LAT upper-limits (Di Mauro et al., 2019) and HAWC observation (Abeysekara et al., 2017). The initial time-scale of pulsar $\tau_{0}=0.1\tau_{c}$. The magnetic-field strength $B=3\mu$G. Top: SDZ size of (30, 50) pc. Bottom: SDZ size of (50, 70) pc.

The synchrotron spectra predicted from 10^∘ region around the PWN is shown in the top panel of Figure 4, together with the corresponding IC emission. Similar to the IC emission, by increasing the size of the SDZ from (30, 50) pc to (50, 70) pc, the synchrotron emission flux is enhanced towards lower energy range (sub-eV to sub-keV). The larger SDZ confines lower energy pairs to cool through emission before diffusing away from the SDZ, thus a flatter spectra profile in the eV–keV range is expected. The emission intensity level and spectra feature are similar at keV energies for SDZ sizes of (30, 50) pc, and (50, 70) pc, as well as for different injection indexes $\gamma_{1}=1.0,2.0$ and $2.2$. If we were going to disentangle the model parameters from observational data, the thermal components in lower than keV energies would make it challenging.

Figure 4: Synchrotron and IC emission spectra of $10^{\circ}$ region around the PWN given different model parameter combinations. The initial time scale of the pulsar is $\tau_{0}=0.1\tau_{c}$ for all cases. Top: multi-wavelength non-thermal emission spectra for the stationary scenario. The magnetic field is fixed at B= 3$\mu$G. Middle: multi-wavelength non-thermal emission spectra for the scenario where pulsar and the SDZ move together. The injection index is $\gamma_{1}=2.0$. Bottom: multi-wavelength non-thermal emission spectra for three scenarios considered. 1) Stationary. 2) Both the pulsar and SDZ move with pulsar proper motion. 3) Only the pulsar moves with the proper motion. The SDZ is fixed. SDZ size is (50,70)pc. The injection index is $\gamma_{1}=2.0$. The magnetic field is $B=3\mu G$. Note: The gap is in between the synchrotron and IC emission is due to the output map setup in the simulation, which coincides with the observational MeV gap.

We present the IC emission for Monogem pulsar assuming varying true ages. We adjusted the parameters $t_{\rm age}$ and $\tau_{0}$ based on Eq. 3 to simulate their emission. Specifically, we compared the emissions from $\tau_{0}$ values of 0.1, 0.45, and 0.9 $\tau_{c}$, which correspond to $t_{\rm age}$ values of 99, 60.5, and 0.11 kyr, respectively. In this simulation, the pulsar is assumed to be stationary. The simulation results are shown in Figure 5, which displays the IC surface brightness and spectra. In the energy range of 8–40 TeV, the surface brightness and the spectra are almost identical for $\tau_{0}=$ 0.1 and 0.45 $\tau_{c}$, but a lower flux is found for $\tau_{0}=\,0.9\,\tau_{c}$. The latter case features a young pulsar at its maximum power output when $t_{\rm{age}}\ll\tau_{0}$. The cooling time of e^± at 100 TeV is 18 kyr if B= 3$\mu$G based on Eq. 5, longer than the injection timescale of the particles. The corresponding diffusion length of these e^± is approximately 2 pc, well below the SDZ size of (50, 70) pc. This means that at such energy, the injection of e^± dominates the system. The SDZ has yet to be saturated with e^± responsible for TeV $\gamma$-rays. Thus it is less bright than the cases with larger pulsar true ages when $\tau_{0}=$ 0.1 and 0.45 $\tau_{c}$. The spectrum at lower energies is enhanced as the pulsar age increases.

3.1.2 Scenario B: Moving PWN and Slow-diffusion Bubble

In this section, we examine the effect of the Monogem pulsar’s proper motion on the anticipated spectra and spatial morphology of the non-thermal radiation. Furthermore, we investigate the possibility of detecting the source of the slow-diffusion bubble around the PWN in the non-thermal emission.

Figure 5: IC emission of varied $\tau_{0}$ of stationary pulsar and SDZ, compared with HAWC observation (Abeysekara et al., 2017) and Fermi-LAT upper-limits (Di Mauro et al., 2019). The size of SDZ is (50,70) pc. The magnetic-field strength is B= 3 $\mu$G. The injection index $\gamma_{1}=2.0$. Top: Surface brightness in 8- 40 TeV for $\tau_{0}=0.1,0.55,0.9$ $\tau_{c}$ fitted to HAWC observation. Bottom: IC spectra of $\tau_{0}=$ 0.1, 0.45, 0.9 $\tau_{c}$ compared to HAWC observation and Fermi-LAT upper-limit.

We begin by discussing the scenario in which the propagation of pairs injected by Monogem self-induces the slow diffusion bubble. In this sense, the SDZ will always be centered at the pulsar location and move with the same proper motion. Figure 6 displays the intensity map of IC and synchrotron emissions near the PWN within a radius of approximately $\sim 30$ pc. The extended halo emission is visible in the GeV-TeV, as well as radio and X-ray range owing to synchrotron emission. As can be seen in this figure, the IC emission at 10 GeV and the synchrotron emission at 100 GHz exhibit roughly spherical emission with the centroid offsetting from the current location of the pulsar. This offset is aligned with the pulsar proper motion direction. Higher energy observations, such as the IC mission at 10 TeV (top-right panel) and synchrotron emission at 5 keV (bottom-rigth panel), show more spherically symmetric profiles at the current pulsar location due to the short cooling time of the TeV electrons/positrons. The TeV-halo around Monogem observed by HAWC has a rather lobbed profile (Abeysekara et al., 2017), suggesting that the asymmetry is not due to the proper motion of the pulsar but possibly caused by a spatial asymmetry in the diffusion coefficient in the ISM, or convection due to a local or Galactic-scale wind. We will discuss the former hypothesis in Section 3.1.3. For the latter case, the expected wind velocity in the environment to modify the emission morphology in the TeV range is greater than 1600 km/s, which is unlikely to be achieved.

(a) IC emission at 10 GeV

(b) IC emission at 10 TeV

(c) Synchrotron emission at 100 GHz

(d) Synchrotron emission at 5 keV

Figure 6: Intensity maps at different wavelengths for the region around Monogem pulsar taking into account the pulsar proper motion. The center of the SDZ moves with the pulsar. The maps are shown for SDZ size of (50,70) pc. The initial time-scale of pulsar $\tau_{0}=0.1\tau_{c}$. Magnetic-field strength $B=3\mu G$, and injection index $\gamma_{1}=2.0$. The lime vector shows the pulsar proper motion, and the blue cross marks the estimated pulsar birthplace. Asymmetry in the emissions due to the proper motion of the PWN and SDZ is significant for 10 GeV IC emission (a), and 100 GHz synchrotron emission (c). While higher symmetry is found in IC emission at 10 TeV (b), and synchrotron emission at 5 keV (d).

The pulsar’s proper motion does not significantly modify the emission morphology in the 8–40 TeV range. The corresponding surface brightness and the spectra of Model B are shown in Fig. 7. In our study, we examined how the strength of the magnetic field affects non-thermal emission and positron flux on Earth. Specifically, we looked at the inverse-Compton emission under three different random magnetic-field strengths (3 $\mu$G, 5 $\mu$G, and 10 $\mu$G) in the solar vicinity , assuming an injection index of $\gamma_{1}=2.0$. As we increased the magnetic-field strength, the surface brightness profile became more concentrated. This is because the stronger magnetic field reduces the cooling time of electrons/positrons, causing them to lose energy more quickly through synchrotron emission near the PWN. Based on our analysis, a magnetic-field strength close to the Galactic average level is favored for the IC spectrum in the nearby region of the PWN. After exploring various parameters, we found that an SDZ size of (50,70) pc and a magnetic-field strength of 3$\mu$G provided the best fit for the observed surface brightness.

Figure 7: IC emission of scenario B that the slow diffusion bubble moves with the pulsar proper motion. The injection spectra power index $\gamma_{1}=2.0$. The initial time-scale of pulsar $\tau_{0}=0.1\tau_{c}$. We show the expected emission of SDZ sizes of (30,50) pc and (50,70) pc, and magnetic-field strength of 3$\mu$G, 5$\mu$G, and 10$\mu$G. Top: Surface brightness in 8-40 TeV in terms of distance from the current location of Monogem pulsar. Bottom: IC spectra of 10^∘ region around Monogem pulsar, compared with Fermi-LAT upper-limits (Di Mauro et al., 2019) and HAWC observation from the region (Abeysekara et al., 2017).

The synchrotron spectra from 10^∘ region around the PWN are shown in the middle panel of Figure 4. We fix the benchmark parameters and compare the synchrotron emission by setting SDZ size of (30, 50) pc and (50, 70) pc, and magnetic-field strength $B=3$, 5, and 10 $\mu$G. As expected, a stronger magnetic field yields a higher synchrotron flux. Note that we fixed the diffusion coefficient in the simulations. If we tune the diffusion coefficient with the magnetic-field strength to make sure the TeV emission morphology always fits, the synchrotron emission should be indistinguishable at keV energy.

3.1.3 Scenario C: Moving PWN and Fixed Slow-diffusion Bubble

To evaluate the case in which the slow-diffusion zone is related to the expansion of its parent SNR “the Monogem Ring,” we fixed the SDZ at the Monogem pulsar’s birth location and left it static in our simulation. Since we assume constant proper motion of the pulsar and zero transverse velocity, the center of the SDZ slightly differs from the estimated center of the Monogem Ring (that comes from observations). Despite this, our simulation should capture the system’s evolution well. The SDZ has a current radius of 70 pc, consistent with the radius of the Monogem Ring in X-rays. Assuming an age of 0.9$\tau_{c}$, and no line of sight velocity, we estimate that Monogem traveled $\sim 6$ pc from its birthplace, making it well located inside the SDZ, as suggested by observations (Knies et al., 2018). However, we cannot exclude the possibility that the pulsar has a significant line of sight velocity component. Consequently, while the projected position of the pulsar appears within the Monogem Ring, the actual position of the pulsar may lie outside the SNR.

Figure 8 displays the intensity maps at 10 GeV and 10 TeV from IC emission and 100 GHz and 5 keV from synchrotron emission. Compared to Fig. 6, where the SDZ moves with the pulsar, we can see that the emission morphologies and intensities are comparable in every energy windows. Given the small proper motion of the pulsar perpendicular to the line of sight, the projected birthplace and the current pulsar location are only 1.2^∘ apart in the sky. The impact of anisotropic diffusion due to the Monogem Ring over the pulsar’s lifetime is minimal. In this scenario, the 10 TeV emission is approximately spherically symmetric around the pulsar, as predicted also in scenario B. This suggests that the asymmetric TeV halo observed by HAWC is unlikely to be caused by the spatial asymmetry of the diffusion coefficient introduced by the Monogem Ring. More realistic modeling of the diffusion bubble, taking into account the interaction between the PWN and the SNR, may help explain the observed lobed TeV emission.
We show the comparison of synchrotron and IC spectra for scenarios A, B, and C in the bottom panel of Figure 4. The plot is shown for SDZ size of (50, 70) pc, magnetic field $B=3\,\mu$G, injection spectral index $\gamma_{1}=2.0$, and $\tau_{0}=0.1\tau_{c}$. Same as the emission morphology we discussed above, the emission spectrum is also indistinguishable for scenarios B and C.

(a) IC emission at 10 GeV

(b) IC emission at 10 TeV

(c) Synchrotron emission at 100 GHz

(d) Synchrotron emission at 5 keV

Figure 8: Intensity maps at different wavelengths for the region around Monogem pulsar taking into account the pulsar proper motion. The center of the SDZ stays at the birth location of the pulsar. The SDZ has a size of $(r_{1},r_{2})=(50,70)$ pc at current time. The initial time-scale of pulsar $\tau_{0}=0.1\tau_{c}$. Magnetic-field strength $B=3\mu G$, and injection index $\gamma_{1}=2.0$. The lime vector shows the pulsar proper motion, and the blue cross marks the estimated pulsar birthplace. Asymmetry in the emissions due to the proper motion of the PWN and SDZ is significant for 10 GeV IC emission (a), and 100 GHz synchrotron emission (c). While higher symmetry is found in IC emission at 10 TeV (b), and synchrotron emission at 5 keV (d).

3.2 Positron Flux

The PAMELA telescope observed a surplus in the CR positron-to-electron ratio at energies greater than 100 GeV, compared to the ratio predicted by the secondary production model (Adriani et al., 2009). This observation was later confirmed by Fermi-LAT (Ackermann et al., 2012) and AMS-02 (Aguilar et al., 2013), amongst others. Many authors have proposed various theories to explain this excess in the measured positron fraction, including nearby astrophysical leptonic sources (e.g., Yüksel et al., 2009; Hooper et al., 2009), and dark matter annihilation (e.g., Ibe et al., 2013; Dev et al., 2014). Geminga and Monogem are well-known nearby sources of $e^{\pm}$ that have emerged as potential sources of the positron excess. The TeV halos observed around these two pulsars provide strong evidence of leptonic particles escaping into the interstellar medium from these sources. However, the extension of the TeV halos also indicates that the $e^{\pm}$ are strongly confined, which suggests a reduced flux of positrons from these sources to Earth.

We have conducted a study on the primary positron spectra that are expected at Earth for different combinations of model parameters. In Fig. 9, we show the data collected by AMS-2 over seven years and compare it to our prediction. The yellow and green solid lines on the plot represent the two preferred combinations of model parameters by the IC surface brightness, which are $(r_{1},r_{2})=(50,70)$ pc, $\gamma_{1}=2.0$, and $B=3\mu$G and $5\mu$G, respectively.

Our study has found that the contribution of positrons with energy greater than 100 GeV from Monogem alone is only a few percent. We have also observed that the peak of the positron flux shifts to lower energies as the ambient magnetic-field strength increases because positrons lose their energy faster around the PWN. We have varied the magnetic field on the Galactic scale, and as a result, we expect the cutoff effect shown in the plot to be more dramatic than it would be in the actual case where only the magnetic field around the source is varied. We have also found that increasing the size of SDZ results in a longer confinement time of positrons. Additionally, high-energy positron spectra become softer due to synchrotron and IC losses.

Despite this result, the pulsar hypothesis remains the most compelling one. Recent studies (e.g., Orusa et al., 2021) have conducted sophisticated simulations of pulsar populations in the Galaxy, showing that nearby bright sources could account for the positron excess. These findings are consistent with previous analyses (Profumo et al., 2018; Jóhannesson et al., 2019; Schroer et al., 2023; Martin et al., 2022) of the Geminga pulsar, that show that the Geminga pulsar could contribute a significant proportion of the observed positron excess.

Figure 9: Positron flux from Monogem pulsar expected at Earth in comparison with AMS-02 7 years data. $\gamma_{1}=2.0$, and $\tau_{0}=0.1\tau_{c}$ for all cases.

4 Conclusions

In this article, we utilized the GALPROP (v57) framework to calculate the expected spectrum and spatial distribution of interstellar diffuse gamma-ray, including IC and synchrotron emission, and positron spectra at Earth. We have demonstrated that HAWC’s observations of the TeV halo surrounding PWN B0656+14 (Monogem) can be well explained by a two-zone diffusion model. Our modeling suggests that a diffusion bubble of about $\sim 50$–70 pc in size encompasses Monogem, where the diffusion coefficient experiences a suppression similar to that of the Geminga pulsar. These results imply that the slow diffusion of ultrarelativistic $e^{\pm}$ pairs around Monogem and Geminga may be due to the same mechanism.

The proper motion of the Monogem pulsar introduces asymmetries in the GeV IC emission and the radio synchrotron emission around the pulsar’s current location in both scenarios: whether the slow diffusion bubble is always centered at the pulsar or fixed at its birthplace. Our analysis shows that the non-thermal spectra and morphologies are indistinguishable between these two cases due to the small distance the pulsar has traveled perpendicular to the line of sight over its lifetime. Notably, both scenarios predict spherically symmetric 10 TeV emission, contrary to the asymmetrical emission morphology observed by HAWC. This suggests that the slow diffusion of pairs around Monogem may be due to additional factors beyond magnetic turbulence caused by the Monogem Ring or the propagation of the CRs themselves. Possible reasons could include pre-existing features of the diffusion coefficient in the surrounding ISM. Furthermore, realistic modeling that determines the exact location of the Monogem relative to the Monogem Ring, and their interactions, is necessary to understand the emission region. We note here that the TeV halo is expected to have synchrotron halo counterpart in the radio to X-ray energy range. However, with the current capabilities of radio and X-ray telescopes, detecting such an extended object in the sky remains a challenge.

Based on comparisons between our models and the data, we obtained that the escape efficiency of ultrarelativistic $e^{\pm}$ from the PWN into the ISM is approximately $\sim 5$–$10\%$ for an injection spectrum with $\gamma_{1}=2.0$. In this case, the corresponding GeV emission spectrum is around one order of magnitude lower than the upper limits obtained with Fermi-LAT (Di Mauro et al., 2019). If the injection spectrum is softer ($\gamma_{1}\gtrsim 2.0$), the resulting IC spectrum will, in turn, be softer, leading to a steeper profile in the GeV range. We note that an improved analysis of Fermi-LAT observations could be used to constrain the properties of the SDZ further. Although it might be challenging to distinguish between the size of the diffusion bubble and the shape of the $e^{\pm}$ injection spectra, our predicted asymmetry in the spatial morphology of the GeV gamma-ray emission could help disentangle the Monogem pulsar from the diffuse gamma-ray background in Fermi-LAT observations.

Moreover, we showed that the strength of the magnetic field has a significant impact on our simulations. After conducting several tests, we determined that the HAWC observed gamma-ray spectrum and spatial profile can be appropriately explained if the Galactic magnetic field around Monogem is close to the Galactic average of $\sim 3$–$5\;\mu$G. When we tested higher magnetic fields, we obtained that they shortened the cooling time of $e^{\pm}$, resulting in the need for a smaller diffusion coefficient to describe the observed gamma-ray emission morphology. Specifically, the predicted energy spectrum drops sharply in the TeV range for a magnetic-field strength of 10 $\mu$G in the local ISM. However, this may still match the data depending on whether a smaller diffusion coefficient inside the slow diffusion bubble is assumed.

Recent observations of Monogem (and Geminga) with the HAWC telescope at TeV energies suggest that diffusion of ultrarelativistic $e^{\pm}$ within PWNe is less effective than the rest of the interstellar medium. Initially, it was claimed (Abeysekara et al., 2017) that the $e^{\pm}$ produced by these pulsars may not contribute significantly to the locally measured positron flux if the diffusion coefficient in the local interstellar medium is similar to the value inferred within the nebula. Here, in agreement with various previous studies (Profumo et al., 2018; Jóhannesson et al., 2018, e.g.,), we have shown in detail that if the diffusion is slow around the PWN, then the expected contribution by Monogem is at the $\sim 10\%$ level, which is not insignificant, considering that there are many more objects like it within the inner $\sim 2$ kpc of the Solar system. Indeed, recent studies (e.g., Orusa et al., 2021) have demonstrated that a population of bright nearby pulsars could explain the AMS-02 positron data in detail.

Finally, we have shown that two-zone diffusion models can adequately explain multi-wavelength data from Monogem. A similar conclusion was obtained by Porter et al. (2017) in the case of Geminga. If slow-diffusion bubbles are a common feature of PWNe in our Galaxy, regions of the sky with high concentrations of these sources (e.g., the Galactic center) may have significantly higher cosmic-ray densities than those predicted by the standard diffusion model. These findings motivate further investigation of diffuse non-thermal radiation from the Galactic disk and bulge in which slow-diffusion bubbles wrap at least a fraction of the Galactic CR sources.