Combined Magnetohydrodynamic- Monte Carlo Simulations of Proton Acceleration in Colliding Wind Binaries

TL;DR

This paper combines magnetohydrodynamic and Monte Carlo simulations to study proton acceleration in colliding-wind binary systems, aiming to improve gamma-ray flux predictions by modeling shock acceleration processes.

Contribution

It introduces a self-consistent approach integrating MHD and Monte Carlo simulations to analyze proton acceleration at both shocks in CWBs, considering magnetic field effects.

Findings

Protons are efficiently accelerated at both shocks in CWBs.

Magnetic field topology significantly influences particle acceleration.

Results can refine gamma-ray emission predictions for CWBs.

Abstract

The interaction between the strong winds of the stars in colliding-wind binary (CWB) systems produces two shock fronts, delimiting the wind collision region (WCR). There, particles are expected to be accelerated mainly via diffusive shock acceleration, and to produce $\gamma$-rays, in processes involving relativistic electrons and/or protons. We investigate the injection and the acceleration of protons in typical CWB systems by means of Monte Carlo simulations, with a test-particle approach. We use magnetohydrodynamic simulations to determine the background conditions in the wind collision region. This allows us to consider particle acceleration at both shocks, on either side of the WCR, with a self-consistently determined large-scale magnetic field, which has an impact on the shape of the WCR, and the topology of which plays an important role in particle acceleration at collisionless shocks. Such studies may contribute to improve $\gamma$-ray flux predictions for CWB systems.

Figures

Tables

Table 1. Table 1: Stellar and wind parameters of a typical colliding-wind binary system, as in Kissmann  et al. [ 23 ] . M ∗ subscript 𝑀 M_{*} is the stellar mass, R ∗ subscript 𝑅 R_{*} the stellar radius, T ∗ subscript 𝑇 T_{*} the effective temperature, L ∗ subscript 𝐿 L_{*} the luminosity, M ˙ ˙ 𝑀 \dot{M} the mass loss rate, v ∞ subscript 𝑣 v_{\infty} the terminal velocity of the wind, and B ∗ subscript 𝐵 B_{*} the surface magnetic field.

Star	$M_{\ast }$	$R_{\ast }$	$T_{\ast }$	$L_{\ast }$	$\dot{M}$	$v_{\infty }$	$B_{\ast }$
	[M☉]	[R☉]	[K]	[L☉]	[M☉ yr-1]	[km s-1]	[G]
B	$30$	$20$	$23000$	${10}^5$	${10}^{-6}$	4000	100
WR	$30$	$10$	$40000$	$2.3\times {10}^5$	${10}^{-5}$	4000	100

Table 2. Table 2: Spectral indices and injection efficiencies of protons injected at z = 40 ​ R ☉ 𝑧 40 subscript 𝑅 ☉ z=40\ R_{\astrosun} , z = − 420 ​ R ☉ 𝑧 420 subscript 𝑅 ☉ z=-420\ R_{\astrosun} , and z = 420 ​ R ☉ 𝑧 420 subscript 𝑅 ☉ z=420\ R_{\astrosun} . The errors refer to the fit of Equation 1 to the data.

z	Side of WCR	Spectral index	Injection efficiency
$\left[R_{☉}\right]$		$\sigma$	$\epsilon$
40	B	$1.86\pm 0.02$	$\approx 26\%$
	WR	$1.96\pm 0.02$	$\approx 16\%$
-420	B	$2.19\pm 0.01$	$\approx 9\%$
	WR	$2.18\pm 0.04$	$\approx 8\%$
420	B	$1.93\pm 0.02$	$\approx 16\%$
	WR	$2.08\pm 0.04$	$\approx 13\%$

Equations

$$J(E)=vn(E)\propto {\left[E\left(E+2m_p{c}^2\right)\right]}^{-\frac{\sigma }{2}},$$

$$\epsilon =\frac{n_{NT}}{n_{TOT}},$$

Full text

We investigate the injection and the acceleration of protons in typical CWB systems by means of Monte Carlo simulations, with a test-particle approach. We use magnetohydrodynamic simulations to determine the background conditions in the wind collision region. This allows us to consider particle acceleration at both shocks, on either side of the WCR, with a self-consistently determined large-scale magnetic field, which has an impact on the shape of the WCR, and the topology of which plays an important role in particle acceleration at collisionless shocks. Such studies may contribute to improve $\gamma$-ray flux predictions for CWB systems.

1 INTRODUCTION

Collisionless shock fronts are known to accelerate particles, primarily via first-order Fermi acceleration, resulting in a power law energy spectrum of the non-thermal population (e.g. [1, 2] and references therein).

Colliding-wind binaries (CWBs) - binary systems of massive, hot stars with strong stellar winds - are expected to be sites where particle acceleration occurs. Shocks form at the wind collision region (WCR), thus creating a suitable environment for accelerating particles by means of diffusive shock acceleration (DSA). Besides the detected radio synchrotron radiation from many such systems [3, 4, 5], processes such as inverse Compton (IC) scattering, relativistic Bremsstrahlung, and decay of neutral pions produced in hadronic interactions are expected to provide sufficiently high fluxes of GeV and TeV $\gamma$-rays to allow their detection by instruments such as Fermi-LAT, HESS, MAGIC and VERITAS (e.g. [6, 7, 8]). In contrast to such predictions, $\gamma$-ray emission from CWBs has not been observed up to now, with the exception of $\eta$ Carinae [9, 10, 11], and the still debated WR 11 [12].

In previous studies, the non-thermal photon emission of typical CWB systems was computed using the spectral energy distributions of high-energy particles obtained by combining three-dimensional hydrodynamic numerical simulations with the solution of the transport equation for protons and electrons [13, 14]. In these works, particle acceleration is provided by a very simplified term in the transport equation, which is the analytical result obtained when considering a population of suprathermal particles at a shock of infinite extent. The injection efficiency of thermal particles into the acceleration process, i.e. how many of the particles from the thermal distribution are able to recross the shock from downstream to upstream and be injected into DSA, is necessarily prescribed “by hand”.

For studying the microphysics of shock formation and injection efficiency, (full) particle in cell (PIC), and hybrid PIC simulations are suited best: particles move in the simulation box and are scattered by self-consistently generated magnetic turbulences (e.g. [15, 16] and references therein). Unfortunately, these approaches are computationally very demanding, and both box sizes and time intervals that can be investigated are limited.

Particle acceleration has also been studied by means of Monte Carlo simulations (e.g. [17, 18, 19]). Here, particles move undisturbed on a predefined background until a scattering occurs. The scattering process is necessarily modelled, but injection efficiencies can be obtained, depending on the prescribed scattering laws. Furthermore, this simplification significantly reduces the computational load, and allows for much larger regions and time intervals (and in turn energy ranges) to be considered.

For this reason, we chose the Monte Carlo method for our purpose of simulating particle acceleration in extended regions, and over wide energy and time intervals. In particular, we employ a test-particle approach, resembling methods used for studying particle injection at oblique shocks (e.g. [19, 20, 21]).

In order to obtain more realistic background conditions, we use magnetohydrodynamic (MHD) simulations of a typical CWB system.

In the following, we will first describe the numerical method employed in this work. We will then show the details of the system studied, as well as the results of the simulations. The last section is devoted to the conclusions.

2 NUMERICAL METHOD

In this work, we investigate particle acceleration in the wind collision region of CWBs. We employ the results of MHD simulations of a typical CWB system, carried out using the cronos code [22, 23], to characterize the background (plasma flow velocity, magnetic field, electric field, temperature, and density) that influence the particle motion in the simulation box. The background is initialized at the beginning of each simulation and, since we are applying a test-particle approach, it is not changed at runtime by the accelerated particles.

In the MHD simulations, the shock front is about three cells wide for numerical reasons only. Moreover, if the shock front is not “vertical” in the simulation box, the cell boundary is not parallel to the actual shock front (see Fig. 1). This can lead to unrealistic results, especially if the gyroradius of the particle is smaller than the cell size, and the particle gyrates many times along the shock surface when encountering the shock. After illustrating the Monte Carlo method itself, we will describe how we deal with these issues.

2.1 Monte Carlo method

Protons are injected into the selected cell, with an isotropic Maxwell-Boltzmann distribution of velocities in the frame comoving with the local plasma flow, with the temperature of the local background. Each particle is followed in the frame where the shock front is stationary, and moves on the background electromagnetic field, driven by the Lorentz force, until: (i) a scattering occurs, (ii) it reaches the boundary of a cell, or (iii) it is removed from the system. The time between scatterings is exponentially distributed, with mean value $t_{c}=\eta\tilde{r}_{g}/v$, where $\eta$ is a proportionality factor, $v$ is the speed of the particle, and $\tilde{r}_{g}=p/qB$ [21]. For $\eta=1$ this corresponds to Bohm diffusion. We assume the scattering to be elastic in the local plasma frame. The new $\mu_{p}=\cos{\theta_{p}}$, where $\theta_{p}$ is the pitch angle in the plasma frame, is a random number between 0 and 1, drawn from a uniform distribution. This models large fluctuations in the magnetic field [21].

When a particle crosses the boundary of a cell, the background changes, and its trajectory is computed accordingly. Particles are removed from the system either after they were scattered $10^{6}$ times (usually the case for particles which were not accelerated at the shock front and are being advected away by the plasma) or when they reach the outer boundary of the entire simulation box.

In order to improve statistics at higher energies, the technique of particle splitting is employed. A statistical weight is assigned to each particle at the beginning of the simulation. When the energy of a particle increases by a factor of 10, its statistical weight is halved, and the particle is splitted into two instantaneously identical ones, which will follow different paths from that moment onwards. The spectrum of the flux through a certain surface is recorded by adding to the appropriate energy bin the statistical weight of the particle crossing it. We verified our code by reproducing the results obtained in Ref. [19].

The background in our simulation is divided into two regions: one includes the shock-front cells, identified as being part of the shock front, the other one comprises all the other cells. When particles are in the non-shock-front cells, the background is just the one resulting from the MHD simulations. When, on the other hand, a particle enters the area marked as being part of the shock front, a different setup is used, with two bigger cells -upstream and downstream-, the background of which is obtained as described in the following section. The position of the particle as it enters this new regime is recorded, so that, when leaving the super-cell regime, the new position in the normal background is found by adding the displacement in the appropriate reference system to the recorded coordinates (see Fig. 2 (b)).

2.2 Background

The background is obtained from MHD simulations. At the shock fronts, the shock transition layer is about 3 cells wide, and the cell surface is in general not aligned with the shock front. We therefore divided the computational domain into two regions: “shock front” and “normal” cells. In the shock front domain, a system of “superimposed” cells (super-cells) is used. For technical reasons, these have dimensions $(2\Delta x)^{3}$ upstream, and $3\Delta x\times(2\Delta x)^{2}$ downstream. In order to initialize the super-cells, one first needs to locate the position of the shock front. This is done by setting a threshold for the temperature: at the shock fronts the plasma temperature abruptly rises from $\sim$ 104K to $\sim 10^{7}$-$10^{8}$K. The temperature gradient was found to be a good tracer of the shock front position by Reitberger et al. [13].

Every upstream super-cell is associated with a downstream one (and vice versa). The background of the upstream super-cells is obtained by averaging the fields of the neighbouring upstream shock-front cells. Following Reitberger et al. [13], the background of the associated downstream super-cell is obtained by computing a weighted average of the fields within a distance $\Delta x$ ($\equiv$ cell size) from the point 3$\Delta x$ downstream of the upstream shock-front cell, in the direction normal to the shock front (see Fig. 2 (a)). The vector fields (flow velocity, magnetic field, electric field) are then rotated, so that shock normal and normal to the surface between upstream and downstream super-cells are parallel (reference frame $\bf{S^{\prime}}$ in Fig. 2 (b)).

3 RESULTS

3.1 Parameters of the simulations

In the following we consider a B star and Wolf-Rayet star binary system, with the parameters listed in Table 1.

The stellar separation is $R=1440\mbox{ R}_{\astrosun}$. The region used in the Monte Carlo simulations consists of ($151\times 81\times 151$) cubic cells of dimension (3.9$\mbox{ R}_{\astrosun}$)3. We stress that no analytical prescriptions are necessary concerning the large-scale magnetic field at the WCR, since the magnetic field was evolved dynamically in the MHD background simulations.

Particles are injected upstream of the shock fronts, at different positions along the WCR, on the $x$-$z$ plane, at $z=40\ R_{\astrosun}$, $z=-420\ R_{\astrosun}$, and $z=420\ R_{\astrosun}$. The fluxes are recorded when particles cross the shock fronts. Here, we show the results of simulations carried out with $\eta=1$ (highly turbulent medium) [21]. The impact of varying $\eta$ on the results is subject to future studies.

3.2 Spectral indices and injection efficiencies

The analytical result for the density of particles accelerated via DSA, at a shock with compression ratio $r$, yields the well known dependence $n(p)\propto p^{-\sigma}$, where $p$ is the momentum, $n$ is the differential particle density, and $\sigma=(r+2)/(r-1)$ is the spectral index. The differential current of protons in terms of the kinetic energy $E$ is:

[TABLE]

where $v$ is the speed of the proton, $m_{p}$ is its rest mass, and $c$ is the speed of light.

The spectra resulting from the simulations can be seen in Fig. 3. The magnetic field on the WR-side of the WCR is weaker, which causes a difference of up to two orders of magnitude in the maximal energy reached by particles injected on the WR-side ($E^{WR}_{max}\sim 10^{11}$ eV) and those injected on the B-side ($E^{B}_{max}\sim 10^{12}-10^{13}$ eV). Moreover, the compression ratio is in general higher on the B-side, which results in harder spectra.

The Monte Carlo simulations also allow to estimate injection efficiencies for the system, once the scattering law (i.e. scattering operator and mean free path dependence on different parameters) has been chosen, and under the assumption that the shock front is seen as a sharp transition by the protons. If we define the injection efficiency $\varepsilon$ as:

[TABLE]

where $n_{NT}$ and $n_{TOT}$ are the particle density in the non-thermal tail and in the total particle distribution, respectively, we obtain $8\%\leq\varepsilon_{{}_{WR}}\leq 16\%$ for the WR shock, and $9\%\leq\varepsilon_{{}_{B}}\leq 26\%$ for the B shock. The variation in the injection efficiency seems to depend not only on the side of the WCR where the particles are injected, but also on the distance of the injection position from its apex, as can be seen in Table 2. This is probably due to both different plasma flow velocities and shock obliquities. The spectra obtained when injecting particles at $z=420\ R_{\astrosun}$ are harder than those of particles injected at $z=-420\ R_{\astrosun}$. This is ascribable to different compression ratios at the two injection positions. Although an association of this feature with the asymmetry of the plasma flow is tempting, we note that the compression ratio changes not continuously along the shock front, and further simulations are needed in order to verify if spectra are systematically harder above the plane at $z=0$. We note that efficiencies are rather high, therefore including feedback of the accelerated particles on the plasma would yield more realistic results. Nevertheless, our simulations indicate that the change in injection efficiency along the WCR may indeed be relevant when modelling non-thermal emission from CWB systems.

4 CONCLUSIONS

In this work, we investigated the acceleration of protons in a typical colliding-wind binary system. Our Monte Carlo test-particle simulations employ the results of magnetohydrodynamic simulations, which determine the background on which protons move.

We found a difference in both the spectral index and the highest energy reached by the accelerated particles, depending on the considered side of the WCR, due to different compression ratios and strengths of the magnetic field. Moreover, we found that these two characteristics of the spectra of non-thermal protons also vary moving away from the apex of the WCR. A similar remark can be done concerning the injection efficiencies, which seem to be in general higher on the B-side of the WCR, and to decrease further away from its apex, along the shock fronts. Our results indicate that a variation of the injection efficiency probably needs to be taken into account in models aiming at predicting $\gamma$-ray fluxes produced by non-thermal particles accelerated in CWB systems.

5 ACKNOWLEDGMENTS

E. G. and A. R. acknowledge financial support from the Austrian Science Fund (FWF), project P 24926-N27.