© The Authors 2024

1 Introduction

Over the past two decades, observational studies of the multiplicity of massive O-type stars have led to the conclusion that binarity and multiplicity are far more widespread among these stars than previously thought. As a result, binary interactions could play a predominant role in the evolution of massive stars (e.g. Sana et al. 2012). An interesting object in this context is HD 149404. This system consists of an O7.5 I(f) primary and an ON9.7 I secondary on circular orbits with a period of 9.81 days (Rauw et al. 2001). Raucq et al. (2016) used spectral disentangling to reconstruct the spectra of the individual binary components and subsequently analysed these spectra with the CMFGEN model atmosphere code (Hillier & Miller 1999). The authors found that the ON9.7 I secondary star exhibits a nitrogen overabundance with [N/C] ~ 100 [N/C]_⊙ that can only be explained if the system has gone through an episode of Case A Roche lobe overflow, confirming previous suggestions by Penny et al. (1999) and Thaller et al. (2001). The secondary star of HD 149404 thus displays an evolutionary status intermediate between an O-type star and a WN Wolf-Rayet star.

Rauw et al. (2019) analysed BRIght Target Explorer (BRITE) and Solar Mass Ejection Imager (SMEI) photometry of the binary, finding ellipsoidal variations with a peak-to-peak amplitude near 0.04 mag. Combining the results of these light curves with the Gaia-DR2 parallax allowed us to constrain the orbital inclination of HD 149404 in the range 23°–31° and the Roche lobe filling factor of the secondary star to ≥0.96.

The optical spectrum of HD 149404 displays strong line profile variability in the Hα emission line and in the emission part of the He II λ 4686 line (Rauw et al. 2001; Thaller et al. 2001; Nazé et al. 2002). These variations were interpreted as the signature of a wind-wind interaction occurring between the stars. Evidence of such an interaction was also provided by the BRITE photometry that suggests the presence of a bright spot at the surface of the secondary (Rauw et al. 2019). Such a feature could result from the heating of the secondary photosphere by the wind interaction or from the crash of the primary wind onto the secondary’s photosphere.

In some massive binaries, the shock-heated gas in the windwind collision zone produces significant X-ray emission that manifests itself through an overluminosity and orbital modulation of the observed X-ray emission (Stevens et al. 1992; Pittard & Parkin 2010; Rauw & Nazé 2016). However, the current generation of X-ray satellites also revealed that strong X-ray emissions from the wind interaction zone are the exception, rather than the rule. Most massive binaries actually do not display any significant excess emission above the average canonical level (e.g. Nazé 2009).

From a theoretical viewpoint, HD 149404 is a promising target when it comes to using the X-ray emission from the stellar wind interactions to probe the properties of the winds of the two binary components. HD 149404 is among the 12 X-ray brightest O + O binary systems detected during the ROSAT All-Sky Survey (RASS; Berghöfer et al. 1996) and it should be one of the most straightforward cases to be studied that way. Actually, many of the X-ray brighter O-star binaries in the RASS are higher multiplicity systems (e.g. δ Ori, Oplištilová et al. 2023; τ CMa, Maíz Apellániz & Barbá 2020; HD 150136, Sana et al. 2013), which complicates the interpretation of their X-ray emission. Several other X-ray bright O-star systems have wide eccentric orbits (e.g. ι Ori, Marchenko et al. 2000; 9 Sgr, Rauw et al. 2016a) and display X-ray emission that is notably modulated by the changing orbital separation. On the contrary, for a binary system with a circular orbit, such as HD 149404, the situation should be simpler as putative orbital variations in the observed X-ray emission are expected to arise mostly from occultation effects and from the changing line-of-sight column density towards the shock region as the stars revolve around each other.

To study the X-ray properties of HD 149404, we collected three XMM-Newton observations at key orbital phases to probe the properties of the stellar winds of this system and their interaction. In Sect. 2 we present these new observations and their processing. The X-ray spectra are analysed in Sect. 3, and their variations with orbital phase are modelled under different assumptions in Sec. 4. Our conclusions are presented in Sect. 5.

2 Observations and data processing

HD 149404 was observed with the XMM-Newton satellite (Jansen et al. 2001) on three occasions in August 2018 (see Table 1). The times of the observations were chosen to cover the two conjunction phases and one quadrature phase of the binary. XMM-Newton features three mirror modules that focus X-rays onto three European Photon Imaging Camera (EPIC; Turner et al. 2001; Strüder et al. 2001) instruments and two Reflection Grating Spectrometer (RGS; den Herder et al. 2001) devices (RGS1 and RGS2). The EPIC instruments consist of two Metal Oxide Semi-conductor CCDs (MOS; Turner et al. 2001) and one p-n junction CCD (pn; Strüder et al. 2001).

Our observations of HD 149404 were taken with the EPIC cameras operated in full frame mode. Owing to the optical brightness of the source, we used the thick filters to prevent contamination of the X-ray signal by optical and UV photons. The first 10% of the second observation was affected by a soft proton background flare and was thus discarded. We processed the data with the Science Analysis System (SAS) software version 18.0.0 using the current calibration files available in February 2020. Whilst there are several weak secondary sources in the EPIC field of view, HD 149404 is by far the brightest object in the field and is well isolated from neighbouring X-ray sources.

The EPIC spectra of the source were extracted from a circular source region centred on the SIMBAD coordinates of our target and adopting a radius of 22″.5. The background was extracted from an annular region surrounding the source and with an outer radius of 30″. We also extracted background-corrected light curves of HD 149404 over the 0.5–4.0 keV band for each of the three EPIC instruments and for each of the three observations with temporal bin sizes of 100 and 1000 s. We performed χ² tests to assess the significance of short-term variability adopting different null hypotheses: a constant count rate, a linear trend with time, or a quadratic trend with time. In the majority of the light curves we found no significant variations. During our first observation, the EPIC-pn data showed significant variability consistent with a linear increase of few per cent over the duration of the exposure. However, no significant variation was found in the MOS1 and MOS2 light curves recorded during this same observation. For the second observation, the 100 s bin MOS1 and MOS2 light curves suggested marginally significant variability, which was not confirmed by the 1000 s bin light curves of these instruments nor by the EPIC-pn light curves. Hence, we conclude that in most cases no significant intra-pointing variations were observed beyond the fluctuations consistent with Poisson noise. In addition to the EPIC spectra, we extracted first-and second-order grating spectra of HD 149404 for both RGS instruments.

3 Spectral fitting

The X-ray spectrum of HD 149404 displays emission lines that are characteristic of a thermal plasma (see Fig. 1). Several lines are clearly visible in the combined RGS data. These are Ne X Ly α at 12.13 Å, the He-like fir triplet of Ne IX at 13.70, 13.55, and 13.45 Å, the Fe XVII lines at 15.01 and 17.05 Å, and the O VIII Ly α line at 18.97 Å. Whilst the quality of the RGS data is not sufficient to perform a detailed analysis of the line profiles for individual observations, the combined RGS spectra indicate that the f component of the Ne IX triplet is suppressed, which is expected in the presence of a strong UV radiation field (Porquet et al. 2001) such as that in the circumstellar environment of HD 149404. No line is seen at longer wavelengths. This reflects a strong absorption by circumstellar and interstellar material, which can also be seen by the drop in the EPIC-pn spectra at energies below ~0.7 keV. Above photon energies of 4 keV the EPIC spectra reveal no significant emission, suggesting that the X-ray emitting plasma is rather soft (see the bottom panel of Fig. 1).

To characterise the overall spectral energy distribution of HD 149404’s X-ray spectrum, we used version 12.9.0i of the xspec code (Arnaud 1996). We attempted to simultaneously fit the EPIC (i.e. MOS1, MOS2, and pn) and RGS (both orders of RGS1 and RGS2) spectra, using various combinations of absorbed optically thin thermal plasma models of the kind $$\mathtt{\text{ TBabs * phabs * }}\sum _{i=1}^n\mathtt{\text{vapec}}\left(k{T}_i\right).$$(1)

In expression (1) the interstellar absorption is accounted for by using the TBabs model (Wilms et al. 2000). Based on the results of Diplas & Savage (1994) and Rachford et al. (2009), Gudennavar et al. (2012) quoted a total interstellar neutral hydrogen column density of (3.72 ± 0.94) × 10²¹ cm⁻². This value is in excellent agreement with the total neutral column density estimated from the reddening (3.94 × 10²¹ cm⁻² for E(B − V) = 0.68) using the N(H I)/E(B − V) relation of Bohlin et al. (1978), and agrees reasonably well with the value (4.16 × 10²¹ cm⁻²) inferred from the N(H I)/E(B − V) relation of Gudennavar et al. (2012). We thus set the interstellar neutral hydrogen column density to 3.72 × 10²¹ cm⁻². Additional photoelectric absorption by the stellar winds is included by means of the phabs model component, where the hydrogen column density is considered a free parameter of the model. The emission by the hot plasma is modelled by the sum of up to three collisionally ionised, optically thin vapec thermal plasma components (Smith & Brickhouse 2001). Such models are known to provide a good overall description of the X-ray spectral energy distribution of O-type stars (e.g. Nazé 2009). We explicitly requested the different vapec plasma components to have identical chemical compositions.

Whilst the X-ray spectra of massive stars are frequently described by expressions such as (1), this is certainly not a unique approach. For instance, one could also associate a particular circumstellar column density with each plasma component. However, there obviously exists some model degeneracy between the action of absorbing material and the intrinsic hardness of the various plasma components. Though we try to be as self-consistent as possible, we note that the primary goal of our xspec fitting procedure is to obtain a good description of the overall X-ray spectral energy distribution. The specific parameters of the various model components should thus not be overinterpreted.

Our very first fitting trials indicated that a single plasma component is not sufficient to achieve a decent quality fit of the spectra. Good fitting qualities were instead achieved for models with two or three components. We set the chemical abundances to solar according to Asplund et al. (2009), except for the C, N, and O abundances. Tests with these three abundances allowed to vary indicated that the spectral fits are essentially insensitive to the C abundance, whereas they systematically yielded a significant N overabundance and an O abundance consistent with solar. These results most likely indicate that the X-ray emitting plasma contains some gas from the secondary star. We thus tested several assumptions for the CNO abundances of the X-ray emitting plasma: solar composition, a pure primary wind composition, a 50–50 mix of primary and secondary wind compositions, and a pure secondary wind composition, the last three with abundances as inferred by Raucq et al. (2016).

The quality of the fits significantly improved when allowing for non-solar composition. The best fit was achieved for an equal proportion mix of primary and secondary abundances (χ_ν between 1.05 and 1.15; see Tables 2 and 3). However, the fit quality significantly degraded $(\mathrm{\Delta }{\chi }_v^2\sim +0.20)$ for a 100% secondary composition of the emitting plasma. The reason for this is the low O abundance of the secondary wind inferred by Raucq et al. (2016), which seems inconsistent with the EPIC spectra and with the presence of the O VIII Ly α line in the RGS spectra. A less extreme degradation of the fit quality ($(\mathrm{\Delta }{\chi }_v^2\sim +0.05)$) occurs if we assume a pure primary wind composition. These tests hence suggest that the X-ray emitting plasma in HD 149404 contains a mix of material from both winds. Finally, we also tested models where the proportion of primary and secondary abundances was left as a free parameter. Those models yielded a secondary abundance contribution of (46 ± 12)%. The quality of these fits was essentially the same $(\left|\mathrm{\Delta }{\chi }_{\nu }^2\right|\le 0.01)$ as for the 50–50 mix models. In the remainder of this paper we thus focus on the results obtained for the equal proportion mix.

For a two-temperature (2-T) plasma consisting of a 50–50 mix, we fitted the spectra of the three observations independently keeping the wind column density and the plasma temperatures as free parameters. The temperatures of the two plasma components were found to be very similar in all three observations: kT₁ = 0.20 ± 0.04 keV and kT₂ = 0.65 ± 0.06 keV. We thus repeated the fits, fixing the two plasma temperatures, but allowing their emission measures to vary from one observation to another (see Table 2).

We repeated the same procedure for a three-temperature (3-T) plasma model with chemical abundances consisting of a 50–50 mix. Again, the temperatures of the three plasma components were found to be very similar in all three observations: kT₁ = 0.14 ± 0.04 keV, kT₂ = 0.35 ± 0.07 keV, and kT₃ = 0.73 ± 0.02 keV. The results obtained with these temperatures fixed but allowing the emission measures and the wind column density to vary are quoted in Table 3. An F-test to assess whether including a third vapec component significantly improves the quality of the fit yields false alarm probabilities of 10⁻⁶, 10⁻⁸, and 6 × 10⁻⁷ for the first, second, and third observations, respectively. Hence, we conclude that the third plasma component indeed provides an improvement in the fit quality.

This soft X-ray emission might be surprising at first sight for a putative interacting wind system, but it is not unusual for short period O-star binaries. For instance, XMM-Newton spectra of the O7.5 III-II(f) + O7.5 III-II(f) system HD 152248 (P_orb = 5.8 d) were fitted by Sana et al. (2004) with a dominant plasma component with kT ~ 0.2 keV and a hottest plasma component with kT ~ 0.8 keV. Likewise, for the ASCA-SIS spectra of the highly eccentric O9 III + B1 III system ι Ori (P_orb = 29.1 d, e = 0.76), Pittard et al. (2000) inferred best-fit plasma temperatures of 0.15 keV and 0.61 keV.

Over the 0.5–4.0 keV energy domain, the fluxes vary between 4.9 × 10⁻¹³ erg cm⁻² s⁻¹ and 7.2 × 10⁻¹³ erg cm⁻² s⁻¹. When correcting for the absorption by the interstellar material only, the corresponding fluxes are 1.09 × 10⁻¹² erg cm⁻² s⁻¹ and 1.53 × 10⁻¹² erg cm⁻² s⁻¹. Comparing these numbers to the bolometric fluxes of the system (Raucq et al. 2016; Rauw et al. 2019), we obtain a distance-independent estimate of log $\frac{f_X}{f_{\text{bol }}}=-7.13\pm 0.03$ at minimum and −6.97 ± 0.03 at maximum. This result is very similar to the canonical value of O-type stars (Nazé 2009), indicating that HD 149404 is one of those O-star binaries that do not exhibit any strong X-ray overluminosity. If anything, the system would appear to be underluminous. We come back to this point in Sec. 4.3.

Together with the overall softness of the X-ray emission, these results suggest that HD 149404 does not exhibit any strong X-ray excess emission arising from the shocked winds. This situation probably stems from the fact that the stars that make up this binary system are large compared to their orbital separation. As a result, the stellar winds do not have enough room to accelerate to high velocities before they interact. Therefore, the material in the wind interaction zone may be too cool to produce any strong and hard X-ray emission (see Sec. 4.3).

The variations in the X-ray spectrum and fluxes of HD 149404 as a function of orbital phase are illustrated in Figs. 1 and 2. In Fig. 2 we also include the results of an archival ROSAT-PSPC observation, taken at phase φ = 0.70, which was previously discussed in Rauw et al. (2001). This figure shows that the X-ray fluxes display rather limited variability, except for the drop around phase 0.0, which corresponds to the conjunction with the ON9.7 I secondary being in front. At this phase, the line of sight towards the wind interaction zone mostly samples the wind of the ON9.7 I secondary, whilst it probes the wind of the O7.5 I(f) primary at the other conjunction phase 0.5.

It is worth noting that the time of minimum flux, around phase 0.0, corresponds to the lowest circumstellar column density in our xspec fitting. This behaviour is also observed for the 2-T and 3-T plasma models. We note that the model parameters should not be overinterpreted, and that the formal errors on N_H do not allow us to rule out the possibility of a constant column density. However, our results suggest that the observed flux variation is not due to phase-modulated absorption column density. Instead, the fits in Tables 2 and 3 indicate that the lower flux at phase 0.0 stems from a global reduction of the emission measure of the emitting plasma at that phase (see also Fig. 3). This suggests that occultation of the emitting region plays an important role.

Fig. 1 X-ray spectra of HD 149404. Top panel: combined RGS (RGS1 + RGS2, first and second orders, all three observations) spectrum of HD 149404. The labels indicate the lines that are apparent in this spectrum. The error bar in the top left corner indicates the mean error on the energy bins in the combined RGS spectrum between 10 Å and 20 Å. Bottom panel: comparison of the EPIC-pn spectra of HD 149404 during the three observations. The first observation (φ = 0.227) is shown in blue, the second (φ = 0.992) in red, and the third (φ = 0.487) in green.

Fig. 2 Variations of the X-ray emission of HD 149404 with orbital phase. From top to bottom, the panels illustrate the full EPIC-pn count rates for the different observations and the observed fluxes evaluated over three different energy bands: soft (0.5–1.0 keV), medium (1.0–2.0 keV), and hard (2.0–4.0 keV). The magenta symbols in the second and third panels correspond to the fluxes measured on an archival ROSAT-PSPC spectrum (Rauw et al. 2001).

Fig. 3 Ratio of the observed fluxes evaluated over narrow energy bands. The blue (resp. red) symbols indicate the ratios of the fluxes measured at phase 0.227 (resp. 0.992) to those determined at phase 0.487, the phase at which the emission is strongest.

4 The variations in the X-ray emission

Figure 3 illustrates the ratio of the observed fluxes in various energy bands at the orbital phases sampled by our XMM-Newton observations. A remarkable property of these variations is that they are essentially grey (i.e. nearly independent of photon energy). In this section, we attempt using the observed variations in the X-ray spectrum between different orbital phases to probe the properties of the winds of HD 149404. For this purpose, we consider various scenarios: either the observed X-ray emission is dominated by the wind interaction zone, or it is mostly intrinsic to the stellar winds with only a moderate contribution from the wind interaction zone.

4.1 X-ray opacity of the stellar winds

To estimate the X-ray opacities κ_E of the stellar winds of HD 149404’s components, we used the CMFGEN non-LTE model atmosphere code (Hillier & Miller 1999) together with the stellar parameters and chemical compositions inferred by Raucq et al. (2016) and recalled in Table 4. Overall, the computed opacities of the two winds are rather similar, except for some differences in the region around 0.5 keV, which are mainly due to the enhanced nitrogen and depleted carbon abundances of the secondary star. Over the energy domain between 0.5 and 4.0 keV, which is of interest for our present study, the opacities vary relatively little with position inside the wind. It is only for the innermost positions in the wind that the opacity drops significantly owing to the higher ionisation stage of the wind material in these regions. We thus consider that the opacity at a distance from the stellar centres r = a/2, where a is the semi-major axis of the orbit, is representative of the wind opacity over those regions that are most relevant to our present purpose.

The optical depth along a given line of sight is obtained by integrating the product of κ_E by the wind density ρ $${\tau }_E=\int {\kappa }_E\rho \mathrm{d}s,$$(2)

where the integration starts from the position of the considered point in the wind interaction zone and proceeds along the line of sight towards the observer. To compare the effects of the two winds, Fig. 4 shows the product of κ_E by ρ evaluated at a position r = a/2 in the wind. As is apparent, the CMFGEN models predict the primary wind to be slightly more opaque than the secondary wind at most energies. This would be in line with the marginally higher column density in the xspec models around phase 0.5 compared to phase 0.0. However, our observational results (Figs. 2 and 3) indicate that the observed X-ray fluxes are lower at the conjunction with the secondary star being in front (φ = 0.0), thus suggesting once more that the optical depth variations are not the main cause of the reduced observed flux at phase 0.0.

Whilst the values of κ_E ρ of both stars are remarkably similar (Fig. 4), they vary very strongly with energy: the values at 0.5 and 5.0 keV differ by about a factor of 100. This provides another argument suggesting that the nearly grey attenuation of the observed fluxes seen in Fig. 3 cannot stem from the sole effect of changing wind absorption, at least not if the emission at all energies arises from the same geometrical region. A possibility to produce a grey attenuation would be through a phase-dependent occultation effect of the emitting region by the stellar bodies. Adopting our best estimates of the stellar radii (Table 4), the surface area of the secondary star is about twice as large as that of the primary. If the X-ray emitting region is located very close to the secondary’s surface facing the primary, we would expect a smooth phase-dependent modulation of the observable emission with the minimum flux occurring at phase 0.0. Because of the low orbital inclination, the more extended the emission region, the smaller the expected variations in the X-ray flux between different orbital phases.

Fig. 4 Product of wind X-ray opacity κE times wind density ρ as a function of energy. The opacity and wind density are computed from the CMFGEN models of Raucq et al. (2016) at a distance r = a/2 from the stellar centre. The continuous red line (labelled 1) shows the results for the primary star, whereas the dashed blue line (labelled 2) corresponds to the secondary star of HD 149404.

4.2 X-ray emission from the wind interaction zone

Although we show in Sect. 4.1 that the X-ray luminosity of HD 149404 is rather modest and that the spectrum is soft, we start by considering as a first hypothesis that the observed X-ray emission of HD 149404 is dominated by the wind interaction zone. For this part of the analysis, we followed the same approach as previously done by Rauw et al. (2014) in the study of HDE 228766. We used the formalism of Cantó et al. (1996) to compute the shape of the axisymmetric wind interaction zone (see Fig. 5) and discretise the three-dimensional interaction region into a total of 200 × 360 cells along the θ and φ spherical coordinates.

As a first step we computed the density and plasma temperature in the wind interaction zone and the associated X-ray emission by the shock-heated plasma in the same way as done in Rauw et al. (2016b). This approach assumes that the plasma in the wind interaction zone cools adiabatically. We then computed the optical depth along the line of sight towards a given cell of material in the wind interaction zone (see Appendix A).

Our calculations further accounted for occultations of the plasma cells by either of the two stars, as well as the lines of sight that cross both winds. For the plasma emission, we used emissivities that are representative of prominent spectral lines that are present over the range of temperatures found in the fits of the XMM-Newton spectra (see Table 3). For instance, the emissivity Q_E(T) in the 0.5–0.75 keV energy band was taken to be that of the strong O VIII Ly α line at 0.65 keV.

The relevant parameters in this model are the mass-loss rates ${\dot{M}}_1$ and ${\dot{M}}_2$ of the stars with the stronger and weaker wind respectively, the associated wind velocities υ_∞,1 and υ_∞2, as well as the orbital inclination i. Taking advantage of the constraints on i established from the BRITE and SMEI light curves (Rauw et al. 2019), we restricted our investigation to inclinations in the range 20°-40°. The geometrical shape and the location of the shock region are defined by the wind-momentum ratio $\eta =\frac{{\dot{M}}_1{v}_1}{{\dot{M}}_2{v}_2}$ (Cantó et al. 1996). The distance from the centre of the star with the strongest wind (assumed here to be star 1) to the shock is given by $$x_{\text{stag }}=\frac{\sqrt{\eta }a}{1+\sqrt{\eta }},$$(3)

whilst the asymptotic opening angle θ_∞ of the shock cone (see Fig. 5) is the solution of the equation (Cantó et al. 1996) $${\theta }_{\mathrm{\infty }}-\tan {\theta }_{\mathrm{\infty }}=\frac{\pi \eta }{\eta -1}.$$(4)

In the most extreme case the dominant wind crashes directly onto the surface of the star with the less energetic wind. In the case of HD 149404, this happens for values of η exceeding 5.44, which correspond to θ_∞, values larger than about 120°. In our models, situations where η exceeds this limiting value are treated as cases where observable X-rays are only emitted by the part of the shock region that is detached from the secondary surface. The emission zone along the wind interaction region was assumed to be restricted to cells located within a distance r_max from the centre of the star with the weaker wind (see Fig. 5).

For a given set of parameters, we can then use this model to predict the ratio of the X-ray fluxes that would be observed at phases 0.0 and 0.5 on the one hand, and phases 0.25 and 0.5 on the other hand. We finally compare these predictions to the observed values by evaluating $${\chi }^2={\displaystyle \sum _{{\varphi }_j}\sum _{E_i}}\frac{{[{\left(\frac{f_X({\varphi }_j,E_i)}{f_x(0.487,E_i)}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}-{\left(\frac{f_X({\varphi }_j,E_i)}{f_x(0.487,E_i)}\right)}_{\mathrm{m}\mathrm{o}\mathrm{d}}]}^2}{{\sigma }^2{\left(\frac{f_x({\varphi }_j,E_i)}{f_x(0.487,E_i)}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}}$$(5)

for the two values of φ_j and the six values of E_i.

We computed extensive grids of models in the five-dimensional parameter space (i, ${\dot{M}}_{\text{Prim }},v_{\text{Prim }},{\dot{M}}_{\text{Seco }},v_{\text{Seco }}$). We tested various assumptions on the size of the emitting region, adopting r_max of a/2, 3 a/4, a, 3 a, 5 a, and no limitation on the size of the emitting zone. We further tested the assumptions that the wind momentum of the primary dominates over that of the secondary or that the reverse scenario holds.

Models where the secondary wind is assumed to dominate are significantly worse in reproducing the observed variations than models where the primary wind dominates. However, even for the latter models, the achieved quality of the best fits (i.e. those having the lowest χ²) is comparatively poor. The best results are obtained for r_max = a/2 and values of η close to 14 (see Fig. 6) with i = (23.0 ± 0.9)°, ${\dot{M}}_{\text{Prim }}$ = (1.9 ± 0.4) × 10⁻⁷ M_⊙ yr⁻¹ ${\dot{M}}_{\text{Seco }}$ = (2.7 ± 0.4) × 10⁻⁷ M_⊙ yr⁻¹, υ_Prim = (2539 ± 184) km s⁻¹, and υ_Seco = 125 km s⁻¹. Figure 6 shows that even these best-fit models fail to reproduce correctly the energy dependence of the flux ratios, especially for the phase φ = 0.227 observation. We also note that the best-fit parameters are quite different from the wind properties inferred by Raucq et al. (2016). The low secondary wind velocity is not fully unexpected since in these models the primary wind crashes onto the secondary surface. What is more surprising is the fact that, in these models, both stars have nearly identical mass-loss rates, with ${\dot{M}}_{\text{Seco }}$ slightly exceeding ${\dot{M}}_{\text{Prim }}$, whereas the opposite is expected (see Table 4). Whilst the secondary mass-loss rate would be close to the value expected from the models of Raucq et al. (2016), the primary mass-loss rate would be five times below the expected value for a star of this spectral type.

Owing to the small orbital separation of the stars, the plasma density in the wind interaction region is likely sufficient for radiative cooling to play a significant role (Stevens et al. 1992). Under these circumstances, the shock-heated plasma cools before being advected out of the wind interaction region. This effect is amplified by the radiative inhibition of the wind velocity by the radiation pressure of the companion star (Stevens & Pollock 1994). Radiative inhibition lowers the pre-shock velocities, thereby leading to cooler and denser post-shock plasma. In this case, the contribution of a given cell of gas to the X-ray luminosity is no longer proportional to the square of the plasma density, but instead scales linearly with the kinetic energy flux flowing into the shock (e.g. Stevens et al. 1992; Kee et al. 2014). Following Kee et al. (2014), the emission integrated over both shock cooling layers is proportional to ${\rho }_0{v}_0^3$, where ρ₀ and υ₀ are respectively the pre-shock wind density and the pre-shock wind velocity perpendicular to the shock. The power radiated in X-rays by a given cell is then evaluated as the product of ${\rho }_0{v}_0^3$ times the area of the cell along the shock surface. However, even this model provides only a first approximation to the effect of radiative cooling. As shown by Kee et al. (2014), the occurrence of thin-shell instabilities (Vishniac 1994) not only leads to a significant reduction of the overall X-ray emission compared to the above scaling, but further yields a rather complex shape of the wind interaction zone. As a result, X-ray emission arises from a limited number of discrete patches associated with the tips and troughs of the unstable shear structure.

Another aspect that might be of relevance here is related to the impact of the orbital motion. In a short-period binary system such as HD 149404, the Coriolis pseudo force could lead to a significant curvature and non-axial geometry of the wind interaction zone. Evidence of such a curvature was seen in the variability of the optical emission lines (Rauw et al. 2001; Nazé et al. 2002). It should be emphasised though that these optical emission lines probably arise from relatively distant parts of the wind interaction zone (Rauw et al. 2001; Nazé et al. 2002), whilst the X-ray emission is expected to form in the inner parts.

To determine how these two effects impact the putative X-ray emission from the wind interaction zone of HD 149404, we used the LIFELINE code (Mossoux & Rauw 2021), which combines the treatment of radiative driving and inhibition of Stevens & Pollock (1994) with the formalism of Parkin & Pittard (2008) to treat the Coriolis deflection. LIFELINE uses the steady-state planar post-shock cooling layer model of Antokhin et al. (2004) to describe the density and temperature profiles of the plasma in the post-shock region that cools radiatively. The Antokhin et al. (2004) model is known to predict X-ray luminosities and plasma temperatures that are too high (e.g. De Becker et al. 2004). This is notably because of the above-described impact of thin-shell instabilities (Vishniac 1994; Kee et al. 2014). Hence, the plasma temperatures and X-ray intensities predicted by LIFELINE should be seen as upper limits only.

We performed LIFELINE calculations for three sets of model parameters (see Table 5). Model I corresponds to the parameters of Table 4. Model II uses half mass-loss rates, whilst Model III considers a larger orbital inclination. In all three models, we observed a severe radiative inhibition effect. For instance, the primary wind reaches pre-shock velocities between 272 km s⁻¹ (Model III) and 465 km s⁻¹ (Model II). In Models I and III, the primary wind crashes onto the secondary star despite the inclusion of radiative inhibition (see Fig. 7). The LIFELINE simulations predict that the hottest regions of the wind interaction zone are located in the zone where the shock detaches from the secondary surface. These regions are essentially unaffected by the Coriolis deflection (see Fig. 7). Therefore, the Coriolis deflection is unlikely to play a significant role in the interpretation of the variations in the X-ray emission of HD 149404.

Fig. 5 Schematic illustration of the wind interaction zone in the orbital plane of HD 149404 modelled with the formalism of Cantó et al. (1996).

Fig. 6 Comparison between the observed flux ratios and the best-fit predicted values for different models tested in this work. The red and blue symbols with error bars correspond respectively to the measured flux ratios displayed in Fig. 3. The orange and cyan curves yield the best-fit theoretical models for the phases shown respectively by the red and blue symbols. The continuous lines assume the X-rays arise entirely from the wind interaction zone of HD 149404 and using temperature and density profiles computed with the Cantó et al. (1996) model for an adiabatic wind interaction zone (Sect. 4.2). The long-dashed lines yield the results for X-rays arising entirely from wind-embedded shocks intrinsic to the primary and secondary stars (Sect. 4.3).

4.3 An alternative description of the X-ray emission

Since HD 149404 does not appear to be a very bright and hard X-ray source, a significant part of its X-ray emission could arise from embedded wind shocks resulting from the intrinsic instability of the wind driving mechanism (Lucy 1982; Owocki et al. 1988; Feldmeier et al. 1997). These shocks are expected to be distributed throughout the winds of the individual binary components. This mechanism is considered to produce the intrinsic X-ray emission of massive stars, the X-ray luminosity of which scales with the bolometric luminosity L_bol (e.g. Nazé 2009; Rauw 2022, and references therein).

Unlike a single massive star, the stellar wind of a given component of a close binary system such as HD 149404 cannot deploy equally in all directions. This is especially true in the case of a collision of the dominant wind with the photosphere of the companion. In such a situation the weaker wind only deploys over a solid angle limited by the wind interaction region. Using the notations of the Cantó et al. (1996) model, the solid angle occupied by the weaker wind amounts to 2 π (1 + cos θ_∞). For θ_∞ = 120°, which is the asymptotic shock angle predicted by the Cantó et al. (1996) formalism for a wind-photosphere collision (see above) and which agrees well with the LIFELINE computations (see Fig. 7), the secondary’s wind is able to deploy freely only over about one-quarter of the full 4 π solid angle. Hence, compared to a single star of same L_bol, the embedded wind shocks of the secondary star would contribute at most one-quarter of the X-ray emission expected from the canonical $\frac{L_X}{L_{\text{bol }}}$ relation. Likewise, the primary wind can fully deploy only over a solid angle 2 π (1 − cos θ_∞), that is three-quarters of the full solid angle. However, if the separation between the two stars is sufficient, that is to say, more than a ≃ 2 R_*, embedded wind shocks might also develop in the part of the primary wind located between the two stars, ahead of the wind interaction zone. As we describe above, the wind velocity in that region is strongly impacted by radiative inhibition due to the secondary’s radiation field (Stevens & Pollock 1994), probably leading to weaker shocks and softer X-ray emission from the wind embedded shocks. To be conservative, we consider only the reduction of the intrinsic emission of the secondary star by a factor of four. This effect diminishes the intrinsic X-ray emission to 65% of the value expected from the total bolometric luminosity of the binary system. Any effect on the primary’s intrinsic emission would further reduce this value. In terms of log f_X this corresponds to a reduction in the expected flux due to wind embedded shocks of about 0.19 dex.

Observational studies of high-resolution X-ray spectra of single O-type stars showed that most of the X-ray emission arises from rather deep inside the stellar winds (typically within the inner 5 R_*), though there are clear differences in the emission formation regions according to the plasma temperature and the wind density (Hervé et al. 2013; Rauw et al. 2015). The outer boundary of an emission region is ruled by the competition between the outward decreasing emission measure and the outward decreasing optical depth of the overlaying cool stellar wind. For instance, in the case of ζ Pup (log $\dot{M}$ = −5.46), Hervé et al. (2013) derived emission regions that range between 1.5 and 38 R_* for the softest plasma component (kT = 0.20 keV), 2.7 and 4.0 R_* for the kT = 0.40 keV component and between 3.1 and 4.1 R_* for the hottest plasma component (kT = 0.69 keV). In the case of λ Cep (log $\dot{M}$ = −5.85), most of the X-ray emission was found to arise from regions between 1.1 and 2.5 R_* (Rauw et al. 2015). In both ζ Pup and λ Cep, the volume filling factors of the hot gas were found to be near 0.01 (Hervé et al. 2013; Rauw et al. 2015). The mass-loss rates of the stars in HD 149404 are smaller than those of either ζ Pup or λ Cep, resulting in lower density winds. Hence, we assume the intrinsic X-ray emission of those stars to form dominantly within about 2 R_* from the stellar centre.

To test the hypothesis that the entire observable X-ray emission arises from wind embedded shocks, we computed grids of models where three plasma components at temperatures found in the 3-T plasma models (see Table 3) are distributed across the wind of the primary and the part of the secondary wind contained within a shock cone attached to the secondary surface and with θ_∞ = 120°. Each plasma component was assumed to extend out to a radial distance of r_out. The emission region was discretised into cells of 0.1 R_* in radial extension, and 2° in both latitude and azimuth angles. For each of these emitting cells, we computed the emission measure and the optical depths across the wind using Eqs. (A.3)–(A.6). We accounted for the fact that a given line of sight can intersect the shock cone up to two times, thus implying that for such lines of sight the winds of both stars contribute to the optical depth. We varied r_out in steps of 0.1 R_* between 1.05 and 2.45 R_*, independently for each of the three plasma temperatures. We further considered volume filling factors of the two hotter components, which were varied between 1/16 and 8 times the volume filling factor of the coolest plasma component. Finally, we tested inclinations ranging from 23 to 31°. We adopted the wind properties given in Table 4.

Although we explored a wide range of parameters, we found no combination of model parameters that successfully reproduced the observed variations (see the long-dashed lines in Fig. 6). Instead, the results were much worse than those obtained assuming the X-ray emission to arise from the sole wind interaction zone. The models that came closest to the observations predict fluxes at all energies to be larger at φ = 0.0 (secondary in front) than at φ = 0.5 which is the reverse of what we observe. This situation stems from the combination of a slightly larger primary wind opacity and occultation effects of the secondary wind by the body of the secondary star at φ = 0.5. Actually, the largest flux is predicted at quadrature phases (0.25 or 0.75), whilst observations indicate this to be the case at φ = 0.5. Moreover, the corresponding best-fit models predict a significant energy-dependence of the flux ratios, at odds with the observations. We also tested a grid of models where the primary mass-loss rate was reduced by a factor two, all other parameters remaining unchanged. These models did not improve the agreement with the observations. We thus conclude that the observed variations in the X-ray flux of HD 149404 cannot be explained if the X-rays arise solely from embedded wind shock emission attenuated by an overlying smooth wind.

Fig. 7 Theoretical temperature (top panel) and density (bottom panel) maps of the wind interaction zone in HD 149404 computed with the LIFELINE code (Mossoux & Rauw 2021), adopting the parameters of Model III (see Table 5). For this set of parameters, the wind of the primary star (on the right) crashes onto the surface of the secondary (on the left).

5 Discussion and conclusions

HD 149404 displays variations in its X-ray emission at the ~30% level. The system appears X-ray brightest when the O7.5 I(f) primary is in front, and X-ray faintest when the ON9.7 I secondary is in front. In a system with a circular orbit, such as HD 149404, this behaviour is at first sight indicative of a phase-dependence of the optical depth of circumstellar column towards the X-ray source. This is reminiscent of the case of the O7 + Of⁺/WN8ha 10.7-day period binary HDE 228766 (Rauw et al. 2014). For this latter system, the X-ray flux is severely absorbed at phases when the Of⁺/WN8ha star is in front compared to phases when it is behind, which was consistently interpreted as the result of the higher opacity and higher density of the Of⁺ /WN8ha wind and the X-ray emission arising mostly from a wind-wind interaction zone. However, there are clear differences between the behaviours of HD 149404 and HDE 228766. First, the flux variations observed in HD 149404 are of significantly lower amplitude than those of HDE 228766. This could be due to the lower inclination of the orbit of HD 149404, found to be in the range 23°-31° (Rauw et al. 2019), whilst it was estimated in the range 54°-61° for HDE 228766 (Rauw et al. 2014), and to the lower mass-loss rate of the primary in HD 149404 compared to the more extreme wind of the Of⁺/WN8ha star. More surprisingly, the amplitude of the changes in the measured fluxes of HD 149404 appears rather independent of photon energy (see Fig. 3), whereas a clear energy-dependence is expected for an absorption effect and was indeed observed for HDE 228766 (Rauw et al. 2014). The grey energy-dependence of the variations observed for HD 149404 therefore suggests that the bulk of the variations instead stem from phase-dependent occultations of the X-ray emitting zone by the stellar bodies.

In our present study, we considered two possible physical origins for the X-ray emission: pure wind-wind collision and pure intrinsic emission arising from embedded wind shocks in the winds of the individual stars. Reproducing the observed variations as a function of energy is a challenge for both models. Because of their expected strong energy dependence, the models involving intrinsic emission perform significantly worse in reproducing the observed behaviour. Taken at face value, this suggests that the X-ray emission mostly arises from the innermost parts of the wind interaction between the stars and with a location very close to the secondary star. However, with log $\frac{L_X}{L_{\text{bol }}}$ ranging between −7.13 and −6.97, the overall level of the X-ray emission of HD 149404 is comparatively low, unlike what one would expect in the presence of an X-ray bright wind-wind interaction.

The most plausible explanation for the system’s grey variations of the X-ray emission could be that the stellar winds consist of partially optically thick fragments. Feldmeier et al. (2003) and Oskinova et al. (2004) investigated a situation, based on the results of the hydrodynamic models of Feldmeier et al. (1997), where the stellar wind consists of randomly distributed fragments of compressed shells with a small lateral extension. The individual fragments could be optically thick, allowing photons to escape only via the lateral gaps between the fragments. Equation (35) in Oskinova et al. (2004) provides an expression for the ensuing effective optical depth of such a porous wind. For optically thick fragments, this effective optical depth is indeed independent of photon energy, but depends on the average number of absorbing fragments along the radial direction < N_r > and the radial extent of the wind region filled with such fragments. Wind porosity would allow soft X-ray photons arising inside the winds to escape, whereas such photons would be almost totally absorbed by a homogeneous wind. In the absence of a self-consistent formalism to predict the radial extent of the fragmented winds and the values of < N_r >, it is difficult to perform a quantitative comparison with the observations. However, since the highest X-ray emission level is observed around conjunction with the primary star in front, we can conclude that the primary wind should be more porous than the secondary wind, despite the higher mass-loss rate of the primary star (Raucq et al. 2016).

In this context, it is also worth recalling that the BRITE photometry of HD 149404 displayed a substantial red noise stochastic variability (Rauw et al. 2019). This could reflect the presence of a sub-surface convection zone (Cantiello et al. 2009) in at least one of the stars of the binary. Such a sub-surface convection zone could in turn trigger clumping, leading to a highly structured inner part of the wind.

The most likely explanation of the phase behaviour of HD 149404’s X-ray emission thus appears to be that the stellar winds are clumpy and that these clumps are at least partially optically thick in the X-ray domain.

Acknowledgements

The Liège team acknowledges support from the Fonds National de la Recherche Scientifique (Belgium) notably under grant no T.0192.19, and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contract HERMeS). We are greatful to Dr. Y. Nazé for constructive discussions. ADS and CDS were used for this research.

Appendix A Optical depth along the line of sight

To evaluate the optical depth from a given plasma cell towards the observer, we express the position of the plasma cell in the cylindrical (p, z) coordinates of the star that is in front. For a wind velocity law $$v=v_{\mathrm{\infty }}(1-\frac{R_{\ast }}{r}),$$(A.1)

and defining $${\tau }_{E,\ast }=\frac{{\kappa }_E(r=a/2)\dot{M}}{4\pi v_{\mathrm{\infty }}{R}_{\ast }},$$(A.2)

the optical depth along the line of sight towards an observer located at z → ∞ then becomes $${\tau }_E(p,z)={\tau }_{E,\ast }\ln \left(\frac{z/R_{\ast }}{z/R_{\ast }-1}\right)\text{ for }p=0,$$(A.3) $${\tau }_E(p,z)=\frac{{\tau }_{E,\ast }}{\sqrt{1-\frac{p^2}{R_{\ast }^2}}}\ln \left(\frac{\alpha -x_2}{\alpha -x_1}\right)\text{ for }0<p<R_{\ast }\text{, }$$(A.4) $${\tau }_E(p,z)=\frac{2{\tau }_{E,\ast }}{\alpha -1}\text{ for }p=R_{\ast },\text{ and }$$(A.5) $${\tau }_E(p,z)=\frac{2{\tau }_{E,\ast }}{\sqrt{\frac{p^2}{R_{\ast }^2}-1}}[\frac{\pi }{2}-\arctan \left(\frac{\frac{\alpha p}{R_{\ast }}-1}{\sqrt{\frac{p^2}{R_{\ast }^2}-1}}\right)]\text{ for }p>R_{\ast }\text{. }$$(A.6)

In these equations $$\alpha =\tan (\frac{\pi }{2}-\frac{1}{2}\arccos \left(\frac{z}{\sqrt{p^2+z^2}}\right)),$$(A.7) $$x_1=\frac{R_{\ast }}{p}+\sqrt{{\left(\frac{R_{\ast }}{p}\right)}^2-1}\text{, and }$$(A.8) $$x_2=\frac{R_{\ast }}{p}-\sqrt{{\left(\frac{R_{\ast }}{p}\right)}^2-1}.$$(A.9)

For the X-ray emission arising from the wind interaction zone, only the cases p ≥ R_* are relevant, whereas for intrinsic emission arising from wind-embedded shocks, all situations described by Eqs. (A.3) – (A.6) matter.

References

All Tables

All Figures

Fig. 1 X-ray spectra of HD 149404. Top panel: combined RGS (RGS1 + RGS2, first and second orders, all three observations) spectrum of HD 149404. The labels indicate the lines that are apparent in this spectrum. The error bar in the top left corner indicates the mean error on the energy bins in the combined RGS spectrum between 10 Å and 20 Å. Bottom panel: comparison of the EPIC-pn spectra of HD 149404 during the three observations. The first observation (φ = 0.227) is shown in blue, the second (φ = 0.992) in red, and the third (φ = 0.487) in green.

In the text

Fig. 2 Variations of the X-ray emission of HD 149404 with orbital phase. From top to bottom, the panels illustrate the full EPIC-pn count rates for the different observations and the observed fluxes evaluated over three different energy bands: soft (0.5–1.0 keV), medium (1.0–2.0 keV), and hard (2.0–4.0 keV). The magenta symbols in the second and third panels correspond to the fluxes measured on an archival ROSAT-PSPC spectrum (Rauw et al. 2001).

In the text

	Fig. 3 Ratio of the observed fluxes evaluated over narrow energy bands. The blue (resp. red) symbols indicate the ratios of the fluxes measured at phase 0.227 (resp. 0.992) to those determined at phase 0.487, the phase at which the emission is strongest.
In the text

	Fig. 4 Product of wind X-ray opacity κE times wind density ρ as a function of energy. The opacity and wind density are computed from the CMFGEN models of Raucq et al. (2016) at a distance r = a/2 from the stellar centre. The continuous red line (labelled 1) shows the results for the primary star, whereas the dashed blue line (labelled 2) corresponds to the secondary star of HD 149404.
In the text

	Fig. 5 Schematic illustration of the wind interaction zone in the orbital plane of HD 149404 modelled with the formalism of Cantó et al. (1996).
In the text

Fig. 6 Comparison between the observed flux ratios and the best-fit predicted values for different models tested in this work. The red and blue symbols with error bars correspond respectively to the measured flux ratios displayed in Fig. 3. The orange and cyan curves yield the best-fit theoretical models for the phases shown respectively by the red and blue symbols. The continuous lines assume the X-rays arise entirely from the wind interaction zone of HD 149404 and using temperature and density profiles computed with the Cantó et al. (1996) model for an adiabatic wind interaction zone (Sect. 4.2). The long-dashed lines yield the results for X-rays arising entirely from wind-embedded shocks intrinsic to the primary and secondary stars (Sect. 4.3).

In the text

	Fig. 7 Theoretical temperature (top panel) and density (bottom panel) maps of the wind interaction zone in HD 149404 computed with the LIFELINE code (Mossoux & Rauw 2021), adopting the parameters of Model III (see Table 5). For this set of parameters, the wind of the primary star (on the right) crashes onto the surface of the secondary (on the left).
In the text