11institutetext: Department of Space, Earth and Environment, Chalmers
University of Technology, SE-412 96 Gothenburg, Sweden
11email<EMAIL_ADDRESS>22institutetext: Instituto Argentino
de Radioastronomía (CCT La Plata, CONICET; CICPBA; UNLP), C.C.5, (1894) Villa
Elisa, Buenos Aires, Argentina. 33institutetext: Space sciences, Technologies
and Astrophysics Research unit – STAR, University of Liège, Quartier Agora,
19c, Allée du 6 Août, B5c, B-4000 Sart Tilman, Belgium 44institutetext: School
of Physics and Astronomy, University of Southampton Highfield Campus,
Southampton SO17 1PS, UK. 55institutetext: Departament de Física Quàntica i
Astrofísica, Institut de Ciències del Cosmos (ICC), Universitat de Barcelona
(IEEC-UB), Martí i Franquès 1, E08028 Barcelona, Spain. 66institutetext: Joint
Institute for VLBI ERIC, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, the
Netherlands.

# Evidence for non-thermal X-ray emission from the double WR colliding-wind
binary Apep

S. del Palacio 11 F. García 22 M. De Becker 33 D. Altamirano 44 V. Bosch-Ramon
55
P. Benaglia 22 B. Marcote 66 G. E. Romero 22

(Received —; accepted —)

###### Abstract

Context. Massive colliding-wind binaries (CWBs) can be non-thermal sources.
The emission produced in their wind-collision region (WCR) encodes information
of both the shocks properties and the relativistic electrons accelerated in
them. The recently discovered system Apep, a unique massive system hosting two
Wolf-Rayet stars, is the most powerful synchrotron radio emitter among the
known CWBs, being an exciting candidate to investigate the non-thermal
processes associated with stellar wind shocks.

Aims. We intend to break the degeneracy between the relativistic particle
population and the magnetic field strength in the WCR of Apep by probing its
hard X-ray spectrum, where inverse-Compton (IC) emission is expected to
dominate.

Methods. We observe Apep with NuSTAR for 60 ks and combine this with a re-
analysis of a deep archival XMM-Newton observation to better constrain the
X-ray spectrum. We use a non-thermal emission model to derive physical
parameters from the results.

Results. We detect hard X-ray emission consistent with a power-law component
from Apep. This is compatible with IC emission produced in the WCR for a
magnetic field of $\approx 105$–190 mG, corresponding to a magnetic-to-thermal
pressure ratio in the shocks of $\approx 0.007$–0.021, and a fraction of $\sim
1.5\times 10^{-4}$ of the total wind kinetic power being transferred to
relativistic electrons.

Conclusions. This is the first time that the non-thermal emission from a CWB
is detected both in radio and high energies. This allows us to derive the most
robust constraints of the particle acceleration efficiency and magnetic field
intensity in a CWB so far, reducing the typical uncertainty of a few orders of
magnitude to just within a factor of a few. This constitutes an important step
forward in our characterisation of the physical properties of CWBs.

###### Key Words.:

stars: Wolf-Rayet, winds — radiation mechanisms: non-thermal — acceleration of
particles — X-rays: stars

##  1 Introduction

Colliding-wind binaries (CWBs) are binary systems in which the powerful winds
of the massive stars collide. The strong shocks at the wind-collision region
(WCR) produce very hot ($>10^{6}$ K) X-ray emitting plasma. Morever, they can
also accelerate relativistic particles (Eichler & Usov, 1993; Benaglia &
Romero, 2003) and constitute a subset of objects called Particle-Accelerating
Colliding-Wind Binaries (PACWBs; De Becker & Raucq, 2013). The efficiency of
this particle acceleration process is however still poorly constrained both
theoretically and observationally. The usually assumed scenario for particle
acceleration in PACWBs is Diffusive Shock Acceleration (DSA; Drury, 1983).

Relativistic electrons, expected in general to radiate more efficiently than
relativistic protons, can up-scatter stellar optical/ultraviolet photons to
X-ray or $\gamma$-ray emission by the Inverse Compton (IC) process.
Relativistic electrons can also radiate synchrotron emission in the radio band
by interacting with the magnetic fields in the WCR. Many CWBs present non-
thermal radio emission (De Becker & Raucq, 2013), but this is insufficient to
characterise both the relativistic electron population and the magnetic field
intensity in the emitter without severe partitioning assumptions (De Becker,
2018). Meanwhile, detections at hard X-rays and above remain scarce:
$\eta$-Car has been clearly detected in both hard X-rays (Hamaguchi et al.,
2018) and $\gamma$-rays (Tavani et al., 2009; Reitberger et al., 2015; H. E.
S. S. Coll. et al., 2020; Martí-Devesa & Reimer, 2021), while $\gamma^{2}$ Vel
has been recently confirmed as a $\gamma$-ray source (Martí-Devesa et al.,
2020), and a tentative detection of non-thermal hard X-rays ($E\lesssim 18$
keV) has been associated with HD 93129A (del Palacio et al., 2020).

The X-ray spectral energy distribution (SED) of a CWB is determined by the
thermal and the non-thermal radiation components, which depend on the WCR
properties, together with the local wind absorption (Pittard & Parkin, 2010).
The emission from individual stellar winds can only produce soft X-rays at
energies $\lesssim 1$ keV, and the total absorption from most stellar winds
and the interstellar medium (ISM) is not relevant above 2 keV. Thus, the SED
at energies $>3$ keV is determined solely by processes in the WCR. Thermal
processes are likely to dominate the SED up to $\sim 10$ keV given the high
wind velocities and consequent post-shock temperatures, and thus the non-
thermal processes can only be investigated at energies above 10 keV. It is
therefore necessary to have a broadband measurement of the X-ray SED to
disentangle these two components.

The system Apep is a peculiar case of a massive binary made up of two Wolf-
Rayet stars (Callingham et al., 2019). The stars are separated by more than
100 AU, which allows the stellar winds to accelerate to full speed before
collision. Radio observations of this system revealed that it is a very
powerful synchrotron source (Callingham et al., 2019), which also establishes
it as an efficient particle accelerator. In addition, Marcote et al. (2021)
confirmed that this emission raises from the WCR using very long baseline
interferometric observations. Further constraints on the radio spectrum by
Bloot et al. (2022) allowed del Palacio et al. (2022) to model the source
broadband emission in order to infer properties of the stellar winds and
predict the SED of the source at high energies. However, these predictions are
highly degenerate as it is not possible to disentangle the relativistic
particle energy distribution and the magnetic field strength in the WCR,
$B_{\mathrm{WCR}}$, solely from radio data (del Palacio et al., 2022). A
recent analysis of Fermi-LAT data in $\gamma$-rays placed stronger constraints
on the high-energy SED of Apep, but were still unable to detect its emission
(Martí-Devesa et al., 2022).

In addition, Apep was also observed in soft X-rays on different occasions with
XMM-Newton and Chandra. Callingham et al. (2019) analysed this data-set and
concluded that this source: i) is point-like in X-rays; ii) is not variable on
scales of years; and iii) has a predominantly thermal spectrum with a
significant absorption below 2 keV.

Here we aim to investigate the hard X-ray spectrum of Apep and search for
signatures of a non-thermal IC component, as predicted by del Palacio et al.
(2022). Measuring such a component can better constrain both the energy budget
in relativistic particles and the magnetic field strength in the WCR. With
this purpose we conducted observations of Apep with NuSTAR, probing for the
first time its spectrum at energies $>10$ keV. In this work we present the
analysis and interpretation of such observations.

##  2 Observations and data reduction

###  2.1 The system Apep

Table 1: Parameters of the WC+WN system Apep. Parameter | Value | Reference
---|---|---
Distance |  $d=2.4^{+0.2}_{-0.5}$ kpc | Callingham et al. (2019)
Projected system separation |  $D_{\mathrm{proj}}=47\pm 6$ mas | Han et al. (2020)
Projection angle | $\psi=85^{\circ}$ | del Palacio et al. (2022)
Wind momentum rate ratio | $\eta=0.44\pm 0.08$ | Marcote et al. (2021)
Stellar temperature |  $T_{\mathrm{eff,WN}}=65\,000$ K | Typical (e.g. Crowther, 2007; Hamann et al., 2019)
Stellar radius |  $R_{\mathrm{WN}}=6$ $\mathrm{R}_{\odot}$ | Typical (e.g. Hamann et al., 2019)
Wind terminal velocity |  $v_{\infty,\mathrm{WN}}=3500\pm 100$ km s-1 | Callingham et al. (2020)
Wind mass-loss rate |  $\dot{M}_{\mathrm{WN}}=(4\pm 1)\times 10^{-5}$ $\mathrm{M}_{\odot}$ yr-1 | del Palacio et al. (2022)
Wind mean atomic weight | $\mu_{\mathrm{WN}}=2.0$ | Typical (e.g. Leitherer et al., 1995)
Stellar temperature |  $T_{\mathrm{eff,\mathrm{WC}}}=60\,000$ K | Typical (e.g. Crowther, 2007; Sander et al., 2019)
Stellar radius |  $R_{\mathrm{WC}}=6.3$ $\mathrm{R}_{\odot}$ | Typical (e.g. Sander et al., 2019)
Wind terminal velocity |  $v_{\infty,\mathrm{WC}}=2100\pm 200$ km s-1 | Callingham et al. (2020)
Wind mass-loss rate |  $\dot{M}_{\mathrm{WC}}=(2.9\pm 0.7)\times 10^{-5}$ $\mathrm{M}_{\odot}$ yr-1 | del Palacio et al. (2022)
Wind mean atomic weight | $\mu_{\mathrm{WC}}=4.0$ | Typical (e.g. Cappa et al., 2004)

The massive binary Apep (2XMM J160050.7$-$514245) is located at
$RA=10^{\mathrm{h}}\,43^{\mathrm{m}}\,57.5^{\mathrm{s}}$, $DEC=-59\degr
32\arcmin 51.4\arcsec$ (J2000). Its orbit is wide, with a separation between
the stars of tens of AU (Han et al., 2020). The primary star is a WN star
while the secondary is a WC star. These stars have very massive and fast winds
with kinetic powers of $L_{\mathrm{WN}}\approx 1.5\times 10^{38}$ erg s-1 and
$L_{\mathrm{WC}}\approx 4.1\times 10^{37}$ erg s-1 (with an uncertainty of
$\approx 30$%; see Table 1). This constitutes an abundant energy reservoir to
feed emission processes at the WCR. In fact, the WCR in Apep is exceptionally
luminous, being the brightest PACWB detected at radio wavelengths (Callingham
et al., 2019). A more detailed list of the relevant system parameters is given
in Table 1.

This source has been observed in X-rays with both XMM-Newton and Chandra
during 2015–2021. Most of these observations have been previously analysed by
Callingham et al. (2019), who showed that the source does not present
significant variability in X-rays. We observed this system for the first time
with NuSTAR in 2022. In addition, we complemented this with a re-analysis of a
deep XMM-Newton observation in order to characterise better the X-ray SED of
Apep. We summarise the analysed observations in Table 2 and describe them in
more detail in the following subsections.

Table 2: Summary of the X-ray observations analysed. In the case of NuSTAR, A and B refer to FPMA and FPMB, respectively. Instrument | Obs. ID | Date (start) | Exposure time (ks) | Effective time (ks) | Offset (’)
---|---|---|---|---|---
NuSTAR | 30402001002 | 2022-06-17 | 31.0 | 31.0 (A), 30.8 (B) | 1.36 (A), 1.99 (B)
NuSTAR | 30402001004 | 2022-06-18 | 30.1 | 29.3 (A), 28.8 (B) | 1.43 (A), 2.06 (B)
XMM-Newton | 0742050101 | 2015-03-08 | 105 (PN), 137 (MOS2) | 79.9 (PN), 106.3 (MOS2) | 8.8

###  2.2 XMM-Newton

Apep is in the field of view of several archival XMM-Newton observations. Of
these, Obs. 0742050101 is the deepest one ($>100$ ks) and therefore the one we
chose to analyse. The major drawback of this observation is the large offset
from on-axis to the position of Apep, which is $8.8\arcmin$. The source
appears in both PN and MOS2 cameras, and the observation was carried out in
full frame mode. Other details of this observation are summarised in Table 2.

Data processing was performed using the Science Analysis Software SAS v.20.0.0
and the calibration files (CCF) available in August 2022. We used the
metatasks emproc and epproc to reduce the data. We then filtered periods of
high background or soft proton flares. Standard screening criteria were
adopted, namely pattern $\leq 12$ for MOS and pattern $\leq 4$ for PN. We
determined good time intervals by selecting events with PI$>$10000 and
PATTERN$==$0 and adopting the standard rejection thresholds RATE$\leq$0.35 for
MOS2 and RATE$\leq$0.4 for PN. The effective time after filtering is reported
in Table 2. These values are 25–30 ks shorter than those of Callingham et al.
(2019), which suggests that our selection of GTIs was more conservative.

To produce the spectra, the radius of the extraction region was set to 50″, as
the source is $8.8\arcmin$ off axis and therefore the PSF is larger than for
on-axis sources (for which typically 10–30″ is used). The background spectrum
was extracted in an elliptical region located on the same chip, in an area
devoid of point sources, selected using the ebkreg task. However, the adopted
background region has a negligible impact on the results given that the source
is very bright. Adequate response matrix files (RMF) and ancillary response
files (ARF) were produced using the dedicated tasks (rmfgen and arfgen,
respectively) for all spectra. On this last point, we note that the standard
psfmodel ELLBETA does not work well for MOS2 due to Apep being a very bright
source and quite off-axis ($8.8\arcmin$), so we used the psfmodel=EXTENDED
option (XMM support, priv. comm.). This correction improves the match between
the MOS2 and the PN spectra (a mismatch around 6 keV can be seen in
Supplementary Information Fig. 4 from Callingham et al. 2019, where the
standard ELLBETA model was used for MOS2). Other parameters adopted for the
ARF file were extendedsource=no, detmaptype=psf, and applyabsfluxcorr=yes,
being the last one used to improve the cross-calibration between XMM-Newton
and NuSTAR. We finally grouped the spectra using the task ftgrouppha with
grouptype=opt.

###  2.3 NuSTAR

The NuSTAR X-ray observatory was launched in 2012 and its major asset is its
unique imaging capacity in hard X-rays. The observatory includes two co-
aligned X-ray grazing incidence telescopes, known as FPMA and FPMB for their
focal plane modules, which are comprised of four rectangular solid state
CdZnTe detectors. NuSTAR is capable of observing in the 3–79 keV energy range
with an angular resolution of 18″(half power diameter of 58″; Harrison et al.,
2013).

We observed the massive binary Apep in June 2022 with NuSTAR under program
8020 (PI: del Palacio). The observations were carried out in two 30-ks visits,
adding up to roughly 60 ks of exposure time with both cameras. We refer to the
observations in each epoch as 2022a and 2022b. We summarise the relevant
details of these observations in Table 2.

We reduced the data using Heasoft 6.30.1 and the latest calibration files
available in June 2022 (CALDB 4.9.7-0). We used the nupipeline task to create
level 2 data products with the options saacalc=2, saamode=optimized and
tentacle=yes to filter high background epochs. This led to negligible data
loss in the 2022a observations and $<3\%$ data loss in the 2022b
observation111http://www.srl.caltech.edu/NuSTAR_Public/NuSTAROperationSite/SAA_Filtering/SAA_Filter.php..
We then used the nuproducts task to create level 3 data products. We extracted
the source spectrum from a 55″ region centred in Apep, while the background
was extracted from an ellipse located in the same chip, sufficiently far from
the source as to avoid contamination. The selected background region also
avoids contamination from the supernova remnant G330.2+1.0, as shown in Fig.
1. Further analysis of the influence of the selected background region is
presented in Appendix B. Finally, we binned the spectra using the task
ftgrouppha with the option grouptype=opt.

###  2.4 XSPEC spectral model

Once the spectra have been obtained, one needs to fit them with a spectral
model in order to extract physical information. Any model we adopt should be
both physically-motivated and simple enough as to reproduce the data without
requiring a very large number of parameters. A spectral model for a PACWB
should include an absorption component –that can take into account both
internal absorption in the stellar winds and external absorption in the ISM–,
a thermal component –dominated by the WCR emission–, and a non-thermal
component –relevant only at $E>10$ keV, also produced by the WCR–. The thermal
emission from CWBs is usually approximated using an apec model (Smith et al.,
2001), because it is simple and can reproduce well the emission from an
optically thin plasma. Multiple apec components can be used to emulate the
temperature gradient along the WCR (e.g. Pittard & Parkin, 2010). However, the
apec model assumes that electrons and ions are in equilibrium, and the shocks
in the WCR can be collisionless under certain conditions, leading to an
ionisation that can be out-of-equilibrium. This depends on the relation
between two timescales: the dynamical timescale, $t_{\mathrm{dyn}}$, and the
electron-ion temperature equalisation timescale, $t_{\mathrm{eq}}$. The
timescale $t_{\mathrm{dyn}}$ depends on the characteristic size of the WCR
($\sim D$) and on the velocity at which the material is advected away ($\sim
v_{\infty}$), while $t_{\mathrm{eq}}$ depends on the post-shock temperature
and density, and therefore on $v_{\infty}$ and $\dot{M}$. If
$t_{\mathrm{dyn}}<t_{\mathrm{eq}}$, the electrons and ions cannot get in
thermal equilibrium through Coulomb interactions before the post-shock plasma
is advected away. This condition can be summarised through the parameter
$\zeta_{\mathrm{eq}}=\dfrac{t_{\mathrm{dyn}}}{t_{\mathrm{eq}}}\approx\dfrac{13.36}{\bar{\mu}\,\mu^{1.5}}\left(\dfrac{\dot{M}}{10^{-6}\,M_{\odot}}\right)\left(\dfrac{V_{\infty}}{1000\,\mathrm{km\,s^{-1}}}\right)^{5}\left(\dfrac{10^{14}\,\mathrm{cm}}{D}\right)$,
such that if $\zeta_{\mathrm{eq}}<1$ the difference between electron and ion
temperatures should be taken into account (Zhekov & Skinner, 2000). For long
period binaries this condition is more likely to be fulfilled, such as in the
case of Apep. Actually, for the conditions in the shocks of the primary and
secondary and using the parameters given in Table 1, we obtain $\zeta\sim
0.01-0.1$. Thus, the use of a non-equilibrium model, such as pshock, is
justified in this case and we explored this possibility as well.

Regarding the non-thermal emission, the simplest way to parameterise it is as
a power-law component. We can constrain its spectral index considering that
this component is expected to be IC radiation emitted by the same relativistic
electron population in the WCR that is responsible for the synchrotron
emission observed in the radio band (del Palacio et al., 2022). For a flux
density $S_{\nu}\propto\nu^{\alpha}$ in the radio band, the specific photon
flux density $F\propto E^{-\Gamma}$ in X-rays has a spectral index
$\Gamma=-\alpha+1$ (del Palacio et al., 2022). In the case of Apep, a value of
$\alpha=-0.72$ was reported by Callingham et al. (2019), which leads to
$\Gamma=1.72$. We note that this value is slightly steeper than the canonical
$\alpha=-0.5$ ($\Gamma=1.5$) expected to arise from electrons accelerated by
DSA in high Mach number shocks under the test-particle assumption.
Nonetheless, this trend of steeper spectra is also seen in supernova remnants
and can be related to the back-reaction of cosmic rays in the shocks (e.g.
Drury, 1983; Gabici et al., 2019).

Finally, the emitted X-ray radiation can be absorbed intrinsically in the
source and externally in the ISM. In XSPEC, the standard model used to
calculate the ISM absorption is TBabs. The value of the $N_{\mathrm{H}}$
column can be taken from HI4PI Collaboration et al.
(2016)222https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl. For Apep,
the value retrieved is $N_{\mathrm{H}}=1.63\times 10^{22}$ cm-2. Additional
intrinsic absorption (mostly by the stellar winds) can be included using a
phabs model. Throughout this work the confidence intervals are obtained using
the error command in XSPEC and given at a 1-$\sigma$ level unless stated
otherwise.

##  3 Results

We now focus on the analysis and model fitting of the spectra obtained with
the NuSTAR and XMM-Newton observatories.

###  3.1 XMM-Newton

Figure 1: RGB XMM-Newton image of the field of view of Apep (red = 0.3–1.2
keV, green = 1.2–2.5 keV, blue = 2.5–8 keV) from the EPIC-PN detector. We also
mark the position of the supernova remnant G330.2+1.0.

In Fig. 1 we show an RGB exposure-corrected image of the field of view of XMM-
Newton. In the figure we also show the source and background extraction
regions.

To fit the XMM-Newton spectra, we first considered a model of the form
constant*TBabs*apec, with abundances set to Wilms (Wilms et al., 2000) as done
by Callingham et al. (2019). The normalisation constant was set to unity for
PN and then fitted to MOS2, obtaining $C=1.06\pm 0.01$. We show the fitted
spectra in Fig. 2. In general, we found very similar results to Callingham et
al. (2019) ($kT\approx 5.1$ keV, $N_{\rm H}\approx 2.7\times 10^{22}$ cm-2,
abundance $A\approx 0.5$, observed flux $F_{\mathrm{0.3-10\,keV}}\approx
8\times 10^{-12}$ erg s-1 cm-2), despite the differences in the data reduction
(Sect. 2.2). Nonetheless, the goodness of the fit was actually quite poor,
yielding large structured residuals around 1–2.5 keV at energies coincident
with known Si and S transitions (Fig. 2) and a C-Stat 439.0/224. This
motivated us to look for a better model.

As discussed in Sect. 2.4, adopting a pshock model instead of an apec is
physically justified in the case of Apep. The use of a pshock model required
only one additional free parameter, namely the upper limit on ionisation
timescale ($\tau_{\mathrm{u}}$), and it improved significantly the quality of
the fit (C-Stat 311.5/223). We note that this affected only the low energy
portion of the spectrum, in particular by improving the ratios in the Si and S
lines. This suggests that ionisation was indeed out of equilibrium, and is
also consistent with the value obtained of $\tau_{\mathrm{u}}<10^{12}$ s cm3.
However, this had a completely negligible ($\sim 1\%$) impact on the spectra
at energies above 3 keV. Further improvement can be done by setting variable
abundances: although Callingham et al. (2019) claimed that using a vapec model
did not improve the fit significantly, we found that a vpshock model can
indeed introduce a significant improvement, reaching C-Stat 270.9/220. This
can also be appreciated in Fig. 2, which highlights the smaller residuals
retrieved in the 1–2.5 keV energy range with this model. When fitting the
individual abundances, we obtained for Fe a similar value to that of $A$,
indicating that this element dominates the value of $A$ for fixed relative
abundances. A few elements presented different abundances (mainly Ne), while
others like C and N were not constrained by the data and were left fixed to
one. At last, as the stellar winds are expected to contribute significantly to
the (photoelectric) absorption, we changed the absorption model to
TBabs*vphabs, fixing $N_{\mathrm{H}}=1.63\times 10^{22}$ cm-2 for the TBabs
component (Sect. 2.4) and fixing the abundances of the vphabs to those of the
vpshock model. This improved the fit slightly (C-Stat 268.4/220) without
adding extra free parameters, and is therefore our preferred model for the
XMM-Newton data. In Table 3 we present in detail all the fitted parameters for
the most relevant models, including the data from NuSTAR as discussed in the
next section. Taking as a reference the PN camera, the observed flux in the
0.3–10 keV energy range is $F_{\mathrm{0.3-10\,keV,obs}}=(7.92\pm 0.10)\times
10^{-12}$ erg s-1 cm-2, the ISM-unabsorbed flux is
$F_{\mathrm{0.3-10\,keV,unabs}}=(9.69\pm 0.05)\times 10^{-12}$ erg s-1 cm-2,
and the unabsorbed flux from the vpshock component only is
$F_{\mathrm{0.3-10\,keV,vpshock}}=(1.89\pm 0.04)\times 10^{-11}$ erg s-1 cm-2.

Figure 2: Apep XMM-Newton spectra in the 0.3–10 keV energy range. In the left
panel the fitted model is TBabs*apec, while in the right panel it is
TBabs*vphabs*vpshock. The latter model is preferred as it leads to smaller
residuals between 1–2.5 keV.

###  3.2 NuSTAR

In Fig. 3 we present an image in the 3–20 keV energy range ($PI$ channels
35–460) with NuSTAR for each observing epoch and camera. Apep is clearly
detected and no other bright sources appear in the field. We note the presence
of straylight in the FPMB observations, although this is not problematic given
that Apep is far away from it. In Fig. 3 we also show the selected source and
background extraction regions.

Figure 3: NuSTAR image in the 3–20 keV energy range for the 2022a (left
panels) and 2022b (right panels) observations. For each epoch and camera we
show the selected source and background extraction regions. Straylight can be
seen at the bottom right corner of the FPMB observations.

The obtained spectra for each epoch and camera are shown together in Fig. 7
for comparison. All observations reveal a very similar spectrum in which the
source is detected above the background up to $\gtrsim 20$ keV. We calculated
the integrated source flux in the 3–10 keV and 10–25 keV energy ranges for
each observation to make a quantitative comparison333For this test we adopted
the models described in Sect. 2.4 and used later in Sect. 3, and for all of
them we reached the same conclusion.. In all cases the fluxes differ in less
than 10% and are compatible within 1-$\sigma$ level. We checked whether it was
possible to combine the data from both observations to increase the signal-to-
noise. For this we co-added the spectra of both cameras for each epoch using
the addspec task with the options qaddrmf=yes qsubback=yes and the default
value bexpscale=1000. We repeated the calculation of the integrated fluxes and
obtained variations below 2% ($<1\sigma$ difference), indicating that they are
perfectly compatible, which is to be expected considering that the
observations were taken very close in time to each other.

We then co-added the spectra of both epochs for each camera to compare if
there was any systematic difference between FPMA and FPMB. We obtained that
the FPMA fluxes in the 3–10 keV were slightly ($\approx 5\%$) higher than for
FPMB (with a 1-$\sigma$ significance), while in the 10–25 keV range the fluxes
between both cameras match up to 2% and this difference is less significant
($<1\sigma$). This is within the calibration uncertainties of the instrument
(Madsen et al., 2015, 2022). From these tests we concluded that both
observations are compatible but with a small difference between the cameras.
We therefore decided to work with both NuSTAR observations co-added for each
camera separately to improve the signal-to-noise at the highest energies while
preventing cross-calibration errors to introduce further errors to the fitting
of the spectra. In Fig. 7 we show the spectrum obtained for the combined
observations. In this case the statistics are improved, as expected, and the
source is brighter than the background up to $\sim$25 keV. In the fitting we
include data up to 35 keV as there is valuable spectral information between
25–35 keV, and we make use of Cash statistics in XSPEC to deal properly with
the low number of counts in this energy range.

We tried fitting different models to the NuSTAR spectra. The main conclusion
is that these spectra are not sensitive to the adopted absorption model (as
the absorption at energies $>3$ keV is negligible) nor to the specifics of the
emission model (as the information from emission lines is poor). Thus, we
simply adopted the same model used to fit the XMM-Newton spectra,
constant*TBabs*vphabs*vpshock, leaving the abundances and absorption fixed
(but allowing temperature and normalisation to vary). Re-fitting the spectrum
allowed us to obtain the normalisation constant between FPMA (taken as unity)
and FPMB, $C=0.958\pm 0.015$. A high-temperature ($k\,T\sim 5$ keV) thermal
component naturally extends to energies above 10 keV and can explain most of
the emission detected with NuSTAR. In Fig. 4 we show the spectra and the
residuals. However, the residuals at energies $>20$ keV become increasingly
large, with deviations between 3–7 $\sigma$. Such deviations can be attributed
to the putative non-thermal component. We therefore included an additional
power-law (po) component. This way we obtained a significantly improved fit,
as can be seen in the residuals in the right panel of Fig. 4, as well as in
the lower C-stat value (which decreased from 237.5/182 to 203.1/181). For
completeness, we note that an additional high-temperature vpshock component
with $kT>10$ keV also leads to a similar improvement in the fit (although with
higher residuals at energies $\sim 30$ keV). However, such high temperatures
are not expected in CWB shocks, while a non-thermal component should arise
naturally given the (already established in the radio band) presence of non-
thermal particles. Thus, the main conclusion here is that the NuSTAR data by
itself strongly supports the existence of an additional high-energy component,
which we interpret hereon as a (non-thermal) power-law component. More robust
results can be derived by including the data from XMM-Newton consistently,
which we address in the next section.

Figure 4: Apep NuSTAR spectra in the 3–35 keV energy range from co-adding the
observations for each camera separately. In the left panel the fitted model is
TBabs*vphabs*vpshock, while in the right panel an additional power-law
component is added as TBabs*vphabs*(vpshock+po).

###  3.3 Joint analysis

The previous fitting to each instrument separately allowed us to understand
the behaviour of the X-ray spectra and which observations are more sensitive
to each component. Namely, the XMM-Newton data allows us to better constrain
the absorption and thermal emission models, whereas the NuSTAR data is
sensitive to the putative high-energy non-thermal component. We now present a
joint analysis of the whole data set for the same models as before. In
general, we allowed for a different normalisation constant between the
different instruments444Differences in absolute flux calibration between XMM-
Newton and NuSTAR can be of 5–15% (Madsen et al., 2017). and tied the
remaining physical parameters.

In Fig. 5 we show the spectra and the combined fitting for two different
models, while in Table 3 we detail the fitted parameters for the most relevant
models. The TBabs*apec model, used previously by Callingham et al. (2019),
fails to reproduce the spectra measured by both XMM-Newton and NuSTAR. The
TBabs*vphabs*vpshock model satisfactorily fits the XMM-Newton spectra, but it
struggles with the NuSTAR spectra at energies above 20 keV (Fig. 5). Finally,
the TBabs*vphabs*(vpshock+po) model can fit all the spectra simultaneously. In
this case the C-stat diminished from 510.5/402 to 491.2/401, being the
improvement more significant for NuSTAR at the expense of a slightly worse fit
for XMM-Newton (Table 3). We further quantified the significance of the power-
law component using the task simftest in XSPEC. We ran 11 000 simulations and
obtained a probability $<0.01\%$ of the data being consistent with a model
without the power-law component, which corresponds to a significance
$>3.91\sigma$. We also note that the inclusion of the power-law component has
little effect on the overall fit, mainly by diminishing slightly the
temperature of the thermal component (from $\approx$5.3 to $\approx$4.9 keV;
Table 3). At last, we introduced a cflux component to calculate both the total
flux in the 10–30 keV band, $F_{\mathrm{10-30\,keV}}=(1.99\pm 0.11)\times
10^{-12}$ erg s-1 cm-2, and the one coming only from the power-law component,
$F_{\mathrm{10-30\,keV}}=4.8^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2.

We note that the flux of the power-law component is susceptible to the
background extraction region chosen for NuSTAR, although its presence is
always statistically favoured. A detailed exploration of different background
regions and the caveats in their selection is given in Appendix B, together
with a complementary analysis of the background using nuskybgd Wik et al.
(2014). The latter yields a flux of $F_{\mathrm{10-30\,keV}}\sim
3.9^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2 for the power-law component,
which is consistent with the previous result within $1\sigma$.

We also corroborated whether the value adopted for the power-law index had a
significant impact on the results. For $\Gamma$ in the range 1.6–1.8,
$F_{\mathrm{10-30\,keV}}$ varied only slightly ($\sim 2\%$), much less than
1$\sigma$. Thus, the value adopted for $\Gamma$ does not affect the integrated
flux in the hard X-ray band. We conclude that the presence of a power-law
component is robust, and that its flux can only be measured with a rather high
uncertainty.

Figure 5: Apep unfolded X-ray spectra in the 1–35 keV energy range with XMM-Newton and NuSTAR. The fitted models are constant*TBabs*vphabs*vpshock (left panel) and constant*TBabs*vphabs*(vpshock+po) (right panel). The spectra has been rebinned for clarity. The fitted parameters are given in Table 3. Table 3: Results of the fitting of NuSTAR and XMM-Newton spectra of Apep using different models. C-statistics were used. The errors at 1-$\sigma$ level are specified for all parameters. A multiplicative constant was also added with a value set to 1 for FPMA and fitted for the other instruments to $C_{\mathrm{FPMB}}=0.957\pm 0.015$, $C_{\mathrm{PN}}=0.849\pm 0.010$, and $C_{\mathrm{MOS2}}=0.902\pm 0.011$. The best model is highlighted in boldface. Parameter | Units | TBabs*apec | TBabs*vphabs*vpshock | TBabs*vphabs*(vpshock+po)
---|---|---|---|---
$N_{\rm H}$ |  $10^{22}$ cm-2 | $2.67\pm 0.02$ | 1.63 (fixed), $1.23^{+0.06}_{-0.06}$ | 1.63 (fixed), $1.20^{+0.07}_{-0.05}$
$k\,T$ | keV | $5.15\pm 0.07$ | $5.27^{+0.07}_{-0.06}$ | $4.91\pm 0.11$
norm1 |  $10^{-2}$ cm-5 | $1.102^{+0.020}_{-0.019}$ | $1.023^{+0.018}_{-0.017}$ | $1.014^{+0.017}_{-0.019}$
$A$ |  | $0.52\pm 0.02$ | — | —
Fe |  | — | $0.54\pm 0.02$ | $0.58\pm 0.02$
Ne |  | — | $2.88^{+0.38}_{-0.34}$ | $3.36^{+0.50}_{-0.44}$
S |  | — | $0.47\pm 0.07$ | $0.43^{+0.08}_{-0.06}$
Ca |  | — | $0.58^{+0.23}_{-0.21}$ | $0.57^{+0.20}_{-0.22}$
$\tau_{\mathrm{u}}$ |  $10^{11}$ s cm3 | — | $9.01^{+0.96}_{-0.86}$ | $9.36^{+1.31}_{-0.90}$
$\Gamma$ |  | — | — | 1.72 (fixed)
norm2 | 10-4 keV-1 cm-2 s-1 | — | — | $1.30^{+0.26}_{-0.33}$
$C_{\mathrm{stat}}$/bins (PN) |  | 234.3/112 | 136.6/112 | 137.4/112
$C_{\mathrm{stat}}$/bins (MOS2) |  | 206.9/115 | 134.6/115 | 141.0/115
$C_{\mathrm{stat}}$/bins (FPMA) |  | 125.7/91 | 124.8/91 | 112.1/91
$C_{\mathrm{stat}}$/bins (FPMB) |  | 122.6/95 | 114.5/95 | 101.4/95
$\mathbf{C_{\mathrm{stat}}}$/d.o.f. (total) |  | 689.5/406 | 510.5/402 | 491.2/401

##  4 Discussion

The main result from our spectral analysis is the detection of hard X-ray
emission consistent with a power-law component with a flux of
$F_{\mathrm{10-30\,keV}}=4.8^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2. We
also constrain the thermal emission from the WCR much better. Previous
estimates of the plasma temperature by Callingham et al. (2019) were highly
uncertain, spanning a range $kT\sim 4.7$–$6.3$ keV (for the particular XMM-
Newton observation that we also analysed, they obtained $kT\sim 4.9$–5.3 keV).
By means of a more careful analysis of the XMM-Newton observations, combined
with the unique information provided by NuSTAR above 10 keV, we constrain it
to $kT\approx 4.85$–$5.07$ keV.

We now focus on the interpretation of the non-thermal component. To be more
conservative, in what follows we consider additional sources of errors in the
flux values, both statistical and systematic. For example, a 10% systematic
error due to absolute calibration uncertainties are estimated for NuSTAR
(Harrison et al., 2013). There is also an additional dependency of the
retrieved flux with the chosen background region, as shown in Appendix B.
Based on this, we adopt a less constrained flux of
$F_{\mathrm{10-30\,keV}}=(2$–$6)\times 10^{-13}$ erg s-1 cm-2 in our analysis,
which is slightly broader than the 90% confidence interval obtained from the
analysis with nuskybgd in Appendix B
($F_{\mathrm{10-30\,keV}}\sim(2.3$–$5.6)\times 10^{-13}$ erg s-1 cm-2).

###  4.1 Modelling the hard X-ray emission

According to del Palacio et al. (2022), the hard X-ray emission should arise
from IC scattering of stellar photons by relativistic electrons in the WCR.
These same electrons should also be responsible for the non-thermal
(synchrotron) emission observed in the radio band (Callingham et al., 2019;
Marcote et al., 2021; Bloot et al., 2022). In order to relate the observed
fluxes with the particle acceleration in the shocks, we need to take into
account that only a fraction of the wind kinetic power can be converted into
relativistic particles in the shocks, that this energy is distributed into
electrons and protons, and that each of this particle species radiate only a
fraction of their energy at any given frequency range (e.g. De Becker & Raucq,
2013).

Thus, to model this emission, we use a code based on the non-thermal emission
model presented in del Palacio et al. (2016) with the system parameters listed
in Table 1. The model solves for the acceleration and transport of
relativistic particles for both shocks in the WCR (one for each stellar wind).
These relativistic particles radiate by different processes (synchrotron, IC,
p-p collisions), and this radiation is mitigated by absorption processes in
the stellar winds or radiation fields. The model has two free parameters that
determine the leptonic emission: the ratio between the magnetic field pressure
to thermal pressure in the WCR, $\eta_{B}$, and the fraction of the available
power at the shocks that is converted into relativistic electrons,
$f_{\mathrm{NT,e}}$. The available power for particle acceleration is the wind
kinetic power injected perpendicularly into the WCR shocks. Denoting this
power by $L_{\mathrm{inj,\perp}}$ and the total wind kinetic power of a star
by $L_{\mathrm{w}}=0.5\dot{M}v_{\mathrm{w}}^{2}$, we can write
$L_{\mathrm{inj,\perp}}=\epsilon L_{\mathrm{w}}$, with
$\epsilon_{\mathrm{WN}}=8\%$ and $\epsilon_{\mathrm{WC}}=18\%$ for the system
Apep. Here, the value of $\epsilon$ depends on the geometry of the WCR, which
is governed by the value of the wind-momentum rate ratio $\eta$, and is
calculated numerically in the model. Further details of the model are
described in Appendix C.

Figure 6: Modelled broadband non-thermal SED of Apep. We show SEDs for the two
extreme combinations of parameters that are compatible with the hard X-ray
emission component detected with NuSTAR (this work) and with the ATCA radio
data (Callingham et al., 2019). Dashed lines show the synchrotron component,
dot-dashed lines the IC component, dotted lines the p-p component, and solid
lines the total SED. We also show the Fermi-LAT upper-limits from Martí-Devesa
et al. (2022) and the sensitivity curve for CTA (Funk et al., 2013).

It is possible to tie the two free parameters, $f_{\mathrm{NT,e}}$ and
$\eta_{B}$, by modelling the observed synchrotron component (del Palacio et
al., 2016). In this case the relation is $f_{\mathrm{NT,e}}\eta_{B}=constant$
(del Palacio et al., 2020). However, it is not possible to break the
degeneracy between these two parameters from radio data alone. Fortunately,
this can be solved when observations in hard X-rays measure the flux from the
IC component, which depends only on $f_{\mathrm{NT,e}}$ as
$F_{\mathrm{IC}}\propto f_{\mathrm{NT,e}}$ (del Palacio et al., 2020). Once
$f_{\mathrm{NT,e}}$ is derived from the measured hard X-ray flux, we can
obtain $\eta_{B}$ by fitting the synchrotron emission to the observed radio
flux of $\approx 120$ mJy at 2 GHz (Callingham et al., 2019). In Fig. 6 we
show the SEDs fitted using this procedure. In addition, for a given value of
$\eta_{B}$ we can calculate the magnetic field in the apex of the WCR,
$B_{\mathrm{WCR}}$.

Previous estimates by del Palacio et al. (2022) based only on the synchrotron
emission from the source had more than one order of magnitude of uncertainty
in $f_{\mathrm{NT,e}}$, namely $f_{\mathrm{NT,e}}\approx(0.11$–$2.7)\times
10^{-3}$. Our new estimates based on the NuSTAR detection yield a very well-
constrained value with less than a factor two uncertainty,
$f_{\mathrm{NT,e}}\approx(0.7$–$2)\times 10^{-3}$. This corresponds to roughly
$1.5\times 10^{-4}$ of the total wind kinetic power being converted into
relativistic electron acceleration. Moreover, the magnetic field in the WCR
was also poorly constrained by del Palacio et al. (2022) to
$B_{\mathrm{WCR}}\approx 70$–400 mG ($\eta_{B}=(3$–$100)\times 10^{-3}$),
while now we constrained it to $B_{\mathrm{WCR}}\approx 105$–190 mG
($\eta_{B}=0.007$–$0.021$). This translates into a ratio between the energy
density in relativistic electrons and the magnetic field of
$U_{\mathrm{e}}/U_{B}\approx 0.02$–0.2. For reasonable values of a ratio
between power injected in electrons and protons of $K_{\mathrm{e,p}}<0.1$,
this leads to a magnetic field in subequipartition with the non-thermal
particles, supporting the possibility of relativistic protons driving the
magnetic field amplification (Bell, 2004).

These values can be compared with those found for the O+O binary HD 93129A
during its periastron passage (del Palacio et al., 2020),
$f_{\mathrm{NT,e}}\approx 6\times 10^{-3}$ and $\eta_{B}\sim 0.02$
($B_{\mathrm{WCR}}\approx 0.5$ G). Comparisons with other systems are
complicated given the uniqueness of the detection of a PACWB in both radio and
high energies. Nonetheless, we can comment on the PACWB $\eta$-Car, for which
non-thermal hard X-rays were also detected with a power-law index $\Gamma\sim
1.65$ (although poorly constrained; Hamaguchi et al., 2018). Another systems
studied by De Becker (2018), based on radio observations and equipartition
assumptions between relativistic particles and magnetic fields, yielded that
the fraction of the wind kinetic power converted into relativistic electrons
is $\sim 10^{-4}$–$10^{-6}$ for Cyg OB2 #8a, $\sim 10^{-7}$–$10^{-9}$ for WR
140, and $\sim 10^{-5}$–$10^{-7}$ for HD 167971. Compared with the value
obtained here for Apep ($1.5\times 10^{-4}$), it is clear that this system is
a much more efficient electron accelerator. This is consistent with the fact
that this binary is the brightest synchrotron-emitting PACWB. We also tried to
compare the values of $\eta_{B}$ with those derived from De Becker (2018).
These values are $10^{-7}$ for WR 140, $2\times 10^{-4}$ for Cyg OB2 #8a, and
$5\times 10^{-5}$ for HD 167971, but these have an uncertainty that spans 2–3
orders of magnitude, so all we can say is that our value of $\eta_{B}\approx
10^{-2}$ is exceptionally well-constrained. In addition, Pittard et al. (2021)
fitted the radio SED of the system WR 146 to derive a magnetic field
compatible with $\eta_{B}\approx 10^{-3}$, although these authors required a
very large particle efficiency in return ($f_{\mathrm{NT}}\approx 0.3$).

One last parameter we could derive from our results is the surface magnetic
field of the stars. For this we assumed a toroidal stellar magnetic field that
drops as $r^{-1}$ and is adiabatically compressed in the WCR shocks, and a
stellar rotation velocity of $V_{\mathrm{rot}}\sim 0.1v_{\infty}$ (del Palacio
et al., 2016, and references therein). Under these assumptions we obtain
values of the surface stellar magnetic fields in the ranges
$B_{\mathrm{WN}}=650$–1100 G and $B_{\mathrm{WC}}=280$–490 G. Nonetheless, it
is possible that magnetic field amplification processes take place in the WCR
shocks (e.g. Bell, 2004; Pittard et al., 2021). In this case the
aforementioned values should actually be interpreted as upper limits to the
stellar magnetic fields.

###  4.2 Predictions of $\gamma$-ray emission

The previous estimates of the power in relativistic electrons also allowed us
to compute the expected IC luminosity in the $\gamma$-ray domain. In Fig. 6 we
show the modelled broadband SED extending to $\gamma$-ray energies, together
with characteristic sensitivity thresholds of $\gamma$-ray observatories. We
first focus on the 0.1–100 GeV energy range, which can be tested with
observations with the Fermi-LAT instrument. The $\gamma$-ray luminosity in
this case is $F_{\mathrm{0.1-100\,GeV}}=(1.9\pm 0.9)\times 10^{-12}$ erg s-1
cm-2, though it can be larger if a hadronic component is included (for
example, $K_{\mathrm{e,p}}=0.04$ yields to a total flux of
$F_{\mathrm{0.1-100\,GeV}}=(2.0\pm 1.5)\times 10^{-12}$ erg s-1 cm-2). These
values are mostly consistent with a non-detection of this source with Fermi-
LAT at a level of $F_{\mathrm{0.1-100\,GeV}}\sim(1$–$2)\times 10^{-12}$ erg
s-1 cm-2 (Martí-Devesa et al., 2022). The small tension between the higher
fluxes predicted for the lower magnetic field scenarios might suggest that the
higher magnetic field scenarios are to be preferred (Fig. 4). Nonetheless,
this tension can also be attributed to even small uncertainties in the
particle energy distribution that lead to a significant difference in the
predicted $\gamma$-ray fluxes. Assuming that the injected electron energy
distribution is slightly harder, $p=2.3$ (equivalently, $\Gamma=1.65$), we
obtain $F_{\mathrm{0.1-100\,GeV}}\sim 4.2\times 10^{-12}$ erg s-1 cm-2, while
a slightly steeper distribution, $p=2.55$ ($\Gamma=1.77$), yields an IC
emission of $F_{\mathrm{0.1-100\,GeV}}\sim(1.0\pm 0.7)\times 10^{-12}$ erg s-1
cm-2. Thus, a hardening of the electron energy distribution is strongly
disfavoured as it would overpredict the $\gamma$-ray luminosity. Moreover, we
conclude that the source is either close to be detected by Fermi, or that a
non-detection with deeper sensitivity ($F_{\mathrm{0.1-100\,GeV}}<10^{-12}$
erg s-1 cm-2) would mean that the SED softens at energies above hard X-rays,
which in turn would require a softening in the electron energy distribution at
energies $E_{\mathrm{e}}>100$ MeV.

Finally, we address the prospects for detection of TeV emission from Apep. We
predict an IC flux of $F_{\mathrm{TeV}}\sim 1.8\times 10^{-14}$ erg s-1 cm-2
in the 0.1–10 TeV energy range, although the poorly-constrained hadronic
component is likely dominant (del Palacio et al., 2022): assuming
$K_{\mathrm{e,p}}=0.04$, the predicted total flux (p-p + IC) is
$F_{\mathrm{TeV}}\sim 8\times 10^{-14}$ erg s-1 cm-2. Moreover, small
variations in the spectral index of the particle energy distribution
($p=2.3$–2.5) can lead to significantly different TeV fluxes,
$F_{\mathrm{TeV}}\sim(0.3$–$18)\times 10^{-14}$ erg s-1 cm-2. Only the higher
fluxes would potentially be detectable by the Cherenkov Telescope Array (CTA;
Funk et al., 2013), but these are already disfavoured in view of the lack of
detections of GeV emission. Thus, the TeV emission from Apep seems too faint
to be detected with current and upcoming TeV observatories.

##  5 Conclusions

We presented the first hard X-ray view of the PACWB Apep. This system is the
brightest synchrotron source among the known PACWBs. The NuSTAR spectrum
revealed strong evidence of a power-law component, consistent with the
predicted IC emission produced by relativistic electrons in the WCR. The
detection of this non-thermal high-energy emission from a system that also
presents non-thermal emission in the radio band represents an observational
breakthrough in the study of PACWBs. In particular, it has allowed us to place
the tightest constraints to the magnetic field and electron acceleration
efficiency in the WCR of a PACWB. We also predict that Apep is close to be
detected at $\gamma$-rays with Fermi unless the electron energy distribution
softens at energies $>100$ MeV.

We highlight the importance of multi-wavelength observations for improving our
understanding of PACWBs. Unfortunately, the high-energy emission from these
systems is rather weak and difficult to detect, but at least for the brightest
sources observations in the hard X-ray and high energy $\gamma$-rays bands
have proven to be successful, paving the way for moving forward in the
research of PACWBs.

###### Acknowledgements.

We would like to thank the referee for providing useful feedback that helped
us to improve the manuscript. This work was carried out in the framework of
the PANTERA-Stars555https://www.astro.uliege.be/~debecker/pantera/ initiative.
FG is CONICET Researcher and was supported by PIP 0113 (CONICET),
PICT-2017-2865 (ANPCyT), and PIBAA 1275 (CONICET). FG was also supported by
grant PID2019-105510GB-C32/AEI/10.13039/501100011033 from the Agencia Estatal
de Investigación of the Spanish Ministerio de Ciencia, Innovación y
Universidades, and by Consejería de Economía, Innovación, Ciencia y Empleo of
Junta de Andalucía as research group FQM-322, as well as FEDER funds. D.A.
acknowledges support from the Royal Society. V.B-R. is Correspondent
Researcher of CONICET, Argentina, at the IAR. GER was supported by grant PIP
11220200100554CO (CONICET). This work received financial support from the
State Agency for Research of the Spanish Ministry of Science and Innovation
under grant PID2019-105510GB-C31/AEI/10.13039/501100011033 and through the
”Unit of Excellence María de Maeztu 2020-2023” award to the Institute of
Cosmos Sciences (CEX2019-000918-M). This work also made use of the softwares
ds9 (Joye & Mandel, 2003) and Matplotlib (Hunter, 2007).

## References

  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Benaglia & Romero (2003) Benaglia, P. & Romero, G. E. 2003, A&A, 399, 1121
  * Bloot et al. (2022) Bloot, S., Callingham, J. R., & Marcote, B. 2022, MNRAS, 509, 475
  * Callingham et al. (2020) Callingham, J. R., Crowther, P. A., Williams, P. M., et al. 2020, MNRAS, 495, 3323
  * Callingham et al. (2019) Callingham, J. R., Tuthill, P. G., Pope, B. J. S., et al. 2019, Nat. Ast., 3, 82
  * Cappa et al. (2004) Cappa, C., Goss, W. M., & van der Hucht, K. A. 2004, AJ, 127, 2885
  * Crowther (2007) Crowther, P. A. 2007, ARA&A, 45, 177
  * De Becker (2018) De Becker, M. 2018, A&A, 620, A144
  * De Becker & Raucq (2013) De Becker, M. & Raucq, F. 2013, A&A, 558, A28
  * del Palacio et al. (2022) del Palacio, S., Benaglia, P., De Becker, M., Bosch-Ramon, V., & Romero, G. E. 2022, PASA, 39, e004
  * del Palacio et al. (2016) del Palacio, S., Bosch-Ramon, V., Romero, G. E., & Benaglia, P. 2016, A&A, 591, A139
  * del Palacio et al. (2020) del Palacio, S., García, F., Altamirano, D., et al. 2020, MNRAS, 494, 6043
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973
  * Eichler & Usov (1993) Eichler, D. & Usov, V. 1993, ApJ, 402, 271
  * Funk et al. (2013) Funk, S., Hinton, J. A., & CTA Consortium. 2013, Astroparticle Physics, 43, 348
  * Gabici et al. (2019) Gabici, S., Evoli, C., Gaggero, D., et al. 2019, International Journal of Modern Physics D, 28, 1930022
  * H. E. S. S. Coll. et al. (2020) H. E. S. S. Coll., Abdalla, H., Adam, R., et al. 2020, A&A, 635, A167
  * Hamaguchi et al. (2018) Hamaguchi, K., Corcoran, M. F., Pittard, J. M., et al. 2018, Nat. Ast., 2, 731
  * Hamann et al. (2019) Hamann, W. R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57
  * Han et al. (2020) Han, Y., Tuthill, P. G., Lau, R. M., et al. 2020, MNRAS, 498, 5604
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Joye & Mandel (2003) Joye, W. A. & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Leitherer et al. (1995) Leitherer, C., Chapman, J. M., & Koribalski, B. 1995, ApJ, 450, 289
  * Madsen et al. (2017) Madsen, K. K., Beardmore, A. P., Forster, K., et al. 2017, AJ, 153, 2
  * Madsen et al. (2022) Madsen, K. K., Forster, K., Grefenstette, B., Harrison, F. A., & Miyasaka, H. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 034003
  * Madsen et al. (2015) Madsen, K. K., Harrison, F. A., Markwardt, C. B., et al. 2015, ApJS, 220, 8
  * Marcote et al. (2021) Marcote, B., Callingham, J. R., De Becker, M., et al. 2021, MNRAS, 501, 2478
  * Martí-Devesa & Reimer (2021) Martí-Devesa, G. & Reimer, O. 2021, A&A, 654, A44
  * Martí-Devesa et al. (2020) Martí-Devesa, G., Reimer, O., Li, J., & Torres, D. F. 2020, A&A, 635, A141
  * Martí-Devesa et al. (2022) Martí-Devesa, G., Reimer, O., & Reimer, A. 2022, arXiv e-prints, arXiv:2212.10146
  * Martinez et al. (2022) Martinez, J. R., del Palacio, S., Bosch-Ramon, V., & Romero, G. E. 2022, A&A, 661, A102
  * Merten et al. (2017) Merten, L., Becker Tjus, J., Eichmann, B., & Dettmar, R.-J. 2017, Astroparticle Physics, 90, 75
  * Park et al. (2009) Park, S., Kargaltsev, O., Pavlov, G. G., et al. 2009, ApJ, 695, 431
  * Pittard & Parkin (2010) Pittard, J. M. & Parkin, E. R. 2010, MNRAS, 403, 1657
  * Pittard et al. (2021) Pittard, J. M., Romero, G. E., & Vila, G. S. 2021, MNRAS, 504, 4204
  * Reitberger et al. (2015) Reitberger, K., Reimer, A., Reimer, O., & Takahashi, H. 2015, A&A, 577, A100
  * Sander et al. (2019) Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91
  * Tavani et al. (2009) Tavani, M., Sabatini, S., Pian, E., et al. 2009, ApJ, 698, L142
  * Wik et al. (2014) Wik, D. R., Hornstrup, A., Molendi, S., et al. 2014, ApJ, 792, 48
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Zhekov & Skinner (2000) Zhekov, S. A. & Skinner, S. L. 2000, ApJ, 538, 808

##  Appendix A NuSTAR co-added data

In the top panel of Fig. 7 we show the NuSTAR spectra for both NuSTAR
observations and both cameras independently. As detailed in Sect. 2.3, we co-
added the spectra of both epochs for each camera to compare if there was any
systematic difference between FPMA and FPMB. The result is shown in the middle
panel of Fig. 7666We note that the apparently very small background errorbars
returned by the addspec task are appropriate (FTOOLS helpdesk, priv. comm.)..
Even though both cameras are compatible, there is a small difference between
them at energies $<10$ keV, and for that reason they were not combined in the
detailed analysis carried out in Sect. 3. Nonetheless, we also explored the
result of co-adding both cameras in order to improve the signal-to-noise at
the highest energies (neglecting the error introduced at lower energies). The
result is shown in the bottom panel of Fig. 7. This last spectrum has improved
statistics at energies $>20$ keV and allows us to better see by eye the shape
of the SED. Additionally, the fitting of this co-added dataset retrieved fully
compatible results on the best-fitting parameters to the ones obtained in
Sect. 2.3 without co-adding the cameras.

Figure 7: Apep NuSTAR spectra in the 3–35 keV energy range. The spectra was
further rebinned in XSPEC for clarity. The source spectra is shown with
crosses and the background with open circles. Top panel: spectra for the 2022a
and 2022b observations showing each camera separately. Middle panel: spectrum
from the co-added observations from both epochs (for each camera separately).
Bottom panel: spectrum from the co-added observations (both epochs, both
cameras).

##  Appendix B NuSTAR background

As discussed in Sect. 3, the selection of the background extraction region can
have a mild impact in the high-energy spectrum of Apep derived with NuSTAR. We
explore this in more detail to make sure that our results are robust.

In Fig. 8 (top panel) we show an image in the 20–35 keV range of the Apep
field of view for the first observation with FPMB, in which Apep is clearly
detected above the local background. We note that other high-intensity
background regions appear in this image, but these are not consistent in both
observations and FPM cameras; on the contrary, the count excess at the
position of Apep is confirmed in all cases, strongly supporting its detection.
In this image we also mark the different background regions analysed here,
labelled as bkg1–bkg5. The corresponding background-subtracted spectra of Apep
and the background spectra are also shown in Fig. 8 (bottom panel).

  * •

bkg1: This elliptical region is our preferred choice as it is within the same
chip as the source, which is the usual recommendation for analysing NuSTAR
data, and it is sufficiently far from Apep and other sources in the field.
This background region leads to a power-law component with a significance of
$>3.91\sigma$ ($<0.01\%$ of by-chance detection) with a flux of
$F_{\mathrm{10-30\,keV}}=4.8^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2.

  * •

bkg2: An annulus region centred at the position of Apep and with inner radius
of 135″and outer radius of 170″. We note that $\approx 90\%$ of the photons
from Apep should be contained within a 130″from the source; nonetheless, the
presence of weak Fe line emission in the bkg spectrum (seen in Fig. 8, bottom
panel) suggests a non-negligible contamination from the source (at least at
energies below 10 keV). This background region leads to a power-law component
with a significance of $2.12\sigma$ ($3.4\%$ of by-chance detection) with a
flux of $F_{\mathrm{10-30\,keV}}=2.3^{+1.1}_{-1.2}\times 10^{-13}$ erg s-1
cm-2.

  * •

bkg3: This is an elliptical region within the chip with the highest background
level (also appreciable in Fig. 3). This background region leads to a power-
law component with a significance of $1.66\sigma$ ($9.7\%$ of by-chance
detection) with a flux of $F_{\mathrm{10-30\,keV}}=1.3^{+1.1}_{-1.0}\times
10^{-13}$ erg s-1 cm-2.

  * •

bkg4: This is a circular region in another . This background region leads to a
power-law component with a significance of $3.4\sigma$ ($0.07\%$ of by-chance
detection) with a flux of $F_{\mathrm{10-30\,keV}}=3.5^{+1.0}_{-1.1}\times
10^{-13}$ erg s-1 cm-2.

  * •

bkg5: This elliptical region is also within the same chip as the source, but
we note that it is likely contaminated by the hard-spectra supernova remnant
G330.2+1.0 (Park et al. 2009) shown in Fig. 8. This background region leads to
a power-law component with a significance of $>3.85\sigma$ ($0.012\%$ of by-
chance detection) with a flux of
$F_{\mathrm{10-30\,keV}}=4.3^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2.

  * •

At last, we considered all of the previous regions altogether in order to have
an estimate of an average background. This leads to a power-law component with
a significance of $2.53\sigma$ ($1.1\%$ of by-chance detection) with a flux of
$F_{\mathrm{10-30\,keV}}=2.8^{+1.1}_{-1.2}\times 10^{-13}$ erg s-1 cm-2.

We emphasise that the alternative backgrounds bkg2–bkg4 analysed here do not
lie (at least completely) within the same chip as the source, and so they
should be taken with caution. Nonetheless, from this experiment we conclude
that the presence of a power-law component is statistically supported for all
background regions considered. We highlight that the significance of the
detection is lowest for bkg3 ($\approx 1.66\sigma$), while for the remaining
regions the significance is $>2\sigma$. In particular, the significance of the
detection for the backgrounds within the same chip as the source (bkg1 and
bkg5) is $>3.8\sigma$. At last, the flux of the power-law component fluctuates
between $F_{\mathrm{10-30\,keV}}\sim(1.3$–$4.8)\times 10^{-13}$ erg s-1 cm-2,
depending on the background.

Figure 8: Top: NuSTAR image in the 20–35 keV energy range for FPMB (2022a). We
mark the position of Apep with a 30″-circle, the background extraction region
used in the analysis (bkg1), and the alternative backgrounds used as
independent checks (bkg2–bkg5). We also show in cyan contours the emission
from the supernova remnant G330.2+1.0 obtained in the 2.5–8 keV range with
XMM-Newton. The image is in linear scale (instead of logarithmic scale as in
Fig. 3) and it was smoothed in ds9 using a Gaussian with Radius=5. Bottom:
Background-subtracted spectra (crosses) and background spectra (open circles)
of Apep for the different background regions considered in this Appendix B. In
all cases the whole dataset was co-added and rebinned to increase the S/N and
ease the visual comparison.

Complementarily, we followed a different approach to deal with uncertainties
in the background. In this case, a background model was introduced and fitted
using nuskybgd777https://github.com/achronal/nuskybgd-py (Wik et al. 2014).
The advantage of this approach is that the underlying background spectral
model is well-known, thus limiting the uncertainties to statistical ones. The
total spectrum for the source region is extracted and then fitted with a fixed
background model for that location. This is repeated for each observation and
camera separately, using the five background regions discussed previously to
have a good sampling of the background. The result is shown in Fig. 9. In this
case the presence of a power-law component is again favoured, with a C-stat
that diminishes from 682.4/589 to 666.3/588 when a power-law component is
introduced. Similarly as discussed in Sect. 3, the overall decrease in C-stat
is due to the better fit of the NuSTAR data, as the C-stat for the XMM-Newton
cameras actually increases slightly. This power-law component has a flux of
$F_{\mathrm{10-30\,keV}}=3.9^{+1.0}_{-1.2}\times 10^{-13}$ erg s-1 cm-2. The
90% confidence interval for the flux is
$F_{\mathrm{10-30\,keV}}=(2.3$–$5.6)\times 10^{-13}$ erg s-1 cm-2, probing
consistent results with those derived from fitting the background-subtracted
spectra.

Figure 9: NuSTAR spectra from the region of Apep for both observations and
cameras. This is not background-subtracted, but instead a background model was
fitted using nuskybgd; the multiple components of this background model are
shown in thin lines. The source model is TBabs*vphabs*(vpshock+po).

##  Appendix C Non-thermal emission model

Here we present a review of the multi-zone model used to calculate the non-
thermal radiation from Apep. This model is suitable for this CWB as the stars
are separated by tens of AU, which means that the shocks in the WCR are
essentially adiabatic and quasi-stationary (del Palacio et al. 2022). The WCR
structure is treated as an axi-symmetric surface under a thin shock
approximation (i.e., the fluid is considered homogeneous along the direction
perpendicular to the shock normal). The thermodynamical properties along the
WCR (density, magnetic field intensity, etc.) are calculated with semi-
analytical prescriptions based on mass and energy conservation of the fluid
elements (Martinez et al. 2022), which is slightly more precise than the
original assumption of Rankine-Hugoniot jump conditions in del Palacio et al.
(2016). The only free parameter in the hydrodynamical model is the ratio
between the thermal pressure (set by the properties of the stellar wind) and
the magnetic pressure in the WCR, $\eta_{B}$.

Relativistic particles accelerate when a fluid line from a stellar wind enters
the WCR. The relativistic particle distribution injected at a given position
in the WCR is a power law with the spectral index given by the radio
observations ($p=2.4$). This distribution is normalised such that the injected
power is a fraction $f_{\mathrm{NT}}$ of the total power available for
particle acceleration ($L_{\mathrm{inj,\perp}}$). This power is distributed in
electrons and protons as
$f_{\mathrm{NT}}=f_{\mathrm{NT,e}}+f_{\mathrm{NT,p}}$. One common
parameterisation is $f_{\mathrm{NT,e}}=K_{\mathrm{e,p}}f_{\mathrm{NT}}$, with
$K_{\mathrm{e,p}}\sim 0.01-0.1$ (see Merten et al. 2017, for a discussion of
uncertainties in this value). Upon injection, particles are attached to the
fluid lines via the magnetic fields and flow together with the shocked fluid.
As they stream, particles cool down due to different processes. The
approximate solution used for the transport equation is given in Appendix A of
del Palacio et al. (2022).

Finally, the non-thermal particles produce broadband radiation. The most
relevant emission mechanisms are synchrotron and IC for electrons, which
dominate in the radio and high-energy domain (hard X-rays and $\gamma$ rays),
and proton-proton collision for protons, which can contribute to the
$\gamma$-ray flux (though this hadronic contribution is subdominant if
$K_{\mathrm{ep}}>0.01$). This emission is then corrected for absorption,
namely free-free absorption in the ionised stellar winds for low-frequency
radio emission, and $\gamma$-$\gamma$ absorption in the stellar radiation
fields for high-energy $\gamma$ rays.