Probing the Properties of the Pulsar Wind in the Gamma-Ray Binary HESS J0632+057 with NuSTAR and VERITAS Observations

Gamma-ray binaries are systems composed of a massive star and a compact object in which the nonthermal luminosity is dominated by the gamma-ray emission (>1 MeV) (Dubus 2013). Objects of this class are extremely rare, and currently only a handful of sources have been unambiguously identified as gamma-ray binaries: PSR B1259-63, LS 5039, LS I+61 303, HESS J0632+057, 1FGL J1018.6–5856, LMC P3, PSR J2032+4127 (Paredes & Bordas 2019), and 3FGL J1405.4-6119 (Corbet et al. 2019). All of these systems host compact objects, either neutron stars or black holes, orbiting a massive star of O or B type. Their orbital periods vary substantially, spanning from 3.9 days (LS 5039) to ∼50 yr (PSR J2032+4127). Only PSR B1259-63 and PSR J2032+4127 contain known radio pulsars, while the nature of the compact object in the remaining systems is still unknown.

Two distinct scenarios are commonly invoked to explain the origin of the nonthermal emission in gamma-ray binaries. In objects known as microquasars, the accretion of stellar material onto the compact object produces relativistic jets where particles are accelerated to high energies and can radiate in the gamma-ray band (Mirabel & Rodríguez 1998). Where the system hosts a pulsar, the spin-down power of a young pulsar is transferred to high-energy particles through the acceleration of electron pairs from the pulsar wind in the termination shock formed by the interaction of the pulsar and the stellar winds (Tavani & Arons 1997).

HESS J0632+057 was first detected as an unidentified point-like source during H.E.S.S. observations of the Monoceros region (Aharonian et al. 2007). Its binary nature in association with the Be star MWC 148 was suggested at first and supported later by XMM-Newton observations in X-rays (Hinton et al. 2009). Subsequent observations by VERITAS did not lead to a detection (Acciari et al. 2009), indicating a substantial flux variability, which is characteristic of gamma-ray binaries. Since then, extensive observations have been performed in soft X-rays (Falcone et al. 2010) and the TeV gamma-ray band (Aleksić et al. 2012; Aliu et al. 2014). Both X-ray and gamma-ray light curves present two clear periodic peaks around ϕ ≈ 0.3–0.4 and ϕ ≈ 0.6–0.8 when folded to an orbital period of around 315–320 days (Aliu et al. 2014). The system is uncommonly faint in the GeV band compared with other gamma-ray binaries, having been detected only recently in Fermi-LAT data (Li et al. 2017).

Aragona et al. (2010) estimated the distance of the system to be 1.1–1.7 kpc through optical spectroscopy of the companion star MWC 148. The orbital period of ${P}_{\mathrm{orb}}$ = 321 days was first derived by Bongiorno et al. (2011) using X-ray light-curve data, and was later refined to be ${P}_{\mathrm{orb}}={315}_{-4}^{+6}$ by Aliu et al. (2014). Two distinct orbital solutions have been proposed by Casares et al. (2012) and, more recently, by Moritani et al. (2018). Both studies follow similar methodologies based on optical spectroscopic data; however, Moritani et al. (2018) uses a more extensive data set, which results in a set of orbital parameters entirely different from those by Casares et al. (2012).

The nature of the compact object in HESS J0632+057 is still uncertain, as the implications of the two orbital solutions point in the opposite directions. Moritani et al. (2018) suggest that the mass of the compact object must be consistent with the pulsar scenario (${M}_{\mathrm{co}}\lt 2.5$ ${M}_{\odot }$) unless the inclination of the system is extremely small (<3°). On the other hand, the solution by Casares et al. (2012) favors the black hole scenario with ${M}_{\mathrm{co}}\gt 2.1\,{M}_{\odot }$ (Zamanov et al. 2017), being only marginally compatible with the pulsar scenario. In addition to the mass functions, indirect evidence of the pulsar scenario can only be claimed by comparing the observational features of HESS J0632+057 to other known sources, especially PSR B1259-63, in which the compact object is a known pulsar.

In this paper, we present simultaneous observations in X-ray by NuSTAR and TeV gamma-ray by VERITAS during 2017 November and December. These observations correspond to ϕ ≈ 0.22 and 0.30, respectively, when the X-ray and TeV fluxes are rising toward the first flux peak (Aliu et al. 2014). The pulsar scenario is probed by fitting the spectral energy distribution (SED) to a model in which the nonthermal emission is assumed to be produced by electrons from the pulsar wind that are accelerated at the termination shock formed by the collision of the stellar and pulsar wind. The X-ray and TeV gamma-ray spectra are assumed to be produced through synchrotron radiation and inverse Compton scattering (ICS) of stellar photons, respectively. Because the distance between the Be star and the compact object at the time of the observations is substantially larger than the estimated size of the circumstellar disk, the model relies on a few assumptions about the properties of the polar stellar wind, neglecting the disk component. As the result of the model fitting, we obtain a relation between the pulsar spin-down luminosity (${L}_{\mathrm{sd}}$) and the pulsar-wind magnetization (σ) that is consistent with theoretical expectations. This result shows that our data can be satisfactorily described within the pulsar hypothesis.

The observations and data analysis are first described in Sections 2 and 3. In Section 4 we describe the system parameters, including both orbital solutions, and define the sets of parameters adopted in this work. The model based on the pulsar-wind scenario is described in Section 5, and the results of the SED model fitting in Section 6. Finally, we summarize and discuss the results in Section 7.

2.1. NuSTAR

2.2. VERITAS

3.1. NuSTAR Spectral Analysis

Source photons were extracted from a r = 30'' circle centered on the source position at R.A. = 982472 and decl. = 58015 (J2000), yielding a total of 2534/2468 net counts (FPMA and FPMB combined) for the 2017 November/December observations, respectively, for a net count rate of ∼0.02 cts s⁻¹ for each detector module. Background spectra were extracted from a rectangular source-free region on the same detector chip as the source. We used nuproducts to generate the response matrix and ancillary response files.

NuSTAR spectra were grouped to a minimum significance of 5σ in each bin and fitted with the XSPEC (v12.9.1) package (Arnaud 1996) using χ² statistics. We fit NuSTAR module A and B spectra jointly in the 3–30 keV energy band, above which the background dominates. Given the previously measured column density values (N_H ∼ (2.1–4.7) × 10²¹ cm⁻²) (Moritani et al. 2018), we find that the ISM absorption is negligible above 3 keV. Therefore, NuSTAR 3–30 keV spectra allow us to determine the intrinsic continuum spectral index independently, without degeneracy with N_H.

We fit a single power-law model to the spectra. For the November observation, the best-fit photon index was Γ = 1.77 ± 0.05 (χ/dof = 0.92 for 105 dof); the 3–30 keV flux for the unabsorbed power law is (2.42 ± 0.13) × 10⁻¹² erg cm⁻² s⁻¹, and the corresponding luminosity is (5.67 ± 0.30) × 10³² erg s⁻¹. The December spectrum is significantly harder, with photon index Γ = 1.56 ± 0.05 (χ/dof = 0.71 for 102 dof), while 3–30 keV flux is (2.45 ± 0.13) × 10⁻¹² erg cm⁻² s⁻¹, for a luminosity of (5.75 ± 0.30) × 10³² erg s⁻¹. All uncertainties were calculated at 1σ intervals. Luminosity values assume a distance of 1.4 kpc.

We find no significant spectral break or cutoff in the NuSTAR spectra. The NuSTAR SED derived from the spectral analysis is shown in Figure 1 (upper panels).

Zoom In Zoom Out Reset image size

3.2. NuSTAR Timing Analysis

3.3. VERITAS Analysis

Accurate knowledge of the system's geometry and the properties of the companion star are critical ingredients for any attempt at modeling gamma-ray binaries. In this section, we summarize the available orbital solutions and the parameters related to the companion Be star, including the properties of the stellar wind. We also describe the sets of parameters that will be assumed for the model fitting presented in the next sections.

The properties of the companion Be star MWC 148 (=HD 259440) were first derived through optical spectroscopy by Aragona et al. (2010). The effective temperature T was found to be 27.5–30.0 kK, the mass 13.2–19.0 ${M}_{\odot }$ and the radius 6.0–9.6 ${R}_{\odot }$. The distance of the system was also derived to be 1.1–1.7 kpc. Based on these ranges of values, we assume ${T}_{\mathrm{Be}}=30$ kK, ${R}_{\mathrm{Be}}=7.8$ ${R}_{\odot }$ and $d=1.4\,\mathrm{kpc}$. The mass of the star ${M}_{\mathrm{Be}}$, however, will be allowed to vary within the derived range according to the mass function f(M) obtained in the orbital solutions and the system's inclination i.

The first orbital solution was obtained by Casares et al. (2012) through spectroscopic studies of Hα emission lines. Recently, however, Moritani et al. (2018) proposed a distinct solution, obtained with the same methodology but with a larger data set. For the present study, both solutions will be considered and two sets of parameters will be defined. The orbital period (${P}_{\mathrm{orb}}$) of the system was assumed to be 321 days by Casares et al. (2012) and derived to be 313 days by Moritani et al. (2018).³⁰ Although ${P}_{\mathrm{orb}}$ is correlated to the more recent orbital parameters (Ho et al. 2017), we assume here that the impact of varying it by a few days is negligible and we will set ${P}_{\mathrm{orb}}=315$ days for both sets of parameters, following the results obtained by Aliu et al. (2014). To estimate the impact of the ${P}_{\mathrm{orb}}$ uncertainties on the results of our studies, we will consider an arbitrary ±2 days uncertainty range. The sets of parameters defined for our study are shown in Table 1 and the orbit of the compact object is illustrated in Figure 2. The uncertainties shown in Table 1 were accounted for to probe the robustness of the model fitting. The uncertainties on ${R}_{\mathrm{Be}}$ and ${T}_{\mathrm{Be}}$ are not relevant and will be neglected.

Zoom In Zoom Out Reset image size

Table 1.  System Parameters of HESS J0632+057 as Used in This Work

	Casares et al. (2012)	Moritani et al. (2018)
e	0.83 ± 0.08	0.64 ± 0.29
ω (°)	129 ± 17	271 ± 29
f (${M}_{\odot }$)	0.01	0.0024
i (°)	69.5 ± 10.5	37 ± 5
a2 (au)	${3.90}_{-0.22}^{+0.13}$	${2.13}_{-0.17}^{+0.14}$
${P}_{\mathrm{orb}}$ (day)	315 ± 2
${M}_{\mathrm{psr}}$ (${M}_{\odot }$)	1.4
${M}_{\mathrm{Be}}$ (${M}_{\odot }$)	13.2–19.0
${R}_{\mathrm{Be}}$ (${R}_{\odot }$)	7.8
${T}_{\mathrm{Be}}$ (K)	30000
${v}_{{\rm{w}}}$ (km s−1)	1500
${\dot{M}}_{{\rm{w}}}$ (${M}_{\odot }{\mathrm{yr}}^{-1}$)	${10}^{8.5\pm 0.5}$

Note. See the text for the description and references. The MJD of ϕ = 0 is 54857.5 as defined by Falcone et al. (2010).

Download table as:  ASCIITypeset image

The compact object is assumed to be a pulsar with ${M}_{\mathrm{psr}}=1.4$ ${M}_{\odot }$ in the model considered here. With this assumption, the mass function f(M) from a given orbital solution constrains the relation between ${M}_{\mathrm{Be}}$ and the inclination i. Thus, the allowed range of i can be defined by imposing that ${M}_{\mathrm{Be}}$ is consistent with the range derived by Aragona et al. (2010). ${M}_{\mathrm{Be}}\,\mathrm{versus}\,i$ is shown in Figure 3, from which we find that (a) $32^\circ \lt i\lt 42^\circ$ for the Moritani et al. (2018) solution and (b) there is no value of i allowed for the Casares et al. (2012) solution assuming the nominal value of f(M) = 0.06, as can be seen by the solid red line. To define a set of orbital parameters based on the Casares et al. (2012) solution, we redefine f(M) to its lower bound within the quoted uncertainties, given by f(M) = 0.06–0.05 = 0.01. By doing this, we find that $i\gt 59^\circ$ becomes allowed (see the dashed red line in Figure 3). The upper bound of i will be set to $80^\circ$ following the argument given by Casares et al. (2012) based on the lack of observed shell lines, which would have indicated obscuration of the star by the circumstellar disk.

Zoom In Zoom Out Reset image size

In the orbital parameters defined for our study, the inclination will be taken as the center of the allowed range, with the full range being used to probe the robustness of the results. As ${M}_{\mathrm{Be}}$ is linked to the inclination through the mass function, i and ${M}_{\mathrm{Be}}$ are strongly correlated. A similar relation occurs with the semimajor axis of the pulsar orbit (a₂), as the orbital solutions directly provide a₁ sin i, where a₁/a₂ = ${M}_{\mathrm{psr}}$/${M}_{\mathrm{Be}}$. The uncertainties of a₂ indicated in Table 1 are due to the effect of varying i within its allowed range.

In shocked wind models, the properties of the stellar winds are relevant ingredients because the location of the intrabinary shock is determined by the dynamical pressure balance between the stellar and pulsar winds. For Be stars, the stellar wind is commonly described as being composed of a fast low-density polar wind and a slow dense equatorial disk wind (Waters et al. 1988). The disk properties of MWC 148 were studied by Moritani et al. (2015), Zamanov et al. (2016) through optical spectroscopy. It was found that the disk size, estimated as the radius of the Hα emitting region, is 0.85–1.4 au. The average disk radius of 1.12 au is indicated by dotted lines in Figure 2, assuming that the disk and the orbit of the compact object are co-planar.

For the sets of observations reported here, both orbital solutions imply that the distance between the Be star and the pulsar is substantially larger than the disk size. Therefore, the influence of the disk on the shock formation in the colliding wind scenario is expected to be relevant only if the shock is located relatively close to the companion star. We will show, by the results of the model fitting, that the pulsar spin-down luminosity ${L}_{\mathrm{sd}}$ required for the latter condition is not supported by our observations. Therefore, we will neglect the interaction with the disk in our model and only the polar-wind parameters will be defined.

The velocity profile of the polar wind is approximately described by ${v}_{{\rm{w}}}(r)={v}_{0}+({v}_{\infty }-{v}_{0})(1-{R}_{\mathrm{Be}}/r)$ (Waters et al. 1988), where typically v₀ ∼ 20 km s⁻¹ and ${v}_{\infty }\sim 1000\mbox{--}2000$ km s⁻¹. Because $r\gg {R}_{\mathrm{Be}}$ in the context of our study, we use the approximation ${v}_{{\rm{w}}}={v}_{\infty }$, assuming ${v}_{{\rm{w}}}=1500$ km s⁻¹. The mass-loss rate ${\dot{M}}_{{\rm{w}}}$ of the wind in Be stars is only poorly constrained and it is usually assumed to be in the range 10⁻⁹–10⁻⁸ ${M}_{\odot }{\mathrm{yr}}^{-1}$ based on Snow (1981); Waters et al. (1988). The impact of the ${\dot{M}}_{{\rm{w}}}$ uncertainties on the results of the model fitting is substantial, because the relative location of the intrabinary shock is strongly affected by this parameter. Thus, we will adopt ${\dot{M}}_{{\rm{w}}}={10}^{-8.5}$ ${M}_{\odot }{\mathrm{yr}}^{-1}$ as a reference value and the range 10⁻⁹–10⁻⁸ ${M}_{\odot }{\mathrm{yr}}^{-1}$ will be used to provide an estimation of its uncertainties.

Our model is based on the assumption that the compact object is a pulsar of mass ${M}_{\mathrm{psr}}=1.4$ ${M}_{\odot }$. The collision between the pulsar and the stellar wind creates a termination shock which is assumed to accelerate cold electron pairs from the pulsar wind to high energies. Hence, nonthermal radiation emitted by the accelerated electrons located around the apex of the shock is believed to produce the observed X-ray and gamma-ray photons through synchrotron and ICS, respectively. The pulsar spin-down luminosity, ${L}_{\mathrm{sd}}$, is a central parameter of the model, as it drives the pulsar-wind pressure as well as the B-field strength. It is important to note that the following fitting approach is only meant to describe the observations presented in this paper with a minimum number of assumptions. The application of the same ideas to observations taken at other points in the orbit would require to take into account further processes, in particular, the effect of the circumstellar disk on the shock formation.

5.1. Shock Formation and the Pulsar-wind Magnetization

By imposing hydrodynamic balance between the pulsar and the stellar wind, the distance of the shock apex to the pulsar is given by (Harding & Gaisser 1990; Tavani & Arons 1997; Ball & Kirk 2000)

$\begin{eqnarray}&&{R}_{\mathrm{sh}}=\displaystyle \frac{\sqrt{\eta }}{1+\sqrt{\eta }}D,\end{eqnarray} \tag{ 1 }$

where

$\begin{eqnarray}&&\eta =\displaystyle \frac{{L}_{\mathrm{sd}}}{\dot{M}{v}_{{\rm{w}}}c},\end{eqnarray} \tag{ 2 }$

and D is the distance between both objects. $\dot{M}$ and v_w are the mass-loss rate and the wind velocity assumed for the polar component of the stellar wind. The circumstellar disk component of the wind is neglected here because the distance between the objects is much larger than the size of the disk as estimated by Moritani et al. (2015), Zamanov et al. (2016). The influence of the disk, if present, would move the shock position closer to the pulsar, changing the properties of the nonthermal emission. The same effect would be observed by assuming different values of $\dot{M}$, while still neglecting the disk component of the wind. Thus, in this model, the effect of disk interactions on the emission features are estimated by changing $\dot{M}$ within its uncertainties.

The particle acceleration, and consequently, nonthermal emission, is assumed to occur in a relatively small region around the shock apex, and therefore, the shock morphology will be neglected. This assumption is supported by the results of hydrodynamical simulations of the gamma-ray binary LS 5039 (Dubus et al. 2015).

The relativistic electron pairs from the pulsar wind are accelerated at the termination shock, forming the high-energy electron population, upstream to the shock that is responsible for the nonthermal radiation. Following Kennel & Coroniti (1984a, 1984b), the B-field upstream of the shock is given by

$\begin{eqnarray}&&B=\sqrt{\displaystyle \frac{{L}_{\mathrm{sd}}\sigma }{{R}_{\mathrm{sh}}c(1+\sigma )}\left(1+\displaystyle \frac{1}{{u}^{2}}\right)},\end{eqnarray} \tag{ 3 }$

where u is the radial four-velocity of the wind downstream from the shock, given by

$\begin{eqnarray}\begin{array}{rcl}{u}^{2} & = & \displaystyle \frac{8{\sigma }^{2}+10\sigma +1}{16(\sigma +1)}\\ & & +\displaystyle \frac{{\left[64{\sigma }^{2}{\left(\sigma +1\right)}^{2}+20\sigma (\sigma +1)+1\right]}^{1/2}}{16(\sigma +1)}.\end{array}\end{eqnarray} \tag{ 4 }$

The magnetization of the wind σ is commonly assumed to depend on the distance from the pulsar, as $\sigma \propto {R}_{\mathrm{sh}}^{-\alpha }$ with $\alpha \sim 0.5\mbox{--}2.0$ (Kong et al. 2012; Takata et al. 2017).

5.2. Energy Spectrum of the High-energy Electron Population

5.3. Radiative Processes

The population of high-energy electrons accelerated at the termination shock is assumed to radiate through synchrotron and ICS, producing the X-ray and gamma-ray photons observed at Earth, respectively. The B-field responsible for the synchrotron emission is calculated following Equation (3). The seed photon field for the ICS is composed of the thermal photons radiated by the companion star, with a density given by

$\begin{eqnarray}&&U=\displaystyle \frac{{\sigma }_{\mathrm{SB}}{T}_{\mathrm{Be}}^{4}{R}_{\mathrm{Be}}^{2}}{c{\left(D-{R}_{\mathrm{sh}}\right)}^{2}},\end{eqnarray} \tag{ 5 }$

where σ_SB is the Stefan–Boltzmann constant. We also account for the anisotropic nature of the ICS by calculating the photon scattering angle (${\theta }_{\mathrm{ICS}}$) from the geometry given by the orbital solutions.

Because of the relatively dense field provided by the stellar photons, the impact of pair-production absorption of gamma-rays (Gould & Schréder 1967) becomes relevant. For each evaluation of the model, we calculate the optical depth (τ_γγ) as a function of the gamma-ray photon energy, accounting for the corresponding position of the termination shock and geometry of the system (see e.g., Dubus 2006; Sierpowska-Bartosik & Torres 2008; Sushch & van Soelen 2017 for the details of the τ_γγ calculation). The gamma-ray spectrum expected to be observed is then attenuated by a factor of ${e}^{-{\tau }_{\gamma \gamma }}$. The impact of pair-production absorption is the strongest at the low end of the observed gamma-ray spectrum (∼0.2 TeV), and it is larger for the Casares et al. (2012) orbital solution because of the longer path of the gamma-ray photons within the stellar photon field.

5.4. Model Fitting

The model described above was fitted to the SED derived from both NuSTAR and VERITAS observations. In what follows, the two observation sets will be labeled by the lower indices 0 and 1, for the 2017 November and December observations, respectively. From the orbital solutions described in Section 4, we calculated the distance between the objects (D) and the ICS angle (${\theta }_{\mathrm{ICS}}$), which are shown in Figure 4 as a function of the system's phase. Assumptions are made regarding the stellar wind properties ${\dot{M}}_{{\rm{w}}}$ and ${v}_{{\rm{w}}}$. The impact of the uncertainties on ${\dot{M}}_{{\rm{w}}}$ will be estimated by repeating the procedure with different values.

Zoom In Zoom Out Reset image size

The remaining parameters of the model are the pulsar spin-down luminosity (${L}_{\mathrm{sd}}$), the pulsar-wind magnetization at termination shock for both periods (σ₀ and σ₁), and the parameters of the electron spectrum for both periods (${N}_{e,0}$, ${{\rm{\Gamma }}}_{0}$, ${N}_{e,1}$ and ${{\rm{\Gamma }}}_{1}$). The dependence of σ on ${R}_{\mathrm{sh}}$ is assumed to constrain σ₀ and σ₁ and reduce them to a single parameter. Here, we maintained σ₀ as a free parameter of the model and we calculated ${\sigma }_{1}={\sigma }_{0}{\left(\tfrac{{R}_{\mathrm{sh},0}}{{R}_{\mathrm{sh},1}}\right)}^{\alpha }$, with α = 1. As the difference between ${R}_{\mathrm{sh},0}$ and ${R}_{\mathrm{sh},1}$ is relatively small for our observations, the impact of the choice of α = 1 is minimal.

Regarding the electron spectrum, only ${N}_{e,0}$ and ${N}_{e,1}$ were treated as free parameters, while ${{\rm{\Gamma }}}_{0}$ and ${{\rm{\Gamma }}}_{1}$ were set to the values derived from the single power-law fit of the X-ray spectrum (see Section 3). The values are ${{\rm{\Gamma }}}_{0}=1.77\times 2\mbox{--}1=2.54$ and ${{\rm{\Gamma }}}_{1}=1.56\times 2\mbox{--}1=2.12$. In the end, the model contains four free parameters: ${L}_{\mathrm{sd}}$, σ₀, ${N}_{e,0}$ and ${N}_{e,1}$.

The model evaluation starts by inserting ${L}_{\mathrm{sd}}$ into Equations (1) and (2) to calculate ${R}_{\mathrm{sh},0}$ and ${R}_{\mathrm{sh},1}$ (using the assumptions of ${\dot{M}}_{{\rm{w}}}$ and ${v}_{{\rm{w}}}$ and D₀/D₁ given by the orbital solutions). With ${R}_{\mathrm{sh}}$, the photon density U of the ICS seed field and the optical depth τ_γγ are calculated. In the following, the values of the B-field at the shock, B₀ and B₁, are calculated by inserting the σ and ${R}_{\mathrm{sh}}$ values in Equation (3). Finally, with the electron spectrum determined by ${N}_{e,0}$ and ${N}_{e,1}$, the synchrotron and ICS spectra are calculated, using ${\theta }_{\mathrm{ICS}}$ from the orbital solutions for the ICS emission.

The synchrotron and ICS emission were computed by following Blumenthal & Gould (1970) using the Naima package (Zabalza 2015). The SED fit is performed by means of a χ² method. While ${L}_{\mathrm{sd}}$ and σ₀ are scanned over a predefined grid, ${N}_{e,0}$ and ${N}_{e,1}$ are fitted by minimizing χ² using Minuit framework (James & Roos 1975).

We show in Figure 5 the results of the fit for the ${L}_{\mathrm{sd}}$–σ₀ plane, indicating the 1 and 2σ regions for both orbital solutions. The best solutions are indicated by stars and the corresponding χ²/dof is 0.786 for both sets of orbital parameters. The dashed lines indicate the σ₀ that minimizes the χ² for each value of ${L}_{\mathrm{sd}}$. These curves will be labeled as best curves in what follows. The SED model–data comparisons are shown in Figure 6, in which the best solution of the fit is used to evaluate the model. To illustrate the SED throughout all the energy bands, the electron spectrum was assumed to be a power law, starting at E_min, with an exponential cutoff characterized by E_cut. Although different values of E_min and E_cut are illustrated in Figure 6, the model fitting does not depend on these parameters, as explained in Section 5.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The degeneracy observed between ${L}_{\mathrm{sd}}$ and σ₀ implies that neither of these parameters can be individually constrained with our approach. Nevertheless, our results demonstrate that the observations can be consistently described by a pulsar wind model and indicate the region of the ${L}_{\mathrm{sd}}$-σ₀ space that allows for it. To further characterize the allowed solutions, we show in Figure 7 the B-field, ${R}_{\mathrm{sh}}$ and τ_γγ for the November observation (B₀, ${R}_{\mathrm{sh},0}$ and ${\tau }_{\gamma \gamma ,0}$) as a function of ${L}_{\mathrm{sd}}$. The lines in Figure 7 correspond to the best curves shown in Figure 5.

Zoom In Zoom Out Reset image size

The estimated total luminosity of the full SED shown in Figure 6 is in the range 10³⁴–10³⁵ erg s⁻¹. Estimations of the efficiency of nonthermal radiation in gamma-ray binaries are highly uncertain. For instance, Fermi-LAT observations show that nearly all the ${L}_{\mathrm{sd}}$ is converted into nonthermal radiation in PSR B1259-63 around the periastron passage (Abdo et al. 2011). Sierpowska-Bartosik & Torres (2008), on the other hand, assume an efficiency of 0.01 for LS 5039. Therefore, the truncation of our results at ${L}_{\mathrm{sd}}={10}^{34}$ erg s⁻¹ represents a conservative lower limit.

The impact of the different system geometries provided by the two orbital solutions is minimal. The observations are equally well described by both sets of orbital parameters, with no indication of one of them being favored. While σ₀ for ${L}_{\mathrm{sd}}\lt {10}^{37}$ erg s⁻¹ is two to three times larger for the Casares et al. (2012) solution, the difference in B₀ amounts to only a factor of 0.6 at the highest ${L}_{\mathrm{sd}}$ regime.

At the high ${L}_{\mathrm{sd}}$ regime of the solutions, the upper limit within 1σ is ${L}_{{\rm{sd}}}\lt 7\times {10}^{37}$ erg s⁻¹ for both sets of orbital parameters, which is consistent with expectations for very young pulsars. In this regime, the relatively strong pulsar wind moves the termination shock closer to the companion star (see Figure 7, middle). For both orbital solutions, the termination shock is approximately halfway between the stars for the highest ${L}_{\mathrm{sd}}$ values, although the absolute distances are different. The proximity of the shocked region to the Be star implies a larger density for the ICS seed-photon field, and consequently, enhanced ICS emission, but also stronger pair-production absorption. As the ratio between ICS and synchrotron emission is constrained by the data, the B-field has to change accordingly to compensate for the ICS behavior, which is done by changing ${\sigma }_{0}$. In this regime, B₀ ≈ 0.25 and ≈0.5 G for the Casares et al. (2012) and Moritani et al. (2018) solutions, respectively.

As the ${L}_{\mathrm{sd}}$ of the solutions decreases, the termination shock moves closer to the pulsar, within a small fraction of the distance between the objects. At this regime, the B-field also decreases and reaches B₀ ≈ 0.06 and 0.07 G for the Casares et al. (2012) and Moritani et al. (2018) solutions, respectively.

The impact of the uncertainties on the system's parameters was evaluated by varying the assumed values and repeating the fitting procedure. We found that the most relevant uncertainty is that of the mass-loss rate of the stellar wind, ${\dot{M}}_{{\rm{w}}}$. The results of the fit for ${\dot{M}}_{{\rm{w}}}={10}^{-9.0},{10}^{-8.5}$ and 10^−8.0 ${M}_{\odot }{\mathrm{yr}}^{-1}$ are shown in Figure 8. We also show the results in Figure 9 by accounting for the uncertainties in the system's inclination. Note that as described in Section 4, the variations of the inclination also result in variations in the semimajor axis, as well as in the mass of the companion star. Apart from the ${\dot{M}}_{{\rm{w}}}$ and inclination, we also tested the robustness of the results considering the uncertainties in the eccentricity, ω, ${P}_{\mathrm{orb}}$, and α. In all cases the effect is similar to or slightly larger than that observed for the inclination. The uncertainty range considered for each of these parameters is indicated in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The impact of the poorly known properties of the stellar wind, represented by the mass-loss rate ${\dot{M}}_{{\rm{w}}}$, is the most relevant source of uncertainty on the constraining power of the pulsar-wind magnetization from our approach. The large range of ${\dot{M}}_{{\rm{w}}}$ assumed here is sufficient to incorporate any effect due to the presence of the circumstellar disk, which has been neglected by the model assumptions.

We presented the results of two sets of combined observations of the gamma-ray binary HESS J0632+057, by NuSTAR in the hard X-ray band, and by VERITAS, in the TeV gamma-ray band. These observations correspond to the rise of the first peak observed in the X-ray and TeV light curves, at ϕ ≈ 0.22 and 0.30. The data provided by both instruments allows us to derive the combined SED for these two periods. The spectral and timing analysis performed on the NuSTAR observations show that: (a) the spectra are well described by a single power-law model with a spectral hardening observed between the two observations (Γ going from 1.77 ± 0.05 to 1.56 ± 0.05), and (b) no evidence of time variability, red noise, or pulsation has been found, consistent with observations of other gamma-ray binaries.

The SED data is used to probe a model based on the pulsar-wind scenario, in which the nonthermal emission is produced by high-energy electrons accelerated at the termination shock formed by the interaction between pulsar and stellar wind. The applied model relies on a minimum number of assumptions which are sufficient for the description of the observations presented in this paper. The description of further observations would require to expand the number of assumptions and the complexity of the model.

The model fitting benefited from two characteristics of our observations. First, for both orbital solutions, the observations correspond to periods in which the distance between the Be star and the compact object is sufficiently large so that the influence of the circumstellar disk is negligible, reducing the number of assumptions required to describe the properties of the stellar wind. second, in systems like gamma-ray binaries, there is a large overlap between the energy range of the electrons responsible for producing hard X-rays and very-high-energy gamma-rays (through synchrotron and ICS, respectively), which allows us to simplify the assumptions about the shape of the electron spectrum to a single power-law description, with only two parameters.

The results of the model fitting show the regions of the ${L}_{\mathrm{sd}}$–σ plane that are allowed by our data. We find an upper limit of ${L}_{\mathrm{sd}}\lt 7\,\times \,{10}^{37}$ erg s⁻¹ within 1σ, independent of the orbital solution, which is consistent with expectations of young pulsars. The σ parameter is constrained to be 0.003–0.03 at the location of the shock. Constraints on σ are particularly important for understanding the physical process behind the transport of energy from the rotation powered pulsar to the surrounding medium, which is still a subject of intense discussions (Arons 2002; Kirk et al. 2009). While theoretical models predict that at the light cylinder the pulsar wind is dominated by Poynting energy (σ_L ≫ 1), observations of the Crab Nebula constrain σ at much larger distances to be kinetic particle dominated (σ_N ≪ 1). The transition between these two regimes is not well described within the current theoretical framework, originating the so-called "σ problem." In gamma-ray binaries, the pulsar-wind termination is typically located at intermediate distances between the light cylinder and the termination shock in pulsar-wind nebulae, which makes these systems interesting in this context.

In Figure 10 we compare the σ versus ${R}_{\mathrm{sh}}$ of the best curves obtained from our fitting with the following selection of results: (a) at the light cylinder (${R}_{\mathrm{sh}}\approx {10}^{8}$ cm), Kong et al. (2012) estimates theoretically that $\sigma \approx 1\,\times \,{10}^{5}$ and $8\,\times \,{10}^{3}$ for the Crab Pulsar and for PSR J1259–63, respectively, (b) Tavani & Arons (1997) assume σ = 0.02 at the termination shock to describe observations of PSR J1259–63 around the periastron, (c) Dubus et al. (2015) constrains σ ≈ 1 at the termination shock to describe observations of the system LS 5039 using a numerical hydrodynamical model and (d) the classical result by Kennel & Coroniti (1984b) that constrains σ ≈ 0.003 at ${R}_{\mathrm{sh}}\approx {10}^{17}$ cm for the Crab Nebula. The comparison indicates that all systems follow a similar trend of σ decreasing with ${R}_{\mathrm{sh}}$. The consistency observed between our results and other systems provides further support to the pulsar scenario of HESS J0632+057.

Zoom In Zoom Out Reset image size

This work used data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. We made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA).

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, and by NSERC in Canada. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the VERITAS instrument.

Facilities: NuSTAR - The NuSTAR (Nuclear Spectroscopic Telescope Array) mission, VERITAS. -

Software: HEAsoft (v6.22), Stingray (Huppenkothen et al. 2016), NuSTARDAS (v1.8.0), XSPEC (v12.9.1; Arnaud 1996), Naima (Zabalza 2015).