Pulsar wind nebulae with bow shocks: non-thermal radiation and cosmic ray leptons

TL;DR

This paper explores bow shock pulsar wind nebulae, focusing on their non-thermal radiation, particle reacceleration, and potential contribution to galactic cosmic ray positrons, supported by multiwavelength observations and theoretical analysis.

Contribution

It introduces the concept of particle reacceleration in BSPWNe, explaining hard spectra and their role in cosmic ray positron production, with detailed observational case studies.

Findings

Reacceleration produces hard particle spectra with photon indices < 1.5.

BSPWNe can significantly contribute to galactic positron flux.

High-resolution imaging reveals pulsar wind and ISM interactions.

Abstract

Pulsars with high spin-down power produce relativistic winds radiating a fraction of the power in the range from radio to gamma-rays in the pulsar wind nebulae (PWNe). The rest of the power is dissipated in the interactions of the PWNe with the interstellar medium (ISM). Some of the PWNe are moving relative to the ISM with supersonic speeds producing bow shocks. In this case, the ultrarelativistic particles accelerated at the termination surface of the pulsar wind may be reaccelerated in the converging flow system formed by the outflow from the wind termination shock and the inflow from the bow shock. An outcome of this reacceleration is the creation of particle distributions with hard spectra, such as required to explain the observed synchrotron spectra with photon indices Gamma <~ 1.5. The presence of this hard component is specific to PWNe with bow shocks (BSPWNe). The accelerated particles may end up containing a substantial fraction of the shock ram pressure. For typical ISM and pulsar parameters, the positrons released by these systems are numerous enough to contribute a substantial fraction of the positrons detected as galactic cosmic ray particles above few tens and up to several hundred GeV. The escape of ultrarelativistic particles from a BSPWN and its appearance in the far-UV and X-ray bands is determined by the directions of the interstellar magnetic field, the velocity of the astrosphere and the pulsar rotation axis. In this respect we review the observed appearance and multiwavelength spectra of three different types of BSPWNe: PSR J0437-4715, the Guitar and Lighthouse nebulae, and Vela-like objects. We argue that high resolution imaging of such objects provides unique information on pulsar winds and on the ISM. We discuss the interpretation of imaging observations and estimate the BSPWN contribution to the positron flux observed at the Earth. (ABRIDGED)

Figures

Equations

$$T_{\mathrm{esc}}={\left(\frac{R_{\mathrm{ST}}}{V_{\mathrm{PSR}}{T}_{\mathrm{ST}}^{2/5}}\right)}^{\frac{5}{3}}\approx 45\mathrm{kyr}{\left(\frac{E_{\mathrm{SN}}}{10^{51}\mathrm{erg}}\right)}^{\frac{1}{3}}{\left(\frac{n_H}{{\mathrm{cm}}^{-3}}\right)}^{-\frac{1}{3}}{\left(\frac{V_{\mathrm{PSR}}}{400km/s}\right)}^{-\frac{5}{3}},$$

$$I\Omega \dot{\Omega }=\frac{4}{9}\frac{B_{\ast }^2{R}_{\ast }^6{\Omega }^4}{c^3},$$

$$\Omega (t)=\frac{{\Omega }_0}{{\left(1+t/{\tau }_b\right)}^{\frac{1}{2}}},\mathrm{with}{\tau }_b=\frac{9I{c}^3}{8B_{\ast }^2{R}_{\ast }^6{\Omega }_0^2},$$

$$E_{\mathrm{res}}=\frac{1}{2}I{\Omega }^2(T_{\mathrm{esc}})=\frac{1}{2}\frac{I{\Omega }_0^2}{\left(1+t_{\mathrm{esc}}/{\tau }_b\right)}.$$

$${\mathcal{M}}_{\mathrm{a}}=v_{\mathrm{sh}}(4\pi {\rho }_{\mathrm{a}}{)}^{1/2}/B_0\approx 460v_8{n}_{\mathrm{a}}^{1/2}/B_{-6},$$

$${\mathcal{M}}_{\mathrm{s}}\approx 85v_8\cdot [T_4\cdot (1+f_{\mathrm{ei}}){]}^{-1/2},$$

$$R_{\mathrm{cd}}\approx \sqrt{\frac{{\zeta }_{\mathrm{K}}\dot{E}}{6\pi {\rho }_{\mathrm{ism}}{V}_{\mathrm{psr}}^{2}c}}\approx 2\times 10^{16}{\dot{E}}_{35}^{1/2}{n}_{\mathrm{ism}}^{-1/2}{V}_{200}^{-1}\mathrm{cm},$$

$${\mathcal{L}}_{\mathcal{M}}\approx 3\times 10^{33}\cdot \frac{{\mathrm{\Gamma }}_{\mathrm{flow}}^2}{{\beta }_{\mathrm{flow}}}{\left(\frac{E/Z}{10^{14}\mathrm{eV}}\right)}^2\mathrm{erg}{\mathrm{s}}^{-1},$$

$$E_{\max }\approx Z\times 10^{14}\cdot \frac{{\beta }_{\mathrm{flow}}^{1/2}}{{\mathrm{\Gamma }}_{\mathrm{flow}}}{\left(\frac{\dot{E}}{3\times 10^{33}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}\mathrm{eV},$$

$${\gamma }_{\max }\approx 1.9\times 10^7\cdot {\left(\frac{\dot{E}}{10^{38}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{3/4}{\left(\frac{\dot{N}}{10^{40}{\mathrm{s}}^{-1}}\right)}^{-1/2},$$

$$D(E_{\max })\sim v_{\mathrm{sh}}{R}_{\mathrm{sh}},$$

$$\frac{\delta B_{\mathrm{w}}^2}{8\pi }\approx \frac{v_{\mathrm{sh}}}{2c}{U}_{\mathrm{CR}},$$

$$\delta B_{\mathrm{w}}\approx \sqrt{2\pi {\zeta }_{\mathrm{CR}}\rho v_{\mathrm{sh}}^3/c}\sim 20{\zeta }_{\mathrm{CR}}^{1/2}{v}_{sh8}^{3/2}{n}^{1/2}\mu G,$$

$$\delta B_{\mathrm{w}}\approx \sqrt{4\pi {\zeta }_{\mathrm{B}}\rho }{v}_{\mathrm{sh}}\sim 50{\left[\frac{{\zeta }_{\mathrm{B}}}{0.01}\right]}^{1/2}{v}_{sh8}{n}^{1/2}\mu G.$$

$$E_{\max }$$

$$\approx \frac{0.2Z({\zeta }_{\mathrm{K}}{\zeta }_{\mathrm{CR}}{)}^{1/2}}{\eta }{\dot{E}}_{35}^{1/2}{v}_{sh8}^{3/2}\mathrm{TeV}$$

$$E_{\max }$$

$$\approx \frac{4.5Z({\zeta }_{\mathrm{K}}{\zeta }_{\mathrm{B}}{)}^{1/2}}{\eta }{\dot{E}}_{35}^{1/2}{v}_{sh8}\mathrm{TeV}.$$

$$\Lambda (p)=\frac{4}{\pi }\frac{R_{\mathrm{st}}^2(p)}{L_{\mathrm{cor}}}\propto p^2,$$

$$D_{\perp }\approx D_{\parallel }(p)\frac{1}{1+{\zeta }^2},D_{\parallel }(p)=\frac{v\Lambda (p)}{3},\zeta =\Lambda (p)/R_{g0}.$$

$$N_i^{\pm }(z,\gamma ,t)=\frac{N_0({\gamma }_0,{\gamma }_m)}{\gamma }\cdot H(\gamma -{\gamma }_0)\cdot H(t-{\tau }_a(\gamma ))\cdot \exp \left(-u_{\mathrm{i}}\cdot |z|/D_{\mathrm{i}}(\gamma )\right),$$

$${\tau }_a=\frac{3}{\left(u_1+u_2\right)}\underset{{\gamma }_0}{\overset{\gamma }{\int }}\left(\frac{D_1(\gamma )}{u_1}+\frac{D_2(\gamma )}{u_2}\right)\frac{d\gamma }{\gamma }.$$

$${\tau }_a\approx \frac{9}{c}\underset{{\gamma }_0}{\overset{\gamma }{\int }}\frac{D_2(\gamma )}{u_2}\frac{d\gamma }{\gamma }.$$

$$\frac{9D({\gamma }_{m1})}{c{u}_{\mathrm{sh}}}=\frac{R_{\mathrm{cd}}^2}{4D({\gamma }_{m1})}.$$

$$R_g({\gamma }_{m1})\approx \frac{1}{2}\sqrt{\frac{u_{\mathrm{sh}}}{c}}{R}_{\mathrm{cd}}.$$

$${\tau }_s(\gamma )\approx 7\times 10^{19}{\gamma }^{-1}{\left[1+{\left(\frac{B}{3\mu G}\right)}^2\right]}^{-1}s,$$

$${\gamma }_{m2}^2\approx 4\times 10^{17}{\eta }^{-1}\left(\frac{u_2}{100\mathrm{km}{\mathrm{s}}^{-1}}\right)\left(\frac{B}{3\mu G}\right){\left[1+{\left(\frac{B}{3\mu G}\right)}^2\right]}^{-1}.$$

$$r\left(\theta \right)=\{\begin{array}{rl}& a_{in}\text{mbox},if\theta \le \pi /2\\ & a_{in}/\sin \theta \text{mbox},if\theta >\pi /2.\end{array}$$

$$r\left(\theta ,\sigma \right)=\sigma R_{cd}\csc \theta \sqrt{3\left(1-\theta \cot \theta \right)}.$$

$${\lambda }_1\left(B,\eta ,\gamma \right)=\eta R_g\left(B,\gamma \right)=\eta \frac{m_e{c}^2\gamma }{eB}.$$

$${\lambda }_{\mathrm{am}}\left(\gamma \right)=3\times 10^{18}{\left(\frac{E}{1\text{mbox}GeV}\right)}^{1/3}\text{mbox}cm.$$

Full text

††journal: SSRv

1 Introduction

Pulsars are magnetized fast rotating neutron stars, whose rotation energy is carried away by relativistic winds blowing out of their magnetic poles and mainly consisting of electron-positron (e±) pairs. Multiband emission observed from these winds allows one to study the complex physics of energy transport in space (Arons, 2002; Gaensler and Slane, 2006; Kirk et al., 2009; Arons, 2012; Bucciantini, 2014; Kargaltsev et al., 2015). Multi-wavelength observations have revealed a vast variety of radiative efficiencies of pulsar winds (PWs) and pulsar wind nebulae (PWNe) produced by the winds, ranging from some tens percent in the Crab nebula (see e.g. Bühler and Blandford, 2014) to much lower values in a number of other PWNe (see, e.g., Li et al., 2008; Kargaltsev and Pavlov, 2008). Relativistic particles leaving the observable PWNe interact with the ambient medium on larger scales (upto tens of parsec) and thus provide a tool to reveal the structure of magnetic field and some other parameters of the interstellar medium (ISM), as discussed below in §6. This is a probable case for the fast enough, supersonically moving pulsars.

An analysis of radio pulsar velocity distribution based on the large scale survey at 0.4 GHz by Arzoumanian et al. (2002) suggested a two-component velocity distribution with characteristic velocities of 90 and 500$\rm~{}km~{}s^{-1}$ and inferred that $\sim$15% of the pulsars have natal velocities greater than 1000$\rm~{}km~{}s^{-1}$. The authors also concluded that under some simplified assumptions on the supernova remnant (SNR) evolution, about 10% of the pulsars would leave the host remnant before they are 20 kyr old.

The catalog of 233 pulsars (both young and recycled) compiled by Hobbs et al. (2005) contains derived pulsar speeds measured in one coordinate (1D) and the transverse velocities (2D). No evidence for a bimodal velocity distribution was found. The authors derived the mean 1D pulsar speed of 152$\pm$10$\rm~{}km~{}s^{-1}$ for the normal pulsars and 54$\pm$6$\rm~{}km~{}s^{-1}$ for the recycled ones. The corresponding mean 2D velocities in the catalog are 246$\pm$22$\rm~{}km~{}s^{-1}$ and 87$\pm$13$\rm~{}km~{}s^{-1}$, with the highest values $\sim 1,600\rm~{}km~{}s^{-1}$ inferred for the pulsars B2224+64 and B2011+38. The mean 3D natal pulsar velocity was derived to be 400$\pm$40$\rm~{}km~{}s^{-1}$.

For typical pulsar and ISM parameters, most pulsars are expected to leave their parent SNRs during the Sedov-Taylor phase of remnant expansion. The time $T_{\rm esc}$ at which this occurs can be estimated as $V_{\rm PSR}T_{\rm esc}=R_{\rm ST}(T_{\rm esc}/t_{\rm ST})^{2/5}$, where $V_{\rm PSR}$ is the velocity of pulsar proper motion, $t_{\rm ST}$ and $R_{\rm ST}$ are the time and remnant radius at the beginning of the Sedov-Taylor phase. With an average estimated velocity of about 400$\rm~{}km~{}s^{-1}$, an estimate of the time at which a pulsar escapes from its SNR can be found as

[TABLE]

A pulsar aged a few tens of kyr still has a non-negligible fraction of its original spin-down energy to spend. In general, one can write $\dot{\Omega}\propto\Omega^{n}$, with $n$ – the pulsar braking index, $\Omega$ and $\dot{\Omega}$ – the star rotation frequency and its time derivative, respectively. For dipole spin-down one expects $n=3$, since $\dot{E}=I\Omega\dot{\Omega}\propto B_{*}^{2}\Omega^{4}$, where $I\approx 10^{45}$ g cm2 is the star momentum of inertia and $B_{*}$ – its surface magnetic field. While it is true that for the few cases, where the braking index has been actually measured, its value has always been found to be less than 3, with an average of 2.5 (see, e.g., Faucher-Giguère and Kaspi (2006) and references therein), in the following we will assume the standard dipole spin-down law. Hence, one can compute the rotational energy the pulsar is left with at the moment $T_{\rm esc}$ when it leaves the remnant. From the equation:

[TABLE]

where $R_{*}\approx 10$ km is the neutron star radius, we find:

[TABLE]

and $\Omega_{0}=2\pi/P_{0}$ – the initial star rotation frequency ($P_{0}$ is the initial star spin period). Finally, the residual rotation energy of the star at $T_{\rm esc}$ is:

[TABLE]

The average value of $E_{\rm res}$ can be estimated by taking the average values of the pulsar surface magnetic field and initial period. From the work of Faucher-Giguère and Kaspi (2006) we have: $\langle B_{*}\rangle\approx 4.5\times 10^{12}$ G and $\langle P_{0}\rangle=300$ ms. Using the estimate of $T_{\rm esc}$ from Eq. (1) one obtains $E_{\rm res}\approx 10^{47}$ erg and $\dot{E}_{\rm res}\approx 5\times 10^{34}$ erg/s.

The fast average proper motion means that a non-negligible fraction of the pulsars propagating through the ISM may form a cometary-like bow shock structure. Optical and radio observations discovered a few bow shocks associated with fast moving pulsars (see, e.g., Kulkarni and Hester, 1988; Cordes et al., 1993; Frail et al., 1996; Ng et al., 2012; Brownsberger and Romani, 2014). High resolution X-ray observations of PWNe revealed a rich collection of their morphologies (Kargaltsev and Pavlov, 2008) including a number of cometary type objects, some of which harbour bow shocks.

Bow-shock nebulae are known to originate from the interaction of winds from moving stars with the ambient medium and appear in association with very different types of stars: of solar type (see, e.g., Baranov et al., 1971), young massive stars (see, e.g., Weaver et al., 1977), and fast rotating magnetized neutron stars (see, e.g., Bucciantini and Bandiera, 2001; Blondin et al., 2001). The dynamics of astrophysical bow shocks has been modelled both analytically, within the thin shell approximation, and by means of numerical hydrodynamics schemes (see, e.g., Bandiera, 1993; Wilkin, 1996).

The rich set of observed morphological and spectral appearances of PWNe traces the differences in the intrinsic properties of the parent pulsars (spin-down power and inclination), in their proper velocities, and in the evolutionary stages during the interaction of PWNe with the ISM (Gaensler and Slane, 2006).

The structure of PWNe can be heavily affected by the supersonic motion of the medium in which the nebula is expanding. This can be either associated with the fast proper motion of the pulsar through the normal ISM or with the reverse shock of the supernova propagating back towards the origin of the explosion. The effect of the external ram pressure of the fast moving ambient medium inside the supernova remnant on the PWN structure was studied in the frame of magnetohydrodynamical models by van der Swaluw (2003); van der Swaluw et al. (2004). The authors applied the model to three SNRs: N157B, G327.1-1.1, and W44, and concluded that the head of the PWN is not bounded by a bow shock in the case of N157B, G327.1-1.1, while in the case of W44 they argued for the supersonic scenario with a bow shock PWN (BSPWN) due to the fast motion of the pulsar. Another case of a bow shock type structure occurs when the PW located inside an SNR is interacting with a mildly supersonic flow of Mach number  >$$\sim  1 produced in the transition phase when the reverse shock reaches the center of the remnant (Chevalier and Reynolds, 2011). A detailed one-zone model of spectral evolution of PWNe inside SNRs was presented by Bucciantini et al. (2011). The authors found that the one-zone model provides a satisfactory description of a number of PWNe assuming that relativistic pairs are injected in these nebulae with a spectral distribution in the form of a broken power law, rather flat at low energy (energy spectral index close to 1) and steepening (spectral index larger than 2) above some break energy. They found that the intrinsic spectral break turns out to occur at similar energies for PWNe with rather different characteristics. A more sophisticated modeling of BSPWN dynamics was undertaken by Bucciantini et al. (2005), who investigated it in the framework of axisymmetric relativistic MHD, highlighting the dependence of the dynamics on the PW magnetization.

The Alfvénic Mach number of a bow shock of speed $v_{\rm sh}$ (at the bow apex it is about the pulsar speed) is

[TABLE]

where $n_{\rm a}$ is the ionized ambient gas number density measured in $\rm~{}cm^{-3}$, $B_{-6}$ is the local magnetic field in the shock upstream measured in $\mu$G and $v_{8}$ is the shock velocity in 10${}^{8}~{}\rm~{}cm~{}s^{-1}$. The sonic Mach number for a shock propagating in the interstellar plasma of the standard (Solar) abundance is

[TABLE]

where $T_{4}$ is the plasma ion temperature measured in 104 K and $f_{\rm ei}=T_{\rm e}/T_{\rm i}$ is the ratio between the electron and ion temperatures.

Consider a pulsar moving supersonically with a velocity $V_{\rm psr}$ relative to the local rest frame of the ambient medium. The standoff distance $R_{\rm cd}$ of the PW termination shock in this frame can be estimated as

[TABLE]

where $\dot{E}$ is the pulsar spin-down power, $\rho_{\rm ism}$ is the ambient ISM density. In fact, $R_{\rm cd}$ is rather the location of the contact discontinuity, than either of the PWN termination shock, or of the bow shock. The parameter $\zeta_{\rm K}$ depends on the magnetic inclination angle of the pulsar; $\zeta_{\rm K}=$1 if the PW is spherically symmetric (e.g., van der Swaluw et al., 2003). However, observations of the Crab nebula as well as particle-in-cell and MHD modeling of the relativistic pulsar wind indicate that the wind power is likely highly anisotropic. Such anisotropy affects the value of $\zeta_{\rm K}$ and, therefore, $R_{\rm cd}$, which scales $\propto\zeta^{1/2}$. This means that if the direction of the pulsar proper motion is close to the equatorial plane of a low-inclination rotator (i.e., the direction of the maximal wind power) the position of the contact discontinuity will be at an angular distance \geq 6\hbox{{}^{\prime\prime}}\zeta_{\rm K}^{1/2}d_{250}^{-1} for the most likely pulsar and ISM parameters. The images like the one presented in Fig. 8 may allow us to constrain the wind power anisotropy.

The model of a PWN evolving inside the remnant of the parent supernova developed by van der Swaluw et al. (2004) employed a hydrodynamical approach to simulate the evolution of such a system when the pulsar’s velocity is high. The authors considered four different stages of PWN evolution: the supersonic expansion stage, the stage of reverse shock interaction followed by the subsonic expansion stage and, finally, the bow shock stage. Within this model the bow shock stage occurs at roughly half the crossing time, when the pulsar is located at about 0.68 times the radius of the SNR forward shock. The authors suggested that the pulsar in the W44 SNR has a bow shock because of its supersonic proper motion. However, Chevalier and Reynolds (2011) argued that even in the case of the Vela pulsar, which is located closer than 0.68 of the SNR radius to the center of the remnant, and whose measured proper velocity is below 70 $\rm~{}km~{}s^{-1}$ (Dodson et al., 2003b), a bow shock still can be formed. It is likely that the SNR reverse shock has recently passed over the Vela pulsar and a mildly supersonic flow of ${\cal M}_{\rm s}\sim$ 1.3 is naturally produced. Chevalier and Reynolds (2011) noted that van der Swaluw et al. (2004) had reached their conclusion having assumed that the SNR flow can be described by the Sedov solution, while there is a transition solution preceding the Sedov stage, which is likely to describe the current evolutionary stage of Vela more accurately. Within such a refined description, Chevalier and Reynolds (2011) found that for Vela-like pulsars the flow velocity downstream of the reverse shock of speed $v_{\rm sh}$ is characterized by a velocity 0.75$\,v_{\rm sh}$ which appears to be about 1.3 times the speed of sound.

In what follows, we discuss in some detail how particles are accelerated in PWNe (§ 2), with a specific attention to acceleration in the PWNe driving bow shocks. Namely, efficient acceleration of particles in colliding shock flows (CSFs) between the pulsar wind termination surfaces and bow shocks is considered in § 3 and a Monte-Carlo model is discussed in § 4, which is further employed to simulate acceleration, transport, and synchrotron emission of the high energy e*±* pairs in BSPWNe. In § 5 we review far-UV and X-ray observations of emission structures in the vicinities of some fast moving pulsars with BSPWNe and apply the model of § 4 to simulate the synchrotron images and spectra of the nebula around PSR J0437-4715. § 6 is devoted to peculiar extended X-ray structures associated with the Geminga PWN, the “Guitar nebula”, and the “Lighthouse nebula”. The PWNe are moving supersonically and have hard photon indexes of X-ray emission in the vicinity of the bow shocks. The fascinating extended structures observed in X-rays can be associated with multi-TeV regime particles escaping from the bow shock PWN astrospheres to the ISM producing the magnetized ballistic beams of accelerated particles along the interstellar magnetic field lines. In §7 we discuss the X-ray observations of several BSPWNe possibly interacting with the remnants of their natal supernovae, and discuss the Vela X nebula in §8 as a generic case, which also includes a broadband modeling of its spectrum. Finally, in §9, we discuss the role of Bow Shock PWNe as plausible sources of the positron excess recently highlighted by PAMELA and AMS-02 measurements.

2 Particle Acceleration in Pulsar Wind Nebulae

PWNe are powered by highly magnetized relativistic winds which are terminated at some distance from the star. At the termination surface a considerable fraction of the pulsar spin-down power is transferred to accelerated particles producing non-thermal radiation (see Rees and Gunn, 1974; Kennel and Coroniti, 1984; Chevalier, 2000; Arons, 2012; Sironi et al., 2015).

The maximal energy of accelerated particles is limited according to a given magnetic luminosity of the source ${\cal L_{M}}$ (see, e.g., Lemoine and Waxman, 2009). Such energy can be estimated by the condition on the particle acceleration time scale in the comoving frame of the flow, which should be shorter than the dynamical time of the flow. Within the validity limits of relativistic MHD the acceleration time exceeds the particle gyroperiod (this condition does not apply to the region where magnetic reconnection occurs).

Let us consider acceleration of particles of charge $Z$ in a Poynting flux dominated source of magnetic luminosity ${\cal L_{M}}$. According to the analysis of Norman et al. (1995); Waxman (1995); Lemoine and Waxman (2009), the luminosity required to reach the energy $E$ (in the absence of energy losses) is

[TABLE]

where $\beta_{\rm flow}$ is the flow speed divided by the speed of light $c$, and ${\rm\Gamma}_{\rm flow}$ the flow Lorentz factor. Therefore the maximal energy of a particle which is achievable in the MHD outflow without the account for energy losses is limited by:

[TABLE]

This approximate relation indicates that the most favourable sites to achieve the highest energy are indeed in the trans- or sub-relativistic flows of \beta_{\rm flow}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;1 located between the ultrarelativistic termination shock and the non-relativistic bow shock which we will discuss in §3 below.

Sironi et al. (2013) provided a detailed analysis of the maximal energies of the particles accelerated at relativistic shocks. They found that the relativistic perpendicular shocks propagating in e*±* pair plasma (which is the likely case in PWNe) can efficiently accelerate the pairs only if the magnetization in the upstream plasma is  <$$\sim  10*-3* (and even lower,  <$$\sim  3$\times$10*-5*, for the electron-ion plasmas). The maximal Lorentz factor of the pairs accelerated at the PW termination shock was estimated by Sironi et al. (2013) from the constraint that the acceleration time of the highest energy pairs be shorter than the termination shock crossing time. These authors found that

[TABLE]

where $\dot{N}$ is the total particle flux entering the nebula, which is somewhat uncertain and is likely within the range 3$\times$ 10{}^{38}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;\dot{N}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;10${}^{40}~{}~{}{\rm s}^{-1}$. The estimated maximal energy is rather low even for the Crab nebula in spite of its high spin-down power. Therefore Sironi et al. (2013) concluded that Fermi acceleration at the termination shock of PWNe is not a likely mechanism for producing the synchrotron emitting pairs and some other mechanisms should be considered.

A number of pulsars with observed bow shocks, including PSR B2224+65, which powers the Guitar nebula,  the old recycled pulsar J0437-4715,   PSR J0633+1746 (Geminga) and others, have a spin-down power \dot{E}_{35}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\; 0.3. Several BSPWNe harbor regions of extended X-ray emission, where the magnetic field is estimated to be a few tens of $\mu$G. In order to produce X-ray synchrotron photons in such a field, pairs with a Lorentz factor \gamma\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;108 are required. But according to Eq. (8) such high energies are compatible with the amount of spin-down power mentioned above only if {\rm\Gamma^{2}_{\rm flow}}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;3 while \beta_{\rm flow}\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\; 0.1, i.e., the outflow is transrelativistic in the particle acceleration zone. The Vela pulsar as well as the pulsar IGR J11014-6103 (which powers the Lighthouse nebula) have much higher spin-down power and the condition for acceleration of pairs above the Lorentz factor \gamma\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;108 is relaxed to {\rm\beta_{\rm flow}}\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;10*-3*.

In the light of these estimates, we will discuss below the acceleration of particles in fast flows around PW termination shocks and bow shocks.

2.1 Particle Acceleration and Magnetic Field Amplification at Bow Shocks

Diffusive shock acceleration (Axford et al., 1977; Krymskii, 1977; Bell, 1978; Blandford and Ostriker, 1978) has been demonstrated to be a very attractive mechanism of particle acceleration in supersonic flows of a moderate magnetization. There is a growing body of observational evidence in favor of this mechanism acting in SNRs, AGN jets and some other objects (see, e.g., Bell, 2013; Helder et al., 2012; Marcowith et al., 2016). Particle acceleration in fast MHD shocks is strongly coupled with magnetic field amplification (MFA) which is due to the growth of seed fluctuations by the cosmic ray (CR) driven instabilities. A number of such instabilities resulting in MFA can be excited during diffusive shock acceleration (see, e.g., Bell, 2004; Amato and Blasi, 2009; Bykov et al., 2012; Schure et al., 2012; Bykov et al., 2014).

The maximal energy $E_{\rm max}$ of a particle of charge $Z$ accelerated by a strong quasi-steady bow shock of velocity $v_{\rm sh}$ and characteristic radius $R_{\rm sh}$ is limited by the shock curvature. It can be roughly estimated from

[TABLE]

where $D(E)$ is the particle diffusion coefficient in the shock upstream. Particle propagation at a pulsar bow shock depends on the shock obliquity, i.e., on the angle between the proper velocity of the pulsar and the direction of the local mean interstellar magnetic field. If the shock is quasi-parallel, it is determined by the magnetic fluctuations amplified by the CR-driven instabilities, as shown below in § 2.2. On the other hand, a transverse shock geometry may reduce the radial diffusion coefficients and therefore drastically change the conditions of particle confinement.

2.2 Quasi-parallel Bow Shocks

The amplitude $\delta B_{\rm w}$ of the fluctuating magnetic field amplified by CR-driven instabilities in the upstream of a shock can be estimated from the saturation condition of Bell (2004):

[TABLE]

where $U_{\rm CR}\approx\zeta_{\rm CR}\,\rho v_{\rm sh}^{2}/2$ is the energy density of the CRs accelerated at the shock. Then the amplitude of the fluctuating magnetic field is

[TABLE]

where $\zeta_{\rm CR}$ is the fraction of shock ram pressure converted into the energy of the accelerated particles, $v_{\rm sh8}$ is the shock velocity in 108 cm s*-1* and $n$ is the number density of the ISM. This amplitude may reach \delta B_{\rm w}\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;B_{0}, where $B_{0}$ is the mean magnetic field in the ISM. In any case, in the vicinity of a bow shock, $\delta B_{\rm w}$ is likely well above the estimated amplitude of the field fluctuations in the background ISM turbulence at the scales below 1015 cm, which are important for the scattering of TeV range particles.

The mean free path $\Lambda(E)$ of a relativistic proton in the quiet ISM away from SNRs or other active space objects can be estimated from the galactic CR propagation models as $\Lambda(E)\sim\left(3\mbox{--}5\right)\times 10^{18}\,(E/\mbox{GeV})^{a}$ cm, where the index $a=$0.3–0.6 (see, e.g., Strong et al., 2007). Therefore, the background magnetic turbulence in the quiet ISM can not confine the relativistic $e^{\pm}$ pairs accelerated in the PWNe within a region of scale size $\sim R_{\rm cd}$ as given by Eq. (7), for a bow shock speed above 100 $\rm~{}km~{}s^{-1}$. Depending on the bow shock velocity, strong MFA by CR-driven instabilities at the bow shock would affect the propagation and confinement of relativistic $e^{\pm}$ pairs up to Lorentz factors  >$$\sim 107.

Nonlinear Monte-Carlo modeling of the diffusive shock acceleration process in the presence of MFA (Bykov et al., 2014) has revealed that the regime consistent with (12) holds for shock velocities below 5,000 $\rm~{}km~{}s^{-1}$. At higher shock velocities the energy density of the fluctuating magnetic field scales as a fraction $\zeta_{\rm B}$ of the shock ram pressure. The MFA efficiency is illustrated in Fig. 1 where the ratio $\zeta_{B}$ of the amplified magnetic field energy density to the shock ram pressure is shown as a function of the shock velocity. In the high shock speed regime one can estimate the amplitude of the enhanced magnetic field as

[TABLE]

The validity of the Bohm diffusion approach requires strong scattering of particles by quasi-resonant fluctuations, i.e., $\delta B_{\rm w}\sim B_{0}$. In such a case the diffusion coefficient $D(E)\approx\eta\,cE/(Ze\ \delta B_{\rm w}),$ where \eta\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;1 is a numerical factor (see, e.g., Reville et al., 2008).

Models of extended fast shocks demonstrated that diffusive shock acceleration is rather efficient in this framework with $\zeta_{\rm CR}>$ 0.2 (see, e.g., Helder et al., 2012). The energy density of the CR-driven magnetic fluctuations in the shock vicinity is dominated by the fluctuations with typical wavenumbers $k_{\ast}$ which satisfy the condition k_{\ast}R_{\rm g}(E_{\rm max})\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;1, where $R_{\rm g}$ is the gyroradius of a particle.

Using Eqs. (11), (7) one can obtain the maximal energy of accelerated ions as

[TABLE]

in the regime which is described by Eq.(13), while in the regime of high shock velocity described by Eq.(14) it can be approximated by

[TABLE]

Under the approximations described above, the maximal energies of the accelerated nuclei do not depend on the ambient ISM density but depend on the bow shock velocity in the ISM rest frame.

2.3 Propagation of Pulsar Wind-Accelerated e*±* Pairs at a Bow Shock

The protons accelerated at the bow shock produce magnetic fluctuations at wavelengths comparable to their gyroradii in the vicinity of the bow shock. Moreover, magnetic fluctuations of wavelengths up to a scale-length $L_{\rm cor}$ which is much longer than the gyroradii of protons of energy $E_{\rm max}$ can be produced via the inverse turbulent cascading and the non-resonant CR-driven instabilities (see, e.g., Bykov et al., 2013b, and references therein). These CR-driven fluctuations will determine the diffusive transport of $e^{\pm}$ pairs accelerated at the PW termination surface at energies above 10 GeV in the region between the contact discontinuity and the bow shock. The strong long-wavelength fluctuations would provide a Bohm-type diffusion coefficient $D(E)\approx\eta\,D_{B}(E)$ with $\eta\sim 1$ for e*±* pairs up to a Lorentz factor $\gamma_{1}$ determined by the condition $R_{\rm g0}(mc^{2}\gamma_{1})\approx L_{\rm cor}$, where $R_{\rm g0}(E)=E/eB_{0}$ is the particle gyroradius in the mean magnetic field $B_{0}$.

Accurate estimations of the scale-length $L_{\rm cor}$ are not available yet because of the limited dynamical range of numerical simulations of the CR-driven turbulence. Therefore, we consider $L_{\rm cor}$ as a free parameter which determines the Lorentz factor $\gamma_{1}$ where a transition from the Bohm-like diffusion to the small-scale scattering regime given by Eq. (17) takes place.

For particles with gyroradii larger than the scale-length of the strongly amplified magnetic turbulence $L_{\rm cor}$, the mean free path produced by the short-scale modes can be approximated as

[TABLE]

where $R_{\rm st}(p)=cp/(e\ \delta B_{\rm w})$. In section 4 we will apply the diffusion model to simulate distinctive features of synchrotron images of some BSPWNe. Apart from the effect on propagation of the $e^{\pm}$ pairs, the MHD flow between a fast bow shock and the PW termination shock (a colliding shock flow) may result in a specific way of acceleration of the highest energy $e^{\pm}$ pairs. The hard energy spectra of the pairs produced in the colliding shock flow (CSF) can be responsible for the hard spectra of synchrotron and inverse Compton photons (of photon indexes ${\rm\Gamma}\sim$ 1.0–1.5) recently observed from the heads of some of BSPWNe in the X-rays and gamma-rays, respectively.

2.4 Quasi-Transverse Bow-Shocks

When a shock is highly oblique in respect to the direction of the local mean magnetic field, the transverse diffusion coefficient plays an important role in the diffusive shock acceleration. When the particle mean free path along the magnetic field $\Lambda\left(p\right)$ is large compared to the gyroradius $R_{\rm g0}$, the transverse diffusion coefficient $D_{\perp}$ is much smaller than that along the mean field $D_{\parallel}$:

[TABLE]

When $\zeta\gg 1$ one obtains $D_{\perp}=[R_{\rm g0}(p)/\Lambda(p)]^{2}D_{\parallel}(p)$. Thus, the orientation of the local magnetic field can significantly affect the transport of particles across the shock front, and, consequently, the efficiency of particle acceleration (see, e.g., Jokipii, 1987; Takamoto and Kirk, 2015).

The different bow shock obliquities may explain the substantial differences in the luminosities and observed spectra of BSPWNe with otherwise similar properties, in particular, the nebulae associated with PSR J0437-4715 and PSR J1741-2054.

3 Particle Acceleration in Colliding Shock Flows

Diffusive shock acceleration with MFA is rather an efficient way to produce relativistic particles. However, an even more efficient realization of the Fermi acceleration mechanism is acceleration of particles in colliding MHD flows (colliding shock flows, CSFs) carrying magnetic inhomogeneities, which may scatter relativistic particles (Bykov, 2005; Bykov et al., 2013). This kind of MHD flows is likely realized in a number of astrophysical objects, in particular, when an expanding shell of a supernova is closely approaching or interacting with a fast powerful wind of a nearby young massive star. Converging MHD flows occur in the vicinity of a contact discontinuity of an astropause, downstream of a fast bow shock produced by a supersonically moving astrosphere driven by a stellar wind (the relativistic PW in the case of a PWN).

Consider a simple analytical test particle kinetic model of particle acceleration by the Fermi mechanism in a CSF (Bykov, 2005). Two colliding plane-parallel flows of velocities $u_{\rm 1}$ and $u_{\rm 2}$ aligned with the positive and negative directions of axis $z$ and occupying half-spaces $z<0$ and $z>0$, respectively, are separated by a discontinuity in the plane $z=0$. The flows carry fluctuating magnetic fields amplified by anisotropic energetic particle distributions (e.g., the CR current driven instability). The particles propagate in a diffusive medium described by the fast particle diffusion coefficients $D_{\rm i}(\gamma)$ which depend on the particle Lorentz factor $\gamma$.

This simplified model allows us to illustrate the Fermi acceleration in a CSF of a BSPWN and obtain some important estimates. The collision of the flow emerging from the PW termination region with that coming from the bow shock can be viewed within the discussed simple scheme, but in the vicinity of the contact discontinuity the flows have a complex multi-dimensional structure. Note that the high CR pressure in this region may be dynamically important for the flow structure if the MHD turbulence provides a strong coupling of the CRs with the flow. However, the highest energy e*±* pairs accelerated in the PW termination region may have so high Lorentz factors that their mean free paths become larger than the characteristic size of the multi-dimensional bypassing flows area in the vicinity of the contact discontinuity. Let us denote the minimal Lorentz factor corresponding to this regime as $\gamma_{0}$.

Then the accelerated pair distribution function above $\gamma_{0}$ in between the two shocks can be approximated by the solution obtained from the kinetic equations in the discussed simple plane-parallel model:

[TABLE]

where $H(x)$ is the standard Heaviside step function, $\gamma_{m}$ is the maximal Lorentz factor of an accelerated particle, and $N_{0}(\gamma_{0},\gamma_{m})$ is a normalization constant. The acceleration time of the highest energy e*±* pair is

[TABLE]

We take $u_{1}\approx c/3$, $u_{2}\ll c$ in accordance with the typical flow velocities downstream of the PW termination shock and at the bow shock, respectively, while $D_{2}(\gamma)=\eta cR_{g}(\gamma)/3$ ($\eta\geq 1$) is the Bohm diffusion coefficient in the half-space $z>0$. Here $R_{g}(\gamma)$ is the particle gyroradius in the magnetic field amplified by the CR-driven instabilities. In the limit $u_{2}\ll c$ Eq. (20) simplifies to

[TABLE]

As the mean magnetic field behind the termination shock of a PWN is azimuthal, the radial transport of the pairs at high energies across the mean field is due to the transverse diffusion. The pair propagation at the bow shock depends on the shock obliquity. The transverse shock geometry reduces the radial diffusion coefficients.

3.1 The Maximal Energies of e*±* Pairs Accelerated in CSFs

The maximal Lorentz factors $\gamma_{m}$ of the e*±* pairs accelerated in CSFs are mainly limited either (i) by particle escape from the accelerator or (ii) by the particle energy losses.

In the first case, the particle confinement condition in a CSF requires that the acceleration time Eq.(21) does not exceed the time of diffusion through the region within the bow shock. Thus, the maximal Lorentz factor should not exceed $\gamma_{m1}$ satisfying the equation

[TABLE]

Therefore, assuming the Bohm diffusion in the accelerator, the upper limit on $\gamma_{m1}$ can be obtained from

[TABLE]

However, the estimations based on the analytic one-dimensional solutions given by Eqs. (19)–(20) are rather rough. Hence in §4 we will discuss test particle Monte-Carlo simulations of the CR acceleration and transport in BSPWNe.

In the second case, the synchrotron-Compton energy loss time of the pairs, $\tau_{s}$, should be longer than the acceleration time. This provides another constraint on the maximal Lorentz factor $\gamma_{m}\leq\gamma_{m2}$ given by the equation $\tau_{s}(\gamma_{m2})=\tau_{a}(\gamma_{m2})$. Using a simplified approximation for $\tau_{s}$ which accounts for inverse Compton losses on the cosmic microwave background photons (as appropriate for very high energy pairs) and synchrotron losses in the magnetic field $B$:

[TABLE]

one obtains

[TABLE]

Then the maximal Lorentz factor $\gamma_{\rm m}$ of the pairs accelerated in a CSF can be estimated as $\gamma_{\rm m}={\rm min}(\gamma_{m1},\gamma_{m2})$.

The spectral upturn at the high energies, namely, the fact that a large fraction of energy is carried by the particles with the longest mean free path, can result in an important non-linear modification of the flow in the region of convergence. Such an effect was studied in the context of CR-modified shocks by Bykov et al. (2013) within a plane-parallel time-dependent model. In the case of interest here, i.e., for the diffusive transport and re-acceleration of the high energy end of the spectrum of the $e^{\pm}$ pairs which were produced at the PW termination shock, the non-planarity is very important. Therefore, in § 4 we discuss a Monte-Carlo model of acceleration and propagation of high energy $e^{\pm}$ pairs in the astrosphere produced by the relativistic PW which moves supersonically through the ambient medium.

Acceleration in a CSF can provide very fast and efficient creation of a nonthermal particle population in between the PW termination region and the bow shock. A very hard energy spectrum of the e*±* pairs is formed, $\propto\gamma^{-b}$ with $b\sim$1, which may contain a substantial fraction of the kinetic power released by the pulsar in the highest energy particles. The synchrotron X-ray emission from the extended region between the PW termination shock and the bow shock where acceleration of pairs occurs would have a hard photon index, 1$\leq{\rm\Gamma}\leq$1.5. The associated inverse Compton emission may reach the TeV regime and extend to large distances downstream of the bow shock. Such an extended structure of hard emission may be treated as a distinctive feature of BSPWNe.

4 The Monte-Carlo Particle Acceleration and Transport Model

The modeling of the structure of PWNe interacting with supersonic flows either inside or outside of supernova remnants with an adequate account for the magnetic field structure is a complicated problem (see, e.g., Blondin et al., 2001; Bucciantini and Bandiera, 2001; van der Swaluw et al., 2004; Bucciantini, 2014). The study of particle transport in these systems faces the complexities related to a non-trivial behaviour of magnetized flows in the vicinity of the contact discontinuity. Fortunately, for electrons and positrons emitting synchrotron radiation in the UV and X-ray range this problem is alleviated due to their large mean free paths. Along with the deceleration of the colliding flows and dissipation of the transported magnetic inhomogeneities it significantly simplifies the problem. The transport of high energy e*±* pairs is not very sensitive to the complicated structure of the multidimensional flow in the vicinity of the contact discontinuity. However, the flow advects the low energy pairs along the nebula trail. We discuss here a simplified Monte-Carlo model of the high energy pair transport in BSPWNe.

4.1 The Model Approach

Simulation of the relativistic PW particle acceleration and transport through a PWN interacting with a supersonic flow is performed via a stationary Monte-Carlo approach. The e*±* pairs injected into the system at the PW termination surface are treated as test particles and propagated through the PWN with a fixed set of diffusion parameters, flow velocities, and geometry. An important feature of the approach is the account of the main geometrical properties of the modelled system.

4.2 The Model Geometry

The modelled system is assumed to be axisymmetic and treated as a composition of axisymmetric regions with spatially homogeneous but energy-dependent particle diffusion coefficients. A sketch of the modelled system geometry is presented in Fig. 2. The spherical coordinates $r$, $\theta$, and $\phi$ are centered at the pulsar position and the zenith direction is aligned with the pulsar velocity vector. Also the cylindrical coordinates ($\rho$, $\alpha$, $x$) orientated along the symmetry axis of the system ($Ox$) are used.

The innermost spherical region with radius $r_{ts}$ corresponds to the cold pulsar wind. Through its boundary – the termination shock – the particles are injected into the system. The region labeled “1” corresponds to the shocked pulsar wind (the PWN). Its outer boundary has the following shape:

[TABLE]

The scaling factor $a_{in}$ is a model parameter. The shocked PW flow velocity is directed radially and its value is $u_{w}=\frac{c}{3}r_{ts}^{2}/r^{2}$ at distance $r$ from the pulsar. The region labeled “3” corresponds to the postshock flow of the ambient ISM. Its boundaries are described by the equation for the bow shock shape derived by Wilkin (1996) in the thin shell approximation:

[TABLE]

Here $R_{\rm cd}$ is the standoff distance defined by Eq. (7). The positions of inner and outer boundaries of region “3” are determined by the model parameter $\sigma$. The apexes of these surfaces are denoted as $a_{\rm out}$ and $a_{\rm sys}$, respectively. The velocity of the flow through the bow shock in the pulsar rest frame is fixed and directed oppositely to the pulsar velocity in both regions “3” and “4”. In region “2”, located near the contact discontinuity, advection is neglected. Region “4” corresponds to the unperturbed ambient ISM. The system is bounded by a cylinder whose axis coincides with the system symmetry axis and whose height is $h_{\rm sys}=x_{\rm FEB}+a_{\rm sys}$, where $x=-x_{\rm FEB}$ is the coordinate of the tail-oriented base of the cylinder. The cylinder radius $\rho_{\rm sys}$ matches the radius of section of region “3” outer boundary by the $x=-x_{\rm FEB}$ plane. At the cylinder, free escape boundary conditions are imposed.

4.3 Particle Propagation

In each of the regions “1”-“4” a spatially uniform mean value of magnetic induction $B$ is specified. The parameters of particle diffusion are chosen as follows. In the PWN we adopt the mean free path $\lambda_{1}$ proportional to the particle gyroradius $R_{g}$:

[TABLE]

For the ambient medium the mean free path is taken in the form:

[TABLE]

In regions ”2” and ”3” the mean free path is chosen in the form:

[TABLE]

In Eqs. (27)–(29) $e$ and $m_{e}$ are the charge and mass of an electron or a positron, and $E=m_{e}c^{2}\gamma$ is its energy. The values of $\eta\geq 1$ and $\gamma_{1}$ are free parameters of the scattering model, while $\gamma_{2}$ is derived from the continuity of the mean free path value.

At the beginning of a simulation the initial spectrum of particles, with a given shape, is generated. Then, in the course of Monte-Carlo modeling of particle acceleration and transport, the generated particles are injected into the system one by one. In the developed numerical code the injection is implemented as follows. Two random numbers are generated for every given particle to define its initial position on the termination shock. Two more are generated to specify the initial velocity direction. The particle then propagates for a distance equal to its mean free path in region “1” and after that is scattered for the first time.

The scatterings of the particles are isotropic in the local plasma rest frame. At each scattering the value of particle momentum in the local plasma rest frame is calculated. Then two random numbers which define two scattering angles are generated and the particle momentum is rotated according to the paper of Ostrowski (1991). After that the new particle momentum in the pulsar rest frame is calculated. The value of $\delta\Omega_{\rm max}$, introduced by Ostrowski (1991) is taken equal to $\pi$ (due to that, $\eta$ actually should not be less than $2$).

After injection, the particle is propagated through the system. The particle moves along a straight line during the time interval $dt_{\rm mfp}=\Lambda/v$ (where $\Lambda$ is its mean free path at the current location, and $v\approx c$ is its velocity), and then is scattered. This cycle repeats until the particle finally crosses the free escape boundary, and, therefore, leaves the system. Then the next particle is injected and its propagation starts. When the particle crosses the border between two regions with different scattering parameters, the crossing point is determined and the particle is moved from this point in the same direction for a distance equal to $\lambda_{\rm new}(1-\Delta x/\lambda_{\rm old})$, where $\lambda_{\rm old}$ is the mean free path in the region it has left, $\lambda_{\rm new}$ is the mean free path in the region it has entered and $\Delta x$ is the path of the particle in the region it has left after the last scattering. If the particle enters the cold PW region, it is reflected from its boundary in the case when its gyroradius in the magnetic field downstream of the PW is smaller than $r_{ts}$ and crosses it without scattering otherwise.

The phase space of the system is divided into bins – small volumes of phase space, which are used to ’detect’ the particle spatial and momentum distribution. The coordinate space of the system is permeated by a cylindrical ’grid’ of bins. The cylinder containing the system is splitted into $N_{x}$ equal cylinders along the $Ox$ axis. Each of these smaller cylinders is divided into $N_{\rho}$ coaxial cylindrical layers, whose outer radii grow linearly as $\rho_{n}=n\rho_{\rm sys}/N_{\rho}$, where $n=1,2,\dots,N_{\rho}$. Each cylindrical layer is divided into $N_{\alpha}$ equal segments by semi-planes $\alpha=2\pi k/N_{\alpha}$, $k=1,\dots,N_{\alpha}$. To simulate the particle energy distribution, recording of the absolute values of particle momenta at each detector is required. In order to do this, we introduce a momentum range which is somewhat wider than the range of particle momenta in the initial (injected) particle distribution. This ’detector’ range is logarithmically binned; hence $N_{m}$ momentum bins are introduced. In case of considerable anisotropy of particle energy distribution function detection of particle angular distribution is also of some interest. Thus, the ranges of $\mu=\cos\theta$ and $\phi$ values are splitted into $N_{\mu}$ and $N_{\phi}$ bins, respectively. The particle is detected by the bin with a given set of numbers ($n_{x}$, $n_{\rho}$, $n_{\alpha}$, $n_{m}$, $n_{\mu}$, $n_{\phi}$) if its energy corresponds to the $n_{m}$-th momentum bin, direction of velocity to $n_{\mu}$-th and $n_{\phi}$-th angular bins, and it crosses the left boundary of the spatial bin ($n_{x}$, $n_{\rho}$, $n_{\alpha}$) corresponding to $x=-x_{\rm FEB}+(n_{x}-1)h_{\rm sys}/N_{x}$. To obtain the particle distribution function from the ’detection’ data, we employ the formulae analogous to those presented in §3.1.4 of (Vladimirov, 2009). In order to compensate for the growth of the bin size with the increasing distance from the axis, the values measured by the detectors are divided by the ratio of the current bin volume to the volume of the bin at $n=$1.

To take into account the effect of energy losses due to synchrotron radiation, after each straight line piece of trajectory, the particle momentum value is re-evaluated as $p_{2}=p_{1}/(1+p_{1}Sdt)$. Here $p_{1}$ is the momentum in the local plasma rest frame at the beginning of the current piece of trajectory, $p_{2}$ is the new momentum value in the same frame, $dt$ is the propagation time, and $S$ is the synchrotron loss coefficient:

[TABLE]

Since for all the explored sets of model parameters the particles do not suffer strong energy losses over a diffusion length, the particle energy is effectively decreased due to these losses only after the particle has crossed all the ’detectors’ along a mean free path. If a particle crosses the boundary between regions with different values of $B$, the decrease of its energy is calculated consequently for pieces of trajectory in each region.

4.4 The Model Spectra

The Monte-Carlo procedure described in §4.3 allows us to simulate particle spectra at certain locations inside the BSPWN system as well as spectra of synchrotron emission radiated by the accelerated particles. The intensity of synchrotron emission at various photon energies can be calculated for the lines of sight crossing these locations, thus allowing us to simulate emission maps in the sky plane. The locations chosen to present the particle spectrum modeling results are specified as follows. We choose a certain axial section plane of the axisymmetric BSPWN system. A right-handed Cartesian coordinate system $(xyz)$ centered at the point $P$, the pulsar position, is introduced. The $x$ axis matches the symmetry axis and the plane $Pxz$ matches the fixed axial section. We choose four groups of points in the $Pxz$ plane. Namely, we take four rays starting at $P$, lying in $z\geq 0$ and oriented at angles 0, 30, 60 and 90 degrees with respect to the symmetry axis. Along each ray we select three locations, each in the center of the interval delimited by the intersections of the ray with the surfaces separating regions with different diffusion properties. Thus, we select 3$\times$4 = 12 points in regions “1”–“3”. We include five additional points in region “4”. Five rays parallel to the $Pz$ axis and lying in the $Pxz$ plane ($z\geq 0$) are chosen at $x$ coordinates $x_{j}=-x_{\rm FEB}+j\left(a_{sys}+x_{\rm FEB}\right)/6,~{}~{}j=1..5$. The locations are taken at the centers of the intervals between the points of intersection of these rays with the outer boundary of region “3” and with the outer boundary of the whole system.

The orientation of the considered system with respect to the direction to the observer is defined as follows. We introduce an axis $Y$, directed from the observer towards the pulsar (point $P_{0}$) and the plane of the sky (PoS), orthogonal to this axis and containing $P_{0}$. The orientation of the source is determined by three angles $\Phi$, $\Theta$, and $\Psi$. $\Phi$ is the angle between the $Py$ axis and the line $l_{Y}$ of intersection of the planes $Pzy$ and PoS (chosen in the positive direction: counterclockwise when looking at the $Pzy$ plane against the $Px$ axis direction). $\Theta$ is the angle between $Y$ and $Px$ axes. $\Psi$ is the angle between $l_{Y}$ and the North direction in the PoS. Thus, one may introduce a coordinate system $(XYZ)$ centered at $P_{0}$, $X$ and $Z$ axes in the PoS, where the $Z$ axis is aligned with the direction to the North. Then the angles $\Phi$, $\Theta$, and $\Psi$ would be the three Euler angles defining the orientation of the $XYZ$ system with respect to the $xyz$ system. In the following and in Fig. 2 we use $\Theta=90^{\circ}$.

5 Observations and Modeling

of X-ray and Far-UV Emission from BSPWNe

Multi-wavelength observations in the radio, optical and X-ray bands have revealed the presence of bow-shock type structures in a number of PWNe (see, e.g., Pellizzoni et al., 2004; Gaensler and Slane, 2006).

An example of an apparent bow shock PWN (BSPWN) is the Mouse nebula powered by the radio pulsar J1747-2958, which is moving with transverse velocity (306 $\pm$ 43)$\times({d/5~{}{\rm kpc}})\rm~{}km~{}s^{-1}$ (Gaensler et al., 2004; Hales et al., 2009). The spin-down power of J1747-2958 is about 2.5$\times 10^{36}\rm~{}erg~{}s^{-1}$, while the X-ray luminosity of the nebula in the 0.5-8 keV range is about 5$\times 10^{34}\rm~{}erg~{}s^{-1}$, and the observed PWN structure shows a narrow 45*′′*-long X-ray tail (the length is $\sim$ 1.1 pc at the 5 kpc distance). The best-fit power law photon index of the X-ray spectrum ${\rm\Gamma}$ = $2.0\pm 0.3$ was derived from Chandra observations by Gaensler et al. (2004), who also noticed some steepening of the spectrum from the halo to the tail of the PWN.

High resolution observations of PWNe in Hα (see, e.g., Kulkarni and Hester, 1988; Cordes et al., 1993) have provided images of cometary-like bow shock structures ahead of fast moving pulsars. A recent survey of such bow shocks (Brownsberger and Romani, 2014) listed 9 pulsars with resolved Hα emission at the apex and a number of pulsars with upper limits set on their Hα flux.

5.1 Synchrotron Spectra and Images of the Far-UV Bow Shock

of PSR J0437-4715

One of the most interesting and best studied objects in the compilation of Brownsberger and Romani (2014) is the bow shock of the nearest millisecond pulsar J0437-4715. The pulsar is located at a reliably estimated distance of about 160 pc in a binary system with measured parallax and transverse velocity $\sim$ 104$\rm~{}km~{}s^{-1}$ derived from the motion of the pulsar companion and the bow shock apex measured over 17 yrs. The estimated spin-down power of J0437-4715 is $\dot{E}~{}\sim$ 6$\times$10${}^{33}\rm~{}erg~{}s^{-1}$.

Extended X-ray emission from the \sim 5\hbox{{}^{\prime\prime}} vicinity of J0437-4715 was detected with Chandra by Rangelov et al. (2016). This emission indicates the presence of a faint ($L\sim$ 3$\times$10${}^{28}\rm~{}erg~{}s^{-1}$) PWN produced by J0437-4715. The X-ray photon index of the PWN is ${\rm\Gamma}=1.8\pm 0.4$.

Far-ultraviolet (FUV) imaging of the vicinity of J0437-4715 performed by Rangelov et al. (2016) with the Hubble Space Telescope has revealed a bow-like structure which stands 10*′′* ahead of the pulsar and is coincident with the apex part of the Hα bow shock earlier observed by Brownsberger and Romani (2014). The observed 1250–2000 Å luminosity of the FUV bow is about 5$\times 10^{28}\rm~{}erg~{}s^{-1}$ – an order of magnitude higher than the Hα luminosity of the same structure. The FUV observations also revealed a likely flux enhancement in the extended (about 3*′′*) region at the limb of the bow, which is not seen in the wide-band optical/near-UV images.

Both Hα and FUV emission at the bow of J0437-4715 can be associated with line and continuum radiation of shocked interstellar plasma. In particular, even a simple one-dimensional model of gas flowing through a fast enough interstellar shock (see, e.g., Bykov et al., 2013a) can be used to explain the fluxes of broadband FUV and Hα emission observed by Rangelov et al. (2016). Within such a model the 1250–2000 Å band emission is generated in the hot downstream of the shock and dominated by several individual lines such as C IV 1549 Å, O IV 1403 Å, Si IV 1397 Å, C II 1335 Å, and He II 1640 Å rather than by continuum UV radiation.

Alternatively, synchrotron radiation of accelerated electrons and positrons in the magnetic field amplified at the bow shock could also provide the observed optical and FUV fluxes. As illustrated in Fig. 4, even the observed morphology of the FUV and X-ray emission around J0437-4715 can be understood within such an interpretation.

Modeling of broadband synchrotron emission from PWNe with a single power law particle injection spectrum fails to reproduce the observed spectral energy distributions (e.g., Reynolds and Chevalier, 1984; Atoyan and Aharonian, 1996; Amato et al., 2000). Hence, a spectral break is inferred from the data of radio and X-ray observations of young PWNe in supernova remnants (see, e.g., Gelfand et al., 2009; Bucciantini et al., 2011). Namely, Bucciantini et al. (2011) found that injection of a particle spectrum in the form of a broken power law with 3\times 10^{4}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;\gamma\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;10^{6} and high energy injection index $\geq$ 2.14 can yield a satisfactory description of the available observational data for a number of PWNe. The maximal Lorentz factors of e*±* pairs estimated from the voltage of the total magnetospheric potential of a pulsar are $\sim 10^{8}$ for the spin-down power $\dot{E}\sim$ 6$\times$10${}^{33}\rm~{}erg~{}s^{-1}$.

To model the observed optical, FUV and X-ray emission of the J0437-4715, a power law spectrum of e*±* pairs is injected at the PW termination surface with the index 2.2–2.3 (for the Lorentz factors $\gamma_{0}\leq\gamma\leq\gamma_{\rm m}$). This corresponds to the common models of particle acceleration at relativistic shocks (see, e.g., Achterberg et al., 2001; Ellison and Double, 2004; Keshet and Waxman, 2005; Pelletier et al., 2009; Bykov et al., 2012). While it is not clear that the standard Fermi acceleration might be effective at the PW termination shock, unless highly efficient dissipation of the PW magnetic field occurs before (or at) the termination shock surface, it is true, nevertheless, that this spectral slope is what is inferred from X-ray observations for the majority of PWNe. The shape of the particle spectrum below $\gamma_{0}$ does not need to be specified in the present modeling. The maximal Lorentz factor of the e*±* pairs accelerated at the PW termination surface is likely scaled with the pulsar spin-down power as $\gamma_{\rm m}\propto\dot{E}^{0.5}$ (as the size of the termination shock) and thus may reach (3–5)$\times 10^{7}$ for J0437-4715, i.e., it is just below the value achievable in the magnetospheric potential gap. Pair (re-)acceleration in the colliding shock flow (CSF) of the PW termination shock and the bow shock would result in substantial spectral hardening up to the maximal energy given by Eq. (23). For PSR J0437-4715, with $u_{\rm sh}\sim 100\rm~{}km~{}s^{-1}$ and $R_{\rm cd}\sim 2\times 10^{16}$ cm, from the simplified estimation of Eq. (23) one can obtain $\gamma_{m1}\approx 3\times 10^{6}$. Indeed, the modeled spectra of e*±* pairs re-accelerated in the CSF demonstrate a break at $\gamma\sim\left(2\mbox{--}3\right)\times 10^{6}$ in all the three considered spatial zones (Fig. 3).

In Fig. 3 the accelerated e*±* pair spectra in different zones of a BSPWN for parameters similar to those of the source associated with PSR J0437-4715 are shown. They have been simulated within a Monte-Carlo model, where a power law spectrum $\propto\gamma^{-s}$ with $s$ = 2.3 upto $\gamma_{\rm max}$ = 2$\times$108 (which is \;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;\sqrt{3\dot{E}/2c}) is injected in the immediate downstream of the PW termination surface. As described in §4, the injected pairs propagate through the colliding flow where the efficient re-acceleration reshapes the initial spectrum. The spectra are presented for a few bins, whose locations are shown in Fig. 2. Namely, the red curves correspond to the PW termination shock downstream: the dashed one refers to the $A2$ bin, while the solid one refers to the bin located at the same distance from the pulsar as the $A4$ bin, but lying on the ray directed from the pulsar to observer. The blue curves correspond to the bow shock region (bin $C1$ – the solid curve, bin $C3$ – the dashed curve). The green curves correspond to the region between the two shocks (bin $B1$ – the solid curve, bin $B3$ – the dashed curve). Worse statistics for bins $B1$ and $C1$ is related to the smaller squares of ’detectors’ near the system symmetry axis (see §4.3).

To confront the modeling with the observational data, one needs to construct images and spectra of synchrotron emission produced by the modelled distributions of accelerated e*±*. The intensity of the synchrotron radiation averaged over the magnetic field orientations was obtained using the approximations given by Zirakashvili and Aharonian (2007). The modifications of the synchrotron emission spectra related to the Doppler beaming were taken into account.

Simulated synchrotron images of an object similar to J0437-4715 at $E=$ 8 eV and 1 keV are shown in Fig. 4. On the left panel of Fig. 4 a bright bow-like structure can be seen, which corresponds to the inner part of region “3” in Fig. 2. A very faint counterpart of the bow is hardly seen in this zone on the right panel, where the 1 keV emission is shown. Some details of the simulated source morphology can be affected by varying the parameters of modeling, such as the pulsar axes inclination angles and others, over the allowed parameter space. However, the presence of the bright bow-like structure in the FUV regime and its faint appearance in the X-rays remains.

The spectra of synchrotron emission produced by the re-accelerated pairs in the modelled BSPWN are shown in Fig. 5. The simulated synchrotron spectra are shown for a source similar to J0437-4715. The orientation of the source is given by angles $\Phi=90^{\circ}$, $\Theta=90^{\circ}$ and $\Psi=149.5^{\circ}$ (see §4.4), i.e. the pulsar proper velocity vector lies in the plane of sky perpendicular to the line of sight crossing the pulsar (plane $PoS$) and the position angle $\theta_{pos}=270^{\circ}-\Psi=120.5^{\circ}$. The curves show the spectral energy distribution ($\nu F_{\nu}$) calculated by integration of the spectral emissivity along the lines of sight crossing the bins which were chosen to show the local particle spectra. The colors and styles of the curves are the same as in Fig. 3.

The ratio of the modelled FUV luminosity in the bow-like structure to the modelled X-ray luminosity in the vicinity of the TS is about 2. The modelled luminosity of the J0437-4715 bow in the Hα band is lower than the measured value, so most of the observed Hα emission in the WIYN filter W012 (Brownsberger and Romani, 2014) cannot be attribited to the synchrotron radiation. In general, the simulated synchrotron spectra are consistent with the data of multiband observations of J0437-4715 mentioned above.

The extended (about 3*′′*) FUV flux enhancement (dubbed ”blob” ) which is apparent at the limb of the bow, but not seen in the wide-band optical/near-UV images of Rangelov et al. (2016), can be explained in the frame of the synchrotron model. The optical-UV continuum spectrum is steep at the bow as it is clearly seen in Fig. 5. This means that even relatively modest fluctuations of the magnetic field may produce a sharp structure at higher photon energies, which is not prominent at lower energies. The intermittent clumpy images are a characteristic feature of the systems with steep synchrotron spectra, as it was shown by Bykov et al. (2008).

Another case of a bow shock driven by a supersonically moving pulsar is the nebula of the gamma-ray pulsar J1741-2054, whose estimated velocity is $\sim 150\rm~{}km~{}s^{-1}$ (assuming the uncertain distance of about 400 pc). The Hα bow shock standing at 1\mbox{.!!^{\prime\prime}}5 from the pulsar position was reported by Romani et al. (2010). The pulsar spin-down power is about 9.5$\times 10^{33}\rm~{}erg~{}s^{-1}$. Deep Chandra observations of PSR J1741-2054 performed by Auchettl et al. (2015) have revealed a PWN with an extended (\sim 2\hbox{{}^{\prime}}) trail, whose power law photon index is ${\rm\Gamma}$ = 1.67$\pm$0.06 and the total 0.5–10 keV luminosity is about 2.4$\times$ 10${}^{30}\rm~{}erg~{}s^{-1}$. The trail shows no for synchrotron cooling. Note that the total X-ray luminosity of PSR J1741-2054 is two orders of magnitude higher than that of PSR J0437-4715, though their estimated spin-down powers and proper velocities are comparable. The measured X-ray photon index of J1741-2054 is consistent with the model discussed above, which predicts harder synchrotron spectra in the optical-UV regime followed by asymptotical softening towards the X-ray band. The shape of the simulated energy distribution of the high-energy X-ray emitting particles reproduces the shape of the spectra of e*±* pairs accelerated at a relativistic termination shock. For the power law spectra $\propto\gamma^{-s},$ where $s\sim 2.2$ it gives the photon indices close to 1.6. As the spin-down powers and velocities of the two PSRs are similar, the break energy between the optical-UV component with the hard spectrum produced in the CSF between the termination and bow shocks and that produced at the termination shock of PSR J1741-2054 is expected to be at a few eV. Therefore, the FUV synchrotron luminosity of the bow shock of PSR J1741-2054 is expected to be about a few times 10${}^{28}\rm~{}erg~{}s^{-1}$.

6 Extended X-ray Structures Associated with Bow Shock Pulsars

High resolution observations with the Chandra X-ray observatory (Tananbaum et al., 2014) have revealed a great diversity of X-ray appearances of PWNe (Gaensler and Slane, 2006; Kargaltsev and Pavlov, 2008). For the subset of PWNe with bow shocks (BSPWNe), born due to fast supersonic motion of pulsars, the reasons for diversity are even more numerous than for the rest of the class, as in addition to the pulsar and PW properties (such as pulsar obliquity and wind magnetization), which determine the wind anisotropy (e.g., Bühler and Giomi, 2016), the angle between the star proper velocity and its rotation axis as well as the direction of the local interstellar magnetic field are also important, and greatly extend the parameter space. Some particular cases and configurations will be discussed below.

Peculiar extended X-ray structures associated with PWNe — the “Guitar nebula” and the “Lighthouse nebula” — were observed around two apparently fast supersonically moving pulsars PSR 2224+65 (see, e.g., Romani et al., 1997; Hui and Becker, 2007; Johnson and Wang, 2010; Hui et al., 2012) and IGR J11014-6103 (see, e.g., Pavan et al., 2011, 2014, 2016). Also, Chandra observations revealed a highly extended X-ray feature likely associated with the nearby $\gamma$-ray pulsar PSR J0357+3205 (see, e.g., De Luca et al., 2013), whose space velocity is estimated to be about 390$\rm~{}km~{}s^{-1}$ (assuming the $\sim$ 500 pc distance) and is well aligned with the X-ray feature.

The observed structures of the Guitar and Lighthouse nebulae, as well as the appearance of extended X-ray structures around the Geminga pulsar may be understood within the two jets scenario combined with re-acceleration of the hard-spectrum e*±* component in colliding shock flows (CSFs).

6.1 The Geminga PWN

A bright nearby gamma-ray pulsar, PSR J0633+1746 (“Geminga”) , is located at the estimated distance of 250${}^{+230}_{-80}$pc (Verbiest et al., 2012). Its age is about 340 kyr and its spin-down power is about 3.3$\times 10^{34}\rm~{}erg~{}s^{-1}$. The estimated transverse velocity of Geminga is about 211$\rm~{}km~{}s^{-1}$, which is supersonic for the most widespread ISM phases, and thus it should make a BSPWN. Pulsed X-ray emission up to 20 keV was recently detected from PSR J0633+1746 with the Nuclear Spectroscopic Telescope Array (NuSTAR) by Mori et al. (2014). These authors found that two models can fit the observed spectrum. One of the models consists of a power law component with photon index about 1.7 and two blackbody components of temperatures $kT\sim$ 44 eV and 195 eV. An equally good model consists of a blackbody component of $kT\sim$ 42 eV and a broken power law with photon indexes ${\rm\Gamma}\sim$ 2.0 and 1.4 with a break at 3.4 keV. Much harder photon indexes are needed to model the X-ray spectrum of the Geminga PWN.

With a deep X-ray Multi-Mirror Mission-Newton (XMM-Newton) observation of the PWN, Caraveo et al. (2003) revealed an extended structure with two $\sim~{}2^{\prime}$ elongated, nearly parallel, X-ray tails trailing behind the fast moving pulsar PSR J0633+1746. However, no Hα emission from the nebula was found in a 5-hour exposure of the ESO VLT/FORS1 spectrograph. In addition to the two lateral $\sim$ 0.2 pc long tails, an axial tail of $\sim$ 0.05 pc length resolved into a few individual segments was found in deep Chandra observations by Posselt et al. (2017). The X-ray spectra of the lateral tails, which can be modelled with power law photon indexes ${\rm\Gamma}\approx$ 1, are much harder than that of the axial tail, where ${\rm\Gamma}\approx$ 1.6. The lateral tails are connected to the pulsar, Posselt et al. (2017) found some indications of apparent motions of their footpoints.

The Milagro collaboration (Abdo et al., 2009) reported a detection of a several degrees size TeV emission region in the vicinity of the Geminga PWN. Recently the High Altitude Water Cherenkov Observatory (HAWC) confirmed the TeV source to be about 2 degrees large and found that its TeV spectral index ${\rm\Gamma}\sim$2 (Abeysekara et al., 2017). Contribution of TeV range positrons from the Geminga pulsar is considered as a potential source of the excess over the standard predictions of secondary production in the interstellar matter (see, e.g., Hooper et al., 2017), which will be discussed in detail in section 9.

The observed lack of Hα emission from the PWN may indicate that PSR J0633+1746 is moving through a rarefied ambient medium: an analysis of the geometry of the PWN tails allowed Posselt et al. (2017) to conclude that the ambient number density around the nebula is below 0.01$\rm~{}cm^{-3}$. For the bow shock moving at about 210$\rm~{}km~{}s^{-1}$ such an extremely low density would indicate that the surrounding ISM is hot and thus the Mach number of the shock is moderate. The shock velocity is somewhat higher than that of PSR J0437-4715 and, more importantly, the estimated spin-down power of PSR J0633+1746 is about an order of magnitude higher than that of PSR J0437-4715. Hence, the hard spectrum of e*±* pairs re-accelerated in the colliding shock flow between the bow shock and the PSR wind termination shock of Geminga can be extended to sufficiently high energy to produce synchrotron emission in the X-ray band (according to the maximal energies of the e*±* pairs estimated in §3.1). The X-ray emitting accelerated pairs may flow away through the two magnetic jets and form the apparent lateral tails bent behind due to the pulsar proper motion. Within such a scenario the axial tail may be filled with the pairs accelerated at the PW termination shock and therefore their emission spectrum would have a steeper spectrum, possibly with a photon index close to the observed ${\rm\Gamma}\approx$ 1.6.

An important constraint for the colliding flow model is the spatial distribution of spectral hardness of the observed X-ray emission. The spectra of the “ring” region (Posselt et al., 2017) and of the axial tail are softer than those of the lateral tails. The observed photon indexes of the lateral tails (${\rm\Gamma}\sim$ 1) are substantially harder than those from the other regions of the PWN (${\rm\Gamma}>$ 1.5) as well as that of the pulsar. Assuming a synchrotron origin of the hard emission, one may locate the sources of radiating electrons in the region where the relativistic outflow behind the termination shock is colliding with the bow shock downstream flow. In this region the multi-TeV electrons emitting X-rays can have mean free paths long enough to be accelerated via the Fermi mechanism between the converging flows. This would result in the spectrum of particles $\propto\gamma^{-s}$ with s\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\; 1.

If the pulsar velocity is close to the normal to the “wind plane”, the model of particle acceleration in the colliding flows of the anisotropic wind and of the bow shock downstream could explain the hard X-ray spectral indexes in the lateral tails of the Geminga PWN: here the accelerated electrons are advected through the region of the contact discontinuity between the shocked PW and the interstellar plasma. Numerous simulations of planetary and heliospheric bow shocks predict the enhanced plasma density and magnetic field magnitude in the vicinity of the extended contact discontinuity. The accelerated X-ray radiating electrons would produce the “two tails” configuration as a projection of the rotation figure with a thin radiating shell. The ratio of X-ray surface brightness in the middle of the nebula, at a distance $z$ from the pulsar, to that in the tails would be then $\propto[\delta(z)/R_{\rm cd}(z)]^{1/2},$ where $\delta(z)$ is the radiating shell thickness.

In the frame of the discussed model, the axial tail of the Geminga PWN may be attributed to the synchrotron X-rays from the magnetic tail where the magnetic field reconnection regions and plasma sheets are similar to planetary magnetotails (e.g., Eastwood et al., 2016) An important difference, however, is that, contrary to the planetary and heliospheric magnetotails, reconnecting magnetic flows in PWNe are expected to be composed of moderately magnetized relativistic plasma. Recent two-dimensional particle-in-cell (PIC) simulations of relativistic reconnection in such objects (Sironi et al., 2016) have revealed a phenomenon of formation of quasi-spherical plasmoids, which may show up as blobs apparent in the X-ray images of the the axial tail of Geminga. To be consistent with the observed spectra of the blobs, the simulated spectra of the plasmoids derived by Sironi et al. (2016) must have moderate magnetization parameter 3 $\leq\sigma\leq$ 10.

6.2 Extended X-ray Structures in the Guitar Nebula

The “Guitar nebula” is associated with the bow shock from PSR B2224+65 (see, e.g., Cordes et al., 1993; Romani et al., 1997; Chatterjee and Cordes, 2004). The pulsar is likely one of the fastest pulsars observed with a radio measured proper motion. Harrison et al. (1993) estimated its transverse velocity to be about 986 $\rm~{}km~{}s^{-1}$ at the distance of 1.1 kpc, derived from the dispersion measure. The multi-wavelength study of the nebula by Chatterjee and Cordes (2004) resulted in the distance estimation between 1 and 2 kpc suggesting even faster velocity of PSR B2224+65. The pulsar has a modest spin-down power $\sim$10${}^{33}\rm~{}erg~{}s^{-1}$. Its observed X-ray spectrum is well described by a photon index ${\rm\Gamma}=$1.58${}^{+0.43}_{-0.33}$ (see Hui and Becker, 2007).

The X-ray observations of PSR B2224+65 revealed a jet-like feature extended to \sim 2.4\hbox{{}^{\prime}} (i.e. of about a parsec length) away to the north-west from the pulsar and offset by about 118*∘* from the pulsar direction of motion (see, e.g., Hui and Becker, 2007; Johnson and Wang, 2010; Hui et al., 2012, and the references therein). The observed spectrum of the extended X-ray filament is much harder than of that of the pulsar. Johnson and Wang (2010) derived a photon index of ${\rm\Gamma}=$ 1.00${}^{+0.53}_{-0.47}$ in the extended X-ray filament, while a softer spectrum, with ${\rm\Gamma}=1.69^{+0.39}_{-0.34}$, was observed from the head segment of the PWN. Hui et al. (2012) discussed the possible progenitor star of PSR B2224+65, which is running away from the Cygnus OB9 association. The authors also suggested that the inverse Compton scattering of e*±* pairs of $\gamma\sim 10^{4}$ could produce the extended X-ray feature. The total magnetospheric potential of PSR B2224+65 can provide a maximum particle energy \sqrt{3\dot{E}/2c}\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;7.5$\times 10^{13}$ eV, which is likely below the energy of the X-ray emitting pairs in the observed extended filament, assuming synchrotron emission in the typical ISM magnetic field $\sim$ 3 $\mu G$. We will discuss further the acceleration of the e*±* pairs in the PWN associated with PSR B2224+65.

A scenario to explain the complex phenomenology of the extended X-ray filament was proposed by Bandiera (2008). The author suggested that the highest energy electrons of Lorentz factor $\sim 10^{8}$ are accelerated at the termination shock and then escape from the bow shock into the ISM. They emit synchrotron X-rays in a magnetic field of strength $\sim 45~{}\mu G$ directed along the filament. Such magnetic field is higher than typical ISM values and must have been amplified in some way (e.g., by a streaming instability). The idea behind this scenario is to attribute the observed X-ray emission to collective radiation of accelerated pairs diffusing into the ISM, rather than interpreting the filament as a particle beam. It seems indeed difficult to keep a parsec size approximately linear structure of the beam without a strong bending, despite the high ram pressure due to the fast motion of the pulsar. For example, the X-ray emitting tails of the Geminga pulsar are highly bent, while the estimated proper velocity of PSR B2224+65 is much larger than that of Geminga.

For viability of this scenario a hard spectrum of accelerated pairs escaping the PWN with Lorentz factors $\sim 10^{8}$ is needed. Indeed, the X-ray luminosity of the filament is rather high, corresponding to about a percent of the pulsar spin-down power (if the distance to PSR B2224+65 is about 1 kpc, or even higher for a 2 kpc distance). The photon index of the X-ray filament, estimated from Chandra observations, is also hard: ${\rm\Gamma}=1.00^{+0.53}_{-0.47}$.

Diffusive shock acceleration at the highly relativistic termination shock would provide a pair spectrum of index 2.2-2.3 and corresponding photon index of synchrotron radiation softer than 1.6, which is consistent with that obtained by Johnson and Wang (2010) for the head segment. With the moderate spin-down power derived for PSR B2224+65, the maximal Lorentz factor of e*±* estimated from Eq. (8) is $\sim 10^{8}$. This maximum energy can only be achieved in a transrelativistic outflow with ${\rm\Gamma}_{{\rm flow}}^{2}/{\rm\beta_{\rm flow}}\sim 1$. This excludes as the acceleration sites both the highly relativistic wind before the termination shock and the region around the non-relativistic bow shock. On the other hand, acceleration of pairs up to $\gamma\sim 10^{8}$ may occur in the transrelativistic colliding flow behind the termination shock and the bow shock.

A model of the long linear X-ray filament in PSR B2224+65 also has to explain why it is seen just in the north-west direction from the pulsar. In this context it is important that a similar X-ray filament was reported for the pulsar IGR J11014-6103 by Pavan et al. (2011, 2014, 2016). Recently, Yoon and Heinz (2017) presented hydrodynamic simulations illustrating that the assumption of a density profile of the ISM with inhomogeneities in the form of a series of density discontinuities ahead of the fast moving PSR B2224+65 can reproduce the form of the bow shock and the multiple bubbles seen in Hα observations of the Guitar nebula. Such a model may also help to understand the peculiar appearance of the Hα bow shock of the recycled millisecond pulsar J2124-3358.

6.3 Extended X-ray Structures in the Lighthouse Nebula

IGR J11014-6103 is a source discovered with the IBIS/ISGRI camera aboard the gamma-ray telescope INTEGRAL (Bird et al., 2010). With follow up Chandra observations, Pavan et al. (2011) revealed, in addition to the point-like source, an extended X-ray structure and a helical type tail of \sim 4\hbox{{}^{\prime}} extension. The tail is nearly transverse to the system proper motion (see Fig. 6) and somewhat resembles a similar structure in PSR B2224+65. The authors associated IGR J11014-6103 with the PWN (the ”Lighthouse nebula”) which is produced by the high-velocity ($\sim$1,000 $\rm~{}km~{}s^{-1}$) pulsar of estimated spin-down power $\dot{E}\sim 10^{37}\rm~{}erg~{}s^{-1}$. The pulsar might originate from the supernova remnant MSH 11-61A at an estimated distance of 7$\pm$1 kpc. At this distance the physical extension of the tail would be about 10 pc. The spectral map of the extended filament shows some patchy patterns with spectrum softer than that of PSR B2224+65. If the tail extension is indeed well above 1 pc, the synchrotron burn-off effects may play a role. Deep Chandra observations indicated large deviations from a simple helical model at small and large distances, and the possible presence of an apparent brightness dip at about 50*′′* distance from the pulsar (Pavan et al., 2016). Pavan et al. (2016) concluded that both a ballistic jet scenario and the scenario of Bandiera (2008), which considers synchrotron radiation by high energy pairs propagating along pre-existing interstellar magnetic field lines, can explain only some of the observed features.

Apart from the Guitar and Lighthouse nebulae, several other elongated structures associated with pulsars are known. A very extended (\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;9\hbox{{}^{\prime}}) X-ray filament in the trail of the radio-quiet middle-aged gamma-ray pulsar J0357+3205 was revealed by De Luca et al. (2013). The space velocity of the pulsar was estimated as $\sim 390\rm~{}km~{}s^{-1}$ if it is at a distance of about 500 pc, though no Hα emission was detected. Marelli et al. (2016) found two tail-type nebulae associated with the old pulsar PSR J2055+2539. The pulsar has a modest spin-down power, $\dot{E}\sim 5\times 10^{33}\rm~{}erg~{}s^{-1}$, a characteristic age of about 1.2 Myr, and the estimated distance is $\sim$ 600 pc. The bright tail has an angular size of 12\hbox{{}^{\prime}}\times20*′′* corresponding to $\sim(2.1\times 0.05)$ pc, while the faint one is of 250\hbox{{}^{\prime\prime}}\times 30\hbox{{}^{\prime\prime}} size, i.e. $(0.7\times 0.09)$ pc. The X-ray luminosity of the bright tail in the 0.3-10 keV range is $\sim 2\times 10^{-3}\dot{E}$ (assuming a 600 pc distance). No information on the proper motion of PSR J2055+2539 is available yet, but Marelli et al. (2016) suggested that it is likely to produce a bow-shock.

6.4 Magnetic Jets in BSPWNe

The structure of an astrosphere produced by a relativistic pulsar wind confined by a counter propagating supersonic flow has much in common with that of the heliosphere. Recent progress in the heliosphere observations with Voyager 1 and 2, Cassini, and Interstellar Boundary Explorer satellites, have initiated a new model of the structure of the heliosphere developed by Opher et al. (2015) and Drake et al. (2015). The model suggests that a two-jet structure might provide a better description of the system than the standard comet-like shape. The oppositely directed jets are eventually deflected into the tail region by the motion of the Sun through the ISM. In the models of Drake et al. (2015) and Golikov et al. (2017) the heliosphere is axisymmetric and the structure of the subsonic flows downstream of the termination shock, in the heliosheath and heliopause, is governed by the solar magnetic field. Tension of the solar magnetic field produces a drop in the total pressure between the termination shock and the heliopause. The pressure drop accelerates the plasma flows downstream of the termination shock into the north and south directions to form the collimated jets.

The fact that the magnetic pressure of the azimuthal field may dominate the structure of the subsonic flow in the inner heliosheath was revealed in the early models discussed by Axford (1972). Later, Begelman and Li (1992) proposed that the pinching effect of the toroidal magnetic field in the subsonic flow downstream of the termination shock can be responsible for the observed elongation of the Crab nebula. Similar axisymmetric structures extended in the polar direction were advocated by Chevalier and Luo (1994) for the magnetic shaping of planetary nebulae. Finally, a number of studies (see, e.g., Lyubarsky, 2002; Komissarov and Lyubarsky, 2004; Del Zanna et al., 2004; Porth et al., 2014; Olmi et al., 2016) demonstrated that magnetic collimation of the plasma downstream of a PW termination shock can be a plausible mechanism behind formation of jets in PWNe.

The double jet structure proposed for the heliosphere has not been studied in the regime relevant for supersonically moving PWNe. However, there are reasons to speculate that the tension of the azimuthal magnetic field between the PWN termination shock and the bow shock may drive formation of a double jet structure in this system as well. There is no reason to believe that the orientation of the interstellar magnetic field with respect to the pulsar rotation axis is the same for all BSPWNe. In the case when the plane defined by the directions of the interstellar field and the proper velocity of the pulsar is nearly perpendicular to the pulsar rotation axis, one may expect that the jets would be confined to the region between the apex of the bow shock and the downstream of the termination shock. This constraint may be fulfilled for a narrow subset of bow shock PWNe. For instance, in the Vela PWN (which is likely a BSPWN), shown in the left panel of Fig. 8 the two apparent X-ray jets are connected to the pulsar position. This may be not too surprising given that the bow shock formed in the Vela PWN due to the mildly supersonic flow behind the supernova reverse shock (Chevalier and Reynolds, 2011) is likely quasi-parallel. The available PWN jet models are axially symmetric relative to the pulsar rotation axis, so the scenario described above needs to be confirmed with more complex simulations.

The magnetic jets would channel away some fraction of the high energy particles accelerated by the colliding shock flow in the equatorial wind plane. Magnetically dominated force-free structures are observed in the solar atmosphere and in terrestrial laboratories in the form of magnetic flux ropes (see, e.g., Russell et al., 1990; Birn and Priest, 2007; Daughton et al., 2011; Vinogradov et al., 2016). Magnetic flux ropes typically have helical field structure with the maximum of the (untwisted) field at the rope axis. The magnetic field lines in the jets may reconnect with the interstellar magnetic field lines providing a way for the magnetized ultrarelativistic particles of energies  >$$\sim  10 TeV to escape into the ISM. If the particles are (re)accelerated in CSF, they would have hard spectra of indexes below 2 and therefore contain most of the energy in the highest energy end of their distribution. This is an important point in regard of explaining the observed high efficiencies of the X-ray emission from the extended structures in the Guitar nebula (see, e.g., Johnson and Wang, 2010) and in the Lighthouse nebula (see, e.g, Pavan et al., 2016). The X-ray efficiency may reach a percent of the pulsar spin-down power, a fact that requires (i) hard spectra of accelerated pairs and (ii) an amplified magnetic field in the extended X-ray emitting filaments. We shall discuss below, in §6.5, a possible mechanism of magnetic field amplification by a CR beam propagating along the local magnetic field.

Reconnection of magnetic field lines in the jet with the oppositely directed interstellar field would happen for only one of the two jets associated with the bow shock PWN. This may explain the apparent asymmetry of the X-ray emitting filaments in both the Guitar and Lighthouse nebulae. The X-ray appearance would depend on the orientation of the pulsar rotation axis (providing the jet directions) as well as the direction of the pulsar proper motion. In both the Guitar and the Lighthouse nebulae the pulsar proper motion is directed nearly transverse to the interstellar magnetic field.

6.5 Magnetic Field Amplification in Extended Ballistic Beams

Relativistic particles leaving the accelerator may produce streaming instabilities which would amplify fluctuations of the interstellar magnetic fields to high magnitudes (e.g., Wentzel, 1974; Cesarsky, 1980; Blandford and Eichler, 1987; Berezinskii et al., 1990; Lucek and Bell, 2000).

A fast, non-resonant, almost purely-growing instability driven by CR current was discovered by Bell (2004, 2005). The growth rate for the Bell’s mode propagating along the mean field is

[TABLE]

where

[TABLE]

Here $J_{\rm CR}$ is the electric current carried by the ultrarelativistic particles in the beam escaping the accelerator. Note that in the case of CRs leaving a BSPWN the CR beam is composed of ultrarelativistic e*±* pairs and protons. The highest energy CRs can be accelerated in colliding shock flows (see §3) which energize the highest energy end of the CR spectrum. While the pairs are likely injected into CSF from the PW, the relativistic protons can be injected due to the diffusive shock acceleration at the bow shock. The exact CR composition in the PWNe is somewhat uncertain: both injectures are subject of the intensive studies. However, the CR current $J_{\rm CR}$ exists in the beam leaving the accelerator if the net charge is nonzero in the CR beam.

The saturation level $B_{\rm b}$ of the short scale magnetic field in the beam amplified by Bell’s instability can be estimated from the condition $k_{\rm c}R_{g}\sim$ 1 (Bell, 2004). Then using Eq. (31) one can get

[TABLE]

where $\eta_{\rm CR}$ is the fraction of pulsar spin-down power converted into the accelerated high energy particles forming the beam and $R_{\rm b}$ is the radius of the beam which is R_{g}(\gamma_{\rm max})\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;R_{\rm b}<R_{\rm cd}.

The scale length $l_{\rm b}$ of the beam $\sim c/V_{\rm psr}\times R_{g}$, unless the length is limited by the particle energy losses which is the case for a small $V_{\rm psr}$. The leading edge of the beam is located at a distance $\delta_{\rm bmf}\sim c{\rm\Gamma_{\rm B}}^{-1}$ ahead of the region filled with the amplified magnetic field. The maximal growth rate is ${\rm\Gamma_{\rm B}^{max}}\sim v_{\rm a}k_{\rm c}/2$, and therefore

[TABLE]

The condition $\delta_{\rm bmf}<l_{\rm b}$ must be fulfilled to fill the beam with the magnetic field amplified by the CR driven short scale Bell’s instability. Bell’s instability amplifies highly only the magnetic field fluctuations of short scales $kR_{g}~{}>~{}1$, which are inefficient for scattering the highest energy pairs. The known long-wavelength instabilities (see, e.g., Bykov et al., 2013b) are too slow to be able to fill the beam with long wavelength turbulence. Therefore the high energy pair mean free paths are determined by the background interstellar magnetic fluctuations. The synchrotron-Compton losses of  >$$\sim  10 TeV pairs in the amplified magnetic field can limit the X-ray filament length not to exceed greatly $\sim$ 10 pc.

The mean free path $\Lambda(p)$ of a 10 TeV energy particle along the interstellar magnetic field is a few parsecs (e.g., Strong et al., 2007), while its gyroradius $R_{g}$ is about 1016 cm. Therefore the particles are highly magnetized and their diffusion coefficient transverse to the mean magnetic field estimated from Eq. (18) is $D_{\perp}\approx vR_{g}/3\times(R_{g}/\Lambda)$. The transverse diffusion coefficient is small enough to keep the transverse size of the beam $l_{\rm b\perp}$ to be about a few $R_{g}$.

The main X-ray filament in the Lighthouse nebula extends for about 11 pc at the estimated distance of about 7 kpc (see, e.g., Pavan et al., 2016), while the scale size of the X-ray filament in the Guitar nebula is of about a parsec length (see, e.g., Johnson and Wang, 2010). The X-ray filament length $l_{\rm b}$ in this simplified model scales as $\sqrt{\dot{E}}/V_{\rm psr}$. In the case of the Lighthouse nebula the synchrotron-Compton losses of the X-ray emitting pairs are likely to limit the filament scale size, while in the Guitar nebula, the filament size is roughly consistent with the discussed scaling.

The essential feature of the model described above is formation of a hard spectrum of particles thanks to acceleration in the CSF between the PW termination shock and the bow shock: these particles may power the extended X-ray filaments. The mechanism efficiently converts a sizeable fraction of the wind power into CR luminosity dominated by the highest energy pairs. The accelerated pairs can escape the system through the jets which can be formed in the magnetized region between the shocks. These jets, bent by the ram pressure due to the supersonic motion of the pulsar, can be an explanation of the hard tails observed in the Geminga PWN. If a fast pulsar moves transverse to the interstellar magnetic field, then the reconnection of the ISM magnetic field with the field in a pulsar jet can launch a beam of ultrarelativistic pairs along the interstellar field, which may explain the extended X-ray features observed around PSR B2224+65 and IGR J11014-6103.

In Fig. 7 (right panel) a simulated synchrotron image of the magnetized ballistic beam is presented. A simplified kinetic modeling is performed for a set of parameters likely relevant for the elongated structure of the Guitar PWN. In the discussed scenario only the particles at the highest energies achieved in the colliding shock flow can escape the accelerator (see Fig. 2 in Bykov et al., 2013) to be injected into the beam. The angular distribution of these particles is anisotropic. While the particle pitch-angle distribution is nearly isotropic inside the acceleration region, only particles with pitch-angles $\theta\leq\theta_{c}<\pi/2$ are injected into the beam. The magnetic field amplification by the CR-current driven instability provided $B_{aism}=$ 30 $\mu$G along the beam. An apparent bright region in the central part of the observed structure is related to an inhomogeneity of the field in the magnetic rope. In the presented simulation a circularly polarized field variation with the wavelength 1.3 pc and amplitude $\delta B_{aism}=$ 15 $\mu$G propagating along the beam was imposed onto the background field. The beam was supposed to be directed towards the observer at an angle of 45*∘* relative to the line of sight. The beam had a cylindrical shape and its radius was $3\times 10^{17}$ cm. The distance to the source was $d_{\rm beam}=1.7$ kpc. The particle pitch angle evolution along the beam conserves the adiabatic invariant $p_{\perp}^{2}/B$, since the mean free path of the X-ray emitting pairs is larger than the beam length.

7 BSPWNe Possibly Associated with SNRs

Some of the fast moving young pulsars with PWNe reviewed by Kargaltsev and Pavlov (2008) are located within young supernova remnants (SNRs). Radio observations of the young active pulsar PSR B1853+01 by Frail et al. (1996) revealed that the pulsar, of the spin-down power about 4.3$\times 10^{35}\rm~{}erg~{}s^{-1}$, is located at the apex of an extended feature which has a cometary morphology. The pulsar has a derived transverse velocity of about 375 $\rm~{}km~{}s^{-1}$. Petre et al. (2002) observed in X-rays the pulsar PSR B1853+01 and its associated PWN, embedded into the supernova remnant W44. The authors found the photon index of the extended X-ray nebula to be about 2.2, much steeper than that of PSR B1853+01 which is $\sim$1.4.

A similar cometary shape PWN around the putative pulsar CXOU J061705.3+222127 was observed in the young supernova remnant IC 443 (see, e.g., Bocchino and Bykov, 2001; Swartz et al., 2015). Deep X-ray observations reported by Swartz et al. (2015) revealed a ring-like structure surrounding the pulsar and a jet-like feature. The authors found no evidence for a bow shock or contact discontinuity neither spectrally nor morphologically, since the nearly circular ring is not distorted by motion through the ambient medium. They placed an upper limit on the proper motion of the pulsar of about 310$\rm~{}km~{}s^{-1}$ at a distance of 1.5 kpc. Earlier Bocchino and Bykov (2001) revealed a significant spectral softening along the extended cometary shape X-ray structure, with the power law photon index of about 1.6 in the head and about 2.2 in the tail of the nebula. This is consistent with the expectations based on the simple one zone model of PWNe (see e.g. Chevalier, 2000). However, the softening is not typical for many of the bow shock PWNe which we discuss here.

The PWN energized by PSR B1951+32 is located in the field of the supernova remnant CTB 80 (G69.0+2.7). This is another interesting case of a BSPWN. PSR B1951+32 is an energetic pulsar at a 2 kpc distance with a spin-down power $\dot{E}=3.7\times 10^{36}\rm~{}erg~{}s^{-1}$, an estimated age $\sim 10^{5}$ yrs, and a moderate proper velocity of $240\pm 40\rm~{}km~{}s^{-1}$ (Moon et al., 2004; Li et al., 2005). Radio and optical Hα observations revealed a bow shock and Moon et al. (2004) estimated that the contact discontinuity between the shocked PW and the shocked ambient ISM is at a distance about 0.05 pc from the pulsar.

Spatially resolved X-ray spectra of the main components of the PWN were studied with Chandra by Li et al. (2005). They found that the photon index of the emission from the bow shock region is $1.6^{+0.1}_{-0.2}$ and it is harder than that of the entire nebula, $1.73\pm 0.03$. Moreover, Mattana et al. (2009) reported a significant detection of hard X-ray emission from CTB 80/PSR B1951+32 up to 70 keV with the IBIS/ISGRI camera aboard the INTEGRAL mission.

Li et al. (2005) pointed out that the observed hard emission spectrum may be due either to the acceleration of PW particles at the termination shock or to a contribution from a second $e^{\pm}$ component. They proposed that the new component of emitting $e^{\pm}$ may originate at the bow shock, which was suggested by Hester and Kulkarni (1988) to explain the optical filaments in CTB 80. PSR B1951+32 has spin-down power and proper velocity not very different from those of the Vela pulsar and we will argue in §8 that the hard component may naturally arise due to the acceleration of relativistic $e^{\pm}$ pairs in a CSF between the termination and bow shocks.

An extended (arcminute size) faint X-ray nebula to the south of PSR J0633+0632 was detected by Ray et al. (2011) with Chandra. The pulsar, which is a gamma-ray active Fermi source (see Abdo et al., 2013), was resolved in the X-ray observation, while its optical counterpart is not yet found (Mignani et al., 2016). The X-ray spectrum of the large scale extended PWN has a very hard photon index of 0.74–1.29 for a fixed $N_{H}=0.15\times 10^{22}\rm~{}cm^{-2}$ (0.9${}^{+0.5}_{-0.4}$ for a varied column density). Such a hard spectrum is very difficult to explain as synchrotron emission of the relativistic $e^{\pm}$ accelerated only at the PW termination shock. The radio-quiet pulsar has an estimated spin-down power $\dot{E}\sim 10^{35}\rm~{}erg~{}s^{-1}$ and the pulsar position suggests the Rosette nebula as its possible birth place (see for a discussion Danilenko et al., 2015). In this case the distance to the pulsar would be $\sim 1.5$ kpc and then its proper velocity is $\sim 1,000\rm~{}km~{}s^{-1}$.

Now we will discuss a representative case of a PWN interacting with its SNR. The Vela PWN is a bright nearby object, well studied in multi-wavelength observations. The Vela pulsar is located within its maternal SNR and thus the PWN is interacting with the SNR.

8 The Vela PWN and Vela-like Objects

The bright radio and gamma-ray pulsar PSR B0833–45 (the Vela pulsar) has a period about 89.3 $ms$, characteristic spin-down age of 11,000 yr, and a spin-down power of $\sim 7\times 10^{36}\rm~{}erg~{}s^{-1}$. It is located at a distance of 287${}^{+19}_{-17}$ pc, derived from the parallax measured through the very long baseline interferometry observation by Dodson et al. (2003b). The authors derived the proper transverse space velocity of pulsar PSR B0833–45 to be 61 $\pm 2\rm~{}km~{}s^{-1}$. This velocity is likely subsonic being the pulsar located inside the Vela SNR, where Bocchino et al. (1999) derived from ROSAT observations the presence of a multi-temperature inhomogeneous plasma with a temperature range from about 0.1 to $\sim$ 1 keV and a typical number density 0.1-0.5 $\rm~{}cm^{-3}$.

The nearby Vela supernova remnant appears in multi-wavelength images as a nearly circular outer shell of about 8*∘* in diameter. The bright Vela X region, of about $2^{\circ}\times 3^{\circ}$ scale size, is located near the center of the remnant to the south from the pulsar PSR B0833–45. This region, which has a high degree of radio polarization and a radio spectral index (0.08$\pm$0.10) much flatter than the main SNR shell, is likely powered by the pulsar (see e.g. Weiler and Panagia, 1980; Frail et al., 1997). VLA radio observations of the Vela X region by Frail et al. (1997) revealed a network of fine, extended overlapping linear filaments over the region. The authors also estimated an equipartition magnetic field $\sim 20\mu G$ in the region.

Radio observations of the vicinity of the Vela pulsar by Dodson et al. (2003a) with Australia Telescope Compact Array revealed a highly polarized structure of a few arcminutes extent with two emission lobes to the south and north of the Vela pulsar. While the radio structure surrounds the pulsar, no excess at the position of the X-ray bright nebula was detected. The observed polarization vectors show a symmetry with respect to the spin axis of the Vela pulsar. Dodson et al. (2003a) suggested that the complex emission from the degree scale size Vela X region is the result of the integrated evolutionary history of the Vela pulsar, which makes it different from the few arcminutes vicinity of the pulsar, where emission is dominated by the currently injected particle population (see also the high frequency radio observations of Vela X by Hales et al. (2004)). The observed polarized radio structure was modeled by Chevalier and Reynolds (2011) as an initially thick toroidal structure produced by the PW which is interacting with a mildly supersonic inward flow from the supernova remnant. The flow can be produced in a transient phase of the SNR evolution when the reverse shock reaches the center of the remnant.

The high resolution Chandra X-ray observations have revealed a fascinating synchrotron nebula around the Vela pulsar with a bright structured torus and two jets (Helfand et al., 2001; Pavlov et al., 2001a, b, 2003). The total 1-8 keV luminosity of the bright nebula within 42*′′* from the pulsar is about 6$\times 10^{32}\rm~{}erg~{}s^{-1}$ (see e.g. Pavlov et al., 2003) (assuming a distance of 300 pc) which is a very small fraction ($\sim 10^{-4}$) of the pulsar spin-down power. The photon indexes of the power law fits vary across the nebula being about $1.3\pm 0.1$ for the outer jet of X-ray luminosity about 3$\times 10^{30}\rm~{}erg~{}s^{-1}$, close to 1.1-1.2 for the inner jet and between 1.2 and 1.5 within the bright torii (e.g. Pavlov et al., 2003). A two arcminute scale size northern jet demonstrated helical type dynamical behavior (Durant et al., 2013) and a number of fine scale structures in the nebula were identified by Levenfish et al. (2013) in deep Chandra images of the Vela PWN. The few arcminutes size nebula with polarized radio lobes and bright synchrotron X-ray structures is most likely a result of the recent activity of the Vela pulsar over a few decades time scale after the reverse shock passed through the nebula (see e.g. Chevalier and Reynolds, 2011). A few degrees size extended Vela X region in this scenario is produced by the integrated past activity of the Vela pulsar with account of the relativistic particle escape, their transport and magnetic field structure (see e.g. Aharonian et al., 2006; de Jager et al., 2008; Hinton et al., 2011; Grondin et al., 2013).

Hard X-ray and gamma-ray observations are very useful to study the radiation of the highest energy particles escaped from the Vela PWN. Suzaku/XIS observations just outside the few degrees size radio emitting Vela X region were reported by Katsuda et al. (2011) who found a large scale hard X-ray emission of constant power law photon index $\approx$ 2.4 extending throughout the field of view. While the photon index remains constant, the observed flux is decreasing with increasing distance from the Vela pulsar. Therefore, the authors suggested that the hard X-ray emission is associated with the Vela pulsar wind nebula and that distinct populations of relativistic e*±* pairs are needed to explain both the X-ray PWN and Vela X region.

A wide field hard X-ray image of the Vela X region was obtained by Mattana et al. (2011) with the IBIS/ISGRI camera aboard the INTEGRAL satellite. The image revealed a significant extended emission of $\sim$ a degree scale size both in the northern and southern sides of the region around the Vela pulsar. The spectrum of the Vela nebula in the 18-400 keV energy range was measured by IBIS/ISGRI and SPI on board of the INTEGRAL satellite. Given the point spread function of the coded mask instrument IBIS/ISGRI it was possible to obtain the spectrum within 6\hbox{{}^{\prime}} from the Vela pulsar. The fit to the combined INTEGRAL and Suzaku XIS data yielded a broken power law with the break energy 27$\pm 3$ keV and the photon indexes 1.642${}^{+0.005}_{-0.006}$ and 2.07$\pm 0.05$ below and above the break respectively. The plasma temperature in the hot Vela supernova remnant environment was derived as $kT=0.214^{+0.003}_{-0.005}$ keV. The authors pointed out that the observed change of the slope at $\sim$ 25 keV is compatible with the standard value of 0.5 due to the cooling break in models with a continuous particle injection and radiative losses (Chevalier, 2000). The slope before the break is consistent with the value of photon index 1.6 which is expected in the model of diffusive particle acceleration at relativistic shock. However, this slope is steeper than the photon index 1.3-1.4 observed with Chandra within the bright $\sim$ arcminute size torus-jet nebula around the Vela pulsar. This very hard photon spectrum is indeed reminiscent of those observed at keV photon energies from a number of bow shock PWNe (Kargaltsev and Pavlov, 2008).

The AGILE satellite studied the region and the particle population there in the energy range from 100 MeV to 3 GeV (Pellizzoni et al., 2010). The source AGL J0834-4539 was detected at the significance level about 5$\sigma$ and the spectrum sampled in three energy bands (0.1-0.5 GeV, 0.5-1 GeV, 1-3 GeV) yields a power law fit with a photon index $1.67\pm 0.25$.

Gamma-ray emission from the Vela X region was detected by the Fermi Large Area Telescope, and the spectrum in the 200 MeV to 20 GeV range was measured in the first year of observations (Abdo et al., 2010). With a detailed morphological and spectral analysis of the Vela X region based on four years of observations in the 0.3 - 100 GeV range, Grondin et al. (2013) detected the emission from the northern wing (from the pulsar) of the Vela X region.

The total gamma-ray luminosity of the Vela X region detected by Fermi LAT above 0.2 GeV is 2.4$\times 10^{33}\rm~{}erg~{}s^{-1}$ which is 0.03% of the pulsar spin-down power (Grondin et al., 2013). These authors studied separately the spectra of the northern wing (to the Vela pulsar position) of the Fermi LAT image which can be fitted with a single power law of photon index $2.25\pm 0.07\pm 0.20$. Fitting the southern wing required instead two power laws, with a break at an energy of $2.1\pm 0.5\pm 0.6$ GeV: the photon index before the break was found to be $1.81\pm 0.10\pm 0.24$ while after the break it would steepen to $2.90\pm 0.25\pm 0.07$. Grondin et al. (2013) pointed out that while the energy break at 2.1 GeV is significant, the spectral difference between the northern and southern wings is only marginal. Fermi LAT data analysis by Abdo et al. (2010) favored a scenario with two distinct electron populations to understand the multiwavelength spectrum of Vela X. Fast diffusion of particles accelerated in the vicinity of the pulsar was discussed by Hinton et al. (2011) and Grondin et al. (2013) to explain the origin of the steep Fermi LAT spectrum. Moreover, Hinton et al. (2011) suggested that if the accelerated particles being released early in the evolution of the nebula then travel with a large diffusion coefficient to avoid severe adiabatic losses, this could help to understand the rising positron fraction in the local CR spectrum, an issue which we shall discuss in some detail below.

Very high energy gamma-ray emission from 0.75 to 70 TeV was detected with the groundbased Cherenkov telescope H.E.S.S. from the Vela X region (Aharonian et al., 2006; Abramowski et al., 2012). The source HESS J0835-455 is extended to the south from the pulsar PSR B0833–45 and the morphology of the source was characterized as an ellipse with a Gaussian width along the major axis of $0.52^{\circ}\pm 0.02^{\circ}$, while $0.12^{\circ}\pm 0.02^{\circ}$ is the width along the minor axis. The source shape is different from that obtained by Grondin et al. (2013) with the four years of Fermi LAT data, where the northern wing of the emission was detected. However the source is positionally coincident with the source AGL J0834-4539 imaged by the AGILE telescope in the energy range from 0.1 to 3 GeV. The spectral model of the entire Vela X region (about $1.2^{\circ}$ circle around the central position) in range from 0.75 to 70 TeV is the power law with a photon index $1.32\pm 0.06_{\rm stat}\pm 0.11_{\rm sys}$ and exponential cut off at $14.0\pm 1.6_{\rm stat}\pm 2.6_{\rm sys}$ TeV, while the index is even somewhat harder for the external ring from 0.8 to 1.2 degrees (Abramowski et al., 2012).

Recently Marubini et al. (2015) discussed optical observations of the Vela X region using Digital Sky Survey and SAAO 1.0 m telescope observations in an attempt to constrain the multi-component models of the high energy lepton injection suggested by multiwavelength studies of Vela X.

A PWN that is rather similar to Vela (Chevalier and Reynolds, 2011) is that powered by the pulsar J2229+6114 of spin-down power $\sim 2.2\times 10^{37}\rm~{}erg~{}s^{-1}$ and associated with SNR G106.6+2.9 (also known as the “Boomerang”). The comet-shaped structure of G106.6+2.9 appears as a highly polarized incomplete radio shell with a flat spectrum and an X-ray source of photon index $\sim 1.5$ (Halpern et al., 2001a, b). The authors suggested that the X-ray source may be a wind-blown bubble or bow-shock nebula produced by a young pulsar. A harder X-ray photon index ${\rm\Gamma}=1.3\pm 0.1$ and an X-ray luminosity of $8.2\times 10^{32}\rm~{}erg~{}s^{-1}$ were estimated for the PWN (assumed at a distance of 3 kpc as estimated from X-ray absorption) by Kargaltsev and Pavlov (2008).

A study by Kothes et al. (2001) of the kinematics of the neutral hydrogen and molecular material associated with the SNR has suggested a distance of $\sim$ 0.8 kpc to the SNR and a size of about 14$\times$6 pc. The J2229+6114 PWN would have a diameter of 0.8 pc. This estimate of the distance is supported by the highly polarized emission. On the other hand the distance derived from X-ray absorption (about 3 kpc) and the pulsar dispersion measure ($\sim$ 7 kpc) are much larger, making the distance still very uncertain (Chevalier and Reynolds, 2011).

8.1 A Model of X-ray Emission from the Vela BSPWN

An analysis of the ROSAT HRI X-ray image by Markwardt and Ögelman (1998) suggested that the compact $\sim$ 2 arcmin ”kidney bean” shaped X-ray nebula is due to the bow shock produced by the fast moving Vela pulsar, with an inferred velocity of 260$\rm~{}km~{}s^{-1}$. However, as we referred above, the derived proper transverse space velocity of pulsar PSR B0833–45 is only 61$\pm 2\rm~{}km~{}s^{-1}$ (Dodson et al., 2003b), while the estimated plasma temperature inside the Vela SNR is kT \;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;0.2 keV (though inhomogeneous). Therefore the proper velocity of PSR B0833–45 is too slow to launch a bow shock in such a hot medium. On the other hand the flow behind the SNR reverse shock (see e.g. Truelove and McKee, 1999) may produce the bow shock. Chevalier and Reynolds (2011) noted that there is a transition period before the Sedov solution is reached and this would allow to produce the mildly supersonic flow for the centrally located pulsar PSR B0833–45. They suggested that the reverse shock has recently passed over the Vela pulsar and the plasma flow is mildly supersonic with sonic Mach number ${\cal M}_{\rm s}\sim$ 1.3 at the pulsar position.

As we discussed in §3, the colliding shock flow is expected to occur between the termination surface of the highly relativistic PW and the pulsar bow shock. The magnetized PW is likely highly anisotropic (see e.g. Komissarov and Lyubarsky, 2003; Del Zanna et al., 2004; Kirk et al., 2009; Porth et al., 2014). Relativistic magnetohydrodynamic simulations, supplied with some recipe for the transport and spectral evolution of the accelerated particles, have been extensively used to assess the relation between the wind properties and nebular emission. This has been done through synthetic maps of PWNe, representing their synchrotron emissivity, spectral index and polarization (see e.g. Del Zanna et al., 2006; Volpi et al., 2008; Olmi et al., 2016; Porth et al., 2016).

Recently Bühler and Giomi (2016) performed axisymmetric relativistic magnetohydrodynamic simulations to test the dependence of the PWN morphology appearance on the PW magnetization and obliquity of the pulsar. The simulations suggested that in order to reproduce qualitatively the double rings appearance of the Vela PWN, the pulsar obliquity angle should be about $\pi/4$ and the wind magnetization should be rather high, with $\sigma_{0}\sim 3$ before annihilation of the magnetic field in the striped region, and $\sigma\sim 0.12$ after that (see e.g. Sironi and Spitkovsky, 2011). The authors simulated the PWN synchrotron emission following the evolution of the PW and assuming that the photon emissivity depends on the Doppler factor of the local flow and magnetic field transverse to the line of sight. The origin and spectral evolution of the accelerated particles were not directly modeled in these RMHD simulations: particles were just assumed to be there, due to acceleration either at the termination shock or further downstream in the nebula. The relatively large residual magnetization of the flow in the simulation by Bühler and Giomi (2016) even after the field annihilation region makes diffusive shock acceleration of the pairs at the ultrarelativistic termination shock rather problematic (Sironi et al., 2013), and the maximal energies estimated according to the Eq.(10) are not enough to provide the extended hard X-ray emission observed by Katsuda et al. (2011); Mattana et al. (2011).

Since we aim to discuss here mostly the effect of the bow shock on the Vela PWN high energy emission spectrum, let us assume that a power law distribution of accelerated pairs $\propto Q_{0}\gamma^{-s}$, with $s=2.2$ up to $\gamma_{\rm max}=10^{9}$ (which is below the value estimated from the magnetospheric potential $\sim\sqrt{3\dot{E}/2c}$), is injected at the pulsar wind termination surface. Then we model the particle propagation and re-acceleration in the bow shock PWN through the Monte-Carlo model which described in §4.

Several studies have addressed the propagation of accelerated pairs across a PWN and the consequent variations of spectral index. Tang and Chevalier (2012) concluded that a purely diffusive model can well reproduce the spectral index maps as well as the nebular size for PWNe without an apparent bow shock signature, like the Crab nebula, G21.5-0.9 and 3C 58. This is because in these PWNe one-dimensional advection of the pairs dominates only close to the pulsar, while on larger scales the streamlines become so twisted that the motion can be described as turbulent diffusion. Porth et al. (2016) recently modeled the transport of high-energy accelerated pairs in PWNe using a three-dimensional magnetohydrodynamic simulation of the nebular dynamics and a test-particle propagation model for the accelerated pairs. The authors found that strong fluctuations of the flow velocity resulted in energy independent diffusion in PWNe with P$\grave{e}$clet number of the order of unit. The bow shock flows were not modeled in these papers.

As extensively discussed above, in the case of bow shock PWNe, the MHD flows are important further away from the pulsar and the colliding flows may result in pair re-acceleration causing a spectral hardening. This is clearly seen in Fig. 10 where the local spectra in the different parts of the bow shock PWN are presented. The simulation run is for a magnetic field strength of 200 $\mu G$ immediately downstream of the PW termination shock (section 1 in Fig.2), 120 $\mu G$ at the bow shock (section 3 in Fig.2), 20 $\mu G$ in the region between the two shocks (section 2 in Fig.2). The flow velocity is taken to be at 500 $\rm~{}km~{}s^{-1}$.

In Fig. 10 we illustrate the simulated local spectra of e*±* pairs in different zones of a BSPWN (see §4) for parameters similar to those estimated for the Vela PWN. Spectra are shown for a few bins of the Monte-Carlo model, whose locations are given in Fig. 2. The red curves correspond to the PW termination shock downstream: the solid one refers to the $A2$ bin, the dashed one — to the $A4$ bin. The blue curves correspond to the bow shock region (bin $C2$ – the solid curve, bin $C4$ – the dashed curve), while the green ones — to the region between the two shocks (bin $B2$ – the solid curve, bin $B4$ – the dashed curve). Because of the efficient $e^{\pm}$ pairs re-acceleration in the colliding flows the pairs spectra in the latter two regions are much harder than what was injected at the wind termination surface.

The radiation spectra in the multi-zone model of a Vela-like PWN are shown in Fig. 11. The orientation of the source is given by the angles $\Phi=90^{\circ}$, $\Theta=90^{\circ}$, and $\Psi=270^{\circ}$ (see §4.4), i.e., the pulsar velocity in the local ambient matter rest frame lies in the plane of sky perpendicular to the line of sight crossing the pulsar (plane $PoS$) and is directed to the north. The curves show the spectral energy distribution ($\nu F_{\nu}$) calculated by integration of the spectral emissivity along the lines of sight crossing the same locations. The colors and styles of the curves are the same as in Fig. 10

It should be noted here that the radiation spectra presented in Fig. 11 are simulated for the outer regions of the PWN where the flow velocity typically does not exceed 0.3 c. The effect of particle acceleration in the colliding plasma flows between the downstream of the relativistic termination shock and the bow shock results in a very hard spectrum of the diffuse X-ray emission in the vicinity of the bow shock. The bright X-ray arcs, rings, and knots dominating the emission from the inner part of the nebula in Fig. 8 are likely highly enhanced by the Doppler boosting in the relativistic outflow in the vicinity of the termination shock (see, e.g., Komissarov and Lyubarsky, 2004; Del Zanna et al., 2006; Bühler and Giomi, 2016). Hence, to derive the flux of the e*±* pairs injected into the CSF acceleration from the observed X-ray image shown in Fig. 8 one has to use a quantitative model of the Doppler boosted synchrotron X-ray radiation from the transrelativistic outflow at the wind termination surface. Having this in mind, we show in Fig. 10 the spectra of e*±* pairs simulated for a given arbitrarily normalized spectrum of the pairs, which were injected into the downstream of the termination shock. The absolute flux of the injected pairs (the red curve in Fig. 10) determines the fluxes of the pairs accelerated at the CSF (shown as blue and green curves in Fig. 10). The spectra were simulated within the test particle approximation, which assumes that the pressure of the accelerated particles is below the magnetic pressure at the bow shock.

The spectral shape above a few keV is sensitive to the position of the apex and the magnetic field strength at the bow shock. The spectra presented in the Fig. 11 illustrating the modeling of a Vela-like PWN correspond to a rather large value of the bow shock apex.

The acceleration time to reach in CSF the maximal $e^{\pm}$ energies shown in Fig. 10 (blue and green curves) according to Eq. (21) is about 20–30 yrs for the bow shock upstream flow velocity $\sim 500\rm~{}km~{}s^{-1}$ in a Vela-like PWN.

In Fig. 12 the power law photon spectral indexes obtained by fitting the simulated synchrotron spectra in a wide energy range are presented. In the Chandra ACIS energy range, at both the bow shock and in the region between the termination and bow shocks, we find photon indexes ${\rm\Gamma}<1.5$ consistent with what observed in the arcminute-size X-ray nebula (Helfand et al., 2001; Pavlov et al., 2001a). These indexes are substantially harder than what one would expect for accelerated $e^{\pm}$ pairs injected at the wind termination surface. The magenta line in Fig. 12 corresponds to a position outside the PWN, before the bow shock (the northern region), but inside the Vela SNR where large scale X-ray emission of power law photon index $\approx$ 2.4 was found in the 2-10 keV Suzaku/XIS observations reported by Katsuda et al. (2011).

The red curve in Fig. 11 shows the spectrum of the synchrotron radiation on the lines of sight which are crossing the PW termination shock downstream. The spectrum was simulated assuming a regular magnetic field in the emitting volume. It is likely, however, that the magnetic fields in PWNe are fluctuating. This would result in a particularly strong variability in the steep gamma-ray tail of the spectrum, where modest fluctuations of the magnetic field may lead to strong flares of spectral flux (Bykov et al., 2012). eASTROGAM, a space mission dedicated to the gamma-ray observations with sensitivity $\sim 10^{-12}\rm~{}erg~{}cm^{-2}~{}s^{-1}$ in the energy range 0.2 – 100 MeV, extending up to GeV energies (see, e.g., Tatischeff et al., 2016), can be used to detect or constrain the flux variations in MeV energy regime from the Vela PWN.

It is worth to note that the Vela SNR, being a nearby young SNR hosting a PWN is expected to be a major contributor to the observed leptonic cosmic ray component (Hinton et al., 2011). The spectrum of $e^{\pm}$ injected by the Vela SNR consists of both particles accelerated at the SNR forward shock and pairs accelerated in the Vela PWN which escape from the nebula and its bow shock. The spectral index of the $e^{\pm}$ escaped the BSPWN is likely harder than that produced at the SNR forward shock and extends up to $\sim$ 50 TeV. Therefore the spectrum of $e^{\pm}$ injected from the Vela SNR would have a spectral break at some energy below 50 TeV with a softer spectrum above the break energy. A break at much lower energy might actually be required not to overpredict the flux of CR leptons, as constrained from AMS-02 measurements (Aguilar et al., 2014) and the latest analysis of Fermi/LAT (Di Mauro et al., 2017). We shall discuss more extensively the role of PWNe as the CR positron sources in the Galaxy in §9.

The ultrarelativistic $e^{\pm}$ pair spectra shown in Fig. 10 were produced in the colliding flow behind the bow shock. If the Bohm diffusion model is a good approximation for the transport of pairs in region “3” of the Monte-Carlo model (see §4.2), then, from Eq.(21) the time needed to accelerate the pairs to \gamma\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;3\times 10^{8} in a magnetic field of \;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;10\muG is below 30 yrs. Therefore, within this model the high energy pair spectrum of the Vela PWN is indeed multi-component because of the colliding flow produced by the bow shock. This study supports the interpretation of the Vela PWN as the result of the interaction of the PW with a mildly supersonic MHD inflow (Chevalier and Reynolds, 2011). The flow bypasses the pulsar with speed \;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;500\rm~{}km~{}s^{-1}, as might happen in a transient phase while the reverse shock approaches the center of the progenitor supernova remnant.

The multi-component model outlined above may help to understand the observed X-ray morphology and spectra on the few arcminutes size nebula. We discussed in the §6 the problem of the escape of $e^{\pm}$ pairs as well as protons accelerated at a bow shock PWN. The structure of the large scale magnetic field around the PWN can determine the escape of some of the accelerated particles. In the case of the Guitar nebula and the Lighthouse nebula, as we have argued earlier, the interstellar magnetic field is likely nearly transverse to the pulsar proper motion.

The local magnetic field in the degree scale size region around the Vela PWN is determined by the structure of the MHD flow inside the Vela SNR. If the magnetic field is nearly radial and aligned with the flow velocity behind the SNR reverse shock, then one may understand what observed in the very high energy gamma-rays by H.E.S.S. and shown in Fig. 9. Indeed, pairs with very hard spectra, with particle spectrum index $\sim$ 1.1–1.5, accelerated up to \gamma\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;10^{8} at the colliding flow in the region between the wind termination surface and the bow shock (see the blue curve in Fig. 10) can escape the source. The escaped particles would propagate along the mean field to a few parsec scale distances without appreciable losses (assuming the large scale magnetic field to be below $\sim 50\mu$G). They would produce the TeV range emission by inverse Compton scattering of the background photons. Note that earlier X-ray studies by Bamba et al. (2010) of PWNe of different ages, up to $\sim$100 kyr, but all characterised by low surface brightness, suggested that the accelerated electrons up to $\sim$ 80 TeV can escape from the PWNe without relevant energy losses. Many authors (see, e.g., de Jager et al., 2008; Aharonian et al., 2006; Abdo et al., 2010; Pellizzoni et al., 2010; Grondin et al., 2013; Abramowski et al., 2012), who studied different aspects of the extended Vela X region discussed above, have concluded that multiple populations of accelerated particles are needed to explain the object. In the case of bow shock PWNe, acceleration of pairs in the colliding flow structure can naturally produce an extra particle population, in addition to acceleration at the termination shock.

Some Vela-like objects show hard indexes of the X-ray emission. The pulsar wind nebula G106.6+2.9, with a thick toroidal structure, may have much in common with the Vela PWN, and the scenario discussed above can explain the hard photon index ${\rm\Gamma}=1.3\pm 0.1$ observed in X-rays (Kargaltsev and Pavlov, 2008). Also the PWN of the young pulsar PSR B1823-13 showing properties similar to the Vela pulsar has a rather hard photon index, ${\rm\Gamma}=1.3\pm 0.4$ (see e.g. Kargaltsev and Pavlov, 2008), while the progenitor SNR has not been detected yet (Gaensler et al., 2003). The radio pulsar B1046-58 is associated with an asymmetric PWN, which may imply a supersonic proper motion (Gonzalez et al., 2006). The X-ray emission from this source also shows a hard photon index ${\rm\Gamma}=1\pm 0.2$ (Kargaltsev and Pavlov, 2008). The Vela-like object, PWN G284.0-1.8, associated with pulsar J1016-5857, is characterized by an X-ray spectrum which can be fitted by a power law with index ${\rm\Gamma}=1.32\pm 0.25$ (Camilo et al., 2004). The PWN around the fast-spinning energetic pulsar PSR J0855-4644 likely produced an elongated axisymmetric X-ray nebula with a longitudinal extent of 10*′′*, with a very hard spectrum of photon index ${\rm\Gamma}\sim 1.1$ (Maitra et al., 2016). Moreover, H. E. S. S. Collaboration et al. (2016) showed that at least a part of the shell emission observed by H.E.S.S. from the supernova remnant RX J0852.0-4622 (also dubbed Vela Junior) is likely due to the presence of a PWN around PSR J0855-4644 located near the south-eastern rim of the SNR. The origin and age of PSR J0855-4644 are not yet established (see e.g. Acero et al., 2013a): it can be interacting with the ejecta of RX J0852.0-4622 SNR, likely located at a distance of about 1 kpc and may be associated with the Vel OB1 association (0.8 kpc) (see e.g. Allen et al., 2015). Among the PWNe sources which demonstrated a prominent GeV-TeV emission listed in Kargaltsev et al. (2013); Acero et al. (2013b) only very few have photon indexes of the TeV emission below 2.0. All of these very interesting objects need further multi-wavelength studies to see if there is indeed a subclass of Vela-like PWNe.

9 Positrons from BSPWNe

As already mentioned in previous sections, BSPWNe are likely to play an important role in explaining one of the most intriguing anomalies recently discovered in the field of cosmic ray physics, the so-called positron excess. In the following we briefly review what this excess is, why it is both mysterious and potentially very important, and try to assess what role BSPWNe might play in explaining it.

9.1 The Positron Excess

Galactic Cosmic Rays (CRs) are known to be mainly made of protons (90% by number), with a minor fraction of heavier nuclei (about 10% of Helium in the GeV range and minor percentages of more massive species), a tiny fraction, about 1%, of electrons, and an even smaller fraction of antimatter, positrons and antiprotons. In recent years a few experiments have been devoted to precise measurements of the composition and spectrum of CRs, finding intriguing results on both their hadronic and leptonic components. The most relevant results for the discussion in this section concern the latter and have been first obtained with sufficiently high significance by PAMELA (Adriani et al., 2009) and more recently confirmed with even higher statistics by AMS02 (Aguilar et al., 2013). What these two experiments found is that the ratio between the positron flux ($\Phi_{e^{+}}$) and the total electron ($\Phi_{e^{-}}$) plus positron flux in CRs, $\chi=\Phi_{e^{+}}/(\Phi_{e^{+}}+\Phi_{e^{-}})$, increases with energy for energies above $\sim$ 6 GeV. The results published by AMS-02 (Aguilar et al., 2013) are reported in Fig. 13.

A growing positron fraction is in plain contrast with the standard scenario of CR origin and propagation and indeed, its unequivocal detection by PAMELA (following earlier, but statistically less significant claims by CAPRICE (Boezio et al., 2001), HEAT (Beatty et al., 2004) and AMS-01(AMS-01 Collaboration et al., 2007)), raised considerable interest in the scientific community, one of the main reasons being the fact that the excess positrons might be revealing dark matter related processes.

Before discussing how this excess could be related to the class of sources being discussed in this article, let us briefly review why it appears so puzzling and interesting.

In the standard model for the origin and propagation of Galactic CRs, these particles are accelerated in the blast waves of Supernova explosions and propagate diffusively throughout the Galaxy interacting with the interstellar gas. The SN blast wave picks up the particles to be accelerated from the Interstellar Medium (ISM), and hence only ordinary matter enters the acceleration process. Electrons are accelerated together with protons and heavier nuclei, but positrons (and antiprotons) present in CRs have long been considered to be of purely secondary origin, namely deriving from the nuclear collisions of CR protons with the ISM during their journey through the Galaxy. The most relevant chain of reactions for positron production involves p-p scattering leading to the generation of charge pions. These subsequently decay into muons that finally produce positrons and neutrinos:

[TABLE]

A similar reaction involving negatively charged pions lead to electron production. In the standard scenario, all these secondary leptons are much less numerous than primary electrons. In addition, their ratio to primaries is expected to monotonically decrease with increasing energy. The reason for this expectation is more easily understood in the case of the ratio between secondaries and primaries that do not suffer energy losses, while the actual demonstration in the case of leptons that do suffer appreciable energy losses in the relevant energy range is more cumbersome.

Let us therefore consider, just for illustrative purposes, the ratio between secondary and primary nuclei that do not suffer relevant energy losses during propagation. In relation to CR confinement, let us describe the Galaxy as a cylinder of radius $R_{d}\approx 15$ kpc (the star and gas disk radius) and height $H\approx 3$ kpc (as the most recent propagation models suggest for the extent of the magnetized halo of the Galaxy, see e.g. Trotta et al. (2011)). Let us further assume that CRs diffuse through the volume of the Galaxy with a spatially uniform diffusion coefficient $D(E)\propto E^{\delta}$, where $E$ is the CR energy and $\delta>0$. Finally, particles are considered free to leave the Galaxy through the vertical bounds at $\pm H$. This is the so-called Leaky Box Scenario for CR propagation.

In the absence of energy losses, the equilibrium spectrum of primary CRs, $N_{\rm prim}(E)$, will then simply result from balance between injection and escape from the Galaxy

[TABLE]

where we have taken the number of particles injected by the sources per unit time and energy interval to scale with energy as $Q_{\rm prim}(E)\propto E^{-\gamma_{\rm inj}}$.

Following the same reasoning we can also derive the equilibrium spectrum of secondaries. Now the injection relies on the spallation reactions of primary nuclei and is then described by

[TABLE]

where $\sigma_{\rm sp}$ is the spallation cross-section (assumed to be independent of energy, which is a good approximation in the GeV-TeV energy range), $c$ is the speed of light and $n_{H}$ is the density of target material.

As a result, the equilibrium spectrum of secondaries will be:

[TABLE]

and the secondary-to-primary flux ratio will scale as:

[TABLE]

While absence of energy losses is a good approximation for protons and nuclei of all energies between a few GeV and the knee, the same is not true for leptons, which lose energy through synchrotron and Inverse Compton radiation while propagating through the Galaxy. If we take the standard value of magnetic and radiation energy density in the Galaxy, $U_{B}\approx U_{\rm rad}\approx 1$ eV cm*-3*, the radiative loss time, $\tau_{\rm loss}(E)\approx 2.5\times 10^{8}$ yr $(E/{\rm GeV})^{-1}$, is shorter than the propagation time $\tau_{\rm esc}(E)\approx 2\times 10^{8}$ yr $(E/{\rm GeV})^{\delta}$ already at $E>$ few GeV, for $0.2<\delta<0.8$, as inferred from recent studies of CR propagation (Trotta et al., 2011). In this case, in the energy range of our interest, the equilibrium spectrum of electrons, which we approximate as all primaries, is described by:

[TABLE]

where the first term on the rhs comes from injection and the second from propagation. The form of the latter is determined by the combination of diffusive propagation and radiative energy losses that limit the number of contributing sources with respect to the case of nuclei, imposing constraints in terms of distance in space and time.

Analogously, for positrons, that are all secondaries, we find:

[TABLE]

where what appears in the first term of the rhs is the equilibrium proton spectrum ($\propto E^{-\gamma_{p}}$). The result for the positron fraction is then:

[TABLE]

and in the standard scenario, with $\gamma_{\rm inj,p}=\gamma_{\rm inj,el}$, $\chi\propto E^{-\delta}$ is expected to monotonically decrease with energy. This is why the observed increase was found to be puzzling and stimulated a plethora of attempts at explaining it.

Before discussing the possible explanations and how BSPWNe enter the scene, an important remark is in order. The prediction of a decreasing positron fraction is a strong one only in a scenario that assumes a spatially uniform distribution of sources. This is indeed the assumption underlying Eqs. 39 and 40, but we know that this assumption does not apply to SNRs, the putative sources of primary electrons: these are discrete sources in space and time and their distribution in the Galaxy is not uniform at all. The effects of the discreteness of the sources and their rapidly decreasing density with Galactocentric distance were analysed for CR nuclei by Blasi and Amato (2012a), and in that case they were found to be poor: the resulting spectrum of all different species still respected on average the predictions of the simplified leaky box model (LBM hereafter). With electrons, however, due to radiative losses, some difference might appear, since especially at the highest energies the number of sources that might contribute to the flux at Earth is likely to be very small, and strongly dependent on the recent SN explosion history in the solar neighborhood. A possibility would then be that it is actually the electron flux that steepens with energy with respect to the LBM prediction, rather than a flattening occurring in the positron flux. This possibility has been ruled out by more recent analysis of the PAMELA (Adriani et al., 2011) and AMS-02 data (Aguilar et al., 2014), which clearly show a flattening of the positron flux, a real positron excess. The spectra of electrons and positrons obtained by AMS-02 (Aguilar et al., 2014) are shown in Fig. 14.

9.2 Qualitative Evaluation of the PWN Contribution

An intriguing suggestion immediately made after this discovery was that these extra positrons could come from dark matter annihilation in the Galaxy (see e.g. Serpico (2012); Di Mauro et al. (2016) for a review). While a dark matter related origin of the rise of $\chi$ is difficult to rule out completely, astrophysical explanations seem favored at the current time and the most natural place where to look for a source of extra positrons is the magnetosphere of a pulsar, where a large number of electron-positron pairs are created by every electron that leaves the star surface.

These $e^{-}-e^{+}$ pairs end up in the PW and are then accelerated at its termination shock. The efficiency of conversion of the pulsar spin-down power into particle acceleration is usually thought to be very high, and measured to reach about 30% in the case of the Crab nebula, the prototype PWN. In addition, we know from observations of many such objects that the particle spectrum is usually a broken power law: from radio observations we infer a particle spectral index $\gamma_{R}$ in the range $1<\gamma_{R}<2$, while the X-rays come from particles with spectral index $\gamma_{X}>2$. The boundary between the two regimes is inferred to occur in the several hundreds GeV range. It is apparent then that these are extremely promising candidates to explain the excess at least from the spectral point of view: the excess is observed in the energy region between few and 300 GeV (no data are available at higher energies) and in this interval PWNe could well enter the scene with their flat positron spectrum.

Of course as long as the particles are part of the PWN they are confined within the remnant of the parent supernova explosion. Their Larmor radii in the nebular magnetic field (of typical strength around a few $\times$100 $\mu G$) are very small and diffusive escape is not effective. The situation is different when the pulsar breaks out of the remnant and finds itself in the ISM. The time at which this happens and the energy that is left in the system at that time were derived in § 1, where we found, for typical values of the pulsar and ISM parameters, $T_{\rm esc}\approx few\times 10^{4}$ yr and $E_{\rm res}\approx 10^{47}$ erg.

After time $T_{\rm esc}$ the supersonic motion of the pulsar through the ISM transforms its PWN in a Bow Shock PWN. As we extensively discussed in earlier sections, the plasma in front of the pulsar is confined by the bow shock, but particles flowing through the tail are free to leave the system and become part of the cosmic ray pool.

These particles are the interesting ones in view of explaining the positron excess. Their spectrum, at the energies of interest, is shown to be flat, similar to what described above for their younger counterparts. In particular, radio emission again shows a flat particle spectrum in the up to hundreds of GeV. In addition, numerical simulations support the idea that these low energy particles can easily escape from the system along the tail with very little energy losses (Bucciantini et al., 2005).

These particles are easily shown to be numerous enough to explain the excess. This can be done through the following back-of-the-envelope calculation.

Let us go back to the rough description of the electron spectrum in Eq. (39), and let us assume that these particles are accelerated in SNRs with the same spectrum as protons: given the protons’ spectrum measured at Earth (Eq. 35) $\gamma_{\rm inj,p}+\delta\approx 2.7$, one can promptly derive the spectrum that must be injected by the sources as soon as the energy dependence of the diffusion coefficient is known. For the latter we assume a slope $\delta=1/3$, in agreement with a Kolmogorov’s model for the turbulence in the Galaxy (which is preferred to steeper energy dependencies by considerations concerning the CR anisotropy (Blasi and Amato, 2012b)). We then obtain that electrons must be injected with $\gamma_{\rm inj,el}\approx 2.4$. The fraction of SN explosion energy that is usually assumed to be converted into CR acceleration is $\xi_{\rm CR}=0.1$ and most of it is carried by protons. The ratio between accelerated electrons and protons that is needed to interpret the CR measurements at Earth is usually $f_{\rm e,p}\approx 0.01$. We can use this information to find the number of electrons that SNRs inject in the Galaxy per unit time and unit energy interval:

[TABLE]

with $\cal R$ the SN rate in the Galaxy and $E_{0}=1$ GeV corresponding to the particles that carry most of the energy in the distribution. If we now follow a similar line of reasoning for the positrons, assuming that they are injected by BSPWNe, which convert some fraction $\eta$ of $E_{\rm res}$ in a distribution of electrons and positrons characterized by a flat spectrum $E^{-\gamma_{pos}}$ with $\gamma_{pos}<2$, extending up to $E_{\rm max}\approx 500$ GeV, we can write:

[TABLE]

where $E_{\rm res}$ is given in Eq. (4) and $f_{\rm II}$ takes into account the fact that only $\approx$ 80% of the SN explosions give birth to a pulsar. Since propagation affects electrons and positrons in the same way, the ratio $\chi$ will simply scale as the injection, and for energies large enough that the secondary positrons can be neglected, will read

[TABLE]

where the last approximate equality holds for the typical values of the parameters mentioned above. At 30 GeV this value of $\chi$ is within a factor of 2 of the AMS-02 measurement reported in Fig. 13, which strongly suggest that the contribution of PWNe must definitely be taken into account when wondering about the positron excess. In addition the energy dependence in Eq. (44) also approximately agrees with the data.

9.3 Positrons from BSPWNe

A quantitative calculation of the flux of electrons and positrons expected from BSPWNe was carried out by Blasi and Amato (2011), who showed that BSPWNe are likely to contribute a large fraction of the needed extra positrons, if not all of them. That work was adopting some simplifying assumptions about the spatial distribution of the SN explosions in the Galaxy (the presence of arms was not considered) and was associating to every explosion a universal pulsar with typical values of the parameters, rather than accounting for the spread of parameter across the pulsar population. In the following we show the results of a more refined version of the same exercise, including the above mentioned missing ingredients and the most recent information available from AMS-02 on the spectrum of primary protons (Aguilar et al., 2015) and on the diffusion coefficient (Aguilar et al., 2016). We first provide the analytical expression for the contribution at Earth of a single source and then sum up the contribution of all the potentially relevant sources in the Galaxy, generated at random from an appropriate distribution. Some other studies along the same line were performed by Di Mauro et al. (2014), who considered the contribution of putative nebulae associated with pulsars from the ATNF catalogue, rather than a generic pulsar population.

The equation describing the transport of particles in the presence of spatial diffusion, energy losses and catastrophic losses reads:

[TABLE]

where $f$ is the number of particles of energy $E$ per unit volume and energy interval at position $\vec{r}$ and time $t$, $D(E)$ is the particle diffusion coefficient, $B(E)=dE/dt$ represents the particles’ energy losses, $T$ is the timescale for catastrophic losses (particle destruction) and $Q(E,\vec{r},t)$ is a source term (number of particles with energy $E$ injected in the system per unit time and energy interval at position $\vec{r}$ and time $t$ ). The solution of Eq. (45) can be found through the Green function method (Syrovatskii, 1959), and allows one to compute the particle distribution function at any location once $D(E)$, $B(E)$, $T$ and $Q(E,\vec{r},t)$ are known. In the absence of catastrophic losses, the Green function reads:

[TABLE]

where $\delta$ is the Dirac delta function, $\tau$ is the timescale of energy losses, and $\lambda$ is the particle diffusion length. In formulae we have:

[TABLE]

The full solution of Eq. (45) can be found through convolution of the Green function in Eq. (46) with the appropriate source term. Sources such as SNRs and BSPWNe are localised in space and time. The injection by one such source, located at position $\vec{r}_{s}$ at time $t_{s}$, can be described as

[TABLE]

Convolution of the Green function in Eq. (46) with the source term in Eq. (48) provides an expression for the particle density per unit volume and unit energy interval contributed at $\vec{r}=(x,y,z)$ by a single source of space-time coordinates $(x_{s},y_{s},z_{s};t_{s})$:

[TABLE]

where

[TABLE]

In Eq. (49) $z_{n}^{\prime}=(-1)^{n}z_{s}+2nH$ and the sum on the rhs is introduced to take into account the boundary condition of particle escape from $\pm H$ (see e.g. Blasi and Amato (2012a) for a discussion).

Eq. (49) can be used with any distribution of sources. In particular, when a spatially uniform distribution of sources is considered, it is easy to see that its integration over the source volume leads to a particle spectrum as described by Eq. (39), which we used for our qualitative discussion.

However for a quantitative evaluation of the electron and positron spectrum contributed by SNRs and PWNe the procedure to follow is different. SN explosions must be generated at random places in the Galaxy and at random times. In order to be sure to correctly reproduce the spectrum of low energy CRs, one needs to include all the explosions occurring in a time span corresponding to the escape-time of these particles. To do this, we assume (as we always will, hereafter) a Kolmogorov diffusion coefficient, in agreement with the the estimate by the AMS-02 team based on the B/C ratio (Aguilar et al., 2016). The normalisation is taken according to the most recent GALPROP analysis (Trotta et al., 2011), resulting into:

[TABLE]

where $H_{0}=3$ kpc is the halo size. The escape time of 10 GeV particles is then $\tau_{\rm esc}(10GeV)=H_{0}^{2}/D({\rm 10GeV})$=45 Myr. With an explosion rate as described by Faucher-Giguère and Kaspi (2006), namely a gaussian distribution around an average rate of $2.8$ explosions/century and a standard deviation of $0.1$ events/century, the number of explosions that have to be taken into account is $\approx 1.2\times 10^{6}$. The explosions are assumed to occur at random places in the Galaxy, spatially distributed as:

[TABLE]

where a cylindrical coordinate system has been adopted, with $\rho_{s}$ the radial coordinate of the explosion and $z_{s}$ its height above or below the Galactic plane. The parameters appearing in Eq. (52) have the following values: the sun radial distance from the Galactic centre is $R_{\odot}=8.5$ kpc, $h_{\rm SN}=50$ pc (form Faucher-Giguère and Kaspi (2006)), $\alpha=1.09$ and $\beta=3.87$ from Green (2015). The presence of Galactic arms has been included on top of this profile, according to the description given by Faucher-Giguère and Kaspi (2006). The resulting distribution of SN explosions can be seen in Fig. 15, where the coordinates of about 30000 events are reported in a Galactocentric frame.

All these events are considered as actual contributors to the primary electron flux at Earth, while only 80% of them, randomly selected, are assumed to be associated to type II explosions giving origin to a pulsar and hence included in the count of positron contributors. The flux of electrons and positrons is obtained as the sum over all relevant sources of the single source contribution described in Eq. (49), where $N_{s}(E)=Q_{s}(E)/{\cal R}$ with $Q_{s}(E)$ given by Eqs. 42 and 43 for electrons and positrons respectively. The electrons’ injection spectral index is assumed to be the same as that of protons, and the spectrum of the latter is taken equal to the fit to the newest data provided by the AMS-02 collaboration (Aguilar et al., 2015). The SN explosion energy that enters the normalisation is assumed to be $E_{\rm SN}=10^{51}$ erg for all events. As far as the positrons are concerned, observations of PWNe suggest that $\gamma_{pos}$ (and accompanying $\gamma_{el}$) is different from source to source but without an obvious relation to any other observable. The results that are shown in the following refer to the choice $\gamma_{pos}=1.2$ up to a maximal energy $E_{\rm max}=500$ GeV. A steeper power law with index $2.3$ is assumed at higher energies.

$E_{\rm res}$ will also vary from source to source depending both on the pulsar properties and on the environment in which the explosion occurs. In the absence of a good and easy recipe for how to vary the latter, previous work can be improved in realism by taking into account variations of the former. The surface magnetic field is assumed to be the same for all pulsars and corresponding to the average value proposed by Faucher-Giguère and Kaspi (2006), $B_{*}=4.5\times 10^{12}$, but the pulsar’s spin period at birth and the proper motion vary among different objects according to the distributions provided by the same authors.

In addition to primary positrons from BSPWNe and primary electrons from the same sources plus SNRs, we also account for secondary electrons, generated from CR interactions with the ISM. The density of the latter is taken to correspond to 1 particle/cm3 and the incident proton spectrum is taken again from AMS-02 (Aguilar et al., 2015).

In Fig. 16 the spectrum of electrons and positrons resulting from this calculation is shown together with the AMS-02 data (Aguilar et al., 2014). Some of the individual curves, corresponding to different realisations, are also shown in order to give an idea of how different the electrons’ and positrons’ spectra can be at high energy depending on the source configuration. A conversion efficiency of $\eta=$22% of the pulsar spin-down energy into acceleration of pairs has been adopted.

The difference between the different realisations at high energies depends on the strong decrease with increasing electron energy of the number of contributing sources. Fig. 17 serves the purpose of highlighting this effect. The two panels of the figure show the sources that in this calculation actually make a contribution to the flux of positrons (and electrons) of different energies at Earth. It is apparent that while at 10 GeV the number of relevant sources (green diamonds) is large (more than 2000 in the particular realisation of the Galaxy to which the plots refer), it strongly decreases with increasing energy: at 100 GeV about 300 sources are relevant (black triangles in the plots), and only a handful at 1 TeV (blue asterisks). The panel on the left shows the projected position in the Galactic plane. The fact that the Earth is in an inter-arm region, poorly populated with sources is evident from the plot and plays an important role: this is likely the main reason for the need of a larger value of $\eta$ with respect to the calculation by Blasi and Amato (2011).

From the left panel of Fig. 17 one sees instead what a small range of distances in space and time can contribute to positrons at the highest energies: already at 100 GeV all the leptons that reach Earth come from sources closer than 2 kpc and younger than 2 Myr.

In Fig. 18, finally, the rise of the positron fraction with energy that this model provides is compared with AMS-02 data. The different curves are drawn following the same notation as in Fig. 16. A generally good agreement is seen between the model and the data.

Further refinement of this kind of calculation are possible and desirable, possibly with inclusion of some model for the density distribution of the ISM. However the main point this section aims at making is that bow shock PWNe are expected to play an important role as contributors of cosmic ray positrons in the Galaxy. That their contribution to the positron excess be non-negligible seems unavoidable and this must be properly assessed and taken into account before using positrons as a diagnostics of other phenomena. In this sense, achieving a better understanding of the physics of these sources is a task with potential far reaching implications for cosmic ray and dark matter physics.

Acknowledgements.

The authors thank the referee for constructive comments. A.M. Bykov and E. Amato thank the staff of ISSI for their generous hospitality and assistance. A.M. Bykov and A.E. Petrov were supported by the RSF grant 16-12-10225. A.B. thanks A.V. Artemyev, S.M. Osipov, and G.G. Pavlov for discussions. Some of the modeling was performed at the “Tornado” subsystem of the St. Petersburg Polytechnic University supercomputing center.