Table Of Content

MHD Models of Gamma-ray Emission in WR 11 K. Reitberger1,a), R. Kissmann1, A. Reimer2and O. Reimer1 1Institutfu¨rAstro-undTeilchenphysik,Leopold-Franzens-Universita¨tInnsbruck,A-6020Innsbruck,Austria 7 1 2Institutfu¨rTheoretischePhysik,Leopold-Franzens-Universita¨tInnsbruck,A-6020Innsbruck,Austria 0 2 a)Correspondingauthor:[email protected] n a Abstract.Recentreportsclaimingtentativeassociationofthemassivestarbinarysystemγ2Velorum(WR11)withahigh-energy J γ-raysourceobservedbyFermi-LATcontrasttheso-farexclusiveroleofηCarinaeasthehithertoonlydetectedγ-rayemitterinthe 5 sourceclassofparticle-acceleratingcolliding-windbinarysystems.Weaimtoshedlightonthisclaimofassociationbyproviding 2 dedicatedmodelpredictionsforthenonthermalphotonemissionspectrumofWR11. Weusethree-dimensionalmagneto-hydrodynamicmodelingtotracethestructureandconditionsofthewind-collisionregionof ] WR11throughoutits78.5dayorbit,includingtheimportanteffectofradiativebrakinginthestellarwinds.Atransportequation E isthensolvedinthewind-collisionregiontodeterminethepopulationofrelativisticelectronsandprotonswhicharesubsequently H usedtocomputenonthermalphotonemissioncomponents. . Wefindthat–ifWR11beindeedconfirmedastheresponsibleobjectfortheobservedγ-rayemission–itsradiationwillunavoid- h ablybeofhadronicoriginowingtothestrongradiationfieldsinthebinarysystemwhichinhibittheaccelerationofelectronsto p energiessufficientlyhighforobservableinverseComptonradiation.Differentconditionsinwind-collisionregionneartheapastron - o andperiastronconfigurationleadtosignificantvariabilityonorbitaltimescales.Thebulkofthehadronicγ-rayemissionoriginates r ata∼400R(cid:12)wideregionattheapex. t s a [ INTRODUCTION 1 v Thesurprisingdiscrepancybetweenthepredictionofhigh-energyγ-rayemissionfromcollidingwindbinaries[1,2,3] 4 andtheirongoingnon-detection[4]isstillunresolved.Thehithertoonlyexceptiontotherule,thebrightγ-raysource 8 ηCarinae[5,6],rendersthelackofemissionfromsimilarsourcesevenmoreinteresting. 2 Therecentlyreportedanalysisof7promisingcolliding-windbinarysystemswith7yearsofdatafromtheFermi- 7 LATbyPshirkov[7]evenaddstotheproblembyfurtherloweringtheupperfluxlimitofWR140–oncepredictedto 0 beapromisingγ-source–tolessthan1.1×10−9 phcm−2s−1 inthe0.1−100GeVenergyrange.Therespectivevalue . 1 forηCarinaeis184±30×10−9phcm−2s−1[8].Takingintoaccountthatthetwosystemshaveasimilardistancefrom 0 Earth,similareccentricitiesandstellarseparations,andeventhatthetotalkineticwindenergyavailableforparticle 7 accelerationandultimatelyforγ-rayemissioninWR140isafactorof2largerthanforηCarinae[9]-howcanthis 1 : differenceinfluxofmorethantwoordersofmagnitudebeexplained?Theanswermightbeconnectedwiththemany v peculiaritiesoftheηCarinaesystem.Itissingledoutamongstothercolliding-windbinarysystemsbythenatureof i X itsprimarystar,presumablyaLuminousBlueVariablewithastellarwindthatis 1.5ordersofmagnitudedenserthen theprimarywindofe.g.WR140. r a Pshirkov [7] now reports the detection of WR 11, a short-period WC8 + O7.5 binary system very unlike the systemsdiscussedabove[seee.g.,9].However,itsclosenesstoEarthwithadistanceofmerely340pc[10]andthe knowledge we have of its stellar and stellar wind parameters make WR 11 a prime target for numerical modeling. Insightsgainedbyreproducingtheobserveddatawithourmodelsmaybeusedtofurtherconstrainmodelparameters neededforthemodelingofmorecomplexsystems.Theabilitytosuccessfullymodeltheobservedγ-rayemissionwill supporttheclaimofassociation. Byperformingthree-dimensionalmagneto-hydrodynamic(MHD)simulationofWR11,WR140,andηCarinae, we strive to obtain a better understanding of the physics of colliding-wind binary systems, ultimately seeking an explanation for the apparent lack of γ-ray emission in some of these systems. In this work we present our model resultsforWR11. CHALLENGESOFMODELLINGCOLLIDING-WINDBINARIES Modeloverview Following a procedure similar to the one detailed in Reitberger et al. [11] and Kissmann et al. [12], we use the Cronos-code [13] to perform a threedimensional MHD simulation of the colliding-wind binary system. The stellar windsareimplementedusingamodifiedCastor-Abbott-Klein(CAK)approximation[basedon14].Inadditionalto theeightMHDvariables(beingdensityρ,velocity(cid:126)v,temperatureT,andmagneticfieldstrength B(cid:126)),weinclude200 additionalscalarfieldsrepresentingelectronsandprotonsat100energybinsintherangeof1MeVto10TeV.The transportequationissolvedinruntimealongwiththeMHDequationsandincludestherelevantenergylossandgain mechanismsfortheparticles,aswellasenergy-dependentdiffusion.Simulationsalsoconsidertheorbitalmotionof thestars. This method provides spectra of high-energy electrons and protons for every region of the simulated volume. FollowingtheproceduredetailedinReitbergeretal.[15]thisdistributionofparticlesisthenusedtocomputeleptonic andhadroniccomponentsofhigh-energyγ-rayemissiondependingontheparticles’interactionwiththesurrounding radiationfields,magneticfields,andwindplasmadensity. Moredetailsonthecurrentmodeofoperationandtheperformanceofourcode–asitisusedtomodelWR11 andothercolliding-windbinarysystems–willbegiveninasubsequentpublicationonmodelresultsfortheWR11 system.Here,werestrictourselvestolistanumberofmajorchallengesinvolvedinthesimulation. Radiativebraking Radiativebrakingoccurswhenthestellarwindinabinarysystemisefficientlysloweddownbytheradiativeinfluence ofthesecondstarbeforeitreachesthecollisionzone.ThereisevidencethatthiseffectisimportantinWR11[16],as wellasduringtheperiastronpassageofηCarinae[e.g.,17].Theproperimplementationofradiativebrakingisnon- trivialasthecomputationoftheforcecausedbyonestar’sradiationfieldonthesecondstar’swindisambiguously treatedinliterature.Therearethetwopossibilitiesofeitherusingstar-specificCAK-parameters(weakcoupling)or wind-specific CAK-parameters (strong coupling). Numerical simulations of CWB-systems have hitherto predomi- nantly used the former method [e.g., 18, 19, 20, 11], whereas theoretical papers on radiative inhibition and braking havegenerallyappliedthelattermethod[e.g.,21,22,23].WefindthatinthecaseofWR11sufficientagreementwith observationsfromX-rayspectroscopycanonlybeachievedbyusingstrongcoupling.Thereforeweconsistentlyuse thismethod. Themagneticfield Care must be taken in the implementation of the magnetic field. The effect that a strong dipole field has on the wind plasma of a single has star has been studied by ud-Doula and Owocki [24] in great detail. Its impact on the shape of and conditions at the wind collision region in a binary system were recently investigated by Kissmann et al. [12]. It is shown that the wind-collision region can be significantly distorted by a strong magnetic field (e.g., B =100G).Windvelocitiestowardsthepolesofthemagneticfieldaregreatlyincreased.Inordernottoencounter surface theaforementioneddistortions,wechoosearelativelylowsurfacemagneticfieldof10G.Althoughsuchafielddoes nothaveanysignificantinfluenceonthewindstructure,wefindthatitdoesleaveasignificantimprintontheelectron distribution.Wherethepolesofthemagneticdipoleareclosetothewind-collisionregion,synchrotronlossesarehigh andelectronsreachlowerenergiesthanalongthemagneticequatorwherethefieldisweak. Energy-dependentdiffusion WhereasmaximumelectronenergiesaredeterminedbysynchrotronandinverseComptonlosses,whichusuallydom- inate the acceleration at higher energies, the maximum energies of the protons are determined by diffusion. An ap- proximationoftenusedistheBohmcutoff,whichsimplyimposesamaximumenergyatwhichtheprotonsleavethe shockfrontsbeforebeingfurtheraccelerated.Suchanapproximationisnolongernecessaryifanenergy-dependent diffusioncoefficientisused[25].Inourmodel,weapplyadiffusioncoefficientoftheformD(E)= D Eδ.According 0 toliterature,theexponentδisgenerallysettoδ=0.3foraKolmogorovtypespectrumandtoδ=0.5foraKraichnan typeturbulencespectrum[26].Wetrybothandcompareittothedata.D remainsafree-parameter. 0 600 (cid:65)(cid:65)(cid:65)(cid:65) 6060003e+06 (cid:3)(cid:3)(cid:3)(cid:3) 3e3+e0+606 (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) 400 404000 (cid:65)(cid:65)(cid:65)(cid:65) 2.4e+06 (cid:3)(cid:3)(cid:3)(cid:3) 2.42e.4+e0+606 (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) 200 202000 (cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3) (cid:65)(cid:65)(cid:65)(cid:65) 1.8e+06 (cid:3)(cid:3)(cid:3)(cid:3) 1.81e.8+e0+606 ](cid:12) θ(cid:65)(cid:65)(cid:65)(cid:65) (cid:3)(cid:3)(cid:3)(cid:3)θ [R 0 00 z 1.2e+06 1.21e.2+e0+606 200 202000 6e+05 6e6+e0+505 400 404000 600 6060000 0v0[ms−1] 800 600 400 200 0 200 400 404000 202000 00 202000 404000 606000 808000 x[R(cid:12)] x[R(cid:12)] FIGURE1.Absolutevelocityofwindplasmaforapastron(left)andperiastron(right)configuration.Theapproximateopening angleθisindicatedbytheyellowline. Orbitalmotion AlthoughWR11hasaneccentricityof0.3whichiscomparablylowtothevalueof∼0.9forηCarinaeandWR140,the conditionsatperiastronandapastronarefarfromsimilar.Whereastheflowdownstreamoftheshockisfairlylaminar aroundapastron,turbulenceemergesasthesystemclosesinonitsperiastronpassage.Ithasalsobeenstudiedtowhat degreethedistortionofthewind-collisionregioncausedbytheorbitalmotioninfluencestheparticledistributionand γ-rayemission.Wefindnosignificantdifferenceinthefinalγ-rayemissionofasimulationthatfullyconsidersorbital motioncomparedtoasimulationshowingthesteadystateatap-orperiastron. RESULTSOFMODELLINGWR11 Shapeofwindcollisionregioninagreementwithobservations In their analysis of high-resolution X-ray spectroscopy of Chandra data, Henley, Stevens, and Pittard [16] find evi- denceforalargeshock-coneopeningangle.Theyremarkthatsimulationshithertofailedtoreproducesuchafeature whichmightbeduetosignificantradiativebraking.Inoursimulation,includingtheeffectsofradiativebrakingwith strongcouplingoftheCAKparameters,wecanindeedreproducethislarge-shockconeopeningangle,whichismost evidentclosetotheperiastronpassage. Figure1showsabsolutevelocityofthewindplasmafortheapastronandperiastronconfiguration.Theopening angleisalsoindicated.Thestrongeffectofradiativebrakingisclearlyvisiblebytheblueareaupstreamofthewind collision where the secondary wind is effectively slowed down by the primary star’s radiation. Also, the periastron plot displays the effect of shadowing where the secondary wind is accelerated less in the region where the primary star is eclipsed by the disk of the secondary. As mentioned above, the shock is fairly laminar at apastron and more turbulentatperiastron.Theshock-coneopeningangleis∼65◦forapastronand∼76◦forperiastron. Hadronicdominanceinhigh-energyγ-rayemission Inashort-periodbinarysystemlikeWR11,thestellarseparationandthereforethedistancebetweenthestarsandthe windcollisionregionarefairlysmall.Stellarradiationfieldsandmagneticfieldsatthesiteoftheparticleacceleration are therefore high compared to long-period binary systems. This leads to severe energy losses by inverse Compton emissionandsynchrotronemission.WefindthattheelectronsinWR11donotreachenergieshigherthan100MeV. ThisisshownintheleftandcenterplotofFigure2whichshowsthemaximumenergiesoftheelectronsinthe x−yand x−zplane.Theasymmetryisduetothemagneticdipolefieldwhichisalignedalongthez-axis.Electrons 606000 1e1+e0+202606000 1e1+e0+202606000 1e1+e0+505 404000 404000 404000 202000 202000 202000 R](cid:12) R](cid:12) R](cid:12) [00 1010 [00 1010 [00 1e1+e0+404 y z z 202000 202000 202000 404000 404000 404000 606008008000 606000 404000 x2[02R000(cid:12)] 00 202000 40400E0max11[Me6V060080]08000 606000 404000 x2[02R000(cid:12)] 00 202000 40400E0max11[Me6V060080]08000 606000 404000 x2[02R000(cid:12)] 00 202000 404E000max1[e1+Me0+30e3V] FIGURE2.Maximumparticleenergiesforelectronsinthey−zplane(left),electronsinthe x−zplane(middle)andprotons (right). 0.002.0020 1e-10e9-09 0.001.5015 1e-11e0-10 0.001.0010 1e-11e1-11 0.000.5005 1e-11e2-12 0.005arcsec 0.000.0000 1e-11e3-13 0.000.0000 0.000.5005 0.001.0010 0.001.5015 0.002.00F2l0ux[m−2s−1] FIGURE3.Left:Projectedfluxabove100MeVforneutralpiondecay.Right:Schematicviewofthetwostarswithinthecompu- tationaldomain.Thelineofcentersisrepresentedbythehorizontaldashedorangeline.Variousviewingangles(lines-of-sight)are indicated,includingtheonethatisusedinoursimulation:i=65◦,Φ=68◦. reachhigherenergiesalongtheequatorofthedipolewheresynchrotronlossesarelow.Highestenergiesarereached intheouterwingsofthethecollisionregionwherethedistanceofthestarsincreases. Asprotonsarenotaffectedfromsimilarlossesrelatedtothemagneticandradiationfields,theyreachenergiesup to100GeV.Attheseenergiestheenergy-dependentdiffusionisstrongandallowstheacceleratedparticlestotravel upstreamoutsidethecollisionregion.Thehighestenergyprotonsarefoundattheapexofthewindcollisionregion. ThisisshownintherightplotofFigure2. Duetothelowmaximumenergyofelectrons,WR11isaclearcaseofhadronicdominanceinthehigh-energy γ-rayemission. Confinementofemissionregiontoapexofwind-collision The inclination i of the orbital plan of the γ2 Velorum binary system and its argument of apastron ω are well WR constrained. Schmutz et al. [27] find i = 65◦ ± 8◦ and ω = 68◦ ± 4◦. Using this information as well as the WR distributionofhigh-energyprotonsandthewind-plasmadensitywecomputeaprojectedfluxmapasshownintheleft part of Figure 3. The right part illustrates the orientation of the orbital plane and the stars at apastron configuration withregardtothelineofsight.Theprojectedfluxmapisshownasitwouldbeseenalongtheredarrowatadistance of d =342 pc. It indicates that most of the γ-ray emission by neutral pion decay has its origin at the apex of the wind-collisioninaregionroughly5milliarcsecor 400R wide. (cid:12) ) ) ])1 −1s]1 −1s]0 −3m2 −2m2 −2m1 V V V2 e3 e e [M [M3 [M3 N4 F F 2E 2E4 2E4 (5 ( ( g g g5 o o o l l5 l 6 6 70 1 log2(Energ3y[Me4V]) 5 6 64 2 log(E0nergy[M2 eV]) 4 6 70 1 log2(Ener3gy[Me4V]) 5 6 FIGURE4.Left:spectralenergydistribution(SED)forneutralpiondecayatperiastron(redsolid)andapastron(greendashed). Middle:SEDforinverseComptonemission(redsolid),bremsstrahlung(greendashed),andneutralpiondecay(bluedotted)at apastron.Right:SEDforneutralpiondecaywithadiffusionexponentδ = 0.3(redsolid)andδ = 0.5(greendashed)alongwith thedatapointsasofPshirkov[7]. Variabilityonorbitaltimescales ThecomparisonoftheabsolutevelocitiesinFigure1alreadysuggeststhatthedifferentconditionsatthewindcollision region for periastron and apastron configuration will influence particle acceleration and γ-ray emission. Indeed, we find that the maximum γ-ray energies as well as the maximum differential flux level differ by almost 2 orders of magnitude.ThisisshownintheleftplotofFigure4.Apastronconditionsareclearlymorefavourableforhigh-energy γ-rayemission.However,theverylowstatisticswillmakesuchanexpectedvariabilitydifficulttoobserve. Broadbandspectra Althoughelectronenergiesarefartoolowtoproducehigh-energyγ-rayemissionwecanstillusetheirdistribution to compute nonthermal photon emission from inverse Compton scattering and bremsstrahlung emission at lower energies.Thuswecomputebroadbandspectra,asshowninthecenterplotofFigure4,whichcanalsohelptoconstrain importantparameters,e.g.limitsonthemagneticfieldstrengths. Fittingthediffusioncoefficient As mentioned in the above discussion of energy-dependent diffusion, the diffusion index δ can be assumed to have eitherthevalue0.3(Kolmogorov)or0.5(Kraichnan).WeapplyourmodeltothedatabyPshirkov[7]andfindthat thedataagreesfarbetterwithmodellingresultsforδ = 0.3whichisexpectedforathreedimensionalmodel.Model resultsforthetwovaluesofδ-allotherparametersequal-areshownintherightplotofFigure4alongwiththedata. SUMMARYANDOUTLOOK Weareusinganimprovedversionofapreviouslypresented[11,15,12]MHDcodetosimulatethehigh-energyγ-ray emissionofmassive-starbinarysystemswithcollidingwinds.Threesystems,WR11,ηCarinae,andWR140,areof particularinterestastheparametersnecessaryformodellingaresufficientlyconstrainedandobservationseitheryield detection or suprisingly low uper limits. In this work, we present our results for WR 11, showing a) agreement of thewindstructurewithobservationsfromX-rayspectroscopy,b)hadronicdominanceofhigh-energyγ-rayemission, c) a well confined emission region close to the apex of the collision region, and d) expected variability in apastron to periastron flux. Our model results support the claim of association of the observed γ-ray source with the WR 11 system. InthenearfuturewewillpresentourresultsonηCarinaeandWR140,focusingonfindinganexplanationfor the on-going non-detection of the WR 140 which might be connected to the peculiar role of the primary wind in ηCarinaeandtheinfluenceofwindcompositiononthesystem’scapabilityforefficientparticleacceleration. ACKNOWLEDGMENTS The computational results presented have been achieved (in part) using the HPC infrastructure of the University of Innsbruck.A.R.acknowledgesfinancialsupportfromtheAustrianScienceFund(FWF),projectP24926-N27. REFERENCES [1] P.BenagliaandG.E.Romero,A&A399,1121–1134,March(2003). [2] A.Reimer,M.Pohl, andO.Reimer,ApJ644,1118–1144,June(2006). [3] J.M.Pittardetal.,A&A446,1001–1019,February(2006). [4] M.Werneretal.,A&A555,p.A102,July(2013). [5] M.Tavanietal.,ApJLetters698,L142–L146,June(2009). [6] K.Reitberger,A.Reimer,O.Reimer, andH.Takahashi,A&A577,p.A100,May(2015). [7] M.S.Pshirkov,MNRAS457,L99–L102,March(2016). [8] F.Aceroetal.,ApJS218,p.23,June(2015). [9] M.DeBeckerandF.Raucq,A&A558,p.A28,October(2013). [10] J.R.North,P.G.Tuthill,W.J.Tango, andJ.Davis,MNRAS377,415–424,May(2007). [11] K.Reitberger,R.Kissmann,A.Reimer,O.Reimer, andG.Dubus,ApJ782,p.96,February(2014). [12] R.Kissmann,K.Reitberger,O.Reimer,A.Reimer, andE.Grimaldo,arXiv:1609.01130,acceptedforpub- licationinApJ,September(2016). [13] J.Kleimann,A.Kopp,H.Fichtner, andR.Grauer,ANNALESGEOPHYSICAE27,989–1004(2009). [14] A.Pauldrach,J.Puls, andR.P.Kudritzki,A&A164,86–100,August(1986). [15] K.Reitberger,R.Kissmann,A.Reimer, andO.Reimer,ApJ789,p.87,July(2014). [16] D.B.Henley,I.R.Stevens, andJ.M.Pittard,MNRAS356,1308–1326,February(2005). [17] E.R.Parkinetal.,MNRAS394,1758–1774,April(2009). [18] J.M.Pittard,MNRAS396,1743–1763,July(2009). [19] E.R.Parkin,J.M.Pittard,M.F.Corcoran, andK.Hamaguchi,ApJ726,p.105,January(2011). [20] T.I.Maduraetal.,MNRAS436,3820–3855,December(2013). [21] I.R.StevensandA.M.T.Pollock,MNRAS269,p.226,July(1994). [22] S.P.OwockiandK.G.Gayley,ApJ454,p.L145,December(1995). [23] K.G.Gayley,S.P.Owocki, andS.R.Cranmer,ApJ475,p.786,February(1997). [24] A.ud-DoulaandS.P.Owocki,ApJ576,413–428,September(2002). [25] J.G.Kirk,F.M.Rieger, andA.Mastichiadis,A&A333,452–458,May(1998). [26] A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin, Annual Review of Nuclear and Particle Science 57, 285–327,November(2007). [27] W.Schmutzetal.,A&A328,219–228,December(1997).