A model for high-energy emission of the Intrabinary shock in pulsar binaries

We present our studies of intrabinary shock emission for astrophysical binary systems with a neutron star. We construct a model for the shock emission and compare the model calculation with the light curve and the spectral energy distribution of the gamma-ray binary 1FGL J1018.6−5856. The model assumes a slow and a fast population of particles accelerated in the shock, and computes the high-energy emission spectra and orbital light curves produced by synchrotron, self-Compton and external Compton processes of the high-energy particles in the shock. The model allows one to study plasma properties and to constrain the binary geometry, most importantly the inclination angle (i). We discuss potential use of this model for other pulsar binaries to determine the inclination angle of the binary hence the mass of the neutron star.


Introduction
In some astrophysical binary objects, it is believed that winds of the two objects interact strongly and form shocks, so called intrabinary shocks (IBSs).These shocks can accelerate particles to very high energies, and the particles in the shock radiate high-frequency electromagnetic waves [1].If one of the stars is a neutron star, the wind of the neutron star is likely composed of cold energetic electrons and positrons.When they are accelerated in the shock, they radiate their energy very efficiently.This radiation can be observed and help to study IBS.
IBS is important for many reasons.Particle acceleration mechanisms in astrophysical shocks are not yet well understood; diffusive shock acceleration (DSA) [2] is a commonly used one but there are others as well.These mechanisms have been studied with particle-in-cell (PIC) simulations [3]; the resulting energy distribution of particles is relatively sensitive to the shock environments which are not well known.Because particles are believed to be accelerated stochastically (e.g., via DSA) in IBS, observations of IBS emission can be used for inferring properties of the shock.IBS emission can also provide us with an independent way of estimating the binary orbital parameters which are important to infer the mass of the neutron star in the binary.In particular, this can supplement the current method of inferring the neutron star mass with optical spectroscopic and photometric observations and modeling.
Formation of shocks in pulsar binaries is similar to that in pulsar wind nebulae [4,5]: by the interaction of pulsar wind with the ambient medium (stellar wind in pulsar binaries and interstellar gas or supernova remnants in pulsar wind nebulae).A contact discontinuity (CD) is formed when the wind of the companion and the relativistic pulsar wind interact in the orbit.Around the CD, shocks are � e-mail: hjan@chungbuk.ac.kr formed of which the relativistic pulsar wind shock is of interest in regard to the high-energy emission in these systems [6].In the pulsar wind shock, particles are stochastically accelerated and are believed to be bulk accelerated as they flow [7,8]; Lorentz beaming is expected which may be observed as spikes in the orbital light curve [9,10].
Using the IBS scenario, we develop a phenomenological model for the high-energy emission [10] in pulsar binaries.Below we show details of the model components and present the results of application to the gamma-ray binary 1FGL J1018.6−5856.We then discuss plasma properties in the IBS of the source and potential use of the model for other pulsar binaries.Note that similar models have been used for other sources such as the gamma-ray binary LS 5039 [8] and the pulsar binary PSR J2215+5135 [9].

The IBS emission Model
Emission in an IBS is very complex.There are relativistic electrons and positrons stochastically accelerated in the shock, hence synchrotron radiation is expected.Then the particles will inverse-Compton upscatter the synchrotron photons (synchrotron-self Compton, SSC) and the stellar optical/UV photons external to the shock (external Compton, EC); these emissions will be in the ∼TeV band [6,10].Due to these emissions, the energy distribution of the particles evolves in time and space, so does the emitted photon spectrum.Furthermore, the IBS emission will be Lorentz boosted because some of the particles are bulk accelerated in the flow as shown in magnetohydrodynamic simulations [7,8], and so the Lorentz beaming also needs to be considered.Of course, the stellar and the pulsar emission can contaminate the IBS emission.Therefore, studying IBS requires detailed numerical simulations and comparison of simulations and observations.However, these simulations are yet to be improved for the relativistic pulsar  wind, so instead, we use a relatively simple phenomenological model for IBS emission in order to infer the flow properties and the orbital geometry using high-energy observations [10].

The Shock Geometry
As noted above, IBS emission is highly anisotropic, so properly estimating the shape of an IBS is very important.The geometry of an IBS is determined by the shape of the CD.The CD produced due to isotropic wind interaction is in principle paraboloid like, and the shape can be calculated [11]: r(θ p ) = Dsinθ s /sin(θ p + θ s ) and θ s cotθ s = 1 + β(θ p cotθ p − 1), where r(θ p ) is the distance from the compact object to a point in the CD, D is the orbital separation, θ s and θ p are angles with respect to the line of centers from the companion star and the compact object, respectively, and β is the wind momentum flux ratio (Figure 1).Once β and D are given, the full locus of the CD can be calculated.
We note that the above calculation is for isotropic winds.Realistically, the wind of a neutron star (and the stellar companion) may not be isotropic, then the CD shape may change.However, isotropic winds may be a good approximation, and the model with the approximation can capture most of the important features of the emission as we show below.
Of course, the shock flow can extend to a large distance and may make other structures [8].We note that the IBS emission becomes rapidly weaker with distance s due to gradually reducing B and seed photon density (for SSC and EC).So we calculate the IBS emission out to 3× orbital separation (s ≤ 3D).Emission at larger s may explain the low-energy radio emission of pulsar binaries.

High-energy emission in IBS
In an IBS, particles are stochastically accelerated, perhaps via DSA.Of course, magnetic reconnection may play a role.The energy distribution of the particles can be diverse depending on the shock environment, but generally a power law with an index of p = 2 is expected (i.e., dN/dγ e ∝ γ −p e ).Some of the particles are further bulk accelerated in the IBS flow to have a Lorentz factor Γ > 1, hence there can be multiple populations with different Γ [7,8].In hydrodynamic simulations, it is shown that Γ can be modest, but the exact value and the distribution of Γ for relativistic pulsar winds are not known.
As we mentioned earlier, in the IBS scenario e − /e + plasma in the shock makes synchrotron radiation under the influence of B field supplied by the pulsar.Although it is impossible to measure B in the shock directly, we can scale the neutron star's field: the ∼ 10 12 G field of the neutron star will decrease following the dipole formula (∝ r −3 ) out to the light cylinder radius (R LC = cP/2π, where c is the speed of light and P is the rotation period of the neutron star), and then the stripped pulsar wind will carry the field (∝ r −1 ) to the shock at r(θ p ) [6].So once the timing solution (spin period P and the derivative Ṗ) of the pulsar is obtained, dipolar B at the neutron star's pole can be inferred with B = 3.2 × 10 19  √ P Ṗ G, and hence B in the shock can be estimated.Note that this estimation is rough and serves only as a guide.The synchrotron loss is dE/dt ∝ −γ 2 e , and the emitted photon energy distribution follows a power law with an index (p − 1)/2 (i.e., F ν ∝ ν −(p−1)/2 ) and lies in the X-ray to the low-energy gamma-ray bands for typical values of shock parameters.
The dominant photon populations in the shock are synchrotron photons emitted by the particles themselves; there are also stellar blackbody photons supplied by the stellar companion.All these photons interact with the particles in the shock via inverse-Compton scattering (i.e., SSC and EC).In these processes, the synchrotron photons gain energy by a factor of ∼ γ 2 e by SSC; when observed, these are further Lorentz-boosted by δ , where θ V is the viewing angle of the observer (Figure 1).The stellar photons, however, are external to the shock and so there is another ∼ Γ 2 energy gain (with the scattering angle dependence), so the gain factor is ∼ Γ 2 γ 2 e for EC.These SSC and EC photons become high-energy gamma-ray photons (>10 GeV) to be observed.
The observed spectral energy distribution (SED) of the IBS-emitted (+reprocessed) photons is characterized with a "double-hump" continuum structure, the low-energy hump being produced by the synchrotron process and the high-energy hump by SSC and EC.Of course, there are contaminations from the stellar blackbody photons at ∼ eV and pulsar magnetospheric emission at ∼ GeV.

The IBS model
The model we develop [10] is basically a combination of the above.Using the masses of the two stars and the orbital EPJ Web of Conferences 168, 04013 (2018) https://doi.org/10.1051/epjconf/201816804013Joint International Conference of ICGAC-XIII and IK-15 on Gravitation, Astrophysics and Cosmology period measurement, we can determine the orbital separation (D) using the Kepler's law.Of course, orbital eccentricity (e) needs to be assumed, and later can be verified because change in the separation due to e will cause the magnetic field strength to vary over the orbit.This change may be visible in the synchrotron light curve (orbital variation).
At each orbital phase, distance to a point on the CD (r(θ p )) is calculated using the formula in Section 2.1 for given β.We then assume the shock-accelerated electron distribution (power law), magnetic field strength B, and the bulk Lorentz factor (Γ = 1) at the apex of the shock.As the particles flow along the shock to a distance s(θ p ) from the apex (see Figure 1), these quantities vary: N e [s(θ p )], B[s(θ p )], and Γ[s(θ p )].Although the exact functional forms need to be found using magnetohydrodynamic simulations, simple analytic forms can be used for our phenomenological model.For example, N e ∝ 1−cosθ p can be used because pulsar wind particles are injected over the whole shock and at θ p particles injected from θ p =0 to θ p are accumulated, and B ∝ r(θ p ) −1 since r(θ p ) >> R LC for typical millisecond pulsars.The functional form for Γ may be more complex [7,8], but a simple linear trend (Γ(s) ∝ s) may be a good starting point.
At a given point s 0 on the shock, there are multiple particle populations injected: a fresh injected one from the pulsar and another that flows from the apex (s < s 0 ).The slope (power-law index) of the distribution is the same between the populations, but those from the apex lose more energy as they already emitted synchrotron/SSC/EC photons (cooling) while they flow to s 0 from the apex at s = 0.Because the energy loss is proportional to ∼ γ 2 e for the synchrotron and the Compton processes, the highest-energy particles lose energy quickly and become low-energy particles.Hence, the flow from the apex lacks high-energy particles.This cooling results in a break in the particle distribution (cooling break at γ e,b ), so the electron energy distribution will be a broken power law: dN/dγ e ∝ γ −p 1 e for γ e ≤ γ e,b and dN/dγ e ∝ γ −p 2 e for γ e > γ e,b .This will be visible as a break in the observed photon spectrum as well.In addition, the minimum and the maximum electron energies need to be specified.Because we are concerned with high-energy emission, we use a reasonably small γ e for the minimum such that the emitted photon energy is well below the X-ray band.The maximum γ e depends on the acceleration process, but in DSA the maximum can be as high as ∼ 100 MeV, the radiation reaction limit.This results in a cut off in the photon spectrum at the highest frequency.
With all these prescriptions, we compute the synchrotron and the inverse-Compton radiation at each point of the shock and Lorentz-boost them by δ (for synchrotron and SSC) to get the observed spectrum (see [12] for example).We repeat this for each orbital phase of interest, calculate fluxes, and construct the model light curves.Also, we construct the phase-summed model SED by summing up the spectra over the orbital phases.
3 Application to the gamma-ray binary 1FGL J1018.6−58561FGL J1018.6−5856 is identified as a gamma-ray binary based on the gamma-ray modulation [13].The orbital period is measured to be 16.544 days [14].Importantly, phase 0 in Figure 2 is identified as the inferior conjunction of the compact object [15].Because pulsations are not detected from this source yet, it is not clear whether the compact object is a neutron star or a black hole, but previous studies favor a neutron star over a black hole [15][16][17][18], so we assume a neutron star in this work.
Figure 2 shows the broadband SED (left) and the highenergy orbital light curves (middle and right) of the source [10].In the SED, the lowest-frequency radio photons (ν ∼ 10 10 Hz) are assumed to be emitted in a different region from that for the high-energy photons, so we do not attempt to explain these points.In the ∼ 10 15 Hz band, stellar blackbody emission dominates, hence these data allow an accurate determination of the seeds for EC.The Xray to low-energy gamma-ray (ν ∼ 10 17 −10 22 Hz) photons are produced by the synchrotron process, hence provide us with information on the particle distribution.The slope of the observed SED in this band in Figure 2 left is ∼ −0.5, which implies p ≈ 2 for the electron distribution, and the maximum energy of the synchrotron emitted photons is ∼1 GeV which implies that δ ≈ 10 (hν max = δ×100 MeV).Although the spectral break point is not detected by any observatory (no observatory operating in that band), the model expects it to be in the MeV range to match both the X-ray and the ∼100 MeV spectra; then B is inferred to be ∼a few Gauss.For this magnetic field, the maximum energy of electrons required to emit synchrotron photons at the radiation reaction limit is estimated to be γ e ≈ 10 8 (see Table 1).
At ∼ 1 GeV, the emission is mainly from the pulsar magnetosphere and contaminates the IBS emission (Figure 2 right).So we model the pulsar emission with a power-law with an exponential cut-off spectrum to separate out the IBS emission (see [10] for the pulsar spectral parameters).SSC and EC photons lie in the ν ≥ 10 26 Hz band; these are also sensitive to the properties of the flow, the electron distribution and Γ.The flow properties we inferred above explain the SSC+EC SED well (see Figure 2  left).
Moving on to the light curve modeling, the X-ray and low-energy gamma-ray light curves are explained by the beamed emission aligned with the observer's line of sight (LoS) for the spike at phase 0 and isotropic emission (broad hump peaking at phase ∼0.4) modulating orbitally due to the change in separation between the stars (and so B) with an addition of the orbitally constant pulsar magnetosphere emission in the gamma-ray light curve (Figure 2 right).The models are shown in lines and the parameters are presented in Table 1 (see [10] for more details).Because the model is very complicated and does not allow us to determine a unique set of parameters, we do not formally fit the data and so do not report the uncertainties for the parameter values.Broadband SED (left), the X-ray (middle) and the gamma-ray (right) light curves of 1FGL J1018.6−5856.The X-ray data are taken with Swift, XMM-Newton and NuSTAR observatories, and the gamma-ray data are taken with the Fermi gamma-ray space telescope.The IBS models for this source described in Section 3 are shown in lines.Various emission components are denoted in different colors in the left figure.In the low-energy gamma-ray light curve (right), we also plotted the pulsar emission (green dotted) to isolate the IBS emission (red).The figures are taken from [10].a These vary along the orbit and the shock distance s.

Discussion
We developed an IBS emission model for binary systems in which winds of the two stars interact, producing IBS.
In particular, we are concerned with high-energy (>Xray) emission in astrophysical binary systems which have a relativistic wind and so can produce a relativistic wind shock.The characteristics of the model are the doublehump structure in the SED, and spikes (for proper inclinations) and broad hump (for eccentric orbits) in the highenergy light curve.We applied the model to the gamma-ray binary 1FGL J1018.6−5856 because it has features that can be easily explained by our IBS model [10].Although the model parameters are not uniquely determined, there are several things we learned about the shock properties and the orbit.(1) the electron distribution in the shock is very similar to that expected by diffusive shock acceleration theories, p ≈ 2, (2) in the shock, electrons are maximally accelerated to the radiation reaction limit of DSA, (3) the observed data require two populations of electrons, a slow Γ = 1 population and a fast bulk-accelerated population Γ ≤ 7, the latter being only a small fraction (a few %) of the former, (4) our LoS is aligned well with the shock tail (asymptotic tangent).
(1) and ( 2) suggest that particles can be accelerated to very high energies (to tens of TeV) in the IBS which may have interesting implications for the cosmic-ray studies.
(3) qualitatively agrees with the limited numerical magnetohydrodynamic simulations [7,8].Although these simulations were not done using a realistic pulsar wind with bulk Lorentz factor ≥ 10 4 , they seem to predict some features of the flow: multiple populations with a smaller particle fraction for higher Γ.In our model, we used two populations, but the fast population has multiple values of Γ (linearly growing with s in our simple model) so in principle, we have multiple populations.In addition, our modeling requires that the fast population is only a small fraction of the slow population, which is qualitatively similar to the prediction of the numerical studies [8].We find that the result is not very sensitive to the choice of the function Γ(s), hence more realistic Γ distribution in an IBS flow obtained with the simulations may work; this can be done in the future when both simulation and the IBS model are further improved.
Finally, (4) has very interesting implications for studying pulsar binaries because it can give us a new estimation of the inclination of the systems [9,19].Even without detailed modeling, (4) can be concluded based on the amplitude and the width of the spike.This can be done in the following way.In the IBS model, the spike in the light curve is produced because of beaming which is maximum when our LoS is aligned with the accelerated flow, i.e., the shock tail.Suppose an observer is looking at the binary near the normal to the orbit.Then the shock tail will not cross the LoS anywhere in the orbit; the observer sees only the side of the shock.In this case, beaming (strongly dependent on the viewing angle θ V ) will be very weak and the observer will not see any spike in the light curve.On the other hand, if an observer is near the orbital plane, the observer will see two edges of the shock as the binary orbits and hence will detect two spikes.In between, if an observer is near the tail of the shock, the observer will see only one spike.So simply the number of spikes in the light curve can already tell us much about the inclination, and the light curve modeling can give us more detailed infor-EPJ Web of Conferences 168, 04013 (2018) https://doi.org/10.1051/epjconf/201816804013Joint International Conference of ICGAC-XIII and IK-15 on Gravitation, Astrophysics and Cosmology mation.Of course, the situation may be more complex if the orbital speed of the stars is comparable to the wind speed [9].
Determination of the orbital inclination is particularly important for inferring the mass of the neutron star in pulsar binaries.The mass of the neutron star in a pulsar binary can be estimated with the mass function if the mass ratio of the two stars and the orbital inclination are determined.The former can be determined by the radial velocity measurement, and the latter by the optical light curve modeling (see [9,20] for example).The mass determined in this way is subject to relatively large uncertainties (to the third power in sini) due to the correction for the mismatch between the centers of light and mass, and an independent determination will be very helpful.Of course, for relativistic pulsar binaries measuring the relativistic delay of pulses can be used for accurate determination of the neutron star mass (see [21] for example).
As we show above, the IBS emission can tell us about the orbital inclination of pulsar binaries, especially when there are spikes in the high-energy light curve [9,19].In particular, the orbital inclination can be determined with the synchrotron light curve modeling only which can run reasonably quickly [9] compared to the one with the full SSC+EC emission models.Of course, there is a weakness in this method; measuring the crucial parameter β, the wind momentum flux ratio, for determining the inclination is hard.The way the shock bends can tell us about the relative strength of the two winds, but accurate determination of β (and so i) requires measurements of the momentum flux of both winds.Furthermore, the pulsar wind can be anisotropic, so the angular profile of the pulsar wind momentum needs to be considered.Also note that for relativistic shocks [7] the opening angle may differ from the analytic derivation of [11], hence more realistic magnetohydrodynamic simulations need to be performed.These will benefit from the high-energy modeling presented here for inferring the shock properties which can be compared with the simulation results, helping to improve the simulations further. The

Figure 1 .
Figure1.Vertical cross section of the binary system and the IBS.Two stars (pink) and the CD (blue) are shown.The CD should look like a paraboloid in 3D.This plot is generated for the mass ratio q = 17 and the wind momentum flux ratio β = 25.Various orbital parameters are denoted.The figure is taken from[10].

Figure 2 .
Figure2.Broadband SED (left), the X-ray (middle) and the gamma-ray (right) light curves of 1FGL J1018.6−5856.The X-ray data are taken with Swift, XMM-Newton and NuSTAR observatories, and the gamma-ray data are taken with the Fermi gamma-ray space telescope.The IBS models for this source described in Section 3 are shown in lines.Various emission components are denoted in different colors in the left figure.In the low-energy gamma-ray light curve (right), we also plotted the pulsar emission (green dotted) to isolate the IBS emission (red).The figures are taken from[10].
Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden).Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged.This work performed in part under DOE Contract DE-AC02-76SF00515.This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2017R1C1B2004566).

Table 1 .
The orbital parameters and the shock properties obtained with the modeling