High-energy particle transport in 3D hydrodynamic models of colliding-wind binaries

8.1 MB·

by  K. Reitberger

High-energyparticletransport
in3Dhydrodynamicmodels
ofcolliding-windbinaries
K.Reitberger,R.Kissmann,A.Reimer,O.Reimer
Institutfu¨rAstro-undTeilchenphysikandInstitutfu¨rTheoretischePhysik,Leopold-Franzens-Universita¨t
4 Innsbruck,A-6020Innsbruck,Austria
1
0
and
2

n
G.Dubus
a
J
7 UJF-Grenoble1/CNRS-INSU,InstitutdePlane´tologieetd’AstrophysiquedeGrenoble(IPAG),UMR5274,
Grenoble,France
]
E
H [email protected]
.
h
p ABSTRACT
-
o
r
t
s Massivestarsinbinarysystems(asWR140,WR147orηCarinae)havelongbeenregarded
a
aspotentialsourcesofhigh-energyγ-rays. Theemissionisthoughttoariseintheregionwhere
[

the stellar winds collide and produce relativistic particles which subsequently might be able
1
v to emit γ-rays. Detailed numerical hydrodynamic simulations have already offered insight in
3
thecomplexdynamicsofthewindcollisionregion(WCR),whileindependentanalyticalstud-
2
3 ies, albeit with simplified descriptions of the WCR, have shed light on the spectra of charged
1
particles. Inthispaper, wedescribeacombinationofthesetwoapproaches. Wepresenta3D-
.
1
hydrodynamicalmodelforcollidingstellarwindsandcomputespectralenergydistributionsof
0
4 relativisticparticlesfortheresultingstructureoftheWCR.Thehydrodynamicpartofourmodel
1
incorporatestheline-drivenaccelerationofthewinds, gravity, orbitalmotionandtheradiative
:
v cooling of the shocked plasma. In our treatment of charged particles we consider diffusive
i
X shockaccelerationintheWCRandthesubsequentcoolingviainverseComptonlosses(includ-
r ingKlein-Nishinaeffects),bremsstrahlung,collisionsandotherenergylossmechanisms.
a
Subjectheadings: accelerationofparticles–binaries: general–gammarays: stars–hydrody-
namics–stars: winds,outflows
1. INTRODUCTION
Inrecentyears,severalmodelsofmostlyanalyticalnaturehaveaddressedthequestionwhetherbinary
systems of massive stars without compact objects are liable to produce high-energy γ-ray emission (e.g.,
–2–
Reimeretal.2006;Pittardetal.2006;Benaglia&Romero2003). Thesestudiesarguethatsuchobjects(i.e.
systemscontainingaWolf-Rayet(WR)starandanOB-typestar)providesuitableenvironmentsforefficient
particle acceleration and subsequent γ-ray emission. Electrons and protons are thought to be accelerated
at the shock fronts tracing the edge of the regions where stellar winds collide with supersonic velocities.
Several γ-ray emission mechanisms – inverse Compton (IC) scattering, relativistic bremsstrahlung and π0-
decay – can produce γ-rays at GeV and TeV energies with sufficient fluxes to allow for detection with
instrumentsastheFermi-LargeAreaTelescope(LAT)orperhapsevenHESS,MAGICandVERITAS.
In contrast to expectations, no such detection of γ-ray emission linked to colliding-wind binaries
(CWBs) has been reported so far, with one notable exception: the highly unusual object η Carinae is the
onlyCWBsystemunambiguouslylinkedtohigh-energyγ-rayemission(Reitbergeretal.2012). Thissys-
tem most likely consists of a WR star and a high-mass Luminous Blue Variable (LBV). Both stars are
envelopedinahugedustandgascloud,theHomunculus-nebula,whichoriginatedinamassiveoutburstof
theLBVintheyearof1843. Althoughthissourceexhibitsanumberofuniquecharacteristics,noexplana-
tionforitshighγ-rayflux–comparedtothenon-detectionofotherCWBs–canbegivensofar. Dedicated
observations of WR-OB systems as WR 147 or WR 140 (for which models predicted fluxes above LAT
detectionthresholds)haveyieldedupperlimits(Werneretal.2013).
Detailed 3D hydrodynamical (HD) simulations (Pittard 2009) have recently explored the highly dy-
namical nature of the WCR in CWBs and its strong dependence on stellar and orbital parameters. The
complexdensity,velocityandtemperaturestructureofthecollidingwindshavefurtherbeenusedtomodel
thethermalradioandX-rayemissioninsuchsystems(Pittard2010;Pittard&Parkin2010).
This work aims at a numerical computation of the spectral energy distribution of charged particles
within a numerical HD model of the WCR. By solving a transport equation including spatial convection,
diffusive shock acceleration (DSA) and various cooling processes for electrons and protons at every grid
point of the HDsimulation, we simulate the time-dependent 3D spatial distributionof particles at different
energies. AscoolingprocesseswetakeICemission,bremsstrahlung,nucleon-nucleoninteraction,adiabatic
cooling and others into account. The particles are injected at the shock fronts of the wind collision region
andsubsequentlygainenergybyDSA.
We developed a code which – by taking stellar, stellar wind and orbital parameters of a given binary
system as input – solves the 3D distribution of density, velocity and temperature fields of the wind plasma
as well as the energy spectra for electrons and protons for every point on a numerical grid. The resulting
spatiallyvaryingparticlespectrawillserveasinputforcomplementarystudies,computingthecomponents
ofensuingγ-rayemission.
In Section 2 we introduce the numerical and HD setup that we use to simulate a 3D distribution of
radiativelydrivenwindsinabinarysystem. Ourmethodofdealingwiththespatialandenergeticevolution
of particle species via solving a transport equation is thoroughly discussed in Section 3. In Section 4 we
present results obtained for a typical CWB system in terms of spatial and spectral distribution of high-
energyelectronsandprotons. Section5providesasummaryofourfindingsaswellasanoutlookonfuture
developments.
–3–
2. HYDRODYNAMICS
2.1. NumericalSetup
To simulate the hydrodynamics of the stellar winds we use the Cronos code (Kissmann et al. 2008;
Kleimann et al. 2009) which is a finite-volume magneto-hydrodynamic (MHD) code optimized for the
simulation of compressible astrophysical plasmas. The code is second-order accurate in space and allows
for Cartesian, cylindrical, and spherical grid layouts. Cronos uses approximate Riemann solvers for the
timeintegrationoftheHDandMHDequations,wherethechoiceofRiemannsolversincludesHarten-Lax-
vanLeer(hll, usedinthisstudy), hll-Contact(hllc)andhll-Discontinuities(hlld). Timeintegrationisdone
via a second- or third-order Runge-Kutta integrator, where the time discretisation utilises the semidiscrete
approach.
Inthepresentcase,CronosisusedtosolvetheHDequations,
∂ρ
+∇(cid:126) ·(ρ(cid:126)v) = 0 (1)
∂t
∂ρ(cid:126)v
+∇(cid:126) ·(ρ(cid:126)v(cid:126)v+Pˆ) = ρf(cid:126) (2)
∂t
(cid:32) (cid:33)2
∂(cid:15) ρ
+∇(cid:126) ·[((cid:15) +P)(cid:126)v] = Λ(T)+ρf(cid:126)·(cid:126)v (3)
∂t m
H
where the left-hand side of the equations is solved by Cronos intrinsically, whereas the right-hand side
(i.e. theforceterm f(cid:126)andtheradiativecoolingtermΛ(T)aredealtwithindedicatedadditionalmodulesas
outlinedbelow.
ρ is the mass density,(cid:126)v is the velocity vector, P the scalar pressure, f(cid:126)the sum of all external forces,
(cid:15) = ρv2+eisthetotalenergypervolume,etheinternalenergypervolumeandTthetemperature. Weuse
2
theidealgasequationofstateP = (γ−1)ewiththeadiabaticindexforamonoatomicgasγ = 5.
3
2.2. Theforceterm f(cid:126)
Theforcedensityterm f(cid:126)consistsofthreecomponents: gravity,radiativelineaccelerationduetheions
inthewindandtheradiativeaccelerationduetophotonsscatteringoffelectrons. Itcanbewrittenas
f(cid:126)= (cid:88)2 −GM∗,ir(cid:126)r3i +(cid:126)grLad,i+(cid:126)gerad,i (4)
i=1 i
where the index i indicates each star. The vector (cid:126)r is given relative to the star i. We assume radiative
i
accelerationtobedirectedradiallyfromthestars. Therefore,
(cid:126)r σ L
(cid:126)ge = ge i with ge = e ∗,i (5)
rad,i rad,iri rad,i 4πr2c
i
–4–
whereσ isthespecificelectronopacityduetoThomsonscatteringandL istheluminosityofstari.
e ∗,i
Determining the line-acceleration requires an integration over the finite stellar disk (as a point source
basedapproachwouldyieldsignificantlyerroneousresultsclosetothestar(seeLamers&Cassinelli1999).
This integral can be simplified by assuming azimuthal symmetry of the stellar disk. To approximate the
contribution of the wide spectrum of optically thick and thin spectral lines of the ions in the wind, we rely
on the standard Castor-Abbott-Klein (CAK) formalism first introduced by Castor et al. (1975) and later
modified and improved by Pauldrach et al. (1986). This approach replaces the required summation over
all lines by a simple parametrization with two parameters α and k. Taking all this into account, the line-
accelerationterm(asgiveninGayleyetal.1997)can–withsomemodification–beexpressedasfollows:
(cid:126)grLad,i = grLad,i(cid:126)rrii with grLad = 2kσc1e−α4πLR∗,i2∗,i (cid:90)0θ∗,in(cid:126)i·∇(cid:126)ρv(nt(cid:126)hi·(cid:126)v)αsinθcosθdθ (6)
withk andαbeingtheCAK-parametersmentionedabove. ρisthemassdensityofthewindinthecellfor
whichtheaccelerationiscalculated. v isthethermalvelocityofionizedhydrogen,theangleθ marksthe
th ∗,i
edgeofthestellardiskrelativetothepointforwhichwederivetheacceleration. Itisdefinedbysinθ = R∗,i.
∗,i ri
Theprojectedvelocitygradientn(cid:126) ·∇(cid:126)(n(cid:126) ·(cid:126)v)dependsontheunitvectortowardsapointonthestellarsurface
i i
n(cid:126) dependingagainonθ. InEuclideancoordinatesitcanbeexpressedas
i
(cid:32) (cid:33)
∂v ∂v ∂v ∂v ∂v
(cid:126)n·∇(cid:126)((cid:126)n·(cid:126)v) = n2 x +n2 y +n2 z +n n y + x +
x ∂x y ∂y z ∂z x y ∂x ∂y
(cid:32) (cid:33) (cid:32) (cid:33) (7)
∂v ∂v ∂v ∂v
+n n z + x +n n z + y
x z y z
∂x ∂z ∂y ∂z
In evaluating the integral in Eq. 6 we use(cid:126)n = sinθP(cid:126) +cosθN(cid:126), where N(cid:126) is the unit vector pointing
i
fromthegridpointtowardsthestellarcenter. P(cid:126) isaperpendicularvectorchosensuchthatitisnonzeroand
liesinaplaneofthenumericalgrid. IntegrationisthenperformednumericallyusingasimpleSimpson-rule
withfivestepsintheinterval[0,cosθ ].
∗,i
Duetotheionizationoftheplasma,linedrivingissettozeroincellswithtemperatureabove106K.
2.3. Radiativecooling
DuetohighdensitiesandtemperaturesintheWCR,radiativecoolingbecomesimportantandthushas
tobeconsidered,seeEq. (3). HereweusethecoolingfunctionΛ(T)fromSchureetal.(2009)whopresent
anewradiativecoolingcurvebasedonancontemporaryplasmaemissioncodeprovidingdetailedlogT-logΛ
tables. Weusethetabulateddataandinterpolatebetweenindividualdata-points.
–5–
2.4. Geometricalsetup&initialconditions
We use a Euclidean coordinate system defined such that the x-axis coincides with the semimajor axis
ofthebinarysystem,thusconnectingperiastronandapastron. Theoriginisatthecenterofmass. Thez-axis
isperpendiculartotheorbitalplane.
Initially we place the stars along the x-axis (choosing between periastron and apastron passage) and
initializethestellarwindsbyaβ-lawapproximation.
(cid:32) ξR (cid:33)βi
v(r) = v 1− ∗,i (8)
i ∞,i
r
FromthecontinuityequationEq. (1)itfollowsthat
M˙
ρ(r) = i (9)
i
4πr2v(r)
i
fortheindividualstars. ThetemperatureT isinitializedwiththevalue104K.Itisnotallowedtodropbelow
thisvalue,simulatingtheeffectsofphotoionizationheating. (cf.Pittard2009).
Asitiscomputationallyexpensivetoself-consistentlysimulatethewindupfromthestellarsurface,it
has been common practice to prescribe a fixed solution in at least 3 cells above the stellar surface at every
time step (e.g. Pittard 2009). A classical β-law solution (ξ = 1, only accurate in the point source limit)
usually results in an non-physical kink between fixed cells and those dynamically solved. Multiple trials
with1Dsimulationsclearlyshowthatasmoothtransitionisbestachievedbyinitiatingthefixedcellswith
ξ ∼ 0.9983–whichisalsothevaluesuggestedbyLamers&Cassinelli(1999). Toachievesmootherwind-
velocity profiles even at lower spatial resolutions we prescribe the wind with this modified β-law at every
time step up to 5 cells above the stellar surface. For a typical resolution of ∼15.6 R per cell, this gives
(cid:12)
satisfactoryvelocityanddensityprofiles.
For given stellar parameters, we determine the best fit value of β, as well as the values of the CAK-
parametersαandkbyaniterativesequenceof1DsimulationsbasedontheknownobservablesR , L , M ,
∗ ∗ ∗
thewind’sterminalvelocityv andthemasslossrate M˙.
∞
To simulate the orbital motion of the stars, we use the standard equations for a Keplerian orbit. The
differential equation for the eccentric anomaly Ψ(t) (being ωt = Ψ−esinΨ with eccentricity e) is solved
by a ten step Newton-Raphson method. The wind velocity that is prescribed at each step above the stellar
surfaceismodifiedbythestellarvelocityduetoitsorbitalmotion.
–6–
3. PARTICLESPECTRA
3.1. Thetransportequation
IninjectingelectronsandprotonsinsidetheWCRwhichthensufferfromvariousenergy-lossmecha-
nisms,wewidelyrelyontheworkbyReimeretal.(2006)inwhichananalyticalapproachforasimplified
windset-upiscarriedoutingreatdetail.
We distinguish between acceleration cells which are located at the shock front of the WCR, and all
othercells. Inordertodiscriminatebetweenthesetworegions,weusethetemperaturestructureofthewind
plasma. Due to the imposed lower limit of 104 K, the wind is isothermal until it reaches the WCR. There,
the temperature increases by three to four orders of magnitude within a few grid cells (see Fig. 1(b)). In
our approach, we declare a cell to be an acceleration cell if three conditions are met: i) the divergence of
thevelocityvectorbenegative(thewindiseffectivelyslowingdown),ii)thetemperaturebehigherthanin
theunshockedwind, andiii)thetemperatureofatleastoneofthesixneighbouringcellsbeatthevalueof
the unshocked wind. This yields a one-cell-thick skin of acceleration cells that envelop the WCR (see Fig.
1(d)).
Thetime-dependenttransportequation-asitissolvedforeachgridpoint-reads:
∂N(E) ∂ (cid:104) (cid:105) N(E)
+∇·[(cid:126)vN(E)]+ E˙N(E) + = Q δ(E−E ) (10)
0 0
∂t ∂E τ
(1) (4)
(2) (3)
where N is the differential number density of particles at energy E in a grid cell at position (cid:126)r. We now
discusseachtermindetail:
(1) Thespatialconvectionterm(where(cid:126)visthevelocityvectorofthewindmaterial)handlesthetransport
ofchargedparticlesdownstreamalongtheWCR.Itisusedforallgridcellsexceptaccelerationcells
whichistreatedinleakyboxapproachanalogoustotheaccelerationregioninReimeretal.(2006).
(2) Thistermhandlesallenergygainandlossprocesseswhichdependontheconsideredparticlespecies
andonwhetherthecellislocatedattheshockfrontornot. Foraccelerationcellstheterm E˙ includes
DSA (not active outside the shock front) and radiative losses. We also consider adiabatic cooling
whichbecomesimportantastheparticlesacceleratetravellingdownstream. Foraccelerationcellsand
all cells where ∇·(cid:126)v < 0 (true for cells just behind the shock front) the adiabatic cooling term has to
be switched off as it would allow additional acceleration which is explicitly already taken care of by
DSA.

E˙ = EEE˙˙˙rDaSdiAati+veE+˙raEd˙iative feeollssreeawwchhceeelrreeeriiaffti∇∇on··(cid:126)(cid:126)vvce<>lls00 (11)
radiative adiabatic
TheindividualenergygainandenergylosstermsarediscussedindetailinSections3.2and3.3.
(3) By considering the acceleration cells similarly to a leaky box model (see e.g. Protheroe & Stanev
1999), the escape time τ describes the rate at which particles are lost by diffusing out of the system.
–7–
Theescapetimecanbeapproximatedviathediffusioncoefficientandthewindvelocityperpendicular
totheshockV . ForWCRcellsoutsidetheshockfront,thediffusionalleakageoutofthesystem
Shock
issettozerobychoosinganinfiniteescapetime(seeMartin&Dubus2013).

τ = ∞VcS2r,hDock, elsefworhaecrecelerationcells (12)
withthecompressionratioc andtheenergy-independentdiffusioncoefficientD,whichisanapprox-
r
imation for the sum of the upstream and downstream component of the diffusion coefficients at the
shock (D ≈ D + c D ) (e.g., Schure et al. 2010). Kirk et al. (1998) provide a comparison to the
1 r 2
alternative assumption of an energy-dependent diffusion coefficient. In the present work, we choose
theenergy-independentapproachinordertoallowdirectcomparisonwithReimeretal.(2006). Thus,
D is constant throughout this work. Note that an energy-independent diffusion coefficient demands
additional consideration of the Bohm limit, at which further acceleration of particles is inhibited as
their gyroradii become comparable to the characteristic size of the shock. We include this in our
simulations by setting τ = 0 as soon as the Bohm diffusion coefficient (∝ E) exceeds the chosen
energy-independentdiffusioncoefficientD.
τ = 0, forE > E = 3eBD (13)
Bohm
wherethemagneticfieldstrength Bisapproximatedasdescribedbelow.
(4) Particles are injected into the system at an energy E (usually chosen to be 1 MeV). The choice of
0
theinjectionrateQ islimitedtotheconstraintsofparticlenumberandenergyconservation. Agiven
0
gridcellintheshockregioncannotinjectmoreparticlesthanitinitiallycarries,neithercanitdeposit
more energy in the injectedparticles than it has. Following Martin& Dubus (2013) we approximate
theamountofinjectedparticlesbyaconstantfractionηofthemassdensity.
Assuming a wind that predominantly consists of ionized hydrogen, partly ionized helium and elec-
trons,itfollowsthat
ρ 1+I ζ ρ
n = He He and n = (14)
e m 1+4ζ p m (1+4ζ )
H He H He
with n and n being the number densities of free electrons and protons in a given cell, ρ being the
p e
local wind density, ζ = nHe, and I being the number of electrons provided per helium nucleus.
He nH He
Hydrogen is assumed to be completely ionized. (We assume ζ = 0.1 and I = 2 for wind tem-
He He
peratures higher than 30 000 K and I = 1 below.) This imposes a limit on the maximum number
He
densityforelectronsQe andprotonsQp thatcanbeacceleratedpertimestepdt,ofwhichwetakethe
0 0
fractionsη andη . OtherparticlespeciesliableforsignificantaccelerationintheWCR(asHeions)
e p
arenotconsideredinthepresentstudy.

Q = ηdetmρH 11++IH4ζeζHHee or ηdpt mH(1+ρ4ζHe), foraccelerationcells(electrons,protons) (15)
0 0, elsewhere
–8–
The injection fractions are used as free parameters. Typically, we take η = 10−3 and η = 10−5.
p e
The ensuing mass density decrease in the wind plasma is small enough to be negligible. As the
energydepositedintheshockisalsoproportionaltoη (andthusverysmall),thereisnosignificant
e,p
reductionofthekineticenergyofthewind. Energyconservationthereforeholds. Wedonotconsider
alternativesourcesofparticlesenteringtheaccelerationprocess(e.g.,viaγ−γpairproduction.)
3.2. Particleacceleration
Ofthevariousaccelerationmechanismsfornon-thermalparticlesthatarediscussedinliterature,DSA
(as a variant of the first-order Fermi acceleration principle in which particles gain energy by repeatedly
traversing a shock) appears to be the most feasible process to be taken into account for our models. Al-
ternative mechanisms as various turbulent processes (e.g., the second-order Fermi processes) or magnetic
reconnection are assumed to be either less efficient or incapable of supplying particles that reach energies
sufficientforgamma-rayemission(seePittard&Dougherty2006). Thus, werestrictourselvestoconsider
DSAasthesoleaccelerationprocessofelectronsandprotons,fortheremainderofthiswork. Theaverage
rateofmomentumgainbyDSA(seee.g.;Schureetal.2010)is
(cid:32)c −1(cid:33) V2
E˙ = r ShockE (16)
DSA
3c D
r
whereV istheshockvelocity,c thecompressionratioand Dtheenergy-independentspatialdiffusion
Shock r
coefficientasmotivatedabove.
We define the shock velocity V as the upstream velocity normal to the shock and determine it by
Shock
computingthevelocitycomponentoftheplasmaperpendiculartotheWCR.Asatracerfortheorientation
ofthecollisionregion,wetakethegradientofthetemperaturefield∇(cid:126)T andcompute
(cid:126)v(cid:5)∇(cid:126)T
V = (17)
Shock
| ∇(cid:126)T |
Thecompressionratioc attheshockisalsodirectlydeterminedfromthewindstructurebycomputingthe
r
ratioofpost-shockmassdensityandpre-shockmassdensity. Theformerisdeterminedbyinterpolatingthe
massdensityatadistanceofthreecellwidthsalongtheshocknormaltowardstheWCR.Thechoiceofthree
cell widths is motivated by the fact that the shock front in the simulations is in general three cells wide. A
choice of two cell widths would yield too low postshock mass densities and result in a lower compression
ratio. Likewise, a choice of four or more cell widths probes the mass density too far inside the WCR. As
themassdensityminimumofthewindplasmapriortoreachingtheWCRisgenerallylocatedattheshock
positionforwhichthecompressionratioiscomputed,wesetthepreshockmassdensitytothelocalvalue.
Onceaccelerated,theelectronsandprotonsaresubjectedtoseverallossmechanismwhichwediscuss
below.
–9–
3.3. Energylosses
Inverse Compton emission (electrons only) IC cooling occurs as relativistic electrons scatter on stellar
radiation fields within the binary system. It is a major energy-loss mechanism for non-thermal electrons.
HereweusethefullKlein-Nishinacrosssectionresultinginthelossterm:
E˙ = −b E2f (E) (18)
IC IC KN
with
4
b = σ u (19)
IC Th ph
3m2c3
e
where σ is the Thomson cross section and f (E) is the correction factor for the Klein-Nishina regime
Th KN
(followingModerskietal.2005)being
(cid:34)(cid:32) (cid:33) (cid:32) (cid:33) (cid:35)
9 1 6 11 1
f (b˜) = b˜ +6+ ln(1+b˜)− b˜3+6b˜2+9b˜ +4 −2+2Li (−b˜) (20)
KN b˜3 2 b˜ 12 (1+b˜)2 2
with the dilogarithm-function Li and the dimensionless variable b˜ = 4(cid:15) E with (cid:15) being the energy of
2 Tm2c4 T
e
thetargetphoton. Forb˜ < 10−3 itissafetoset f = 1.
KN
Theradiationenergydensityu (takingintoaccountbothstars)is
ph
 
uph = 41πc Lr∗2,1 + Lr∗2,2 (21)
1 2
where L are the stellar luminosities and r the distance to the stars. At present, we do not consider addi-
∗,i i
tionalradiationfields(e.g.,photonsfromsynchrotronemission)astargetsforICscattering.
Synchrotronemission(electronsonly) Thelosstermforsynchrotronemissionis
E˙ = −b E2 (22)
syn syn
with
4
b = σ u (23)
syn Th B
3m2c3
e
wherethemagneticenergydensityu isdefinedas
B
B2
u = (24)
B
2µ
0
Duetothelackofknowledgeofthemagneticfieldstrength Binmostcolliding-windbinaries,werelyhere
ontheapproximationsinUsov&Melrose(1992)whodescribethemagneticfieldaseitheraclassicaldipole
field,aradiallydominatedfieldandoratoroidallydominatedfield,dependingonthedistancefromthestar.
B≈ BBB∗∗∗(cid:18)(cid:18)(cid:16)RrvvrAR∞r∗or3∗rt(cid:17)2RA3(cid:19)2∗r,(cid:19),, fffooorrrrrr <>> rrrAAA aannddrr <> RR∗∗vvvvrr∞∞oott (25)
–10–
where B isthemagneticfieldatthestellarsurfaceandv isthesurfacerotationvelocityofthestar, typi-
∗ rot
callyapproximatedbyv ∼ 0.1v . FortheAlfve´nradiusr wetake
rot ∞ A

R∗(1+ξ), forξ (cid:28) 1
r ≈ (26)
A R∗ξ1/4, forξ (cid:29) 1
withξ = B2∗R2∗ . Asareasonablevalueforthesurfacemagneticfield B = 0.01Tisassumedthroughoutthis
M˙v∞ ∗
workforbothstars(seee.g.,Reimeretal.2006).
Thus, we have an estimate for the magnetic field depending on the distance from each star. Since
we merely approximate the absolute value of the magnetic field vector without having knowledge of its
individual components, we cannot compute its vector sum for both stars. As an approximation, we merely
considerthedominantcomponentatagivengridcelltocalculateu andb .
B syn
InordertopreserveconsistencywithReimeretal.(2006),wedonotacknowledgethecompressionof
themagneticfieldinproportiontothecompressionofthewindinsidetheWCR.Thus, thefieldstrength B
islikelytobeunderestimatedbyafactorof∼4whichhasnoeffectatlowerenergies,butcaninfluencethe
maximumparticleenergyduetohighersynchrotronlossesandahigherBohm-energy. Infuturework,afull
MHDdescriptionofthewindwillallowgreaterprecisioninthetreatmentofthemagneticfield.
Losses by thermal bremsstrahlung (electrons only) As charged particles interact with the Coulomb
fields of ions in the wind plasma, bremsstrahlung emission occurs. To account for the ensuing energy loss
weuse
E˙ = −b E (27)
br br
with
2
b = ασ cN (28)
br Th H
π
where α is the fine-structure constant and N is the number density of the thermal ions in the wind which
H
weapproximatefromthemassdensityinagivengridcell–dividingitbym (1+4ζ ).
H He
ρ
N ≈ (29)
H m (1+4ζ )
H He
Coulomblosses(electrons) Aslosstermfortheelectronsweconsidertheenergy-independentCoulomb
lossesgivenby
E˙ = −b = −55.725cσ N m c2 (30)
coul coul Th H e
Coulomb losses (protons) In the dense winds Coulomb losses can also become significant for protons.
Theyareexpressedby
3cσ m c2Z2lnλ β2
E˙ = − Th e N (31)
coul 2 Hx3 +β3
m