Observation of the Gamma-Ray Binary HESS J0632+057 with the H.E.S.S., MAGIC, and VERITAS Telescopes

The H.E.S.S. (High Energy Stereoscopic System) experiment is an array of five imaging atmospheric Cherenkov telescopes located in the Southern Hemisphere in the Khomas Highland of Namibia ($23^\circ {16}^{{\prime} }$ S, $16^\circ {30}^{{\prime} }$ E), at an elevation of 1800 m above sea level. Four telescopes of the array (CT 1–4) have 12 m diameter mirror dishes. A larger, 28 m diameter telescope (CT 5) was deployed at the site in 2012 July during the H.E.S.S. phase-II upgrade.

The telescopes CT 1–4 have been in operation since 2004 (Hinton & the HESS Collaboration 2004) and have a 5° diameter field of view (FoV). CT 5 is equipped with a camera with a 3.5° diameter FoV. Using only CT 1–4, as is done in this work, H.E.S.S. is sensitive to gamma-rays with energies from ∼0.2 TeV up to tens of TeV for observations at zenith.

The H.E.S.S. observations of HESS J0632 + 057 cover almost 14 yr (2004 March–2018 February), totalling 120 h of data after quality cuts, dead-time corrected and taken under dark sky conditions (Table 1). Analysis of 73 h of data obtained prior to 2012 February (MJD 55,951) was reported in Aliu et al. (2014) and these data are reanalyzed in this work, applying improved calibration and reconstruction methods for this paper. All available H.E.S.S. data were processed with the Image Pixel-wise fit for Atmospheric Cherenkov Telescope (ImPACT) analysis (Parsons & Hinton 2014) with std_ImPACT cut configuration, which requires a minimum of 60 photoelectrons per image.

The signal was extracted from a circular region of radius 0.07° centered on the position of HESS J0632 + 057. For the background estimation, the reflected background method (Fomin et al. 1994; Berge et al. 2007), which ensures that the regions used for signal (ON) and background (OFF) extraction have the same acceptance in the FoV of the camera, was applied. A cross-check of the main analysis was performed using the Model++ analysis (de Naurois & Rolland 2009), yielding compatible results.

The MAGIC (Major Atmospheric Gamma-Ray Imaging Cherenkov Telescopes) experiment is a stereoscopic system of two 17 m diameter imaging atmospheric Cherenkov telescopes. It is located on the Canary Island of La Palma, Spain, at the observatory El Roque de Los Muchachos ($28^\circ {45}^{{\prime} }$ N, $17^\circ {53}^{{\prime} }$ W, 2200 m above sea level). It is equipped with pixelized cameras covering a field of view with a diameter of ∼35. The MAGIC telescopes reach a trigger threshold down to 50 GeV and a sensitivity of ∼ 0.6% of the Crab Nebula flux above 250 GeV in 50 hr of observations at zenith angles <30° (Aleksić et al. 2016).

MAGIC observed HESS J0632 + 057 from 2010 October until 2017 November for a total of 68 hr (Table 1). The data set taken during an X-ray outburst occurring from 2010 October to 2011 March was published in Aleksić et al. (2012). The most recent data cover 2015 December to 2017 November. After quality cuts and correction for dead-time, 57.3 hr of good data remain. These observations were taken both under dark conditions (30.9 hr) and under moderate-to-strong moonlight (26.4 hr). They were carried out in wobble mode, pointing at two regions symmetrically 0.4° away from the source position to simultaneously evaluate the background (Fomin et al. 1994). The data have been analyzed using the standard MAGIC procedure (Zanin et al. 2013). The recorded images were calibrated, cleaned, and parameterized (Hillas 1985; Albert et al. 2008a), while applying different cleaning levels and size cuts for the data in moderate-to-strong moonlight according to Ahnen et al. (2017). The background rejection and the estimation of the gamma-ray arrival direction were performed using the random forest method (Albert et al. 2008b).

The VERITAS (Very Energetic Radiation Imaging Telescope Array System) telescope array is composed of four imaging atmospheric Cherenkov telescopes located at the Fred Lawrence Whipple Observatory in southern Arizona (1268 m above sea level, 31°40'N, 110°57'W; Weekes et al. 2002). Each telescope has a 12 m diameter segmented mirror of ≈100 m² area. Cherenkov light emitted in extensive air showers is reflected onto photomultiplier cameras covering a field of view of 35 diameter.

The VERITAS observations of HESS J0632 + 057 (VER J0633+057) have been obtained over 11 yr, from 2008 December to 2019 January. For a detailed discussion of the evolution of the performance of the VERITAS instrument with time, see Park & VERITAS Collaboration (2015).

VERITAS observed HESS J0632 + 057 for a total dead-time corrected exposure time of 260 hr, after the application of standard data selection cuts to account for nonoptimal weather or hardware conditions (see Table 1 for details). HESS J0632 + 057 was observed at an average elevation angle of 60° and at a fixed offset of 05 from the center of the camera's field of view. Observations were taken typically in dark sky conditions, with the exception of a small fraction (∼10 hr) of the data set which was obtained during 30%–65% moonlight. The latter results in a higher energy threshold and reduced sensitivity compared to observations taken under nominal sky conditions (Archambault et al. 2017).

Data were analyzed using VERITAS standard analysis software tools utilizing boosted decision trees for gamma–hadron separation based on image parameters (see, e.g., Krause et al. 2017; Maier & Holder 2017). The source region was defined as a circle with radius 009 around the nominal position of the X-ray source XMMU J063259.3 + 054801. Background event rates are estimated using the reflected background method (Fomin et al. 1994), assuming typically six background regions located symmetrically with respect to the signal region around the camera center. Light-intensity calibration factors obtained from regular monitoring of the optical throughput and of the detector performance are applied to take time-dependent changes in the instrument response into account (Nievas & VERITAS Collaboration 2021). The observations up to 2017 March (MJD 57,822) have been presented before (Aliu et al. 2014; Maier & VERITAS Collaboration 2017). These data have been reanalyzed applying improved calibration and background suppression methods.

This study includes data obtained with the Swift-XRT, XMM-Newton, Chandra, and NuSTAR instruments. All data presented in Malyshev et al. (2019) were used, adding also Swift-XRT and Chandra data obtained after the publication of the aforementioned paper. Here, only details of analysis for these additional data are provided; for a description of the analysis of NuSTAR and Suzaku data, see Sections 2.4 and 2.5 of Malyshev et al. (2019).

The Swift-XRT data taken between 2018 December 5 and 14 were analyzed with the heasoft v.6.22 software package and reprocessed with xrtpipeline v.0.13.4, as suggested by the Swift-XRT team. ¹¹⁹ The spectra were extracted with xselect, using a $36^{\prime\prime}$ circle for source counts and an annulus centered at the source position with inner/outer radii of $60^{\prime\prime}$/$300^{\prime\prime}$ for the background. Chandra data were analyzed using CIAO v.4.9 software and CALDB 4.7.6. The data were reprocessed with the chandra_repro utility, source and background spectra with corresponding redistribution matrix files (RMF), and ancillary response files (ARF) were extracted from circular regions of radii $11^{\prime\prime}$ (source) and $50^{\prime\prime}$ (background) with the specextract tool. Two Chandra observations (obs. IDs 20974 and 20975; total exposure ∼70 ks) were performed over 2018 February 19–21, about two weeks after the last Chandra observation presented in Malyshev et al. (2019).

For each analyzed X-ray energy spectrum, neighboring energy bins were merged until each bin contained at least one count. Obtained spectra were then fitted in the 0.3–10 keV energy range. The reported fluxes were extracted from the fit of the spectra with an absorbed power-law model (cflux*phabs*po) with XSPEC v.12.9.1m. The fit was performed with Cash statistics, suitable for the analysis of low-statistics data (Cash 1979).

High-dispersion Hα profiles of the Be star in HESS J0632 + 057 were obtained using the HIDES (High-Dispersion Echelle Spectrograph), a fiber-fed system (Kambe et al. 2013) deployed on a 188 cm telescope at Okayama Astrophysical Observatory (OAO). ¹²⁰ Furthermore, spectra using ESPaDOnS (Echelle SpectroPolarimetric Device for the Observation of Stars) at the Canada–France–Hawaii Telescope (Manset & Donati 2003) from 2013 October 31 to 2017 April 13 were obtained; see Table 2. The resolving power, R, of HIDES and ESPaDOnS is ∼50,000 and ∼68,000, respectively.

The observations and the data reduction steps are described in detail in Moritani et al. (2018); a short summary is provided here. The HIDES data were reduced using IRAF. The ESPaDOnS data were reduced using the Libre-Esprit/Upena pipeline, provided by the instrument team. The continuum level around the Hα line was refitted and the spectrum was normalized in order to compare with the HIDES data.

The radial velocity was measured as the bisector velocity of the wing of the Hα profile (Moritani et al. 2018) and the FWHM was determined by fitting the middle and bottom of the profiles with a single Gaussian function. The equivalent width (EW) was estimated by simply integrating the profiles. The observations are summarized in Table 2.

The long-term light curves (Figure 2) confirm that HESS J0632 + 057 is highly variable at X-ray and gamma-ray energies (see also Acciari et al. 2009; Bongiorno et al. 2011; Aliu et al. 2014). Above 350 GeV, the flux varies between 5% and 10% with respect to the flux of the Crab Nebula (above the same energy), dropping below the detection limits of the different ground-based instruments for some periods.

Zoom In Zoom Out Reset image size

The large data sets allow us to search for periodic variability patterns on different timescales. The period was first determined to be 321 ± 5 days by Bongiorno et al. (2011) using approximately two years of Swift-XRT measurements (MJD 54,857 to 55,647). Subsequent updates to this analysis, based on extended data sets from Swift-XRT, yielded comparable results of P_orb = 313–321 days (Aliu et al. 2014; Moritani et al. 2018; Malyshev et al. 2019). Flux modulations of the X-ray light curve are interpreted in the literature as being due to the orbital period of the binary system. The only non-X-ray measurement of the orbital period, based on optical observations of Hα radial velocity profiles, yielded ${P}_{\mathrm{orb}}={308}_{-23}^{+26}$ days (Moritani et al. 2018), again compatible with previous results (see Table 3 for a summary of the measurements of orbital period). Although Casares et al. (2012) did not estimate the orbital period, they showed that the main parameters of the Hα emission line (equivalent width, FWHM, and centroid velocity) are modulated within the 321 day X-ray period and calculated an orbital solution of the binary.

Table 3. Results of Derivation of the Orbital Period from Literature in Comparison with the Outcome of This Work

Publication	Energy Range	MJD Range	Method	Orbital Period
				(days)
Bongiorno et al. (2011)	Swift-XRT	54,857–55,647	Peak fitting	321 ± 5
Aliu et al. (2014)	Swift-XRT	54,857–55,972	ZDCF	${315}_{-4}^{+6}$
Moritani et al. (2018)	Swift-XRT	54,857–57,052	ZDCF	${313}_{-8}^{+11}$
Moritani et al. (2018)	Hα	56,597–57,856	Fourier analysis	${308}_{-23}^{+26}$
Malyshev et al. (2019)	Swift-XRT	54,857–57,860	PCC	${316.2}_{-2.0}^{+1.8}$
This work	X-ray	54,857–58,466	PCC	317.3 ± 0.7(stat) ± 1.5(sys)
This work	X-ray	54,857–58,466	PDM	316.7 ± 1.5(stat) ± 1.5(sys)
This work	X-ray	54,857–58,466	DCF	318.9 ± 2.2(stat) ± 1.5(sys)
This work	Gamma-ray	53,087–58,490	PCC	316.7 ± 4.4(stat) ± 2.5(sys)
This work	Gamma-ray	53,087–58,490	PDM	319.0 ± 3.4(stat) ± 2.5(sys)

Note. The applied methods are: peak fitting (Bongiorno et al. 2011), Z-transformed discrete correlation functions (Alexander 1997, ZDCF), phase dispersion minimization (Stellingwerf 1978, PDM), discrete correlation functions (Edelson & Krolik 1988, DCF), and correlation analysis comparing the light curves with a binned-average light curve (Malyshev et al. 2019, PCC). Errors indicate 1σ statistical uncertainties.

Download table as:  ASCIITypeset image

The analysis of the most recent available Swift-XRT data set including the reanalysis of older data (264 data points for MJD 54,857–58,466; see Figure 2, top), considered in this study using a method based on Pearson's correlation coefficient (PCC method; see Appendix A), results in an orbital period of P_orb = 317.3 ± 0.7(stat) ± 1.5(sys) days (see Table 3). Systematic uncertainties are derived from the application of the different methods for the period derivation applied to 1000 Monte Carlo (MC)-generated data sets. The description of the application of additional methods, such as discrete correlation functions (DCF) or phase dispersion minimization (PDM), together with details of the analysis of the period and derivation of the associated systematic uncertainties, can be found in Appendix A. Periodic signals on timescales different from the orbital period have been observed in high-mass X-ray binaries, e.g., the superorbital period of 1667 days detected in radio, optical, X-ray and gamma-ray observations of the binary LS I +61° 303 (Gregory 2002; Li et al. 2011; Ackermann et al. 2013; Paredes-Fortuny et al. 2015b; Ahnen et al. 2016; Kravtsov et al. 2020). A signal at periods other than the orbital period has not been found for periods between 50 and 1000 days in an analysis of the X-ray observations of HESS J0632 + 057 using the PCC method (see Appendix A and Figure 12(d)). The analysis of MC-generated light curves with similar time structures shows that the PCC periodogram obtained from the analysis of the measured light curve is consistent with expectations from statistical fluctuations.

The large gamma-ray data set, spanning almost 15 yr of gamma-ray observations (Figure 2, bottom), allows the orbital period of this system to be determined from gamma-ray data alone for the first time. The approach of the analysis is equivalent to that of the X-ray analysis discussed above (PCC method), resulting in an orbital period P_γ = 316.7 ± 4.4(stat) ± 2.5(sys) days, consistent with the value obtained at X-ray energies (see again Appendix A for details on the analysis). Given the known strong correlation of X-ray and VHE gamma-ray fluxes observed in the past (Aliu et al. 2014, see also Subsection 3.3), this result is not unexpected.

Table 3 summarizes the outcome of the different analyses of orbital period in comparison with results from the literature. For period-folding, the result from the PCC method was used in this paper. It was applied to the Swift-XRT data with the following values: MJD₀ = 54,857 (arbitrarily set to the date of the first observations of the source with Swift) and period P = 317.3 ± 0.7(stat) ± 1.5(sys) days. It should be stressed that none of the conclusions on the physical properties of the binary system depend on this particular choice (see also Appendix B and Figures 13 and 14).

The light curves at X-ray and gamma-ray energies were both folded with the adopted orbital period of 317.3 days using ϕ = 0 as MJD=54,857.0, arbitrarily set to the date of the first Swift-XRT observations (Bongiorno et al. 2011). These folded light curves reveal several prominent features (see Figure 3): a maximum around phases 0.2–0.4, a minimum at phases 0.4–0.6, a second (lower) peak around phases 0.6–0.8, and a broader plateau around 0.8–0.2. HESS J0632 + 057 is now detected at energies >350 GeV with high statistical significance at all orbital phases with the exception of the region around the first minimum, at phases 0.4–0.6 (Table 4). According to the orbital solution of Casares et al. (2012), the periastron is located at phase ϕ = 0.967, coincident with the gamma- and X-ray flux plateau, while apastron takes place during the dip/minimum, at phase ϕ = 0.467. For the orbital solution proposed by Moritani et al. (2018), periastron takes place at ϕ = 0.663, during the second peak, while apastron happens at ϕ = 0.163 (see Figure 1 for a representation of the suggested solutions). The phase-folded parameters of the Hα measurement results are shown in Figure 4. No significant modulation is seen in EW, whereas the radial velocity and FWHM show small variations as a function of orbital phase. Note that Moritani et al. (2018) proposed an orbital solution using the radial velocity of Hα.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HESS J0632 + 057 has been observed over 14 orbits at X-ray and 18 orbits at gamma-ray energies (see Appendix C and especially Figures 15 and 16). Despite more than 270 X-ray pointings and ∼450 hr of gamma-ray observations, these data cover most orbits relatively sparsely. Only five orbits (orbits 7–11 in Figure 15) have been regularly monitored at X-ray energies over more than 50% of an individual orbit. Ground-based gamma-ray instruments covered an even smaller range, as they cannot observe the binary for a large fraction of its orbit due to visibility constraints. Four orbits have been observed in the phase range 0.2–0.4 with notable coverage in X-ray and gamma-ray observations (Figure 9) and are discussed in more detail in Section 3.4.1.

Correlations between X-ray and gamma-ray flux measurements are shown in Figure 5 (top). Gamma-ray observations are averaged over several nights in order to achieve sufficient statistical detection significance, while X-ray observations are nightly averages. Measurements are considered to be contemporaneous in this context for X-ray observations taking place during the time period of the grouped gamma-ray observations (the maximum interval length for the grouped gamma-ray observations is ∼10% of the orbital period). Average X-ray fluxes are used for periods of gamma-ray observations including several X-ray data points. It is worth stating that both the X-ray and gamma-ray fluxes are measured to be variable on timescales as short as hours, and the total X-ray exposure associated with each point is typically just a small fraction of the gamma-ray exposure. This might contribute to some of the scattering between fluxes in both energy bands as observed in Figure 5 (top). The DCF method (Edelson & Krolik 1988) is used to quantify the correlation and to search for possible time lags between X-ray and gamma-ray emission patterns. Figure 5 (bottom) shows the DCF using a binning of 5 days (different binnings have been tried with no significant changes in the results) and with confidence limits calculated using the toy Monte Carlo technique as recommended in Uttley et al. (2003). The correlation coefficient is 0.82 at time lag τ = 0 and the p-value for noncorrelation is 4 × 10⁻¹⁵, calculated from 59 flux pairs. The ratio of gamma-ray (>350 GeV) to X-ray flux (0.3–10 keV) is on average 2.9 ± 0.3, estimated from a linear fit as shown in Figure 5 (top) for the chosen integration ranges in energy. This underlines the equality or even dominance of the gamma-ray energy range for the emission of the binary system with respect to the X-ray regime (see also Figure 8). The fit also indicates a nonzero X-ray flux (F_{0.3−10 keV} = (6.1 ± 1.5) × 10⁻¹³ erg cm⁻² s⁻¹) for a vanishing gamma-ray component, suggesting an X-ray source partially unrelated to the gamma-ray emission. No indication of a significant time lag between X-ray and gamma-ray emission is observed (see Figure 5, bottom), although this measurement is not sensitive to time lags τ ⪅ 10 days due to the required grouping of the gamma-ray observations. In general, this confirms the strong correlation between gamma-ray and X-ray emission, as already shown in Aliu et al. (2014). Ranked and unranked correlation analyses result in similar correlation coefficients between 0.75 and 0.82. This points toward a close-to-linear correlation, although statistical uncertainties of the photon fluxes observed in both energy bands prevent us from making any strong statements regarding possible nonlinearities between gamma-ray and X-ray emission.

Zoom In Zoom Out Reset image size

No correlation is found between any parameters obtained from the optical Hα measurements and X-ray or gamma-ray fluxes (Figure 6). The Spearman correlation coefficients are between –0.4 and 0.33 for time lags τ = 0, corresponding to p-values for noncorrelation between 0.02 and 0.85 (calculated from 18 gamma-ray/Hα and 31 X-ray/Hα pairs). A relaxed definition for contemporaneous data is used for the Hα correlation analysis with a time span of ±5 days, well below the shortest timescales of ≈50 days observed for Hα profile variability (Moritani et al. 2015). This relaxation of the contemporaneity criteria is necessary to obtain a reasonable number of pairs of observations, but might obscure variability on shorter timescales. The low gamma-ray fluxes and large uncertainties might hide possible correlations. It should also be noted that the optical Hα data were obtained during five orbital periods (MJD 56,597–57,857; orbits 12–16 in Figures 15 and 16) with generally sparse coverage of X-ray and gamma-ray observations. This results in no contemporaneous observation pairs at TeV/Hα around orbital phases 0.3–0.6, which correspond to the first maximum and the dip in the folded light curve. In the case of the contemporaneous observation pairs at X-ray/Hα, there is no coverage at orbital phases 0.0–0.2 and 0.5–0.6. This imperfect coverage of orbital phases might also obscure possible correlations.

Zoom In Zoom Out Reset image size

The data sets are divided into four different bins in orbital phase for the spectral analysis: the two higher flux states with orbital phases 0.2–0.4 and 0.6–0.8, the low state at orbital phases 0.4–0.6, and the medium (plateau) state for the remaining phases 0.8–0.2. This approach ignores the variability of the source on timescales shorter than these bins, as well as possible orbit-to-orbit variability, but allows a sufficient number of excess events to be collected for high-statistics spectra.

Results from the spectral analysis are summarized in Figure 7 and Table 4. The differential gamma-ray energy spectra in the four phase ranges are consistent within statistical uncertainties with simple power laws ${dN}/{dE}={N}_{0}{(E/(0.5\,\mathrm{TeV}))}^{-{\rm{\Gamma }}}$ with the exception of the long-exposure data set obtained with VERITAS for phase bin 0.2–0.4. Here, a spectral fit assuming a power law with exponential cutoff (${dN}/{dE}={N}_{0}{(E/(0.5\,\mathrm{TeV}))}^{-{\rm{\Gamma }}}{e}^{-E/{E}_{c}}$) with E_c = 1.75 ± 0.38 TeV provides a significantly better reduced χ² than a power-law fit without cutoff. Shorter exposure times at other phase bins of H.E.S.S. and MAGIC observations prevent us from measuring or refuting similar changes in the spectral shapes with those data sets. All measurements of the differential energy spectra by H.E.S.S., MAGIC, and VERITAS are compatible, within statistical and systematic uncertainties. No variability of the spectral index Γ is observed among spectra in different orbital phase bins, considering the statistical uncertainties. The spectral energy distributions for the four different phase bins are shown in Figure 8, including the X-ray spectra, which have been obtained in the same way by averaging over the phase bins discussed above. HESS J0632 + 057 is only weakly detected in the energy range covered by the Fermi-LAT, therefore only a phase-averaged spectrum is shown, taken from Li et al. (2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4.1. Detailed Examination of Orbits 9, 10, 16, and 17

Four orbits in the data sets presented provide reasonably good X-ray and gamma-ray coverage around orbital phase values of 0.2–0.45, which includes the first maximum and subsequent flux decay. Figure 9 shows the light curves for the four orbits with gamma-ray and X-ray coverage: 9, 10, 16, and 17. The shapes of the light curves of these four orbits differ notably and suggest the existence of orbit-to-orbit variability. The short-timescale variability on timescales of 20 days and less has been seen at X-ray energies (see also Figure 2). The observed differences in fluxes between orbits might be a biased view due to observations taking place at slightly different orbital phases. Observations during orbits 9 (2011 Jan) and 17 (2018 Jan) reveal gamma-ray and X-ray emission in bright states, comparable to the highest ever observed from this binary. The X-ray flux in orbits 10 and 16 is a factor of ∼1.5 lower than in orbits 9 and 17 at orbital phases of 0.35. The maximum gamma-ray flux also appears to be lower, although the low cadence of observations in this energy range does not allow a firm statement to be made.

Zoom In Zoom Out Reset image size

The bright state in orbit 9 has been reported in previous work by the MAGIC and VERITAS collaborations (Aleksić et al. 2012; Aliu et al. 2014). The gamma-ray flux for orbit 17 reached a level of (7.0 ± 1.3) × 10⁻¹² photons cm⁻² s⁻¹ above 350 GeV (corresponding to ≈7% of the flux of the Crab Nebula above the same energy). In this same orbit, HESS J0632 + 057 was detected on three individual nights within one week of observations with >7σ statistical significance per night (see Table 5).

Table 5. Summary of Observations and Spectral Analyses at Gamma-Ray Energies for Orbits 9 and 17

MJD Range	Exposure	Significance a	Flux Normalization	Photon Index	χ2/N
	(hr)	(σ)	Constant at 0.5 TeV
		(pre-trial)	[cm−2 s−1 TeV−1]
VERITAS observations of orbit 9 (2011 Jan)
55,585–55,600	11.3	6.1	(3.8 ± 0.9) × 10−12	2.71 ± 0.30	2.7/3
55,600–55,603	8.7	11.6	(8.3 ± 1.2) × 10−12	2.69 ± 0.17	4.5/4
MAGIC observations of orbit 9 (2011 Jan)
55,585–55,600	4.5	5.9	(4.9 ± 0.9) × 10−12	2.48 ± 0.20	4.2/3
VERITAS observations of orbit 17 (2018 Jan)
58,136	1.8	7.5	(11.9 ± 2.6) × 10−12	2.22 ± 0.24	3.3/4
58,141	2.5	7.7	(9.2 ± 2.2) × 10−12	2.20 ± 0.24	2.0/4
58,142	3.6	7.8	(11.6 ± 2.4) × 10−12	2.56 ± 0.28	1.5/3
58,143	3.1	5.8	(7.3 ± 2.0) × 10−12	2.15 ± 0.22	3.2/2

Note. The table lists the results of the power-law fits to the differential energy spectra obtained for VERITAS and MAGIC. Errors are 1σ statistical errors only. The systematic errors are 11%–20% on the flux constant and 0.10–0.15 on the spectral index.

^a Significances are calculated following the maximum likelihood ratio test method proposed in Li & Ma (1983).

Download table as:  ASCIITypeset image

The determination of variability timescales is limited by the cadence of the observations. For orbit 17, the flux changes from the highest state on MJD 58,136 to a flux below the detection limit on MJD 58,153 (flux upper limit for MJD 58,153: F_99%CL(>350 GeV) < 2.6 × 10⁻¹² photons cm⁻² s⁻¹), indicating a flux decay time shorter than 17 days. A similar timescale of roughly 20 days or less, again limited by the cadence and detection statistics of the observations, is observed for orbit 9.

The spectral energy distributions for orbits 9 and 17 are shown in Figure 10 and results for spectral fits to the gamma-ray data, assuming power-law functions, are given in Table 5. Despite the dramatic changes of the overall flux levels, no evidence for variability of the spectral index is detected, within statistical errors, for all six periods of gamma-ray observations.

Zoom In Zoom Out Reset image size

HESS J0632 + 057 has been observed at VHE with irregular coverage across 18 orbits, spanning an observational period from 2004 to 2019. These observations result in a statistically significant detection of nonthermal emission at all orbital phases, except phases 0.4–0.6. The phase-averaged luminosity of HESS J0632 + 057 above 1 TeV is ≈10³² erg s⁻¹ (assuming a distance of 1.4 kpc), making it one of the faintest gamma-ray binaries known to date (see Dubus 2013). The VHE gamma-ray emission extends beyond several TeV (see Figure 7), indicating the existence of radiating particles with at least similar energies. The energy range for the detected emission and the spectral shapes of individual spectra are, within statistical uncertainties, independent of the orbital phase of the measurement. The spectra follow approximately power-law shapes with spectral indices in the range 2.3–2.6 (Table 4), in line with generic expectations from diffusive shock acceleration.

In general, and independent of the nature of the compact object (black hole or neutron star), the VHE emission can be produced through interactions of accelerated hadrons or leptons with surrounding matter, magnetic, or radiation fields. The literature (e.g., Khangulyan et al. 2007) provides detailed discussions of mechanisms for acceleration and losses (e.g., synchrotron radiation, anisotropic inverse Compton scattering, or adiabatic losses) as functions of density and magnetic field strengths in binary systems. The discussion of the nonthermal processes at work in HESS J0632 + 057, including particle acceleration and X- and gamma-ray emission and absorption processes, is limited by the uncertainties present in some of the fundamental characteristics of the system. In particular, the unknown nature of the compact object allows different scenarios for the required power source for particle acceleration to be considered (e.g., the interaction of the stellar wind with an accretion-driven relativistic jet or with a pulsar wind). Equally important are the difficulties in determining the orbital geometry, given that two discrepant solutions with large uncertainties are found in the literature (Figure 1). Further uncertainties include the size and orientation of the equatorial disk of the Be star, the mass-loss rates and velocities of both equatorial and isotropic stellar winds, the properties of the accreting object (in the case of a black hole) or the pulsar wind (in the case of a neutron star), and the magnetic field configuration in the acceleration region.

The 15 yr of H.E.S.S., MAGIC, and VERITAS observations of HESS J0632 + 057, combined with measurements at optical and X-ray energies, suggest multiple patterns of variability on different timescales in the binary system. The most prominent modulation, interpreted as the orbital period, is consistently seen, within statistical uncertainties, in X-ray and gamma-ray measurements (P_X =317.3 ± 0.7(stat) ± 1.5(sys); P_γ = 316.7 ± 4.4(stat) ± 2.5(sys), see Section 3.1). In both bands, the orbital-phase-folded light curves (Figure 3) show two maxima at phases Φ ≈ 0.3 and Φ ≈ 0.6–0.8, separated by a deep flux minimum at orbital phase Φ ≈ 0.4. In the VHE band, the first maximum is more pronounced and exhibits, at least for some orbits, a sharp peak with a duration of a few days. The locations of the extrema in the X-ray and gamma-ray light curves along the binary orbit are not known, due to the large uncertainties in the orbital solutions, as discussed above. Assuming orbital parameters as presented in Casares et al. (2012) (Figure 1, top), the first maximum and the second maximum occur before and after apastron, where the environment around the compact object could likely be the least disturbed by the winds of the massive star. For the orbital solution proposed by Moritani et al. (2018), instead, the first maximum is after apastron and the second maximum shortly after periastron. A direct comparison with other gamma-ray binaries is affected by these uncertainties. However, it is worth pointing out that the two other long-period systems known (PSR B1259-63/LS 2883 and PSR J2032 + 4127/MT91 213) show their emission maxima close to or around the periastron passages.

The analysis of a total of 59 X-ray–VHE flux pairs reveals a strong correlation between these two energy bands (Figure 5, top). No time lag is observed (Figure 5, bottom), although the grouping of the gamma-ray observations (to gain statistics) prevents us from probing time lags τ ≤ 10 days. The X-ray–VHE correlation suggests that emission in both bands is produced by the same population of accelerated particles. One-zone leptonic models can describe the X-ray to TeV spectral distribution and the correlation between these two energy bands in HESS J0632 + 057, as seen in Hinton et al. (2009), Aleksić et al. (2012), and Aliu et al. (2014). In this scenario, a population of electrons produces X-rays via synchrotron emission and VHE gamma-rays through inverse Compton scattering off thermal photons from the companion Be star. It has been shown in Archer et al. (2020) that a one-zone leptonic model can describe well a small set of contemporaneous SED data in hard X-rays and VHE gamma-rays observed with NuSTAR and with VERITAS, respectively.

One of the main difficulties in modeling HESS J0632 +057's SED data lies in accounting for the observed orbit-to-orbit variability. Detailed hydrodynamical simulations of HESS J0632 + 057 by Bosch-Ramon et al. (2017) and Barkov & Bosch-Ramon (2018) propose a trapping of parts of the shocked pulsar wind by the massive star as a cause of orbit-to-orbit variability. These authors suggest that the peaks of the X-ray and VHE gamma-ray emission are a consequence of the accumulation of the hot shocked plasma injected by the shocked pulsar wind and the consequent disruption of the spiral arm in the periastron–apastron direction. It is proposed that the observed drop in the X-ray and gamma-ray fluxes is caused by the disruption of the stellar wind by the pulsar wind. In the simpler model of of Archer et al. (2020), this variability is assumed to be due to fluctuations of the electron density in the emission zone, which is treated as a free parameter of the fit. This kind of approach can only be applied to contemporaneous X-ray and VHE gamma-ray observations and it is only possible because of the large overlap between the electron energies responsible for the X-rays and VHE gamma-rays in this leptonic scenario.

For the first time, contemporaneous X-ray and gamma-ray observations have been presented for orbital phases 0.2–0.45 (around first maximum and minimum) for four different orbits (Figures 9 and 10). The measurements suggest that flux levels are different between these periods, indicating the existence of orbit-to-orbit variability in the emission pattern. No evidence for a superorbital period, as observed, e.g., in the gamma-ray binary LS I +61° 303 (Ahnen et al. 2016), is found. On shorter integration times, strong variability of the VHE flux on timescales of a few days was seen during observations in 2011 (orbit 9) and 2018 (orbit 17) at orbital phase ∼0.35 (prior to the local flux minimum; see Figure 9). Similar day-scale variability has been detected during VERITAS observations of the apastron passage of the compact object in the LS I +61° 303 binary system (Archambault et al. 2016).

In the GeV band, the weak flux of HESS J0632 + 057 hinders the detection of orbital variability with the Fermi-LAT (Li et al. 2017). No variability pattern is seen in Hα equivalent width, whereas the radial velocity and the Hα FWHM increase in the orbital phase range 0.3–0.6 and decrease between 0.7 and 0.3. The data demonstrate (with low statistical significance) a somewhat higher scatter at ϕ > 0.5, which may indicate a perturbed state of the stellar disk at these orbital phases. The correlation coefficients obtained for the different Hα parameters (EW, FWHM, and centroid velocity) and X-ray or VHE emission do not show significant values for the observed time lags (Figure 6). Optical spectroscopic measurements of MWC 148 during and after the time of the high state around 2018 January 25 with the RCC telescope at Rozhen NAO, Bulgaria, indicate changes in the structure or size of the equatorial disk of the Be star (Stoyanov et al. 2018), although comparable measurements have been carried out only during different orbital phases (2018 January 25: phase 0.36 and 2014 January 23: phase 0.71). Measurements of equivalent width by the Southern African Large Telescope (SALT Transient Program, L. Townsend 2020 personal communication) confirm this: pointings with a measured ∣EW∣ of >40 Å in 2018 February indicate a larger circumstellar disk than at the same phase in earlier cycles (e.g., ∣EW∣ of 2017 March 28 was ≈30 Å; see also measurements by Casares et al. 2012; Moritani et al. 2015).

Similarities of the emission pattern observed for HESS J0632 + 057 with those of other gamma-ray binaries are notable. Similar correlations between X-rays and VHE emission have been found during an outburst of LS I +61° 303 (Anderhub et al. 2009), but not consistently across different orbits (Acciari et al. 2011). However, many short-period systems like LS I +61° 303 are characterized by single-maximum light curves (Dubus 2013), while long-period gamma-ray binaries (e.g., PSR B1259-63/LS 2883, PSR J2032 + 4127/MT91 213) show a similar double-peaked structure in X- and gamma-ray emission (see, e.g., Dubus 2013; Abeysekara et al. 2018; H.E.S.S. Collaboration et al. 2020) as observed in HESS J0632 + 057. Moritani et al. (2018) and Malyshev et al. (2019) compare the emission patterns of HESS J0632 + 057 with those of PSR B1259-63/LS 2883 and relate the double-peak structure of the orbital-phase-folded light curves in the keV–VHE bands in both systems to the inclined disk geometry and the change in environmental parameters during the first and second crossings of the circumstellar disk by the compact object. The higher ambient density and increased magnetic field strengths in the stellar equatorial wind zone lead to enhanced particle acceleration via wind–wind interaction and possibly increased radiation efficiency due to an increased energy loss timescale for the electron population. Close to periastron, when the compact object leaves the disk, acceleration and radiation efficiency decrease, leading to the dip in the keV and VHE light curves. An & Romani (2017) obtain a quantitative description of the orbital variability pattern in X-rays of the binary 1FGL J1018.6-5856, which shows a similar double-peaked structure to that observed in HESS J0632 + 057. The model used a simple leptonic model based on intrabinary shock emission as discussed in the previous section, but with the following modifications: two electron populations are assumed, where one moves slowly in the shock and the other is accelerated along the shock flow; beaming effects, as proposed by Dubus et al. (2010), are added to the accelerated component in the cone-shaped shock front. Applying this model to HESS J0632 + 057, these components would lead to the enhanced flux state around orbital phases of 0.7–0.9 (when the flow direction and the line of sight align). The stellar disk of the Be star in HESS J0632 + 057 is the main difference compared to the diskless O-star in 1FGL J1018.6-5856 and might lead to enhanced absorption of the X-ray and gamma-ray emission. Additional changes in the magnetic fields at the intrabinary shock position due to the eccentric orbit and the resulting change in acceleration efficiency might be important to provide favourable conditions for the high emission state around phases 0.3–0.4 (Tokayer et al. 2021).

An alternative scenario to explain the variability pattern of HESS J0632 + 057 could be provided by flip-flop models. Torres et al. (2012) and Papitto et al. (2012) proposed this scenario to explain the emission from LS I +61° 303. In this scenario, the compact object is a pulsar. Different authors support the pulsar nature of the compact object in HESS J0632 + 057 (see Mochol & Kirk 2014; Moritani et al. 2015), as might be the case for the binaries LS I +61° 303 and LS 5039 (Yoneda et al. 2020), or as in the confirmed pulsar binaries PSR B1259-63/LS 2883 (Johnston et al. 1992; Kirk et al. 1999) and PSR J2032 + 4127/MT91 213 (Abdo et al. 2009; Camilo et al. 2009). In the model proposed by Torres et al. (2012) for LS I +61° 303, the pulsar switches from a propeller regime at periastron to an ejector regime at apastron, which is triggered by changes in the mass-loss rate of the companion star. Hence, for periods when the disk is at its largest, the propeller regime could be active also at apastron, reducing or impeding VHE emission and explaining the long-term variability in the binary. It should be noted that a recent Astronomer's Telegram (Weng et al. 2021) reported on the detection of a periodicity from LS I +61° 303 at a period 269.196 ms. This period is almost exactly that marked as the preferred one for the flip-flop model to work (see Figure 4 of Papitto et al. 2012). In this model, optical and VHE emissions are predicted to be anticorrelated. Moritani et al. (2015) argued that a similar mechanism could work in HESS J0632 + 057, although the orbital phases in which the gamma-ray emission occurs are different due to the geometry of the system. In the case of HESS J0632 + 057, the strong gas pressure will overcome the pulsar-wind pressure at periastron, suppressing VHE emission (assuming the ephemerides of Casares et al. 2012). The second minimum will take place at apastron, where the photon field and magnetic field density are low. In the case of LS I +61° 303, Ahnen et al. (2016) found the VHE periodicity as predicted, but no hints for (anti)correlation between optical and VHE were observed. The latter could be the result of the strong, short-timescale intraday variation displayed by the Hα fluxes. In HESS J0632 + 057, a similar situation can be at work, since both the contemporaneity of the observations (5 day time span) and the integration times are large, even more than for the case of LS I +61° 303, and were already problematic. The lack of (anti)correlation between optical and high-energy emission (X- and gamma-rays) observed in HESS J0632 + 057 can also be explained by the difference in the integration times required at these frequencies, which is in agreement with the results found by Ahnen et al. (2016) for the case of LS I +61° 303.

Regardless of the assumed model and the nature of the compact object, two additional processes can affect the orbital modulation of VHE emission, namely gamma–gamma absorption processes and the anisotropic nature of inverse Compton radiation. Close to periastron, the increased soft photon density can produce severe gamma–gamma absorption losses, which are potentially able to explain dips in the light curves of several gamma-ray binary systems (see, e.g., Böttcher & Dermer 2005; Dermer & Böttcher 2006; Sushch & van Soelen 2017). Anisotropic inverse Compton effects are important for orientations close to edge-on, as suggested by Bednarek (2006). In the context of LS 5039, Khangulyan et al. (2008) and Dubus et al. (2008) have shown that a model taking into account both effects can reproduce well the observed orbital light curve, producing somewhat different results to the spectrum in the isotropic inverse Compton case at inferior and superior conjunctions. The observed cutoff in the energy spectrum for the orbital phase bin 0.2–0.4 at E_c = 1.75 ± 0.38 TeV can have multiple, possibly correlated, origins: limiting properties of the shock determining the acceleration efficiency and the maximum energy of accelerated particles; photon fields determining the absorption of multi-TeV gamma-rays; or the effective temperature distribution of the low-energy target photons (e.g., due to substantial contribution from the stellar disk), which impact the energy of the gamma-rays obtained from inverse Compton processes.

An additional channel of information that can allow the proposed models to be discriminated against each other is the observation of flux variability patterns and outbursts at shorter timescales than the orbital period, as shown in Figure 9. Such variability suggests a rapid change in the environmental properties leading to particle acceleration, or in the surrounding medium impacting the VHE emission, possibly because of a peculiar geometry in the system (e.g., beaming, as proposed for PSR B1259-63/LS 2883 in, e.g., Khangulyan et al. 2011; Chernyakova et al. 2020). In terms of the discussed models, however, short-term variability can be accommodated by several theoretical explanations. In the flip-flop scenario these timescales could be interpreted as the timescales on which the pulsar wind is quenched due to increased circumstellar disk density. In the PSR B1259-63/LS 2883-like model, the same timescales would correspond to the times at which the compact object intersects the disk's region of influence and exits from it close to periastron. Additionally, stellar winds of massive stars can exhibit strong density variations. Flux variability can be connected to the clumpiness of the equatorial or isotropic stellar wind, which can lead to keV–VHE flux changes; see, e.g., Paredes-Fortuny et al. (2015a) and Archambault et al. (2016). In this case the variability time can allow us to estimate the characteristic size of the corresponding high-density structure. Deeper, simultaneous multiwavelength (keV–GeV–VHE) observations of HESS J0632 + 057 with current or future facilities can help to clarify these points. Such observations may allow one to find the shortest timescales of variability in each band and more firmly establish the similarities and differences in the variability patterns. Additional future Hα observations used as a tracer for the radius of the circumstellar disk are also helpful, since the data set shown in this paper is rather sparse.