# Colliding-Wind Binaries as a Source of TeV Cosmic Rays

- Escola de Artes, Ciências e Humanidades, Universidade de São Paulo, São Paulo, Brazil

In addition to gamma-ray binaries which contain a compact object, high-energy and very high–energy gamma rays have also been detected from colliding-wind binaries. The collision of the winds produces two strong shock fronts, one for each wind, both surrounding a shock region of compressed and heated plasma, where particles are accelerated to very high energies. Magnetic field is also amplified in the shocked region on which the acceleration of particles greatly depends. In this work, we performed full three-dimensional magnetohydrodynamic simulations of colliding winds coupled to a code that evolves the kinematics of passive charged test particles subject to the plasma fluctuations. After the run of a large ensemble of test particles with initial thermal distributions, we show that such shocks produce a nonthermal population (nearly 1% of total particles) of few tens of GeVs up to few TeVs, depending on the initial magnetization level of the stellar winds. We were able to determine the loci of fastest acceleration, in the range of MeV/s to GeV/s, to be related to the turbulent plasma with amplified magnetic field of the shock. These results show that colliding-wind binaries are indeed able to produce a significant population of high-energy particles, in relatively short timescales, compared to the dynamical and diffusion timescales.

## 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^{4}–10^{5} 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 ^{2} 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 *Γ* is the adiabatic index set to 1.1. In our model, we set

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

where ρ, *p*,

The MHD equations were solved numerically using a high-order shock-capturing adaptive refinement Godunov-type code AMUN^{1} 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^{5} (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^{6} K.

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

where *q* and *m* are particle charge and mass, respectively, *c* is the speed of light, and

We assumed the particle velocity *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 *B*_{*} = 10^{−5} to 10^{0} 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^{5} 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^{2} 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

## 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^{−7} c.u. (white, in the middle of the presented color bar). However,

**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 ^{–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**. 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

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

### 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

As an example, we provide a movie demonstrating the evolution of protons during the time of integration for the model with *B*_{*} = 10^{−1} 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 *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**. 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 *h* being the effective cell size *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^{3} s and 10^{4} s, and eventually reaching nearly the speed of light between *t* = 10^{5} s and 10^{6} s (see Figure. 3c). At this stage, that is, after *t* = 10^{5} 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^{4} s, a significant fraction of protons can increase their kinetic energies, defined as

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

**FIGURE 4**. Left: Distribution of the particle largest kinetic energy for integration done with different values of

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):

with

where

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 ^{−3}, 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 *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

**FIGURE 5**. Left: Particle kinetic energy distribution at moments

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

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

### 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

**FIGURE 6**. Distribution of the maximum acceleration rate

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

**FIGURE 7**. Probability distribution function of the maximum acceleration rate ^{−5} to 10^{0} 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

## 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

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,

• The highest available energies achieved by the accelerated protons depend on the strength of magnetic field. For

• 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

• The regions of the highest acceleration rates

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

## Author’s Note

Visualization of the protons’ motion in the model with

## 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

^{1}AMUN code is open source and freely available from https://gitlab.com/gkowal/amun-code.

^{2}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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Keywords: stars: Wolf–Rayet, stellar winds, gamma-ray sources, particle astrophysics, cosmic rays, turbulence, magnetohydrodynamics, methods: numerical

Citation: Kowal G and Falceta-Gonçalves DA (2021) Colliding-Wind Binaries as a Source of TeV Cosmic Rays. *Front. Astron. Space Sci.* 8:667805. doi: 10.3389/fspas.2021.667805

Received: 14 February 2021; Accepted: 28 April 2021;

Published: 28 May 2021.

Edited by:

Luca Comisso, Columbia University, United StatesReviewed by:

Alberto Martin Gago, Pontifical Catholic University of Peru, PeruZhenbin Wu, University of Illinois at Chicago, United States

Copyright © 2021 Kowal and Falceta-Gonçalves. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Grzegorz Kowal, grzegorz.kowal@usp.br

^{†}These authors have contributed equally to this work and share first authorship