The high-energy emission from HD 93129A near periastron

ABSTRACT

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 Chandra and NuSTAR and the γ-ray satellite AGILE. High-energy emission coincident with HD 93129A was detected in the X-ray band up to ∼18 keV, whereas in the γ-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_NT,e < 0.02. In addition, we set a lower limit for the magnetic field in the wind-collision region as B_WCR > 0.3 G. We also argue a putative interpretation of the emission from which we estimate f_NT,e ≈ 0.006 and B_WCR ≈ 0.5 G. We conclude that multiwavelength, dedicated observing campaigns during carefully selected epochs are a powerful tool for characterizing the relativistic particle content and magnetic field intensity in colliding wind binaries.

acceleration of particles, radiation mechanisms: non-thermal, gamma-rays: stars, X-rays: stars

1 INTRODUCTION

The role of massive stars as accelerators of Galactic cosmic rays has been gaining more interest in the recent years (e.g. Seo, Kang & Ryu 2018; Aharonian, Yang & de Oña Wilhelmi 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 spectrum, typically a power law with additional features such as a spectral break or an exponential cut-off. 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 has 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 wind-collision 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 characterization 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 2020), although this fascinating object is hardly representative for CWBs in general. Another CWB, γ² 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 yr, 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 include 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 |$\sim 5{{\ \rm per\ cent}}$| 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 characterization 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 non-thermal 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 broad-band 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 multiwavelength 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 organized as follows. The description of the observations is given in Section 2, and the procedure adopted to reduce and analyse each data set in Section 3. We present the results from this analysis in Section 4. In Section 5, we discuss the constraints of our results in the context of theoretical models, and we present a summary and concluding remarks in Section 6.

2 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.

Table 1.

Summary of the X-ray observations analysed.

Instrument	Obs. ID	Date (MJD)	Exposure time (ks)	Effective time (ks)
Chandra	4495	53269.7	58.06	57.33
Chandra	20152	58288.4	25.06	23.84
Chandra	20153	58343.1	25.06	23.84
NuSTAR	30402001002	58290.76	27.59	27.40
NuSTAR	30402001004	58350.03	33.07	32.20

2.1 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. 2002). 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 view of HD 93129A is shown in Fig. 1.

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.

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 arcsec) to estimate contamination by nearby sources, weight was set to yes. We found little variation between both observations from 2018 (less than 2 per cent 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).

2.2 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 arcsec (half power diameter of 58 arcsec; 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 time-scale. 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 |$\lt 3{{\ \rm per\ cent}}$| data loss.¹ 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.

2.3 AGILE observations

We analysed the data collected between 2018 January 1 and 2018 December 31 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 and 400 MeV is 4.2° and 1.2° (68 per cent 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⁻⁶ photons cm⁻² s⁻¹ (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 analysed 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 multisource 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 maximum-likelihood 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 per cent 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.

3 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, Castor & Rybicki 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² 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), normalization (N), and abundance (A). This model is adequate to represent the X-ray emission from a hot plasma such as those produced by hot stellar winds of massive stars (Feldmeier, Puls & Pauldrach 1997; Owocki et al. 2013) and CWBs (Stevens, Blondin & Pollock 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 high-temperature 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 characterize the spectrum of the source (see Section 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 minimization 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.

4 RESULTS

4.1 Chandra

We aim to characterize 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 arcsec on axis). Several sources are detected within 10 arcsec from HD 93129A. We therefore extract the spectrum from HD 93129A using a 1 arcsec circular region so that the background/foreground contamination is minimal. The PSF is such that it ensures that almost 100 per cent of the photons with E < 1 keV, and |$\sim 85{{\ \rm per\ cent}}$| for those with E ∼ 10 keV, are captured (Weisskopf et al. 2002).

We compare the spectra from the 2004 and 2018 observations in Fig. 2. The source was detected up to 8 keV. We fitted both spectra simultaneously with a model [phabs*(vapec₁ + vapec₂)]. 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 minimize 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₁, was also tied between epochs. The values of N_H, T₂, and normalizations norm₁ and norm₂, were fitted independently for each epoch. The best-fitting values are shown in Table 2. We found that between 2004 and 2018 the observed flux increased by more than 50 per cent in the 0.5–3 keV range and 100 per cent in the 3–8 keV range.

Figure 2.

Spectra of HD 93129A as observed with Chandra in 2018 (black) and 2004 (red) extracted from a 1 arcsec circular region. The best-fitting model ([phabs*(vapec₁ + vapec₂)]) is also shown, with the parameters detailed in Table 2.

Table 2.

Results of the fitting of Chandra spectra of HD 93129A using a [phabs*(vapec + vapec)] thermal model. C-statistics were used. The errors at 1σ level are specified for all parameters and for the observed flux.

	1 arcsec extraction region		30 arcsec extraction region
Parameter	2004	2018	2004	2018
NH (1022 cm−2)	|$0.69^{+0.06}_{-0.05}$|	|$0.73^{+0.06}_{-0.05}$|	|$0.64^{+0.04}_{-0.03}$|	|$0.60^{+0.04}_{-0.03}$|
A	|$0.17^{+0.04}_{-0.03}$|		|$0.13^{+0.02}_{-0.02}$|
Si	|$0.48^{+0.09}_{-0.08}$|	|$0.52^{+0.09}_{-0.07}$|	|$0.26^{+0.05}_{-0.04}$|	|$0.39^{+0.11}_{-0.10}$|
S	|$1.13^{+0.26}_{-0.21}$|	|$0.69^{+0.17}_{-0.15}$|	|$0.58^{+0.14}_{-0.13}$|	|$0.36^{+0.11}_{-0.10}$|
Fe	|$0.07^{+0.02}_{-0.02}$|		|$0.06^{+0.01}_{-0.01}$|
|$k\, T_1$| (keV)	|$0.35^{+0.04}_{-0.03}$|		|$0.35^{+0.01}_{-0.02}$|
norm1 (10−2 cm−5)	|$1.30^{+0.52}_{-0.30}$|	|$2.56^{+1.02}_{-0.58}$|	|$2.86^{+0.64}_{-0.49}$|	|$3.32^{+0.75}_{-0.55}$|
|$k\, T_2$| (keV)	|$1.80^{+0.13}_{-0.10}$|	|$2.35^{+0.16}_{-0.14}$|	|$2.41^{+0.14}_{-0.12}$|	|$2.68^{+0.05}_{-0.04}$|
norm2 (10−3 cm−5)	|$1.66^{+0.19}_{-0.18}$|	|$2.22^{+0.20}_{-0.19}$|	|$3.48^{+0.23}_{-0.24}$|	|$3.98^{+0.22}_{-0.19}$|
F0.5–3 (10−13 erg cm−2 s−1)	|$8.39^{+0.01}_{-0.46}$|	|$13.0^{+0.1}_{-0.8}$|	|$16.9^{+0.1}_{-0.5}$|	|$21.3^{+0.0}_{-0.6}$|
F3–8 (10−13 erg cm−2 s−1)	|$1.85^{+0.12}_{-0.14}$|	|$3.79^{+0.17}_{-0.26}$|	|$5.88^{+0.18}_{-0.26}$|	|$7.82^{+0.13}_{-0.31}$|
Cstat/d.o.f. (⁠|$\chi ^2_\mathrm{red}$|⁠)	140.7/110 (1.40)		174.3/124 (1.40)

In order to compare the Chandra and NuSTAR spectra, we need to use a 30 arcsec extraction region for both. For this reason, we compared how selecting a 30 arcsec 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 arcsec. 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 arcsec region with respect to the one from a 1 arcsec region.

In addition, we also checked the 2004–2018 variability considering a 30 arcsec extraction region. The spectral fitting parameters are summarized in Table 2. In this case, the flux between the two epochs increased only in |$26{{\ \rm per\ cent}}$| in the 0.5–3 keV band and |$33{{\ \rm per\ cent}}$| in the 3–8 keV band; similar increase factors were obtained for a 10 arcsec 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 arcsec extraction regions. We obtain that the flux variations for the nearby (background) sources between 2004 and 2018 are |$\lt 10{{\ \rm per\ cent}}$| in both of these energy bands (this can be appreciated more clearly in Fig. A3). Therefore, the observed flux variability within the 30 arcsec 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).

4.2 NuSTAR

The 3–18 keV fluxes between both FPMA and FPMB cameras in each observation were consistent within less than 5 per cent. 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{{\ \rm per\ cent}}$| depending on the selected background. Moreover, the fluxes between the two observations differ only by |$12{{\ \rm per\ cent}}$|⁠, which is almost within 1σ level and can be attributed to calibration uncertainties (up to 10 per cent; 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.

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 underpredicts the emission above 8 keV, which can be explained by the presence of either a hotter component or a non-thermal component.

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 arcsec 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.

In Fig. 4, we show a fit to the Chandra and NuSTAR spectra using an [apec + power-law] model. This fit yields kT = 2.1 ± 0.2 keV, A = 0.21 and |$\Gamma \approx 2.0^{+0.5}_{-0.7}$| (C = 92.6, with 87 d.o.f., and χ²_red = 1.11). We use cflux to obtain the unabsorbed flux from the model and obtain |$F_{3-8} = 7.84^{+0.35}_{-2.61} \times 10^{-13}$| erg s⁻¹ cm⁻² and |$F_{8-18} = 2.35^{+0.20}_{-0.21} \times 10^{-13}$| erg s⁻¹ cm⁻². We also calculate the individual unabsorbed flux of the power-law component and obtain |$F_{8-18} = 1.93^{+0.26}_{-0.36} \times 10^{-13}$| erg s⁻¹ cm⁻².

Figure 4.

Chandra (red) and NuSTAR (black) spectra for 2018. The fitted model is a [apec + power-law], with kT ≈ 2.1 keV and Γ ≈ 1.9 (see text for details).

We also considered the possibility of having only thermal emission. We fitted the Chandra + NuSTAR data with an [apec + apec] model, which yields A = 0.24 (tied for both components), |$kT_1 = 1.7^{+0.4}_{-0.2}$| keV and |$kT_2 = 8.6^{+13}_{-1.5}$| keV (C = 89.8, with 80 d.o.f., and χ²_red = 1.13). The unabsorbed fluxes in this case are |$F_{3-8} = 7.7^{+0.26}_{-0.24} \times 10^{-13}$| erg s⁻¹ cm⁻²and |$F_{8-18} = 2.30^{+0.18}_{-0.18} \times 10^{-13}$| erg s⁻¹ cm⁻².

4.3 AGILE

The result from the analysis of the AGILE data is consistent with an upper limit to the flux in 50 MeV to 100 GeV of 2 × 10⁻⁷ cm⁻² s⁻¹. 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 Section 5.

5 DISCUSSION

In Section 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.

5.1 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 arcsec extraction region is considered. For larger extraction regions with size ≥10 arcsec, 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 arcsec 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. None the less, 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 per cent 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.

• 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⁻¹³ erg s⁻¹ cm⁻² (considering the 1σ upper limit).

• 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⁻¹³ erg s⁻¹ cm⁻² (again assuming the 1σ upper limit).

• The analysis in Section 4.1 suggests that roughly |$\sim 50{{\ \rm per\ cent}}$| 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 F_8–18 ≈ 1.1 × 10⁻¹³ erg s⁻¹ cm⁻².

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.17 v_\mathrm{w,8}^2$| keV (Stevens et al. 1992), with v_w,8 the wind velocity in units of 10⁸ cm s⁻¹. 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 Section 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. 2020). Thus, it is not likely that the detected flux of 8 × 10⁻⁷ cm⁻² s⁻¹ is due to η-Car. None the less, 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⁻¹² cm⁻² s⁻¹, which is more than four orders of magnitude below the detected γ-ray flux (more details in Section 5.2). We also considered a possible hardening in the non-thermal 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.

5.2 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 |$\dot{M}_1=10^{-5}$| M_⊙ yr⁻¹ and |$\dot{M}_2 = 6\times 10^{-6}$| M_⊙ yr⁻¹, and wind terminal velocities v_∞,1 = 3200 km s⁻¹ and v_∞,2 = 2800 km s⁻¹. 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 Benaglia et al. (2015). In this case the relation is f_NT,eB² = 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 Section 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,eB² = constant and F_X ∝ f_NT,e, we get the relation |$B \propto F_X^{-1/2}$|⁠. 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.

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 Section 5.1).

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⁻¹³ erg s⁻¹ cm⁻², 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 periastron 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.

Figure 6.

Modelled SED adopting f_NT,e = 0.006, f_NT,p = 10 f_NT,e, and B_WCR ≈ 0.5 G. Solid lines are calculated for an injected particle energy distribution with constant spectral index of p = 3.2, whereas dash–dotted lines consider a possible hardening of p = 2 in the high-energy spectrum. Dotted lines show the intrinsic emission components before correcting for absorption effects (see Appendix  B for details of the model). We also show the NuSTAR data in the 8–18 keV energy range, the AGILE upper limit assuming a photon index Γ = 2, and the expected sensitivity for 4-yr Fermi-LAT and for 100-h CTA (extracted from Funk, Hinton & CTA Consortium 2013).

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⁻¹. Following del Palacio et al. (2016, and references therein), we can express |$B_\star \approx 2.5 B_\mathrm{WCR} \, (r/R_\star)$|⁠. 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, Owocki & Townsend 2009).

We clarify that magnetic field amplification processes in the WCR (e.g. Bell 2004; Drury & Downes 2012; del Valle, Lazarian & Santos-Lima 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.

6 CONCLUSIONS

We report the results from an observational campaign on the extreme CWB 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 X-ray 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 to 100 GeV range, F < 3 × 10⁻⁷ cm⁻² s⁻¹, 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 multizone 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 wind-collision region as B_WCR > 0.3 G. Neglecting possible magnetic-field amplification in the wind-collision region, we derive a lower limit for the surface stellar magnetic field of B_⋆ > 130 G.

• We consider that |$\sim 50{{\ \rm per\ cent}}$| 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 multiwavelength, dedicated observing campaigns during carefully selected epochs is a powerful tool for characterizing the relativistic particle content and magnetic field intensity in CWDs. 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 angular-resolution observations in the γ-ray band in order to better characterize the non-thermal particle population in colliding-wind binaries.

ACKNOWLEDGEMENTS

We thank the anonymous referee for his/her comments that helped to improve the manuscript. SdP thanks the Committee on Space Research (COSPAR) for support via its Capacity Building Program and its Capacity Building Fellowship Program. SdP, DA, and FG acknowledge support from the Royal Society (RS) International Exchanges ‘The first step for High-Energy Astrophysics relations between Argentina and UK’. FG acknowledges support from Athena project number 184.034.002, partly financed by the Dutch Research Council (NWO). DA acknowledges support from the RS. JMP, VBR, and GER acknowledge support by the Spanish Ministerio de Economía, Industria y Competitividad (MINEICO/FEDER, UE) under grants AYA2016-76012-C3-1- P, MDM-2014-0369 of ICCUB (Unidad de Excelencia ‘María de Maeztu’) and the Catalan DEC grant 2017 SGR 643. GER was also supported by PIP 0338 from the Argentinian Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET). JMA acknowledges support from the Spanish Government Ministerio de Ciencia, Innovación y Universidades through grant PGC2018-095 049-B-C22. RHB acknowledges DIDULS Regular Project No. 18143. KH is supported by the Chandra grant GO8-19010A.

Footnotes

1

http://www.srl.caltech.edu/NuSTAR_Public/NuSTAROperationSite/SAA_Filtering/SAA_Filter.php.

2

http://www.atomdb.org/

3

In del Palacio et al. (2016) this angle was referred as i, which can be confused with the orbital inclination.

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.

Figure A1.

NuSTAR image with FPMB in the 3–11 keV energy range. The source and background extraction regions have a 30 arcsec and 2 arcmin radius, respectively. Left is the first observation, right is the second one.

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 arcsec. For the 1 arcsec source extraction region the background is negligible, whereas for radii ≥10 arcsec 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 arcsec and an inner radius of 1 arcsec 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 |$\lt 5{{\ \rm per\ cent}}$| in the 0.5–3 keV range and |$\lt 10{{\ \rm per\ cent}}$| in the 3–8 keV energy range.

Figure A2.

Comparison of the Chandra spectra extracted from a 1 arcsec (black), 10 arcsec (red), and 30 arcsec (blue) region centred at the position of HD 93129A. The crosses represent the source emission and the circles (with error bars) the background. The spectra are rebinned (setpl rebin 4 8 in xspec) for clarity.

Figure A3.

Comparison of the Chandra spectra of the sources within an annulus of outer radius 30 arcsec and inner radius 1 arcmin centred at the position of HD 93129A; black is 2004 and red 2018. The crosses represent the source emission and the circles the background.

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.

Figure A4.

AGILE light curve in the 50 MeV to 100 GeV energy range during 2018. The top panel does not consider emission from η-Car, while the bottom panel considers it as an extra background source.

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 normalization 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,pf_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⁻³, but with a large uncertainty.

The emission in the radio band is produced by synchrotron emission. This radiation can be significantly attenuated by free–free absorption in the ionized 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 = Dcos ψ is no longer a free parameter of the model.³ A value of ψ ∼ 32° is used in accordance with the most recent orbital ephemeris.

• We considered an increased free–free opacity at radio-frequencies 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, Aharonian & Kelner (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 (equation 16 from del Palacio et al. 2018).