3D Magnetohydrodynamic Models of Nonthermal Photon Emission in the Binary System γ² Velorum

Implementing stellar winds with a fixed terminal velocity does not reflect reality and is disfavored by observations showing that the winds in ${\gamma }^{2}$ Velorum hit each other long before reaching terminal velocity (Henley et al. 2005). A more sophisticated approximation of radiative wind acceleration becomes necessary. Its details of implementation are a crucial point of our model as they significantly alter the outcome. In previous work we have relied on the standard Castor–Abbott–Klein (CAK) approximation (Castor et al. 1975) in its modified form proposed by Pauldrach et al. (1986). So far, we have used its variant of "weak coupling." The current version of the code can also use the case of "strong coupling." In the following we motivate the existence of these two approaches and argue about our choice of method for the ${\gamma }^{2}$ Velorum system.

A commonly used expression for the line acceleration of a single star is

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{\mathrm{rad},i}^{l}={g}_{\mathrm{rad},i}^{l}\displaystyle \frac{{{\boldsymbol{r}}}_{i}}{{r}_{i}},\end{eqnarray} \tag{ 1 }$$

with the radial acceleration term

$$\begin{eqnarray}&&{g}_{\mathrm{rad},i}^{l}=\displaystyle \frac{{\sigma }_{e}}{c}\displaystyle \frac{{L}_{\star ,i}}{4\pi {r}^{2}}{{kt}}^{-\alpha }{I}_{\mathrm{FD}},\end{eqnarray} \tag{ 2 }$$

the optical depth parameter

$$\begin{eqnarray}&&t=({\sigma }_{e}\rho {v}_{{th}})\,/\left|\displaystyle \frac{\partial u}{\partial r}\right|,\end{eqnarray} \tag{ 3 }$$

the finite-disk correction

$$\begin{eqnarray}{I}_{\mathrm{FD}} & = & \displaystyle \frac{1}{1-\cos {\theta }_{\star ,i}^{2}}\displaystyle \frac{1}{\pi }{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{{\theta }_{\star ,i}}{\left(({{\boldsymbol{n}}}_{i}\cdot {\rm{\nabla }}({{\boldsymbol{n}}}_{i}\cdot {\boldsymbol{u}}))/\left|\displaystyle \frac{\partial u}{\partial r}\right|\right)}^{\alpha }\\ & & \times \cos {\theta }_{i}d{{\rm{\Omega }}}_{i},\end{eqnarray} \tag{ 4 }$$

and the angle

$$\begin{eqnarray}&&\sin {\theta }_{\star ,i}=\displaystyle \frac{{R}_{\star ,i}}{{r}_{i}}.\end{eqnarray} \tag{ 5 }$$

The radial acceleration component ${g}_{\mathrm{rad},i}^{l}$ is therefore a function of the distance to the star r_i, its bolometric luminosity ${L}_{\star ,i}$, the velocity of the wind ${\boldsymbol{u}}$, the stellar radius ${R}_{\star ,i}$, and the standard CAK parameters α and k,

$$\begin{eqnarray}&&{g}_{\mathrm{rad},i}^{l}=f({\boldsymbol{u}},{L}_{\star ,i},{R}_{\star ,i},{r}_{i},\alpha ,k).\end{eqnarray} \tag{ 6 }$$

If applied to a system of only one star, these equations can be readily applied. However, if we consider two stars, their interpretation is more difficult. The effect of the companion star radiation field on the primary star wind, and vice versa, will lead to radiative inhibition (the star B radiation field reduces the wind acceleration in the atmosphere of star A) and radiative braking (rapid deceleration of the star A wind—shortly before reaching the WCR—due to the increasing radiative force of star B). To address both effects correctly, one has to know how the radiation of star B acts on the wind of star A and vice versa. This, however, is not entirely clear because Equation (6) leaves room for two different interpretations.

The key question is whether the two CAK parameters, α and k, are related to the properties of the wind plasma or the stellar radiation field. Do they characterize the radiation of the star (method A) or the capability of the wind material to react to radiation (method B)? In mathematical terms,

$$\begin{eqnarray}\mathrm{method}\ {\rm{A}}\qquad {g}_{\mathrm{rad}}^{l} & = & f({\boldsymbol{u}},{L}_{\star ,1},{R}_{\star ,1},{r}_{1},{\alpha }_{1},{k}_{1})\\ & & +\,f({\boldsymbol{u}},{L}_{\star ,2},{R}_{\star ,2},{r}_{2},{\alpha }_{2},{k}_{2}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathrm{method}\ {\rm{B}}\\ \quad {g}_{\mathrm{rad}}^{l}=\left\{\begin{array}{ll}f({\boldsymbol{u}},{L}_{\star ,1},{R}_{\star ,1},{r}_{1},{\alpha }_{1},{k}_{1}) & \\ \quad +\,f({\boldsymbol{u}},{L}_{\star ,2},{R}_{\star ,2},{r}_{2},{\alpha }_{1},{k}_{1}) & {\rm{if}}\,{\rm{in}}\,{\rm{wind}}\,{\rm{1}}\\ f({\boldsymbol{u}},{L}_{\star ,1},{R}_{\star ,1},{r}_{1},{\alpha }_{2},{k}_{2}) & \\ \quad +\,f({\boldsymbol{u}},{L}_{\star ,2},{R}_{\star ,2},{r}_{2},{\alpha }_{2},{k}_{2}) & {\rm{if}}\,{\rm{in}}\,{\rm{wind}}\,{\rm{2}}\end{array}\right..\end{eqnarray} \tag{ 8 }$$

Both methods yield roughly the same result for CWB systems in which the two stars and the location of the WCR are far apart. However, the differences are significant for systems where this is not the case (e.g., η Carinae during periastron, and ${\gamma }^{2}$ Velorum during its entire orbit).

The CAK approximation has initially been devised to describe the line acceleration of a single star and thus is to be used with great caution in the case of CWB systems. While the physical interpretation of the parameter α is clear—it represents the ratio of optically thin and optically thick lines in the wind plasma (Cranmer & Owocki 1995)—, the meaning of the parameter k is rather elusive as it depends on the wind thermal temperature while determining the strength of the line acceleration component (Pittard 2009). In general, the two CAK parameters represent a fit to the total line driving force depending on the properties of the ions in the wind as well as on the radiation field of the star. Accordingly, methods A and B both constitute borderline cases of the problem, which, depending on the actual composition of a binary systems, should be applicable to a varying degree at the same time.

Past numerical simulations of CWB system have predominately used method A (e.g., Pittard 2009; Parkin et al. 2011; Madura et al. 2013; Reitberger et al. 2014b), while theoretical works on radiative inhibition and braking generally apply method B (e.g., Stevens & Pollock 1994; Owocki & Gayley 1995; Gayley et al. 1997). The latter (also known as "strong coupling") is slightly more complicated and costly to implement in 3D hydrodynamical simulations, as it requires the introduction of numerical means that allow to discriminate between the two wind components.

The use of method A (also known as "weak coupling") in studies on the η Carinae system (e.g., Parkin et al. 2011) is also understandably motivated: if method B were used for this system, the turbulent WCR close to the periastron passage would be a factor of ∼2 closer to the surface of the WR star, thus making it even more complicated to simulate than it already is. This problem, however, is specific only to the η Carinae system and its extreme values of bolometric luminosity and mass-loss rate of its primary component. For other CWB systems it is usually method B that leads to a physically more meaningful representation of the wind, and as we show for the case of ${\gamma }^{2}$ Velorum, to improved agreement with observations linked to the shape of the WCR.

Following our extension of CWB simulations from HD to MHD (as detailed in Kissmann et al. 2016), we no longer rely on any magnetic field approximation as in previous studies. In computing the synchrotron losses for the electrons accelerated at the shock, we previously approximated the magnetic field following Usov & Melrose (1992). Now, we use the three-dimensional magnetic field components as determined by the MHD simulations. The additional information of the local direction of the magnetic field makes it possible to consider the contributions of both stars, not just the star with the dominant field.

As the surface magnetic field of high-mass stars remains poorly constrained and as this study is not focused on the significant distortions that the presence of strong stellar dipole fields (${B}_{\mathrm{surface}}\gtrsim {10}^{-2}\,T$) have on the structure of the WCR, we choose a moderate field strength of 10⁻³ T for each star. Such a field does not produce a significant effect on the density and velocity distribution of the wind plasma. However, because of its strong influence via synchrotron losses, it still has a significant effect on the electron population. As we show in Section 3, electrons can reach higher energies along the magnetic equator because of the lower field strength and therefore lower synchrotron losses. Along the stellar dipole of the magnetic field, the losses are higher.

The surface magnetic field strength of the two stars are important free parameters in determining the shape of and conditions at the WCR, thus influencing not only the MHD properties of the system, but also the arising particle populations and nonthermal emission components. While a choice of ${B}_{\mathrm{surface}}={10}^{-3}T$ differs only insignificantly from the case of no surface magnetic field at all, higher surface magnetic field strengths lead to severe distortions in the WCR, which becomes narrower, more turbulent, develops stronger curvature, and also features a nose-like structure, as discussed in Kissmann et al. (2016). Accordingly, other free parameters such as the diffusion coefficient and the injection rate of high-energy protons (as discussed and identified below) have to be chosen differently for higher magnetic field strengths in order to find agreement between modeling and observations.

In our earlier studies (Reitberger et al. 2014a, 2014b) the computation of diffusive shock acceleration (DSA) of charged particles at the shock fronts was made via the energy-gain rate

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{DSA}}=\displaystyle \frac{{c}_{r}-1}{3{c}_{r}\kappa }{V}_{s}^{2}E,\end{eqnarray} \tag{ 9 }$$

where c_r is the compression ratio, κ the energy-independent diffusion coefficient, V_s the shock velocity, and E the energy of the particle. This method required the following steps: 1. identification of "acceleration cells" at the shock front and determination of the shock orientation therein, 2. computation of V_s as the velocity component of the wind perpendicular to the shock front, 3. computation of c_r via interpolation of the density structure perpendicular to the shock front, and 4. solving Equation (9) for every acceleration cell. Steps 1 to 4 become unnecessary when using our new method to describe particle acceleration. In the Appendix, we show that Equation (9) corresponds to the expression for the adiabatic cooling as used in the transport equation,

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{adiab}}=-\displaystyle \frac{E}{3}{\rm{\nabla }}\cdot {\boldsymbol{u}}.\end{eqnarray} \tag{ 10 }$$

At the shock front, where the divergence of the velocity is negative, this term describes the acceleration of the particles.

This new method greatly simplifies the code and at the same time allows for computation of particle acceleration for turbulent shock fronts where the exact location of the acceleration region and the compression ratio are difficult to define.

While maximum electron energies are governed by synchrotron and inverse-Compton losses, which usually suppress the acceleration at higher energies, the maximum energies of the protons are governed by the diffusion and convection of the particles. An often-used approximation is Bohm diffusion, which imposes a cutoff at maximum energies where the protons leave the shock fronts before being further accelerated. Such an approximation is no longer necessary if an energy-dependent diffusion coefficient is used (Kirk et al. 1998). In earlier simulations (Reimer et al. 2006; Reitberger et al. 2014a, 2014b) the effects of spatial diffusion were approximated via a fundamental loss time in the transport equation—similar to the treatment in a leaky-box model. Now, the transport equation includes a classical diffusion term $D(E){{\rm{\nabla }}}^{2}j$. The full equation (further motivated in Appendix 1) is

$$\begin{eqnarray}&&\displaystyle \frac{\partial j}{\partial t}-\mathop{D(E){{\rm{\nabla }}}^{2}j}\limits_{(1)}+\mathop{{\rm{\nabla }}\cdot ({\boldsymbol{uj}})}\limits_{(2)}-\mathop{\displaystyle \frac{\partial }{\partial E}\left[\left(\displaystyle \frac{E}{3}{\rm{\nabla }}\cdot {\boldsymbol{u}}+{\dot{E}}_{\mathrm{loss}}\right)j(E)\right]}\limits_{(3)}\\ &&\quad =\,\mathop{{Q}_{0}\delta (E-{E}_{0})}\limits_{(4)},\end{eqnarray} \tag{ 11 }$$

where j is the differential number density of particles at energy E at position ${\boldsymbol{r}}$, ${\boldsymbol{u}}$ is the velocity of the wind plasma, and Q₀ is the injection rate of particles at energy E₀. The term ${\dot{E}}_{\mathrm{loss}}$ includes losses by inverse-Compton, synchrotron, and bremsstrahlung emission, Coulomb cooling, as well as losses by nucleon–nucleon collisions (details in Reitberger et al. 2014b). While the spatial convection term (2) is solved along with the MHD equations via the treatment of high-energy particles as advected scalar fields, the diffusion term (1) and the energy-loss and energy-gain term (3) are solved in separate routines, both of which are capable of subcycling if a high-diffusion coefficient D(E) (or a very high energy-loss rate) should require it.

As term (1) in Equation (11) indicates, diffusion is treated isotropically. The modifications now allow considering an energy-dependent diffusion coefficient of the form $D(E)\,={D}_{0}{\left(\tfrac{E}{1\mathrm{MeV}}\right)}^{\delta }$. The exponent δ can be set to $\delta =0.3$ for a Kolmogorov-type spectrum (recommended for 3D MHD simulations) and to $\delta =0.5$ for a Kraichnan-type turbulence spectrum (Strong et al. 2007). We consistently use $\delta =0.3$ in this study (although it remains a free parameter in principle). The diffusion coefficient D₀ remains a free parameter.

In Section 3 we study the effects of changing the values of the parameters δ and D₀ on the resulting particle spectra. Spatial diffusion as specified above determines the cutoff in the proton spectrum in most cases.

The hot shocked gas of the WCR is expected and observed to be a strong thermal X-ray emitter. In our model, the significant energy loss of the shocked plasma due to its thermal emission is accounted for by the radiative cooling term Λ appearing in the energy equation (see Reitberger et al. 2014b). As in the earlier versions of our code, our cooling term is based on the work of Schure et al. (2009). However, our previous implementation of radiative cooling was only valid in the case of a wind that predominantly consists of fully ionized hydrogen (${\mu }_{i}=1,{\mu }_{e}=1$). While this is safe to assume for OB type stars, WR stars demand a different assumption. For the corrected cooling terms considering metallic wind composition, we find

$$\begin{eqnarray}&&{\rm{\Lambda }}\to \displaystyle \frac{1}{{\mu }_{e}}\displaystyle \frac{1}{{\mu }_{i}}{\rm{\Lambda }},\end{eqnarray} \tag{ 12 }$$

where ${\mu }_{i}$ and ${\mu }_{e}$ depend on whether the location for which the cooling term is computed is situated within the O wind or the WR wind. To discriminate between the two wind components, we make use of the same numerical means that were used in regard to the radiative wind acceleration. Assuming the WR star to consist of fully ionized pure He (${\mu }_{i}=4,{\mu }_{e}=2$), a disregard of this correction due to composition would overestimate Λ by a factor of 8. This leads to overefficient cooling, resulting in a very dense and turbulent WCR. When the above correction is applied, we find a less turbulent WCR.

Stellar and stellar wind parameters used in our model of ${\gamma }^{2}$ Velorum along with the corresponding references are given in Table 1. We used the most up-to-date and best-constrained values found in literature, mostly relying on the fundamental parameter determination by North et al. (2007). The CAK parameters α and k where fitted to obtain the desired wind parameters in a 1D simulation that is consequently adapted to 3D. Details of this procedure are found in Kissmann et al. (2016). For the WR star we decided to run simulations for both values of the mass-loss rate that have been determined either by polarimetric measurements ($\dot{M}=8\times {10}^{-5}\,{M}_{\odot }$ yr⁻¹) or via radio-continuum observations ($\dot{M}=3\times {10}^{-6}\,{M}_{\odot }$ yr⁻¹). The figures and values presented in this publication show the results determined with the latter, which we deem more likely due to reasons given below. Note that in Reitberger et al. (2017) we used the former mass-loss rate, which led to different results.

Table 1.  Stellar and Stellar Wind Parameters of ${\gamma }^{2}$ Velorum Used in this Study

Star	M*	R*	T*	L*	$\dot{M}$	${v}_{\infty }$	α	k	Bsurface
	(M⊙)	(R⊙)	(K)	(L⊙)	(M⊙ yr−1)	(km s−1)			T
O7.5	28.5 [1, 2]	17 [2]	35,000 [1]	2.8 × 105 [2]	1.78 × 10−7 [1]	2500 [1]	0.613	0.055	10−2
WC8	9 [2]	6 [2]	56,000 [2]	1.7 × 105 [2]	8 × 10−6 [2, 3]	1450 [2, 3]	0.526	0.747	10−2
					3 × 10−5 [1, 2]			1.498

Note.  Taken from [1] De Marco et al. (2000), [2] North et al. (2007), and [3] Schild et al. (2004).

Download table as:  ASCIITypeset image

The orbital parameters are shown in Table 2. Our choice reflects the best-constrained distance determination so far and reliable values for stellar separation and eccentricity. The inclination i of the orbital plane of the ${\gamma }^{2}$ Velorum binary system and its argument of apastron ${\omega }_{\mathrm{WR}}$ are well constrained. Schmutz et al. (1997) find $i=65^\circ \pm 8^\circ$ and ${\omega }_{\mathrm{WR}}=68^\circ \pm 4^\circ$. The exact definition of the angle ${\omega }_{\mathrm{WR}}$ becomes clear by the consideration that if the inclination were 90° (edge-on view), the WR star would eclipse the O star at the apastron passage if ${\omega }_{\mathrm{WR}}=0$, and ∼9 days after apastron if ${\omega }_{\mathrm{WR}}=68^\circ$. Thus, the meaning of the angle ${\omega }_{\mathrm{WR}}$ (as in Schmutz et al. 1997) is consistent with the definition of the angle Φ in Reitberger et al. (2014a).

Table 2.  Orbital Parameters Used in this Study

Distance d	Semimajor Axis a	Period P	Eccentricity e	Inclination i	${\omega }_{\mathrm{WR}}={\rm{\Phi }}$
(pc)	(${R}_{\odot }$)	(d)		(°)	(°)
336 [1]	258 [1]	78.53 [2]	0.334 [2]	65 [2]	68 [2]

Note.  Taken from [1] North et al. (2007) and [2] Schmutz et al. (1997).

Download table as:  ASCIITypeset image

All other parameters relevant for particle acceleration or nonthermal γ-ray emission are either the same as in Reitberger et al. (2014a, 2014b) or are discussed individually in the sections below.

In a short-period binary system like ${\gamma }^{2}$ Velorum, the shape and plasma conditions of the WCR strongly depend on the choice of weak or strong coupling. Figure 1 shows the possible components of the total line acceleration ${g}_{\mathrm{rad}}^{l}$ as in Equations (7) and (8) in the primary wind along the line of centers connecting the two stars. The acceleration acting on the primary wind by the primary star radiation field $f({\boldsymbol{u}},{L}_{\star ,1},{R}_{\star ,1},{r}_{1},{\alpha }_{1},{k}_{1})$ is shown in solid red, the acceleration acting on the primary wind by the secondary star radiation field is shown in dashed green for the case of weak coupling ($f({\boldsymbol{u}},{L}_{\star ,2},{R}_{\star ,2},{r}_{2},{\alpha }_{2},{k}_{2})$, method A) and in dotted blue for the case of strong coupling ($f({\boldsymbol{u}},{L}_{\star ,2},{R}_{\star ,2},{r}_{2},{\alpha }_{1},{k}_{1})$, method B).

Zoom In Zoom Out Reset image size

The computation of ${g}_{\mathrm{rad}}^{l}$ for all three cases is based on the primary wind initial velocity profile neglecting the presence of the secondary. Figure 1 shows the system at a state before the influence of the secondary acts. It indicates that—even without the influence of the ram pressure of the secondary wind—the primary wind will be pushed back toward the primary star merely by the acceleration of the secondary's radiation field. For the case of weak coupling, the distance where the two oppositely oriented acceleration components cancel out is merely $\sim 35\,{R}_{\odot }$ above the stellar surface for periastron and $\sim 85\,{R}_{\odot }$ for apastron. After including the presence of the secondary's wind (and not just the radiative influence of the secondary star) the wind ram pressure will push the WCR even closer toward the O star. For the case of strong coupling the distance of equal radiative acceleration components is $\sim 83\,{R}_{\odot }$ from the stellar surface for periastron and $\sim 180\,{R}_{\odot }$ for apastron configuration. There is a difference of more than a factor of 2 between the respective values for strong and weak coupling.

This significant difference between weak and strong coupling is shown in Figure 2—now including the presence of the secondary's wind. In both the apastron and the periastron state, weak coupling leads to a narrow shock-cone of the primary wind that is engulfed by a dominant secondary wind from the WR star. However, by using strong coupling, we find a much larger opening angle. The wind collision happens also farther away from the primary with strong turbulence in the apastron configuration. The blue-shaded regions of the secondary wind between the two stars show effects of radiative braking where the velocity of the secondary wind is significantly lowered by the radiation field of the primary star. The narrow blue region on the far side of the WR star results from the effect of shadowing in this region. The primary star is obscured by the secondary. Therefore the wind plasma in the shadow of the secondary does not feel any radiative acceleration component caused by the primary.

Zoom In Zoom Out Reset image size

The choice of weak or strong coupling not only influences the distance of the WCR from the primary star, but also severely alters the opening half-angle of the shock cone: $\sim 24^\circ$ in the case of weak coupling, $\sim 72^\circ$ in the case of strong coupling. As already mentioned in Secton 1, there is evidence from high-resolution X-ray spectroscopy of Chandra data (Henley et al. 2005) that ${\gamma }^{2}$ Velorum exhibits a large shock-cone opening half-angle of $\sim 85^\circ$. In the same study, the authors remark that simulations hitherto failed to reproduce such a feature, which might be due to significant radiative braking. In our simulation, including the effects of radiative braking with strong coupling of the CAK parameters, we can indeed reproduce this large shock-cone opening angle, which is most evident close to the periastron passage (see Figure 2, lower panel). Consequently, sufficient agreement with observations from X-ray spectroscopy can only be achieved by using strong coupling to the wind. We therefore consistently use this method in the following analysis.

The remaining fluid properties of density, temperature, and magnetic field strength are shown in Figure 3. The density of the wind plasma has a similar structure as the velocity: a thin laminar WCR for periastron and a turbulent WCR for apastron configuration. Maximum densities inside the WCR are $\sim 5\times {10}^{16}$ m⁻³ at the apex of both cases. The wind plasma reaches temperatures up to  5 × 10⁷ K in the turbulent WCR for the apastron configuration. Periastron temperatures are slightly lower and reach a maximum of just below 10⁶ K in the inner regions of the laminar WCR. The apex of the periastron region remains cool due to severe radiative cooling in the high-density environment of the thin WCR. Regarding the magnetic field, the high field strengths of the O star dipole strongly influence the conditions at the WCR, where field strengths of 10⁻⁵ T are reached. Electrons accelerated close to the apex of the WCR will suffer from severe synchrotron losses due to the influence of the O star dipole field. This is also valid for the periastron configuration, where the more laminar WCR is also close to regions with high magnetic field strength. While inverse-Compton losses are clearly dominant in the x–y plane where magnetic fields are low, the order reverses in regions of higher magnetic field strengths, where losses by synchrotron emission dominate all others.

Zoom In Zoom Out Reset image size

By solving the transport equation (along with the diffusion equation and the MHD equations including the spatial convection) for 100 logarithmically spaced energy bins between 1 MeV and 10 TeV for electrons and protons each, we obtain populations of high-energy particles at every location throughout the simulated volume. Details are given in Reitberger et al. (2014a, 2014b). As for the updated treatment of diffusion, we first only study the effects of the two parameters δ and D₀ for a single grid cell.

Figure 4 shows the resulting proton and electron spectra at the apex of the WCR in ${\gamma }^{2}$ Velorum for different values of D₀. For the results shown in Figure 4(a), diffusion is constant ($\delta =0,D={D}_{0}$). In Figures 4(b) and (c) the diffusion follows a power law in energy with index 0.3 ($D(E)={E}^{0.3}$). Figure 4(c) is the only figure where the advection of the particles with the flow of the wind is taken into account. It has been deactivated for Figures 4(a) and (b) in order to make the effects of different diffusion setups visible. Below these three particle density plots, two auxiliary plots show the various energy-loss and energy-gain rates that are taken into account at the location for which the spectra are shown.

Zoom In Zoom Out Reset image size

The proton spectra of Figure 4(a) indicate that—without the implementation of either Bohm diffusion or energy-dependent diffusion—the spectra have no cutoff. This is also apparent from the relative strengths of the energy-loss and energy-gain rates in Figure 4(e). The acceleration (solid red) dominates the combined losses for all energies. The spectral index can be controlled by the normalization of the diffusion coefficient D₀. In our model we find that for ${D}_{0}\geqslant {10}^{16}$ m² s⁻¹ and $\delta =0.3$, protons will not reach energies above 10 GeV. At this energy they simply diffuse out of the region in which further acceleration would be possible. The change in spectral index at ∼100 MeV—most visible in the two spectra for stronger diffusion in Figure 4(a)—can be readily explained by the change in Coulomb losses for protons, which roughly become constant above this energy. As the difference between energy gains and losses in Figure 4(e) increases, the proton spectra become harder.

For the electrons in Figure 4(a), the situation is very different. The cutoff reached just above 10 MeV for the case of low diffusion is caused by the strong inverse-Compton losses in the vicinity of the O star. Plot d) illustrates that acceleration is dominant only in a small energy range below 10 MeV. For the cases of higher diffusion, we find—as for protons—a softening of the spectra. Together with the severe energy losses, this leads to an earlier cutoff.

Note that the situation for the electrons changes somewhat as larger distances from the stars are reached in the outer wings of the WCR. As the effect of inverse-Compton losses subsides, electron energies of several 100 MeV can be reached. However—considering the magnetic dipole field—synchrotron losses will become dominant in the direction of the poles of the stellar dipoles and lead to an even lower energetic cutoff than caused by inverse-Compton losses. While an accurate representation of the magnetic field (as provided by the full-MHD simulation) is not of vital importance near the apex, where losses linked to the stellar radiation fields dominate, it is of high relevance at other regions of the WCR.

In Figure 4(b) the diffusion is not constant anymore but grows with ${E}^{0.3}$. This produces a cutoff—even if the diffusion is low at low energies. The graph illustrates that an implementation of a Bohm-diffusion related cutoff is no longer needed for the case of energy-dependent diffusion. It also emphasizes the role of D₀ and δ as important parameters to determine the maximum energy reached by protons. For electrons these parameters would only become relevant in case of very low loss-rates, as might be the case for systems where the WCR is far from the stars, but certainly not in ${\gamma }^{2}$ Velorum.

Figure 4(c) includes the effect of advection, while in the previous discussion, the particles have been allowed to diffuse freely away from the shock front without being affected by the flow of the wind. It has to be stressed that the WCR in a CWB system cannot be described by the properties of a typical 1D shock. An important difference is the significant velocity component perpendicular to the shock normal in the post-shock wind, which leads to particle transport in the downstream direction. The effects are seen in the figure. Including advection, we find a significant softening of all spectra. However—the maximum energies that are reached are still very close to the case without advection in it. Although the spectra are much softer at low energies (compared to Figure 4(b)), the spectral indices are similar at higher energies where diffusion dominates and the effects of advection disappear. The apparent irregularities for the case of very low diffusion (wiggles in the red solid spectra in Figure 4(c)) are explained by the local turbulence of the plasma, which leaves an imprint on the spectra. This effect can be seen for low diffusion and is enhanced if advection is dominant. For the spectra with ${D}_{0}={10}^{14}$ m² s⁻¹ in Figure 4(c), the effect of the vanishing influence of Compton losses and the corresponding hardening of the spectrum above ∼100 MeV becomes apparent once again.

Figure 5 explores the spatial distribution of electrons and protons at different energies. To simplify matters, only apastron conditions are shown here. Apparently, the strong stellar radiation field close to the apex of the WCR leads to severe inverse-Compton losses, which prevent electrons from reaching higher energies than ∼10 MeV at the apex and ∼100 MeV in the outer wings. The proton densities illustrate the effect of energy-dependent diffusion. The higher the particle energy, the larger the populated region. The turbulent structure of the apastron WCR disappears at higher energies as small-scale variations are smoothed out by the dominant diffusion.

Zoom In Zoom Out Reset image size

Figure 6 shows the highest particle energies reached. A clear difference can be seen between the electron energies in the x–y (left) and x–z (center) plane. This is due to the effect of synchrotron losses caused by the dipolar shape of the magnetic field. Around the dipole equator in the x–y plane, the field strengths are low and electrons therefore reach higher energies. In the center plot, large parts of the WCR are in close vicinity to the lobes of the dipole field, which causes severe losses by synchrotron emission. Electron energies decrease accordingly. For protons there is no such effect. Diffusion allows for a population of protons up to ∼1 TeV around the apex of the WCR where the acceleration of particles is most efficient. In regions farther downstream, the maximum energies are generally lower. This is a combined effect of less acceleration due to lower velocity gradients at the shock and of energy loss by collisions and Coulomb losses as the protons move downstream.

Zoom In Zoom Out Reset image size

Figure 8 (left) shows the resulting particle spectra when integrating over the simulated volume. The agreement with the previously discussed single-cell spectra and contour maps of particle densities is apparent.

Processing the obtained particle spectra (for ${D}_{0}=8\,\times {10}^{-14}$ m² s⁻¹, $\delta =0.3$ and proton injection ratio ${\eta }_{p}={10}^{-3}$) and the MHD variables with the schemes presented in Reitberger et al. (2014a), we obtain 2D projection maps, spectral energy distributions (SEDs), and integrated flux values for high-energy emission via inverse-Compton scattering, bremsstrahlung emission, and neutral-pion decay. Owing to its anisotropic nature, the inverse-Compton component is sensitive to the viewing angle which is determined by the parameters as listed in Section 3.1. A schematic view of the topology of the computational domain is shown in Figure 7 (left), which depicts the orbital plane with the stars in apastron configuration along with various viewing angles. The red arrow marks the viewing angle that is suggested by observations (Schmutz et al. 1997). Figure 7 (right) shows a projected emission map for the neutral-pion component if seen along the red line by an observer located on Earth. It becomes clear that the bulk of the emission stems from the apex of the WCR where high-energy proton densities as well as wind plasma densities are highest. We estimate the diameter of the region where >99% of the emission above 100 MeV occurs to be ∼700 ${R}_{\odot }$. It is smaller for periastron. The analogous plot for bremsstrahlung and inverse-Compton components would be empty. Owing to the low maximum energy of the electrons, the leptonic channels do not emit any γ-rays above 100 MeV. This can be seen in the SEDs in Figure 8 (middle), which represent the case of apastron. Even if the orbital motion of the system is turned off (as throughout this study), the turbulence of the WCR will lead to minor fluctuations in the SED. This is indicated by the gray-shaded area in Figure 8 (right), which shows the range of fluctuation for the neutral-pion component. We also modeled the system of ${\gamma }^{2}$ Velorum using the alternative WR mass-loss rate of ${\dot{M}}_{\mathrm{WR}}=8\times {10}^{-6}$ ${M}_{\odot }$ yr⁻¹. This 3.75 times lower mass-loss rate leads to a significant reduction of the maximum plasma densities on the WR side of the WCR. Consequently, the high-energy proton densities become lower as well. A model result that comes close to the measured data points demands an injection ratio of ${\eta }_{p}\sim 1$, which is highly problematic. However, for ${\dot{M}}_{\mathrm{WR}}=3\times {10}^{-5}$ ${M}_{\odot }$ yr⁻¹, a far more realistic injection rate of ${\eta }_{p}={10}^{-3}$ is sufficient to reproduce the measured spectrum. The total integrated γ-ray flux emitted above 100 MeV via the channel of neutral-pion decay is $\sim 5.6\times {10}^{-5}$ m⁻² s⁻¹, above 10 GeV it is $\sim 1.4\times {10}^{-7}$ m⁻² s⁻¹. The bremsstrahlung and inverse-Compton channel are only relevant at lower energies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The process of photon-photon absorption in the stellar radiation fields has been taken into account throughout the analysis. This has been done by integrating the optical depths $d\tau$ inside the simulated volume along the line of sight (as detailed in Reitberger et al. 2014a). The emissivities are then lowered by a factor of ${e}^{-\tau }$. As the difference to a blackbody spectrum is small, we assume the stellar photons to be monochromatic. From the simple estimate ${E}_{T}{E}_{\gamma }\geqslant {({m}_{e}{c}^{2})}^{2}$, one deduces that the process becomes relevant for γ-ray photons at ${E}_{\gamma }\gtrsim 50$ GeV (assuming monochromatic seed photons at ${E}_{T}/{k}_{B}\,=$ 56,000 K). Some of the spectra from neutral-pion decay reach maximum energies of ∼100 GeV. The effect of photon-photon absorption is clearly visible in the SED of Figure 8 (center) as a softening above $\sim 50\,$ GeV. It has no effect in the energy regime of the Fermi-LAT data points.

The data presented by Pshirkov (2016) suggests a γ-ray emissivity of $\sim 5\times {10}^{31}$ erg s⁻¹. This is consistent with our model results, which yield a total γ-ray emissivity of $\sim 5.6\times {10}^{31}$ erg s⁻¹. The presented projected emission maps for neutral-pion decay suggest that the γ-rays are primarily emitted on the WR side of the shock where plasma densities as well as shock velocities are high. The bolometric luminosity of the WR star is $\sim 6.5\times {10}^{38}$ erg s⁻¹. From the given mass-loss rate and terminal velocity, one can compute that $\sim 2\times {10}^{37}$ erg s⁻¹ or roughly 3% of L_bol are given as the kinetic energy of the wind. However, due to radiative braking, ${v}_{\infty }$ is not reached in the region relevant for the emission. Extracting density and velocity values from the simulation output, we compute the total kinetic energy available in this region. The sum yields $\sim 5\times {10}^{35}$ erg s⁻¹ or roughly 0.1% of L_bol. If only one ten-thousandth of this energy budget were used for the emitted spectrum, we could account for the signal we see.

The analysis presented in the previous sections focused on the orbital state of apastron (stellar separation $d\sim 344{R}_{\odot }$), which is not necessarily representative for the whole orbit. In order to compare matching predictions of our simulations with the measured data, other orbital states have to be taken into account. We therefore selected five orbital phases ($\phi =0,0.1,0.2,0.3,0.4$) for which the same analysis steps as described above for the case of apastron have been carried out.

We find a significant decrease of γ-ray flux toward the periastron state (ϕ = 0, $d=172\,{R}_{\odot }$). This is due to several factors, such as lower shock velocities, a reduced volume of the WCR, and increased Coulomb losses due to higher densities as the stars draw nearer toward each other.

For the periastron phase, the measured flux levels cannot be reached by tuning the free parameters of injection rate and normalization of diffusion coefficient (within reasonable ranges). However, a spectrum that lies within the gray-shaded area shown in Figure 8 (right) and is thus consistent with the data can easily be found for the orbital state $\phi =0.3$ (with the parameters of ${D}_{0}=3\times {10}^{-14}$ m² s⁻¹, δ = 0.3, and ${\eta }_{p}=6\times {10}^{-3}$). At this phase, the stellar separation nicely corresponds to the average distance obtained by integrating over the whole orbit. For the same set of parameters, the apastron flux lies considerably above the measured data, the periastron flux considerably below.

After obtaining the magnetic field and the electron spectra throughout the simulated volume, we apply the formalisms detailed by Blumenthal & Gould (1970) to reach a prediction for the synchrotron emission at radio wavelengths. At 0.2 GHz the resulting spectrum lies ∼8 orders of magnitude below the ATCA observations of ${\gamma }^{2}$ Velorum. This is due to the rather low chosen surface magnetic field strength of 10⁻³ T and the early cutoff of the electron spectra. For a higher surface magnetic field of 10⁻² T (while keeping D₀, δ, and ${\eta }_{p}$ constant at the values stated above for apastron passage), the predicted synchrotron emission is merely one order of magnitude below the ATCA observations. However, the difference between model results and the measured data points increases for higher wavelengths as the simulated synchrotron spectrum decreases much more steeply than the observed signal.

An alternative set of parameters in conjunction with smaller inverse-Compton losses may bring the predicted nonthermal radio emission in principal agreement with observations. Still, this fails to address the composite nature of the radio intensity between thermal and nonthermal radio emission, with the former considered being dominant (Chapman et al. 1999). A strong limit is only imposed in a way that our synchrotron prediction cannot supersede the total measured radio intensity, which clearly is not the case.

The transport Equations (13), (15) have been investigated by several authors (see Axford et al. 1977; Krymskii 1977; Blandford & Ostriker 1978) in the context of particle acceleration at an infinitely extended shock perpendicular to the magnetic field in the plasma. By assuming constant upstream and downstream velocity, matching the respective solutions of the transport equation at the position of the shock shows that particles injected at the shock obtain a power-law spectrum.

For the same scenario, it is also possible to investigate the temporal change of the energy-dependent particle flux. This is, e.g., extensively discussed in Drury (1983), where the acceleration rate of the energetic particles can be found by applying a Laplace transform to the transport equation. Here, we recapture the discussion from that study for the specific case of our numerical setup.

To investigate the acceleration rate, we consider an infinitely extended shock with homogeneous upstream and downstream plasmas. Additionally assuming spatially constant diffusion, Equation (15) can be used in the form (see Drury 1983)

$$\begin{eqnarray}&&\displaystyle \frac{\partial j}{\partial t}-D\displaystyle \frac{{\partial }^{2}j}{\partial x}+\displaystyle \frac{\partial }{\partial x}({uj})-\displaystyle \frac{1}{3}\displaystyle \frac{\partial u}{\partial x}\displaystyle \frac{\partial }{\partial p}({pj})=q(x,p).\end{eqnarray} \tag{ 16 }$$

When we assume injection to be localized at the shock at a given momentum ($q(x,p)={q}_{0}\delta (p-{p}_{0})\delta (x-{x}_{\mathrm{shock}})$, this equation becomes considerably simpler in the homogeneous upstream or downstream medium:

$$\begin{eqnarray}&&\displaystyle \frac{\partial j}{\partial t}-D\displaystyle \frac{{\partial }^{2}j}{\partial x}+{u}_{i}\displaystyle \frac{\partial j}{\partial x}=0,\end{eqnarray} \tag{ 17 }$$

where we explicitly used the constant velocity in the upstream u_i = u_u and the downstream u_i = u_d regions.

The Laplace transform

$$\begin{eqnarray}&&J(x,p,s)=\underset{0}{\overset{\infty }{\int }}{e}^{-{st}}j(t,x,p)dt\end{eqnarray} \tag{ 18 }$$

applied to Equation (17) gives

$$\begin{eqnarray}&&{sJ}+{u}_{i}\displaystyle \frac{\partial J}{\partial x}-D\displaystyle \frac{{\partial }^{2}J}{\partial {x}^{2}}=0,\end{eqnarray} \tag{ 19 }$$

where $j(t=0,x,p)=0$. Since the particles originate and are accelerated at the shock, J = 0 for $x\to \pm \infty$. Thus, the solution of Equation (19) is

$$\begin{eqnarray}&&J\propto {e}^{{\beta }_{i}x}\qquad \mathrm{with}\qquad {\beta }_{i}=\displaystyle \frac{{u}_{i}}{2D}\left(1\pm {\left(1+\displaystyle \frac{4{Ds}}{{u}_{i}^{2}}\right)}^{1/2}\right),\end{eqnarray} \tag{ 20 }$$

with $(+)$ for the upstream and $(-)$ for the downstream case.

Now, the downstream and upstream solutions need to be matched at the shock position. First, the flux needs to be continuous there (see Drury 1983), implying

$$\begin{eqnarray}&&{J}_{u}(x\to {x}_{\mathrm{shock}})={J}_{d}(x\to {x}_{\mathrm{shock}}).\end{eqnarray} \tag{ 21 }$$

Second, an additional constraint is found by integrating Equation (16) from a position just upstream to one just downstream of the shock. When we consider that j is also continuous at the shock, this leads to

$$\begin{eqnarray}&&-D{\left.\displaystyle \frac{\partial j}{\partial x}\right|}_{d}+D{\left.\displaystyle \frac{\partial j}{\partial x}\right|}_{u}+{u}_{d}j-{u}_{u}j-\displaystyle \frac{1}{3}\displaystyle \frac{\partial ({pj})}{\partial p}({u}_{d}-{u}_{u})\\ &&\quad =\,-D{\left.\displaystyle \frac{\partial j}{\partial x}\right|}_{d}+D{\left.\displaystyle \frac{\partial j}{\partial x}\right|}_{u}+\displaystyle \frac{2}{3}({u}_{d}j-{u}_{u}j)\\ &&\qquad -\,\displaystyle \frac{1}{3}p\displaystyle \frac{\partial j}{\partial p}({u}_{d}-{u}_{u})={{pq}}_{0}\delta (p-{p}_{0}).\end{eqnarray} \tag{ 22 }$$

Taking the Laplace transform of this at the position of the shock leads to

$$\begin{eqnarray}&&D({\beta }_{u}-{\beta }_{d}){J}_{0}-\displaystyle \frac{2}{3}({u}_{u}-{u}_{d}){J}_{0}+\displaystyle \frac{1}{3}p\displaystyle \frac{\partial {J}_{0}}{\partial p}({u}_{u}-{u}_{d})\\ &&\quad =\,{\displaystyle \int }_{0}^{\infty }{q}_{0}\delta (p-{p}_{0}){e}^{-{st}}{dt}=\displaystyle \frac{{{pq}}_{0}\delta (p-{p}_{0})}{s},\end{eqnarray} \tag{ 23 }$$

where ${J}_{0}=J(s,0,p)$ is the Laplace transform at the shock. In analogy to the derivation in Drury (1983), we introduce

$$\begin{eqnarray}&&{A}_{i}(s)=\left({\left(1+\displaystyle \frac{4{Ds}}{{u}_{i}^{2}}\right)}^{1/2}-1\right).\end{eqnarray} \tag{ 24 }$$

With this, it can be shown that Equation (23) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}({u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)){J}_{0}+\displaystyle \frac{1}{3}({u}_{u}+2{u}_{d}){J}_{0}\\ &&\quad +\,\displaystyle \frac{1}{3}p\displaystyle \frac{\partial {J}_{0}}{\partial p}({u}_{u}-{u}_{d})=\displaystyle \frac{{{pq}}_{0}\delta (p-{p}_{0})}{s}.\end{eqnarray} \tag{ 25 }$$

Considering that the inhomogeneity is only non zero at $p={p}_{0}$ yields the solution

$$\begin{eqnarray}{J}_{0}(s,{p}_{1}) & = & \displaystyle \frac{3{q}_{0}}{s({u}_{u}-{u}_{d})}{\left(\displaystyle \frac{{p}_{1}}{{p}_{0}}\right)}^{-\tfrac{{u}_{u}+2{u}_{d}}{{u}_{u}-{u}_{d}}}\\ & & \times \exp \left(-\underset{{p}_{0}}{\overset{{p}_{1}}{\displaystyle \int }}\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)}{{u}_{u}-{u}_{d}}\displaystyle \frac{{dp}}{p}\right).\end{eqnarray} \tag{ 26 }$$

From this the time-dependent particle flux at the shock follows from the related Bromwich integral:

$$\begin{eqnarray}&&{j}_{0}(t,p)=j(t,x=0,p)=\displaystyle \frac{1}{2\pi {\imath }}\underset{\eta -{\imath }\infty }{\overset{\eta +{\imath }\infty }{\int }}{J}_{0}(s,p){e}^{{ts}}{ds},\end{eqnarray} \tag{ 27 }$$

where $\eta \to 0$ since the highest value of s for any singularity is s = 0, here. The asymptotic behavior is found from the contribution of the residual with the largest real part, which is the simple pole at s = 0, leading to

$$\begin{eqnarray}&&{j}_{0}(t\to \infty ,p)={j}_{0}(\infty ,p)=\displaystyle \frac{3{q}_{0}}{{u}_{u}-{u}_{d}}{\left(\displaystyle \frac{{p}_{1}}{{p}_{0}}\right)}^{-\tfrac{{u}_{u}+2{u}_{d}}{{u}_{u}-{u}_{d}}}\propto {p}_{1}^{-a}\\ &&\qquad \mathrm{with}\qquad a=\displaystyle \frac{{c}_{r}+2}{{c}_{r}-1}.\end{eqnarray} \tag{ 28 }$$

This leads to the usual result that for a compression ratio of ${c}_{r}\to 4$, the spectral index becomes $a\to 2$. Additional consideration of the spatial dependence upstream of the shock in the limit $t\to \infty$ shows that Equation (20) leads to

$$\begin{eqnarray}&&j\propto {e}^{{\beta }_{u}x}\qquad \,{\beta }_{u}=\displaystyle \frac{{u}_{u}}{D},\end{eqnarray} \tag{ 29 }$$

representing the usual exponential decrease in the upstream direction.

In the general case, the inverse transform reads

$$\begin{eqnarray}{j}_{0}(t,p) & = & \displaystyle \frac{1}{2\pi {\imath }}{\displaystyle \int }_{\eta -{\imath }\infty }^{\eta +{\imath }\infty }\displaystyle \frac{3{q}_{0}}{s({u}_{u}-{u}_{d})}{\left(\displaystyle \frac{{p}_{1}}{{p}_{0}}\right)}^{-\tfrac{{u}_{u}+2{u}_{d}}{{u}_{u}-{u}_{d}}}\\ & & \times \exp \left(-\underset{{p}_{0}}{\overset{{p}_{1}}{\displaystyle \int }}\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)}{{u}_{u}-{u}_{d}}\displaystyle \frac{{dp}}{p}\right){e}^{{ts}}{ds}\\ & = & \displaystyle \frac{1}{2\pi {\imath }}\underset{0}{\overset{t}{\displaystyle \int }}{dt}^{\prime} {\displaystyle \int }_{\eta -{\imath }\infty }^{\eta +{\imath }\infty }\displaystyle \frac{3{q}_{0}}{{u}_{u}-{u}_{d}}{\left(\displaystyle \frac{{p}_{1}}{{p}_{0}}\right)}^{-\tfrac{{u}_{u}+2{u}_{d}}{{u}_{u}-{u}_{d}}}\\ & & \times \exp \left(-\underset{{p}_{0}}{\overset{{p}_{1}}{\displaystyle \int }}\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)}{{u}_{u}-{u}_{d}}\displaystyle \frac{{dp}}{p}\right){e}^{t^{\prime} s}{ds}\\ & = & \displaystyle \frac{1}{2\pi {\imath }}\underset{0}{\overset{t}{\displaystyle \int }}{dt}^{\prime} {\displaystyle \int }_{\eta -{\imath }\infty }^{\eta +{\imath }\infty }{j}_{0}(\infty ,p)\\ & & \times \exp \left(t^{\prime} s-\underset{{p}_{0}}{\overset{{p}_{1}}{\displaystyle \int }}\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)}{{u}_{u}-{u}_{d}}\displaystyle \frac{{dp}}{p}\right){ds}.\end{eqnarray} \tag{ 30 }$$

Apparently, this can be expressed as

$$\begin{eqnarray}&&{j}_{0}(t,p)={j}_{0}(\infty ,p)\underset{0}{\overset{t}{\int }}\phi (t^{\prime} ){dt}^{\prime} ,\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&\phi (t)=\displaystyle \frac{1}{2\pi {\imath }}\underset{\eta -{\imath }\infty }{\overset{\eta +{\imath }\infty }{\int }}{e}^{{ts}-h(s)}ds\end{eqnarray} \tag{ 32 }$$

and

$$\begin{eqnarray}&&h(s)=\underset{{p}_{0}}{\overset{{p}_{1}}{\int }}\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{u}{A}_{u}(s)+{u}_{d}{A}_{d}(s)}{{u}_{u}-{u}_{d}}\displaystyle \frac{{dp}}{p}.\end{eqnarray} \tag{ 33 }$$

From the definition of $\phi (t)$, it follows that

$$\begin{eqnarray}&&\underset{0}{\overset{\infty }{\int }}\phi (t){e}^{-{ts}}{dt}={e}^{-h(s)}.\end{eqnarray} \tag{ 34 }$$

Thus, $\phi (t)$ can be interpreted as the acceleration time distribution for a given p₀ and p₁ for h = 0, where the distribution is normalized (see Drury 1983). From the definition of h in Equation (33), this represents the case s = 0. Thus, differentiating Equation (34) with respect to s and setting s = 0 gives

$$\begin{eqnarray}&&-\underset{0}{\overset{\infty }{\int }}{tdt}=-{\left.\displaystyle \frac{\partial }{\partial s}h(s)\right|}_{s=0}.\end{eqnarray} \tag{ 35 }$$

This leads to the following expression for the mean acceleration time:

$$\begin{eqnarray}&&\langle t\rangle ={\left.\displaystyle \frac{\partial }{\partial s}h(s)\right|}_{s=0}=\displaystyle \frac{3}{{u}_{u}-{u}_{d}}\underset{{p}_{0}}{\overset{{p}_{1}}{\int }}D\left(\displaystyle \frac{1}{{u}_{u}}+\displaystyle \frac{1}{{u}_{d}}\right)\displaystyle \frac{{dp}}{p}.\end{eqnarray} \tag{ 36 }$$

This finally shows that the rate of momentum gain is given as

$$\begin{eqnarray}\dot{p} & = & \displaystyle \frac{p}{{t}_{{acc}}}\qquad \mathrm{with}\\ {t}_{{acc}} & = & \displaystyle \frac{3}{{u}_{u}-{u}_{d}}D\left(\displaystyle \frac{1}{{u}_{u}}+\displaystyle \frac{1}{{u}_{d}}\right)=\displaystyle \frac{3D}{{u}_{u}^{2}}\displaystyle \frac{(r+1)r}{r-1},\end{eqnarray} \tag{ 37 }$$

which is the usual expression used in many studies of diffusive shock acceleration.

The acceleration rate can also be found by investigating the term reflecting the adiabatic energy changes in the transport equation. When we consider only the first and the fourth term of Equation (15), this part of the transport equation can be viewed as an advection equation in momentum with an advection velocity:

$$\begin{eqnarray}&&{v}_{p}=\displaystyle \frac{1}{3}(({\rm{\nabla }}\cdot {\boldsymbol{u}})p),\end{eqnarray} \tag{ 38 }$$

representing the rate of momentum gain. When we evaluate this expression at the position of the shock, the velocity divergence is

$$\begin{eqnarray}{\rm{\nabla }}\cdot {\boldsymbol{u}} & = & \displaystyle \frac{{u}_{d}-{u}_{u}}{{\rm{\Delta }}x}=\displaystyle \frac{1}{{\rm{\Delta }}x}\displaystyle \frac{{u}_{d}-{u}_{u}}{{u}_{d}}{u}_{d}\\ & = & \displaystyle \frac{1}{{\rm{\Delta }}x}(1-r){u}_{d}\displaystyle \frac{{u}_{u}}{{u}_{u}}=\displaystyle \frac{1}{{\rm{\Delta }}x}\displaystyle \frac{1-r}{r}{u}_{u}.\end{eqnarray} \tag{ 39 }$$

The relevant length-scale can, additionally, be found from physical arguments. First, we noted the exponential decrease of the flux in the upstream direction. According to Blasi (2004), the point where the flux has decreased by a factor e can be interpreted as the position where the particles on average change direction. According to Equation (29), this happens at a distance

$$\begin{eqnarray}&&{x}_{u}={({\beta }_{u})}^{-1}=\displaystyle \frac{D}{{u}_{u}}\end{eqnarray} \tag{ 40 }$$

upstream of the shock. On the downstream side, such an estimate is no longer trivial since the distribution becomes spatially constant in the limit $t\to \infty$. Drury (1983), however, shows that the mean residence time upstream and downstream of the shock follow the same analytical expression, thus motivating

$$\begin{eqnarray}&&{x}_{d}=\displaystyle \frac{D}{{u}_{d}}.\end{eqnarray} \tag{ 41 }$$

Using ${\rm{\Delta }}x={x}_{u}+{x}_{d}$ then yields

$$\begin{eqnarray}{\rm{\nabla }}\cdot {\boldsymbol{u}} & = & \displaystyle \frac{1}{D}\displaystyle \frac{1-r}{r}{u}_{u}\displaystyle \frac{1}{\tfrac{1}{{u}_{u}}+\tfrac{1}{{u}_{d}}}=\displaystyle \frac{1}{D}\displaystyle \frac{1-r}{r}{u}_{u}\displaystyle \frac{{u}_{u}{u}_{d}}{{u}_{u}+{u}_{d}}\\ & = & \displaystyle \frac{1}{D}\displaystyle \frac{1-r}{r}{u}_{u}^{2}\displaystyle \frac{1}{r+1}.\end{eqnarray} \tag{ 42 }$$

Thus, we find for the momentum rate of change:

$$\begin{eqnarray}&&{v}_{p}=\displaystyle \frac{\partial p}{\partial t}=\displaystyle \frac{{u}_{u}^{2}}{3D}\displaystyle \frac{1-r}{r(r+1)}p,\end{eqnarray} \tag{ 43 }$$

with a resulting timescale that is identical to the one found in Equation (37).

Additionally, this shows the limits of the numerical representation of the transport equation. A shock computed in a numerical simulation has a finite thickness on the order of a few numerical cells. As long as this thickness is below ${\rm{\Delta }}x$ as computed above, the accelerated particles can be expected to experience the shock as a discontinuity, leading to the same acceleration rate as for an actual discontinuity. Deviations are only to be expected when the particle diffusion is so low that ${\rm{\Delta }}x$ becomes smaller than the simulated shock thickness.