Colliding-Wind Binaries as a Source of TeV Cosmic Rays

Introduction

Star-forming regions are among the sources of high-energy (HE) photons in our galaxy. The reason for that is still a subject of investigation, given that these regions are populated by a mix of magnetized turbulent plasma, as well as a number of stellar objects at different evolutionary stages, given their initial masses, such as collapsing gas cores of stars yet to form, stars of main sequence, as well as more evolved massive ones, which are already at the edge of being subjects to energetic explosive events that end their evolution. Shocks and magnetic fields, basic ingredients for particle acceleration, are present in all the mentioned scenarios and give a hint of possible generation of HE photons from HE particles by synchrotron, inverse Compton, and pion decay processes. While the magnetized plasma of the interstellar medium is understood as important for the particle diffusion and scattering, stellar sources and compact objects have been mentioned as main drivers of particle acceleration in star-forming regions (SFRs).

Supernova (SNe) remnants have been well associated as sources of HE photons and particles for many decades. Among these, the type-II, Ib, and Ic SNe are strongly related to the loci of star formation, given their progenitor’s high mass. Accretion—and associated disks and jets—is also commonly related to the observed HE photons from SFRs, for example, the X-ray binaries. While low-mass X-ray binaries (LMXB) are spatially distributed over the whole galaxy, the high-mass counterpart (HMXB) is associated to SFRs (see, e.g., Grimm et al., 2002).

Not only massive binaries with compact objects emit HE photons, strong shocks generated by high-velocity winds of both stars in massive binary systems can also accelerate particles to very high energies, being observed in hard X-rays and gamma-rays bands, the so-called colliding-wind binaries (CWBs). These objects are also strongly related to SFRs. Particle acceleration to high energies from CWBs has been suggested by Eichler and Usov (1993). Indeed, nonthermal emission—especially in radio frequencies—associated to high-energy particles interacting with local magnetic fields is often detected from these objects (De Becker and Raucq, 2013).

The theory of diffusive shock acceleration (DSA) is based on the scattering and change of particle momenta as they interact with plasma waves in a second-order Fermi acceleration fashion (e.g., Fermi, 1949; Bell, 1978). As pointed by Eichler and Usov (1993), there is the need for a relatively strong magnetic field at the shock region for the particles to be efficiently accelerated. The minimum strength of magnetic field can be estimated in order for the collisionless regime to occur; that is, the gyrofrequency must be larger than the thermal collision frequency. There is a number of recent works attempting to estimate HE particle populations and the related emission of radiation for CWBs that do not properly consider the distribution of magnetic fields in the shocks (e.g., Dougherty et al., 2003; Reimer et al., 2006; Reitberger et al., 2014b, a; del Palacio et al., 2016; Pittard et al., 2020).

These previous works, although very useful to generally describe spectral features of nonthermal and high-energy sources associated to CWBs, have not been able to provide a proper answer to two major open issues in this field: 1) how the particles are accelerated in CWBs—which includes how fast this acceleration occurs and how they diffuse, and 2) where these particles are mainly accelerated. To answer these two questions, a more complex analysis is needed, based on full magnetohydrodynamics (MHD) combined with a full treatment of particle acceleration either by evolution of momentum distributions or, even better, by studying the kinematics of individual particles.

Full MHD simulations were first employed for the study of nonthermal emission from CWBs by Falceta-Gonçalves and Abraham (2012), which showed that parallel shocks are unrealistic, and magnetic field amplification at the downstream is much larger than expected from adiabatic modeling. Full MHD simulations of CWBs were further presented by Kissmann et al. (2016). Geometry complexity and amplification are products of strong cooling of the post-shocked plasma, which brings a scenario much different from the ones typically assumed for semi-analytical approximations. In Falceta-Gonçalves and Abraham (2012), however, the particle acceleration was treated in a simplified model, and the energy distribution of particles could not be properly determined. Falceta-Gonçalves (2015) studied the first- and second-order Fermi acceleration processes finding that second-order effects are dominant for their MHD numerical simulations of eta Carinae. Energies as high as 1–10 TeV, at timescales of 10⁴–10⁵ s, were given as parametric estimates for that particular MHD setup of eta Carinae, although the energy distribution of particles was not determined.

Given the context described above, a natural improvement of the theoretical modeling of particle acceleration in CWBs would be the development of a full kinetic description of particles as subject to magnetized plasma fluctuations derived from a full MHD simulation of the shocks. This is the main objective of the present work. Here, we provide the first model to integrate a population of thermal hadronic particles, subject to the local distributions of plasma velocity, density, and magnetic fields derived from a full MHD numerical simulation. The organization of the manuscript is given as follows.

In Numerical Modeling, we describe the magnetohydrodynamic (MHD) model of the CWB system we simulated and the method to evolve acceleration and diffusion of test particles in this system. In Results, we describe the results of the MHD model and the most important characteristics of nonthermal particles accelerated in the studied system. We also discuss the spatial distribution of the acceleration rate. In Discussion and Conclusions, we briefly discuss our results and draw our main conclusion.

Numerical Modeling

Magnetohydrodynamic Model of Colliding-Wind Binary

We model a binary system consisting of primary and secondary stars placed at a distance of approximately 7.2 AU. The stars are described numerically by spherical regions of the radius of about 0.12 AU at which the density, pressure, outflow velocity, and magnetic field are maintained during the whole simulation, and determined by emitted stellar winds. The wind density of the primary star is set to ${\rho }_m={10}^7$ in the code units (c.u.), while in the secondary star, it is 10² times smaller. We assumed the wind speeds of the primary and secondary stars to be 500 km/s and 2500 km/s, respectively. The wind pressure at the star surfaces is determined by the formula $p_{wind}={\left(v_{wind}/M_{wind}\right)}^2{\rho }_{wind}/\Gamma$, where $M_{wind}$ is the sonic Mach number of the wind at the stellar surface and Γ is the adiabatic index set to 1.1. In our model, we set $M_{wind}$ to 3 and 15 for the primary and secondary stars, respectively. Although magnetized, the massive stellar winds are super-Alfvénic, that is, the magnetic pressure is smaller than the kinetic one. As a consequence, near the star, the magnetic field lines follow the velocity field profile, which is radial. The radially oriented fields follow a $r^{-2}$ scaling, anchored at the stellar surface foot points. The foot points are constant in time and are determined considering an internal magnetic field of ${10}^{-3}$, in the code units, oriented along the Z direction. The surface and internal magnetic fields are kept constant during the simulation, while the plasma outside of the stars is evolved according to the MHD equations. For the sake of setup simplicity, we also set initially the magnetic field of the interstellar gas as uniform, with the same intensity. This, however, is irrelevant for the results presented here as the whole plasma of the interstellar medium is pushed out of the computational domain as the winds expand and reach the quasi–steady-state condition. It is important to notice that the interstellar plasma parameters analyzed further in this work are therefore independent of this external initial setup, depending only on the stellar surface wind parameters.

The above setup is evolved by the ideal adiabatic compressible MHD system of equations:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathit{u}\right)=0,(1)$

$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla \cdot \left[\rho \mathit{u}\mathit{u}+\left(p+\frac{B^2}{2}\right)I-\mathit{B}\mathit{B}\right]=0,(2)$

$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times \left(\mathbf{u}\times \mathbf{B}\right)=0,\nabla \cdot \mathbf{B}=0,(3)$

$\frac{\partial E}{\partial t}+\nabla \cdot [\left(E+p\right)\mathbf{u}-\left(\mathbf{u}\cdot \mathbf{B}\right)\mathbf{B}]=0,(4)$

where ρ, p, $\mathbf{u}$, and $\mathbf{B}$ are the plasma density, thermal pressure, velocity, and magnetic field, respectively, $E=(\rho u^2/2)+p/\left(\Gamma -1\right)+\left(B^2/2\right)$ is the total energy, and I is the identity matrix in a 3D rectangular domain described on the Cartesian coordinate system with physical dimensions 9L×6L×6L (assuming L ≈ 1.2 AU) corresponding to the X, Y, and Z directions, respectively, the base resolution of 48×32×32, and the local mesh refinement up to the fifth level resulting in the effective resolution of 768×512×512 (the effective grid size being h ≈ 0.012 c. u. ≈ 0.014 AU along all directions). The refinement criterion was based on the normalized value of vorticity and current density to follow the shock, turbulence, and magnetic reconnection development. We also set the outflow boundaries along all directions.

The MHD equations were solved numerically using a high-order shock-capturing adaptive refinement Godunov-type code AMUN¹ with 5^th order optimized compact monotonicity preserving reconstruction of Riemann states (Ahn and Lee, 2020), the HLLD Riemann flux solver (Mignone, 2007), and the 3^rd order 4-step Strong Stability Preserving Runge–Kutta (SSPRK) method (Gottlieb et al., 2011) for time integration. The system was evolved until t = 21 c. u. (corresponding to roughly 44 days considering the assumed length and velocity units), at which the interaction surface between the colliding winds stopped moving through the domain, and a turbulent region produced by their interaction became statistically stationary. Once such a state is reached, the system can be evolved for an arbitrary period of time maintaining its statistical properties.

Test Particle Integration

In order to study the behavior of test particles in the binary system evolved as described in the previous subsection, we took the last available domain snapshot (velocity and magnetic field spatial distributions) at t = 21 and injected 10⁵ (if not specified otherwise) protons in random locations near the turbulent region created by colliding winds with random orientations of velocity and the initial particle velocity distribution corresponding to the thermal distribution at temperature T = 10⁶ K.

The trajectories of injected particles were integrated using the relativistic equation of motion for a charged particle

$\frac{d}{dt}\left(\gamma \mathit{v}\right)=\frac{q}{m}\left(\mathit{E}+\mathit{v}\times \mathit{B}\right),\frac{d\mathit{x}}{dt}=\mathit{v},(5)$

where $\mathbf{x}$ and $\mathbf{v}$ are the particle position and velocity, respectively, $\gamma ={\left(1-v^2/c^2\right)}^{1/2}$ is the Lorentz factor, q and m are particle charge and mass, respectively, c is the speed of light, and $\mathbf{E}=-\mathbf{u}\times \mathbf{B}$ is the plasma-induced electric field, with $\mathbf{u}$ and $\mathbf{B}$ being the plasma velocity and magnetic field taken from the grid simulation, respectively.

We assumed the particle velocity $\mathbf{v}$ in Equation. 5 to be expressed in units of the speed of light c and its position in the coordinates of the grid simulation. Since the grid simulations use dimensionless units, one needs to assume the unit of plasma velocity $\mathbf{u}$ and magnetic field $\mathbf{B}$. In this work, we assumed that the plasma velocity is in units of 1000 km/s, and for the magnetization level of the winds, we set the magnetic field anchored at the stellar surfaces to range from B_* = 10⁻⁵ to 10⁰ G. Both stars are assumed to maintain the same intensity of surface magnetic fields during the whole simulation time. This value is a model parameter and should not be confused with the real stellar magnetic strength. Due to the numerical modeling limitations, it simply represents the magnetic field which is ejected from the star together with the wind. We should also note that in the integration, the unit of time was chosen as 50 h = 1.8×10⁵ s which, assuming the plasma velocity unit of 1000 km/s, results in the unit of length equal to L ≈ 1.2 AU.

The integration of Equation. 5 was done with the help of the code GAccel² using the 8^th order embedded Dormand–Prince (DOP853) method (Hairer et al., 2008) with an adaptive time step based on an error estimator. In order to determine the values of $\mathbf{u}$ and $\mathbf{B}$ at the current particle positions, we used tri-cubic interpolation (a cubic interpolation applied dimension by dimension), enforcing the monotonicity by limiting the tangents at the interpolation nodes with the minmod TVD limiter (van Leer, 1979). In the integration of the particle trajectories, we set the absolute and relative tolerance to atol = rtol = 1e-10 and the initial time step dt = 1e-6, if not specified otherwise.

Results

Characteristics of the Colliding-Wind System

The model of colliding-wind binary described in the MHD Model of Colliding-wind Binary was evolved up to the equilibrium position of the contact discontinuity of the shock apex (t ≈ 21 c. u.). At this time, the region at which the winds interact reach stochastic saturation. It is characterized by the presence of strong turbulence produced by the interacting winds and magnetohydrodynamic instabilities. This region is also characterized by strong compression, allowing for the increase of the magnetic field energy to the levels above the magnetic energy maintained within the stars.

In Figure 1, we show a 3D visualization of the magnetic energy distribution at time t = 21. The magnetic energy maintained in both stars is equal to 5×10⁻⁷ c.u. (white, in the middle of the presented color bar). However, $E_{mag}$ quickly decays as the magnetic field is transferred by the wind from the stellar surfaces and distributed into the surroundings (see blueish envelopes around the stars). Since the winds are super-Alfvénic, the field decays radially, that is, $r^{-2}$. It is again compressed in the region of the wind interaction forming strongly magnetized structures (reddish clumps in Figure. 1). The maximum magnetic field strength in these structures could reach values in an order of magnitude higher than those set at the stellar surfaces. The distribution of these magnetic clumps is clearly related to the turbulent region produced by the interacting winds.

FIGURE 1

FIGURE 1. Distribution of magnetic energy in the computational domain of the studied colliding-wind binary system at time t = 21. On the right, we indicate the primary star, and on the left, the secondary one. The turbulent region produced by the interacting winds is characterized by a number of clumps of increased magnetic energy. Regions of magnetic energy larger than the energy within the stars are shown in red.

In Figure 2, we show the distribution of the divergence of velocity $\left(\nabla \cdot \mathbf{u}\right)$ in the same projection as the previous figure. For clarity, the values have inverted sign, so, for example, the value 100 indicates $\left(\nabla \cdot \mathbf{u}=-100\right)$. It is evident that the negative part of the divergence of plasma velocity is a good tracer for determining the extension of the volume affected by the interaction of the colliding winds. Expressing the largest values seen in Figure 2 in the units of u = 1000 km/s and L ≈ 1.2 AU gives the maximum compression rate of around 3.46 h^–1, indicating high compression efficiency of these interactions. The regions of compressed gas relax the increased pressure converting it into the mechanical energy and producing turbulence and pushing away the turbulence structure. This compression–relaxation process can be maintained in such a system as long as the stars emit their powerful winds.

FIGURE 2

FIGURE 2. Distribution of the divergence of velocity in the computational domain at time t = 21 in the same projection as Figure. 1. Only the positive values of $-\nabla \cdot \mathbf{u}$ are shown. The volume where the gas is affected by the wind interaction can be clearly seen. It is characterized by the development of strong compressing shocks and turbulence.

In order to allow for comparison of observations with our results on the particle acceleration and maximum energies reached, which are presented in the following sections, we have estimated the mean magnetic field ${\overline{B}}_{shock}$ within the shock region of the performed numerical modeling of CWBs. The region was determined by the condition $\left(\nabla \cdot \mathbf{u}\le -50\right)$. We have obtained the mean intensity of the magnetic field in the shock region to scale as ${\overline{B}}_{shock}\sim {10}^{-1}{B}_{\mathrm{*}}$, where, as mentioned earlier, $B_{\mathrm{*}}$ is the model parameter determining the strength of the magnetic field at the distance of 0.12 AU from the star center. The variance of the magnetic field strength in the shock region was estimated to be of the same order as the mean value.

Cosmic Rays Produced by the Colliding-Wind System

As we explained in Test Particle Integration, the protons are injected within the turbulent region of the interacting winds. More precisely, we injected ${10}^5$ protons into the region $[-1.5,-0.5]\times [-1.0,1.0]\times [-1.0,1.0]$ of the computational domain corresponding to region limits along the x, y, and z coordinates, respectively. The protons are integrated up to $t_{\max }=5\times {10}^4\text{hours}$ or until they reach the boundary of the domain leaving the system, whichever occurs first.

As an example, we provide a movie demonstrating the evolution of protons during the time of integration for the model with B_* = 10⁻¹ in the Supplementary Material. In the movie protons are represented as color spheres, with the color corresponding to the kinetic energy level E_k. The integration time elapses in the logarithmic scale.

In Figure 3, we show the evolution of distributions of a few most important properties characterizing the protons: 1) the evolution of particle kinetic energy distribution, 2) the evolution of particle gyroradius in the code units, 3) and 4) the evolution of parallel and perpendicular, with respect to the local magnetic field, components of velocity in the units of c. Figure 3a demonstrates that the protons suffer initially a very quick, within a few seconds, pre-acceleration, most probably due to the electromagnetic plasma drift or magnetic field curvature. At this stage, they can reach the kinetic energies up to a few hundreds of MeV. In our analysis, this value was insensitive to the value of $B_{\mathrm{*}}$ used in the integration. The distribution of gyroradia does not change significantly, although some sort of a reorganization is seen characterized by a population of protons for which the gyroradius increases (see the period of time before t = 100 s in Figure. 3b). The increase of gyroradia seems to be related to the increase of perpendicular velocities, since similar patterns are observed in Figure 3d during this period. Contrary to the perpendicular component, the parallel component of particle velocity is essentially unaltered.

FIGURE 3

FIGURE 3. Evolution of distributions of the particle kinetic energy (A), the gyroradius (B), and the parallel and perpendicular particle velocities (C) and (D) for the model with $B_{\mathrm{*}}={10}^{-1}\text{G}$. In the energy plot, the gray area indicates the range of the maximum kinetic energy reached by 1% of the most energetic protons. In the gyroradius plot, the gray area indicates the scales permitted by the MHD simulation, with h being the effective cell size $h=1/512$ and $L_y=6L\approx 7.2 \text{AU}$ the size of the domain in the Y and Z directions.

The situation significantly changes at around t = 200 s when the parallel component starts to increase, initially to around one percent of c between t = 10³ s and 10⁴ s, and eventually reaching nearly the speed of light between t = 10⁵ s and 10⁶ s (see Figure. 3c). At this stage, that is, after t = 10⁵ s, one can notice a significant increase of the perpendicular velocity component, reaching the speed of light, as well. Clearly, the increase of the perpendicular velocity results in the increase of the gyroradius, observed in Figure 3b. It is interesting to see, however, that even due to relativistic particle speed, only a small fraction of protons has gyroradia larger than the grid size h, barely reaching the specified length unit of L ≈ 1.2 AU.

As seen in Figure 3a, after t = 10⁴ s, a significant fraction of protons can increase their kinetic energies, defined as $E_k\equiv \left(\gamma -1\right)m_p{c}^2$, where $m_p$ is the proton mass, by several orders of magnitude, from levels below $E_k\sim 0.1\text{MeV}$ to tens and hundreds of GeV. The gray area seen in this plot indicates the range of kinetic energies reached by 1% of the most energetic protons. We should clarify that this range does not represent the energies at which protons left the computational domain but the maximum energy they reached during the whole integration. As we explained above, the energetic level after the pre-acceleration at $t<{10}^4\text{\hspace{0.17em}s}$ is essentially insensitive to the used strength of the magnetic field $B_{\mathrm{*}}$. As we will see below, the energy band of the most energetic protons, on the contrary, is strongly sensitive to $B_{\mathrm{*}}$.

In the left plot of Figure 4, we show distributions of the maximum energies reached by protons for all studied particle integration cases, that is, for $B_{\mathrm{*}}={10}^{-5}\text{G}$ to $1\text{G}$. As one can see, the low-energy population ${\text{E}}_k<{10}^4\text{eV}$ is very similar for all cases. Also, the second more energetic population of protons is relatively insensitive to the magnetic field strength. Nevertheless, from around the energy level of ${10}^6\text{eV}$ up, one can identify that the high-energy tail of these distribution strongly depends on $B_{\mathrm{*}}$. The vertical dashed lines indicate the maximum energy reached in each case. The relation between these energies and the values of ${\overline{B}}_{shock}$, the mean field in the shock region, with its dispersion as horizontal error bars, is shown in the right plot of Figure 4. It is interesting to see that the maximum energy correlates nearly linearly with the magnetic field strength. Moreover, from this plot, we can see that around 0.1% or more of very energetic protons in the range of 100 GeV to one TeV can be produced for ${\overline{B}}_{shock}={10}^{-3}\text{G}$ or stronger. This dependence also allows to extrapolate the maximum energies for colliding-wind systems characterized by even stronger magnetic fields, indicating that at least one particle of the energy over 10 TeV could be produced once every around 4 months (most protons typically leave the system within ${10}^7\text{s}$, as seen in Figure 3).

FIGURE 4

FIGURE 4. Left: Distribution of the particle largest kinetic energy for integration done with different values of $B_{\mathrm{*}}$, from ${10}^{-5}$ to ${10}^0\text{G}$. Vertical gray-dashed lines indicate the maximum kinetic energies reached in each integration. Right: Dependence of the maximum kinetic energy (indicated by the vertical lines in the left panel) on the mean magnetic field in the shock region ${\overline{B}}_{shock}$.

The maximum energies shown here were obtained neglecting any process of energy loss. We are basically interested in the acceleration of hadrons, that is, the typical electron energy losses such as synchrotron emission, inverse Compton, Coulomb collisions, and bremsstrahlung emission losses are unimportant and may be disregarded. The protons are expected to suffer energy loss mainly due to hadronic collisions. The process of pion generation, and further decay, is the main energy loss phenomenon at the energy range obtained from our simulations. In order to validate our results, we may estimate the loss rate for proton–proton collisions, as given by Kelner et al. (2006) and Pittard et al. (2020):

$\frac{d\gamma }{dt}\simeq \frac{1}{2}{n}_p{\gamma }_{p}c{\sigma }_{pp},(6)$

with ${\sigma }_{pp}$being

${\sigma }_{pp}=\left(34.3+1.88\lambda +0.25{\lambda }^2\right){\left[1-{\left(\frac{1.22GeV}{E_p}\right)}^2\right]}^2\times \left({10}^{-27}c{m}^2\right),(7)$

where $\lambda =\ln \left(E_p/1\text{TeV}\right)$, $E_p\equiv \gamma m_p{c}^2$ is the particle total energy, and $n_p$ is the proton number density. Equation 6 was solved for each time step of the particle integration procedure. The deceleration rates $d\gamma /dt$ for the nonthermal energetic particles result in maximum energy loss scaling as $\sim {10}^{-10}{n}_p\left({\overline{B}}_{shock}/1\text{G}\right)$ GeV/s, where ${\overline{B}}_{shock}$ is the mean magnetic intensity in the shock region. The dependence on the thermal proton number density $n_p$ has been kept in order to allow us to estimate under what conditions the energy loss rate would be relevant compared to the acceleration process.

As discussed in the next subsection, the acceleration rates for these particles lie in the range of few to several GeV/s. In order to halt the acceleration process, the hadronic collisions will be dominant, in a mG (miliGauss) shocked region, only for proton number densities above ${10}^{13}$ cm⁻³, which is well above any typical post-shock condition expected for the massive systems known. Therefore, the estimates of acceleration and energy distributions obtained here are likely to be independent on these loss terms.

The analysis of the maximum energies of protons produced in the colliding-wind system gives an insight into the possible potential of the system in the production of high or very high energetic protons. In order to analyze how the thermal protons reach higher energies, we have to analyze the particle distributions at given times. In the left plot of Figure 5, we demonstrate the distributions of the kinetic energy for the integration done with $B_{\mathrm{*}}={10}^{-2}\text{G}$ for three different moments: the initial thermal distribution of $T={10}^6\text{K}$ (gray line), at the intermediate stage at t = 5 h when the high-energy power-law tail starts to develop, and the distribution at t = 50 h, just before the protons start leaving the computational domain. We can notice that the purely thermal protons are initially heated toward slightly higher temperatures and develop a nonthermal energetic tail. At t = 50.0, the distribution is characterized by two thermal populations of protons (fitted by green and red dashed curves) with temperatures estimated to around $4.5\times {10}^6\text{K}$ and $5.9\times {10}^7\text{K}$ and a nonthermal power-law tail, $dN\propto E^{\alpha }dE$, with the index estimated to be around $\alpha \approx -0.41$ (gray dashed line). In the right plot of Figure 5, we show the dependence of the fitted power-law index on mean magnetic field strength of the turbulent region, ${\overline{B}}_{shock}$. The index depends on the value of ${\overline{B}}_{shock}$ decaying from around 0.66 for ${\overline{B}}_{shock}={10}^{-6}\text{G}$ to around 0.18 for ${\overline{B}}_{shock}=0.1\text{G}$. The flattening of the power-law tail is well justified, considering that in all cases the initial thermal distribution has the same temperature and the energy band of the tail, as well as the energy gain efficiency, both increase with the value of ${\overline{B}}_{shock}$.

FIGURE 5

FIGURE 5. Left: Particle kinetic energy distribution at moments $t=0$, $t=5$, and $t=50$ hours (gray, blue, and orange, respectively) for the integration done with $B_{\mathrm{*}}={10}^{-2}\text{G}$. At the initial time, the distribution is purely thermal with temperature $T={10}^6\text{K}$. At the final time shown, the distribution is characterized by two thermal populations with $T\approx 4.5\times {10}^6\text{K}$ (green dashed line) and $T\approx 5.9\times {10}^7\text{K}$ (red dashed line) and the high-energy power-law tail (gray dashed line) with the fitted index of around $-0.4$. Right: The dependence of the power-law index fitted to the nonthermal tail on mean magnetic field in the shock region ${\overline{B}}_{shock}$ at the same moment $t=50$.

Observations of nonthermal emission from WR-O binary systems have been very successful in the radio wavelengths, typically associated to the synchrotron emission from relativistic electrons. The nonthermal leptonic emission at radio wavelengths presents spectral indices as steep as $\alpha \sim -0.62$ with a few systems, where α is in the range from −0.4 to −0.2 (Dougherty and Williams, 2000). Radio emission observations of the galactic super star cluster Westerlund 1 revealed some nonthermal radio sources with the spectral index α in the range from −0.6 to −0.06 (Dougherty et al., 2010). Since we integrated hadron trajectories in our simulations, we cannot overinterpret our results for electrons and related synchrotron emission spectra.

Emission from relativistic hadrons however, such as nonthermal bremsstrahlung and γ-rays from pion decay, has only been systematically observed for the system of η Carinae. Fermi-LAT detected very energetic emission $>10$ GeV close to the periastron passage (Abdo et al., 2010). Emission from ηCar, in the range of 100 GeV up to ∼400 GeV, has also been detected by H.E.S.S. (H. E. S. S. Collaboration et al., 2020). The combined data in the range from few up to hundreds of GeVs show a spectral index of −0.7 to −0.5. Since the proton-proton interactions generate γ-ray spectrum with a slope similar to that of the energy distribution of the parent protons, we may compare these observation estimates to our results shown in Figure 5b. It suggests that by varying the typical magnetic field intensity of the stellar surface $B_{\mathrm{*}}$, we can also obtain the spectral indices in a compatible range from −0.2 to −0.7. Clearly, we do not take into account the particle energy losses for obtaining the particle energy spectra, but still the fact that we could obtain similar spectral indices is encouraging and suggests future studies taking into account the radiative energy losses and transfer, or more complex magnetic field stellar geometries.

Acceleration Rate Distribution

Besides the basic properties of the acceleration mechanism such as the maximum kinetic energy reached by the protons, their distribution, and the index of the nonthermal tail, it is also interesting to determine what the range of the acceleration rates available in the given system is and the spatial distribution of the sites where the efficient acceleration is favorable.

In Figure 6, we show the spatial distribution of such acceleration sites in our computational domain for the integration of protons done with $B_{\mathrm{*}}={10}^{-1}\text{G}$. There is a clear correlation of these sites with the turbulent region produced by the interacting winds. Protons crossing the sites shown with the red color could gain the energies at rates as large as one TeV/h, with the maximum acceleration rate determined for this case reaching nearly 34 TeV/h $(\sim 9.3\text{GeV/s})$.

FIGURE 6

FIGURE 6. Distribution of the maximum acceleration rate $d{E}_k/dt$ in the computational domain for the integration done with $B_{\mathrm{*}}={10}^{-1} \text{G}$. Regions indicated by red color correspond to sites where the detected maximum energy gain rate was above one TeV/h. The maximum acceleration rate detected in this model was nearly 34 TeV/h.

In Figure 7, we show the distributions of the acceleration rate for all studied particle models. There are a few important characteristics seen in these distributions. First of all, the maximums are shifted toward higher values of the acceleration rates for stronger magnetic fields. Moreover, the bands where the most common rates are located depend in the magnetic field strength $B_{\mathrm{*}}$. With the increasing $B_{\mathrm{*}}$, the acceleration rates become more efficient, reaching orders of hundreds of GeV per second for the most magnetized case $B_{\mathrm{*}}=1\text{G}$. We should note that in this case, we were able to integrate only 4,000 protons, which explains the visible deviation of this distribution from other cases.

FIGURE 7

FIGURE 7. Probability distribution function of the maximum acceleration rate $d{E}_k/dt$ for the integration done with $B_{\mathrm{*}}$ from 10⁻⁵ to 10⁰ G. Relatively high fraction of the computational volume is characterized by higher values of the acceleration rates, from keVs to GeVs per second. The maximum acceleration rate is as high as nearly 80 GeV/s for the model with $B_{\mathrm{*}}=1\text{G}$.

Discussion and Conclusion

In this work, we have investigated a colliding-wind binary system as a possible site of high-energy and very high–energy cosmic ray production. We have performed the magnetohydrodynamic simulations of interacting high velocity and high mass loss rate stellar winds from a typical massive binary system. The winds are characterized by supersonic speeds, 2500 km/s and 500 km/s, for the primary and secondary stars, respectively. Collision of such strong magnetized winds causes the production of shocks, turbulence, and magnetic reconnection. These constitute fundamental ingredients for efficient diffusive acceleration. The magnetic reconnection has been proposed as an alternative acceleration mechanism, also represented by the first-order process, gaining more attention recently (see, e.g., de Gouveia dal Pino and Lazarian, 2005; Lyubarsky and Liverts, 2008; Kowal et al., 2011, 2012; de Gouveia Dal Pino and Kowal, 2015; del Valle et al., 2016; Comisso and Sironi, 2019).

In order to study the energy gain of nonthermal protons in such system, we injected test particles and integrated their trajectories using the relativistic equation of the charged-particle motion. We present the evolution of distributions of the particle kinetic energy, gyroradia, and parallel and perpendicular, with respect to the local magnetic field, components of velocity. In order to get a deeper insight in the possible energies, the protons can reach in such a system where we estimate how the maximum kinetic energy and nonthermal power-law index depend on the average magnetic field strength in the shock region ${\overline{B}}_{shock}$. Finally, we determine the spatial distribution of the maximum acceleration rate and its typical and maximum values.

Based on our results, we draw the following conclusions:

• The colliding-wind binary studied in this manuscript is able to produce strongly magnetized turbulent region with shock being able to compress the magnetic field to values higher than the magnetic field maintained within the stars. For the limited resolution used in our modeling, the maximum magnetic field strength in the system is about an order of magnitude higher than the magnetic strength, $B_{\mathrm{*}}$, at a distance 0.12 AU from the star. Such a turbulent region surrounded by incoming plasma favors the production of nonthermal protons.

• The highest available energies achieved by the accelerated protons depend on the strength of magnetic field. For ${\overline{B}}_{shock}\ge {10}^{-2}\text{G}$, 1% of most energetic protons, a significant amount, reached energies between tens and hundreds of GeV. For stronger magnetic field, the level between hundreds of GeV and a few TeV was observed. We estimated that the relation between the maximum energy particle and magnetic field scales nearly as $E_{\max }\propto {\overline{B}}_{shock}$. This result confirms that acceleration of thermal protons up to the range of 1–10 TeV is feasible in the colliding-wind binary systems.

• For all studied values of magnetic field intensity, we observe the production of nonthermal tail in the energy distribution with the index decreasing from around 0.7 for the weakest ${\overline{B}}_{shock}={10}^{-6}\text{G}$ to around 0.2 for ${\overline{B}}_{shock}=0.1\text{G}$.

• The regions of the highest acceleration rates $d{E}_k/dt$ are essentially limited to the turbulent region produced by the colliding winds. In the case of the strongest magnetic field studied here, $B_{\mathrm{*}}=1\text{G}$, the PDF of $d{E}_k/dt$ shows the more common range between a few tens of MeV/s and a few GeV/s, with the maximum value reaching nearly 80 GeV/s.

Our results may be compared to those of previous works as follows. There is no record in the literature of full kinematic evolution of thermal test particles interacting with the shocked MHD plasma of a CWB with the objective of understanding the acceleration process of high-energy particles. Estimates of maximum energies and slopes, based on the first-order Fermi acceleration process, as well as DSA, have been posed in the past (e.g., Eichler and Usov, 1993; Dougherty et al., 2003; Reimer et al., 2006; Falceta-Gonçalves and Abraham, 2012; Reitberger et al., 2014b, a; Falceta-Gonçalves, 2015; del Palacio et al., 2016; Pittard et al., 2020; and others, see references therein).

Our findings of energies up to tens of TeV are in agreement with the upper limit DSA estimates of Falceta-Gonçalves (2015), based on MHD simulations, as well as to the approximate solutions of the transport equation, based on HD simulations with assumed magnetization levels, given by Reitberger et al. (2014b). However, the latter have assumed much stronger surface magnetic fields, in the order of 0.01 T. Some discrepancies are expected, given that our models contain self-consistent turbulent magnetic fields generated in the shocks. With respect to the particle energy spectrum of the protons, Reitberger et al. (2014b) obtained values in the range of −0.2 to 0.0, in good agreement with our models with stronger magnetic fields.

In the forthcoming studies, it is important to test the particle acceleration in CWBs for particles other than proton. Therefore, in the first step, we plan to determine the dependence of the maximum energy limit on the ratio $q/m$, extending it toward heavier particles, such as carbon, as well as toward lighter ones, like electrons. We expect that the processes responsible for energy losses would be more important in the case of lighter particles, especially electrons. Therefore, we plan to include some of the most important processes, such as synchrotron, inverse Compton, and bremsstrahlung losses, directly in the integration equations. In future works, it would be mandatory to perform these models with a larger ensemble of particles in order to confirm the obtained results, and more importantly to benchmark upper limit estimates, as well as semi-analytic transport equation calculations for the simulation of nonthermal emission from CWBs. Our results may also be directly compared to the observations in the near future, for example, the Cherenkov Telescope Array (CTA) (Chernyakova et al., 2019).

Author’s Note

Visualization of the protons’ motion in the model with $B_{\mathrm{*}}={10}^{-1}$. 100,000 protons are represented as color spheres, with the color corresponding to the kinetic energy level $E_k$. Please note, that the time elapses in the logarithmic scale.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

DF-G formulated the scientific problem of particle acceleration in CWB. GK prepared and performed numerical simulations of CWB model and the integration of test particle trajectories. Both authors analyzed and interpreted the results and contributed equally to the writing of the article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

GK acknowledges support from the Brazilian National Council for Scientific and Technological Development (CNPq no. 304891/2016–9) and FAPESP (grants 2013/10559–5 and 2019/03301–8). DF-G thanks the Brazilian agencies CNPq (no. 311128/2017–3) and FAPESP (no. 2013/10559–5) for their financial support. This work has made use of the computing facilities: the cluster of Núcleo de Astrofísica Teórica, Universidade Cruzeiro do Sul, Brazil.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2021.667805/full#supplementary-material

Footnotes

¹AMUN code is open source and freely available from https://gitlab.com/gkowal/amun-code.

²GAccel code is open source and freely available from https://gitlab.com/gkowal/gaccel.

References

Abdo, A. A., Ackermann, M., Ajello, M., Allafort, A., Baldini, L., Ballet, J., et al. (2010). Search for Gamma-Ray Emission from Magnetars with the Fermi Large Area Telescope. Astrophysical J. Lett. 725, L73–L78. doi:10.1088/2041-8205/725/1/L73

CrossRef Full Text | Google Scholar

Ahn, M.-H., and Lee, D.-J. (2020). Modified Monotonicity Preserving Constraints for High-Resolution Optimized Compact Scheme. J. Sci. Comput. 83, 34. doi:10.1007/s10915-020-01221-0

CrossRef Full Text | Google Scholar

Bell, A. R. (1978). The Acceleration of Cosmic Rays in Shock Fronts - I. Monthly Notices R. Astronomical Soc. 182, 147–156. doi:10.1093/mnras/182.2.147

CrossRef Full Text | Google Scholar

Chernyakova, M., Malyshev, D., Paizis, A., La Palombara, N., Balbo, M., Walter, R., et al. (2019). Overview of Non-transient γ-ray Binaries and Prospects for the Cherenkov Telescope Array. Astron. Astrophysics 631, A177. doi:10.1051/0004-6361/201936501

CrossRef Full Text | Google Scholar

Collaboration, H. E. S. S., Abdalla, H., Adam, R., Aharonian, F., Ait Benkhali, F., Angüner, E. O., et al. (2020). Detection of Very-High-Energy γ-ray Emission from the Colliding Wind Binary η Car with H.E.S.S. Astron. Astrophysics 635, A167. doi:10.1051/0004-6361/201936761

CrossRef Full Text | Google Scholar

Comisso, L., and Sironi, L. (2019). The Interplay of Magnetically Dominated Turbulence and Magnetic Reconnection in Producing Nonthermal Particles. Astrophysical J. 886, 122. doi:10.3847/1538-4357/ab4c33

CrossRef Full Text | Google Scholar

De Becker, M., and Raucq, F. (2013). Catalogue of Particle-Accelerating Colliding-Wind Binaries. Astron. Astrophysics 558, A28. doi:10.1051/0004-6361/201322074

CrossRef Full Text | Google Scholar

de Gouveia Dal Pino, E. M., and Kowal, G. (2015). Part. Acceleration by Magn. Reconnection 407, 373. doi:10.1007/978-3-662-44625-6_13

CrossRef Full Text

de Gouveia dal Pino, E. M., and Lazarian, A. (2005). Production of the Large Scale Superluminal Ejections of the Microquasar GRS 1915+105 by Violent Magnetic Reconnection. Astron. Astrophysics 441, 845–853. doi:10.1051/0004-6361:20042590

CrossRef Full Text | Google Scholar

del Palacio, S., Bosch-Ramon, V., Romero, G. E., and Benaglia, P. (2016). A Model for the Non-thermal Emission of the Very Massive Colliding-Wind Binary HD 93129A. Astron. Astrophysics 591, A139. doi:10.1051/0004-6361/201628264

CrossRef Full Text | Google Scholar

del Valle, M. V., de Gouveia Dal Pino, E. M., and Kowal, G. (2016). Properties of the First-Order Fermi Acceleration in Fast Magnetic Reconnection Driven by Turbulence in Collisional Magnetohydrodynamical Flows. Monthly Notices R. Astronomical Soc. 463, 4331–4343. doi:10.1093/mnras/stw2276

CrossRef Full Text | Google Scholar

Dougherty, S. M., Clark, J. S., Negueruela, I., Johnson, T., and Chapman, J. M. (2010). Radio Emission from the Massive Stars in the Galactic Super Star Cluster Westerlund 1. Astron. Astrophysics 511, A58. doi:10.1051/0004-6361/200913505

CrossRef Full Text | Google Scholar

Dougherty, S. M., Pittard, J. M., Kasian, L., Coker, R. F., Williams, P. M., and Lloyd, H. M. (2003). Radio Emission Models of Colliding-Wind Binary Systems. Astron. Astrophysics 409, 217–233. doi:10.1051/0004-6361:20031048

CrossRef Full Text | Google Scholar

Dougherty, S. M., and Williams, P. M. (2000). Non-thermal Emission in Wolf-Rayet Stars: Are Massive Companions Required?. Monthly Notices R. Astronomical Soc. 319, 1005–1010. doi:10.1046/j.1365-8711.2000.03837.x

CrossRef Full Text | Google Scholar

Eichler, D., and Usov, V. (1993). Particle Acceleration and Nonthermal Radio Emission in Binaries of Early-Type Stars. Astrophysical J. 402, 271. doi:10.1086/172130

CrossRef Full Text | Google Scholar

Falceta-Gonçalves, D. (2015). “Magnetic Fields, Non-thermal Radiation and Particle Acceleration in Colliding Winds of WR-O Stars,” in Wolf-Rayet Stars. Editors W.-R. Hamann, A. Sander, and H. Todt, 289–292.

Google Scholar

Falceta-Gonçalves, D., and Abraham, Z. (2012). MHD Numerical Simulations of Colliding Winds in Massive Binary Systems - I. Thermal versus Non-thermal Radio Emission. Monthly Notices R. Astronomical Soc. 423, 1562–1570. doi:10.1111/j.1365-2966.2012.20978.x

CrossRef Full Text | Google Scholar

Fermi, E. (1949). On the Origin of the Cosmic Radiation. Phys. Rev. 75, 1169–1174. doi:10.1103/PhysRev.75.1169

CrossRef Full Text | Google Scholar

Gottlieb, S., Ketcheson, D., and Shu, C.-W. (2011). Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations. World Scientific. doi:10.1142/7498

CrossRef Full Text

Grimm, H. J., Gilfanov, M., and Sunyaev, R. (2002). The Milky Way in X-Rays for an outside Observer. Log(N)-Log(S) and Luminosity Function of X-Ray Binaries from RXTE/ASM Data. Astron. Astrophysics 391, 923–944. doi:10.1051/0004-6361:20020826

CrossRef Full Text | Google Scholar

Hairer, E., Nørsett, S., and Wanner, G. (2008). Solving Ordinary Differential Equations I: Nonstiff Problems. in Springer Series in Computational Mathematics. Springer Berlin Heidelberg).

Kelner, S. R., Aharonian, F. A., and Bugayov, V. V. (2006). Energy Spectra of Gamma Rays, Electrons, and Neutrinos Produced at Proton-Proton Interactions in the Very High Energy Regime. Phys. Rev. D 74, 034018. doi:10.1103/PhysRevD.74.034018

CrossRef Full Text | Google Scholar

Kissmann, R., Reitberger, K., Reimer, O., Reimer, A., and Grimaldo, E. (2016). Colliding-wind Binaries with Strong Magnetic Fields. Astrophysical J. 831, 121. doi:10.3847/0004-637X/831/2/121

CrossRef Full Text | Google Scholar

Kowal, G., de Gouveia Dal Pino, E. M., and Lazarian, A. (2011). Magnetohydrodynamic Simulations of Reconnection and Particle Acceleration: Three-Dimensional Effects. Astrophysical J. 735, 102. doi:10.1088/0004-637X/735/2/102

CrossRef Full Text | Google Scholar

Kowal, G., de Gouveia Dal Pino, E. M., and Lazarian, A. (2012). Particle Acceleration in Turbulence and Weakly Stochastic Reconnection. Phys. Rev. Lett. 108, 241102. doi:10.1103/PhysRevLett.108.241102

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyubarsky, Y., and Liverts, M. (2008). Particle Acceleration in the Driven Relativistic Reconnection. Astrophysical J. 682, 1436–1442. doi:10.1086/589640

CrossRef Full Text | Google Scholar

Mignone, A. (2007). A Simple and Accurate Riemann Solver for Isothermal MHD. J. Comput. Phys. 225, 1427–1441. doi:10.1016/j.jcp.2007.01.033

CrossRef Full Text | Google Scholar

Pittard, J. M., Vila, G. S., and Romero, G. E. (2020). Colliding-wind Binary Systems: Diffusive Shock Acceleration and Non-thermal Emission. Monthly Notices R. Astronomical Soc. 495, 2205–2221. doi:10.1093/mnras/staa1099

CrossRef Full Text | Google Scholar

Reimer, A., Pohl, M., and Reimer, O. (2006). Nonthermal High-Energy Emission from Colliding Winds of Massive Stars. Astrophysical J. 644, 1118–1144. doi:10.1086/503598

CrossRef Full Text | Google Scholar

Reitberger, K., Kissmann, R., Reimer, A., Reimer, O., and Dubus, G. (2014b). High-energy Particle Transport in Three-Dimensional Hydrodynamic Models of Colliding-Wind Binaries. Astrophysical J. 782, 96. doi:10.1088/0004-637X/782/2/96

CrossRef Full Text | Google Scholar

Reitberger, K., Kissmann, R., Reimer, A., and Reimer, O. (2014a). Simulating Three-Dimensional Nonthermal High-Energy Photon Emission in Colliding-Wind Binaries. Astrophysical J. 789, 87. doi:10.1088/0004-637X/789/1/87

CrossRef Full Text | Google Scholar

van Leer, B. (1979). Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov’s Method. J. Comput. Phys. 32, 101–136. doi:10.1016/0021-9991(79)90145-1

CrossRef Full Text | Google Scholar