Key Words.:

Abstract

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.

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

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 (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.

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 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 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.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 ((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)):

η=Lsd˙Mwvwc≈0.06Lsd37˙M−1w−6.5v−1w8.5,

(1)

where ˙Mw−6.5=(˙Mw/3×10−7M⊙yr−1) is the stellar
mass-loss rate, vw8.5=(vw/3×108cms−1) is the stellar wind velocity, and Lsd37=(Lsd/1037ergs−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,

Rp=√ηRorb(1+√η),

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

α≈28.6∘(4−η2/5)η1/3.

(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:

LsdLw=2Lsd˙Mwv2w≈12Lsd37˙M−1w−6.5v−2w8.3.

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

η≈(2GM(Rorbv2w)2)∼2×10−6R−2orb12.5v−4w8.5,

(5)

where M is the pulsar mass, and Rorb12.5=(Rorb/3×1012cm).

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 η<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).

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. ˙Mw>2.5×10−8M⊙ 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.

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Ωvwt 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
P⊥≈ρwv2⊥≈(2Ωt)2ρwv2w≈(2Ωt)2Pw (Pw=ρwv2w).
From this, the location where the
pulsar wind total pressure is balanced by the stellar wind can be
derived with

Lsd4πcx2=P⊥=ρw(4πT)2x2,

(6)

where ρw=ρw0/(1+x/Rorb)2, with
ρw0=˙Mw/4πR2orbvw
the stellar wind density at the pulsar location. For x≪Rorb,
Eq. (6) can be simplified to

x=(LsdvwT2R2orb(4π)2c˙Mw)1/4;

(7)

and for x≫Rorb,

x=(LsdvwT2(4π)2c˙Mw)1/2.

(8)

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

x0≈7×1012L1/2sd37v1/2w8.5T6˙M−1/2w−6.5cm,

(9)

where T6=(T/106s).
This distance is much less than in the purely ballistic case, for which
x0∼cT=3×1016T6 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

vexp∼√2Lsd˙Mw≈109L1/2sd37˙M−1/2w−6.5 cm s−1.

(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).

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 ∼vwT. 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.

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

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/1051erg) and nism1=(nism/10cm−3).
In this context, the binary system will leave the SNR after a time

tsp≈(2.8ESNRR2cρism)1/5v−7/5p;

(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.

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

RSNR≈1.15(ESNR/ρism)1/5t2/5p,

(12)

where tp is the pulsar age, velocity

vSNR≈(2/5)RSNR/tp,

(13)

and the inner pressure is

PSNR≈ESNR2πR3SNR,

(14)

assumed here to be constant.
Equilibrium between
bubble and SNR pressures
determines the location of the boundary between these two regions:
Rbcd≈(3Lsdtp/10πρismv2b)1/3,
and the equilibrium
between bubble preshock and postshock total pressures determines the
location of the bubble reverse
shock:

Rbrs≈√Lsd2πPSNRvexp.

(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

RSNR≈(2.8ESNRR2ct2ρism)1/7,

(16)

a velocity of

vSNR=(2/7)RSNR/tp,

(17)

and an inner pressure of

PSNR≈12πϵESNRR2cR5SNR,

(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

Rb≈0.76(Lsdρism)1/5t3/5p,

(19)

moving with velocity

vb≈(3/5)Rb/tp.

(20)

The inner pressure is now

Pb∼ρismv2b,

(21)

and again the preshock and postshock total
pressure balance gives the location of the bubble reverse shock: Rbrs≈√Lsd/2πPbvexp. 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:

Rbrs∼√Lsd2πvexpρismv2p≈0.3 pc.

(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.

Parameter [units]

value

Star mass-loss rate [M⊙ yr−1]

3×10−7

Stellar wind velocity [cm s−1]

2×108

Pulsar spin-down luminosity [erg s−1]

1037

SNR energy [erg]

1051

ISM number density [cm−3]

10

Magnetic-to-total energy density ratio

0.01, 0.1

Nonthermal-to-total luminosity ratio

0.01

Star luminosity [erg s−1]

1039

Stellar photon energy [eV]

10

Reverse shock velocity [cm s−1]

109

Source distance [kpc]

3

Table 1: Adopted parameter valuesFigure 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).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).

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=(1/2π)(vexp/c)2qBc
(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 ξBPSNR|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, Emax≈√qc/asB, with as=1.6×10−3, and the
diffusive-escape one, Emax≈√3/4π(vb/c)qBRbrs. 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

R∼Rbrs(1+3vexptcool/4Rbrs)1/3∼10Rbrs,

(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.01Lsd. 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 (∝BE2) decreases with time, and it
quenches the X-ray flux for tp\ga105 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=104yr<tc.

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

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\la104 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\la104 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 arcminutes 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 ∼vexpT. 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.

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.

Acknowledgements.

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.

Footnotes

offprints: V. Bosch-Ramon;

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.

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.