Thermal X-ray Emission from the Shocked Stellar Wind of Pulsar Gamma-ray Binaries

Gamma-ray loud X-ray binaries are binary systems that show non-thermal broadband emission from radio to gamma rays. If the system comprises a massive star and a young non-accreting pulsar, their winds will collide producing broadband non-thermal emission, most likely originated in the shocked pulsar wind. Thermal X-ray emission is expected from the shocked stellar wind, but until now it has neither been detected nor studied in the context of gamma-ray binaries. We present a semi-analytic model of the thermal X-ray emission from the shocked stellar wind in pulsar gamma-ray binaries, and find that the thermal X-ray emission increases monotonically with the pulsar spin-down luminosity, reaching luminosities of the order of 10^33 erg/s. The lack of thermal features in the X-ray spectrum of gamma-ray binaries can then be used to constrain the properties of the pulsar and stellar winds. By fitting the observed X-ray spectra of gamma-ray binaries with a source model composed of an absorbed non-thermal power law and the computed thermal X-ray emission, we are able to derive upper limits on the spin-down luminosity of the putative pulsar. We applied this method to LS 5039, the only gamma-ray binary with a radial, powerful wind, and obtain an upper limit on the pulsar spin-down luminosity of ~6x10^36 erg/s. Given the energetic constraints from its high-energy gamma-ray emission, a non-thermal to spin-down luminosity ratio very close to unity may be required.


INTRODUCTION
In the past few years, gamma-ray binaries have emerged as a new category of very high-energy (VHE, E > 100 GeV) gamma-ray sources. In most of the cases, they comprise a stellar mass compact object orbiting a young, massive star. They generally differ from typical X-ray binaries in that they exhibit luminosities at gamma-ray energies at levels similar or above their X-ray luminosity. Five such sources have been already detected: LS 5039 (Aharonian et al. 2005b), PSR B1259−63 (Aharonian et al. 2005a), LS I +61 303 (Albert et al. 2006), HESS J0632+057 (Skilton et al. 2009;Bongiorno et al. 2011;Moldón et al. 2011), and 1FGL J1018.6−5856 (Corbet et al. 2011;Pavlov et al. 2011). Evidence of VHE flaring emission has been found in Cygnus X-1 (Albert et al. 2007), but at present lacks confirmation.
Whereas the pulsar nature of the compact object in PSR B1259−63 is determined thanks to the detection of radio pulses, this is not possible in the case of LS 5039 and LS I +61 303, in which the radio pulses would be likely freefree absorbed in the stellar wind. Broadband radiation modeling in the microquasar scenario has been able to reproduce the spectral energy distribution of the source (e.g., Paredes et al. 2006;Bosch-Ramon et al. 2006), but the lack of accretion signatures, and the similarities with PSR B1259−63, prompted other authors to consider the pulsar wind shock (PWS) scenario for these sources (e.g., Martocchia et al. 2005;Dubus 2006;Chernyakova et al. 2006).
The emission from the non-thermal particle population accelerated in the PWS scenario has been widely studied in the radio, X-ray, and high energy (HE) and VHE gamma-ray bands (Maraschi & Treves 1981;Tavani & Arons 1997;Kirk et al. 1999;Dubus 2006;Khangulyan et al. 2007;Neronov & Chernyakova 2007;. Furthermore, the possibility of detecting the free pulsar wind from its spec-tral signature has also been discussed (e.g., Kirk et al. 1999;Khangulyan et al. 2007;Cerutti et al. 2009). However, one of the sources of emission that could provide valuable information on the characteristics of gamma-ray binaries has been hitherto overlooked: the thermal X-ray emission from the shocked stellar wind. Massive stellar binaries (of types O+O or O+WR) are known to produce powerful thermal X-ray emission in the wind collision region (WCR), where the supersonic winds of the two young stars interact (see, e.g., Pittard et al. 2005, for a review). Even though it has been previously considered that the detection of X-ray emission lines would be an indication of accretion in gamma-ray binaries, similar lines can arise from the thermal X-ray emission of the shocked stellar wind in the WCR. The complexity of the stellar winds from the Be companion stars in LS I +61 303 and PSR B1259−63 makes the analysis of the thermal emission from the shocked stellar wind difficult. The presence of both a polar wind and a decretion disk adds a number of unknown free parameters. On the other hand, the radial wind from the O star in LS 5039, also more powerful than the one from Be stars, turns this source into an ideal candidate to study the observational impact of the thermal wind emission.
In this paper we present a model for the thermal X-rays from pulsar gamma-ray binaries inspired on semi-analytical models of the X-ray emission from massive binary systems (e.g., Antokhin et al. 2004;Parkin & Pittard 2008). This model, which allows the study of the dependence of the thermal Xray emission for different stellar and pulsar wind properties, will be applied to the binary LS 5039 in order to constrain the spin-down luminosity of its putative pulsar. In Section 2, we characterize the shocked stellar wind region, approximated here as the contact discontinuity between the pulsar and the stellar winds, and present a model of its thermal X-ray emission in Section 3. In Section 4 we apply the model to LS 5039 2 Zabalza, Bosch-Ramon, & Paredes and present spin-down luminosity upper limits in Section 5. Finally, we discuss the implications of the results in Section 6.

THE DYNAMICAL MODEL
In the PWS scenario, the pulsar wind and the stellar wind collide forming an interaction region bounded on either side by the reverse shocks of the winds. The shocked winds are separated by a surface called contact discontinuity (CD), located where the ram pressure components perpendicular to this surface are in equilibrium: p ⊥ = p p⊥ . Both winds are assumed here to be radial.
The ram pressure of the relativistic pulsar wind can be represented at any point at a distance r p from the pulsar as p p (r p ) = L sd /(4πcr 2 p ), where L sd is the spin-down luminosity of the pulsar. In the case of the stellar wind, we will assume a β-velocity law in the radial direction, v where r is the distance to the center of the star, R is the radius of the star and v ∞ is the terminal velocity of the stellar wind. The azimuthal component of the wind will be neglected. We have considered β = 0.8, which has been shown to be a good approximation of the outer structure of the winds of hot massive stars (Pauldrach et al. 1986). Assuming a constant mass-loss rate from the starṀ, the density of the stellar wind is ρ =Ṁ/(4πr 2 v(r )). The ram pressure of the stellar wind is given by p (r ) =Ṁv(r )/(4πr 2 ).
The balance of perpendicular components of the ram pressure that defines the position of the CḊ can be rearranged to obtain the dimensionless parameter η where θ and θ p are the angles between the CD surface and the directions toward the center of the star and the pulsar, respectively. The coordinates x and y correspond to the position along the line of centers and the distance from it, respectively (see Figure 1). Following Antokhin et al. (2004), we can obtain the differential equation that describes the shape of the CD as a function of η from Equation (2), where x(y) is the function describing the shape of the CD and D is the binary separation distance. The boundary condition (dx/dy y→0 → 0) is x 0 = D/(1 + √ η), the location of pressure equilibrium along the symmetry axis, also known as stagnation radius. For the case of constant wind velocities, η will be constant, and x 0 and x(y) can be computed analytically. However, given the β-velocity law we adopted for the stellar wind, a numerical integration is required to obtain the shape of the CD since the stellar wind may not yet have reached its terminal velocity at the contact surface. Figure 1 shows the resulting shape of the CD in comparison to the star-pulsar system as a function of η ∞ ≡ (L sd /c)/(Ṁv ∞ ). At the scales of interest, and given the very high speed of the stellar wind of massive stars, the role of orbital motion in the shape of the CD is not expected to be significant. The aberration or skew angle µ of the shock cap is defined as the angle between the axis of symmetry of the shock cap and the C Star Pulsar P y x r p r ⋆ θ ⋆ θ p D Figure 1. Illustration of the dynamical model for a close pulsar gamma-ray binary. The ∼10 R star (but not the pulsar) is at scale. The thick black line indicates the shape of the CD for η ∞ = 0.04. The line CP is tangent to the CD at the point P, for which the angles and radii used in the derivation of the CD are shown. The gray lines show the shape of the CD for values of η ∞ spaced logarithmically in the range 0.003-0.3.
line of centers of the star and compact object. This angle is given by tan µ = v orb /v w , where v orb is the orbital velocity and v w is the speed of the slowest wind (i.e., the wind of the star) at the contact surface. For the orbital parameters of LS 5039, for example, the skew angle is below 23 • at all moments along the orbit.
3. THERMAL X-RAY EMISSION 3.1. Cooling regime of the shocked stellar wind On either side of the CD, the stellar and pulsar winds develop shocks. The shocked pulsar wind is the candidate location for particle acceleration, which would give rise to the observed broadband non-thermal emission (Tavani & Arons 1997). On the other side of the CD, the stellar wind develops a strong shock where the wind material is slowed down and heated. The cooling efficiency of the gas determines the width of the cooling layer. On one extreme, if radiative cooling is very efficient, the gas will be rapidly cooled after the shock and will collapse into a thin dense layer on the CD. The rapid cooling assures that all the wind kinetic energy flux perpendicular to the shock is transformed into thermal emission. On the other extreme, if radiative cooling is slow enough, the hot gas will flow along the CD and cool adiabatically through expansion. In the adiabatic regime, a smaller fraction of the kinetic energy of the stellar wind is emitted in the X-ray band. Stevens et al. (1992) proposed the parameter χ = t rad /t esc as a measure of the cooling regime of the shocked material in colliding wind binaries. χ values significantly below (above) unity indicate that the shocked stellar wind cools radiatively (adiabatically), whereas χ values around unity indicate that a fraction of the incoming kinetic energy will be radiated away and the rest will be dissipated through adiabatical expansion.

Estimation of the expected X-Ray luminosity
Before we enter into the details of the calculation of the thermal X-ray spectrum and luminosity, we will roughly derive the expected X-ray luminosity as a function of the parameters of the system. As mentioned above, we expect the X-ray luminosity to be a fraction f of the wind kinetic luminosity crossing the shock: L X = f L sh kin . The value of f is strongly dependent on the value of χ at the shock, and as a first approximation can be taken as f ≈ 1/(1 + χ) (see Equation (7)). Additionally, THERMAL X-RAY EMISSION FROM THE SHOCKED STELLAR WIND OF PULSAR GAMMA-RAY BINARIES 3 the kinetic luminosity deposited on the shock, L sh kin , will be a fraction of the total kinetic luminosity of the stellar wind, L tot kin , determined by the fraction of the total solid angle of the stellar wind subtended by the shock. For a cone with a half-opening angle of θ, we obtain L X ≈ f 1 2 (1−cos θ)L tot kin , which can be simplified up to O(θ 4 ) as L X ≈ f 1 4 θ 2 L tot kin . To obtain the estimate of the luminosity, we will here consider that the wind reaches the shock front at terminal velocity. The angle θ can then be derived from Equation (3) as πη/(1 + η), which on a first-order approximation for small values of η (i.e., dominance of the stellar wind) will be θ ≈ πη. Using these approximations we obtain L X ≈ f 1 4 π 2 η 2 L tot kin . As seen in Equation (2), η is proportional to L sd . Therefore, we expect the emitted thermal X-ray luminosity from the shocked stellar wind to behave roughly as: where we have taken the wind terminal velocity as v ∞ = 2400 km s −1 , and assumed that half of L sh kin is converted into X-ray luminosity, i.e., f = 0.5. This latter assumption is consistent with a cooling regime between the adiabatic and radiative regimes, i.e., χ ≈ 1.
The above approximations will only hold for small values of η, i.e., for the case of a dominant stellar wind. When the ram pressure of the pulsar wind is similar to that of the stellar wind and η approaches unity, the half-opening angle of the CD will approach π/2. At this limit, the X-ray luminosity will be even higher that for the previous case, and of the order of L X ≈ 3 × 10 34 Ṁ 10 −7 M yr −1 erg s −1 .
Gamma-ray binaries are non-thermal X-ray sources with luminosities of the order of several 10 33 erg s −1 . Therefore, as Equations (4) and (5) show, pulsar spin-down luminosities moderately higher than 10 36 erg s −1 would imply a significant thermal component in their X-ray spectrum.

X-Ray emission model details
For a strong, steady, standing shock with a compression ratio ρ 0 /ρ w = v w⊥ /v 0⊥ = 4 (where the subindices w and 0 indicate the properties of the gas immediately before and after the shock, respectively), the gas will have an immediate postshock temperature of where v w⊥ indicates the preshock speed of the wind component perpendicular to the shock front and µ is the average atomic weight. A value of µ = 0.62 is appropriate for a fully ionized medium with solar abundances (Anders & Grevesse 1989).
The X-ray emission of the shocked stellar gas is computed assuming that, in the emitting plasma, conditions are equal to those at the immediate postshock position, which is true as long as the gas remains subsonic. This will provide an estimate of how much of the kinetic energy of the incoming wind is radiated away. Assuming the shocked gas acts as a calorimeter, in which all the incoming perpendicular kinetic energy is either radiated away as X-ray emission or dissipated through adiabatic expansion, the fraction of the kinetic energy converted into X-ray emission will be where t rad is the radiative cooling timescale, t esc is the escape or flow dynamics timescale and L sh kin is the incoming kinetic energy. Given thatĖ ∝ t −1 , the fraction of times may be interpreted as the fraction of the radiated energy over the total energy lost owing to radiation and escape.
For a gas cooling with a volumetric emission rate Λ(T ), the radiative cooling timescale is given approximately by where n 0 is the postshock number density, and in a strong shock will be four times the preshock wind number density. Instead of using a pre-calculated table for Λ(T ), it can be calculated using the mekal code used to compute the emitted spectrum (see below). The escape timescale is the time the heated gas takes to flow away from the postshock location. Assuming that the stellar wind overpowers the pulsar wind, we can use the distance from the CD to the pulsar as the radius of curvature of the shock front, and thus as a typical flow dynamics distance scale. The escape timescale can then be approximated as t esc = r p /v 0 . We assume that the wind velocity component parallel to the shock is unaffected after traversing it, so that the immediate postshock speed may be found as v 2 0 = (v w⊥ /4) 2 + v 2 w . We used the mekal code as implemented in the spex software package (Kaastra et al. 1996) to generate a look-up table of emission spectra for optically thin hot gas in collisional ionization equilibrium as a function of the gas temperature T . We populated the look-up table with emission spectra in the 0.05-15 keV range, taking 250 logarithmically spaced temperatures from 10 4 to 10 9 K. Additionally, for each of the temperatures, the volumetric emissivity rate Λ(T ) was computed by integrating the emission between photon energies of 0.001 and 50 keV.
In order to compute the total cumulative emission from the shocked stellar wind we assumed that the X-ray emitter is infinitely thin and located on the CD defined through Equation (3). We defined a grid of 2000 bins in the y-direction, away from the line of centers, and 30 azimuthal bins on the shock front. For each segment of the grid the immediate postshock temperature T 0 and X-ray luminosity was obtained through Eqs. (6) and (7). Using the X-ray emission look-up table, we obtained the intrinsic X-ray emission spectrum from each segment. Finally, we computed the neutral hydrogen column density owing to the non-shocked stellar wind, from each of the positions of the segments on the shock front to the observer. We assumed that neither the enhanced density of the shocked stellar wind nor the free pulsar wind contribute significantly to the total column density. The column density of the former can be estimated to be two orders of magnitude below the wind column density (Szostek & Dubus 2011), and the latter is too rarefied to absorb X-rays. We then computed the energy-dependent attenuation factor corresponding to the found column density through the wabs subroutine provided with Xspec v12. In addition, the emission from those areal segments eclipsed by the star was discarded. The absorbed emission from each of the areal segments was added together to obtain the cumulative emitted spectrum of the shocked stellar wind.
From the calculations described in this section, we can see what the general characteristics of the shocked stellar wind emission will be. Along the line of centers of the system we will find both the highest temperatures and the highest densities owing to the perpendicular impact of the wind on the shock front. The hardest X-rays will come from this location, but its volume is relatively small and its luminosity may not dominate the final emission. Away from this point temperatures will drop rapidly as the incoming wind will not be perpendicular to the shock, whereas the volume increases, and density decreases. This region will produce a softer X-ray spectrum. The final spectrum will be a so-called multi-temperature X-ray spectrum, resulting from the emission of the stratified temperature structure of the shocked gas.

Caveats of the model
The semi-analytic model detailed in the previous sections is robust given that the assumptions made are fairly conservative. However, there are some factors that affect the WCR that cannot be accounted for in our semi-analytic framework. Detailed hydrodynamical simulations would be required to take into account effects such as the following.
• In the radiative cooling regime the CD may be disrupted by thin-shell instabilities (Stevens et al. 1992;Pittard 2009). Additionally, stellar wind clumping would likely disturb the shocked wind region to some extent. If this perturbation is propagated to the shocked pulsar wind, it would have implications for both the thermal and nonthermal emission.
• We modeled the radiative output of the postshock gas through the ratio of radiative to total cooling timescales in Equation (7), and assumed an infinitely thin cooling layer. Both the intrinsic X-ray luminosity and cooling layer width could be calculated self-consistently through hydrodynamical simulations. Whereas the results are likely to be very similar in the radiative cooling regime limit, they will probably differ in the adiabatic limit. In the adiabatic limit of the postshock flow, the flow may become supersonic at a distance smaller than r p , cooling the gas and reducing the radiative outcome. This effect would not be important at scales close to the apex, but could affect the contribution produced at the periphery of the system.

Method of spin-down luminosity upper limit derivation
To accurately compare the output of the model described in Sections 2 and 3 with the observed X-ray spectra of gammaray binaries, the computed X-ray spectra can be formatted as a FITS table model in order to load it into the spectral analysis software package Sherpa (Freeman et al. 2001). The FITS table model is generated from the cumulative spectra for 25 pulsar spin-down luminosities spaced logarithmically between 10 36 erg s −1 and 5 × 10 37 erg s −1 , taking a constant value of the mass-loss rate and orbital inclination.
The X-ray spectra observed at a certain orbital phase can be assigned a source model corresponding to an interstellar medium neutral hydrogen absorption model (A), affecting both the non-thermal power law (P) and the computed shocked stellar wind thermal emission (Th): A ISM × (P + Th). The goal is then to find for which value of the spin-down luminosity the thermal component begins to affect the shape of the spectrum and the fit worsens significantly. We used the confidence Sherpa command to evaluate the 99.7% (3σ) confidence level range for acceptable values of the spin-down luminosity. If no thermal features are present in the observed spectrum, the lower bound will be consistent with zero, i.e., lack of thermal component. The upper bound, on the other hand, provides us with the upper limit to the spin-down luminosity of the pulsar.

APPLICATION OF THE MODEL TO LS 5039
The system LS 5039 is located at 2.5±0.1 kpc and contains a compact object with a mass between 1.4 and 5 M , orbiting an O6.5V((f)) donor star every 3.90603 ± 0.00017 days (Casares et al. 2005) in a mildly eccentric orbit (e = 0.24 ± 0.08; Sarty et al. 2011). The detection of elongated asymmetric emission in high-resolution radio images was interpreted as mildly relativistic ejections from a microquasar jet and prompted its identification with an EGRET HE gamma-ray source (Paredes et al. 2000(Paredes et al. , 2002. However, recent Very Long Baseline Array observations by Ribó et al. (2008) show morphological changes on short timescales that might be consistent with a pulsar binary scenario. It shows variable periodic emission in the X-ray, and HE and VHE gamma-ray bands. The X-ray lightcurve has an orbital modulation ), superposed by short timescale features that are quite stable over the years ). The maximum and minimum emission phases in the periodic X-ray lightcurve are located around the inferior and superior conjunctions, so a geometrical effect for the modulation has been suggested. However, inferior and superior conjunctions are located close to apastron and periastron, and the broadness of the peaks, well fitted by a sinusoidal function, makes it difficult to discern between a physical or geometric origin of the modulation. Whereas the VHE gamma-ray lightcurve has a modulation similar to the X-ray one (Aharonian et al. 2006), the HE gamma-ray flux is lower at the inferior conjunction (0.45 < φ < 0.9) and higher at superior conjunction (0.9 < φ < 0.45). In addition, at superior conjunction it shows a spectral cutoff at E = 1.9 ± 0.5 GeV, not statistically required for inferior conjunction data (Abdo et al. 2009).
The radial and powerful wind of LS 5039 makes it an ideal system to apply the model presented in this work. The shocked gas cooling regime (see Section 3.1) will be closer to the radiative than to the adiabatic limit, with χ values around unity. Therefore, it seems safe to apply the thermal emission model presented in Section 3.3 to LS 5039. 4.1. Stellar mass-loss rate and pulsar spin-down luminosity in LS 5039 As explained in Section 2, the shape of the CD depends only on the parameter η ∞ , which in turn depends only on the ram pressure of the winds on either side of the CD. Figure 2 shows the shape of the CD in comparison to the orbit of LS 5039 at periastron and apastron for different values of η ∞ . Casares et al. (2005) derived that the O6.5V((f)) star in LS 5039 has a stellar radius of R = 9.3 R , its wind has a terminal velocity of v ∞ = 2400 km s −1 and the orbital inclination is in the range 20 • -70 • . Therefore, the only parameters still unknown that affect the shape of the CD, and thus the thermal X-ray emission, are the stellar mass-loss rate,Ṁ, and the pulsar spin-down luminosity, L sd . Sarty et al. (2011) performed an optical photometric and spectroscopic campaign on LS 5039, and determined the massloss rate of the star to be 3.7 to 4.8 × 10 −7 M yr −1 from Hα THERMAL X-RAY EMISSION FROM THE SHOCKED STELLAR WIND OF PULSAR GAMMA-RAY BINARIES 5 Star to observer P A Figure 2. Sketch of the orbit of the compact object around the companion star in LS 5039. The star (but not the pulsar) is at scale. The two marked positions correspond to periastron (P) and apastron (A). For each of these positions the shape of the CD is shown for η ∞ = 0.0025, 0.025 and 0.25 (blue, red, and green lines, respectively). line measurements. However, Hα line fitting is a diagnostic known to be affected by wind clumping, and owing to its ρ 2 dependence could easily overpredict mass-loss rates by a factor ∼2 (Markova et al. 2004). Bosch-Ramon (2010) studied the effect of wind clumping in the absorption of X-ray emission from LS 5039, but it must be noted that Sarty et al. (2011) did not detect any evidence of dense wind clumps in the analysis of stellar wind emission lines. The lack of variability of the column density measured in X-ray observations along the orbit may also be used to constrain the maximum value of mass-loss rate. Bosch-Ramon et al. (2007) concluded that, considering a non-thermal X-ray emitter inside the binary system, the mass-loss rate could not exceed a few 10 −8 M yr −1 . Adopting a model in which the non-thermal emission is extended, as expected in the PWS scenario, this upper limit is relaxed up to 1.5 × 10 −7 M yr −1 (Szostek & Dubus 2011). In the following calculations we considered the whole range of mass-loss rates discussed here, i.e., from 5×10 −8 M yr −1 to 4.8×10 −7 M yr −1 . For some cases, in which a single representative mass-loss rate is used to test other properties of the thermal emission, we chose the midpoint of the range:Ṁ = 2.65 × 10 −7 M yr −1 .
The spin-down luminosity of the putative pulsar in LS 5039 is not known as it has never been measured directly. A lower limit may be derived from efficiency considerations taking into account the HE gamma-ray emission detected by Fermi/LAT (Abdo et al. 2009). If the HE gamma-ray emission has a magnetospheric origin, the detected flux above 100 MeV of 2.55 × 10 −10 erg cm −2 s −1 implies a spin-down luminosity between 0.2 and 30 × 10 36 erg s −1 , depending on the beaming model, with higher values favored by the high energy spectral cutoff (see Zabalza et al. 2011, for an explanation of the derivation of these values for the case of LS I +61 303). However, the observed orbital modulation of the HE gamma-ray emission, and the lack of pulsations, indicate that its origin might not be magnetospheric but rather IC.
The minimum non-thermal luminosity required for the Fermi emission of LS 5039 can also be estimated in the context of IC radiation. To obtain a conservative estimate, one can assume that electrons are injected only in the relevant energy range, ∼1-100 GeV, their velocities are isotropically distributed, and have time to cool down only through IC. Stellar IC leading to photons in the Fermi band occurs in the Thomson regime. Since this process is strongly anisotropic (L IC ∝ (1 − cos θ IC ) (p+1)/2 with p being the electron power-law index; e.g., Dermer et al. 1992), one needs to account for the IC interaction angle (θ IC ) between the observer line of sight and the seed photon direction. Under these conditions, the injection luminosity of electrons can be written as where the factor 1.7 accounts for the luminosity radiated below 100 MeV, α = (p+1)/2, and θ IC ≤ π/2 for phases around the inferior conjunction. Taking L GeV ≈ 2×10 35 erg s −1 during those orbital phases (Abdo et al. 2009), p ∼ 3 to explain the Fermi spectrum, and a rather conservative θ IC ≈ 60 • , one obtains L e ≈ 2 × 10 36 erg s −1 . However, escape, adiabatic, and perhaps synchrotron losses could easily increase L e by a factor of several. In addition, not all the L sd will go to non-thermal particles 1 , and the injection electron distribution is expected to be broader than ∼1-100 GeV. From all this, L sd 2 × 10 37 erg s −1 seems more realistic. Doppler boosting has not been accounted for but, although it can reduce the required energy budget in certain orbital phases, it will increase it in others. Note that for θ IC ≈ 90 • and phases between inferior and superior conjunction, i.e., with the flow bulk motion unlikely pointing to us, still L e 10 36 erg s −1 , and thus L sd 10 37 erg s −1 . The GeV emitter may be extended, thus relaxing the θ IC dependence, but it should not change significantly our conclusion.

4.2.
Comparison to X-Ray data LS 5039 has been observed at different positions along its orbit by many X-ray observatories (see Zabalza et al. 2008, for a summary). Recently, a very long observation with Suzaku provided an uninterrupted lightcurve of the source during an orbit and a half , and a comparison with previous observations indicated that it was quite stable over periods of up to a decade ). The Suzaku observation provides a very good orbital coverage, but we have chosen, to compare the output of our numerical model, two XMM-Newton observations performed in 2005 (Bosch-Ramon et al. 2007). The superior effective area of XMM-Newton gives us the possibility to extract very good soft X-ray spectra in relatively short observations of duration ∆φ ∼ 0.04. These short observations assure that the ambient conditions do not vary significantly and the measured spectra can be accurately compared to spectra calculated for a given phase. These two ∼15 ks observations were performed at periastron (ObsId 0202950301, orbital phases 0.02-0.05) and apastron (ObsId 0202950201, orbital phases 0.49-0.53) during the same orbital period, thus probing the two extreme separations between the star and the compact object.
For each of the two observations we processed the data through the standard procedure with version 10.0 of the XMM-Newton Science Analysis Software. We filtered out the periods 6 Zabalza, Bosch-Ramon, & Paredes a Unabsorbed X-ray fluxes of the non-thermal power law (P) and the thermal shocked stellar wind (Th) components in the range 0.3-10 keV are given in units of 10 −12 erg cm −2 s −1 . b Interstellar column density is in units of 10 22 cm −2 .
of high particle background for each of the three detectors aboard XMM-Newton, and extracted source and background spectra from the clean event files, as well as the corresponding redistribution matrix and ancillary response files. We binned the extracted source spectra to a minimum of 20 counts per bin, avoiding an oversampling of the detector energy resolution by more than a factor three. Before introducing our model of the thermal X-ray emission of the shocked stellar wind, we performed a fit of the data to an absorbed power law using the spectral analysis package Sherpa bundled with CIAO 4.3. The fit results for this strictly non-thermal source model may be found in the third column of Table 1.
We assigned each of the two groups of spectra (periastron (p) and apastron (a)) the source model described in Sections 3.5: A ISM × (P p,a + Th p,a ). Whereas the normalization and photon index parameters of the power-law component were left free for the spectral fitting, the interstellar absorption and the spindown luminosity of the pulsar were kept constant between the periastron and apastron source models. We note that such a choice for the source model implies the assumption that the non-thermal emission is not significantly affected by circumstellar absorption, whereas it affects the thermal emission and is already computed in the component Th. This would happen for either extended non-thermal emitters or low circumstellar absorption (see Section 6.1 below for a discussion).
As expected from the previous estimate (see Section 3.2), the X-ray luminosity of the thermal component increased monotonically with the value of the spin-down luminosity of the pulsar. In order to derive the spin-down luminosity upper limit, we performed a simultaneous fit of the six spectra (three detectors for each of the two orbital phases). In all cases, the best fit corresponded to the lowest level of emission from the thermal component, i.e., the non-thermal component alone provided the best fit to the data. The upper bound, as explained in Sec. 3.5, corresponds to the upper limit to the pulsar spin-down luminosity.
5. RESULTS 5.1. Thermal X-Ray spectra We used the dynamical, radiative, and absorption models described in Sections. 2 and 3 to calculate the thermal X-ray spectra emitted by the shocked stellar wind of LS 5039. For the properties of this system we find that the cooling regime xx xx xx  (4) is shown as a gray line. The range of pulsar spin-down luminosities excluded by the thermal emission is shown as a hatched region. In all cases, a mass-loss rate of 2.65 × 10 −7 M yr −1 was assumed. of the shocked gas is an intermediate regime, with χ around unity. As an example, forṀ = 2.65 × 10 −7 M yr −1 and L sd = 3×10 36 erg s −1 (η ∞ = 0.025), we obtain radiative cooling and escape timescales at the CD apex during periastron of 9.1× 10 3 s and 7.7 × 10 3 s, respectively, resulting in χ 1.18. The maximum postshock temperature reached is T 0 = 2.43 keV at the apex, with a postshock number density of n 0 2 × 10 10 cm −3 .
In accordance with the estimates of Section 3.2, we found that in general the computed thermal X-ray luminosities increase monotonically with the spin-down luminosity of the pulsar. In Figure 3 we show the 0.3-10 keV X-ray luminosity of the thermal emission from the shocked stellar wind as a function of the pulsar spin-down luminosity, takinġ M = 2.65 × 10 −7 M yr −1 . For scenarios with η ∞ 1 the behavior is generally in good agreement with Equation (4), and the thermal X-ray luminosity can be acceptably explained by a function of the form L X ∝ L α sd with α = 1.381 ± 0.004. This value is significantly lower than that predicted by the estimate of Equation (4). This could be related to the lack of correction for the obliquity of the incoming flow with respect to the shock in the estimate, which could result in an underestimation for lower η ∞ and an overestimation at higher η ∞ , where a larger fraction of the incoming wind impacts obliquely. At higher η ∞ , the proximity of the CD to the stellar surface at periastron provokes a decrease in the velocity of the wind crossing the shock, and therefore a decrease in the total X-ray luminosity. Figure 4 shows the obtained intrinsic and absorbed spectra for different values of the spin-down luminosity from 3 × 10 35 erg s −1 to 3 × 10 37 erg s −1 . The mass-loss rate of 2.65 × 10 −7 M yr −1 assumed for these spectra provides a dense matter field in the binary. Its effect can be seen as a very strong absorption at energies below 1 keV for the periastron spectrum. The column density at this phase is of the order of 10 22 cm −2 , even higher that the measured interstellar one. In contrast, the apastron spectrum for the highest spin-down luminosity shows negligible circumstellar absorption. At apastron, the pulsar would be close to inferior conjunction, and the high opening angle of the CD would mean that most of the thermal X-ray emission only has to travel through the unshocked pulsar wind to reach the observer (see Figure 2), thus being barely absorbed. As discussed in Section 3.3, the hardest X-rays will come from the region around the apex of the CD. For high values of the spin-down luminosity, the apex of the CD at periastron will be close to the stellar surface, and the velocity of the incoming wind will be low owing to the wind β-velocity law. Therefore, the X-ray emission from the apex will be softer, an effect that can be seen in the high energy part of the topmost spectrum of Figure 4. For the parameters used to compute this spectrum, the stagnation point of the CD is at only 0.39 stellar radii from the surface of the star, and the wind velocity at the shock is about one third of the terminal velocity. Note how this effect is not apparent in the corresponding apastron spectra, for which the stagnation point is much farther from the stellar surface.

Spin-down luminosity upper limits
To derive the upper limits of the pulsar spin-down luminosity, we chose to explore the stellar mass-loss rate values mentioned in Section 4.1, corresponding to the estimates of Bosch-Ramon et al. (2007), Szostek &Sarty et al. (2011): 5 × 10 −8 M yr −1 , 1.5 × 10 −7 M yr −1 , and 4.25 × 10 −7 M yr −1 , respectively. The orbital inclination (i = 20 • -70 • ; Casares et al. 2005) will affect the importance of circumstellar absorption at periastron in the emitted spectrum. The result of the pulsar spin-down luminosity upper limit determination, for the whole ranges of mass-loss rate and orbital inclination, is shown in Figure 5. A summary of the results for the extremes and midpoints of the ranges can be seen in Table 2.
We found that pulsar spin-down luminosities above 6 × 10 36 erg s −1 are excluded by the lack of spectral thermal fea-  tures in the observed periastron and apastron spectra, for any combination of mass-loss rate and orbital inclination. This upper limit corresponds to the extreme inclination of 70 • , and for moderate orbital inclinations the upper limit is at all times below 5.5 × 10 36 erg s −1 . Although the wind kinetic luminosity available forṀ = 5 × 10 −8 M yr −1 is lower than for the other cases, the upper limits derived are at least as strict as for higher mass-loss rates. We found that the reduced absorption owing to the lower density of the stellar wind results in brighter emission below 1 keV. In the last column of Table 1, we show the best-fit results of the thermal plus non-thermal source model when fixing the spin-down luminosity at the upper limit we found. The thermal component was computed taking the midpoint values of mass-loss rate and orbital inclination. As compared to the non-thermal source model, the best-fit power law is slightly harder and its unabsorbed flux is ∼5% lower. An illustration of the thermal and non-thermal contributions to the observed spectrum is shown in Figure 6.
6. DISCUSSION 6.1. On the stellar mass-loss rate The stellar mass-loss rate plays a double role in the PWS scenario. On one hand it determines the shape of the CD and the kinetic luminosity that is shocked onto the CD. On the other hand, it determines the importance of photoelectric absorption in the X-ray spectra of the thermal and non-thermal radiation sources. The lack of X-ray absorption variability in LS 5039 along the orbit indicates that, if the non-thermal emitter is inside the orbital system, the stellar wind density must be low. However, a very low value will result in the pulsar wind overpowering the stellar wind for reasonable values of the pulsar spin-down luminosity. The maximum mass-loss rate compatible with a variation of the neutral hydrogen column density along the orbit below 10 21 cm −2 (Bosch- Ramon et al. 2007) is of ∼ 5 × 10 −8 M yr −1 . For such a low mass-loss rate, the pulsar wind could impact on the stellar surface at periastron for reasonable values of spin-down luminosity. However, the lack of orbital variations of the UV spectra of LS 5039 (McSwain et al. 2004) indicate that the stellar surface is not perturbed at any point along the orbit. This places an upper limit on the spin-down luminosity of ∼ 7 × 10 36 erg s −1 . However, if the star has a strong enough surface magnetic field, then magnetic pressure could be enough to balance the pulsar wind away from the surface (Harding & Gaisser 1990). For a surface magnetic field of B = 50 G and assuming a dipole field, the spin-down luminosity could be as high as 3 × 10 37 erg s −1 without the CD collapsing onto the stellar surface. Note, however, thaṫ M = 5 × 10 −8 M yr −1 and B = 50 G imply a wind magnetic to kinetic energy ratio close to unity.
If a significant fraction of the non-thermal emitter is not located deep inside the binary system, the photoelectric absorption of the thermal and non-thermal X-ray emission could then be decoupled (as we have done in Section 3.5), since their locations would be different. For this scenario, the stellar mass-loss rate could be as high as a few times 10 −7 M yr −1 , as measured by Sarty et al. (2011). Additionally, such a location would favor the acceleration of particles up to the very high energies required to account for the VHE gamma-ray spectrum, as well as avoid strong pair creation absorption of VHE gamma rays Bosch-Ramon et al. 2008).
6.2. On the pulsar spin-down luminosity The found upper limits of the pulsar spin-down luminosity pose strong constraints on further modeling of the PWS scenario for non-thermal emission. As discussed in Section 4.1, the pulsar spin-down luminosity should be at least a few times 10 36 erg s −1 in order to account for the GeV emission. The upper limit we have found from the thermal X-ray emission (L sd < 6 × 10 36 erg s −1 ) indicates that the spin-down luminosity must be very close to this value in order for the PWS scenario to be consistent. A detailed modeling of the GeV emission from LS 5039 is required to accurately check the consistency of the observed GeV luminosity and the PWS scenario. Unfortunately, the structure and location of the GeV emitter is not well known, and IC emission is very sensitive to the θ IC -value, which makes difficult to obtain conclusive results. Nevertheless, the lower limit derived here for L e , in the IC scenario, is still robust enough to show that, likely, L e /L sd → 1. It might appear that such a strong constraint can only be achieved at very low η ∞ , but the effect of the orbital motion on the large-scale structure allows for the pulsar wind to be terminated even in the direction opposite to the star up to η ∞ ∼ 1 (Bosch-Ramon & Barkov 2011). We note that the constraint on the efficiency is slightly less restrictive for high orbital inclinations (i 45 • ) and low stellar mass-loss rates (Ṁ 2 × 10 −7 M yr −1 ).
The recent detection of PSR B1259−63 at GeV energies (Tam et al. 2011;Abdo et al. 2011) provides a perfect comparison of the efficiency of spin-down luminosity to non-thermal emission conversion efficiency. During its highest state of emission the total gamma-ray luminosity was comparable to the pulsar spin-down luminosity, which combined with the recent reclassification of the optical star (Negueruela et al. 2011) results in nearly unbearable energy requirements . Even considering Doppler boosting, the efficiency of the shocked pulsar wind in converting the spin-down luminosity into non-thermal HE gamma-ray emission must be quite high, and would only weaken the constraint for certain phases (see Section 4.1).
The acceleration and emission processes responsible for the high efficiency in the conversion from spin-down luminosity to HE gamma-ray emission are currently unknown. IC emission from particles accelerated beyond the light cylinder in the striped pulsar wind is one of the proposed emission processes (Pétri & Dubus 2011). However, when applied to PSR B1259−63, this process can only account for the lowflux state before periastron, and is unable to explain the high luminosity detected after the periastron passage without the invocation of addition seed photon sources apart from the stellar companion. The extremely high efficiency required for LS 5039, similar to the high-flux state of PSR B1259−63, indicates that this process might not be responsible for its HE gamma-ray emission. THERMAL X-RAY EMISSION FROM THE SHOCKED STELLAR WIND OF PULSAR GAMMA-RAY BINARIES 9 6.3. Concluding remarks The semi-analytic model presented here has allowed us to gain valuable insight into the nature of LS 5039. The assumptions described in Section 3 are chosen such that the resulting level of thermal X-ray emission is fairly conservative and thus the upper limit on the spin-down luminosity robust. We have shown that the study of the thermal emission from the shocked stellar wind can be a powerful diagnostic tool for pulsar gamma-ray binaries. Even for cases where no thermal features are observed in the X-ray spectrum, it can be used to place constraints on important properties of the system such as stellar mass-loss rate and pulsar spin-down luminosity. The application of such a model to LS 5039 allowed us to conclude that, if hosting a non-accreting pulsar, this system has a very efficient non-thermal mechanism, with L sd ∼ (3-6)×10 36 erg s −1 .
We acknowledge support by the Spanish Ministerio de