A On 1{}^{\mathrm{st}}-order Fermi-acceleration

3D magneto-hydrodynamic models of non-thermal photon emission in the binary system  Velorum


Recent reports claiming association of the massive star binary system  Velorum (WR 11) with a high-energy -ray source observed by Fermi-LAT contrast the so-far exclusive role of  Carinae as the hitherto only detected -ray emitter in the source class of particle-accelerating colliding-wind binary (CWB) systems. We offer support to this claim of association by providing dedicated model predictions for the nonthermal photon emission spectrum of  Velorum.

We use three-dimensional magneto-hydrodynamic modeling (MHD) to investigate the structure and conditions of the wind-collision region (WCR) of  Velorum including the important effect of radiative braking in the stellar winds. A transport equation is then solved throughout the computational domain to study the propagation of relativistic electrons and protons. The resulting distributions of particles are subsequently used to compute nonthermal photon emission components.

In agreement with observation in X-ray spectroscopy, our simulations yield a large shock-cone opening angle. We find the nonthermal -ray emission of  Velorum to be of hadronic origin owing to the strong radiation fields in the binary system which inhibit the acceleration of electrons to energies sufficiently high for efficient inverse Compton radiation. We also discuss the strong dependence of a hadronic -ray component on the energy-dependent diffusion used in the simulations. Of two mass-loss rates for the WR star found in literature, only the higher one is able to accommodate the observed -ray spectrum with reasonable values for important simulation parameters such as the injection ratio of high-energy particles within the WCR.

acceleration of particles – binaries: general – gamma rays: stars – hydrodynamics – stars: winds, outflows

1 Introduction

The Velorum (WR 11) binary system is of interest for many different reasons. At a distance of 336 pc (North et al., 2007) it is the closest known binary system containing a WR star (Henley et al., 2005) and contains the brightest example of an WR type star with an O-star companion. The two stars combined reach 1.8 mag and are well visible with the naked eye at the Southern night sky. Velorum, along with its early-type stellar companion Velorum, is part of the Vela OB2 association, a group of 100 early-type stars spread over a diameter of pc at a mean distance of 410 12 pc (Jeffries et al., 2009).

The spectral types of the two stellar components of WR 11 are most probable WC8 and O7.5 (De Marco et al., 2000). Both stars orbit each other with a period of 78.53 0.01 days (Schmutz et al., 1997) with an eccentricity of 0.334 0.003 (North et al., 2007). From periastron to apastron the stellar separation of the two components varies from 172 to 344 R. Velorum is one of very few binary systems where direct determination of the stellar masses is possible and allows for critical tests of (WR and) massive star evolution. It constitutes an excellent laboratory for the study of radiatively-driven winds (Millour et al., 2007).

Velorum is observable at many wavelengths. In the radio regime, it is the strongest known thermal source among the class of hot luminous stars (Purton et al., 1982) and was first detected by Seaquist (1976) using the Parkes 64m telescope. More recent analysis with ATCA gave strong indication for a highly attenuated nonthermal radio component originating at the WCR between the two stars (Chapman et al., 1999). The attenuation of this nonthermal radio component is plausibly explained by the absorption of radio waves in the wind of the WR star (De Becker & Raucq, 2013). The observed indication for nonthermal radio waves via synchrotron emission suggests the presence of high-energy electrons accelerated in the WCR.

Velorum was observed by ISO at 2.3 to 29.6 m (van der Hucht et al., 1996) and by Keck in the K-band during an investigation of several WR systems among which Velorum and WR 140 were the only two objects without any observable dust emission throughout the system and its WCR (Monnier et al., 2002). More recent efforts applying long-baseline interferometry in near-IR with the AMBER/VLTI instrument revealed that the observed infrared emission primarily results from the contribution of the individual stars of the binary system. The stellar separation is too small to allow to spatially resolve the two components. Discrepancies between models and data may be resolved by additional free-free emission originating in the WCR (Millour et al., 2007).

Observations at optical wavelengths are manifold. Most relevant for this study are perhaps the observations with the HEROS spectrograph at the ESO 50 cm telescope (Schmutz et al., 1997; De Marco & Schmutz, 1999; De Marco et al., 2000) and its interpretation regarding system parameters, as well as more recent observations with the Sydney University Stellar Interferometer (North et al., 2007) which set the so far best constraints on many important stellar and stellar wind parameters by a new determination of the orbital parameters of the system. Interestingly, the determination of the WR-star’s mass-loss rate is ambiguous as different methods (polarimetric and radio-emission-based) yield different answers; the former yielding a mass-loss rate of 8.0 Myr, the latter of 3.0 Myr. The difference might plausibly stem from the effect of wind clumping on the radio-emission-based mass-loss rate which is proportional to the square of the wind density. The polarimetrically determined mass-loss rate is proportional to the density and therefore insensitive to clumping (North et al., 2007).

From observations at ultra-violet frequencies some details concerning the structure of WCR in Velorum can be derived. Using data from Copernicus and the International UltraViolet Explorer, St.-Louis et al. (1993) find clear evidence of eclipses of the O-star light caused by the WR-wind, as well as the presence of a wide cavity in the wind that is much closer to the O-star than to its companion. Thus, the first evidence for wind-wind collision in Velorum came for UV data.

The role of wind clumping as well as constraints on the opening-angle of the wind cavity have been further explored by X-ray observations. Velorum is a bright and well observable soft X-ray emitter and has been seen by ROSAT (Willis et al., 1995), ASCA (Rauw et al., 2000), Chandra (Skinner et al., 2001), and XMM-Newton (Schild et al., 2004). The latter study reveals a curious high and low state variation in the 1 to 8 keV regime that may originate from photoelectric absorption in the dense WR-wind whenever it eclipses the hot collisional plasma of the WCR. Whenever the WCR can be seen through a rarefied cavity that builds around and behind the O-star, the X-ray flux is clearly enhanced. Henley et al. (2005) use the measured X-ray variation in data provided by Chandra to estimate that the shock-cone opening half angle of the wind cavity must be rather large ( 85) – a finding which could not be reproduced in their hydrodynamical models which generally show a much narrower WCR. The same study concludes that the winds of both components will not reach terminal velocity before reaching collision. This puts additional constraints on hydrodynamical models as it requires more sophisticated wind implementations (e.g., radiative wind acceleration). The X-ray emission seen from Velorum is attributed to thermal emission reaching up to  keV. Investigations at higher energies towards the -ray regime by INTEGRAL have yielded upper limits (Tatischeff et al., 2004). The same authors also derive an upper limit on a possible nonthermal component in the X-ray data seen by ASCA.

At high-energy rays a recent study of 7 years of Fermi-LAT data shows a weak -ray signal (6.1 ) at 0.1 to 100 GeV at the position of the binary system (Pshirkov, 2016). Variability of the -ray flux consistent with Velorum’s orbital period could not yet be established due to low statistics. If this -ray source should indeed stem from the WCR in Velorum, it would be the second object detected within the source class of particle accelerating CWB systems - next to Carinae (Tavani et al., 2009; Abdo et al., 2010; Reitberger et al., 2015).

Velorum’s closeness to Earth and its well constrained stellar and stellar wind parameters make it a prime target for numerical modeling. The system’s short orbit, low eccentricity and lack of complicated gas and dust structures (like the Homunculus nebula around the Carinae) simplify modeling efforts and increase the reliability of the results.

It is the aim of this study to present a MHD model of Velorum which reproduces the above mentioned large shock-cone opening half angle derived from Chandra data and predicts a -ray signature to be compared with the observations by the Fermi-LAT. Insights gained by our model may be used to further constrain model parameters (injection ratios, diffusion coefficients, etc.) needed for the modeling of further or even more complex systems. The ability to successfully model the observed -ray emission will also support the claim of association.

The numerical modeling framework we use to model Velorum’s wind evolution, its particle populations and nonthermal emission is based on the work presented in Reitberger et al. (2014b), Reitberger et al. (2014a), and Kissmann et al. (2016). As these studies have so far applied the code to generic binary systems, the present study describes the first application of the code to a real astrophysical system. In Section 2 we detail and motivate all additions and alterations to the numerical modeling framework that have been implemented since the above mentioned publications, most notably the inclusion of “strong coupling” in the radiative wind-acceleration, an improved treatment for high-energy diffusion, and also for shock acceleration. In Section 3 we present our modeling results on the system of Velorum and compare it to observations. Finally, Section 4 provides a discussion of the results and an outlook on implications for future modeling of other CWB systems - most notable Carinae and WR 140.

2 Updates of the numerical modeling framework

2.1 Radiative wind-acceleration

Implementing stellar winds with a fixed terminal velocity does not reflect reality and is disfavored by observations showing that the winds in 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 Velorum system.

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


with the radial acceleration term


the optical depth parameter


the finite-disc correction


and the angle


The radial acceleration component is therefore a function of the distance to the star , its bolometric Luminosity , the velocity of the wind , the stellar radius and the standard CAK parameters and .


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’s radiation field on the primary star’s wind, and vice versa, will lead to radiative inhibition (star B’s radiation field reduces the wind acceleration in star A’s atmosphere) and radiative braking (rapid deceleration of star A’s 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 as Eq. 6 leaves room for two different interpretations.

The key question is whether the two CAK parameters, and , are related to the properties of the wind plasma or the star’s radiation field. Are they characterizing the radiation of the star (method A) or the capability of the wind material to react to radiation (method B)? In mathematical terms:


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 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. Whereas the physical interpretation of the parameter is clear – it represents the ratio of optical thin and optical thick lines in the wind plasma (Cranmer & Owocki, 1995) –, the meaning of the parameter is rather elusive, as it depends on the wind’s thermal temperature whilst 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, both method A and B 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), whereas 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 will show for the case of Velorum, to improved agreement with observations linked to the shape of the WCR.

2.2 Magnetohydrodynamics

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 one 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 () 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, due to its strong influence via synchrotron losses, it still has a significant effect on the electron population. As we will 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. Whereas a choice of  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 feature a nose-like structure as discussed in Kissmann et al. (2016). Accordingly, other free parameters 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.

2.3 Particle Acceleration

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


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


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.

2.4 Diffusion

Whereas 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 approximation often used 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 (Reitberger et al., 2014b, a; Reimer et al., 2006) 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 . The full equation (further motivated in Appendix 1) is


where is the differential number density of particles at energy at position , is the velocity of the wind plasma and is the injection rate of particles at energy . The term 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). Whereas 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 gain term (3) are solved in separate routines, both of which are capable of sub-cycling if a high diffusion coefficient (or a very high energy loss rate) should require it.

As term (1) in equation 11 indicates, diffusion is treated isotropically. The modifications allow now to consider an energy-dependent diffusion coefficient of the form . The exponent can be set to for a Kolmogorov type spectrum (recommended for 3D MHD simulations) and to for a Kraichnan type turbulence spectrum (Strong et al., 2007). We consistently use in this study (although it remains a free parameter in principle). The diffusion coefficient remains a free parameter.

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

2.5 Radiative cooling

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 (). Whereas 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


where and 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 (), a disregard of this correction due to composition would overestimate by a factor of 8. This leads to over-efficient cooling, resulting in a very dense and turbulent WCR. Once the above correction is applied we find a less turbulent WCR.

3 Results

3.1 Parameters

Stellar and stellar wind parameters used in our model of  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 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 (= M yr) or via radio-continuum observations (= M 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.

The orbital parameters are shown in Table 2. Our choice reflects the so far best constrained distance determination and reliable values for stellar separation and eccentricity. The inclination of the orbital plane of the  Velorum binary system and its argument of apastron are well constrained. Schmutz et al. (1997) find and . The exact definition of the angle 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 , and  days after apastron if . Thus, the meaning of the angle (as in Schmutz et al., 1997) is consistent with the definition of the angle in Reitberger et al. (2014a).

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

Star () () (K) () ( yr) (km s) - - T O7.5 28.5 [1,2] 17 [2] 35000 [1] 2.8 [2] [1] 2500 [1] 0.613 0.055 10 WC8 9 [2] 6 [2] 56000 [2] [2] [2,3] 1450 [2,3] 0.526 0.747 10 [1,2] 1.498

Table 1: Stellar and stellar wind parameters of  Velorum used in this study, taken from [1] De Marco et al. (2000), [2] North et al. (2007), [3] Schild et al. (2004)

distance semi-major axis period eccentricity inclination (pc) () (d) - () () 336 [1] 258 [1] 78.53 [2] 0.334 [2] 65 [2] 68 [2]

Table 2: Orbital parameters used in this study, taken from [1] North et al. (2007), [2] Schmutz et al. (1997)

3.2 MHD results

In a short-period binary system like  Velorum, the shape and plasma conditions of the WCR strongly depend on the choice of weak or strong coupling. Fig. 1 shows the possible components of the total line acceleration 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’s radiation field is shown in solid red, the acceleration acting on the primary wind by the secondary star’s radiation field is shown in dashed green for the case of weak coupling (, method A) and in dotted blue for the case of strong coupling (, method B).

The computation of for all three cases is based on the primary wind’s initial velocity profile neglecting the presence of the secondary. Fig. 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 towards 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 above the stellar surface for periastron and 85 for apastron. After including the presence of the secondary’s wind (and not just the secondary star’s radiative influence) the wind’s ram pressure will push the WCR even closer towards the O star. For the case of strong coupling the distance of equal radiative acceleration components is from the stellar surface for periastron, and for apastron configuration. There is 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 Fig. 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 which 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 further 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 primary star’s radiation field. 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.

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: in the case of weak coupling, 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  Velorum exhibits a large shock-cone opening half-angle of . 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 Fig. 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.

[]rad. acc. of primary (at position 0)rad. acc. of secondary / weak couplingrad. acc. of secondary / strong coupling(secondary located at -172 )[]rad. acc. of primary (at position 0)rad. acc. of secondary / weak couplingrad. acc. of secondary / strong coupling(secondary located at 344 )
Figure 1: Absolutes of the radiative acceleration in the primary wind close to the primary star (located at position 0) for periastron (left) and apastron (right) configuration. The accleration felt by the primary wind due to the radiation of the secondary star is indicated for the case of weak (green dashed line) and strong (blue dotted line) coupling. Arrows indicate the direction of the various acceleration components either towards or away from the primary. The intersections mark the points of balance where the radiative components are equal. s[][]method Amethod BWRWROOOWROWR
Figure 2: Absolute velocity of the wind plasma for periastron (left) and apastron (right) with weak coupling (upper panel) and strong coupling (lower panel). The approximate opening angle is indicated by the dashed yellow lines.

The remaining fluid properties of density, temperature, and magnetic field strength are shown in Fig. 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 5 m at the apex of both cases. The wind plasma reaches temperatures up to  510 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’s dipole strongly influences 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’s dipole field. This is also valid for periastron configuration where the more laminar WCR is also close to regions with high magnetic field strength. Whereas 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.

Figure 3: Number density of wind particles (left), temperature (center), and magnetic field strength (right) of  Velorum in the x–z plane for apastron (top) and periastron (bottom) configuration.

3.3 High-energy particle results

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. (2014b, a). As for the updated treatment of diffusion we first study the effects of the two parameters and for a single grid cell only.

Fig. 4 shows the resulting proton and electron spectra at the apex of the WCR in Velorum for different values of . For the results shown in Fig. 4 a) diffusion is constant (). In Fig. 4 b) and c) the diffusion follows a power law in energy with index 0.3 (). Fig. 4 c) is the only one where the advection of the particles with the flow of the wind is taken into account. It has been deactivated for Fig. 4 a) and 4 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.

The proton spectra of Fig. 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 gain-rates in Fig. 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 . In our model we find that for m s and , protons will not reach energies above 10 GeV. At this energy they simply diffuse out of the region in which further acceleration were possible. The change in spectral index at 100 MeV – most visible in the two spectra for stronger diffusion in Fig. 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 Fig. 4 e) increases, the proton spectra become harder.

For the electrons in Fig. 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 Fig. 4 b) the diffusion is not constant anymore but grows with . 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 not needed anymore for the case of energy-dependent diffusion. It also emphasizes the role of 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  Velorum.

Fig. 4 c) includes the effect of advection whereas 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 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 Fig. 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 Fig. 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 ms in Fig. 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 4: Upper row: Differential number density of electrons and protons for different values of the diffusion coefficient . Electrons and protons can be distinguished by the different injection ratios (higher for protons) at . The values for are (red solid), (green dashed), (blue dotted), and m s (azur dot–dashed). We used and no advection for plot a), and no advection for plot b), and advection for plot c).
Lower row: Energy loss-rates for electrons in plot d) and for protons in plot e). The colors indicate the acceleration rate (red solid), inverse Compton losses (green dashed), synchrotron losses (blue dotted), bremsstrahlung losses (azur dot–dashed), and Coulomb losses (magenta double-dot–dashed) in plot d) and acceleration rate (red solid), losses by nucleon-nucleon collisions (green dashed), and Coulomb losses (blue dotted) in plot e).

Fig. 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 gets, the larger the populated region becomes. The turbulent structure of the apastron WCR disappears at higher energies as small-scale variations are smoothed out by the dominant diffusion.

Fig. 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 (centre) 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 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 cause 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 further 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.

Fig. 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.

1 MeV10 MeV10 MeV
1 MeVMeV m100 MeVMeV m10 GeVMeV m
Figure 5: Differential number density of electrons (top) and protons (bottom) in MeV m for different values of kinetic particle energy. The countor maps show the x–z plane of a 512 simulation at y=0.
Figure 6: Maximum particle energies for electrons in the plane (left), electrons in the x–z plane (middle) and protons in the x–z plane (right).

3.4 Gamma-ray predictions

Processing the obtained particle spectra (for ms, and proton injection ratio ) 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 Fig. 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 which is suggested by observations (Schmutz et al., 1997). Fig. 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. 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 Fig. 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 grey-shaded area in Fig. 8 (right) which shows the range of fluctuation for the neutral pion component. We also modeled the system of Velorum using the alternative WR mass-loss rate of M 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 too. A model result that comes close to the measured data points demands an injection ratio of 1 which is highly problematic. However, for  M yr, a far more realistic injection rate of =10 is sufficient to reproduce the measured spectrum. The total integrated -ray flux emitted above 100 MeV via the channel of neutral pion decay is  ms, above 10 GeV it is ms. The bremsstrahlung and inverse Compton channel are only relevant at lower energies.

The process of photon-photon absorption in the stellar radiation fields has been taken into account throughout the analysis. This has been done via integrating the optical depths 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 . As the difference to a black-body spectrum is small, we assume the stellar photons to be monochromatic. From the simple estimate one deduces that the process becomes relevant for -ray photons at  50 GeV (assuming monochromatic seed photons at 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 Fig. 8 (centre) as a softening above GeV. It has no effect in the energy regime of the Fermi-LAT data points.

0.005 arcsecFlux [ m s]
Figure 7: Left: view of the two stars within the computational domain. The line of centers is represented by the horizontal dashed orange line. Various viewing angles (lines-of-sight) are indicated, including the one that is used in our simulation: . Right: Projected flux above 100 MeV for neutral pion decay.
(Energy [MeV])(Energy [MeV])(Energy [MeV])
Figure 8: Left: number density of protons (green) and electrons (red) integrated over the computational domain. Middle: resulting SEDs for inverse Compton (solid red), bremsstrahlung (dashed green) and neutral pion decay emission channels, the latter without (dotted blue) and with photon-photon absorption (dot-dashed magenta). Also shown are the Fermi-LAT data points determined for Velorum by Pshirkov (2016) (black) and for Carinae by (Reitberger et al., 2015) (green) as well as the upper limit on the 0.5-10 keV emission of  Velorum derived from ASCA, and the INTEGRAL/IBIS upper limits (for both, see Tatischeff et al., 2004). Right: The gray-shaded area indicates the simulated fluctuation of the neutral pion decay induced SEDs due to the turbulence in the WCR (stellar separation fixed at apastron). The data points are the same as in the center plot.

3.5 Energetics considerations

The data presented by Pshirkov (2016) suggests a -ray emissivity of erg s. This is consistent with our model results which yield a total -ray emissivity of 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 erg s. From the given mass-loss rate and terminal velocity one can compute that erg s or roughly 3% of are given as the kinetic energy of the wind. However due to radiative braking, 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 erg s or roughly 0.1% of . If only one ten-thousandth of this energy budget were used for the emitted spectrum, we could account for the signal we see.

3.6 Orbital considerations

The analysis presented in the previous sections is focusing on the orbital state of apastron (stellar separation ) 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 () 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 towards the periastron state (=0, ). 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 towards each other.
For 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 withing the grey-shaded area shown in Fig. 8 (right) and is thus consistent with the data can easily be found for the orbital state (with the parameters of ms, =0.3, and ). 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, periastron flux considerably below.

3.7 Nonthermal radio emission

Having obtained 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  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 , , and 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 leaves the composite nature of the radio intensity between thermal and nonthermal radio emission unaddressed, 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.

4 Summary and Outlook

We present the first 3D MHD simulation of the Velorum binary system that successfully reproduces a wide opening angle of the WCR as suggested by observations (Henley et al., 2005). This is achieved by implementing the CAK radiative line acceleration in the limit of strong coupling of stellar and wind parameters. Through severe radiative braking the WR wind is slowed efficiently prior to hitting the WCR resulting in a large opening half-angle of 72.

In addition, our model can account for the observed -ray emission from Velorum via diffusive shock acceleration of protons at the WCR and the resulting neutral pion decay emission. Simulating the acceleration of charged particles at the WCR, we obtain proton distributions up to 1 TeV whereas electrons hardly reach 100 MeV. This is mainly due to the significant inverse Compton and synchrotron losses at the WCR. Maximum proton energies and flux levels are determined by our choice of diffusion coefficient, diffusion index and proton injection ratio. To obtain approximate agreement with measured -ray spectra, we choose ms, , and for the apastron phase of the system. Most orbital phases can be modeled to fit the data by changing these parameters within a reasonable range. For phases around periastron this is no longer possible. This suggests the presence of variability on orbital timescales the non-observation of which can be accounted for by the low statistics.

2D projection maps of the -ray emission region show that the bulk of the emission due to neutral pion decay is clearly confined to a region around the apex of the WCR of a diameter of roughly R. A larger computational domain would add little to the total -ray emission.

The observed data by Pshirkov (2016) can only be attained by rays of hadronic origin as the leptonic components are far too weak due to the low energy of the electrons. Of the two mass-loss rates for the WR star found in literature only the higher one is fit to reproduce the observed spectra with a realistic choice for the proton injection ratio. Myr would require which is considered unrealistically large. Owing to the turbulent WCR the output spectra do not converge to a definite flux level but rather fluctuate above and below the observed data.

This work lays the base line for further investigation of the Velorum via numerical simulation. The results on particle acceleration and -ray emission focus on the orbital phase of apastron but also show that similar agreement with observations can be achieved for other orbital states. In the near future we plan to compare numerical simulations of -ray emission in  Velorum with the much brighter -ray source of Carinae and the yet undetected but suspected -ray source WR 140. If we apply the same parameters regarding diffusion and particle injection as for  Velorum – can simulations account for the fact that one source is bright while the other remains dark on the Fermi -ray sky? We hope to provide an answer to this question in the near future.

The computational results presented have been achieved (in part) using the HPC infrastructure of the University of Innsbruck. A.R. acknowledges financial support from the Austrian Science Fund (FWF), project P 24926-N27.

Appendix A On 1-order Fermi-acceleration

In our modelling efforts we use the cosmic-ray transport equation following Parker (1965):


where phase space density of the energetic particles depending on time , position , and momentum where . Additionally, is the magnitude of spatial diffusion, is the advection velocity, and the source term. The measurable quantity is the particle flux . By using:


the transport equation for the flux is found to be:


where the source term relates to the one in Eq. A1 by . Note that this is analogous to Eq. 11 with the minor differences that the latter is expressed in terms of energy , rather than absolute momentum , includes energy losses , assumes isotropic diffusion, and specifies the source term to an injection at a given energy .

a.1 Diffusive shock acceleration

The transport equations (A1), (A3) have been investigated by several authors (see Krymskii, 1977; Axford et al., 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 will recapture the discussion from that study for the specific case of our numerical setup.

a.2 Acceleration Rate

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


Assuming injection to be localised at the shock at a given momentum ( this equation becomes considerably simpler in the homogeneous upstream or downstream medium:


where we explicitly used the constant velocity in the upstream and the downstream regions.

The Laplace transform:


applied to Eq. (A5) gives:


where . Since the particles originate and are accelerated at the shock, for . Thus, the solution of Eq. (A7) is:


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:


Secondly, an additional constraint is found by integrating Eq. (A4) from a position just upstream to one just downstream of the shock. Considering that is also continuous at the shock, this leads to:


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


where is the Laplace transform at the shock. In analogy to the derivation in Drury (1983) we introduce:


With this, it can be shown that Eq. (A.2) can be written as:


Considering that the inhomogeneity is only non zero at yields the solution:


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


where since the largest value of for any singularity is , here. The asymptotic behaviour is found from the contribution of the residual with the largest real part, which is the simple pole at leading to:


This leads to the usual result that for a compression ratio of the spectral index becomes . Additional consideration of the spatial dependence upstream of the shock in the limit shows that Eq. (A8) leads to


representing the usual exponential decrease in the upstream direction.

In the general case the inverse transform reads:


Apparently, this can be expressed as:






From the definition of it follows that:


Thus, can be interpreted as the acceleration time distribution for given and for , where the distribution is normalised (see Drury, 1983). From the definition of in Eq. (A21) this represents the case . Thus, differentiating Eq. A22 with respect to and setting gives:


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


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


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

a.3 A direct estimate of the acceleration rate

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


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


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


upstream of the shock. On the downstream side such an estimate is not trivial anymore, since the distribution becomes spatially constant in the limit . Drury (1983), however, shows that the mean residence time upstream and downstream of the shock follow the same analytical expression, thus motivating:


Using then yields:


Thus, we find for the momentum rate of change:


with a resulting time scale that is identical to the one found in Eq. (A25).

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 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’s diffusion is so low that becomes smaller than the simulated shock thickness.


  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 723, 649
  2. Axford, W. I., Leer, E., & Skadron, G. 1977, International Cosmic Ray Conference, 11, 132
  3. Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29
  4. Blasi, P. 2004, Astroparticle Physics, 21, 45
  5. Blumenthal, G. R. & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237
  6. Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  7. Chapman, J. M., Leitherer, C., Koribalski, B., Bouter, R., & Storey, M. 1999, ApJ, 518, 890
  8. Cranmer, S. R. & Owocki, S. P. 1995, ApJ, 440, 308
  9. De Becker, M., & Raucq, F. 2013, A&A, 558, A28
  10. De Marco, O., & Schmutz, W. 1999, A&A, 345, 163
  11. De Marco, O., Schmutz, W., Crowther, P. A., et al. 2000, A&A, 358, 187
  12. Drury, L. O. 1983, Reports on Progress in Physics, 46, 973
  13. Gayley, K. G., Owocki, S. P., & Cranmer, S. R. 1997, ApJ, 475, 786
  14. Henley, D. B., Stevens, I. R., & Pittard, J. M. 2005, MNRAS, 356, 1308
  15. Jeffries, R. D., Naylor, T., Walter, F. M., Pozzo, M. P., & Devey, C. R. 2009, MNRAS, 393, 538
  16. Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452
  17. Kissmann, R., Reitberger, K., Reimer, O., Reimer, A., & Grimaldo, E. 2016, ApJ, 831, 121
  18. Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, (, 1306
  19. Madura, T. I., et al. 2013, MNRAS, 436, 3820
  20. Millour, F., Petrov, R. G., Chesneau, O., et al. 2007, A&A, 464, 107
  21. Monnier, J. D., Greenhill, L. J., Tuthill, P. G., & Danchi, W. C. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 260, Interacting Winds from Massive Stars, ed. A. F. J. Moffat & N. St-Louis, 331
  22. North, J. R., Tuthill, P. G., Tango, W. J., & Davis, J. 2007, MNRAS, 377, 415
  23. Owocki, S. P., & Gayley, K. G. 1995, ApJ, 454, L145
  24. Parker, E. N. 1965, Planet. Space Sci., 13, 9
  25. Parkin, E. R., Pittard, J. M., Corcoran, M. F., & Hamaguchi, K. 2011, ApJ, 726, 105
  26. Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A& A, 164, 86
  27. Pittard, J. M. 2009, MNRAS, 396, 1743
  28. Pshirkov, M. S. 2016, MNRAS, 457, L99
  29. Purton, C. R., Feldman, P. A., Marsh, K. A., Allen, D. A., & Wright, A. E. 1982, MNRAS, 198, 321
  30. Rauw, G., Stevens, I. R., Pittard, J. M., & Corcoran, M. F. 2000, MNRAS, 316, 129
  31. Reimer, A., Pohl, M., & Reimer, O. 2006, ApJ, 644, 1118
  32. Reitberger, K., Kissmann, R., Reimer, A., & Reimer, O. 2014a, ApJ, 789, 87
  33. Reitberger, K., Kissmann, R., Reimer, A., & Reimer, O. 2017, in American Institute of Physics Conference Series, Vol. 1792
  34. Reitberger, K., Kissmann, R., Reimer, A., Reimer, O., & Dubus, G. 2014b, ApJ, 782, 96
  35. Reitberger, K., Reimer, A., Reimer, O., & Takahashi, H. 2015, A&A, 577, A100
  36. Schild, H., Güdel, M., Mewe, R., et al. 2004, A&A, 422, 177
  37. Schmutz, W., et al. 1997, A&A, 328, 219
  38. Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751
  39. Seaquist, E. R. 1976, ApJ, 203, L35
  40. Skinner, S. L., Güdel, M., Schmutz, W., & Stevens, I. R. 2001, ApJ, 558, L113
  41. St.-Louis, N., Willis, A. J., & Stevens, I. R. 1993, ApJ, 415, 298
  42. Stevens, I. R., & Pollock, A. M. T. 1994, MNRAS, 269, 226
  43. Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle Science, 57, 285
  44. Tatischeff, V., Terrier, R., & Lebrun, F. 2004, in ESA Special Publication, Vol. 552, 5th INTEGRAL Workshop on the INTEGRAL Universe, ed. V. Schoenfelder, G. Lichti, & C. Winkler, 409
  45. Tavani, M., et al. 2009, ApJ Letters, 698, L142
  46. Usov, V. V., & Melrose, D. B. 1992, ApJ, 395, 575
  47. van der Hucht, K. A., Morris, P. W., Williams, P. M., et al. 1996, A&A, 315, L193
  48. Willis, A. J., Schild, H., & Stevens, I. R. 1995, A&A, 298, 549