WIMP diffusion in the solar system including solar depletion and its effect on Earth capture rates Johan Lundberg∗ Division of High Energy Physics, Uppsala University, SE-751 21 Uppsala, Sweden Joakim Edsjo¨† Department of Physics, AlbaNova University Center, Stockholm University, SE-106 91 Stockholm, Sweden (Dated: February 2, 2008) Weakly Interacting Massive Particles (WIMPs) can be captured by the Earth, where they even- tuallysinktothecore,annihilateandproducee.g.neutrinosthatcanbesearchedforwithneutrino telescopes. The Earth is believed to capture WIMPs not dominantly from the Milky Way halo directly, but instead from a distribution of WIMPs that have diffused around in the solar system 4 due to gravitational interactions with the planets in the solar system. Recently, doubts have been 0 raised about thelifetime of theseWIMP orbitsdueto solar capture. Wehereinvestigate this issue 0 bydetailed numerical simulations. 2 Compared to earlier estimates, we find that the WIMP velocity distribution is significantly sup- n pressed below about 70km/swhich results ina suppression of thecapturerates mainly forheavier a WIMPs (above∼ 100 GeV). At 1 TeV and abovethe reduction is almost a factor of 10. J WeapplytheseresultstothecasewheretheWIMPisasupersymmetricneutralinoandfindthat, 8 withintheMinimalSupersymmetricStandardModel(MSSM),theannihilation rates, andthusthe neutrino fluxes, are reduced even more than the capture rates. At high masses (above ∼ 1 TeV), 1 thesuppression is almost two orders of magnitude. v This suppression will make the detection of neutrinos from heavy WIMP annihilations in the 3 Earth much harder compared to earlier estimates. 1 1 1 I. INTRODUCTION [4], Gould refined the calculations of Press and Spergel 0 for the Earthandderivedexactformulaefor the capture 4 rates. In a later paper [5], Gould pointed out that since 0 Thereismountingevidencethatamajorfractionofthe / the Earth is in the gravitational potential of the Sun, h matter inthe Universe isdark. The WMAP Experiment all WIMPs will have gained velocity when they reach o-p wgihveerseasΩaCDbMestisfitthvealrueelitchdaetn[1si]tΩyCoDfMcohl2d=da0r.k11m3±at0te.0r0i9n, tvheeryEsamrtahll.anHdowheevnecre,cGaoputuldrelaotferheraevalyizWedIM[6]Ptshwatoduuldebtoe r unitsofthecriticaldensityandhistheHubbleparameter gravitationalinteractionswith the other planets (mainly t in units of 100 kms−1Mpc−1. One of the main candi- s Jupiter, Venus and Earth), WIMPs will diffuse in the a dates for the dark matter is a Weakly Interacting Mas- solar system both between different bound orbits, but v: siveParticle(WIMP),ofwhichthe supersymmetricneu- also between unbound and bound orbits. Gould showed Xi tralino is a favorite candidate. There are many ongoing that the net result of this is that the velocity distribu- efforts trying to find these dark matter particles, either tion at the Earth will effectively be the same as if the r viadirectdetectionorviaindirectdetectionbydetecting a Earth was in free space. This approximation is widely their annihilation products. used today where one further assumes that the halo ve- Oneoftheproposedsearchstrategiesistosearchfora locitydistributionisGaussian(i.e.aMaxwell-Boltzmann fluxofhigh-energyneutrinosfromthecenteroftheEarth distribution). [2]. ThisideagoesbacktoPressandSpergel[3],whocal- However, Farinella et al. [7] later made simulations of culated the capture rate of heavy particles by the Sun. For the Earth, the idea is that WIMPs can scatter off a NearEarthAsteroids(NEAs)thathadbeenejectedfrom theasteroidbelt. Theyfoundthatmanyofthesehavelife nucleus in the Earth, lose enough energy to be gravita- tionallytrapped,eventuallysinktothecoreduetosubse- timesoflessthantwomillionyears. Afterthattimethey quentscatters,annihilateandproduceneutrinos. Dueto areeither throwninto the sun or thrownoutof the solar purelykinematicalreasons,thecapturerateintheEarth system. If this typical lifetime also applies to WIMPs, depends strongly on the mass and the velocity distribu- this would significantly reduce the number of WIMPs tionof the WIMPs. The heavierthe WIMP is,the lower bound in the solar system, as pointed out in Ref. [8]. the velocity needs to be to facilitate capture. In Ref. This in turn would reduce the expected capture and an- nihilationratesintheEarthandthusreducetheneutrino fluxes. InRef.[8],GouldandAlaminvestigatedwhatthe implications would be if bound WIMPs would actually ∗Electronicaddress: [email protected] be thrown into Sun. They investigated two scenarios: †Electronicaddress: [email protected] an ultra conservative scenario where all bound WIMPs 2 aredepletedandaconservative scenariowhereallbound particlesareincreasedwhen they approachthe potential WIMPs that do not have Jupiter-crossing orbits are de- of the Sun. This reduces capture substantially. On the pleted. In the ultra conservative view solar depletion is other hand, bound solar orbits are allowed. Gould real- assumed to be so efficient that no bound WIMPs exist, ized that particles scattered by the Earth could become whereas in the conservative view, Jupiter is assumed to boundto the solarsystem. This scatteringcanbe oftwo be faster at diffusing WIMPs into the solar system than kinds: gravitational scattering, which is elastic in that solardepletionisatthrowingthemintotheSun. Bothof the velocity with respect to the Earth is conserved or these views significantly reduce the neutrino fluxes from weak scattering off an atom, which can be inelastic and the Earth for heavier WIMPs. either lead to capture by the Earth or make the particle However, the truth probably lies somewhere between bound to the solar system. In this context, an equation the conservative view and the assumption that solar de- for estimating the timescales of weak and gravitational pletion is very inefficient, i.e. some WIMPs on bound scatteringwasdeveloped,followingtraditionsofO¨pik[9]. orbits in the inner solar system will survive, but solar Amongotherthings,Gouldconcludedthat,duetothe capture will diminish their numbers somewhat. The aim differences in total scattering cross section, the gravity of this paper is to investigate the effects of solar capture of the Earth is more effective in changing the orbits of on the distribution of WIMPs in the solar system and bound particles than is weak scattering. For capture by theimplicationthishasonexpectedneutrinofluxesfrom the Earth though, weak scattering is the only process the Earth. We will do this by numerical simulations of at work since gravitationalscattering leaves the velocity WIMPs in the solar system and by reanalyzing the pro- with respect to the Earth unchanged. cess of WIMP diffusion in the solar system. Finally, we In 1991, Gould [6] continued further, moving his at- willapplyourresultsto the casewherethe WIMPis the tentionto thegravitationaldiffusioncausedbythe other neutralino, whicharises naturally in Minimal Supersym- planets. Further, he considered the combined diffusion metric extensions of the Standard Model (MSSM). effect of Jupiter, Venus andthe Earthconcluding that it The layout of this paper is as follows. In section II, will makethe velocity distribution isotropicin the frame we will briefly review the history of WIMP capture cal- of the Earth. Based essentially on Liouville’s theorem, culationsforthe Earth. InsectionIII wewillgothrough this means that the phase space density of unbound and ourassumedhalomodelandtheroleofdiffusioninmore bound particles would be the same. Specifically, for the detail. InsectionIVwewillgothroughtheformalismfor most important parts of velocity space, this would hap- the diffusion caused by one planet and in section V we pen on time scales shorterthan the age of the Solar Sys- addthenewingredient,solardepletion. InsectionVIwe tem. Obviously, such a scenario would substantially en- presentournumericaltreatmentofthediffusionproblem. hance capture of heavy WIMPs by the Earth. Further, All of this will be put together with the dominant plan- heconcludedthatweakcaptureofWIMPstoboundsolar ets for diffusion insectionVII where our mainresults on orbitsisnegligible,andthatonemayusethe”free-space” the velocity distribution at the Earth are presented. In formulae derivedin Ref. [4] for capture, even though the theremainingsectionswewillinvestigatehowthisaffects Earth is deep within the potential well of the Sun. the capture and annihilation rates in the Earth and will Asmentionedintheintroduction,thecalculationstook present results on the expected neutrino-induced muon a new unexpected turn in 1999 when Gould and Alam fluxes in MSSM models in section IX. Finally, we will interpreted[8] the results ofFarinella etal.[7]. Farinella conclude in section X. et al. had numerically calculated the fates of about 50 asteroidsofwhichmostwereconsideredtobenearEarth asteroids (NEAs). They concluded that about a third II. CAPTURE OF WIMPS BY THE EARTH – of the considered asteroids will be ejected to hyperbolic HISTORICAL REMARKS orbits or, more importantly, driven into the sun in less than2millionyears. IftheresultsofFarinellaetal.were Capture of WIMPs by the Sun was first studied by applicable to general Earth crossing orbits of WIMPs, Press and Spergel 1985 [3]. Their calculations were the part of velocity space corresponding to bound solar approximate in nature, especially when applied to the orbits would be effectively empty, since the typical time Earth. This was refined in a series of papers by Gould scales at which such orbits are populated from unbound [4, 5, 6]. In 1987, Gould [4] derived the exact for- orbits are generally muchlonger [6]. The basic results of mulae needed to calculate the capture of WIMPs by a Farinella et al. were later confirmed by Gladman et al. spherically symmetric body. When applied to the cap- [10] and Migliorini et al. [11]. ture by the Sun and the Earth, his approach enhanced To investigate the role of Solar depletion, Gould and the capture by factors of 1.5–3 and 10–300 respectively, Alam [8] analytically investigatedthe difference between compared to the previous approximations by Press and the 1991 case of no solar depletion, and the other ex- Spergel [3]. treme, where there is no dark matter in solar system However, in 1988, Gould [5] refined the analysis, tak- bound orbits at the Earth. In the latter ultra conser- ingintoaccountthattheEarthiswellinsidethegravita- vative view, capture is heavily suppressed. For instance, tionalpotentialoftheSun. Thevelocitiesoftheincoming they found that WIMPs with masses above about 325 3 GeVcouldnotbecapturedbytheEarthatall. Theyalso distributionpassthroughthe solarsystem, the velocities considered a scenario where WIMPs in Jupiter–crossing areboostedandfocusedbytothegravitationalpotential. orbits were not effected by solar capture, the conserva- At the location of the Earth, the solar system escape tive view, which they found allows WIMPs up to about velocityis√2v 42km/s,wherewehaveusedthespeed ⊕≈ 630 GeV to be captured. (See section IIIB below for a of the Earth, v 29.8 km/s. Therefore the velocity at more detailed discussion of these cutoff masses.) Both thelocationoft⊕he≃Earth,w,is,accordingtoconservation ofthesescenariossignificantlysuppressthecapturerates of energy of WIMPs by the Earth and the question to ask is if the results of Farinella et al. can really be applied to all w2 =s2+2v2. (4) ⊕ Earth–crossing WIMP orbits? The orbits of asteroids When a spherically symmetric distribution such as ejected from the asteroid belt are, after all, rather spe- cial as they typically arise from resonances. It is thus Fs(s) is focused by a Coulomb potential such as that of the Sun, the following statement holds: [5] not necessarily so that these results apply to all bound WIMPs. We will in the coming sections go through the F (w)4πw2dw F (s)4πs2ds necessary steps to investigate this question in detail. w = s (5) w s This can be understood as Liouville’s theorem for the III. THE GALACTIC HALO MODEL AND spherically averagedphase space density, since CUTOFF MASSES ds w = F (w)=F (s). (6) A. The galactic halo model w s dw s ⇒ Since the velocity w of the halo particles is always at In order to make the calculations concrete, we use the least equal to the escape velocity, there will be a hole in Maxwell-Boltzmann model [12], where the local velocity velocity space so that distribution of WIMPs is Gaussian in the inertial frame oftheGalaxy. AtthelocationoftheSunthedistribution F (w)=0 when w <√2v . (7) is w ⊕ 3 e−v2/v0 3 Thisis importantsincecapturebythe Earthisverysen- fv(v)d v = π3/2v3 d v, (1) sitive to Fw(w) at low velocities. 0 The distribution F (w) can now be used to calculate w thedistributionasseenfromthemovingEarthwherethe where v0 = 23v¯ with v¯ being the three-dimensional ve- particle velocity is u=w+v locity disperqsion. We will here use the standard value ⊕ of v¯ = 270 km/s corresponding to v0 = 220 km/s. The Fu(u)=Fw(w)=Fw(u v ). (8) − ⊕ distribution is normalized such that This meansthatthe hole isshifted, sothatitis centered fv(v)4πv2dv =1. (2) aartowuonddi−mve⊕n.siTonhaislisslivciesuoafltizheedtbhyrefiegduirmee1nwsihonicahldvieslpolcaiytys Z space. The velocity distribution can be galileo transformed into the frame of the Sun: f (s), where s = v+v , s Sun and vSun=220 km/s, and averaged over all angles. In B. Cutoff masses when low velocity WIMPs are this special case of a Gaussian distribution the transfor- missing mation can be done in closed form[4]. As Gould have pointed out, the angle between the rotation axis of the Inthe absenceofWIMPsgravitationallyboundto the solar system and that of the galaxy is about 60◦ which solar system, the capture by the Earth is totally sup- makes the velocity distribution very close to spherically pressed for WIMP masses larger than a critical value. symmetric,if one considers averages over a galactic year To understand this, consider a particle approaching the 200 million years[5]. The distribution used is mirror ≈ Earthwithvelocityuatinfinitywithrespecttothegrav- symmetric in the galactic plane which means that the itational potential of the Earth. If it is to be captured, time of averageneed only be 100 million years. it must be scattered by an atom to a velocity less than The symbol F (s) will be used to denote the phase s the escape velocity vesc at the atom. By conservationof space number density energy and momentum, the particle must have a veloc- ρ ity less than (assuming iron to be the heaviest relevant χ F (s)= f (s), (3) s s element of the Earth) M χ whereMχ isthe WIMPmass,andρχ isthe WIMPmass ucut =2 MχMFevesc (9) per unit volume in the halo. When the particles of this M M pχ− Fe 4 Here, u is the speed (at infinity) of the approachingpar- m/s) 50 uuuuuuuuuunnnnnnnnnnbbbbbbbbbboooooooooouuuuuuuuuunnnnnnnnnnddddddddddoooooooooorrrrrrrrrrbbbbbbbbbbiiiiiiiiiittttttttttssssssssss ticle inthe frame of the EarthandMχ,cut is the highest k allowedmassofthe particle ifit is to be capturedby the h( 40 Earth. The escape velocity varies from 11.2 km/s at the t r a surface to 15.0 km/s at the center of the Earth (see sec- E 30 e tion VIIIA for more information about the Earthmodel h t 20 bbbbbbbbbboooooooooouuuuuuuuuunnnnnnnnnnddddddddddoooooooooorrrrrrrrrrbbbbbbbbbbiiiiiiiiiittttttttttssssssssss we use), and capture is thus easiest at the center where t a d the escape velocity is higher. Using vesc = 15.0 km/s, pee 10 θ weplotinFig.2the relationbetweenucut andthe cutoff S mass, M . 0 χ,cut 70 60 50 40 30 20 10 0 10 20 30 SpeedattheEarth(km/s) With Eq. (10), we can now relate to the cutoff masses intheconservativeandultraconservativeviewsbyGould FIG. 1: The ecliptic (φ = π/2) slice of the particle velocity and Alam [8]. In the ultra conservative view, we assume space in the frame of the Earth. The dotted curves show the thatonlyunboundhaloparticlesarecaptured. Halopar- velocityrelativetotheEarthandtheindicatedangleθ,isthe ticles cannot be slower than ucut = (√2 1)v 12.3 angle of the particle with respect to the direction of Earth’s km/s at,and inthe frame of, the Earth(th−is is⊕als≃oseen motion. The angle φ determines in which angle we cut the in Fig. 1). This gives a cutoff mass of about 410 GeV velocity sphere. φ = 0 is the south pole of the solar system over which capture by the Earth is impossible. This dif- andφ=π/2(asshownhere)isthesliceradiallyoutwardfrom fers fromthe value of325 GeV for the ultra conservative the Earth. The region inside the solid semicircle represents view in Gould and Alam [8]. The difference is because boundorbits. It’sradiusistheescapevelocityfromtheSolar system at the location of the Earth, but in the frame of the they used an average escape velocity of 13 km/s instead Sun. Inthesameway,theregionoutsidethedash-dottedline of the maximal one of 15 km/s that we have used in (an almost perfect semi-circle) corresponds to particles that Fig. 2. may reach Jupiter. By repeated close encounters with the In the conservative view, we assume that Jupiter– Earth,particlesmaydiffusealongthedottedcircles(actually crossing orbits are filled. This means that all orbits spheres) of constant velocity only, keeping u constant, but outside the dot-dashed curve and the dashed curve in allowing changes in θ and φ,as explained in thetext. Fig.1 arefilled. The lowestvelocityWIMPatthe Earth that is on a Jupiter–crossing orbit is in the lower right- 103 handendofthe dot-dashedcurveandithasavelocityof /s) iuncguti=n tvh⊕e( sa2m/(e1d−irre⊕c/tiroXn)−as1t)h≃e E8.a8rtkhm)./sr(and5is.2mrovis- m p X ≃ ⊕ the radius ofthe Jupiter orbit. This value ofucut givesa k ( cutoff mass of about 712GeV, whereasGould and Alam t102 [8] got a cutoff mass of about 630GeV. The difference is u c againduetothe differentescapevelocitiesused,butalso u forbidden a different velocity to reach Jupiter. We use the value y t ucut 8.8 km/s as indicated above, whereas they used i ≃ c an approximation for more general orbits than the one o l 101 giving the cutoff derived here. So, to conclude, in the e v conservativeandultraconservativeview,wecannotcap- allowed ff ture WIMPs heavier than about 410 GeV and 712 GeV o t respectively. This is in roughagreementwith the results u C of Gould and Alam [8]. 101001 102 103 104 105 If, on the other hand, the solar system is full of gravi- Wimp mass (GeV) tationallybounddarkmatter,thevelocitiescanbemuch lower. AsthelowestallowedvelocityoftheWIMPsucut FIG. 2: Cutoff velocity vcut against WIMP mass M. Only tends to zero, the mass limit Mχ,cut goes to infinity. Typically, most WIMPs in the Galaxy have velocities combinationsofM andvcut totheleftofthelineiskinemat- ically allowed (in the sense that they can lead to capture by much greater than those of Eq. (9), so only a small frac- theEarth). tion of the WIMPs are possible to capture. Aparticleincloseencounterwithaplanet,forinstance the Earth, may get gravitationally scattered into a new direction and a new velocity as seen from the frame of when it approaches the gravitational potential of the theSun. However,byconservationofenergy,thespeedu Earth. Solving for the WIMP mass M gives χ withrespecttotheframeoftheplanet isunchanged. This means that a particle at a particular place in velocity u2+2vesc (vesc+ u2+ve2sc) space may, by repeated close encounters with the Earth, Mχ,cut =MFe u2 . (10) diffuse to any locationonthe sphere ofconstantvelocity p 5 (with respect to the Earth), and nowhere else. Thelocationofaparticleatsuchaspherecanbespeci- fiedbytheanglesatwhichitpassestheEarth. Theangle θ ismeasuredbetweentheforwarddirectionoftheEarth andthe velocityvectorofthe particle,andφisthe angle R l of rotation around the forward direction of the Earth, with φ=0 at the North pole of the solar system. Figure 1 illustrates how spheres (and circles) of con- stantucrossthe limitofwhereparticleshaveboundand unbound orbits. This corresponds to the possibility of ∆ξ gravitationalcaptureandejectionfromthe solarsystem. A single planet can diffuse particles along spheres of constant velocity only. Therefore, it is clear (from e.g. FIG. 3: The angle of perihelion precession ∆ξ, as the orbit Fig. 1) that orbits with velocities u less than 12.3 km/s enters and leaves the disk of the Earth. In this example, the planeoftheorbitisnearlyperpendiculartotheeclipticplane. cannot be populated by the Earthalone. However,since the velocity spheres of different planet’s are not concen- tric (they need not even be spheres when the particles reach another planet), the combined effect may diffuse orbits. Thisismotivatedbythe factthatduetothe pre- particles down to Earth-crossing velocities u less than cessionofperihelion,allsuchorbitellipseswilleventually 12.3 km/s. This will be investigated in detail in section intersect the Earth ring. The small angle the perihelion VII. In order to do so, we must first understand how a sweepsout, asthe orbitellipse entersandleavesthe ring single planet affects the phase space distribution. is given by l ∆ξ tan∆ξ = tanΘ1 (11) ≈ R| | IV. GRAVITATIONAL DIFFUSION IN THE ONE PLANET CASE where Θ1 is the intersection angle between the WIMP ellipseandtheplaneperpendiculartothelocationvector of the Earth R. Since this happens four times during In this section, we investigate the details of what will eachperihelionrevolution,the meanprobabilityforsuch be called gravitational diffusion. We will develop tools a WIMP to intersect the ring of the Earth during each for detailed investigationof the bound orbit phase space WIMP year T is density,takingtheeffectsofsolardepletionintoaccount. χ We will here start by looking at diffusion effects from a 4∆ξ p = . (12) singleplanetonlyandwilltaketheEarthasanexample. Tχ 2π The exact same formalism is then used for Venus and The probability for t(cid:10)he W(cid:11)IMP to come into close en- Jupiter as well. counter with the Earth is therefore p times the cross In this section we assume that when a particle is in Tχ sectionσofsuchanevent,dividedbytheareaoverwhich Earth crossing orbit (perihelion less than the Earth or- theEarthisdistributed. However,thelengthofthepath bit radius R and aphelion greater than R ), long range which is inside the Earth ring during each encounter is ⊕ ⊕ interactions with other planets are less important, and cosΘ2 −1, where Θ2 is the angle between the axis can be ignored. This is not a problem, as we in section ∝of t|he eclip|tic and u, the velocity of the WIMP as seen VII add the effects of other planets (apart from possible fromEarth[5]. Theprobabilityforareactionwithcross– resonances). We will in this section closely follow Gould section σ, can now be calculated, [6], with some small modifications. p(σ) 1 σ 4l 1 = tanΘ1 . (13) T cosΘ2 2πRl2πR| |Tχ ⊕ | | A. The probability of planet collisions. where we have divided by T to get the probability per unit time. The WIMP year⊕can be written in terms of We are interested in calculating the rate at which u(θ,φ,u)[5], the velocity of the particle in the frame of WIMPs with Earth crossing orbits comes into close en- the Earth, counter with the Earth. This will be used to estimate how the Earth affects the WIMP distribution. A close T = 1 2 u cosθ u2 −3/2T , (14) encounter is an event were the particle’s impact param- χ − v − v2 ⊕ eter is smaller than or equal to some value bmax(u). and Θ1(u) and Θ(cid:0)2(u) ca⊕n be expre⊕ss(cid:1)ed in u, θ and φ: Let’s imagine the Earth as being spread out on a flat hri,nagsoifnifingnuerrera3d.iNusowR,coonustiedrerraadipuasrtRic+lelwaitnhdptehriichkenlieosns cosΘ1 = Rˆ ·vˆχ = RR·u(u++vv⊕) | ⊕| less than the planet orbit radius R and aphelion greater usinθsinφ = , (15) thanR. SuchparticleswillbesaidtohaveEarthcrossing (u2+v2+2uv cosθ)1/2 ⊕ ⊕ 6 u (v R) u′ cosΘ2 =sinθcosφ= · ⊕× , (16) uv R ⊕ u cotΘ1 = δ v ⊕ sinθsinφ η (1.7) 1+2(u/v )cosθ+(u2/v2)(1 sin2θsin2φ) 1/2 u ⊕ ⊕ − By s(cid:0)ubstituting and rearranging we conclude th(cid:1)at the yearly reaction probability for an event with cross– section σ is given by p(σ,u) 3 σ v = ⊕γ(u)−1,where (18) T 2πR2 u ⊕ FIG. 4: Scattering off the Earth in velocity space. A fixed γ(u)= impact parameter fixesδ, but η is evenly distributed. 3πsin2θ sinφcosφ (1 2u/v cosθ u2/v2)−3/2 2 1+2(u|/v )cosθ+| (u−2/v2)(1⊕ sin2−θsin2⊕φ) 1(/129.) possiblybeweaklycapturedtoit)withineachperihelion ⊕ ⊕ − precession revolution. (cid:0) (cid:1) The gravitational scattering probability is dependent Eq. (19) was first derived by Gould [5] in a very similar ontheangulardistancebetweenthevelocitiesbeforeand way. Among other things, he used it to calculate the afterscattering: uandu′,suchthatsmalldeflectionsare ”typical timescales” at which particles diffuse between more common. The angle can be related to the impact different velocity space regions in the absence of solar parameter b, depletion. It is also used for calculating the probability of weak scattering of WIMPs at the Earth. bu2 The equations above are derivedunder some (geomet- δ(b)=π 2arctan , (20) − MG rical)approximationswiththe aimofgetting the correct as well as scattering probabilities on average. There are however a few pathologicalcaseswhere the geometricalmodelused δ(uˆ′,uˆ)=arccos(uˆ′ uˆ). (21) above breaks down. This happens when φ=0, φ=π/2, · θ = 0 and θ = π, in which case the probabilities above whereˆdenotes unitvectors. The scatteringangleabove are unphysical. Since this only happens for these few is the one given by Rutherford scattering (see e.g. [13]). special cases we will artificially solve this by adding a Gould used an approximate formula when deriving the small angle (of about 1 degree)to θ and φ when close to typical time scales [4]: δ(b) = R⊕ve2sc/(bu2). The two theseregions. Notethatinprincipal,the problemscould differ at very small impact parameters, and we use the be resolved, by making ∆ξ a function of the full set of full expression in our calculations. orbit parameters, but this is unnecessarily complicated As mentioned before, scattering can only change the for our purposes. For the interested reader, we refer to directionandnotthe velocity,andwearethereforedeal- a detailed investigation of the mathematical properties ing with random walk on spheres of constant u. The of γ as presented in [5] and [6]. To test our solution of directionη ofthescatteringisevenlydistributed,asseen addingasmallangleinthesepathologicalcases,we have in Fig. 4, where the scattering setup is shown. The arc investigatedthe effect of further increasingthe smallan- length is fixed by δ(b), but the scattering direction is gle added and conclude that the actual value chosen is evenly distributed. not important for the final results. This is reasonable, The cross–section for scattering between δ and δ+dδ since orbits in the vicinity of these critical regions are is dσ =2πbdb, so the yearly probability for scattering in quickly deflected into other orbits anyway. this range is, (using Eq. 18) dp(u,b) 3 2πb v = ⊕γ(u)−1. (22) dbT 2πR2 u B. Gravitational scattering on a planet ⊕ This can be rewritten in terms of the scattering angle Now that we have learnt how to calculate the proba- δ(uˆ′,uˆ), bility for particles to come into close encounter with a dp(u,b(δ)) given planet, it is time to apply this to gravitationaldif- = dδT fusion. FortheEarth,weweremainlyinterestedinthose particlescrossingthesphereofoneAUduringeachrevo- 3 ⊕2π v⊕γ(u)−1(M⊕G)2 cos (δ(uˆ′,uˆ)/2) (23) lution,sincetheyhaveachanceofhittingtheEarth(and −2πR2 u u4 2sin3(δ(uˆ′,uˆ)/2) 7 a sphere of constant velocity in a given time must be an integral over the source cell i and the destination cell j, dN dP(β,α) ji = n(α) dt dt ZZα∈ΩiZZβ∈Ωj Inourcase,thedestinationspaceisconvenientlyspanned by the scattering angles δ and η. The density of bound orbits scattered from cell i to cell j evolves with time as dn dη dp(u,δ) ji = dΩ dδ n(u), (26) dt 2π dδT ZZΩi ZZKj ⊕ where K is defined to be the region in δ–η–space cor- FIG. 5: An example of the timescales of particle scattering. j responding to scattering from the i to the j cell. The The color bar indicates the time for which there is a 10% probability ofscattering an angle δ=π/2±π/64, depending scattering probability to the [δ,δ +∆δ] band is evenly onthepresentlocationoftheparticle. Ingeneral,timescales distributed overallcells in that region. Numerically this are shorter at lower velocities. By repeated close encounters is implemented as as loop over δ as measured from the with theEarth,particles may diffusealong thedotted circles center of the source cell. The probability is then dis- (actually spheres) of constant velocity only, keeping u con- tributed over all discrete cells whose centers are inside stant, but allowing for changes in θ and φ. The figure shows the current band. theφ=75◦ slice of theparticle velocity space only. Itisimportanttounderstandthatwhatweareconsid- ering is the movement of the particle orbits, as opposed to the particles themselves. This means that we do not since need to calculate the actual particle trajectories. When M G δ db M G 1 we areinterestedin the actualparticle densities,we pick b= u⊕2 cot2 ; dδ =− u⊕2 2sin2 δ. (24) orbits from the orbit densities. Note that there are two 2 points on the orbit that pass a given radius, but due to The right-hand side of Eq. (23) may look like a negative the perihelion precession, any given particle could show probability density, but this is artificial since integration upatanyoftheseorbitlocationsdependingontheangle should be done for decreasing δ’s. We integrate Eq. (23) ofperihelion. Sincewehaveanticipatedmirrorsymmetry analytically and use that expressionwhenever numerical in the plane ofthe solarsystem, eachparticle is smeared valuesofthescatteringprobabilitytoδ+∆δ areneeded. outatfourindistinguishableorbitlocationsonthesphere To get a feeling for the significance of the diffusion, ofconstantu. Thiscanalwaysbe done,regardlessofthe we solve Eq. (23) to obtain the typical time scales for existence of symmetries in the free distribution. If the scatteringagivenangletooccur. Asanexample,welook free distribution is not mirror symmetric in the ecliptic at the time scales for which the probability of scattering plane,itcanbeforcedtohavethis(inthiscaseartificial) with δ =π/2 π/64is 10%. This is illustratedinFig.5. symmetrybyaveraging,aslongasweareonlyinterested ± in the absolute capture of WIMPs in the Sun or planets. The equations derived above apply only to particles C. The bound orbit density and orbit capture from whicharealreadygravitationallybound to the solarsys- the halo tem. We now turn to the calculation of the bound orbit density capture rate; ∆n /T from the distribution of jf Let us now define the bound orbit density, n(u) to be free particles. We will use the⊕local phase space density the number of bound particle orbits per infinitesimal ve- Ff(u) [Particles/(m3 m/s)]. locity and solid angle on the velocity sphere. The orbit ConsiderthedistributionofparticlesFf(u)passingthe density is thus free from information about the particle Earthwithimpactparameterb. Thenumberofparticles location along its elliptical orbit. The total number of scattered an angle δ(b db/2) in a given period of time ± bound particle orbits in a thin shell of radius u is T, is Tu2πbdb F (u). (27) dN =duu2 dΩn(u) (25) f ZZΩ=bound orbits dV We will now divide each velocity sphere into cells (that According to Eq. (20|) t{hzey}are scattered an angle δ(b). canatthispointbeviewedasinfinitesimallysmall). The Usingtherelations(24)weconcludethattheboundorbit number of particles scattered between two locations on density at the cell j will evolve with time as 8 dn dη (M G)2 cos (δ/2) dtjf =ZZΩfreedΩZZKjdδ 2π (cid:18)− 2πu u⊕4 2sin3(δ/2) Ff(u)(cid:19), (28) caused by gravitational scattering from the halo. We bound orbits and from free orbits. We have one main now have equations for gravitational diffusion as well as piece remaining to be studied, and that is the effects of capture to the solar system. solar depletion, i.e. how much of the bound WIMPs are actuallycapturedbytheSun,thusreducingtheirdensity in the solar system. D. Relating the phase space density F(u) and the We have done this by numerically calculating the ac- bound orbit density n(u) tualmotionfordifferentWIMPorbitsinthesolarsystem over 49 million years. As a measure of the quality of the Theideasofthelastsectioncanbeusedtowritedown numerical methods, we have also calculated the fates of an expression for the phase space density, which is what the 47 asteroids studied by Farinella et al. [7], as pre- we need for the weak capture calculations. The relation sented in appendix A. betweenthe phasespace density F(u)andthe bound or- bit density; n(u) is derived as follows. For a given orbit in the population of bound orbits, A. The numerical methods and conditions we use Eq. (18) to calculate the number of orbits that will pass trough anarea σ eachyear. We now consider a We have numerically integrated the orbits of about volume dV in space with base area σ and height h such thathisparalleltou. Aparticlepassingtroughthearea 2000particlesintypicalEarthcrossingorbitsinorderto estimate the solar depletion. The particles were spread will spend a time h/u in the volume. This means that out on the bound velocity space with random initial po- the fraction of the WIMP year spent in the volume in sitions on the Earth’s orbit. We have mainly used the case of an encounter is Mercury package [14] by Chambers for the integration. h It has the most important numerical algorithms, such . uT as Everhart’s 15:th order Radao [15] with Gauss–Radao χ spacings, and the equally well-known Bulirsch–Stoer [16] ThefractionoforbitspassingthroughσeachWIMPyear algorithm. Both are variable step size algorithms dedi- is catedtomanybodyproblems,andarecommonlyusedin p(σ,u) asteroidresearchforproblemssimilartoours. Thepack- Tχ age also includes a set of symplectic algorithms, which T ⊕ have been used for some tests. By looking at some test Therefore, since F(u) is the number of particles per orbits,we found that the symplectic algorithms (at least du3dV, the relation between F(u) and n(u) is asimplemented inthe Mercury package)wereslowerand less accurate for our setup. The tested symplectic al- h p(σ,u) gorithmswere ”MVS: mixed-variablesymplectic” [17] as F(u)dV =n(u) T . (29) uT T χ well as ”Hybrid symplectic/Bulirsch-Stoer” [14]. χ ⊕ The calculations included the test particles, the Sun, or the Earth, Jupiter and Venus. Other planets were not p(σ,u) 3 dV v included as they are believed to be sub-dominant. The F(u)dV =n(u)h uT = n(u)2πR2u⊕2γ(u)−1 (30) Bulirsch–Stoer algorithmwasusedto calculatethe orbits ⊕ of all test particles, as well as the planets, during a time OneshouldnotethatbyconstructionF(u)aboveisvalid of 49 million years. This took about 35000 CPU hours, in the frame of the planet. However, the right hand side onavarietyofLinux andAlpha machines. Awide range of the equations above presumes the planet to have a ofdifferentaccuracyparameterswereused,from10−14to constant velocity during the encounter so that du3 are 10−8, to evaluate the role this plays. The numericalrep- equal to the velocity volume element in the frame of the resentationoftherealnumberslimitsthebenefitofgoing Sun, dw3. past about 10−12. The final choice of 10−10 is a balance between time and accuracy. A recent publication [18] in the subject of numerical simulations of a special set of V. SOLAR DEPLETION OF BOUND ORBITS Jupiter crossing asteroids, came to a similar conclusion; WhenusingtheBulirsch–Stoeralgorithmfortheircalcu- Intheprevioussection,weinvestigatedtheevolutionof lations, they found accuracy parameters in the range of the bound orbit densities due to scatterings from other 10−9–10−8 to give statistically similar results as 10−12. 9 In the comparisons carried out, this gave results very similar to those with higher accuracy parameters. The comparison with the Radao algorithm gave qualitatively similar, however not identical, results with a similar cal- culation speed. In some occasions however, the Radao algorithmgaveahighersolardepletionforparticleswith very high velocity (relative to the Earth), u & 50 km/s. Thisisnotofmuchconcernforourpurposesthoughaswe are mainly interested in much lower velocities for Earth capture to be efficient. For ordinary asteroid calculations, a point mass ap- proximation combined with collision detection is suffi- cient. Our case is a little more delicate since WIMPs may pass through the planets. To handle this, the grav- itational routines were modified to use the real gravita- tional potentials inside the planets. For Jupiter and Earth, we used ”true” mass distribu- tions [19, 20]. For Venus we rescaled the mass distri- FIG. 6: The time (linear scale) for ejection (blue/dark gray) bution of the Earth and removed the liquid iron core. andcaptureintheSun(yellow/lightgray)ofasetoftestpar- Other planets included in tests where assumed to be ho- ticles. Eachbinrepresentsonlyoneparticle,sothestatistical mogeneous. Theimprovementallowstheparticlestopass error is high. However, this figure is typical for all angles, throughthe planetswithoutbeinginfinitelyscatteredby exceptthattheplateauoffast solardepletionatlarge”back- a point mass, making the calculations more realistic and ward” velocities are raised when φ approaches 90◦. Some numerically stable. For completeness, it would be inter- particles survived in the Solar system for the whole of the esting to add more planets to the simulations, but it is simulation. Those particles are marked with black dots. unfeasible to do as it slows down the calculations too much. We also believe, that we have included the most lion revolution. Considering Fig. 3, it’s evident that this important planets in our simulations. is equivalent to a symmetry in the sign of φ claimed by Gould; that the φ and φ cases are identical. − The particle orbits were evenly distributed in velocity B. The results of the numerical simulations space, but we solve the diffusion equations on spheres of constant u, hence we interpolate our results. What we The solar depletion was mainly calculated for par- needtoextractfromournumericalsimulationsisthede- ticles in eight planes of u space, with the φ values pletion frequency, i.e. the expected depletion probability 0,15,30,45,60,75,90 and 30 degrees (the φ = 30◦ per given time. Since the form of the actual distribu- − − plane was used to investigate the expected radial mirror tion, of which the results of the numerical calculations symmetry ofthe results). Oursolardepletion resultsare are samples, are unknown, the most reasonable way to notas badasGould feared[8]; Mostofthe particlessur- estimate the depletion probability per unit time is vived two million years. Nevertheless, solar capture is too large to be ignored. Figs 6 and 7 show the φ = 75◦ 1 f = . (31) plane,andthetimesafterwhichtheparticleshittheSun. Sun T Sun We note that ejection is much more common at Jupiter- crossing orbits. This is in compliance with the fact that, Figure8 showsthe logarithmof1/f,interpolatedonto a according to the scattering model used here, the proba- sphere of constant u, namely u=40 km/s. bility of scattering for such orbits is high. The fact that there is a large region at 50 km/s where there are no − ejections or sun captures, is in agreement with the qual- VI. THE EVOLUTION EQUATIONS FOR ONE itative results by Gould, presented in his 1991 paper [6] PLANET (see his Fig. 3, where he assumes that the filling times are about the same as the time of ejection). Apart from Intheprevioussectionswehavepresentedtheanalytic thecalculationsshownhere,someextracalculationswere expressions for the scattering of bound orbits to other carriedoutforrelativevelocitieslowerthan15km/s. The bound orbits, Eq. (26), as well as capture from free to results of those calculations were incorporated and used bound orbits, Eq. (28). We have also, by numerical sim- in the same way as the others. ulations, estimated the rate at which orbits are sentinto Another important, however simple, result is that the Sun and thus captured. We are primarily interested thereseemtoexistamirrorsymmetryinthein-outdirec- in how the bound orbit density evolves with time, and tions. Thisisexpected,sinceparticlesmayhittheEarth will here write down the dynamic equations in a form both on its way out and on the way back on its perihe- suitable for numerical work. 10 A. The dynamic equations of the bound orbit density The bound orbit density develops in time in the fol- lowing way, dn dn dn duu2 j =duu2 ji ij dt "i∈bXound(cid:18) dt − dt (cid:19) dn dn dn jf fj sj + , (32) f∈unXbound(cid:18) dt − dt (cid:19)− dt # where n is the number of orbits in the small cell [31] j j of the sphere. The sum over i is the flow from and to the other bound cells. The n and n terms are jf sj representingcapturefromunboundorbitsandcaptureof FIG. 7: The time (log scale) for ejection and capture in the bound orbits by the Sun, while the n term represents Sunofaset oftestparticles. Thisfigureisidenticaltofigure fj the ejection of bound particles. 6, except for the time scale which is logarithmic. In this scale, it’s easier to see that there is a small region at −30 We will now reformulate Eq. (32) in matrix form suit- km/s where the solar depletion occurs directly. This is not able for numerical calculations. Let us first define our surprising,sincethisregioncorrespondstoparticleswithvery state vectors, low velocity in the frame of the Sun. The plateau of direct solar capture extends further in the special case of φ = 90 N s (not shown) which allows extremely elliptic, or radial orbits. X = n (33) i (The plane of start positions is then parallel to the ecliptic F  f plane.)   whereN isthenumberofparticlescapturedbytheSun, s n istheboundorbitdensityandF isthevelocitynum- i f ber density of free (unbound) orbits. If the cells i are small enough, the various densities can be considered constant over each cell. Using this and the fact that the η part of the integration is independent of n(u) and F (u),thismeansthatEqs.(26)and(28)canbewritten f as dn ji = pbbn (u) and (34) dt ji i dn jf = pbfF (u) (35) dt jf f The solar capture can be written in the same way: dn si =pscn (u) (36) dt si i The p:s can be considered as the probability per unit time to transfer particles/orbits from and to the various cells. A positive p means that we transfer to the cell and a negative p that we transfer from the cell. The pbb ii FIG. 8: The solar depletion at the u=40 km/s sphere. The element requires an explanation. This is the probabil- colorbarindicatesthelogarithmofthetypicaldepletiontime ity per unit time that an orbit in cell i is not scattered 1/fSun. Theregion totheright are thefree orbits, for which to another bound or free cell, i.e. this term includes all the solar depletion is irrelevant. At a ”backward” velocity of the scattering out to both other bound orbits, and un- 30km/s, theSun-depletionsis greater, in agreement with the boundorbits. Theprobabilityforsolarcapturethoughis previous figures of this section. In understanding this figure, handled separately by psc. As the various entries in the itmayhelptotakealookattheu=40km/slineoffigure1, si which corresponds to the central horizontal (φ = 90◦) plane statevectorX havedifferentunits(Ns isanumber,ni is the the orbit density and F is the number density), the of this figure. f required conversion factors are also included in the p:s. The cells can be of various size, and these sizes are also