The high-energy emission from HD 93129A near periastron

We conducted an observational campaign towards one of the most massive and luminous colliding wind binaries in the Galaxy, HD~93129A, close to its periastron passage in 2018. During this time the source was predicted to be in its maximum of high-energy emission. Here we present our data analysis from the X-ray satellites \textit{Chandra} and \textit{NuSTAR} and the $\gamma$-ray satellite \textit{AGILE}. High-energy emission coincident with HD~93129A was detected in the X-ray band up to $\sim$18~keV, whereas in the $\gamma$-ray band only upper limits were obtained. We interpret the derived fluxes using a non-thermal radiative model for the wind-collision region. We establish a conservative upper limit for the fraction of the wind kinetic power that is converted into relativistic electron acceleration, $f_\mathrm{NT,e}<0.02$. In addition, we set a lower limit for the magnetic field in the wind-collision region as $B_\mathrm{WCR}>0.3$~G. We also argue a putative interpretation of the emission from which we estimate $f_\mathrm{NT,e} \approx 0.006$ and $B_\mathrm{WCR} \approx 0.5$~G. We conclude that multi-wavelength, dedicated observing campaigns during carefully selected epochs are a powerful tool for characterising the relativistic particle content and magnetic field intensity in colliding wind binaries.


INTRODUCTION
The role of massive stars as accelerators of Galactic cosmic rays has been gaining more interest in the recent years (e.g. Seo et al. 2018;Aharonian et al. 2019;Prajapati et al. 2019). Unfortunately, it is not possible to directly assess the efficiency of cosmic ray acceleration from observations. The study of massive colliding-wind binaries (CWBs) is key in understanding the non-thermal physics taking place in systems harbouring massive stars (De Becker et al. 2017). The signature of relativistic particles is a non-thermal spec-E-mail: sdelpalacio@iar.unlp.edu.ar (IAR) trum, typically a power-law with additional features such as a spectral break or an exponential cutoff. In the radio band, non-thermal emission can be produced by relativistic electrons interacting with the magnetic field (Eichler & Usov 1993). The detection of this synchrotron radiation in radio observations of CWBs have provided conclusive evidence that at least dozens of CWBs in the Galaxy are capable of accelerating cosmic rays (De Becker & Raucq 2013). The most accepted scenario is that relativistic particles accelerate via diffusive shock acceleration in the strong shocks of the windcollision region (WCR; Pittard & Dougherty 2006). However, the cosmic ray acceleration efficiency cannot be proven from radio data alone (e.g. del Palacio et al. 2016, an references therein). Therefore, additional information of the high-energy spectrum of CWBs is required in order to make progress in the characterisation of their non-thermal physics.
The idea that CWBs can produce significant non-thermal emission at high energies (X-rays and γ-rays) is supported by the detection of such a radiation from the exceptional binary η-Carinae (Tavani et al. 2009b;Reitberger et al. 2015;Hamaguchi et al. 2018;H. E. S. S. Collaboration et al. 2020), although this fascinating object is hardly representative for CWBs in general. Another CWB, γ 2 Vel, has also been likely detected at γ-rays (Pshirkov 2016;Martí-Devesa et al. 2020). However, none of these two systems has been detected as a non-thermal radio emitter, most likely due to synchrotron emission being self-absorbed in the system (e.g. Benaglia et al. 2019). If more CWBs are confirmed as non-thermal emitters at high energies, it will be possible to assess the role of CWBs as cosmic ray injectors. Moreover, it will give compelling support to the idea that massive binary systems are frequently γray emitters (Benaglia & Romero 2003;De Becker et al. 2017).
In this work, we focus on the study of the high-energy emission from the extreme binary HD 93129A. This system, located in the core of the young star cluster Trumpler 14, is the most massive un-evolved binary in the solar vicinity. It is made up by two O2 stars that are among the earliest, hottest, most luminous, and with highest mass-loss rates in the Galaxy (Walborn et al. 2002;Maíz Apellániz et al. 2008). The most recent published ephemeris (Maíz Apellániz et al. 2017) yielded a binary period of ∼ 120 years, a periastron passage in 2018.54 +0.54 −0.32 , and an orbital inclination i ∼ 99 • (i.e., nearly edge-on). We have kept obtaining VLTI astrometry of the HD 93 129 Aa,Ab visual pair after the orbital solution that was published in Maíz Apellániz et al. (2017). The new data includes three post-periastron observations that help to significantly reduce the uncertainties in the orbital determination. Regarding the most relevant parameters for this study, the time of periastron passage is now tightly constrained to 2018.70 +0.22 −0.12 (i.e., a 1-σ uncertainty of just 2 months), and the distance at periastron is also well constrained within a ∼ 5% uncertainty, namely as 7.91 ± 0.42 mas (18.6 ± 1.0 AU). In addition, the non-thermal synchrotron radiation from the WCR has been resolved in radio by Benaglia et al. (2015), allowing for a partial characterisation of the relativistic particle energy distribution, and confirming that particle acceleration occurs in this binary.
We conducted our high-energy observational campaign towards this system near its periastron passage. The high-energy nonthermal emission from the WCR, produced by inverse Compton scattering of stellar photons, was expected to peak towards periastron due to the increased number of target stellar photons. The non-thermal emission from this system was predicted by del Palacio et al. (2016) using a broadband radiative model for the WCR that takes into account (i) the evolution of the accelerated particles streaming along the shocked region, (ii) the emission by different radiative processes, and (iii) the attenuation of the emission propagating through the local matter and radiation fields. On the basis of their analysis, del Palacio et al. (2016) suggested that HD 93129A was a promising candidate for detecting high-energy emission close to its periastron passage. A multi-wavelength observational campaign was carried out to get a complete picture of this event. In this work, we present the analysis of high-energy observations performed with Chandra (soft X-rays, 0.3-10 keV), NuSTAR (hard X-rays, 3-79 keV) and AGILE (γ-rays, 50 MeV-300 GeV) close to periastron.
The paper is organised as follows. The description of the observations is given in Sect. 2, and the procedure adopted to reduce and analyse each data set in Sect. 3. We present the results from this analysis in Sect 4. In Sect. 5 we discuss the constraints of our results in the context of theoretical models, and we present a summary and concluding remarks in Sect 6.

OBSERVATIONS
We observed the binary HD 93129A close to its periastron passage in 2018. A summary of the analysed observations is presented in Table 1. In the following subsections we describe them in more detail.

Chandra observations
The Chandra X-Ray Observatory is unique in terms of the high angular resolution that it reaches in the X-ray band. The ACIS instrument on board of the satellite is capable of focusing X-rays with energies in the range 0.3-10 keV with subarcsecond resolution (Weisskopf et al. 2002a). This is particularly useful for our study as it allows us to quantify the contamination from background sources such as the Trumpler 14 association. An image of the field of vew of HD 93129A is shown in Fig. 1.
Since we are interested in studying the continuum emission from HD 93129A (i.e., not the detailed spectral lines), we restrict our analysis to on-axis observations without gratings. We found one archival observation (Obs. id. 4495) from 2004 with 58 ks exposure taken with ACIS-I. In addition, we have observations from June and August 2018 (Obs. id. 20152 and 20153, respectively), both with 25 ks exposure taken with ACIS-S.
We reduced the Chandra data using CIAO v.4.11 (Fruscione et al. 2006). For all the data sets we used the CIAO task chandra_repro to reprocess the data using the latest calibration files available (CALDB v.4.8.4.1).
The source spectra were extracted using the task specextract with the options appropriate for point source analysis (weight=no and correctpsf=yes). When considering a larger region (10 or 30 ) to estimate contamination by nearby sources, weight was set to yes. We found little variation between both observations from 2018 (less than 2% in the integrated flux above 3 keV), so the task combine_spectra was used to combine their spectra. Finally, we used the task ftgrouppha from HEASOFT v.6.26 to group the spectra with the optimal binning scheme by Kaastra & Bleeker (2016).

NuSTAR observations
The NuSTAR X-ray observatory, launched in 2012, counts with two co-aligned hard X-ray grazing incidence telescopes labelled by their focal plane modules FPMA and FPMB. These instruments are 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).
NuSTAR observations with 30 ks exposure were performed within a week of the Chandra observations. This makes both data sets almost simultaneous considering the expected variability timescale. The nupipeline task was used to create level 2 data products. We used the option saacalc=2 saamode=optimized tentacle=yes to filter high background epochs. This leads to

HD 93129B
HD 93129A Figure 1. Chandra RGB image from obs. 20152 (red is 0.5-1.2 keV, green is 1.2-2.5 keV, and blue is 2.5-7 keV). Circular regions centred at the position of HD 93129A and of different radii are shown for reference, as well as the selected background extraction region. The position of the source HD 93129B (not part of the system HD 93129A) is also marked for clarity.
< 3% data loss 1 . We used the nuproducts task to create level 3 data products. We present an image (3-11 keV; PI channels 35-235) in Fig. A1, showing the selected source and background extraction regions for each observation.

AGILE observations
We analysed the data collected between 01/01/2018 and 31/12/2018 by the Gamma-Ray Imaging Detector (GRID; Barbiellini et al. 2002;Prest et al. 2003) on board the AGILE satellite (see Tavani et al. 2009a, for a detailed description of the AGILE payload). The AGILE-GRID is sensitive in the 30 MeV-50 GeV energy band. The point spread function (PSF) at 100 MeV and 400 MeV is 4.2 • and 1.2 • (68% containment radius), respectively (Sabatini et al. 2015). We restricted our analysis to photon energies from 100 MeV to 50 GeV. The angular resolution of AGILE in this range is 0.9 • , so that significant contamination from sources close to HD 93129A is expected. Since 2009 AGILE observes in 'spinning' mode, covering a large fraction of the sky with a controlled rotation of the pointing axis. In this observing mode, typical two-day integration-time sensitivity for sources in the Galactic plane and photon energy above 100 MeV is ∼ 10 −6 photons cm −2 s −1 (3σ).
The analysis of the AGILE-GRID data was carried out with the new Build 23 scientific software, FM3.119 calibrated filter and I0025 response matrices. The consolidated archive, available from the ASI Data Center (ASDCSTDk), was analyzed by applying South Atlantic Anomaly event cuts and 80 • Earth albedo filtering. Only incoming γ-ray events with an off-axis angle lower than 60 • were selected for the analysis. Statistical significance and flux determination of the point sources were calculated by using the AGILE multi-source likelihood analysis (MSLA) software (Bulgarelli et al. 2012) based on the Test Statistic (TS) method as formulated by Mattox et al. (1996). This statistical approach provides a detection significance assessment of a γ-ray source by comparing maximumlikelihood values of the null hypothesis (no source in the model) with the alternative hypothesis (point source in the field model). In this work we report 68% confidence level flux upper limits if TS < 9 (detection significance < 3). To estimate the likelihood of a detection, two different considerations were made to account for multiple nearby sources: i) excluding η-Car, and ii) including η-Car. For each of these background models we constructed weekly and monthly binned light-curves.

SPECTRAL FITTING
The X-ray spectrum of a CWB has many physical components produced by the two individual stellar winds and the two shocked winds. The individual stellar winds produce only low energy X-rays (typically < 1 keV, Owocki et al. 1988), whereas the shocked winds can produce hard X-rays above 3 keV via thermal emission or by IC scattering of stellar photons by relativistic electrons (Hamaguchi et al. 2018). These emission components are affected by absorption by the intertellar medium and within the stellar winds (Cohen et al. 2011). This means that a realistic and physically consistent model requires several components, which unfortunately leads to a degeneracy when fitting real data. In practice, the observed spectrum can be well approximated with a two-temperature apec 2 component together with a phabs (Anders & Grevesse 1989) absorption component (Pittard 2010). The apec is an optically thin plasma thermal emission model that has three parameters: the plasma temperature (kT ), normalisation (N), and abundance (A). This model is adequate to represent the X-ray emission from a hot plasma such as  Table 2. those produced by hot stellar winds of massive stars (Feldmeier et al. 1997;Owocki et al. 2013) and CWBs (Stevens et al. 1992;Pittard & Parkin 2010;Hamaguchi et al. 2018). The low temperature apec includes the thermal emission from the individual stellar winds and the outer (colder) regions of the WCR. The hightemperature apec represents the thermal emission from near the apex of the WCR. In addition, a multiplicative phabs model is used as an effective absorption component that takes into account the attenuation in the ISM and the stellar winds. We note that this model, despite being an oversimplification, is sufficient to qualitatively characterise the spectrum of the source (see Sect. 4.1). An additional power-law component can also be introduced to account for putative non-thermal emission. Finally, we stress that this study is focused on the high-energy emission, so a detailed spectral fitting of low energy photons (< 3 keV) will be addressed in a separate work, together with a detailed analysis of the temporal variability of such thermal emission.
The spectral analysis was done using the software XSPEC v.12.10.1f (Arnaud 1996;Dorman & Arnaud 2001). The C-stat minimisation approach is used to check the adequacy of the fit, given that the source is not very bright in hard X-rays. Error-bars are calculated at 1σ for all cases.

Chandra
We aim to characterise the X-ray spectrum of HD 93129A and its environment using the Chandra data. In Fig. 1 we show the Chandra field of view centred at the position of HD 93129A. The source is resolved thanks to Chandra's high angular resolution (≈ 0.5 on axis). Several sources are detected within 10 from HD 93129A. We therefore extract the spectrum from HD 93129A using a 1 circular region so that the background/foreground contamination is minimal. The PSF is such that it ensures that almost 100% of the photons with E < 1 keV, and ∼ 85% for those with E ∼ 10 keV, are captured (Weisskopf et al. 2002b).
We compare the spectra from the 2004 and 2018 observations in Fig. 2. The source was detected up to 8 keV. We fitted both spec- tra simultaneously with a model [phabs*(vapec 1 +vapec 2 )]. We considered the same element abundances for both vapec components. We defined a common abundance A for all elements except Si and S, which had prominent emission lines in the spectra, and Fe, which had a significantly lower abundance than the rest. We tried to minimise the amount of free parameters to fit while still obtaining a good spectral fit. For that reason, we tied the abundances of A and Fe between epochs, and allowed only S and Si abundances to vary. The temperature of the colder vapec component, T 1 , was also tied between epochs. The values of N H , T 2 , and normalisations norm 1 and norm 2 , were fitted independently for each epoch. The best fit values are shown in Table 2. We found that between 2004 and 2018 the observed flux increased by more than 50% in the 0.5-3 keV range and 100% in the 3-8 keV range. In order to compare the Chandra and NuSTAR spectra, we need to use a 30 extraction region for both. For this reason, we compared how selecting a 30 extraction region affected the Chandra spectrum from 2018. We fitted the Chandra spectrum in the 3.0-7.6 keV energy range (the maximum energy is set by the larger background contamination dominating the spectrum, see Fig. A2) for two different extraction regions of radii 1 and 30 . We tested whether a complex model was required to fit the spectrum in this energy range. For this we checked that in this energy range the contribution from the low temperature vapec and from the absorption phabs were negligible. We therefore fitted the spectra with a single high temperature apec component. This component has kT ≈ 2.35 ± 0.25 keV for both extraction regions. We conclude that there is no significant spectral shape difference in the 3-7.6 keV energy range induced by the increased background contamination in the 2018 Chandra spectra extracted from a 30 region with respect to the one from a 1 region.
In addition, we also checked the 2004-2018 variability considering a 30 extraction region. The spectral fitting parameters are summarised in Table 2. In this case, the flux between the two epochs increased only in 26% in the 0.5-3 keV band and 33% in the 3-8 keV band; similar increase factors were obtained for a 10 region as well (Fig. A2). Interestingly, we can also estimate the average flux of the surrounding sources within an annulus given by the difference between the 30 -and the 1 -extraction regions. We obtain that the flux variations for the nearby (background) sources between 2004 and 2018 are < 10% in both of these energy bands (this can be appreciated more clearly in Fig. A3). Therefore, the observed flux variability within the 30 region seems to be governed by HD 93129A. We took into account this information when interpreting the NuSTAR data in which the binary is not resolved from the surrounding sources (Fig. 1).

NuSTAR
The 3-18 keV fluxes between both FPMA and FPMB cameras in each observation were consistent within less than 5%. We therefore co-added both cameras using addspec. In Fig. 3 we compare the source and background spectra grouped with ftgroup. The source is detected above background with a high significance at ≈ 13 keV (4σ in the 11.3-15.1 keV energy range), and with a lower significance at ≈ 18 keV (1.7σ in the 15.2-17.8 keV energy range). We confirmed this result by selecting different background extraction regions, although the total flux in 3-10 keV can vary within 10% depending on the selected background. Moreover, the fluxes between the two observations differ only by 12%, which is almost within 1-σ level and can be attributed to calibration uncertainties (up to 10%; Madsen et al. 2015) and a different background level (see top panel from Fig. 3). We therefore assume that the whole NuSTAR data set is comparable and co-add both observing epochs with addspec.
In Fig. 3 we plotted together the X-ray data from Chandra and NuSTAR. The spectra seemed to match well up to 6-7 keV.
To check this, we fitted the NuSTAR spectrum below 7 keV with a model, and the Chandra spectrum from the 30 -region with the same model times a constant. The fitted value of the constant is 1.036 +0.056 −0.050 , which is fully consistent with not needing a re-scaling for the data sets (i.e., setting the constant to one). We found that the model fitted using only the Chandra data matches very well the NuSTAR spectrum up to 6-7 keV. However, the Chandra model underestimates the emission above 7 keV. This could be simply due to a poorly constrained high-energy spectrum in the Chandra data resulting in a bad extrapolation of the fitted model. Regardless of its origin, the fact that this model falls below the observed flux in the NuSTAR spectrum suggests the presence of an additional component, either thermal and hotter or non-thermal.

AGILE
The result from the analysis of the AGILE data is consistent with an upper limit to the flux in 50 MeV-100 GeV of 2 × 10 −7 cm −2 s −1 .  Figure 3. Top: NuSTAR source and background spectra for each observation (first observation is black, second one is red) in the 3-40 keV energy range. The cameras at each observation epoch were co-added with addspec, grouped with ftgroup and rebinned in XSPEC with rebin=4. Bottom: NuSTAR source and background spectra obtained by co-adding both observations, shown in the 3-20 keV energy range. The model overlayed is the one that fits the Chandra spectrum in the 3-7.5 keV range, with the residuals shown in the bottom panel. The model under-predicts the emission above 8 keV, which can be explained by the presence of either a hotter component or a non-thermal component. This is confirmed by the light-curve obtained for HD 93129A using both a weekly and a monthly binning and including the background emission from multiple sources nearby. An example light curve is shown in Fig. A4. We note that there is a low significance detection of emission near April 2018 at a 4-σ level if emission from η-Car is not considered. However, it is not likely related with HD 93129A, as discussed in Sect. 5.

DISCUSSION
In Sect. 4 we presented the results from the observations with Chandra, NuSTAR and AGILE. The source HD 93129A was clearly detected in the X-ray domain (0.5 − 18 keV range), whereas in the γ-ray domain no significant emission was detected.

Observations versus predictions
From the analysis of the Chandra data we can conclude that there is flux variability between 2004 and 2018. This variability is more clearly seen when a small 1 extraction region is considered. For larger extraction regions with size ≥ 10 the background level becomes comparable to the emission from HD 93129A, and therefore the estimated variability intrinsic to HD 93129A is unreliable.
The significant long-term variability detected is qualitatively consistent with the expected orbital modulation of the X-ray emission. The theoretical expectation for a thermal X-ray emission produced by the wind-wind interaction in adiabatic conditions is that the unabsorbed flux varies as ∝ 1/D (Stevens et al. 1992). However, testing this from observations requires to fit more complex spectral models, particularly in the soft X-ray band, which is the most affected by absorption. Moreover, converting the date of observation to system separation D requires a precise knowledge of the orbital parameters. For this reason we did not aim to present a quantitative analysis of the flux variability in this work.
The spectrum above 12 keV is poorly constrained. Additional uncertainties arise given that NuSTAR is unable to resolve the emission coming from a 30 region centred at the position of HD 93129A. Thus, it is impossible to determine the actual fraction of the observed flux that corresponds to HD 93129A. Nonetheless, Chandra's high angular resolution is capable of disentangling the different contributions up to 8 keV. We find that approximately half of the emission in the 3-8 keV range comes from the CWB. The spectral shape of the background sources is consistent with the one from HD 93129A. Moreover, the Γ value we obtained from a power-law fit is consistent with the value of Γ 2 found for η-Car by Hamaguchi et al. (2018). This gives further support to the idea that the CWB HD 93129A is (partly) responsible for the observed hard X-ray emission. If we assume this extrapolation remains valid above 8 keV, we expect that roughly 50% of the observed flux by NuSTAR comes from HD 93129A as well.
Whether the measured 8-18 keV flux is an upper limit or a good estimate of the emission from HD 93129A depends on the assumptions made.
(i) First, we can place a hard limit to the maximum flux for the power-law IC component by considering that none of the observed flux in the 8-18 keV range comes from HD 93129A. In this case we obtain F 8−18 ≤ 2.6 × 10 −13 erg s −1 cm −2 (considering the 1σ upper-limit).
(ii) A similar but less restrictive constraint is obtained by stating that the flux from the power-law IC component cannot be higher than the flux from the fitted power-law component. In this case we obtain F 8−18 ≤ 2.2 × 10 −13 erg s −1 cm −2 (again assuming the 1σ upper limit).
(iii) The analysis in Sect. 4.1 suggests that roughly ∼ 50% of the observed hard X-ray flux comes from HD 93129A. In this case we can estimate the flux from the power law component as The most conservative is the first.
Assuming the hard X-ray spectrum is thermal (with solar abundances), the maximum plasma temperature for an adiabatic shock in the WCR is kT = 1.17v 2 w,8 keV (Stevens et al. 1992), with v w,8 the wind velocity in units of 10 8 cm s −1 . For a wind speed of v w,8 ≈ 3 (Cohen et al. 2011), maximum (i.e. close to the apex) shock temperatures kT ≈ 10 keV are expected. This is within the poorly constrained values obtained in Sect. 4.2.
The CWB HD 93129A does not display a significant γ-ray activity as to be detected by AGILE. We note that there is a low significance detection of emission near April 2018 at a 4-σ level if emission from η-Car is not considered. Despite η-Car being a known γ-ray source located at just ∼ 0.2 • from HD 93129A (i.e., within the same beam for AGILE), during this epoch η-Car was expected to be in a low γ-ray emission state (Balbo & Walter 2017;White et al. 2019). Thus, it is not likely that the detected flux of 8 × 10 −7 cm −2 s −1 is due to η-Car. Nonetheless, it is less likely that it comes from HD 93129A, as its γ-ray flux was expected to be increasing at that epoch and no emission was found afterwards. Moreover, under the constraints imposed by not exceeding the observed hard X-ray flux, the γ-ray flux predicted by our non-thermal radiative model is < 5 × 10 −12 cm −2 s −1 , which is more than four orders of magnitude below the detected γ-ray flux (more details in Sect. 5.2). We also considered a possible hardening in the nonthermal electron distribution as in del Palacio et al. (2016), which can enhance the γ-ray emission in a factor ten, but still it is not possible to account for the detected γ-ray flux even considering uncertainties in the system or model parameters. Finally, this emission is detected only at a 2-σ level in the weekly-binned light-curves, and is thus not considered to be significant. We conclude that this emission is either a statistical fluctuation or due to a variable background.

Theoretical modelling
We use an updated version of the non-thermal emission code presented in del Palacio et al. (2016) to calculate the predicted X-ray flux from HD 93129A. The modifications introduced in the model are described in Appendix B. We adopt similar parameters for the stellar winds; specifically, mass-loss ratesṀ 1 = 10 −5 M yr −1 andṀ 2 = 6 × 10 −6 M yr −1 , and wind terminal velocities v ∞,1 = 3200 km s −1 and v ∞,2 = 2800 km s −1 . This model has two free parameters: the ratio between the magnetic field pressure to thermal pressure in the WCR, η B , and the fraction of the available power at the shocks that is converted into relativistic electrons, f NT,e . The available power for particle acceleration is the wind kinetic power injected perpendicularly to the WCR, which is calculated consistently in the model. It is possible to tie these two parameters by modelling the synchrotron component revealed by radio data by f NT,e B WCR Figure 5. Estimated values of f NT,e and B WCR in order to obtain a certain X-ray flux in the 8-18 keV. We mark with a dark grey line the strict (most conservative) upper limit for the non-thermal X-ray flux, which assumes that none of the flux observed with NuSTAR comes from HD 93129A. The red region corresponds to values ruled out by the observations. We also show with a grey line the estimated flux value assuming that the relative contribution of HD 93129A to the X-ray flux with respect to the background remains constant above 8 keV (see Sect. 5.1). Benaglia et al. (2015). In this case the relation is f NT,e B 2 = constant (del Palacio et al. 2016), which we assume holds along the orbit. However, it is not possible to break the degeneracy between these two parameters from radio data alone. Theoretical estimates of the non-thermal X-ray flux combined with the X-ray fluxes observed in the 8-18 keV allow us to constrain f NT,e . This is the only parameter that has a significant impact in determining the non-thermal X-ray flux, which is F X ∝ f NT,e . We adopt a periastron distance of D = 19 AU (see Sect. 1) and compute the model 8-18 keV flux for reasonable values of f NT,e . It is also possible to estimate the magnetic field in the apex of the WCR, B WCR , by assuming that the relation between f NT,e and η B holds along the orbit. A constant value of η B is consistent with B WCR scaling as 1/D. This, in turn, is consistent with a constant (or null) B-field amplification factor in the WCR. Taking into account that f NT,e B 2 = constant and F X ∝ f NT,e , we get the relation B ∝ F −1/2 X . In Fig. 5 we show at the left vertical axis the required value of f NT,e in order to reach a certain 8-18 keV flux, and the corresponding value of B WCR at the right vertical axis.
We can interpret Fig. 5 in two ways: (i) considering the upper limit to the power-law component in X-rays, from which we obtain an upper limit of f NT,e < 0.02 and a lower limit of B WCR > 0.3 G (η B > 0.01); and (ii) considering a detection at the estimated value of F 8−18 = 1.1 × 10 −13 erg s −1 cm −2 , from which we estimate f NT,e ≈ 0.006 and B WCR ≈ 0.5 G (η B ≈ 0.02). For this latter case, we calculate the SED for two scenarios, one in which the injected particle energy distribution has a constant spectral index p = 3.2, and one more favourable for γ-ray production in which the injected distribution hardens at high energies (> 100 MeV for electrons; see del Palacio et al. 2016). In Fig. 6 we show the modelled SED together with the observational data for the periastron passage. Unfortunately, further constraints to the non-thermal particle population cannot be placed using the AGILE upper limit as it is much higher than the γ-ray flux predicted in the most favourable scenario. In addition, as already discussed by del Palacio et al. (2016), it is difficult to observe this system in the radio band close to periastron passage: despite the intrinsic synchrotron flux increases close to pe- riastron as it scales with B WCR , the absorption of the low-frequency photons in the stellar winds also boosts during this epoch, resulting in a reduced flux below 10 GHz. It is possible to relate the inferred magnetic field intensity in the WCR to that in the stellar surface. At large distances from the star (r R ) the stellar magnetic field is toroidal and drops as ∝ r −1 . Following del Palacio et al. (2016, and references therein), we can express B ≈ 2.5B WCR (r/R ). This expression considers an Alfvén radius r A ∼ R , a stellar rotation velocity v rot ≈ 0.1v ∞ , and also takes into account the adiabatic compression of magnetic field lines in the WCR. We set r equal to the periastron separation D and allow, again, for two interpretations: (i) on the one hand, we obtain B > 130 G for the lower limit of B WCR > 0.3 G; (ii) on the other hand, we get B ≈ 200 G for the estimated value of B WCR ≈ 0.5 G. Such a modest B-field strength would not be sufficient to confine the winds (i.e. r A ≤ R , as assumed). This is also consistent with insignificant emission of X-rays above ∼ 1 keV by the individual stellar winds, which would require them to be magnetically confined (e.g. Ud-Doula et al. 2009).
We clarify that magnetic field amplification processes in the WCR (e.g. Bell 2004;Drury & Downes 2012;del Valle et al. 2016, for studies in the context of supernova remnants) were not considered nor computed in the calculations. If the measured magnetic field in the WCR actually results from some B-field amplification, this would translate into a stellar value lower than the one determined from the sole toroidal geometric dilution and adiabatic compression.

CONCLUSIONS
We report the results from an observational campaign on the extreme colliding wind binary HD 93129A close to its periastron passage in 2018. The obtained optical data allowed us to determine the periastron epoch precisely (results will be presented in a forthcoming work), which we used to trigger our campaign in the high-energy domain. We observed HD 93129A with Chandra, NuSTAR, and AGILE, covering a wide range of energies in the Xray and γ-ray domain. We conclude that there is no significant γray emission from this system and the upper limit of the flux in the 50 MeV-100 GeV range, F < 3 × 10 −7 cm −2 s −1 , is unconstraining for non-thermal emission models. Our main observational result is the detection of emission in the hard X-ray band with NuSTAR consistent with being (partially) produced by HD 93129A. It is not possible at this time to assess to what extent background contamination is accountable for this emission. Future observations of the source might reveal a persistent flux of the same value, which would indicate that this emission is produced by nearby sources, or a diminished flux, which would confirm HD 93129A as the responsible for this emission.
We interpret the derived hard X-ray fluxes using the tightly constrained periastron distance from the optical monitoring and the non-thermal radiative model described in del Palacio et al. (2016). This multi-zone model takes into account the relevant physical processes in the wind-collision region. The model has two parameters that can be constrained or estimated by our observations: the fraction of the wind kinetic power injected into the WCR that is converted into relativistic electron acceleration, f NT,e , and the magnetic field in the wind-collision region, B WCR . We present the conclusions for two different interpretations: • Under the very conservative assumption that none of the X-ray flux above 8 keV is produced by HD 93129A, we obtain f NT,e < 0.02. In addition, we estimate the magnetic field in the windcollision region as B WCR > 0.3 G. Neglecting possible magneticfield amplification in the wind-collision region, we derive a lower limit for the surface stellar magnetic field of B > 130 G.
• We consider that ∼ 50% of the 8-18 keV flux produced by a power-law component comes from HD 93129A. In this case we can estimate f NT,e ≈ 0.006. In addition, we get B WCR ≈ 0.5 G, from which we derive an upper limit for the surface stellar magnetic field B ≤ 200 G taking into account possible B-field amplification.
We conclude that multi-wavelength, dedicated observing campaigns during carefully selected epochs is a powerful tool for characterising the relativistic particle content and magnetic field intensity in colliding wind binaries. This, in turn, allows to constrain the value of the magnetic field on the surface of very massive stars. We also highlight the need for more sensitive and higher angularresolution observations in the γ-ray band in order to better characterise the non-thermal particle population in colliding-wind binaries.

APPENDIX A: IMAGES
We present here additional figures that complement our analysis. In Fig. A1 we show a NuSTAR-FPMB image of the field of view for the two observing epochs; FPMA image (not shown) is similar. The background extraction region was selected within the same chip that contained the source. Significant and variable background emission is observed close to the position of HD 93129A.
In Fig. A2 we show the comparison of the Chandra source and background spectra for different source extraction regions of radii 1 , 10 , and 30 . For the 1 source extraction region the background is negligible, whereas for radii ≥ 10 it becomes comparable to the source emission at energies above 7 keV. In addition, in Fig. A3 we compare the combined Chandra spectra of the nearby sources to HD 93129A between 2004 and 2018. The spectra were extracted from an annulus of outer radius 30 and an inner radius of 1 centred at the position of HD 93129A. The combined spectra from these sources are remarkably stable above 2 keV, where the continuum dominates. The apparent discrepancy below 2 keV is probably caused by a difference in the ACIS instrumental response: ACIS-S has a better soft band sensitivity than ACIS-I, but the soft band efficiency of all ACIS sensors has declined recently with contamination on the optical blocking filter. ACIS-S also shows a moderately strong feature around E 2 keV. The flux difference between the two epochs is < 5% in the 0.5-3 keV range and < 10% in the 3-8 keV energy range.
Finally, in Fig. A4 we present the monthly light curves obtained with AGILE for two different background models. When the possible contribution from η-Car is not considered, a detection with a 4-σ significance is obtained in April 2018. However, when η-Car is included as an additional background source the detection is not significant and only upper limits are obtained.

APPENDIX B: NON-THERMAL EMISSION MODEL
Here we present a review of the non-thermal radiation model used. This model is an update to the one presented in del Palacio et al. (2016), suitable for adiabatic and quasi-stationary shocks with a laminar flow. The WCR structure is treated as a 2-dimensional surface under a thin shock approximation. We assume that relativistic particles are accelerated once a fluid line from the stellar wind enters the WCR region. These particles flow together with the shocked fluid which convects the ambient magnetic field. As they stream, particles cool down due to different processes and produce broadband radiation. This emission is corrected for absorption by interaction with the local matter and radiation fields.
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. The normalisation of this distribution is such that the injected power is a fraction f NT of the total power available for particle acceleration (which is only a fraction of the total power of the stellar winds; del Palacio et al. 2016). This power is distributed in electrons and protons as f NT = f NT,e + f NT,p . It is usual to define f NT,e = K e,p f NT , with K e,p ∼ 0.01 − 0.1 (see Merten et al. 2017, for a discussion of uncertainties in this value). For a canonical value of f NT ∼ 0.1, we expect f NT,e ∼ 10 −3 , but with a large uncertainty.
The emission in the radio band is produced by synchrotron emission. This radiation can be significantly attenuated by freefree absorption in the ionised stellar winds. The non-thermal X-ray emission is produced by anisotropic inverse-Compton up-scattering of stellar photons. This process can dominate the γ-ray emission as well, competing with proton-proton inelastic collisions. The γ-ray photons can be absorbed in the stellar radiation field creating secondary electron-positron pairs.
We introduced the following modifications to the model in del Palacio et al. (2016): • The angle ψ such that the observed distance between the stars is D proj = D cos ψ is no longer a free parameter of the model 3 . A value of ψ ∼ 32 • is used in accordance with the most recent orbital ephemeris.
• We considered an increased free-free opacity at radiofrequencies due to clumping in the stellar winds. This effect enhances the opacity by a factor f −1/2 , where f ∼ 0.1 is the volume filling factor of the wind (Muijres et al. 2011). In this case there is no need to adopt a large value of E min,e to reproduce the spectral break at low frequencies as done in del Palacio et al. (2016). Instead, a typical value of E min,e ≈ 1 MeV yields a good fit of the radio spectra.
• The IC spectra is calculated using the expressions by Khangulyan et al. (2014) suitable for black-body-like target photon fields. This reduces the computation time significantly.
• A small correction was introduced in the way particle evolution along stream lines is calculated. This considers variations of cooling times from cell to cell in the emitter (Eq. 16 from del Palacio et al. 2018). This paper has been typeset from a T E X/L A T E X file prepared by the author.   Figure A2. Comparison of the Chandra spectra extracted from a 1 (black), 10 (red), and 30 (blue) region centred at the position of HD 93129A. The crosses represent the source emission and the circles (with errorbars) the background. The spectra are rebinned (setpl rebin 4 8 in XSPEC) for clarity.