Large-scale flow dynamics and radiation in pulsar \gamma-ray binaries

Large-scale flow dynamics and radiation in pulsar -ray binaries

V. Bosch-Ramon Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland;    M.V. Barkov Max Planck Institut für Kernphysik, Saupfercheckweg 1, Heidelberg 69117, Germany Space Research Institute, 84/32 Profsoyuznaya Street, Moscow, Russia
Received ¡date¿ / Accepted ¡date¿
Key Words.:
X-rays: binaries–ISM: jets and outflows–Radiation mechanisms: nonthermal
offprints: V. Bosch-Ramon;

Context:Several gamma-ray binaries show extended X-ray emission that may be associated to interactions of an outflow with the medium. Some of these systems are, or may be, high-mass binaries harboring young nonaccreting pulsars, in which the stellar and the pulsar winds collide, generating a powerful outflow that should terminate at some point in the ambient medium.

Aims:This work studies the evolution and termination, as well as the related radiation, of the shocked-wind flow generated in high-mass binaries hosting powerful pulsars.

Methods:A characterization, based on previous numerical work, is given for the stellar/pulsar wind interaction. Then, an analytical study of the further evolution of the shocked flow and its dynamical impact on the surrounding medium is carried out. Finally, the expected nonthermal emission from the flow termination shock, likely the dominant emitting region, is calculated.

Results:The shocked wind structure, initially strongly asymmetric, becomes a quasi-spherical, supersonically expanding bubble, with its energy coming from the pulsar and mass from the stellar wind. This bubble eventually interacts with the environment on  pc scales, producing a reverse and, sometimes, a forward shock. Nonthermal leptonic radiation can be efficient in the reverse shock. Radio emission is expected to be faint, whereas X-rays can easily reach detectable fluxes. Under very low magnetic fields and large nonthermal luminosities, gamma rays may also be significant.

Conclusions:The complexity of the stellar/pulsar wind interaction is likely to be smoothed out outside the binary system, where the wind-mixed flow accelerates and eventually terminates in a strong reverse shock. This shock may be behind the extended X-rays observed in some binary systems. For very powerful pulsars, part of the unshocked pulsar wind may directly interact with the large-scale environment.

1 Introduction

Gamma-ray binaries are compact sources that present nonthermal emission in the GeV and/or TeV bands and, typically, they are also moderate-to-strong emitters in radio and X-rays. The radiation in gamma-ray binaries is related to transient or persistent outflows, which can interact either with themselves in the form of internal shocks or with the environment on different scales. The number of members of the class of gamma-ray binaries is growing, being at present around ten, including unconfirmed candidates (e.g. Albert et al. alb06 (); alb07 (); Aharonian et al. aha05a (); aha05b (); Abdo et al. abd09a (); abd09b (); abd09c (); abd10a (); abd10b (); abd11 ();Tavani et al. tav09a (); tav09b (); Sabatini et al. sab10 (); Williams et al. wil10 (); Tam et al. tam10 (); Hill et al. hil11 (); Falcone et al. bon11 (); Corbet et al. cor11 ()).

Several gamma-ray binaries present hints or clear evidence of interactions between an outflow and their medium on large scales: Cygnus X-3, LS I +61 303, LS 5039, PSR B125963, and the gamma-ray emitter candidate Cygnus X-1 (Heindl et al. hei03 (); Sánchez-Sutil et al. san08 (); Paredes et al. par07 (); Pavlov et al. pav11a (); Durant et al. dur11 (); Martí et al. mar96 (); Gallo et al. gal05 (); Russell et al. rus07 ()).

Cygnus X-1 and Cygnus X-3 are high-mass microquasars, in which accretion leads to jets that can interact with the ISM (e.g. Velázquez & Raga vel00 (); Heinz hei102 (); Heinz & Sunyaev hei202 (); Bosch-Ramon et al. bos05 (); Zavala et al. zav08 (); Bordas et al. bor09 (); Bosch-Ramon et al. bos11 ()). On the other hand, PSR B125963 is a system formed by an O9.5 V star and a nonaccreting millisecond pulsar (Johnston et al. joh92 (); Negueruela et al. neg11 ()). Instead of coming from a jet, like in microquasars, the nonthermal emission of this source is thought to originate in the region where the star and the pulsar winds collide (Tavani & Arons tav97 ()), with its radio emission extending far away from the binary (Moldón et al. mol11a ()). PSR B125963 also presents extended X-ray emission on scales of (i.e. a projected linear size of  cm at the source distance of 2.3 kpc; Negueruela et al. neg11 ()), with a flux  erg s cm and photon index (Pavlov et al. pav11a ()). The extended radio and X-ray emission would be produced in different regions but still in the shocked stellar-pulsar wind structure, which propagates away from the binary and should eventually terminate in the external medium.

The nature of LS 5039 and LS I +61 303 is still unknown due to the lack of clear pulsar or accretion features, conclusive system dynamical information, or nonthermal emission evidence (e.g. Casares et al. cas05a (), Casares et al. cas05b (), Sarty et al. sar11 (); Rea et al. rea10 (); rea11 (), McSwain et al. mcs11 (), Ribó rib11 (); see also Bosch-Ramon & Khangulyan bos09 ()). In the past, several works have considered LS 5039 and LS I +61 303 good microquasar candidates because of their radio morphology and emission properties (e.g. Taylor et al. tay92 (); Paredes et al. par00 (); Massi et al. mas01 (); mas04 (); Paredes et al. par02 (); Romero et al. rom05 (); Paredes et al. par06 (); Bosch-Ramon et al. bos06 (); Massi & Kaufman Bernadó mas09 ()). However, it has also been proposed that the phenomenology of these sources and, in particular, their milliarcsecond (mas) scale radio morphology may be more compatible with a nonaccreting pulsar scenario (e.g. Maraschi & Treves mar81 (), Martocchia et al. mar05 (), Dubus dub06 (), Chernyakova et al. che06 (), Sierpowska-Bartosik & Torres sie08 (), Sierpowska-Bartosik & Torres sie09 (); see however Romero et al. rom07 ()). Both binaries are high-mass systems, with LS 5039 harboring an O6.5 V star (Clark et al. cla01 ()), and LS I +61 303, an early Be star (Hutchings & Crampton hut81 ()).

Extended X-rays have been found in LS 5039, with angular size ( cm at 2.5 kpc; Casares et al. cas05a ()), X-ray flux  erg s cm and photon index (Durant et al. dur11 ()), and possibly also in LS I +61 303, with ( cm at 2 kpc; Frail & Hjellming fra91 ()), and X-ray flux  erg s cm (Paredes et al. par07 ()111Rea et al. rea10 () did not find evidence of extended X-rays in LS I +61 303 in longer observations than those studied in Paredes et al. par07 (). However, in their observations the X-ray counts were integrated along one spatial dimension in the CCDs, which only allows finding an extension signal in specific directions.). It is noteworthy that the gamma-ray binaries HESS J0632057 (Aharonian et al. aha07 (); Hinton et al. hin09 (); Skilton et al. ski09 (); Falcone et al. fal10 (); bon11 (); Moldón et al. mol11b ()) and 1FGL J1018.65856 (Corbet et al. cor11 (); Pavlov et al. pav11b (); de Ona Wilhelmi et al. deo10 ()) may host a nonaccreting pulsar as well, although a microquasar scenario cannot be discarded.

Summarizing, among the several gamma-ray binaries presenting hints or evidence of interactions with the medium, three of them may (one for sure) host a nonaccreting pulsar. In addition, extended X-ray emission may be eventually discovered as well in HESS J0632057 and 1FGL J1018.65856, another two pulsar binary candidates. At present, the dynamics and radiation of an outflow from a pulsar gamma-ray binary has not been studied far from the binary system. Given the recent findings of extended emission, it seems necessary to analyze the flow evolution beyond the stellar/pulsar wind region and its interaction with the ISM. This is the goal of this work, which is distributed as follows. In Sect. 2, the main properties of the binary generated flow are described. In Sect. 3, the general properties of the interaction between the flow and the external medium, the progenitor supernova remnant (SNR; young sources), or the interstellar medium (ISM; older sources) are characterized analytically, and in Sect. 4, the expected nonthermal emission from radio to gamma rays is computed. Finally, in Sect. 5, the results are discussed in the context of PSR B125963, LS I +61 303 and LS 5039 (assuming that the last two host a powerful pulsar), and some final remarks are made in Sect. 6. All through the paper, cgs units are used.

2 Physical scenario

2.1 The stellar-pulsar wind shock

In Figure 1, the pulsar high-mass binary scenario is sketched. For a general description of the scenario on the binary scales, see for instance Tavani & Arons (tav97 ()). As seen in the figure, the stellar and the pulsar winds collide to form two bow-shaped standing shocks. An important parameter that characterizes the wind interaction region is the pulsar to stellar wind momentum flux ratio (e.g. Bogovalov et al. bog08 ()):


where is the stellar mass-loss rate, is the stellar wind velocity, and is the pulsar spin-down luminosity, taken here as equal to the kinetic luminosity of the pulsar wind. The stellar wind is assumed isotropic. The pulsar orbital velocity, of a few hundred km s, is much lower than , typically  cm s for massive stars (e.g. Puls et al. pul09 ()), so it has been neglected in this section. However, the orbital velocity may play a relevant role in determining in some cases, e.g. in the presence of a stellar decretion disk, as discussed in Sect. 2.2.3. It is also noteworthy that the pulsar wind is simplified here as isotropic, although a more refined treatment should account for anisotropy (see, e.g., Bogovalov & Khangoulian bog02 (), Bogovalov et al. bog11 ()).

For reasonable -values (a scenario with is briefly discussed in Sect. 2.2.2), the shocked wind structure is dominated by the stellar wind ram pressure, and points away from the star on the binary scales. The size of the wind colliding region can be characterized by the distance between the two-shocks contact discontinuity (CD) and the pulsar,


which is  cm for and a distance between the star and the pulsar of  cm. The distance between the CD and the star is then  cm. The CD, which separates the shocked stellar and pulsar winds, has an asymptotic half-opening angle (see Bogovalov et al. bog08 () and references therein):


Whereas the main contribution of momentum an mass generally comes from the stellar wind, the energy is mostly provided by the pulsar wind, as seen from the pulsar to stellar wind luminosity ratio:


We notice that the wind of the pulsar prevents the latter to accrete from the stellar wind for a wide range of -values. As a conservative estimate, we can derive the -value for which is equal to the gravitational capture radius (Bondi & Hoyle bon44 ()):


where is the pulsar mass, and .

Figure 1: a) Sketch of the pulsar high-mass binary scenario on scales larger than the binary system. The star (S) and pulsar (P) shocked winds trace a spiral-shape due to the Coriolis force and the pulsar orbital motion. Eventually, both shocked winds mix due to Kelvin-Helmholtz instabilities, and the spiral arms merge. b) Cartoon of the bubble (case ) driven by the shocked flows formed within the binary system (BS) and made of pulsar- and stellar-wind material. The bubble expands and accelerates, eventually terminating in a shock where its ram pressure is equal to the pressure of the SNR or the shocked ISM for old enough sources. This occurs at the contact discontinuity (CD).

2.2 Shocked flow evolution

After being shocked, the pulsar wind facing the star moves at a speed , but suffers a strong reacceleration along the CD produced by a strong pressure gradient, related to the anisotropy and inhomogeneity of the stellar wind at the pulsar location (Bogovalov et al. bog08 ()). Because of the very low density in the downstream region, the pulsar wind shock is adiabatic (Bogovalov et al. bog08 (); Khangulyan et al. kha08 ()). The stellar wind shock could be either adiabatic or radiative, but in either case the stellar wind velocity will be . Therefore, the two shocked winds have a relative velocity , depending on how much the shocked pulsar wind has reaccelerated. Such a large velocity difference is likely to lead to the development of Kelvin-Helmholtz (KH) instabilities in the CD. This will produce entrainement of shocked stellar wind into the lighter and faster shocked pulsar wind. If the stellar wind shock is radiative, then the matter downstream collapses into a dense and cool thin layer, diminishing even farther the stability of the whole shocked structure (e.g. Stevens et al. ste92 (); Pittard pita09 ()). Perturbations generated in the interface will initially move with the stellar shocked wind, growing on a timescale , where is the size of the perturbation, and the sound speed of the shocked stellar wind. Taking  cm s and  cm, one obtains  s, approximately the lapse needed by the perturbations to enter in the nonlinear regime. This implies that instabilities, turbulence, and mixing will develop on spatial scales that are similar though slightly larger than those of the binary system.

We now discuss the impact of the orbital motion on the shocked winds. To do this, we focus on three cases: the possibly more common , (e.g.  yr for  cm s and  erg s); more briefly, the extreme case ; and under the presence of a stellar decretion disk.

2.2.1 Stellar wind momentum flux dominance

We analyze the dynamics of the gas flow in the coordinate system that corotates with the binary system at an angular velocity , where is the orbital period. The X-axis coincides with the line connecting the centers of the stars, the Y-axis lies in the orbital plane, and the Z-axis is normal to the orbital plane. In such a coordinate system, the Coriolis acceleration is . The stellar wind moves along the X-axis with speed , acquiring a velocity in the Y-direction. In the frame of the shocked-wind flow, , and is the distance, along the X-axis, between the pulsar and a point farther away from the star. The dynamical pressure of the stellar wind in the Y-direction can be estimated as (). From this, the location where the pulsar wind total pressure is balanced by the stellar wind can be derived with


where , with the stellar wind density at the pulsar location. For , Eq. (6) can be simplified to


and for ,


For typical parameters, the solution is closer to the one given by Eq. (8), therefore the bending distance can be written as


where . This distance is much less than in the purely ballistic case, for which  cm.

It is noteworthy that the pulsar wind is shocked even in the direction opposite to the star, within a distance . Therefore, a significant fraction of the pulsar wind luminosity is reprocessed outside the binary system. Interestingly, this region may dominate the very high-energy output in pulsar gamma-ray binaries with high photon-photon absorption (see the discussion in Bosch-Ramon et al. bos08 ()).

If KH instabilities by themselves have not already mixed the shocked winds, isotropizing their particle, momentum, and energy fluxes, together with the Coriolis force they very likely will. Eventually, the shocked-wind flow will probably end up as a trans- or subsonic turbulent bubble, loaded with stellar wind mass and pulsar wind energy. The bubble is not confined because there is still a strong pressure gradient radially outwards. Therefore, and assuming that the whole process is adiabatic, any thermal, turbulent, or magnetic energy will eventually end up in the form of radial bulk motion. The role of the magnetic field should not change this picture much, unless the initial stellar and pulsar winds had a dynamically dominant magnetic field. (For the role of the magnetic field on binary scales; see Bogovalov et al. bog11 ().) The maximum velocity expansion of the bubble will be roughly its initial sound speed, which for a maximum entropy initial state (i.e. right after isotropization) is


Typically, the mas radio emission in pulsar massive binaries will come from scales larger than , and thus this radiation should be strongly affected by the Coriolis force, as well as by shocked-wind mixing. The situation described in this section is pictured in Fig. 1 (left).

2.2.2 Pulsar wind momentum flux dominance

For a dominant pulsar wind, the shocked structure closes towards the star. In this case, the Coriolis force exerted by the light pulsar wind is too low and the shocked structure moves ballistically, with a bending typical distance . Along the orbit, however, the natural bending of the shocked structure directly exposes it to the pulsar wind ram pressure, which pushes and opens the spiral to some degree. As for , for the instabilities in the CD will also take place, and the formation of a mixing layer is expected to occupy part of the region of the free pulsar wind. However, in the direction perpendicular to the orbital plane, the pulsar wind will likely propagate freely. The energetics of this free wind may be significant, eventually terminating as a jet-like structure interacting with the surrounding medium (see, e.g., Bordas et al. bor09 ()). A sketch of the two cases discussed, with and , is shown in Fig. 2.

Figure 2: Sketch of the two possibilities for the outcome of the interaction between the pulsar and the star winds, in the context of a binary system. To the right, the case for a very powerful pulsar, with , is presented. To the left, the case with .

2.2.3 The equatorial flow case

The situation discussed so far is that of two spherical winds interacting with each other. However, some of the systems considered in this work host Be stars, which present decretion disks around the equator. These disks are thought to be much denser than the radiatively driven stellar wind. The disk flow is quasi-Keplerian, and moves subsonically in the radial direction with velocities around  km s (e.g. Porter & Rivinius por03 ()). Given the slow flow velocity, the flow speed determining the flux momentum ratio, in the pulsar reference frame, is near the orbital velocity of the pulsar. Outside the equatorial disk, the radiatively driven wind in Be stars is most likely polar. Such a configuration of the two stellar flows probably makes the geometry of the colliding wind region quite complex, enhancing the development of instabilities in the shocked-flow contact discontinuity (for simulations of this scenario on binary system scales, see Romero et al. rom07 () and Okazaki et al. oka11 ()). The shocked disk material has a much larger density than the light polar wind. This implies that development of instabilities may be slower, although the sound speed of the shocked disk will be similar to the pulsar’s orbital velocity, the dynamical speed in this case. The high density of the shocked equatorial flow implies that it will likely be radiative, thereby potentiating instabilities and fragmentation. Finally, the mass-loss rate of the disk is similar or even smaller than that of the polar wind, and thus once accelerated by the pulsar wind, the flow evolution on the largest scales should not be very different from the isotropic wind case. The dense and fragmented shocked disk material will likely introduce inhomogeneities into the larger scale shocked flow, although eventually the isotropization of the particle, momentum, and energy fluxes of the larger scale structure seems a reasonable guess. However, hydrodynamical simulations on larger scales than the binary are mandatory for validating the assumption.

2.3 Bubble large-scale evolution

The pulsar started its life after a supernova explosion. For an explosion energy  erg and ISM density  g cm ( cm), the SNR becomes radiative (Blinnikov et al. biu82 ()) after  yr when achieving the radius  pc, where and . In this context, the binary system will leave the SNR after a time


i.e.,  yr for  cm s.

Inside the SNR, given the high sound speed of the ambient medium, the binary driven bubble cannot generate a forward shock. Once in the ISM, the bubble forms a slow forward shock. In both situations, the supersonic expanding bubble suffers a strong reverse shock to balance the external pressure, either from the hot SNR ejecta or the shocked ISM. This reverse shock is expected to have a speed in the bubble reference frame. If the pulsar is powerful enough, the reverse shock can power a moderately bright nonthermal emitter, with a potential nonthermal luminosity as high as . This estimate is done by neglecting the radiation losses on the binary system scales, but as long as adiabatic losses are likely to dominate in the colliding wind region (e.g. Khangulyan et al. kha08 ()), this approximation should suffice at this stage. Thermal emission from the reverse shock can be discarded, given the low densities. Some extended thermal emission will be generated in the shocked ISM, but with its peak in the ultraviolet it may be hard to detect. Figure 1 (right) sketches the bubble termination region.

3 Dynamics of the bubble interacting with the medium

The interaction between the bubble and its environment takes place in three steps. At early times, the bubble is still surrounded by the SNR in its adiabatic or Sedov phase (Sedov sed59 ()). Later on, the SNR forward shock in the ISM becomes radiative, entering into the so-called snow-plow phase (Cox cox72 (); Blinnikov et al. biu82 ()). Finally, after the binary leaves the SNR or the latter dissipates, the bubble interacts directly with the ISM.

In the adiabatic regime, the SNR properties can be characterized through the SNR/ISM forward shock radius (Sedov sed59 ())


where is the pulsar age, velocity


and the inner pressure is


assumed here to be constant. Equilibrium between bubble and SNR pressures determines the location of the boundary between these two regions: , and the equilibrium between bubble preshock and postshock total pressures determines the location of the bubble reverse shock:


Since the hot ejecta region has a high sound velocity, the bubble forward shock inside the SNR quickly becomes a sound wave, and the bubble shape becomes spherical.

For , the SNR enters into the radiative phase, so the snow-plow solution has to be used (Blinnikov et al. biu82 ()). The SNR/ISM forward shock has now a radius of


a velocity of


and an inner pressure of


where . As before, the bubble reverse shock radius can be estimated using Eq. (15).

The bubble/ISM direct interaction has a strong resemblance to the interaction between a supersonic stellar wind and the ISM. Therefore, for this case the solution for a supersonic wind with continuous injection has been adopted (Castor et al. cas75 ()). This renders the ISM forward shock location at


moving with velocity


The inner pressure is now


and again the preshock and postshock total pressure balance gives the location of the bubble reverse shock: . Since the bubble external medium has a low sound speed, the precise shape of the bubble is less clear now, although it has been assumed to be spherical, as discussed in Sect. 2.2.

For old enough sources, say  yr, the proper motion of the system can affect the ISM shock, leading to the formation of a bow-shaped ISM shock. In this case, the reverse shock is located approximately at the point where the bubble and ISM ram pressures become equal:


Using the model just described, the relevant distances, velocities, and pressures can be computed for different system ages. Parameter values similar to those of LS 5039 and LS I +61 303 have been adopted, and their values are presented in Table 1. The dynamical, as well as the radiative (see next section), results may also apply to PSR B125963, although the comparison is less straightforward, as discussed in Sect. 5. It is noteworthy that, as long as is much smaller than the largest scales of interaction with the medium, the size difference between PSR B125963 and the other two sources should not introduce strong differences in the bubble propagation and termination. The results are shown in Figs. 3 and 4. The breaks in reverse shock and CD distances, and pressures, apparent in the figures for the SNR case at  yr, are caused by the significant loss of internal energy (a fraction 1-) that takes place in the adiabatic-to-radiative phase transition (Blinnikov et al. biu82 ()).

The approach carried out here should be reasonable at a semi-quantitative level, but there are some caveats. First, these different stages of the bubble evolution are idealizations, and mixed situations could take place. Also, the initial conditions of the bubble were oversimplified, assuming that most of the complexity of the interaction on the binary scales is smoothed out. In particular, the importance of the geometry of the interacting outflows (e.g. polar winds, disk, etc.) is to be studied in detail. On larger scales, the shocked bubble, SNR and ISM media were approximated as homogeneous static gas. Rayleigh-Taylor instabilities in the CD between the denser shocked ISM shell, and its inner region, were neglected, along with anisotropies in the ISM. To account for all this properly, magnetohydrodynamical simulations of the bubble formation, evolution and termination are planned. We also neglected the pulsar spin-down luminosity evolution with time, which may imply an underestimate of the total bubble injected energy by a factor of a few.

Parameter [units] value
Star mass-loss rate [ yr]
Stellar wind velocity [cm s]
Pulsar spin-down luminosity [erg s]
SNR energy [erg]
ISM number density [cm] 10
Magnetic-to-total energy density ratio 0.01, 0.1
Nonthermal-to-total luminosity ratio 0.01
Star luminosity [erg s]
Stellar photon energy [eV] 10
Reverse shock velocity [cm s]
Source distance [kpc] 3
Table 1: Adopted parameter values
Figure 3: a) Evolution with time of the radii of the SNR/ISM forward shock (black solid line), SNR/bubble contact discontinuity (red dashed line), and bubble reverse shock (blue dotted line). The CD radius becomes larger than the SNR shock radius, which indicates that the model adopted here does not apply after a few  yr. b) Evolution with time of the radii of the bubble/ISM forward (black solid line) and reverse shocks (blue dotted line).
Figure 4: a) Evolution with time of the SNR/ (black solid line) and bubble/ISM forward shock velocities (red dashed line). b) Evolution with time of the SNR/ (black solid line) and bubble/ISM pressures (red dashed line). c) Evolution with time of the SNR/ (black solid/blue dashed line) and bubble/ISM magnetic fields (red dotted/green dot-dashed line).

4 Nonthermal radiation estimates

The reverse shock is a good candidate for nonthermal emission, since it is fast and strong, and basically all the bubble energy goes through it. The forward shock in the ambient medium, on the other hand, is either much slower and also less energetic (bubble/ISM case), or just a sound wave (SNR case).

4.1 Emitter characterization

To compute the particle acceleration rate in the reverse shock, we adopted the expression for a non-relativistic hydrodynamical shock in the test particle and Bohm diffusion approximations: (e.g. Drury dru83 ()), where is the magnetic field and the particle charge. The same and diffusion coefficient values were assumed in both sides of the shock. The resulting energy distribution of the particles injected in the emitter depends on energy as . The -energy density was taken equal to , with and 0.01. Smaller values seem less plausible, because some magnetic field is expected to be transported from the binary scales, but cannot be discarded. The maximum energies are the minimum value between the synchrotron-limited case, , with , and the diffusive-escape one, . Although the stellar photon energy density, the dominant radiation field, is higher than the -energy density at , the Klein-Nishina (KN) effect makes IC losses less important than synchrotron (or escape) ones at . For a wide source age range, in the case of electrons is  TeV, but determined by synchrotron cooling for , and by diffusive escape for . For protons, the maximum energies are  TeV for and  TeV for , both limited by diffusive escape. The differences between the maximum energies in the SNR and the bubble/ISM case are small; in particular, under diffusive escape, they are equal. Synchrotron-limited maximum energies are , but diffusive-escape ones are constant, because is also constant.

The conditions downstream of the bubble reverse shock, in particular the low densities, render hadronic processes inefficient. Therefore, we have focused on electrons emitting via synchrotron and inverse Compton (IC) scattering. The most energetic protons may diffuse away from the reverse shock to reach the shocked ISM, but to be efficient, this mechanism requires target densities to be much higher than those considered here. To compute the synchrotron and IC emission (e.g. Blumenthal & Gould blu70 ()), we modeled the emitter as only one zone, taking homogeneous - and radiation fields. This approach is valid for the synchrotron calculations as long as the flow is subsonic, and is the same everywhere, which is admittedly a rough approximation in an ideal MHD flow. The role of diffusion in electron transport is negligible in the present context unless diffusion coefficients are much higher than Bohm. Owing to fast cooling, synchrotron X-rays are radiated close to the reverse shock region, whereas the synchrotron radio emission comes from a more distant location due to advection. The particle flux conservation means that the advection speed downstream of the reverse shock is roughly , where is the distance to the reverse shock. This makes particles accumulate at


with the relevant cooling time. Most of the IC radiation comes also from a distance to the star. This is not true for very high-energy IC photons, but for the adopted -values and due to the KN effect, their fluxes are minor so were neglected here (see however Sect. 5). Given the quasi-spherical shape of the emitter, the target photon field for IC has been taken as isotropic and monoenergetic, given its narrow band. The stellar luminosity and photon energy have been fixed to  erg s and 10 eV. The luminosity injected in the form of nonthermal electrons at the bubble reverse shock was taken as 1% of the pulsar wind luminosity: . Although is hard to determine, changes in this parameter only affect the nonthermal fluxes linearly. Adiabatic cooling has been included in the calculations, with the cooling timescale taken as .

4.2 Spectral energy distributions and lightcurves

After characterizing the emitter as a region with roughly homogeneous properties, the maximum particle energies, and the dominant cooling processes (adiabatic and radiative), it is easy to calculate the electron evolution. For that, a constant particle injection is assumed, and electrons are evolved up to the source age of interest. Once the particle energy distribution is known, the synchrotron and IC spectral energy distributions, as well as the specific and integrated fluxes, are calculated by convolving the synchrotron and IC specific power per electron by the particle energy distribution.

The results of the radiation calculations are shown in Figs. 5 and 6. Figure 5 shows the computed synchrotron and IC spectral energy distributions for a source age of  yr (SNR) and  yr (bubble/ISM), with the remaining parameter values given in Table 1. Figure 6 presents the time evolution of the surface brightness at 5 GHz, and the 1-10 keV and 0.1-10 GeV bolometric fluxes. The radio surface brightness was computed by just dividing the specific flux by the source solid angle in units at 3 kpc, taking as the source size. We point out that our aim is just to provide with a crude estimate of the expected radiation. More detailed calculations of the nonthermal emission will be presented elsewhere.

As seen in the figures, the break of the adiabatic-to-radiative phase transition appears in all the SNR lightcurves, although is more apparent for the radio surface brightness and IC fluxes, which are strongly affected by the jump in . The general evolution of produces the long-term strong decay of the SNR radio and GeV lightcurves, and a smoother decay of the bubble/ISM radio lightcurve, whereas the bubble/ISM GeV lightcurve remains more or less constant. The differences are partially caused by the evolution of the relative importance of the different radiation channels and the adiabatic one, which is not the same in the SNR and the bubble/ISM cases. The synchrotron and IC emissivities have also different time dependencies, the former depending on as well. For , X-ray emitting electrons radiate all their energy through synchrotron, rendering a flat X-ray lightcurve and photon indexes . For , with the electron maximum energies limited by diffusive-escape, the maximum energy of synchrotron photons () decreases with time, and it quenches the X-ray flux for  yr. This also explains the steepening of the synchrotron X-ray spectrum for the bubble/ISM case seen in Fig. 5. The SNR X-ray spectrum does not show this feature because it has been computed for .

Figure 5: Spectral energy distributions of the synchrotron (thick lines) and IC (thin lines) emission computed for the SNR ( yr) and the bubble/ISM ( yr) cases, and .
Figure 6: a) Lightcurve of the radio (5 GHz) surface brightness. b) Lightcurve of X-ray flux in the range 1–10 keV. c) Lightcurve of the IC flux in the range 0.1–10 GeV. Both the SNR and the bubble/ISM cases are shown in the age range  yr and for .

5 Discussion

The sizes of the bubble reverse shock and its downstream region give an idea of the source angular size at different wavelengths. At 3 kpc, the hard X-ray emitter would have when still inside the SNR ( yr), and in the bubble/ISM case ( yr). These values were calculated for the , , and values given in Table 1, although we need to recall their weak inter-dependence: . Extended X-ray fluxes could be detectable by an instrument like Chandra up to  yr. It is worth noting that, if synchrotron X-ray energies are detected from the reverse shock, the diffusion coefficient should be close to Bohm.

Radio emission is possibly detectable for the SNR case with , it may be hard to be found for the case , and its detection seems discarded for the bubble/ISM case for any -value.

The GeV fluxes seem too low to be detectable by any current instrument, and the chances are even lower in TeV due to the Klein-Nishina effect and dominant synchrotron cooling. However, it may still be possible to get higher IC GeV and TeV fluxes increasing ten times , and reducing strongly not to overcome the radio/X-ray constraints. This would lead to GeV-TeV fluxes that could be detectable by present (GeV: Fermi, AGILE; TeV: HESS, MAGIC, VERITAS) and/or forthcoming instruments (e.g. TeV: CTA). Given the spectral break at  GeV introduced by adiabatic losses, the bubble reverse shock would be a good target for a  GeV-threshold Cherenkov instrument.

We now briefly discuss whether the evidence or hints of extended X-rays from LS 5039, LS I +61 303, and PSR B125963 can be understood in the scenario posed here.

5.1 Ls 5039

LS 5039 is likely a young source with an age range  yr, and it has already abandoned its SNR (e.g. Ribó et al. rib02 (); Moldón et al. mol11c ()). The angular size at 2.5 kpc for the bubble/ISM case is consistent with the observed angular size of the emitter, , and the fluxes can be explained in the framework given here. The photon index at distances is , consistent with the SED shown in Fig. 5, and then softens farther out reaching at (Durant et al. dur11 ()). This can be explained by the cooling of the highest energy electrons close to the reverse shock, thereby depleting the highest energy part of photons and softening the spectrum farther from the shock. The adopted  cm is similar to the values derived for the surroundings of the LS 5039 SNR candidate (Ribó et al. rib02 ()), and the pulsar is probably  erg s, which is required to explain the GeV luminosity of the source (see Sect. 4 in Zabalza et al zab11a ()). The morphology of the extended radio emission found by Durant et al. (dur11 ()) is not completely spherical, which may be explained by the only partial isotropization of the momentum flux in the inner regions of the bubble (see Sect. 2.2). A scenario with could be also behind the anisotropy (see Sect. 2.2.2), but it seems less likely.

5.2 Ls I +61 303

Assuming that the hints of extended X-rays in LS I +61 303 are real, we see that the angular size, , is compatible with the bubble/SNR scenario for  yr, but with a low nonthermal efficiency or a very low magnetic field, given the weak X-ray flux (Paredes et al. par07 ()) and the radio nondetection on those scales (Martí et al. mar98 ()).

The SNR progenitor of LS I +61 303 has not yet been found. However, extended 5 GHz emission has been detected centered on the location of LS I +61 303, with a typical size of 6–8 arcminutes and fluxes of tens of mJy (Martí et al. mar98 (); Muñoz-Arjonilla et al. mun09 ()). It has been argued that the lack of extended X-rays on the same scales, and the low surface radio brightness (in a nonthermal scenario), may go against an SNR interpretation (Martí et al. mar98 (); Muñoz-Arjonilla et al. mun09 ()). However, the radio emission surrounding LS I +61 303 may be mostly thermal (Muñoz-Arjonilla et al. mun09 ()), coming from an SNR/ISM shock close to its radiative phase. If the values in the region of LS I +61 303 were a bit higher than those adopted here, an SNR remnant with energy  erg would enter its radiative phase after  yr, with a radius  cm, a few arcminutes at 2 kpc. Thermal X-rays would not be expected in that case, since the shocked ISM material would be too cold. The slow proper motion of the source (Dhawan et al. dha06 ()) is compatible with LS I +61 303 being at the center of the SNR (see however Martí et al. mar98 () for image deformation). As in LS 5039, the -value in LS I +61 303 is expected to be  erg s because of the high GeV luminosity of the source (see Sect. 5 in Zabalza et al zab11b ()).

5.3 Psr B125963

PSR B125963 has an age of  yr (Wex et al. wex98 ()), and is likely now interacting directly with the ISM. With this age, and the pulsar spin-down luminosity  erg s (Manchester et al. man95 ()), Eq. (19) yields  cm ( 1’), assuming  cm. This is about one arcminute at 2.3 kpc and much larger than the observed (Pavlov et al. pav11a ()). Nevertheless, in PSR B125963 it seems likely that the proper motion has already bow-shaped the ISM shock. In that case, for  cm and  cm s, the reverse shock would be located at a distance of  cm from the system, or neglecting projection effects. Unfortunately, the large errors of the proper motion measurements (Zacharias et al. zac09 ()) and the low statistics of the X-ray data make it difficult to compare the X-ray extension and proper motion directions.

There is another possibility behind the extended X-rays found in PSR B125963. The most significant detection comes from , or  cm in linear (projected) size, which is . Although the flow bending distance, , should be less than that, the X-ray emission may be related to the asymmetric shock produced by the Coriolis force, and/or to spiral-arm merging. This radiation could not be resolved in the case of LS I +61 303 and LS 5039, which are much more compact than PSR B125963.

6 Final remarks

Gamma-ray binaries hosting powerful pulsars can produce supersonic flows that originate in the interaction of a pulsar wind with the wind of the companion, which may present strong anisotropies and inhomogeneities. The likely mixing of these flows on larger scales than the binary will render an expanding, rather isotropic, supersonic bubble, which will eventually terminate in the surrounding medium. The interaction of the supersonic bubble with the surrounding medium can take place in different ways depending on the age of the pulsar. For young sources, the environment will be the hot SNR ejecta, whereas it will be the shocked ISM for older objects. Although this interaction is expected to be rather symmetric, for old enough sources, possibly like PSR B125963, the proper motion can bow-shape the interaction structure with the ISM.

The radiation from the shocked bubble interacting with the environment may explain the observed extended X-ray emission from LS 5039, and possibly also from LS I +61 303 (if real). For PSR B125963, the found extended X-rays could come from inner regions of the bubble, triggered perhaps by Coriolis force shocks or some other type of dissipation mechanism. Extended X-ray emission may also eventually be detected in HESS J0632057 and 1FGL J1018.65856, which could also originate in shocked wind outflows. The shape of the interaction region might allow the distinction of the driving flow, i.e. a jet or a pulsar wind, but as shown in Sect. 2.2.1 it may be not straightforward.

The complexity of the flow evolution on the different relevant scales makes a proper analytical characterization difficult, and thus makes magnetohydrodynamical simulations important. The geometry, level of anisotropy, velocity and density of the stellar flow requires careful study, in particular for binaries hosting stars with an equatorial disk, like PSR B125963 and LS I +61 303. Anisotropy in the pulsar wind has been also neglected here, but may play some role, at least up to intermediate scales.

The authors want to thank Dmitry Khangulyan for thoroughly reading the manuscript, and for his useful comments and suggestions. The authors are grateful to the referee, Ignacio Negueruela, for many constructive and useful comments. The research leading to these results has received funding from the European Union Seventh Framework Program (FP7/2007-2013) under grant agreement PIEF-GA-2009-252463. V.B.-R. acknowledges support by the Spanish Ministerio de Ciencia e Innovación (MICINN) under grants AYA2010-21782-C03-01 and FPA2010-22056-C06-02. B.M.V. thanks to State contract 2011-1.4-508-008/9 from FTP of RF Ministry of Education and Science.


  • (2009a) Abdo, A. A., et al. 2009a, 2009, Science, 326, 1512
  • (2009b) Abdo, A. A., et al. 2009b, ApJ, 701, L123
  • (2009c) Abdo, A. A., et al. 2009c, ApJ, 706, L56
  • (2010a) Abdo, A. A., et al. 2010a, Science, 329, 817
  • (2010b) Abdo, A. A., et al. 2010b, ApJ, 723, 649
  • (2011) Abdo, A. A., et al. 2011, ApJ, 736, L11
  • (2005a) Aharonian, F. A. et al. 2005a, A&A, 442, 1
  • (2005b) Aharonian, F. A. et al. 2005b, Science, 309, 746
  • (2007) Aharonian, F. A. et al. 2007, A&A, 469, L1
  • (2006) Albert, J. et al. 2006, Science, 312, 1771
  • (2007) Albert, J. et al. 2007, ApJ, 665, L51
  • (1982) Blinnikov, S. I., Imshennik, V. S., Utrobin, V. P. 1982, SvAL, 8, 361
  • (1970) Blumenthal, G. R., Gould, R. J., 1970, Rev. Mod. Phys., 42, 237
  • (2002) Bogovalov, S. V., & Khangoulian, D. V. 2002, MNRAS, 336, L53
  • (2008) Bogovalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., Aharonian, F. A. 2008, MNRAS, 387, 63
  • (2011) Bogovalov, S., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., Aharonian, F.A. 2011, MNRAS, submited [astro-ph/1107.4831]
  • (1944) Bondi, H. & Hoyle, F.,1944, MNRAS, 104, 273
  • (2011) Bongiorno, S. D., Falcone, A. D., Stroh, M., et al. 2011, ApJ, 737, L11
  • (2009) Bordas, P., Bosch-Ramon, V., Paredes, J. M., Perucho, M. 2009, A&A, 497, 325
  • (2009) Bosch-Ramon, V. & Khangulyan, D. 2009, IJMPD, 18, 347
  • (2005) Bosch-Ramon, V., Aharonian, A. F., Paredes, J. M. 2005, A&A, 432, 609
  • (2006) Bosch-Ramon, V., Paredes, J. M., Romero, G. E., Ribó, M. 2006, A&A, 459, L25
  • (2008) Bosch-Ramon, V., Khangulyan, D., Aharonian, F. A. 2008, A&A, 489, L21
  • (2011) Bosch-Ramon, V., Perucho, M., Bordas, P., 2011, A&A, 528, 89
  • (1975) Castor, J., McCray, R., Weaver, R. 1975, ApJ, 200, L107
  • (2005a) Casares, J., Ribó, M., Ribas, I., et al., 2005a, MNRAS, 364, 899
  • (2005b) Casares, J., Ribas, I., Paredes, J. M., Martí, J., Allende Prieto, C. 2005b, MNRAS, 360, 1105
  • (2006) Chernyakova, M., Neronov, A., Walter, R. 2006, MNRAS, 372, 1585
  • (2001) Clark, J. S., Reig, P., Goodwin, S. P. 2001, A&A, 376, 476
  • (2011) Corbet, R. H. D., Cheung, C. C., Kerr, M., 2011, ATel, 3221, 1
  • (1972) Cox, D.P., 1972, ApJ, 178, 159
  • (2010) de Ona Wilhelmi, E. 2010, 38th COSPAR Scientific Assembly, Bremen, Germany, p.6
  • (2006) Dhawan, V., Mioduszewski, A., Rupen, M. 2006, Proceedings of the VI Microquasar Workshop: Microquasars and Beyond, Como, Italy, 52.1
  • (1983) Drury, L.O’C., 1983, Reports on Progress in Physics, 46, 973
  • (2006) Dubus, G. 2006, A&A, 456, 801
  • (2011) Durant, M., Kargaltsev, O., Pavlov, G. G., et al. 2011, ApJ, 735, 58
  • (2010) Falcone, A. D., Grube, J., Hinton, J., et al. 2010, ApJ, 708, L52
  • (1991) Frail, D. A., & Hjellming, R. M. 1991, AJ, 101, 2126
  • (2005) Gallo, E., Fender, R., Kaiser, C., et al. 2005, Nature, 436, 819
  • (2003) Heindl, W. A., Tomsick, J. A., Wijnands, R., Smith, D. M. 2003, ApJ, 588, L97
  • (2002) Heinz, S., 2002, A&A, 388, L40
  • (2002) Heinz, S. & Sunyaev, R. 2002, A&A, 390, 751
  • (2011) Hill, A. B., Szostek, A., Corbel, S., et al. 2011, MNRAS, 415, 235
  • (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101
  • (1981) Hutchings, J. B., Crampton, D. 1981, PASP, 93, 486
  • (2008) Khangulyan, D. V., Aharonian, F. A., Bogovalov, S. V., Koldoba, A. V., Ustyugova, G. V. 2008, IJMPD, 17, 1909
  • (1992) Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37
  • (2011) McSwain, M. V, Ray, P. S., Ransom, S. M., et al. 2011, ApJ, 738, 105
  • (1995) Manchester, R. N., Johnston, S., Lyne, A. G., D’Amico, N., Bailes, M., Nicastro, L. 1995, ApJ, 445, L1374
  • (1981) Maraschi, L. & Treves, A. 1981, MNRAS, 194, 1
  • (1996) Martí, J., Rodríguez, L. F., Mirabel, I. F., Paredes, J. M. 1996, A&A, 306, 449
  • (1998) Martí, J., Peracaula, M., Paredes, J. M., Massi, M., Estalella, R. 1998, A&A, 329, 951
  • (2005) Martocchia, A., Motch, C., Negueruela, I., 2005, A&A, 430, 245
  • (2001) Massi, M., Ribó, M., Paredes, J. M., Peracaula, M., Estalella, R. 2001, A&A, 376, 217
  • (2004) Massi, M., Ribó, M., Paredes, J. M., et al. 2004, A&A, 414, L1
  • (2009) Massi, M., Kaufman Bernadó, M. 2009, ApJ, 702, 1179
  • (2011a) Moldón, J., Johnston, S., Ribó, M., Paredes, J. M., Deller, A. T. 2011a, ApJ, 732, L10
  • (2011b) Moldón, J., Ribó, M., Paredes, J. M. 2011b, A&A, 533, L7
  • (2011c) Moldón, J., Ribo, M., Paredes, J. M. 2011c, High-Energy Emission from Pulsars and their Systems, Astrophysics and Space Science Proceedings, ISBN 978-3-642-17250-2. Springer-Verlag Berlin Heidelberg, 2011, p. 1-2
  • (2009) Muñoz-Arjonilla, A. J., Martí, J., Combi, J. A. 2009, A&A, 497, 457
  • (2011) Negueruela, I., Ribó, M., Herrero, A., Lorenzo, J., Khangulyan, D., Aharonian, F. A. 2011, ApJ, 732, L11
  • (2011) Okazaki, A. T., Nagataki, S., Naito, T., et al. 2011, PASJ, in press [astro-ph/1105.1481]
  • (2000) Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340
  • (2002) Paredes, J. M., Ribó, M., Martí, J., 2002, A&A, 393, 99
  • (2006) Paredes, J. M., Bosch-Ramon, V., Romero, G. E., 2006, A&A, 451, 259
  • (2007) Paredes, J. M., Ribó, M., Bosch-Ramon, V., et al. 2007, ApJ, 664, L39
  • (2003) Porter, J. M. & Rivinius, T. 2003, PASP, 115, 1153
  • (2011a) Pavlov, G. G., Chang, C., Kargaltsev, O. 2011a, ApJ, 730, 2
  • (2011b) Pavlov, G. G., Misanovic, Z., Kargaltsev, O, Garmire, G. P. 2011b, ATel, 3228
  • (2009) Pittard, J. M. 2009, MNRAS, 396, 1743
  • (2009) Puls, J., Sundqvist, J. O., Najarro, F., Hanson, M. M. 2009, AIPC, 1171, 1234
  • (2010) Rea, N., Torres, D. F., van der Klis, M., Jonker, P. G., Méndez, M., Sierpowska-Bartosik, A. 2010, MNRAS, 405, 2206
  • (2011) Rea, N., Torres, D. F., Caliandro, G. A., et al. 2011, MNRAS, 416, 1514
  • (2002) Ribó, M., Paredes, J. M., Romero, G. E., et al. 2002, A&A, 384, 954
  • (2011) Ribó, M., 2011, invited talk at HEPROIII, held in Barcelona, Spain, in June 2011
  • (2005) Romero, G. E., Christiansen, H. R., Orellana, M. 2005, ApJ, 632, 1093
  • (2007) Romero, G. E., Okazaki, A. T., Orellana, M., Owocki, S. P. 2007, A&A, 474, 15
  • (2007) Russell, D. M., Fender, R. P., Gallo, E., Kaiser, C. R. 2007, MNRAS, 376, 1341
  • (2010) Sabatini, S., et al. 2010, ApJ, 712, L10
  • (2008) Sánchez-Sutil, J. R., Martí, J., Combi, J. A., et al. 2008, A&A, 479, 523
  • (2011) Sarty, G. E., Szalai, T., Kiss, L. L., et al. 2011, MNRAS, 411, 1293
  • (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press)
  • (2008) Sierpowska-Bartosik, A. & Torres, D. F. 2008, Astropart. Phys., 30, 239
  • (2009) Sierpowska-Bartosik, A. & Torres, D. F. 2009, ApJ, 693, 1462
  • (2009) Skilton, J. L., Pandey-Pommier, M., Hinton, J. A., 2009, MNRAS, 399, 317
  • (1992) Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265
  • (2010) Tam, P. H. T., Hui, C. Y., Huang, R. H. H., et al. 2010, ApJ, 724, L207
  • (1992) Taylor, A. R., Kenny, H. T., Spencer, R. E., Tzioumis, A. 1992, ApJ, 395, 268
  • (1997) Tavani, M., Arons, J., 1997, ApJ, 477, 439
  • (2009a) Tavani, M., et al. 2009a, ApJ, 698, L142
  • (2009b) Tavani, M., et al. 2009b, Nature, 462, 620
  • (2000) Velázquez, P. F. & Raga, A. C. 2000, A&A, 362, 780
  • (1998) Wex, N., Johnston, S., Manchester, R. N., Lyne, A. G., Stappers, B. W., Bailes, M. 1998, MNRAS, 298, 997
  • (2010) Williams, S. J., Gies, D. R., Matson, R. A. et al. 2010, ApJ, 723, L93
  • (2011a) Zabalza, V., Bosch-Ramon, V., Paredes, J. M. 2011, ApJ, submitted
  • (2011b) Zabalza, V., Paredes, J. M., Bosch-Ramon, V. 2011, A&A, 527, 9
  • (2009) Zacharias, N., et al. 2009, VizieR Online Data Catalog, 1315, 0
  • (2008) Zavala, J., Velázquez, P. F., Cerqueira, A. H., Dubner, G. M. 2008, MNRAS, 387, 839
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description