A hydrodynamics-informed, radiation model for HESS J0632$+$057 from radio to gamma rays

Relativistic hydrodynamical simulations of the eccentric gamma-ray binary HESS J0632$+$057, show that the energy of a putative pulsar wind should accumulate in the binary surroundings between periastron and apastron, being released by fast advection close to apastron. To assess whether this could lead to a maximum of the non-thermal emission before apastron, we derive simple prescriptions for the non-thermal energy content, the radiation efficiency, and the impact of energy losses on non-thermal particles, in the simulated hydrodynamical flow. These prescriptions are used to estimate the non-thermal emission in radio, X-rays, GeV, and TeV, from the shocked pulsar wind in a binary system simulated using a simplified 3-dimensional scheme for several orbital cycles. Lightcurves at different wavelengths are derived, together with synthetic radio images for different orbital phases. The dominant peak in the computed lightcurves is broad and appears close to, but before, apastron. This peak is followed by a quasi-plateau shape, and a minor peak only in gamma rays right after periastron. The radio maps show ejection of radio blobs before apastron in the periastron-apastron direction. The results show that a scenario with a highly eccentric high-mass binary hosting a young pulsar can explain the general phenomenology of HESS J0632$+$057: despite its simplicity, the adopted approach yields predictions that are robust at a semi-quantitative level and consistent with multiwavelength observations.


INTRODUCTION
HESS J0632+057 is a moderately bright, highly eccentric gammaray binary hosting a Be star and a compact object of unknown nature, with an orbital period slightly below one year (Casares et al. 2012). Discovered for the first time as a point-like TeV source of unclear counterpart (Aharonian et al. 2007), HESS J0632+057 was latter detected in radio, X-rays, GeV and TeV energies (e.g. Hinton et al. 2009;Bongiorno et al. 2011;Moldón et al. 2011;Aliu et al. 2014;Malyshev & Chernyakova 2016;Li et al. 2017). The non-thermal emission from HESS J0632+057 presents a general tendency to peak before apastron, with a broad secondary peak after apastron in X-rays and TeV (e.g. Aliu et al. 2014). The GeV lightcurve is consistent with this behaviour (Li et al. 2017), and interferometric observations also show a displacement of the radio emission in the sky before apastron (Moldón et al. 2011).
The fact that the only confirmed gamma-ray (high-mass) binary Correspondence author: mbarkov@purdue.edu (MVB) hosting a pulsar is another eccentric Be binary, PSR B1259-63 1 , and the similar spectral energy distribution between HESS J0632+057, and PSR B1259-63 (the latter when not far from periastron), have been considered indications that HESS J0632+057 may be also hosting a non-accreting pulsar (see Bosch-Ramon et al. 2017, and references therein). If this were the case, the pulsar and the stellar wind would collide, forming a complex interaction structure that would be behind the gamma-ray emission from HESS J0632+057, as it is thought to be the case in PSR B1259-63 (Aharonian et al. 2005). The possibility that HESS J0632+057 is hosting a pulsar led Bosch-Ramon et al. (2017) to consider, for this source, a scenario similar to that explored by Barkov & Bosch-Ramon (2016) for PSR B1259-63, in which an unstable spiral structure, made of shocked flow, forms due to the star and the pulsar wind collision and the effect of orbital motion. The authors found that, as in PSR B1259-63, the pulsar wind energy accumulated between periastron and apastron passage due to spiral formation. Before apastron, this energy was quickly released, as the spiral was disrupted by the impact of the pulsar wind in the periastron-apastron direction. The kinematics of the disrupted material would be consistent with X-ray observations and the system geometry in the sky (Pavlov et al. 2015;Miller-Jones et al. 2018). This result led Bosch-Ramon et al. (2017) to propose that this behaviour, directly related to the high eccentricity of the system (e ≈ 0.83; see Casares et al. 2012), may explain the multiwavelength observational properties of HESS J0632+057.
In the present work, we study in detail the behaviour of the radiation produced in the shocked pulsar wind along the orbit. For that purpose, we adopt simple but physically consistent prescriptions for the non-thermal energy content, the radiation efficiency, and the impact of energy losses on the non-thermal particles of the flow. These prescriptions, plus the hydrodynamical information obtained by Bosch-Ramon et al. (2017), allow us to estimate the amount of radiation produced in the shocked pulsar wind in different spectral energies and orbital phases. An alternative scenario, partially based on accretion and focused on the GeV emission of HESS J0632+057, was presented by Yi & Cheng (2017).
The paper is distributed as follows: In Sect. 2, we briefly present the main features of the hydrodynamical study from Bosch-Ramon et al. (2017). Then, in Sect. 3, we describe the prescriptions adopted, and present the calculation results. Finally, the results are discussed in Sect. 4.

HYDRODYNAMICAL APPROACH
The computational scheme adopted in Barkov & Bosch-Ramon (2016); Bosch-Ramon et al. (2017) considers the problem at a distance from the binary much larger than the system semi-major axis (a). This allows spherical coordinates, centred at the binary, to reasonably describe the wind interaction on those scales. On these scales and geometry, the winds are well reproduced as a conical, radial and relativistic wind (pulsar), surrounded by a radial slow wind (star), both supersonic, and centred at the binary location. Along the orbit, the pulsar wind cone points in the changing star-pulsar direction. In this way, the cone rotates with time reproducing the orbital motion of period T = 321 days 2 with e = 0.83. If one takes the masses of the Be star and the pulsar equal to 20M and 1.44M , respectively, one obtains a ≈ 3.5 × 10 13 cm. On scales larger than the binary, the simulation leads to an interaction structure very similar to that found in a 3D simulation when the colliding-wind region within the binary is included (Bosch-Ramon et al. 2015). As noted in Barkov & Bosch-Ramon (2016); Bosch-Ramon et al. (2017), the mass of the relevant region of the disc of the Be star is likely much smaller than that of the stellar spherical wind released during one orbit, and thus the dynamics of the former is neglected on the scales of interest.
Given the large range of spatial scales to cover for eccentric binaries, a fully 3-dimensional (3D) simulation is still too costly even if the binary region is simplified in the way explained. One can however assume that the mass, momentum and energy hydrodynamical fluxes get diluted as 1/r 2 (r: distance from the binary), and one can then focus on the equatorial plane region of the simulated system (i.e. a very low resolution in θ). This formally yields a correct description for that region if most of the material there is supersonic, which is a good approximation at the onset of the spiral and in the apastron side of the simulated region, but not everywhere. Nevertheless, the similarity between the results of (Bosch-Ramon et al. 2015) and those of Barkov & Bosch-Ramon (2016); Bosch-Ramon et al. (2017) indicates that the major features of the hydrodynamics at large scales are already captured by the simplified 3D scheme, which also shortens the computing time by several orders of magnitude. Figures 1 and 2 show density maps at periastron (orbital phase 0; right of the grid) and apastron (orbital phase 0.5; left of the grid) for HESS J0632+057 computed following the approach just presented (see also Bosch-Ramon et al. 2017). To obtain these maps, two orbits were simulated to allow the interacting flows to reach a quasi-steady 3 state. There is a shock related to Coriolis forces at the onset of the spiral. This shock is very close to the binary around periastron, where the shocked pulsar wind bends, and is barely visible at the very centre of Fig. 1. One sees in Fig. 2, slightly before apastron, that the scale of this shock is much larger, and the pulsar wind is also disrupting the spiral arm in the periastronapastron direction. gives the amount of pulsar wind energy accumulated along the orbit within 100 a (roughly the fist spiral turn), whose shape is explained by the growth of the shocked flow, one-arm spiral after periastron, and by the sudden disruption of the spiral arm in the periastron-apastron direction around apastron. As described in the next section, Fig. 3 is in fact giving us information on the non-thermal particles present in the flow assuming that they contain a fixed fraction of the total internal energy of the system. Note that the spiral arm evolution along the orbit is characterized by the shocked dense stellar wind, which is slow enough to consider that the shocked fast pulsar wind is stationary at a given orbital phase.
If non-thermal particles are mostly accelerated after the pulsar wind is affected by the spiral turn, the relatively wide nature of the binary, and the shock location, allow these particles to flow adiabatically with the shocked pulsar wind but around periastron, when this shock is the closest to the binary, and radiation losses are the strongest. Therefore, Fig. 3 can be giving us an accurate description of the non-thermal particle energetics for most of the orbit, and the flow internal energy evolution can inform calculations of the non-thermal emission from the shocked pulsar wind structure. Around periastron, some corrections should be introduced due to fast radiation cooling.

RADIATION
The hydrodynamical approach presented above can provide us with semi-quantitative clues of the non-thermal behaviour. The soundness of these clues relies on the fact that the equatorial region is energetically dominant, on the assumption that particles are (mostly) accelerated at the Coriolis shock, and on the validity of the adiabatic approximation for the non-thermal energy evolution for most of the orbit. When radiation losses dominate over adiabatic losses near periastron, some correction has to be applied. We remark that there are important uncertainties that would render a more accurate approach still very approximate, like inaccurate orbital parameters, a not well   constrained physics of stellar and the pulsar wind physics, unknown dominant acceleration mechanism and location, and unclear magnetic field conditions. In this context, simple analytical prescriptions for the non-thermal processes should suffice at this stage to provide robust trends of the non-thermal emission behaviour.
Here, it is assumed that non-thermal particles are present in the shocked pulsar wind only beyond the point at which the Coriolis shock bends the flow. Beyond that distance, the non-thermal energy density of the shocked pulsar wind in a given cell, u cell,NT , is first computed in the laboratory frame in the adiabatic approximation: where 3P is the (relativistic) internal energy density, Γ the bulk Lorentz factor, χ pw the mass fraction of pulsar wind in the cell, and η NT an unknown non-thermal particle fraction parameter. Other regions, like the wind-colliding region within the binary, may contribute to the lightcurve, but we focus here on the scales of the simulation.
We assume that the non-thermal particles are electrons(/positrons) in the shocked pulsar wind, and focus on synchrotron and inverse Compton (IC) as the relevant emission processes at specific particle energies, i.e. those yielding emission at a particular frequency. The corresponding luminosities can be estimated by taking the non-thermal energy per cell, and divide it by the corresponding cooling timescale of a particular radiation process at the energy of interest.
For electrons interacting through IC with stellar photons of ∼ 3 kT ∼ 10 eV (synchrotron photons are negligible as IC targets), the following expression for t IC = γ/ γ IC is used (adapted from Bosch-Ramon & Khangulyan 2009; see also Khangulyan et al. 2014): γ IC = 5.5 × 10 17 T 3 mcc γ log 10 (1 + 0.55γT mcc )× where γ is the electron Lorentz factor, T mcc = kT * /m e c 2 the stellar temperature in electron rest energy units, R * the stellar radius, and r the distance to the star, which yields the typical size of the emitting region. This radiation timescale is valid assuming that the target photon field is isotropic or, loosely, the angle between the observer and the target direction is not small. As the system inclination is not well known, and the orbital elements may not be fully settled, at this stage we neglect angular IC effects. For synchrotron radiation, the adopted timescale is: where B is the magnetic field in G, assumed to be disordered for simplicity, and The value of B in the shocked pulsar wind is taken as B = √ B 8πP, where B ≤ 1 is a non-dimensional normalization parameter. Values B ∼ 1 would render the hydrodynamics assumption invalid.
The IC and synchrotron lightcurves can be calculated integrating the emissivity of each cell, u cell,NT /t IC/sync , over the emitting volume: L IC/sync ∼ ∫ V (u cell,NT /t IC/sync )dV, where dV is the cell volume in the laboratory frame. This is approximately correct for an electron injection distribution in energy ∝ E −2 , which is common for non-thermal astrophysical sources. Much harder (softer) energy distributions would enhance the synchrotron and the IC fluxes at the highest (lowest) energy bands. We note that the electron energy distribution inferred by Malyshev et al. (2017) through model fitting of multiwavelength data goes from ∝ E −1.3 at low energies to ∝ E −2.3 at high energies, although the observed X-ray photon index presents significant changes along the orbit around ∼ 1.5 (see fig. 3 in that work), i.e. ∝ E −2 . An accurate treatment of the electron energy distribution is important, but it is left for future studies. Given that the bulk velocity of the shocked pulsar wind is mildly relativistic at most, the semi-quantitative nature of our approach, and the uncertainties in the binary and the binary-observer geometry, we neglect at this stage Doppler boosting effects.
The expressions given for L IC/sync are valid under dominant adiabatic losses. For a certain region of size ∼ r, one can estimate the non-radiative losses (adiabatic losses and escape) as where υ f is the typical flow dynamical velocity, which is with υ r being the radial velocity, and υ t ≈ υ 2 r + υ 2 θ + υ 2 φ . This prescription is adopted to account for expansion of the spiral arms even in regions in which v r is small or negative.
If t nonrad > t rad , where t rad = 1/(1/t IC + 1/t sync ), the adiabatic approximation breaks. In the cells where it happens, u NT is depleted by radiation losses, and only a fraction f ∼ min(1, t rad /t nonrad ) of the non-thermal energy would be present. To account for this, f is incorporated to the luminosity calculation as: an additional factor a Band ∼ 0.1 is also included to account for the narrow range in energy considered of the broadband radiation. We assume that the non-thermal particles are injected at twice the turnover radius x (Bosch-Ramon & Barkov 2011;Bosch-Ramon et al. 2012), roughly the location at which the shocked pulsar wind strongly bends, where x = (3/2) υ w η 1/2 PW /Ω,, with υ w being the stellar wind velocity, η PW = 0.1 the pulsar-to-stellar wind thrust ratio, and Ω the orbital angular velocity (see Bosch-Ramon et al. 2015). Thus, we integrate the volume in Eq. (7) from 2 x up to the outer computational boundary. Note that even when particles of high enough energy quickly cool down close to the binary around periastron, f does not reduce the emission from farther distances, as it should if particles are accelerated around 2 x, and therefore our calculation overestimates L IC/sync around periastron.
In this work we focus in the energy bands around 5 GHz and 10 keV for the synchrotron emission, and 1 GeV and 1 TeV for IC. The lightcurves of three contiguous simulated orbits are presented in Fig. 4. The three lightcurves are similar, which means that they are not much sensitive to the initial conditions of the first orbit, nor to the orbit-by-orbit variations of the spiral structure. Figure 5 shows the radio maps at 5 GHz for three different orbital phases: 0.3, 0.41, and 0.52, and two different Gaussian filters that mimic a radio interferometric beam: σ = 1 a (∼ 2 mas) and σ = 50 a (∼ 80 mas). In both Figs 4 and 5, B and η NT have been fixed to 0.5 and 0.1, respectively, keeping them constant with time and in the whole computational domain.

SUMMARY AND DISCUSSION
Using hydrodynamical information, we have computed the nonthermal emission in radio, X-rays, and gamma rays for the gammaray binary HESS J0632+057, in the pulsar scenario, and focusing on particles accelerated at the spiral onset. Remarkably, the multiwavelength predicted lightcurves and fluxes, and even morphology (radio), are in rather good agreement with observations (e.g. Moldón et al. 2011;Aliu et al. 2014;Malyshev et al. 2017;Maier & the VERITAS Collaboration 2017), in particular taking into account the approximate nature of our radiation prescriptions.
All the presented lightcurves show a peak around phases 0.3-0.4, which is connected with the accumulation of the hot shocked plasma in the inner spiral arm, and its later quick release as the spiral arm is disrupted in the periastron-apastron direction. In the IC lightcurves (GeV and TeV) a narrow peak is also present right after periastron passage, possibly because of high IC efficiency at r 2 x. However, the IC peak may be an artifact of our simple way of including radiation cooling, but we cannot completely discard its presence without more realistic radiation calculations. Nevertheless, for the way how radiation cooling is treated here (i.e. the f factor in Eq. 7), we are quite confident that the emission around periastron is overestimated rather than underestimated. It is worth noting however that a weak hint of a small TeV peak around phases ∼ 0.0 − 0.1 is seen in Maier & the VERITAS Collaboration (2017), when folding in phase MAGIC data from Aleksić et al. (2012) with a newer orbital period. The computed lightcurves do not reproduce the broad secondary peak in the X-ray and the TeV lightcurves (e.g. Aliu et al. 2014). We recall nevertheless that we do not consider the physics at the scales of the binary, on which the Be disc, among other factors, could play a relevant role. In addition, our simplified approach to hydrodynamics, and its coupling with radiation, may be hiding effects that would arise in more complicated calculations.
The radio maps unveil an extended structure moving in the periastron-apastron direction. The radio flux evolution, and the displacement of the radio structure on the studied spatial scales, are consistent with those found by Moldón et al. (2011). It is predicted that the motion of the observed radio structure should be aligned with the periastron-apastron direction. Moreover, a weak ellipsoidal extension is hinted also in the apastron-periastron direction, formed by the half ring of shocked pulsar wind present at the periastron side seen in Figs. 1 and 2. Therefore, X-ray (Barkov & Bosch-Ramon 2016) and radio (see also Bosch-Ramon et al. 2017) moving structures may appear prominently at the apastron side, with a fainter counterpart at the periastron side.
The scenario presented here would not apply if the binary were less eccentric, as proposed that HESS J0632+057 may be by Moritani et al. (2018), and discussed also in Malyshev et al. (2017) in the context of X-ray data. In a less eccentric binary (our calculations point at e 0.75), the first spiral arm at the apastron side would survive apastron, and a strong lightcurve peak would not be expected before that orbital phase. In principle, the qualitative agreement of our predictions with observations may suggest that the system is indeed very eccentric and hosts a young pulsar, although the uncertainty in the orbital elements has to be reduced to confirm the high eccentricity. Such a high eccentricity would, in the scenario presented here, strengthen the case for a pulsar in the system. Malyshev et al. (2017) and Moritani et al. (2018) considered the possibility of an inclined Be disc playing an important role in the high-energy emission. If the Be disc mass beyond the pulsar location were not much smaller than that injected by the stellar wind during T, then we would tentatively expect an even stronger accumulation of the pulsar wind energy before apastron because of the disc mass, disrupting the spiral faster. This could make the lightcurve peak before apastron somewhat higher, narrower, and earlier. Nevertheless, the role of a rather massive disc, and its inclination, can only be properly assessed through numerical calculations. Interestingly, Moritani et al. (2018) and Malyshev et al. (2017) also proposed that the periastron and apastron orbital phases are significantly different from those in Casares et al. (2012). If this were the case, then the mechanism proposed by us could be hardly behind the dominant lightcurve peak of HESS J0632+057, regardless of the eccentricity. In the particular case the periastron were ∼ 0.3 in phase earlier, the narrow GeV and TeV predicted peaks may actually be the first GeV and TeV peaks of the observed lightcurve, and the broader predicted X-ray and TeV peaks could correspond to the second observed X-ray and TeV peaks of the lightcurve. In that case, the first X-ray peak of the lightcurve may need adding an extra emitting components (e.g. the binary region, Be disc interaction, etc.). In this scenario, the radio behaviour would be more difficult to explain though.
A next step would be to carry out more accurate radiation calculations applying post-processing to the hydrodynamical solution for HESS J0632+057, and compute consistently the particle evolution in space and energy in a way similar to those employed for instance by Dubus et al. (2015) and de la Cita et al. (2017). Another step forward would be to increase the resolution in the θ-coordinate when simulating the wind collision at large scales in this source. Note that including the Be disc in the simulations would require including the binary region and a Cartesian grid, making simulations much more demanding.