Decade-long X-ray observations of HESS J0632+057

We present results of a decade-long X-ray observations of gamma-ray loud binary HESS J0632+057 with Swift , XMM-Newton, Chandra and Suzaku. Basing on all available data we refine the orbital period of the system to be 316.2 −2.0 d (95% c.r.), consistent with the previous studies although measured with significantly better accuracy. We report on a hydrogen column density and spectral slope variability along the orbit. We argue that the observed variability can be explained within a “similar to PSR B1259-63” model in which the orbit of the compact object is inclined to the disk of Be star. We show that the observed X-ray to TeV emission can originate from a broken cut-off powerlaw population of electrons and describe the way in which future X-ray/TeV observations can distinguish between proposed and flip-flop emission model of this system.


INTRODUCTION
HESS J0632+057 is one of a few known gamma-ray loud binaries (GRLBs). All known systems of this class host compact object orbiting either O-or Be-type star characterised by a presence of dense circumstellar disk (see e.g. Dubus 2013, for the review). Peculiarly the spectral energy distributions (SEDs) in GRLB systems are similar and dominated by the emission in GeV-TeV energy band. Although the nature of this similarity is not yet completely understood, it allows to suggest a similar nature of the compact object and physical mechanisms operating in the systems, probably with the corrections to geometry of each individual system. At least in one GRLB, PSR B1259-63, compact object is known to be a young pulsar (Johnston et al. 1992), which allows to suggest the same type of a compact object for several other binaries (Neronov & Chernyakova 2008;Torres et al. 2010;Zdziarski, Neronov & Chernyakova 2010;Durant et al. 2011;Moritani et al. 2015;Bird & for the VERITAS Collaboration 2017), see however Massi et al. (2001); Casares et al. (2005); Williams et al. (2010); Massi, Migliari & Chernyakova (2017) arguing black hole nature of a compact object in other systems.
Contrary to other gamma-ray loud binaries, HESS J0632+057 remained the only system missed in GeV energy band until very recent hints of its detection with Fermi/LAT (Malyshev & Chernyakova 2016;Li et al. 2017). The brightness of the system in TeV band allowed its initial discovery during HESS observations of Monoceros region ) as unidentified point-like source. Its spatial coincidence with Be star MWC 148 allowed to suggest its binary nature ; Hinton et al. 2009). After dedicated observations the system was detected also in radio (Skilton et al. 2009) and soft X-ray (Falcone et al. 2010) energy bands. The system was also detected in TeV band by VERITAS and MAGIC (Aleksić et al. 2012;Aliu et al. 2014) observatories.
During the orbital cycle HESS J0632+057 demonstrates in X-rays two clear peaks of the emission -first at phase φ ∼ 0.2 − 0.4 and second at φ ∼ 0.6 − 0.8 separated by a deep minimum at φ ∼ 0.4 − 0.5 (Bongiorno et al. 2011;Aliu et al. 2014). Recently, the similar structure of the orbital lightcurve in TeV energy range was reported by Maier et al. (2015). Hints of orbital variability in GeV range were reported in Li et al. (2017).
In this work we summarize archival observations of HESS J0632+057 performed by Swift/XRT during the last decade and accomplish them with results of the most recent Swift TOO observations of a first emission peak in 2017. We also use historical XMM-Newton, Chandra and Suzaku observations. Based on these data we significantly improve the determination of the orbital period of the system and report the variability of hydrogen column density nH and the spectral slope along the compact object's orbit. We propose to explain the observed variability within a model in which the orbit of the compact object is inclined to the disk of Be star.  Table 1. Summary of analyzed XMM-Newton,Chandra and Suzaku observations. The data grouped were fitted in 0.3 − 10 keV range with absorbed powerlaw model(cflux*phabs*po) with hydrogen column density n H , slope Γ 0.3−10 and flux in 0.3-10 keV range F 0.3−10 .

DATA ANALYSIS
In this work we re-analyzed all available X-ray archival data on HESS J0632+057 in order to guarantee the usage of the most recent calibration files and data analysis software. Where applicable, the described below data reduction was performed with the most recent available heasoft v.6.22 software package and spectral analysiswith XSPEC v.12.9.1m. All spectra were grouped to a minimum of 1 count/bin and fitted in 0.3-10 keV energy range with the absorbed power law model (cflux*phabs*po) using c-statistics suitable for the analysis of poor-statistics data 1 . The best-fit parameters (flux, slope and hydrogen column density) are summarized in Table 1 and shown in Fig. 2 Within statistical uncertainties the obtained results are consistent with ones reported previously.

Swift/XRT data analysis
For the analysis below we have considered all publicly available Swift/XRT data on HESS J0632+057 taken between 2009, January, 26th and 2017, April, 17th. The data was reprocessed with xrtpipeline v.0.13.4 as suggested by Swift/XRT team 2 . The spectrum was extracted with xselect at nominal coordinates of HESS J0632+057 ((ra;dec) = (98.246894 ; 5.800327) ) from 36" circle and the background -from annulus centred at HESS J0632+057 position and inner(outer) radii of 60"(300"). Determining the best-fit flux for each observation, we have selected for further analysis only observations in which the uncertainty of the flux does not exceed 50% of the flux value. This criterium allowed us to select for further analysis 239 out of 242 total Swift observations of HESS J0632+057.

XMM-Newton data analysis
The analysis of a single XMM-Newton observation of HESS J0632+057 (taken in 2007, Sept. 17 and labelled hereafter "X1") was performed with the latest XMM-Newton Science Analysis software v.15.0.0. Known hot pixels and electronic noise were removed, the data were cleaned 1 See description of used by XSPEC statistics 2 See e.g. Swift/XRT User's Guide from the influence of the soft proton flares, which strongly affected the end of the observation. Applied cleaning resulted in ∼ 30 ksec of flare-clean exposure. During the observation HESS J0632+057 was located close to the border of PN camera chip and the corresponding data were more affected by soft proton flares. In order to keep the analysis conservative in what follows we use only the data from MOS cameras. The spectrum was extracted from 40" radius region centered at the position of HESS J0632+057 and the background was extracted from nearby 80" radius source-free region. The RMFs and ARFs were extracted using RMFGEN and AR-FGEN tools respectively.

Chandra data analysis
Chandra has observed HESS J0632+057 during its high state in 2011, Feb. 13. We analyzed the data using the most recent CIAO v.4.9 software and CALDB 4.7.6. The data was reprocessed with chandra repro utility, source and background spectra with RMF and ARF, were extracted with specextract tool. The observation was performed in asiccc (continuous clocking) mode, in which only 1-dimensional spatial information is available. To extract the source and background spectra we thus utilised box-shaped regions 3 centered on HESS J0632+057 and located in nearby sourcefree region. Total clean exposure for "C1" observation was ∼ 40 ksec

Suzaku data analysis
Suzaku performed two observations of HESS J0632+057 in 2008, Apr. 23 and 2009 Apr. 20. In what follows we will label these observations "S1" and "S2" correspondingly. For the analysis we used XIS 0,1,3 data reprocessed with aepipeline v.1.1.0 tool. Source(background) spectrum as well as corresponding response files were extracted with xselect tool 4 from a circular(annulus) region centered at HESS J0632+057 and radius of 150" (151" and 288" inner and outer radii correspondingly). We find also that the source is not detected with Suzaku HXD/PIN with the upper limits consistent with power law extrapolation of XIS data model. Shaded area illustrates 95% confidence region (statistic and systematic uncertainties) for the orbital period.

Orbital periodicity search
The data described in Sec. 2 covers about 10 orbits of HESS J0632+057 including the most recent Swift/XRT TOO observations taken during high-flux state of the source. This allowed us to perform an improved, in comparison to Bongiorno et al. (2011) and Aliu et al. (2014), search for the orbital periodicity in this system.
We performed autocorrelation analysis for a periods ranging in 250 d -350 d with the step of 0.1 d. For each of these periods at the first step of the analysis we have binned the data in 20 equal orbital-phase bins defining the flux in each bin as the weighted mean flux over all attributed to this bin observations. At this point we explicitly assume, that the orbital profile is smooth on scales of ∼ 0.05 orbit and the average number of observations per selected bin is large enough for uncertainty on mean profile to be significantly smaller than the uncertainty of each individual measurement.
Basing on this binned orbital flux profile we define smoothed orbital flux profile as a linear interpolation between fluxes in the bins with periodic boundary conditions. The considered smoothed orbital flux profile allowed us to define model fluxes at times of each observation of analyzed data. At the second step of the analysis we built the Pearson correlation coefficient between model and (mean) observed fluxes as a function of considered orbital period, see Fig. 1. The correlation coefficient demonstrates clear maximum at orbital period P orb = 316.2 d shown with vertical solid green line.
In order to estimate the uncertainty for P orb we have considered 10 4 random-trial realizations of the fluxdatapoints basing on observed fluxes values and corresponding uncertainties. For each realization r of the data we repeated described above analysis and determine P orb (r) at which maximum of Pearson coefficient is reached. We define 95% confidence range for P orb as a region containing 95% of all P orb (r) (0.025 -0.975 quantiles). This region is shown with shaded green area in Fig. 1 and corresponds to P orb = 316.2 +0.5 −1.6 d. We checked for the systematic effects of the selected number of phase bins on measured orbital period varying the number of bins from 10 to 40 and performing random trial analysis described above. We define systematic-accounted 95% confidence range as an interval which includes all obtained at this step 95% confidence intervals for P orb . The systematic uncertainty is comparable with statistic one and account for it gives P orb = 316.2 +0.5stat+1.3syst −1.6stat−0.4syst d. The observed flux data convolved with best-fit P orb is shown in Fig. 2, left panel. Colour points marked with "O N " stand for corresponding fluxes observed during the orbit N (T0 = MJD 54587.0 as in Bongiorno et al. (2011)). Green points with the stars and right panel show the averaged over all data points flux orbital profile. We would like to note that for the measured value of P orb two high-quality observations of Chandra and XMM-Newton appeared to be taken during the first narrow X-ray maximum(C1) and close to a deep minimum(X1). In what follow we will refer to these observations as high-and low-state observations of HESS J0632+057 correspondingly.

Hydrogen column density and slope orbital profiles
Besides flux orbital profile Swift/XRT data allow to trace orbital variation of the spectral slope Γ and hydrogen column density nH . We split the orbit of HESS J0632+057 over 10 bins and grouped Swift observations according to these bins. For each orbital bin we fit observations with an absorbed power law model fixing the slopes and hydrogen column densities to be the same for all observations, while the flux levels were allowed to vary between observations. The obtained orbital profiles for nH and Γ are shown in Fig 3 with the green points. Black points show the high-quality results from individual observations of XMM-Newton, Chandra and Suzaku. Thin solid blue lines illustrate the best fit of the corresponding orbital profile with a constant (only Swift data taken into account). The constant orbital profiles can be excluded at significance level of ∼ 4.3σ(nH orbital profile) and ∼ 3.7σ(spectral slope orbital profile).
The results of XMM-Newton, Chandra and Suzaku are in a good agreement with Swift mean orbital profile.

DISCUSSION
Several models were proposed in the literature to explain the observable keV -TeV emission of HESS J0632+057. In flip-flop scenario (see e.g. Moritani et al. (2015) and references therein) the compact object is assumed to be a pulsar. Close to the apastron (orbital phases ∼ 0.4 − 0.6, see Fig. 2) the pulsar is in a rotationally powered regime, while switches into propeller regime approaching the periastron(phases 0.1 − 0.4 and 0.6 − 0.85). In a flip-flop system, if the gas pressure of the Be disk overcomes the pulsar-wind ram pressure, the pulsar wind is quenched (phases 0 − 0.1  and 0.85−1). Because the Be disk of the system is estimated to be about three times larger than the binary separation at periastron, the compact object crosses a dense region of the disk near the periastron. In such a situation, the strong gas pressure is likely to quench the pulsar wind and suppress high-energy emissions.
Alternatively, one can consider the model similar to proposed by Chernyakova et al. (2015) for PSR B1259-63 system. Similar two-peaks X-ray orbital behaviour of these systems allows to assume, that the orbital plane of a compact object in HESS J0632+057 is inclined to the Be star's disk, similar to PSR B1259-63. The X-ray/TeV peaks within this model correspond to the first and second crossing of the disk by a compact object. Higher ambient density during these episodes leads to more effective cooling of the relativistic electrons by synchrotron and inverse Compton mechanisms.
Within this interpretation the positions of the peaks al-low to estimate the position of periastron. Due to higher orbital velocities of the compact object at close to periastron phases, the periastron has to be located at shortest orbital separation between the peaks, i.e. at phases 0.4 -0.5. Relative widths of the peaks allow also to estimate the eccentricity of the orbit to be ∼ 0.5, see Fig. 4 for the sketch of the orbit. We would like to note, that these parameters differ significantly from ephemeris of Casares et al. (2012), but are in a good agreement with the orbital parameters obtained with very recent optical monitoring of the system by Moritani et al. (2017). for LS I+61 303. Within this model at each orbital phase observed nH value is given by the integration of the inclined Be star disk density profile along the line of sight to the observer. We assume that the disk consists of non-ionized hydrogen and has an exponential density profile characteristic of isothermal atmosphere. For such a simple case the sum of constant galactic nH and the model one is shown with dashed blue line in Fig. 3

, left.
Within this simple model the double-peaked nH profile can not be expected for non-inclined disk with monotonic density profile usually suggested in "flip-flop" model. Nonmonotonic complex dynamic density distribution can be formed due to strong disk-compact object outflow interactions (Bosch-Ramon et al. 2017, seen in hydrodynamic simulations) and may explain double-peaked nH profile. Note however, that strong compact object outflow at all orbital phases has been suggested by Bosch-Ramon et al. (2017), while in "flip-flop" model at close to periastron phases this outflow is expected to be minimal.
The spectral energy distribution of HESS J0632+057 in keV-TeV band is shown in Fig. 5 for high(observation C1) and low(observation X1) states. It's peculiar GeV-TeV shape suggests the presence of a break or a cutoff located at ∼ 140 − 200 GeV with a lower-energy slope < 1.6 (Malyshev & Chernyakova 2016). The presence of such a break can be interpreted as a corresponding break in the spectrum of the emitting particles arisen from energy dependence of cooling losses. Generally the value of break energy E br can vary along the orbit.
The SEDs of the system in high and low states can be also described in similar to PSR B1259-63 terms. In Fig. 5 we show the model spectrum expected from a broken cutoff powerlaw distribution of electrons calculated with naima v.0.8 python package (Zabalza 2015). X-ray/GeV and TeV parts of the spectrum are produced via synchrotron and inverse Compton mechanisms shown in Fig. 5 with dashed 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 10 11 10 12 10 13 E, eV and dot-dashed lines correspondingly. Very high energy and Chandra observations of HESS J0632+057 in a high state allow to fix the high energy slope of electron spectrum to Γe,2 ≈ 2.3, the position of the break to E e,br ≈ 0.7 TeV (for IC scattering on photons of Be star with the temperature T∼ 3 · 10 4 K) and the high-energy cutoff Ecut ∼ 100 TeV. The magnetic field strength in the emitting region can be estimated as B ≈ 0.1 G. In calculation we assumed the emitting region to be located at a distance d ∼ 5 a.u. from Be star with luminosity L * = 3 · 10 38 erg/s. Fixing low-energy slope of the electron distribution to one matching XMM-Newton low-state spectrum, namely Γ1,e ≈ 1.3, we can satisfy Fermi/LAT upper limits in 100 GeV range. The softening of electron spectrum by ∆Γ 1 is typical for the effect of the synchrotron cooling on the spectrum (see e.g. Chernyakova et al. 2015, for PSR B1259-63 case). The synchrotron cooling modifies the electron spectrum in the energy range above 1 TeV and fails to do this below 1 TeV. The absence of cooling in the energy band below 1 TeV could be attributed to the escape of the sub-TeV electrons from the system. The observation of the hard slope in a low HESS J0632+057 state requires the increase of E e,br as the source transiting from its high to low state which does not allow XMM-Newton to detect the break in soft X-ray band. This can be naturally understood as the low state is observed at phases when the compact object leaves the disk of Be star. Relativistic electrons in this case could escape from the system mush faster and do not suffer from synchrotron losses rather than in high state when the compact object remains in the disk.
Alternatively the spectral break in electron spectrum can occur at transition between adiabatic and IC/synchrotron loses dominant energies, see e.g. Khangulyan et al. (2007) and Takahashi et al. (2009) for PSR B1259-63 and LS 5039 cases. The adiabatic loss time is naturally shortest in sparse regions outside of Be star disk and longest in dense regions inside. This also allow to explain the shift of the spectral energy break to higher energies in a low state.
For both interpretations it is naturally to expect, that at phases when the compact object approaches the dense regions of the disk (i.e. at phases φ ∼ 0.3 and φ ∼ 0.6) the break energy is located between its maximal and minimum values observed in high (deep in the disk) and low(out of the disk) states and can be detected as a break in X-ray photon spectrum. We find that Suzaku data taken at phase ∼ 0.3 (S2 observation) marginally prefers (∆cstat = 5 for 1 d.o.f.) the model with the break prior to a single powerlaw model with the indication of a break present at ∼ 5 keV.
The non-powerlaw shape of the spectrum at these phases can systematically affect spectral slopes measured by Swift. This can explain some discrepancy between mean over several orbits slope measurement by Suzaku and individual measurement of the slope by Chandra at phase ∼ 0.4. The observed discrepancy can be also explained either by strong gradients in slope orbital profile or by a systematically varying from orbit to orbit spectral slopes. Similar broken powerlaw shape of the spectrum can be expected also within "flip-flop" model, since both interpretations of the break origin can be valid for this model. The observed hardening of the spectrum close to apastron(X1 observation) is as well expected in both, flip-flop and "similar to PSR B1259-63" models since it is connected with out-of-the-disk position of compact object identical in both models. These models however can be distinguished by the observations of the position of the spectral slope/break at close to periastron phases and hydrogen density profile.
In flip-flop model it is naturally to expect that the escape time from the system (as well as adiabatic cooling time) increases as the compact object enters dense disk regions and approaches the periastron. This should lead to the shift of the break to lower energies. Contrary, in "similar to PSR B1259-63" model, the compact object during the periastron passage is away or in the sparse regions of the Be star's disk. This allows to expect relatively short escape/adiabatic cooling times and suggests high E e,br . Within this model the longest escape times (lowest E e,br ) are expected during the disk's crossing by a compact object, i.e. coincide with the maxima of X-ray/TeV activity of the system.
One of the key points of "flip-flop" model is switching of the compact object into accreting regime as it approaches the periastron which unavoidably requires the high ambient density at this point. The double-peak nH profile, if confirmed, suggests the low periastron density for thin inclined disk geometry.
Thus, the detailed observations of E e,br position and hydrogen column densities orbital profiles can allow to distinguish between these models. Note, that break-energy observations can be performed either in TeV or X-ray energy band by observation of a corresponding feature in synchrotron or IC component of SED.