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

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

Abstract

Pulsars with high spin-down power produce relativistic winds radiating a non-negligible fraction of this power over the whole electromagnetic 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 ambient interstellar medium (ISM). Some of the PWNe are moving relative to the ambient ISM with supersonic speeds producing bow shocks. In this case, the ultrarelativistic particles accelerated at the termination surface of the pulsar wind may undergo reacceleration in the converging flow system formed by the plasma outflowing from the wind termination shock and the plasma inflowing from the bow shock. The presence of magnetic perturbations in the flow, produced by instabilities induced by the accelerated particles themselves, is essential for the process to work. A generic outcome of this type of reacceleration is the creation of particle distributions with very hard spectra, such as are indeed required to explain the observed spectra of synchrotron radiation with photon indices 1.5. The presence of this hard spectral component is specific to PWNe with bow shocks (BSPWNe). The accelerated particles, mainly electrons and positrons, may end up containing a substantial fraction of the shock ram pressure. In addition, for typical ISM and pulsar parameters, the released by these systems in the Galaxy are numerous enough to contribute a substantial fraction of the positrons detected as cosmic ray (CR) particles above few tens of GeV and up to several hundred GeV. The escape of ultrarelativistic particles from a BSPWN — and hence, its appearance in the far-UV and X-ray bands — is determined by the relative 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 both on pulsar winds and on the ISM. We discuss the interpretation of imaging observations in the context of the model outlined above and estimate the BSPWN contribution to the positron flux observed at the Earth.

11

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 and inferred that 15% of the pulsars have natal velocities greater than 1000. 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 15210 for the normal pulsars and 546 for the recycled ones. The corresponding mean 2D velocities in the catalog are 24622 and 8713, with the highest values inferred for the pulsars B2224+64 and B2011+38. The mean 3D natal pulsar velocity was derived to be 40040.

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 at which this occurs can be estimated as , where is the velocity of pulsar proper motion, and are the time and remnant radius at the beginning of the Sedov-Taylor phase. With an average estimated velocity of about 400, an estimate of the time at which a pulsar escapes from its SNR can be found as

(1)

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 , with – the pulsar braking index, and – the star rotation frequency and its time derivative, respectively. For dipole spin-down one expects , since , where g cm is the star momentum of inertia and – 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 when it leaves the remnant. From the equation:

(2)

where km is the neutron star radius, we find:

(3)

and – the initial star rotation frequency ( is the initial star spin period). Finally, the residual rotation energy of the star at is:

(4)

The average value of 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: G and ms. Using the estimate of from Eq. (1) one obtains  erg and  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  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 (at the bow apex it is about the pulsar speed) is

(5)

where is the ionized ambient gas number density measured in , is the local magnetic field in the shock upstream measured in G and is the shock velocity in 10. The sonic Mach number for a shock propagating in the interstellar plasma of the standard (Solar) abundance is

(6)

where is the plasma ion temperature measured in 10 K and is the ratio between the electron and ion temperatures.

Consider a pulsar moving supersonically with a velocity relative to the local rest frame of the ambient medium. The standoff distance of the PW termination shock in this frame can be estimated as

(7)

where is the pulsar spin-down power, is the ambient ISM density. In fact, is rather the location of the contact discontinuity, than either of the PWN termination shock, or of the bow shock. The parameter depends on the magnetic inclination angle of the pulsar; 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 and, therefore, , which scales . 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 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 (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 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 is characterized by a velocity 0.75 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 (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 in a Poynting flux dominated source of magnetic luminosity . According to the analysis of Norman et al. (1995); Waxman (1995); Lemoine and Waxman (2009), the luminosity required to reach the energy (in the absence of energy losses) is

(8)

where is the flow speed divided by the speed of light , and 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:

(9)

This approximate relation indicates that the most favourable sites to achieve the highest energy are indeed in the trans- or sub-relativistic flows of 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 10 (and even lower, 310, 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

(10)

where is the total particle flux entering the nebula, which is somewhat uncertain and is likely within the range 3 1010. 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 0.3. Several BSPWNe harbor regions of extended X-ray emission, where the magnetic field is estimated to be a few tens of G. In order to produce X-ray synchrotron photons in such a field, pairs with a Lorentz factor 10 are required. But according to Eq. (8) such high energies are compatible with the amount of spin-down power mentioned above only if 3 while 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 10 is relaxed to 10.

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 of a particle of charge accelerated by a strong quasi-steady bow shock of velocity and characteristic radius is limited by the shock curvature. It can be roughly estimated from

(11)

where 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 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):

(12)

where is the energy density of the CRs accelerated at the shock. Then the amplitude of the fluctuating magnetic field is

(13)

where is the fraction of shock ram pressure converted into the energy of the accelerated particles, is the shock velocity in 10 cm s and is the number density of the ISM. This amplitude may reach , where is the mean magnetic field in the ISM. In any case, in the vicinity of a bow shock, is likely well above the estimated amplitude of the field fluctuations in the background ISM turbulence at the scales below 10 cm, which are important for the scattering of TeV range particles.

The mean free path 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 cm, where the index 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 pairs accelerated in the PWNe within a region of scale size as given by Eq. (7), for a bow shock speed above 100 . Depending on the bow shock velocity, strong MFA by CR-driven instabilities at the bow shock would affect the propagation and confinement of relativistic pairs up to Lorentz factors 10.

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 . At higher shock velocities the energy density of the fluctuating magnetic field scales as a fraction of the shock ram pressure. The MFA efficiency is illustrated in Fig. 1 where the ratio 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

(14)

The validity of the Bohm diffusion approach requires strong scattering of particles by quasi-resonant fluctuations, i.e., . In such a case the diffusion coefficient where 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 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 which satisfy the condition 1, where is the gyroradius of a particle.

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

(15)

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

(16)

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 which is much longer than the gyroradii of protons of energy 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 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 with for e pairs up to a Lorentz factor determined by the condition , where is the particle gyroradius in the mean magnetic field .

Accurate estimations of the scale-length are not available yet because of the limited dynamical range of numerical simulations of the CR-driven turbulence. Therefore, we consider as a free parameter which determines the Lorentz factor 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 , the mean free path produced by the short-scale modes can be approximated as

(17)

where . 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 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 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 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 is large compared to the gyroradius , the transverse diffusion coefficient is much smaller than that along the mean field :

(18)

When one obtains . 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.

Figure 1: The magnetic turbulence energy flux in the shock downstream as a function of shock velocity simulated with a nonlinear Monte-Carlo model of diffusive shock acceleration (Bykov et al., 2014).

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 and aligned with the positive and negative directions of axis and occupying half-spaces and , respectively, are separated by a discontinuity in the plane . 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 which depend on the particle Lorentz factor .

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 .

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

(19)

where is the standard Heaviside step function, is the maximal Lorentz factor of an accelerated particle, and is a normalization constant. The acceleration time of the highest energy e pair is

(20)

We take , in accordance with the typical flow velocities downstream of the PW termination shock and at the bow shock, respectively, while () is the Bohm diffusion coefficient in the half-space . Here is the particle gyroradius in the magnetic field amplified by the CR-driven instabilities. In the limit Eq. (20) simplifies to

(21)

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 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 satisfying the equation

(22)

Therefore, assuming the Bohm diffusion in the accelerator, the upper limit on can be obtained from

(23)

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, , should be longer than the acceleration time. This provides another constraint on the maximal Lorentz factor given by the equation . Using a simplified approximation for 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 :

one obtains

(24)

Then the maximal Lorentz factor of the pairs accelerated in a CSF can be estimated as .

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 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 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, with 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, 11.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.

Figure 2: The geometry of the model: half of the axial section of the system by the plane perpendicular to the line of sight, crossing the pulsar (plane ). The numbers 0 – 4 label the regions with different diffusion modes: 0 is the cold PW, 1 (between the orange curves) is the shocked PW, 2 is the zone near the contact discontinuity, 3 (between the red curves) is the postshock flow of the ISM matter, 4 is the unperturbed ISM. The relative sizes of the regions are not to the scale and depend on the particular properties of a BSPWN. The locations (bins) chosen in the following sections to illustrate simulated particle energy distributions are marked by labelled black squares (). The bins are located at the midpoints of the dashed gray segments of the lines connecting the pulsar and the surfaces, which separate regions 0 – 4 (see §4.2). The gray lines are directed at 0, 30, 60, and 90 with respect to the system symmetry axis.

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 , , and are centered at the pulsar position and the zenith direction is aligned with the pulsar velocity vector. Also the cylindrical coordinates (, , ) orientated along the symmetry axis of the system () are used.

The innermost spherical region with radius 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:

(25)

The scaling factor is a model parameter. The shocked PW flow velocity is directed radially and its value is at distance 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:

(26)

Here is the standoff distance defined by Eq. (7). The positions of inner and outer boundaries of region “3” are determined by the model parameter . The apexes of these surfaces are denoted as and , 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 , where is the coordinate of the tail-oriented base of the cylinder. The cylinder radius matches the radius of section of region “3” outer boundary by the 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 is specified. The parameters of particle diffusion are chosen as follows. In the PWN we adopt the mean free path proportional to the particle gyroradius :

(27)

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

(28)

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

(29)

In Eqs. (27)–(29) and are the charge and mass of an electron or a positron, and is its energy. The values of and are free parameters of the scattering model, while 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 , introduced by Ostrowski (1991) is taken equal to (due to that, actually should not be less than ).

After injection, the particle is propagated through the system. The particle moves along a straight line during the time interval (where is its mean free path at the current location, and 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 , where is the mean free path in the region it has left, is the mean free path in the region it has entered and 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 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 equal cylinders along the axis. Each of these smaller cylinders is divided into coaxial cylindrical layers, whose outer radii grow linearly as , where . Each cylindrical layer is divided into equal segments by semi-planes , . 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 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 and values are splitted into and bins, respectively. The particle is detected by the bin with a given set of numbers (, , , , , ) if its energy corresponds to the -th momentum bin, direction of velocity to -th and -th angular bins, and it crosses the left boundary of the spatial bin (, , ) corresponding to . 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 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 . Here is the momentum in the local plasma rest frame at the beginning of the current piece of trajectory, is the new momentum value in the same frame, is the propagation time, and is the synchrotron loss coefficient:

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 , 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 centered at the point , the pulsar position, is introduced. The axis matches the symmetry axis and the plane matches the fixed axial section. We choose four groups of points in the plane. Namely, we take four rays starting at , lying in 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 34 = 12 points in regions “1”–“3”. We include five additional points in region “4”. Five rays parallel to the axis and lying in the plane () are chosen at coordinates . 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 , directed from the observer towards the pulsar (point ) and the plane of the sky (PoS), orthogonal to this axis and containing . The orientation of the source is determined by three angles , , and . is the angle between the axis and the line of intersection of the planes and PoS (chosen in the positive direction: counterclockwise when looking at the plane against the axis direction). is the angle between and axes. is the angle between and the North direction in the PoS. Thus, one may introduce a coordinate system centered at , and axes in the PoS, where the axis is aligned with the direction to the North. Then the angles , , and would be the three Euler angles defining the orientation of the system with respect to the system. In the following and in Fig. 2 we use .

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 43) (Gaensler et al., 2004; Hales et al., 2009). The spin-down power of J1747-2958 is about 2.5, while the X-ray luminosity of the nebula in the 0.5-8 keV range is about 5, and the observed PWN structure shows a narrow 45-long X-ray tail (the length is  1.1 pc at the 5 kpc distance). The best-fit power law photon index of the X-ray spectrum = 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 104 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  610.

Extended X-ray emission from the vicinity of J0437-4715 was detected with Chandra by Rangelov et al. (2016). This emission indicates the presence of a faint ( 310) PWN produced by J0437-4715. The X-ray photon index of the PWN is .

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 – 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 and high energy injection index 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 for the spin-down power 610.

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 ). 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 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 (as the size of the termination shock) and thus may reach (3–5) 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 and cm, from the simplified estimation of Eq. (23) one can obtain . Indeed, the modeled spectra of e pairs re-accelerated in the CSF demonstrate a break at 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 with = 2.3 upto = 210 (which is ) 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 bin, while the solid one refers to the bin located at the same distance from the pulsar as the bin, but lying on the ray directed from the pulsar to observer. The blue curves correspond to the bow shock region (bin – the solid curve, bin – the dashed curve). The green curves correspond to the region between the two shocks (bin – the solid curve, bin – the dashed curve). Worse statistics for bins and is related to the smaller squares of ’detectors’ near the system symmetry axis (see §4.3).

Figure 3: Non-smoothed Monte-Carlo simulated local spectra of e pairs in different zones of a BSPWN (see §4) for parameters similar to those of the source associated with PSR J0437-4715. The red curves correspond to the PW termination shock downstream: the dashed one refers to the bin, the solid one refers to the bin viewed along the ray directed to the observer. The blue curves correspond to the bow shock region (bin – the solid curve, bin – the dashed curve). The green curves correspond to the region between the two shocks (bin – the solid curve, bin – the dashed curve). The locations of all the bins are shown in Fig. 2. Worse statistics for bins and is due 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 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 , and (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 ) and the position angle . The curves show the spectral energy distribution () 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).

Figure 4: Simulated synchrotron images of an object similar to J0437-4715, for 8 eV ( 1550 Å, left) and 1 keV (right). The absolute intensity is given in normalized units to show the contrast change through the image.
Figure 5: Simulated synchrotron spectra of a source similar to J0437-4715. The orientation of the source is given by the angles , , and (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 ) and the position angle . The curves show the spectral energy distribution () calculated by integration of the spectral emissivity along the lines of sight, which cross the bins, indicated in the caption to Fig. 3. The colors and styles of the curves are the same as in Fig. 3.

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 (assuming the uncertain distance of about 400 pc). The H bow shock standing at 5 from the pulsar position was reported by Romani et al. (2010). The pulsar spin-down power is about 9.5. Deep Chandra observations of PSR J1741-2054 performed by Auchettl et al. (2015) have revealed a PWN with an extended () trail, whose power law photon index is = 1.670.06 and the total 0.5–10 keV luminosity is about 2.4 10. 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 where 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.

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 -ray pulsar PSR J0357+3205 (see, e.g., De Luca et al., 2013), whose space velocity is estimated to be about 390 (assuming the  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 250pc (Verbiest et al., 2012). Its age is about 340 kyr and its spin-down power is about 3.3. The estimated transverse velocity of Geminga is about 211, 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 44 eV and 195 eV. An equally good model consists of a blackbody component of 42 eV and a broken power law with photon indexes 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 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  0.2 pc long tails, an axial tail of  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 1, are much harder than that of the axial tail, where  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 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. For the bow shock moving at about 210 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 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 ( 1) are substantially harder than those from the other regions of the PWN ( 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 with 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 from the pulsar, to that in the tails would be then where 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 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 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 10. Its observed X-ray spectrum is well described by a photon index 1.58 (see Hui and Becker, 2007).

The X-ray observations of PSR B2224+65 revealed a jet-like feature extended to (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  1.00 in the extended X-ray filament, while a softer spectrum, with , 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 could produce the extended X-ray feature. The total magnetospheric potential of PSR B2224+65 can provide a maximum particle energy 7.5 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 3 . 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 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 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 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: .

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 . This maximum energy can only be achieved in a transrelativistic outflow with . 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 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 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 (1,000 ) pulsar of estimated spin-down power . The pulsar might originate from the supernova remnant MSH 11-61A at an estimated distance of 71 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 () 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 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, , a characteristic age of about 1.2 Myr, and the estimated distance is 600 pc. The bright tail has an angular size of 1220 corresponding to pc, while the faint one is of 250 size, i.e. pc. The X-ray luminosity of the bright tail in the 0.3-10 keV range is (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 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.

Figure 6: A 0.5-7 keV Chandra image of IGR J11014-6103 (the Lighthouse nebula) likely powered by PSR J1101-6101. The image shows the BSPWN as well as large-scale jet-like features, all launched by IGR J11014-6103 (see Pavan et al., 2016).

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

(30)

where

(31)

Here 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 exists in the beam leaving the accelerator if the net charge is nonzero in the CR beam.

The saturation level of the short scale magnetic field in the beam amplified by Bell’s instability can be estimated from the condition 1 (Bell, 2004). Then using Eq. (31) one can get

(32)

where is the fraction of pulsar spin-down power converted into the accelerated high energy particles forming the beam and is the radius of the beam which is .

The scale length of the beam , unless the length is limited by the particle energy losses which is the case for a small . The leading edge of the beam is located at a distance ahead of the region filled with the amplified magnetic field. The maximal growth rate is