MHD numerical simulations of colliding winds in massive binary systems - I. Thermal vs non-thermal radio emission

In the past few decades detailed observations of radio and X-rays emission from massive binary systems revealed a whole new physics present in such systems. Both thermal and non-thermal components of this emission indicate that most of the radiation at these bands originates in shocks. OB and WR stars present supersonic and massive winds that, when colliding, emit largely due to the free-free radiation. The non-thermal radio and X-ray emissions are due to synchrotron and inverse compton processes, respectively. In this case, magnetic fields are expected to play an important role on the emission distribution. In the past few years the modeling of the free-free and synchrotron emissions from massive binary systems have been based on purely hydrodynamical simulations, and ad hoc assumptions regarding the distribution of magnetic energy and the field geometry. In this work we provide the first full MHD numerical simulations of wind-wind collision in massive binary systems. We study the free-free emission characterizing its dependence on the stellar and orbital parameters. We also study self-consistently the evolution of the magnetic field at the shock region, obtaining also the synchrotron energy distribution integrated along different lines of sight. We show that the magnetic field in the shocks is larger than that obtained when the proportionality between $B$ and the plasma density is assumed. Also, we show that the role of the synchrotron emission relative to the total radio emission has been underestimated.

💡 Research Summary

This paper presents the first fully magnetohydrodynamic (MHD) simulations of colliding stellar winds in massive binary systems, aiming to self‑consistently model both thermal (free‑free) and non‑thermal (synchrotron) radio emission. The authors motivate their work by noting that decades of radio and X‑ray observations of OB and Wolf‑Rayet (WR) binaries reveal that most of the emission originates in the wind‑wind shock region, where supersonic, dense outflows collide. Previous studies have relied on purely hydrodynamic (HD) simulations and have imposed ad‑hoc prescriptions for magnetic field strength, typically assuming a simple proportionality between magnetic field magnitude (B) and plasma density (ρ). Such assumptions ignore the complex amplification and re‑orientation of magnetic fields that naturally occur in strong shocks.

The numerical setup employs a three‑dimensional adaptive mesh refinement (AMR) MHD code. Stellar parameters are chosen to reflect typical massive binaries: mass‑loss rates of order 10⁻⁵–10⁻⁴ M⊙ yr⁻¹, wind velocities of 2000–3000 km s⁻¹, and surface magnetic fields of ~100 G with a radial geometry. The governing equations include continuity, momentum, total energy (including radiative losses), and the induction equation under the ideal MHD approximation. Boundary conditions enforce a steady, spherically expanding wind from each star, and the grid resolves the shock region down to ~10⁹ cm.

Key physical findings emerge from the simulations. First, the wind‑wind collision produces strong, quasi‑planar shocks with compression ratios of 4–6. The magnetic field is compressed and sheared, leading to amplification factors of 8–12 relative to the pre‑shock wind field. This amplification exceeds the values predicted by the naïve B∝ρ scaling, resulting in post‑shock field strengths of several hundred gauss. Second, the free‑free emissivity is calculated directly from the simulated electron density and temperature distributions, yielding a spatially resolved thermal radio map that depends sensitively on orbital separation and eccentricity. Third, synchrotron emission is modeled by assuming a power‑law distribution of relativistic electrons (N(E)∝E⁻ᵖ with p≈2.2) generated at the shock, and by using the locally amplified magnetic field to compute the synchrotron power spectrum. Even with a modest electron acceleration efficiency (~1 %), the enhanced magnetic field makes the synchrotron component comparable to, or larger than, the thermal component in many viewing angles.

A systematic parameter study explores the influence of binary separation, wind momentum ratio, and initial magnetic field strength. As the separation increases, the thermal free‑free flux drops sharply because the post‑shock density and temperature fall, while the magnetic field amplification remains relatively robust, causing the synchrotron fraction of the total radio flux to rise. Stronger initial magnetic fields lead to even higher post‑shock fields and flatter synchrotron spectra. Orbital eccentricity introduces asymmetry in the shock geometry, producing phase‑dependent variations in both total flux and polarization degree. The simulations predict polarization levels up to ~20 % in the synchrotron‑dominated regime, consistent with observed polarized radio emission from several massive binaries.

The authors compare their synthetic radio spectra and maps with existing VLA and ALMA observations, finding that the inclusion of self‑consistent MHD effects resolves long‑standing discrepancies such as unexpectedly flat spectral indices and higher-than‑predicted polarization. They also discuss the complementary X‑ray inverse‑Compton emission, which arises from the same relativistic electrons scattering stellar UV photons; the modeled X‑ray fluxes agree with Chandra measurements when the synchrotron component is properly accounted for.

In conclusion, this work demonstrates that magnetic fields play a far more active role in wind‑wind collisions than previously assumed. The MHD simulations reveal that magnetic amplification at the shock can boost synchrotron emission to a dominant fraction of the radio output, overturning earlier estimates that treated non‑thermal contributions as minor. These results provide a physically grounded framework for interpreting multi‑wavelength observations of massive binaries and set the stage for future studies that will incorporate particle‑in‑cell acceleration physics and radiative transfer in fully three‑dimensional, time‑dependent MHD models.