Proton acceleration in colliding stellar wind binaries

TL;DR

This paper models proton acceleration in colliding stellar wind binaries using Monte Carlo simulations and MHD data, revealing magnetic field amplification effects that influence gamma-ray emissions.

Contribution

It introduces a combined Monte Carlo and MHD approach to study proton acceleration and magnetic field amplification in CWBs, considering non-linear shock modifications.

Findings

Magnetic field amplification significantly reduces non-thermal proton flux.

Proton acceleration affects shock profiles and local magnetic fields.

Gamma-ray emissions are mainly from neutral pion decay, not inverse Compton.

Abstract

The interaction between the strong winds in stellar 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 (DSA). We investigate the injection and the acceleration of protons in typical CWB systems by means of Monte Carlo simulations, with both a test-particle approach and a non-linear method modelling a shock locally modified by the backreaction of the accelerated protons. We use magnetohydrodynamic simulations to determine the background plasma in the WCR and its vicinity. This allows us to consider particle acceleration at both shocks, on either side of the WCR, with a realistic large-scale magnetic field. We highlight the possible effects of particle acceleration on the local shock profiles at the WCR. We include the effect of magnetic field amplification due to resonant streaming instability (RSI), and compare results without and with the backreaction of the accelerated protons. In the latter case we find a lower flux of the non-thermal proton population, and a considerable magnetic field amplification. This would significantly increase the synchrotron losses of relativistic electrons accelerated in CWB systems, lowering the maximal energy they can reach and strongly reducing the inverse Compton fluxes. As a result, $\gamma$-rays from CWBs would be predominantly due to the decay of neutral pions produced in nucleon-nucleon collisions. This might provide a way to explain why, in the vast majority of cases, CWB systems have not been identified as $\gamma$-ray sources, while they emit synchrotron radiation.

Figures

Tables

Table 1. Table 1: Comparison of the solutions of this work, with the semi-analytical results of Caprioli et al. ( 2009 ) .

${\theta }_B$	$r_{\text{sub}}$	$r_{\text{tot}}$	$S_{\text{sub}}$	$S_{\text{tot}}$	$p_{\text{max}}$	$T_2$	$RSI$
[∘]					[${10}^6$ GeV/c]	[${10}^6$ K]
0	3.58	112.1	3.43	108.7	0.24	0.88	No
0	3.84	9.22	3.79	9.12	1.17	126.5	Yes
0	3.5	118			0.25	0.78	No
0	3.8	9.8			1	110	Yes
30	3.6	71			0.25	1.4	No
30	3.7	10.7			1	65	Yes
60	3.6	37			0.25	1.6	No
60	3.6	12.5			1	13	Yes
0	2.4	17			1	5.1	Yes
30	2.7	17			1	7.6	Yes
60	2.7	17			1	2.4	Yes

Table 2. Table 2: Stellar and wind parameters.

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 3. Table 3: Modified shock parameters at different positions along the WCR.

Position	($y,z$)	$r_{\text{sub}}$	$r_{\text{tot}}$	$T_2$	$M_{A0x}$	$P_{w1}$
	[$R_{☉}$]			[${10}^6$ K]
B1	(0, 20)	2.7	9.3	10	76	$1.8\times {10}^{-2}$
B2	(0, -20)	2.7	9.3	6.4	62	$1.7\times {10}^{-2}$
B3	(0, 420)	2.7	6.8	41	30	$2.2\times {10}^{-2}$
B4	(0, -420)	2.7	7.0	44	33	$2.3\times {10}^{-2}$
B5	(500, 420)	2.7	8.3	15	43	$2.9\times {10}^{-2}$
B6	(-500, -420)	2.7	8.4	15	45	$2.8\times {10}^{-2}$
W1	(0, 20)	2.6	80	0.32	$3.7\times {10}^4$	$1.1\times {10}^{-3}$
W2	(0, -20)	2.5	78	0.30	$3.0\times {10}^4$	$1.4\times {10}^{-3}$
W3	(0, 420)	2.7	33	0.80	$1.7\times {10}^3$	$6.1\times {10}^{-3}$
W4	(0, -420)	2.6	32	0.74	$1.5\times {10}^3$	$5.2\times {10}^{-3}$
W5	(500, 420)	2.5	34	0.51	$1.9\times {10}^3$	$4.6\times {10}^{-3}$
W6	(-500, -420)	2.7	33	0.53	$1.6\times {10}^3$	$6.9\times {10}^{-3}$

Equations

$$F_w=\frac{1}{2}\rho (\delta u{)}^2{u}_x+\frac{1}{{\mu }_0}\{(\mathbf{B}\times \mathbf{\delta u}+\mathbf{\delta B}\times \mathbf{u})\times \mathbf{\delta B}\}\cdot \hat{\mathbf{n}},$$

$$a_3{r}^3+a_2{r}^2+a_{1}r+a_0=0,$$

$$a_3$$

$$a_2$$

$$a_1$$

$$a_0$$

$$\aligned n&=\sum_{i}\left|\frac{u_{0}}{v_{x,i}}\right|w\\ \Phi=\sum_{i}v_{x,i}\left|\frac{u_{0}}{v_{x,i}}\right|w\ .$$

$$f(p_k)=\frac{1}{4\pi p_k^2\Delta p_k}n(p_k),$$

$$n(p_k)=\sum _{p_i\in \Delta p_k}\left|\frac{u_0}{v_{x,i}}\right|w.$$

$$u_x(x)\frac{\partial f(x,p)}{\partial x}=\frac{\partial }{\partial x}\left[D(x,p)\frac{\partial f(x,p)}{\partial x}\right]+\frac{\text{d}{u}_x(x)}{\text{d}x}\frac{p}{3}\frac{\partial f(x,p)}{\partial p}+Q(x,p),$$

$$Q(x,p)=\frac{\eta {\rho }_1{u}_1}{4\pi m_p{p}_{\text{inj}}^2}\delta (p-p_{\text{inj}})\delta (x),$$

$$\aligned D(x,p)&=D_{\parallel}(x,p)\cos^{2}\theta_{B}(x)+D_{\perp}(x,p)\sin^{2}\theta_{B}(x)\ ,\\ D_{\parallel}(x,p)=\xi\frac{p}{qB(x)}\frac{v}{3}\ ,\\ D_{\perp}(x,p)=\frac{\xi}{(1+\xi)^{2}}\frac{p}{qB(x)}\frac{v}{3}\ ,$$

$$f(x,p)=f_1(p)\exp \left\{-\frac{q(p)}{3}{\int }_{-\infty }^0\text{d}{x}^{\prime }\frac{u_x(x^{\prime })}{D(x^{\prime },p)}\right\},$$

$$f_{1}(p)=\frac{\eta n_{0}}{4\pi p_{\textrm{inj}}^{3}}\frac{3\ r_{\textrm{tot}}}{r_{\textrm{tot}}U_{px}(p)-1}\exp\left[-\int_{p_{\textrm{inj}}}^{p}\frac{dp^{\prime}}{p^{\prime}}\frac{3r_{\textrm{tot}}U_{px}(p^{\prime})}{r_{\textrm{tot}}U_{px}(p^{\prime})-1}\right]\$$

$$q(p)=-d\log f_1(p)/d\log p.$$

$${\mathbf{u}}_{\mathbf{p}}(p)={\mathbf{u}}_1-\frac{1}{f_1(p)}{\int }_{-\infty }^0\text{d}x\frac{\text{d}\mathbf{u}(x)}{\text{d}x}f(x,p).$$

$${\mathbf{B}}_{\mathbf{p}}(p)={\mathbf{B}}_1-\frac{1}{f_1(p)}{\int }_{-\infty }^0\text{d}x\frac{\text{d}\mathbf{B}(x)}{\text{d}x}f(x,p).$$

$$1+P_{g0}+P_{B0}=U_x(x)+P_g(x)+P_w(x)+P_B(x)+P_c(x).$$

$$P_g(x)=\frac{U_x(x{)}^{-\gamma }}{\gamma M_{0x}^2},$$

$$P_B(x)=\frac{B_z^2(x)}{2{\mu }_0{\rho }_0{u}_{0x}^2},$$

$$P_w(x)=\frac{U_x(x{)}^{-\frac{3}{2}}}{4M_{A0x}}\left[(1-U_x^2(x))\cos {\theta }_{B0}\right].$$

$$P_c(x)=\frac{4\pi }{3{\rho }_0{u}_0^2}{\int }_{p_{\text{inj}}}^{p_{\text{max}}}dp{p}^{3}v(p)f(x,p).$$

$$K=P_{c1}^{\text{eqref}eq:momFluxConsNorm}/P_{c1}^{\text{eqref}eq:PcAlternative},$$

$$f_1^{\ast }(p)=K{f}_1(p),$$

$$\frac{T_2}{T_1}=\frac{1}{r}\frac{p_{g2}}{p_{g1}}=\frac{1}{r[(\gamma +1)-(\gamma -1)r]}\left\{(\gamma +1)r-(\gamma -1)+(\gamma -1)\frac{M_{A1x}^4(r-1{)}^3}{(M_{A1x}^2-r{\cos }^2{\theta }_{B1}{)}^2}\frac{(p_{B1}+p_{w1})}{p_{g1}}\right\}.$$

$$\eta =\frac{4}{3\sqrt{\pi }}(r_{\text{sub}}-1){\psi }^3{e}^{-{\psi }^2}.$$

$$\delta B_p(p)=\delta B_1-\frac{1}{f_1(p)}{\int }_{-\infty }^0\text{d}x\frac{\text{d}\delta B(x)}{\text{d}x}f(x,p),$$

$$\frac{\partial {\mathcal{F}}_w(k,x)}{\partial x}=u_x(x)\frac{\partial {\mathcal{P}}_w(k,x)}{\partial x}+\sigma (k,x){\mathcal{P}}_w(k,x).$$

$$\sigma (k,x)=\frac{{\pi }^2{m}_p^2{\Omega }_0^2{v}_A}{B_0^2}2\pi \iint d\mu dp{p}^2(1-{\mu }^2)v^2\frac{{\hat{\mathbf{n}}}_{\mathbf{B}}\cdot \nabla f(x,p)}{{\nu }_{-}}\delta (kp|\mu |-m_p{\Omega }_0),$$

$$2U_x(x)\frac{\text{d}{P}_w(x)}{\text{d}x}=V_{Ax}(x)\frac{\text{d}{P}_c(x)}{\text{d}x}-3P_w(x)\frac{\text{d}{U}_x(x)}{\text{d}x},$$

$$P_w(x)=U_x(x{)}^{-\frac{3}{2}}\left[\frac{(1-U_x^2(x))\cos {\theta }_{B0}}{4M_{A0x}}\right].$$

Full text

binaries: general — magnetohydrodynamics — methods: numerical — shock waves — acceleration of particles

††journal: ApJ

1 Introduction

Collisionless shocks are known to be sites where charged particles can be accelerated to relativistic energies. A fraction of the thermal particles of the plasma can be scattered by magnetic fluctuations and, after multiple crossings of the shock, gain energy by means of the first-order Fermi acceleration, in this context also known as diffusive shock acceleration (DSA) (e.g. Drury (1983); Jones & Ellison (1991)).

An intriguing class of astrophysical objects, where particles are expected to be accelerated via DSA, are colliding-wind binaries (CWBs). These binary systems consist of hot, massive stars ejecting supersonic stellar winds, which eventually collide and form a (bow-shaped) wind-collision region (WCR), delimited by two shock fronts. Indeed, many CWBs have been identified as particle accelerators, mainly owing to the detection of a non-thermal radio component in their observed spectra (see e.g. De Becker & Raucq (2013) and references therein). Up to now, only two such systems have been identified as $\gamma$-ray sources, namely $\eta$ Carinae (Reitberger et al., 2012; Werner et al., 2013; Reitberger et al., 2015) and $\gamma^{2}$ Velorum (also known as WR 11, Pshirkov (2016)). This is somewhat unexpected: several models (analytical and semi-analytical) predicted the production of $\gamma$-rays with fluxes above the detection threshold of Fermi-LAT, HESS, MAGIC and VERITAS, mainly via inverse Compton scattering of electrons in the radiation field of the stars or due to the decay of neutral pions produced in hadronic interactions (e.g. Eichler & Usov (1993); Benaglia & Romero (2003); Reimer et al. (2006)). These models have been improved recently by means of hydrodynamic (HD) and magnetohydrodynamic (MHD) simulations, combined with the solution of the transport equations for electrons and protons (Reitberger et al., 2014a, b, 2017). Such models allow to study these systems and possible effects of their geometry on the emission of $\gamma$-rays in a more detailed and realistic manner. A limitation of the method is the need to set an “injection parameter”, which determines the fraction of thermal particles which can enter the acceleration process. Aiming at sidestepping this limitation, Grimaldo et al. (2017) combined MHD simulations of an archetypal CWB system, and Monte Carlo (MC) simulations of shock acceleration (similar to the one developed by Ellison et al. (1995)). They found a high acceleration efficiency at the shocks of CWBs, suggesting that the backreaction of the accelerated particles on the shock structure must be taken into account.

A considerable amount of studies was conducted, focussing on the non-linear effects of particle acceleration at collisionless shocks, and several methods have been developed, mainly in the context of supernova remnants. For example, the Monte Carlo approach conserving momenta and energy fluxes was developed, amongst others, by Ellison & Eichler (1984); Vladimirov (2009); Bykov et al. (2014). Semi-analytical methods (Amato & Blasi, 2005; Caprioli et al., 2009), and time-dependent DSA simulations (e.g. Kang et al. (2012)) have also been employed. A comparison between different methods for modelling shocks with non-linear DSA can be found in Caprioli et al. (2010). The most realistic approach for studying the injection mechanism is certainly the particle-in-cell (PIC) method (e.g. Caprioli & Spitkovsky (2014a)), which simulates the shock dynamics from first principles. However, it requires high computational costs. We therefore employ the MC technique to determine the fraction of accelerated particles with respect to the thermal plasma density. We use this injection parameter in the semi-analytical method developed by Caprioli et al. (2009), which includes the effects of the non-linear resonant instability. The equations of Caprioli et al. (2009) are generalized for the case of oblique shocks. This is necessary for considering the variety of shock obliquities along the WCR.

In the next two sections, we will introduce the most important equations describing the numerical method. In Section 4 we will compare the results for strictly parallel shocks to those for oblique shocks, with parameters typical of SNRs (used as a reference). In Section 5 we will apply the non-linear MHD-MC method to an archetypal CWB system, considering different positions along both shocks delimiting the WCR and comparing test-particle to non-linear results. The last section is devoted to the conclusions.

2 Theoretical background

In order to ensure a better comprehension of the problem, we review the most important equations and definitions. Different methods employed for modelling non-linear shocks show some common features concerning the modification of the velocity profile: a shock precursor develops, where the inflowing plasma is progressively slowed down by the pressure of the non-thermal protons, up to the position of the MHD shock, usually called “subshock”. The magnetic turbulence generated by the backstreaming charged particles during the acceleration process can have an effect on the shock structure (Caprioli et al., 2009). Therefore, we will use the MHD equations including the turbulence pressure. For the sake of simplicity, we consider only Alfvén waves generated by the resonant streaming instability. We use here the indices 0, 1 and 2 for quantities far upstream, directly upstream and downstream of the subshock, respectively. Let us suppose that the only spatial dependence of the background fields (flow velocity $\bm{u}$, magnetic field $\bm{B}$, density $\rho$, and temperature $T$) is in the $x$-direction. Following Scholer & Belcher (1971) and Decker (1988), we write the MHD conservation laws in the presence of Alfén waves as: {subequations} {align} [ρu_x]_1^2&=0  ,

[ρ\bmu u_x + (p_g + B22μ0 + P_w ) ^\bmn - Bx\bmBμ0 ]_1^2=0  ,

[u_x { 12 ρu^2 + γ**γ-1 p_g + B2μ0 + F_w}-Bx(\bmB⋅\bmu)μ0 ]_1^2=0  ,

_1^2=0  ,

_1^2=0  .

Herein, $\hat{\bm{n}}$ is the unit vector normal to the shock, $\mu_{0}$ is the permeability constant, $\gamma$ is the adiabatic index, $p_{g}=nk_{B}T$ is the thermal pressure of the plasma, where $T$ is its temperature, $n$ is the particle density, and $k_{B}$ is the Boltzmann constant. A subscript $x$, $y$ or $z$ indicates the $x$-, $y$- or $z$- component of a vector, while the notation $[]^{2}_{1}$ denotes the difference between downstream and upstream quantities. The pressure term associated with the Alfvén waves is $P_{w}=(\delta B)^{2}/(2\mu_{0})$, where $\delta B$ is the magnetic field amplitude of the wave. The energy flux $F_{w}$ in Eq. \eqrefeq:MHDEFluxConservation includes the kinetic energy flux and the $x$-component of the Poynting vector associated to the wave:

[TABLE]

where $\bm{\delta u}$ is the velocity change of the plasma due to the Alfvén waves. By using the transmission and reflection coefficients of Alfvén waves incident on a shock given by McKenzie & Westphal (1969), Scholer & Belcher (1971) obtain a third-order equation for the compression ratio $r\equiv\rho_{2}/\rho_{1}=u_{1x}/u_{2x}$ equivalent to:

[TABLE]

with coefficients

[TABLE]

Herein, $\theta_{B}$ is the angle between the shock normal and the magnetic field, $\lambda=(\delta B_{1}/B_{1})^{2}$, $M_{A1x}=u_{1x}/v_{A}$, with the Alfvén speed $v_{A}=B/\sqrt{\mu_{0}\rho}$, and $\beta=p_{g}2\mu_{0}/B^{2}$. In the limit $\lambda=0$ Eq. \eqrefeq:thirdr reduces to Eqs. (11) of Decker (1988).111This is true after correcting for a typo in Decker (1988): the term $\cos^{2}\delta_{1}$ should always be multiplied by $M_{A1}^{2}$ in their equations (11a)-(11c). Using Eqs. \eqrefeq:MHDConservation, one can find the expressions relating the variations of the magnetic field, the flow velocity and the density along $x$, assuming that the magnetic field lies in the $x$-$z$ plane: {subequations} {align} u_x(x)&=ρ0u0xρ(x) ,

u_y(x)=u_0y ,

u_z(x)=u_0z+( Bz(x)-B0zμ0 ) B0xρ0u0x  ,

{subequations} {align} B_x(x)&=B_0x ,

B_z(x)=( MA0x2-cos2θB0Ux(x) MA0x2- cos2θB0 ) B_0z  .

Here and in the following sections, velocities and pressures indicated with capital letters are normalized by $u_{0x}$ and $\rho_{0}u_{0x}^{2}$, respectively. The electric field is $\bm{E}=-\bm{u}\times\bm{B}$, and it is entirely due to the motion of the plasma in the magnetic field. In all of our set-ups, $\bm{u}$ is (almost) perfectly parallel to $\bm{B}$, and the electric field is therefore small.

For strong turbulence, $\delta B/B>1$, the definition of a “background magnetic field” becomes questionable. In such cases, we consider $\bm{B}(x)$ just as the field determining the direction of propagation of the Alfvén waves.

3 Numerical Methods

Our method combines MHD simulations of the wind plasma, obtained with the Cronos code (Kissmann et al., 2018, 2016), Monte Carlo simulations of shock acceleration, using the technique developed by e.g. Ellison et al. (1995), and a semi-analytical method for computing the local non-linearly modified shock profile (Amato & Blasi, 2005). These components can be run independently from each other. As far as the MHD code is concerned, we refer the reader to the above-cited references. Below we provide details concerning the Monte Carlo part, the semi-analytical part, and how the different components are combined.

3.1 Monte Carlo simulations

The Monte Carlo method is similar to that of Ellison et al. (1995). We let each particle move using the Bulirsch-Stoer algorithm (Press et al., 1992). The particles are acted upon by the Lorentz force given by the background electromagnetic fields. After a time $t_{c}$, exponentially distributed with a mean value $\bar{t}_{c}=\xi r_{g}/v$, a scattering occurs (elastic in the frame of the local plasma flow). Here, $r_{g}=p/(qB)$ is the gyroradius, where $p$ is the magnitude of the momentum of a particle, $q$ is its charge, $v$ is its speed, $B$ is the local magnetic field strength, and $\xi$ is a proportionality factor relating gyroradius and mean free path ($\lambda_{\textrm{mfp}}=\xi r_{g}$). In the following, we set $\xi=1$, corresponding to Bohm diffusion. This choice is supported by PIC simulations (Caprioli & Spitkovsky, 2014b). The new direction of the momentum vector is randomly determined at each scattering, mimicking strong magnetic turbulence.

The protons of the background plasma are assumed to have a Maxwell-Boltzmann distribution in the local plasma frame, and are injected accordingly, close to the shock front, following the prescription of Vladimirov (2009). The densities $n$ and fluxes $\Phi$ at the shock fronts are given by:

[TABLE]

Herein, $u_{0}$ is the flow speed, $w=n_{0}/N_{p}$, where $n_{0}$ is the particle density and $N_{p}$ is the number of particles injected, and $v_{x,i}$ is the $x$-component of the velocity of the particle crossing the measurement surface (i.e. the shock). The quantities $u_{0}$ and $w$ refer to the point where the particles are injected at the beginning of the simulation. The index $i$ runs over all the crossing events of all the simulated particles. As discussed by Vladimirov (2009) (see Eq. (3.9) of that work), if $\Delta p_{k}$ is the width of the $k$-th momentum bin centred at momentum $p_{k}$, the distribution function for particles passing a surface of interest is:

[TABLE]

with

[TABLE]

Here, the index $i$ runs over all the crossing events of the particles having a momentum within the $k$-th momentum bin.

In order to use the MHD results as the background for the Monte Carlo simulations, we first locate the position of the shocks by setting a threshold for the gradient of the temperature, which abruptly rises there from $\sim 10^{4}$ K to $10^{7}-10^{8}$ K. A system of superimposed cells (super-cells) is then initialized, with the purpose of (i) having a sharp shock jump for thermal particles, and (ii) avoiding artefacts in the acceleration of particles due to misalignment of the shock surface and the boundary between upstream and downstream cells. In fact, the simulated protons would not “see” a sharp shock, because the size of the MHD cells is much larger than the mean free path of the thermal particles, and the transition from upstream to downstream is about three cells wide.222In this work, the MHD cells are cubes of edge length $\Delta x=3.9\mbox{ R}_{\sun}$. A more detailed description of the procedure can be found in Appendix A.

We inject the particles close to the shock in one selected upstream super-cell and let them move and scatter from then onwards. We stop the simulation of a particle when (i) it reaches a distance $x_{D}=10D/u_{2}$ downstream of the shock, where $D$ is the diffusion coefficient given by Eq. \eqrefeq:diffOblique, when still in the initial super-cell system, or (ii) it leaves the whole simulated box, or (iii) the number of scatterings experienced by the particle reaches a pre-set value, much greater than the expected mean number of scatterings needed to reach the highest possible energy in the system. Assuming an infinitely extended shock downstream, as it is effectively for low-energy (thermal) particles, the choice $x_{D}=10D/u_{2}$ corresponds to stopping the simulation for a particle when its probability to return to the shock is $\lesssim e^{-10}$(Ostrowski & Schlickeiser, 1993).

In the next section, we will describe the semi-analytical method used for the determination of the local non-linear modifications of the shock.

3.2 Semi-analytical non-linear calculations

In order to obtain a shock profile which conserves energy and momentum fluxes, we adapt the procedure developed by Caprioli et al. (2009) to the case of oblique shocks. It consists of an iterative method solving the diffusion-advection equation for the accelerated particles at the shock. The background conditions are determined by the velocity $\bm u$, the density $\rho$, the temperature $T$, and the magnetic field $\bm B$. We assume that all the considered quantities change locally only in the $x$-direction, and that $v_{A}\ll u$. We can then write:

[TABLE]

where $f(x,p)$ is the isotropic part of the distribution function of the accelerated particles. The source term, which accounts for the injection of particles in the acceleration process, is given by:

[TABLE]

where $m_{p}$ is the proton mass, the $\delta$ is the Dirac delta distribution, $p_{\textrm{inj}}$ is the “injection momentum”, and $\eta$ is the “injection efficiency” (see below for more details concerning these last two terms). The diffusion coefficient is (e.g. Jones & Ellison (1991)):

[TABLE]

$D_{\parallel}$ and $D_{\perp}$ being the diffusion coefficients parallel and perpendicular to the magnetic field lines, respectively. The solution to Eq. \eqrefeq:diff-adv is given by (see, e.g., Amato & Blasi (2005)):

[TABLE]

where

[TABLE]

is the distribution function immediately upstream of the subshock, with

[TABLE]

Here, $r_{\textrm{tot}}=u_{0x}/u_{2x}$ is the total compression ratio of the shock, while the mean velocity of the scattering centres (of the plasma) “seen” by a particle of momentum $p$ is:

[TABLE]

Similarly, for the magnetic field we have:

[TABLE]

Accordingly, the “mean” electric field is: $\bm{E_{p}}=-\bm{u_{p}}(p)\times\bm{B_{p}}(p)$. The momentum flux conservation equation, normalized by the kinetic momentum flux, reads:

[TABLE]

Here, $P_{g}$ is the thermal pressure of the plasma, $P_{w}$ is the pressure associated with the Alfvén waves produced by the resonant streaming instability, $P_{B}$ is the background magnetic field pressure due to the $z$-component of the field, and $P_{c}$ is the pressure of the accelerated protons. Considering only adiabatic heating in the precursor, one has:

[TABLE]

where $M_{0x}^{2}=\rho_{0}u_{0x}^{2}/(\gamma p_{g0})$. The $z$-component of the magnetic field exerts a pressure:

[TABLE]

where the $B_{z}$ can be found using Eq. \eqrefeq:Bz(x). This pressure term is present only if the shock is not strictly parallel, and is usually negligible in the precursor. However, as we will see, it has an influence on the jump conditions at the subshock and on the total compression ratio in the case of efficient particle acceleration, especially in the absence of magnetic field amplification. The term $P_{w}$, the pressure due to Alfvén waves, is given by (see Appendix B for a brief derivation):

[TABLE]

In order to find the non-linearly modified shock profile, we proceed as illustrated by Caprioli et al. (2009), using the equations adapted for the case of oblique shocks. We first set the compression ratio at the subshock, $r_{\textrm{sub}}$, and therefore the total compression ratio $r_{\textrm{tot}}$ (see below for more details concerning this point). A scheme of the algorithm can be seen in Figure 1.

In the first iteration we set $u_{px}(p)=u_{x}(x)=u_{1x}$. Starting with a different value, e.g. $u_{0x}$ does not affect the final result. The magnetic field, the wave pressure and the diffusion coefficient are calculated according to Eqs. \eqrefeq:Bz(x), \eqrefeq:alpha(x) and \eqrefeq:diffOblique. We compute $f_{1}(p)$ using Eq. \eqrefeq:f1, which in turn allows to compute $P_{c1}$, i.e. $P_{c}(x)$ at the subshock, according to:

[TABLE]

At this point, we calculate the pressure of the accelerated particles also with Eq. \eqrefeq:momFluxConsNorm, and we find

[TABLE]

where the bracketed exponent indicates the equation used for the computation. This allows us to normalize $f_{1}(p)$ and obtain

[TABLE]

so that the momentum flux between far upstream and the subshock is conserved. This normalized distribution function is used to calculate $f(x,p)$ by means of Eq. \eqrefeq:fxp, and $P_{c}$ by means of Eq. \eqrefeq:PcAlternative. Finally, the velocity profile $U_{x}(x)$ can be obtained using Eq. \eqrefeq:momFluxConsNorm, which allows to find $\bm B(x)$ and $P_{w}(x)$, and in turn the diffusion coefficient $D(x,p)$ using Eqs. \eqrefeq:Bz(x), \eqrefeq:alpha(x) and \eqrefeq:diffOblique, respectively. A new iteration is then started by calculating $u_{px}(p)$ by means of Eq. \eqrefeq:up. In order to achieve faster convergence, we average the flow profile between iteration $n$ and $n-1$, before computing $\bm B(x)$, $D(x,p)$ and $u_{px}(p)$. A similar solution has also been used by Amato & Blasi (2005) (private communication). We stop the cycle when $K$ does not change more than a specified amount between consecutive iterations. At this point, we check the value of $K$: if it is within a $15\%$ tolerance interval around 1, the solution has been found, otherwise a new $r_{\textrm{sub}}$ is used and the procedure is repeated (see Section 4 for a discussion on the convergence criterion). The compression ratio is increased if $K>1$, while it is decreased if $K<1$.

Eq. \eqrefeq:thirdr allows us to compute the flow velocity and magnetic field downstream, when combined with Eqs. \eqrefeq:Bz(x) and \eqrefeq:uz(x), once the conditions at the subshock are known. The temperature downstream is given by the relation:

[TABLE]

This equation can be obtained with the same approach as used by Vainio & Schlickeiser (1999), except using the conservation laws and transmission and reflection coefficients of Alfvén waves for the case of oblique shocks given by Scholer & Belcher (1971).

In order to find $r_{\textrm{tot}}$, we numerically solve Eq. \eqrefeq:thirdr keeping $r=r_{\textrm{sub}}$ fixed and employing the relation $U_{1x}=r_{\textrm{sub}}/r_{\textrm{tot}}$. In this way, knowing the far upstream conditions and $r_{\textrm{sub}}$, we can determine the background directly upstream and downstream of the subshock. At this point, the only missing ingredient is the fraction $\eta$ of particles being injected into the acceleration process. Caprioli et al. (2009) use the formula:

[TABLE]

They consider a Maxwell-Boltzmann distribution with the temperature of the shocked plasma and assume that only the particles with momentum $p_{\textrm{inj}}\geq\psi p_{\textrm{th,2}}$, with $p_{\textrm{th,2}}=\sqrt{2mk_{B}T_{2}}$, can be injected into the Fermi acceleration. This solution aims at modelling the thickness of the shock, assuming that particles with gyroradii smaller than the shock thickness can only be advected away from the shock and will not contribute to the non-thermal tail of the particle distribution function. Values of $\psi\simeq 2-4$ are usually chosen. Blasi et al. (2005) showed that, for example, $\psi\approx 2$ corresponds to a shock thickness $\lambda_{sh}=r_{g}^{\textrm{th,2}}$, and $\psi\approx 3.25$ corresponds to $\lambda_{sh}=2r_{g}^{\textrm{th,2}}$, where $r_{g}^{\textrm{th,2}}$ is the gyroradius of a particle with momentum $p_{\textrm{th,2}}$.

In this work, we choose a more accurate way to determine $\eta$ by means of Monte Carlo simulations. At the beginning of the first cycle associated to the first guess for the compression ratio $r_{\textrm{sub}}$, we roughly estimate the “injection efficiency” as follows. We initialize the background for MC simulations with the parameters of the subshock (i.e. $\bm u_{1}$, $\bm B_{1}$, etc., upstream and $\bm u_{2}$, $\bm B_{2}$, etc., downstream). We then inject the particles upstream, close to the shock, and determine $\eta$ as $\eta=n_{\textrm{ret}}/n_{\textrm{tot}}$, where $n_{\textrm{ret}}$ is the number of particles recrossing the shock from downstream to upstream, and $n_{\textrm{tot}}$ is the total number of injected particles. Using the injection efficiency so obtained, we calculate a new $\psi$ which satisfies Eq. \eqrefeq:etaCapr, and a new $p_{\textrm{inj}}=\psi p_{\textrm{th,2}}$. We thus find the modified shock solution for the current $r_{\textrm{sub}}$ employing the semi-analytical method described above. Finally, in order to obtain a more accurate estimate for the injection efficiency, and in turn for the density of the non-thermal population, which is essential for obtaining an energy-conserving solution, we run Monte Carlo simulations letting particles reach momenta $p\approx 100p_{\textrm{inj}}$. The particle density obtained from the Monte Carlo runs is then compared, at the momentum $p=10p_{\textrm{inj}}$, to the particle density obtained from the semi-analytical calculations, in order to find a corrected $\eta$ and the respective $K$. The comparison is done at $10p_{\textrm{inj}}$ in order to avoid the low-energy part of the spectrum after the thermal peak, which shows some oscillations, and the cut-off at the end of the distribution. A more accurate estimation would require to simulate the spectra up to higher energies, since there can be a difference in the slopes, up to momenta of $p\approx m_{p}c$, of the non-thermal distributions obtained with the Monte Carlo technique as compared to the ones obtained with the semi-analytical method. This discrepancy, ascribed to different treatment of the transition between thermal and non-thermal particles, was also found in Caprioli et al. (2010). It has been shown in the same work that the spectra at high energies (i.e. above momenta $p\approx m_{p}c$) are in good agreement. For the purposes of this paper, the possible gain in accuracy in the determination of the normalization of the non-thermal distributions does not justify the related increase of the computation times.

For technical reasons, instead of subdividing the computational box upstream into many cells, we use only two cells (upstream and downstream). The shock modification is taken into account by using the momentum-dependent averaged quantities $\bm{u_{p}}(p)$, $\bm{B_{p}}(p)$: when a particle of momentum $p$ is in the upstream cell, we use the background fields $\bm{u_{p}}(p)$, $\bm{B_{p}}(p)$, given by Eqs. \eqrefeq:up and \eqrefeq:Bp, and the corresponding $\bm{E_{p}}(p)$. Caprioli et al. (2009) chose to use the amplified magnetic field for computing the diffusion coefficient, i.e. they apply Eq. \eqrefeq:diffOblique using $B(x)=\delta B(x)\equiv\sqrt{2\mu_{0}P_{w}(x)\rho_{0}u_{0}^{2}}$. In Figure 2 we show an example of how the magnetic field varies with the distance from the shock. Far upstream, the particle density is $n_{0}=0.5\times 10^{6}$ m*-3*, the magnitude of the magnetic field is $B_{0}=5\times 10^{-11}$ T, the temperature is $T_{0}=10^{4}$ K, and the plasma speed is $u_{0}=5.9\times 10^{6}$ m s*-1*. We show the case with shock obliquity $\theta_{B0}=30^{\circ}$ (see also Table 1). The change of $B(x)$ is due to the change in the $z$-component of the magnetic field (Eq. \eqrefeq:Bz(x)), while $\delta B(x)$ changes according to Eq. \eqrefeq:alpha(x). The diffusion coefficient is therefore clearly changing with the distance from the subshock (see Eq. \eqrefeq:diffOblique). In order to obtain the appropriate “mean” diffusion coefficient for a particle of momentum $p$, we employ an effective $\bm{B^{\textrm{eff}}_{p}}(p)$, oriented like the computed $\bm{B_{p}}(p)$, but with a magnitude $|\bm{B_{p}^{\textrm{eff}}}(p)|=\max{(|\bm{B_{p}}(p)|,\delta B_{p}(p))}$, when making use of the magnetic field amplification due to resonant streaming instability. Here,

[TABLE]

and $\bm{B_{p}}$ is given by Eq. \eqrefeq:Bp. In this way, the gyroradius of a particle in the Monte Carlo simulations is the same as the gyroradius in the expression for the diffusion coefficient (Eq. \eqrefeq:diffOblique).

Accordingly, we use the appropriate electric fields, as well as mean scattering times. Choosing such a diffusion coefficient appears to be reasonable, considering results from PIC simulations: Caprioli & Spitkovsky (2014b) found that, for strong shocks, the energetic particles experience Bohm diffusion in the magnetic field amplified predominantly by the non-resonant hybrid instability (Bell, 2004).

4 Validation and comparison

In this section we compare the results of Caprioli et al. (2009) to ours and we further show the results obtained by keeping the parameters of the plasma unchanged, while changing the orientation of the shock. Our treatment differs to some extent from the one in Caprioli et al. (2009). In fact, they account for the possibility that the scattering centres move in the plasma frame with the Alfvén speed calculated using the background magnetic field. Nevertheless, they do not find a strong discrepancy between the “effective compression ratio” and the MHD compression ratio. For the sake of simplicity, we therefore do not include this effect. We accept solutions where the discrepancy between the pressure calculated with Eq. \eqrefeq:PcAlternative and with Eq. \eqrefeq:momFluxConsNorm lies within a $15\%$ tolerance range. Due to the statistical nature of the MC simulations, some fluctuation in the injection efficiency is unavoidable, at every cycle with $r_{\textrm{sub}}$ fixed, and a tighter tolerance range would require very high statistics, which in turn would require unreasonable computation times. We found that there is a dramatic improvement concerning momentum and energy flux conservation when employing the non-linear approach, as compared to the test-particle set-up. A tighter constraint on the tolerance range would be only a minor correction. Moreover, the uncertainties in the microphysics of acceleration in shocks where the ions reach energies well above their rest mass, as well as the possible development of other instabilities (e.g. the non-resonant hybrid instability) would make it unlikely to improve the reliability of the results by strengthening the convergence criteria. In Table 1 we summarize the results of our calculations, with shock obliquities of $0^{\circ}$, $30^{\circ}$ and $60^{\circ}$ and compare to two examples of Caprioli et al. (2009). For all cases, the particle density is $n_{0}=0.5\times 10^{6}$ m*-3*, the magnitude of the magnetic field is $B_{0}=5\times 10^{-11}$ T, the temperature is $T_{0}=10^{4}$ K, the plasma speed upstream, in the shock frame is $u_{0}=5.9\times 10^{6}$ m s*-1*, and $\psi=3.7$, as was used in the reference paper.

The results for the case of a strictly parallel shock are in good agreement, despite the “loose” convergence criteria and the slightly different approach, as mentioned above.

In Figure 3 (a) we show the distribution function of the non-thermal population of protons obtained by means of the semi-analytical approach alone, for the cases without the RSI effect, summarized in Table 1 (rows 3, 5, 7). A trend towards lower densities of accelerated particles is apparent when the shock obliquity increases. This happens despite the (slight) increase of the compression ratio at the subshock, and is an effect of the presence of a non-zero $z$-component of the magnetic field, and the associated pressure: Eqs. \eqrefeq:B(x) and \eqrefeq:PB imply that the decrease in the $x$-component of the plasma velocity results in an increase in $B_{z}$, and in turn in $P_{B}$. At the subshock, Eq. \eqrefeq:thirdr must hold, resulting in a lower $r_{\textrm{tot}}$: the additional magnetic pressure reduces the overall compressibility of the plasma (see Table 1). The downstream temperature is also affected by $P_{B}$: as expected from Eq. \eqrefeq:T2, it increases when passing from shock obliquities of $0^{\circ}$ to $60^{\circ}$.  The second plot of Figure 3 shows the cases with $P_{w}\neq 0$. Recall that $f_{1}$ is multiplied by $[p/(m_{p}c)]^{4}$. The particle density at the injection momentum $p_{\textrm{inj}}$ increases with increasing obliquity, due to the increasing compression ratio $r_{\textrm{tot}}$ and the decreasing downstream temperature. Nevertheless, the decrease of $p_{\textrm{inj}}$ (caused by the lower $T_{2}$), combined with the form of the distribution function $f\propto p^{-q}$, is such that the curves at momenta $p\geq p_{\textrm{inj,0}}$ ($p_{\textrm{inj,0}}$ being the injection momentum for $\theta_{B0}=0^{\circ}$) are lower for higher obliquities, resulting in an overall shift of the curves downwards for more oblique cases (Figure 3 (b)).

In Figure 4 (a) we show the results of the approach combining Monte Carlo simulations and semi-analytical calculations of the shock modifications. The curves of the semi-analytical calculations, obtained after determining $\eta$ with the combined approach as described above, and the distributions resulting from the full Monte Carlo simulations with the calculated modified background, are in good agreement. The discrepancy, higher for the oblique shock with $\theta_{B0}=60^{\circ}$, is lower than a factor of 4. The difference between the spectra of the semi-analytical calculation and the Monte Carlo particle distributions is mainly in the nonrelativistic regime. As already mentioned, a similar feature has also been found in Caprioli et al. (2010). In Figure 4 (b) we compare the semi-analytical non-thermal spectra obtained by fixing $\psi=3.7$ with those employing Monte Carlo simulations for the determination of the injection efficiency. These latter have lower $p_{\textrm{inj}}$ and much higher injection efficiencies. However, the discrepancies between the two approaches rapidly decrease with increasing momentum, and close to $p_{\textrm{max}}$ are less than a factor of $\approx 2$. The spectral indices are also different, due to the stronger shock modification in the combined approach (see also Table 1). The similarity at high energies of the solutions with the different injection efficiency determination suggests that the results and discussion presented in Section 5 are quite robust.

5 Results and discussion

In this section, we use the combined method described above to investigate the acceleration of protons in colliding-wind binaries, including a comparison of the test-particle results with the ones that include the back-reaction of the non-thermal protons. The unmodified background is given by a snap-shot of the simulation of an archetypal CWB (see Figure 5), the parameters of which are given in Table 2. The stars in the MHD simulation do not rotate, and there is no orbital motion.

The stellar separation is $R=1440\mbox{ R}_{\sun}$. The region used in the Monte Carlo simulations consists of $(x\times y\times z)=(201\times 101\times 101)$ cubic cells of dimension (3.9$\mbox{ R}_{\sun}$)3. The system is the same as Model A2 of Kissmann et al. (2016).

After initializing the background as described in Appendix A, we select 12 super-cells where we inject thermal protons. For the test-particle simulations, we do not modify the background any more and we just let the particles move and scatter in the simulated region. For the simulations including the back-reaction of the accelerated protons, we first find a modified background for the super-cells where the protons are injected, as outlined in Section 3. Once a solution is found, we start a simulation embedding the modified super-cells into the same simulated region of the test-particle case. A self-consistent determination of the maximal energies to which the protons are accelerated in the CWB when back-reaction is taken into account would require the modification of the MHD background of the entire WCR, together with the regions upstream of the shocks being modified by the pressure of the non-thermal protons. This is not possible with the approach presented here. Grimaldo et al. (2017) have shown, with test-particle simulations for the same system, that due to the differences in the magnetic field strength on the two sides of the WCR (see Figure 5 (a)), the maximal energies of the accelerated protons can differ by an order of magnitude or more. In fact, on the WR-side, where the magnetic field is weaker, the protons have larger gyroradii and the distribution functions have cutoffs at lower energies ($\approx 10^{2}\ m_{p}c$), as compared to the B-side ($\gtrsim 10^{3}\ m_{p}c$). Therefore, based on results of the test-particle simulations, the maximal momentum of the protons is set to $10^{2}\ m_{p}c$ on the WR-side and $10^{3}\ m_{p}c$ on the B-side of the WCR. Such an approximation does not take into account that some of the most energetic particles are also accelerated in cells different from the initial one, where the shock is modified. This, combined with the effect described in Section 3.2, can lead to an overestimation (more often) or to an underestimation of the proton density at high energies from the semi-analytical calculation, as compared to the final Monte Carlo spectra. However, this does not considerably affect the results of this work.

In Figure 6 (a) we show the distribution functions $f_{1}(p)[p/(m_{p}c)]^{4}$, obtained by injecting the protons close to the apex of the WCR, on the B-side, for both the test-particle and the feedback approach. Note the difference between the test-particle spectra obtained for injection at B1 and at B2. Whereas the spectral indices are very similar at small energies, they become notably different in the relativistic regime, where the B1 spectrum hardens, while the B2 spectrum softens. This is ascribable to different conditions downstream of the respective shocks. The particles which are energetic enough to return to the shock from farther downstream at B1 effectively “see” a higher compression ratio, caused by the slowdown of the plasma approaching the contact discontinuity. The effect seen in the test-particle simulations is entirely due to the interaction of the stellar winds, and the geometry of the WCR. In the investigated scenario the downstream flow is slower than it would be for a shock of infinite extent downstream, calculated with the shock-jump conditions. Therefore, it is easier for the particles, to cross the shock again from downstream to upstream in our setup than it would be in a one-dimensional shock structure. We do not see, however, any appreciable effect of hardening of the spectra due to scattering of the particles between the upstream “colliding shock flows” on the two sides of the WCR, as modelled by Bykov et al. (2013). This is likely due to different background conditions at, and between, the shocks. Bykov et al. (2013) considered a completely symmetric set-up, and obtained a hard spectrum for the distribution function, namely $f_{1}(p)\propto(U_{p}(p)p)^{-3}$. For energetic particles with $\lambda_{\textrm{mfp}}>L_{\textrm{WCR}}$, where $L_{\textrm{WCR}}$ is the width of the WCR, it should be possible to see this effect even without considering any modification due to the back-reaction of the accelerated protons upstream of the shocks (i.e. with $U_{p}(p)=1$). In the system presented here, however, the conditions on the two sides of the contact discontinuity are not equal. Particularly important is the difference in the magnetic field, which is much weaker on the WR-side. As a consequence, the particles accelerated at the B-side shock which manage to cross the contact discontinuity get much larger mean free paths and can therefore more easily escape the system from the WR-side of the WCR. Moreover, the protons accelerated at the WR-side shock do not reach sufficiently high energies that would allow crossing the contact discontinuity and reaching the B-side shock before they are advected out of the simulation box. Therefore, they are not accelerated by scatterings between the two converging flows upstream of the shocks. This situation might change for different parameters of the system, such that the magnetic field strength on both sides of the contact discontinuity is similar. Other factors might also play an important role, e.g. the width and the curvature of the WCR, and the distances between the stars and the shocks. Further studies are required, in order to better understand such effects.

As opposed to the B1 spectrum, the B2 spectrum does not harden, but instead it softens in the relativistic regime. This happens because the particles downstream of the shock are more efficiently advected away from the shock once they reach the equatorial plasma flow of the B star which enters the WCR and flows downwards with relatively high velocities (for a discussion concerning the stellar winds and WCR structure in this and other systems, see Kissmann et al. (2016)).

When comparing the test-particle spectra with those obtained with shock modification shown in Figure 6 (a), the effect of the back-reaction of the accelerated protons is clearly visible. The density of the non-thermal protons is reduced by up to more than three orders of magnitudes. The shock modification is stronger on the WR-side, as can be seen in Figure 6 (b) and in Table 3: the “thermal peak” moves towards lower momenta, while the total compression ratio reaches much higher values on the WR-side. This is caused by the very high Alfvén Mach number on the WR-side. The Alfvén waves produced via the resonant streaming instability lower considerably the value of $r_{\textrm{tot}}$, but not as much as on the B-side of the WCR, where the pressure associated to the waves is up to about one order of magnitude larger (see Table 3). In Figures 7 and 8 we show the same as in Figure 6, but for different positions. In all cases the test-particle results strongly overestimate the acceleration efficiency, as was found in many studies of SNRs (e.g. Vladimirov (2009), Bykov et al. (2014), Amato & Blasi (2005), Kang et al. (2012)).

A self-consistent quantitative modelling of the $\gamma$-ray emission from real systems is beyond the scope of this paper. Instead we conjecture the impact of back-reaction on the estimated $\gamma$-ray fluxes from CWB systems. The lack of detection of $\gamma$-rays from CWBs is well known. For example, Werner et al. (2013) did not find any evidence for $\gamma$-ray emission from various systems (amongst which, WR 11, WR 140, WR 147 have been considered), despite the predictions of several models (e.g. Eichler & Usov (1993); Benaglia & Romero (2003); Reimer et al. (2006)). Besides $\eta$ Carinae, $\gamma^{2}$ Velorum (WR 11) has been recently detected as a $\gamma$-ray source. In the latter, the dominant $\gamma$-ray production channel is probably the decay of neutral pions produced in nucleon-nucleon collisions, because of synchrotron and inverse Compton losses in the magnetic and radiation fields from the stars at the WCR (Reitberger et al., 2017). In the model of Reimer et al. (2006), the dominant production mechanism for $\gamma$-rays in WR 140 and WR 147 was found to be inverse Compton scattering from relativistic electrons in the stellar radiation field. This picture might considerably change if particle acceleration is efficient and magnetic field amplification takes place, leading to a situation similar to what has been found by Reitberger et al. (2017). There, the inverse Compton losses in the strong radiation fields prevent the electrons from reaching sufficiently high energies for $\gamma$-ray production. The $\gamma$-ray emission is therefore likely of hadronic origin. In the same work it is further claimed that a stronger magnetic field and smaller inverse Compton losses would also yield $\gamma$-ray fluxes in agreement with observations.

As far as the normalization of the non-thermal tail of the particle distribution is concerned, an injection parameter is commonly used. In Reitberger et al. (2017), for example, the proton density at 1 MeV is set to $n(E=1\textrm{MeV})=\eta_{1\textrm{MeV}}n_{0}$, with the injection parameter $\eta_{1\textrm{MeV}}=10^{-3}$, and $n_{0}$ the proton density of the background plasma. In Figure 9 we show the particle density for a selection of injection positions obtained with feedback, together with the non-thermal tails of the test-particle case, shifted in order to match the injection parameter of Reitberger et al. (2017). Amongst B1-B6 and W1-W6, we chose the positions corresponding to the higher densities at high energies in the test-particle case, close to the apex and far away from it, for each side of the WCR. We see that the actual particle densities might be even lower than what was obtained with $\eta_{1\textrm{MeV}}=10^{-3}$.

Together with the spectral energy distributions of the particles, the $\gamma$-ray fluxes from a modified WCR are also expected to change, but not necessarily in the same manner: the increased density in the WCR (higher compression ratios) will also have an effect, since hadronic collisions producing neutral pions will be more frequent there. In the following, we roughly estimate the $\gamma$-ray production for the model studied here, in the local vicinity of the apex of the WCR, using as the target proton population the thermal particles downstream of the considered shock. At position B1, the ratio between the non-thermal proton density (Eq. \eqrefeq:defn1) in the model with the modified shock, $n_{1}^{\textrm{NL}}(60\textrm{ GeV})$, and the one in the test-particle approach, with $\eta_{1\textrm{MeV}}=10^{-3}$, $n_{1}^{\textrm{TP}}(60\textrm{ GeV})$, is $n_{1}^{\textrm{NL}}(60\textrm{ GeV})/n_{1}^{\textrm{TP}}(60\textrm{ GeV})\approx 0.02$. At an energy of $10$ GeV (corresponding to a proton energy of $E_{p}\approx E_{\gamma}/\kappa_{\pi^{0}}\approx 60$ GeV, with the inelasticity factor for pion production $\kappa_{\pi^{0}}\simeq 0.17$ (Aharonian & Atoyan, 2000)), the ratio between the $\gamma$-ray emissivities would be $q_{\gamma}^{\textrm{NL}}(10\textrm{ GeV})/q_{\gamma}^{\textrm{TP}}(10\textrm{ GeV})\approx 0.07$. At the WR-side of the WCR, at position W2, the ratios of the densities of the non-thermal population at $\approx 6$ GeV is $n_{1}^{\textrm{NL}}(6\textrm{ GeV})/n_{1}^{\textrm{TP}}(6\textrm{ GeV})\approx 0.06$, while the $\gamma$-ray emissivity in the case of a modified shock, due to the much larger compression ratio, would be higher than in the test-particle case: $q_{\gamma}^{\textrm{NL}}(1\textrm{ GeV})/q_{\gamma}^{\textrm{TP}}(1\textrm{ GeV})\approx 2$. Despite being a limiting case, since the non-linear compression ratio close to the apex on the WR-side reaches the highest values, this estimate highlights the non trivial modifications of the $\gamma$-ray fluxes due to non-linear modifications of the shocks of the WCR in colliding-wind binary systems.

Another intriguing observable effect of proton acceleration in CWBs, even if highly speculative, might be the change of the opening angle of the WCR-cone. Reitberger et al. (2017) found that modelling $\gamma^{2}$ Velorum with a “strong coupling” between the stellar winds, i.e. including the radiative braking of the wind of the O star due to the photon field of the WR star, yields an opening half-angle of the shock-cone ($\approx 72^{\circ}$), i.e. closer to the observed value of $\approx 85^{\circ}$, as compared to the case when radiative braking is neglected ($\approx 24^{\circ}$). High-energy protons escaping the WCR from the O-side towards the WR-side, might slow down even further the wind of the WR star, and therefore contribute to a further shock-cone opening. Indeed, in our simulation we observe that some of the protons with higher energies eventually reach the contact discontinuity at the interface between the B-side and the WR-side downstream regions. There, the magnetic field changes abruptly, being weaker on the WR-side, increasing the mean free path of the particles and enhancing their escape probability on that side.

The discussion above shows that a viable way to reconcile the theoretical predictions for $\gamma$-ray fluxes with the observations of CWBs in the radio and $\gamma$-ray bands, is to employ the widely accepted idea of the existence of a back-reaction of the accelerated particles at collisionless shocks.

6 Conclusions and outlook

In this work we have presented a combined approach, employing magnetohydrodynamic simulations, a semi-analytical method for obtaining non-linear solutions of modified collisionless shocks, and Monte Carlo simulations of proton acceleration. In order to use the methods usually applied to strictly parallel shocks to the broad variety of obliquities found along the wind-collision regions of colliding-wind binaries, we adapted the equations to the case where the shock normal and the magnetic field are not aligned. By applying our method to a model of a typical CWB, we showed that, similarly to what was found in studies of typical SNR shocks, the Monte Carlo test-particle results differ considerably from the non-linear solutions: the test-particle approach greatly overestimates the injection efficiencies, which are dramatically reduced when energy and momentum conservation at the shock is fulfilled. Remarkably, we found indications that the injection and acceleration efficiencies at the shocks of CWBs may be lower than what is often assumed in approaches based on the transport equation for the accelerated particles, neglecting their back-reaction on the shocks. In the test-particle approximation, the maximal energies reached by the protons are different, mainly depending on which side of the WCR they were injected. On the B-side, energies of up to almost $10$ TeV can be reached, while on the WR-side the cut-off is in the range of $10$-$100$ GeV. This is ascribed to the different strengths of the magnetic fields. We note, however, that this difference might be reduced if the shock modifications were globally taken into account, due to an stronger magnetic field amplification on the WR-side, which reduces the difference of the magnetic field strength. The total compression ratios differ systematically from the B-side to the WR-side, the latter being much higher, due to a smaller magnetic turbulence pressure at the subshocks on that side. Also based on the results of Reitberger et al. (2017), we formulated the hypothesis that magnetic field amplification due to the accelerated protons could increase synchrotron losses of electrons accelerated at the shocks delimiting the WCR. This would reduce the maximal energy reached by the relativistic electrons, preventing them from efficiently producing $\gamma$-rays via inverse Compton scattering in the stellar photon fields. This might help to explain why non-thermal synchrotron emission has been observed by many CWB systems, while so far there has been no detection of $\gamma$-rays from those sources, with only one exception so far, that is $\gamma^{2}$ Velorum ($\eta$ Carinae does not show any synchrotron emission in the radio domain, presumably due to synchrotron self-absorption.). Further and deeper investigations of this and other observable effects of non-linear shock modifications in colliding-wind binaries, including the application to observed CWB systems, will be subject of future studies.

E. G. and A. R. acknowledge financial support from the Austrian Science Fund (FWF), project P 24926-N27. E.G. further acknowledges the receipt of a “Doktoratsstipendium aus der Nachwuchsförderung” (2017/3/MIP-3) from the University of Innsbruck.

Appendix A Background treatment

Here we will give details of the approach applied in order to use the MHD results as the background for the particle acceleration simulations. We justify first the need for a modification of the “raw” outputs.

In the MHD modelling, the transition between upstream and downstream of a shock is about three cells wide. When simulating CWBs, this layer is much bigger than the mean free path of the thermal particles. As a consequence, the compression ratio “seen” by the protons injected into the MHD background would be much different from the actual one, which would in turn effectively prevent particles from being accelerated. The Monte Carlo simulations run under the assumption that the scattering centres of the thermal particles are comoving with the background plasma, with a discontinuity at the position of the shock (or the subshock, in the case of non-linear modifications). Therefore, we chose to set up a system of superimposed cells, of size $(2\Delta x)^{3}$ upstream, and $3\Delta x\times 2\Delta x\times 2\Delta x$ downstream, with a sharp jump between the fields of the two cells ($\Delta x$ is the size of a cell in the MHD simulation). In Figure 10 (a), we depict the procedure for the test-particle approach, described in the following.

First of all, we identify the position of the shock as mentioned in Section 3.1, as well as its orientation. We then mark as “upstream cells” those being immediately before the WCR, in the direction pointed by the plasma flow, and as “downstream cells” the first two cells following the shock. We associate each upstream cell to a couple of upstream-downstream super-cells. The background upstream is determined by averaging the fields of the respective MHD upstream cell and the neighbouring upstream-marked ones. The background downstream is a weighted average of the fields of the cells within a radius $\Delta x$ from the point $3\Delta x$ downstream of the upstream MHD cell (Reitberger et al., 2014a) . The super-cell fields are then rotated so that the shock surface an the cell boundary between the upstream super-cell and the downstream one have the same orientation. This is necessary, since that boundary is, for the simulated particles, the shock surface. In the non-linear case, we initialize the super-cell couple where the protons are injected with the fields obtained as described in Section 3.2.

During the simulation, when a particle enters a shock front cell (upstream or downstream), its position is saved and the simulation is performed with a two-cells setup, until the particle leaves the super-cells. At this point, the position of the particle in the whole simulation box is calculated, starting from the previously saved position and adding the (appropriately rotated) total displacement in the super-cells. The simulation is then continued in the “normal” MHD background until the end or until it enters again the shock front domain.

Appendix B Wave pressure

Here, we summarize the derivation of the pressure associated with the Alfvén waves in the case of an oblique shock. Following Caprioli et al. (2009), we start from the stationary equation for growth and transport of magnetic turbulence which, for an oblique shock, reads:

[TABLE]

Herein, $\mathcal{F}_{w}$ is the energy flux, $\mathcal{P}_{w}$ is the pressure, both per unit logarithmic bandwidth, $\sigma$ is the growth rate of the energy in the magnetic turbulence. The latter is given by (Skilling, 1975), and for the case of only backwards propagating waves, as assumed here and in many other works (e.g. McKenzie & Völk (1982), Kang & Jones (2007), Caprioli (2012)) it reads:

[TABLE]

where $\Omega_{0}=qB_{0}/m_{p}$ is the nonrelativistic gyrofrequency, $B_{0}$ is the (mean) background magnetic field, $\bm{\hat{n}_{B}}$ is the unit vector along $B_{0}$, $\nu_{-}=\pi\Omega_{0}\mathcal{P}_{w}/(4\gamma U_{B})$ is the collision frequency of particles against waves moving forward (backwards) in the frame comoving with the plasma. In the latter expression, $\gamma$ is the Lorentz factor of the particle, while $U_{B}=B_{0}^{2}/(2\mu_{0})$ is the magnetic energy density of the background field. Integrating Eq. \eqrefeq:turbulence over $k$, and normalizing by $\rho_{0}u_{0x}^{2}$, yields:

[TABLE]

where $V_{Ax}(x)=v_{A}(x)\cos\theta_{B}(x)/u_{0x}$. As discussed in Caprioli et al. (2009), we can neglect $P_{B}$, $P_{w}$ and $P_{g}$ in Eq. \eqrefeq:momFluxConsNorm in the precursor (but not at the subshock) with respect to the kinetic momentum flux and the pressure of the accelerated particles, if acceleration is efficient, and use $P_{c}(x)\simeq 1-U_{x}(x)$. The solution of Eq. \eqrefeq:diffPw, assuming no wave pressure far upstream is:

[TABLE]