EDP Sciences
Free Access
Issue
A&A
Volume 535, November 2011
Article Number A20
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201117235
Published online 26 October 2011

© ESO, 2011

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. 2006; 2007; Aharonian et al. 2005a,b; Abdo et al. 2009a,b,c, 2010a,b, 2011; Tavani et al. 2009a,b; Sabatini et al. 2010; Williams et al. 2010; Tam et al. 2010; Hill et al. 2011; Falcone et al. 2011; Corbet et al. 2011).

thumbnail Fig. 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 η < 1) 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).

Open with DEXTER

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 B1259 − 63, and the gamma-ray emitter candidate Cygnus X-1 (Heindl et al. 2003; Sánchez-Sutil et al. 2008; Paredes et al. 2007; Pavlov et al. 2011a; Durant et al. 2011; Martí et al. 1996; Gallo et al. 2005; Russell et al. 2007).

Cygnus X-1 and Cygnus X-3 are high-mass microquasars, in which accretion leads to jets that can interact with the interstellar medium (ISM; e.g. Velázquez & Raga 2000; Heinz 2002; Heinz & Sunyaev 2002; Bosch-Ramon et al. 2005; Zavala et al. 2008; Bordas et al. 2009; Bosch-Ramon et al. 2011). On the other hand, PSR B1259 − 63 is a system formed by an O9.5 V star and a nonaccreting millisecond pulsar (Johnston et al. 1992; Negueruela et al. 2011). 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 1997), with its radio emission extending far away from the binary (Moldón et al. 2011a). PSR B1259 − 63 also presents extended X-ray emission on scales of  ~ 4 − 15′′ (i.e. a projected linear size of  ~ 1.5 − 5 × 1017 cm at the source distance of 2.3 kpc; Negueruela et al. 2011), with a flux  ~ 10-13 erg s-1 cm-2 and photon index  ~ 1.6 (Pavlov et al. 2011a). 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. 2005a,b; Sarty et al. 2011; Rea et al. 2010; 2011; McSwain et al. 2011; Ribó 2011; see also Bosch-Ramon & Khangulyan 2009). 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. 1992; Paredes et al. 2000; Massi et al. 2001; 2004; Paredes et al. 2002; Romero et al. 2005; Paredes et al. 2006; Bosch-Ramon et al. 2006; Massi & Kaufman Bernadó 2009). 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 1981; Martocchia et al. 2005; Dubus 2006; Chernyakova et al. 2006; Sierpowska-Bartosik & Torres 2008; Sierpowska-Bartosik & Torres 2009; see however Romero et al. 2007). Both binaries are high-mass systems, with LS 5039 harboring an O6.5 V star (Clark et al. 2001), and LS I +61 303, an early Be star (Hutchings & Crampton 1981).

Extended X-rays have been found in LS 5039, with angular size θ ~ 1′−2′ (~2−4 × 1018 cm at 2.5 kpc; Casares et al. 2005a), X-ray flux  ≈ 9 × 10-14 erg s-1 cm-2 and photon index  ~1.9−3.1 (Durant et al. 2011), and possibly also in LS I +61 303, with θ ~ 10′′ (~3 × 1017 cm at 2 kpc; Frail & Hjellming 1991), and X-ray flux  ≈ 2    ×    10-14 erg s-1 cm-2 (Paredes et al. 20071). It is noteworthy that the gamma-ray binaries HESS J0632 + 057 (Aharonian et al. 2007; Hinton et al. 2009; Skilton et al. 2009; Falcone et al. 2010; 2011; Moldón et al. 2011b) and 1FGL J1018.6 − 5856 (Corbet et al. 2011; Pavlov et al. 2011b; de Ona Wilhelmi et al. 2010) 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 J0632 + 057 and 1FGL J1018.6 − 5856, 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 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 B1259−63, 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 Fig. 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 (1997). 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. 2008): (1)where w − 6.5 = (w/3 × 10-7   M   yr-1) is the stellar mass-loss rate, vw8.5 = (vw/3 × 108   cm   s-1) is the stellar wind velocity, and Lsd37 = (Lsd/1037   erg   s-1) 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-1, is much lower than vw, typically  ~2 × 108 cm s-1 for massive stars (e.g. Puls et al. 2009), 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 2002, Bogovalov et al. 2011).

For reasonable η-values  < 1 (a scenario with η > 1 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, (2)which is  ≈ 7 × 1011 cm for η = 0.1 and a distance between the star and the pulsar of Rorb = 3 × 1012 cm. The distance between the CD and the star is then  ≈ 2.3 × 1012 cm. The CD, which separates the shocked stellar and pulsar winds, has an asymptotic half-opening angle (see Bogovalov et al. 2008 and references therein): (3)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: (4)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 Rp is equal to the gravitational capture radius (Bondi & Hoyle 1944): (5)where M is the pulsar mass, and Rorb12.5 = (Rorb/3 × 1012   cm).

2.2. Shocked flow evolution

After being shocked, the pulsar wind facing the star moves at a speed c/3, 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. 2008). Because of the very low density in the downstream region, the pulsar wind shock is adiabatic (Bogovalov et al. 2008; Khangulyan et al. 2008). The stellar wind shock could be either adiabatic or radiative, but in either case the stellar wind velocity will be  ≪ c. Therefore, the two shocked winds have a relative velocity  ~ c/3 − c, 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. 1992; Pittard 2009). Perturbations generated in the interface will initially move with the stellar shocked wind, growing on a timescale tKH ~ l/csw, where l is the size of the perturbation, and csw the sound speed of the shocked stellar wind. Taking csw ~ vw ~ 108 cm s-1 and l ~ Rp = 1012 cm, one obtains tKH ~ 104 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 η < 1, (e.g. w > 2.5 × 10-8   M yr-1 for vw ~ 2 × 108 cm s-1 and Lsd ~ 1037 erg s-1); more briefly, the extreme case η > 1; 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 Ω = 2π/T, where T 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 ac =  −2   Ω × v. The stellar wind moves along the X-axis with speed  ~ vw, acquiring a velocity v ⊥  = 2   Ω   vw   t in the Y-direction. In the frame of the shocked-wind flow, t = x/vw, and x 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 (6)where ρw = ρw0/(1 + x/Rorb)2, with the stellar wind density at the pulsar location. For x ≪ Rorb, Eq. (6) can be simplified to (7)and for x ≫ Rorb, (8)For typical parameters, the solution is closer to the one given by Eq. (8), therefore the bending distance can be written as (9)where T6 = (T/106   s). This distance is much less than in the purely ballistic case, for which x0 ~ c   T = 3 × 1016   T6 cm.

It is noteworthy that the pulsar wind is shocked even in the direction opposite to the star, within a distance  ~ x0. 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. 2008).

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. 2011.) 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 (10)Typically, the mas radio emission in pulsar massive binaries will come from scales larger than x0, 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  ~ vw   T. 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 η < 1, for η > 1 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. 2009). A sketch of the two cases discussed, with η > 1 and  < 1, is shown in Fig. 2.

thumbnail Fig. 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 η > 1, is presented. To the left, the case with η < 1.

Open with DEXTER

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  ~ 1 km s-1 (e.g. Porter & Rivinius 2003). 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. 2007 and Okazaki et al. 2011). 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 ESNR ~ 1051 erg and ISM density ρism = 1.7 × 10-23 g cm-3 (nism = 10 cm-3), the SNR becomes radiative (Blinnikov et al. 1982) after tc ≈ 2 × 104 yr when achieving the radius Rc ≈ 11   (ESNR51/nism1)1/3 pc, where ESNR51 = (ESNR/1051   erg) and nism1 = (nism/10   cm-3). In this context, the binary system will leave the SNR after a time (11)i.e., tsp ≈ 8 × 104 yr for vp ≈ 2 × 107 cm s-1.

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  ~vexp 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 Lnt as high as Lsd. 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. 2008), 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 1959). Later on, the SNR forward shock in the ISM becomes radiative, entering into the so-called snow-plow phase (Cox 1972; Blinnikov et al. 1982). 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 1959) (12)where tp is the pulsar age, velocity (13)and the inner pressure is (14)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: (15)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 t > tc, the SNR enters into the radiative phase, so the snow-plow solution has to be used (Blinnikov et al. 1982). The SNR/ISM forward shock has now a radius of (16)a velocity of (17)and an inner pressure of (18)where ϵ = 0.24. 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. 1975). This renders the ISM forward shock location at (19)moving with velocity (20)The inner pressure is now (21)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 tp > 105 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: (22)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 B1259 − 63, although the comparison is less straightforward, as discussed in Sect. 5. It is noteworthy that, as long as x0 is much smaller than the largest scales of interaction with the medium, the size difference between PSR B1259 − 63 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 2 × 104 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. 1982).

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.

Table 1

Adopted parameter values.

thumbnail Fig. 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  × 105 yr. b) Evolution with time of the radii of the bubble/ISM forward (black solid line) and reverse shocks (blue dotted line).

Open with DEXTER

thumbnail Fig. 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 lines) and bubble/ISM magnetic fields (red dotted/green dot-dashed lines).

Open with DEXTER

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: Ė = (1/2π)(vexp/c)2   q   B   c (e.g. Drury 1983), where B is the magnetic field and q the particle charge. The same B 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  ∝ E-2. The B-energy density was taken equal to ξB   PSNR|b, with ξB = 0.1 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 as = 1.6 × 10-3, and the diffusive-escape one, . Although the stellar photon energy density, the dominant radiation field, is higher than the B-energy density at Rbrs, the Klein-Nishina (KN) effect makes IC losses less important than synchrotron (or escape) ones at Emax. For a wide source age range, in the case of electrons Emax is  ~ 100 TeV, but determined by synchrotron cooling for ξB = 0.1, and by diffusive escape for ξB = 0.01. For protons, the maximum energies are  ~300 TeV for ξB = 0.1 and 100 TeV for ξB = 0.01, 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  ∝ B − 1/2, but diffusive-escape ones are constant, because B × Rbrs 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 1970), we modeled the emitter as only one zone, taking homogeneous B- and radiation fields. This approach is valid for the synchrotron calculations as long as the flow is subsonic, and ξB 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  ∝ 1/R2, where R is the distance to the reverse shock. This makes particles accumulate at (23)with tcool the relevant cooling time. Most of the IC radiation comes also from a distance  ~R to the star. This is not true for very high-energy IC photons, but for the adopted B-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 1039 erg s-1 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: Lnt = 0.01   Lsd. Although Lnt 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 RSNR|b/vSNR|b.

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 104 yr (SNR) and 3 × 104 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 1′′ × 1′′ units at 3 kpc, taking R 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 R. The general evolution of R 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 ξB as well. For ξB = 0.1, X-ray emitting electrons radiate all their energy through synchrotron, rendering a flat X-ray lightcurve and photon indexes  ~2. For ξB = 0.01, with the electron maximum energies limited by diffusive-escape, the maximum energy of synchrotron photons ( ∝ B   E2) decreases with time, and it quenches the X-ray flux for tp ≳ 105 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 t = 104   yr < tc.

thumbnail Fig. 5

Spectral energy distributions of the synchrotron (thick lines) and IC (thin lines) emission computed for the SNR (104 yr) and the bubble/ISM (3 × 104 yr) cases, and ξB = 0.01,   0.1.

Open with DEXTER

thumbnail Fig. 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 tp = 3 × (103 − 105) yr and for ξB = 0.01,   0.1.

Open with DEXTER

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 θ ≈ 10′′ when still inside the SNR (t ≲ 104 yr), and  ~ 20′′ − 1′ in the bubble/ISM case (t ≈ 104 − 105 yr). These θ values were calculated for the ρism, ESNR, and Lsd values given in Table 1, although we need to recall their weak inter-dependence: θ ∝ (ESNR   |   Lsd/ρism)1/5. Extended X-ray fluxes could be detectable by an instrument like Chandra up to tp ~ 105 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 ξB = 0.1, it may be hard to be found for the case ξB = 0.01, and its detection seems discarded for the bubble/ISM case for any ξB-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 Lnt, and reducing strongly ξB 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  ~ 10 GeV introduced by adiabatic losses, the bubble reverse shock would be a good target for a  ~ 10 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 B1259 − 63 can be understood in the scenario posed here.

5.1. LS 5039

LS 5039 is likely a young source with an age range (4−10) × 104 yr, and it has already abandoned its SNR (e.g. Ribó et al. 2002; Moldón et al. 2011c). The angular size at 2.5 kpc for the bubble/ISM case is consistent with the observed angular size of the emitter, θ ~1′−2′, and the fluxes can be explained in the framework given here. The photon index at distances  ~1′ is  ~1.9 ± 0.7, consistent with the SED shown in Fig. 5, and then softens farther out reaching  ~3.1    ±    0.5 at  ~2′ (Durant et al. 2011). 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 nism = 10 cm-3 is similar to the values derived for the surroundings of the LS 5039 SNR candidate (Ribó et al. 2002), and the pulsar Lsd is probably  ~1037 erg s-1, which is required to explain the GeV luminosity of the source (see Sect. 4 in Zabalza et al. 2011a). The morphology of the extended radio emission found by Durant et al. (2011) 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 η > 1 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, θ ~ 10′′, is compatible with the bubble/SNR scenario for tp ≲ 104 yr, but with a low nonthermal efficiency or a very low magnetic field, given the weak X-ray flux (Paredes et al. 2007) and the radio nondetection on those scales (Martí et al. 1998).

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 arcmin and fluxes of tens of mJy (Martí et al. 1998; Muñoz-Arjonilla et al. 2009). 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. 1998; Muñoz-Arjonilla et al. 2009). However, the radio emission surrounding LS I +61 303 may be mostly thermal (Muñoz-Arjonilla et al. 2009), coming from an SNR/ISM shock close to its radiative phase. If the nism values in the region of LS I +61 303 were a bit higher than those adopted here, an SNR remnant with energy ESNR ~ 1051 erg would enter its radiative phase after  ~104 yr, with a radius  ~1019 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. 2006) is compatible with LS I +61 303 being at the center of the SNR (see however Martí et al. 1998 for image deformation). As in LS 5039, the Lsd-value in LS I +61 303 is expected to be  ~1037 erg s-1 because of the high GeV luminosity of the source (see Sect. 5 in Zabalza et al 2011b).

5.3. PSR B1259 − 63

PSR B1259 − 63 has an age of tp ≈ 3 × 105 yr (Wex et al. 1998), and is likely now interacting directly with the ISM. With this age, and the pulsar spin-down luminosity Lsd ≈ 8 × 1035 erg s-1 (Manchester et al. 1995), Eq. (19) yields Rbrs ~ 2 × 1018 cm (~1′), assuming nism = 10 cm-3. This is about one arcminute at 2.3 kpc and much larger than the observed θ ~ 4′′−15′′ (Pavlov et al. 2011a). Nevertheless, in PSR B1259 − 63 it seems likely that the proper motion has already bow-shaped the ISM shock. In that case, for nism = 10 cm-3 and vp ~ 107 cm s-1, the reverse shock would be located at a distance of  ~3 × 1017 cm from the system, or  ~10′′ neglecting projection effects. Unfortunately, the large errors of the proper motion measurements (Zacharias et al. 2009) 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 B1259 − 63. The most significant detection comes from  ~ 4′′, or  ~1017 cm in linear (projected) size, which is  ~vexp   T. Although the flow bending distance, x0, 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 B1259 − 63.

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 B1259 − 63, 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 B1259 − 63, 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 J0632 + 057 and 1FGL J1018.6 − 5856, 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 B1259 − 63 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.


1

Rea et al. (2010) did not find evidence of extended X-rays in LS I +61 303 in longer observations than those studied in Paredes et al. 2007. 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.

Acknowledgments

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.

References

All Tables

Table 1

Adopted parameter values.

All Figures

thumbnail Fig. 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 η < 1) 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).

Open with DEXTER
In the text
thumbnail Fig. 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 η > 1, is presented. To the left, the case with η < 1.

Open with DEXTER
In the text
thumbnail Fig. 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  × 105 yr. b) Evolution with time of the radii of the bubble/ISM forward (black solid line) and reverse shocks (blue dotted line).

Open with DEXTER
In the text
thumbnail Fig. 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 lines) and bubble/ISM magnetic fields (red dotted/green dot-dashed lines).

Open with DEXTER
In the text
thumbnail Fig. 5

Spectral energy distributions of the synchrotron (thick lines) and IC (thin lines) emission computed for the SNR (104 yr) and the bubble/ISM (3 × 104 yr) cases, and ξB = 0.01,   0.1.

Open with DEXTER
In the text
thumbnail Fig. 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 tp = 3 × (103 − 105) yr and for ξB = 0.01,   0.1.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.