Dark Stars and Boosted Dark Matter Annihilation Rates

Dark Stars (DS) may constitute the first phase of stellar evolution, powered by dark matter (DM) annihilation. We will investigate here the properties of DS assuming the DM particle has the required properties to explain the excess positron and elec- tron signals in the cosmic rays detected by the PAMELA and FERMI satellites. Any possible DM interpretation of these signals requires exotic DM candidates, with an- nihilation cross sections a few orders of magnitude higher than the canonical value required for correct thermal relic abundance for Weakly Interacting Dark Matter can- didates; additionally in most models the annihilation must be preferentially to lep- tons. Secondly, we study the dependence of DS properties on the concentration pa- rameter of the initial DM density profile of the halos where the first stars are formed. We restrict our study to the DM in the star due to simple (vs. extended) adiabatic contraction and minimal (vs. extended) capture; this simple study is sufficient to illustrate dependence on the cross section and concentration parameter. Our basic results are that the final stellar properties, once the star enters the main sequence, are always roughly the same, regardless of the value of boosted annihilation or concentration parameter in the range between c=2 and c=5: stellar mass ~ 1000M\odot, luminosity ~ 10^7 L\odot, lifetime ~ 10^6 yrs (for the minimal DM models considered here; additional DM would lead to more massive dark stars). However, the lifetime, final mass, and final luminosity of the DS show some dependence on boost factor and concentration parameter as discussed in the paper.


Introduction
first considered the effect of dark matter (DM) particles on the first stars during their formation. The first stars formed when the universe was about 200 million years old, at z = 10-50, in 10 6 M haloes consisting of 85% DM and 15% baryons predominantly in the form of H and He from big bang nucleosynthesis. The canonical example of particle DM is weakly interacting massive particles (WIMPs). In many theories, WIMPs are their own antiparticles and annihilate with themselves wherever the DM density is high. The first stars are particularly good sites for annihilation because they form at high redshifts (the density scales as (1 + z) 3 ) and in the high-density centers of DM haloes. Spolyar et al (2008) found that DM annihilation provides a powerful heat source in the first stars and suggested that the very first stellar objects might be dark stars (DSs), a new phase of stellar evolution in which the DMalthough only a negligible fraction of the star's mass-provides the key power source for the star through DM heating. Note that the term 'dark' refers to the power source, not the appearance of the star. DSs are stars made primarily of hydrogen and helium with a smattering of DM (<1% of the mass consists of DM); yet they shine due to DM heating.
In subsequent works,  and Spolyar et al (2009) studied the stellar structure of the DSs and found that these objects grow to be large, puffy (∼10 AU), bright (∼10 7 L ), massive (grow to at least ∼500M ) and yet cool (∼ 6000-10 000 K surface temperature) objects. They last as long as they are fed by DM. In another paper, Freese et al (2010) considered the possibility of an extended period of DM heating and the consequent growth to supermassive DSs >10 5 M . By contrast, in the standard case when DM heating is not included, Population III stars (the standard fusion-powered first stars) 3 form by accretion onto a smaller protostar from ∼10 −3 M (Omukai and Nishi 1998) up to ∼140M and surface 3 temperatures exceeding 5 × 10 5 K. This higher surface temperature in the standard picture inhibits the accretion of the gas, as various feedback mechanisms become effective at those high temperatures (McKee and Tan 2008). For reviews of the formation of the first stars in the standard scenario where DM heating is not included, see e.g. Barkana and Loeb (2001), Yoshida et al (2003), Bromm and Larson (2004), Ripamonti and Abel (2005), Yoshida et al (2006) and Gao et al (2007).
WIMP annihilation produces energy at a rate per unit volumê where n χ is the WIMP number density, m χ is the WIMP mass and ρ χ is the WIMP energy density. The final annihilation products typically are electrons, photons and neutrinos. The neutrinos escape the star, while the other annihilation products are trapped in the DS, thermalize with the star and heat it up. The luminosity from the DM heating is where f Q is the fraction of the annihilation energy deposited in the star (not lost to neutrinos) and dV is the volume element. The DM heating rate in a DS scales as the square of the WIMP density times the annihilation cross section, as can be seen from equation (1). The WIMP density inside the star adjusts itself in response to changes in the star's baryonic mass profile, since the gravitational potential well of the star is determined by the baryons. (The DM profile responds to changes in the gravitational potential due to the conservation of adiabatic invariants.) In this paper, we investigate the dependence of DS properties on these two quantities: (i) the annihilation cross section and (ii) the density of the halo within which the star forms, as characterized by the concentration parameter.
For a short list of papers by various other authors who have continued the work of Spolyar et al (2008) and explored the repercussions of DM heating in the first stars, see , , Taoso et al (2008), Yoon et al (2008), Ripamonti et al (2009), Gondolo et al (2010), Ripamonti et al (2010) and Sivertsson and Gondolo (2010). Their potential observability has been discussed by Freese et al (2010) and Zackrisson et al (2010aZackrisson et al ( , 2010b. Today, there are in orbit several very powerful near-IR telescopes, such as the Spitzer Space Telescope, the AKARI satellite and the Hubble Space Telescope (HST). The results from those instruments are essential for the understanding of the formation of the first stars in the universe and the end of the 'dark ages'. One unique signature of DSs has been discussed by Zackrisson et al (2010a), who have shown that DSs with masses ranging between 700 and 1000M , 'which contribute only 0.7% of the stellar mass in this galaxy, give rise to a conspicuous red bump in the spectrum at rest-frame wavelengths longward of 0.36 µm (this corresponds to wavelengths longer than 3.96 µm at z = 10). Because of this, galaxies that contain many cool dark stars are expected to display anomalously red colours. A feature like this is very difficult to produce through other means.' The upcoming James Webb Space Telescope (JWST) can even detect this signature, assuming that DSs are not exceedingly rare. Another technique for detection that is extensively used in high-redshift object searches is to look for dropouts, J, H or even K band dropouts, in deep-field surveys of the sky. The Atacama Cosmology Telescope (ACT), a ground-based instrument, could also be used to detect IR backgrounds at high redshifts. In order to actually determine whether a detection with a near-IR instrument is indeed due to a DS, one would need 4 a confirmation from a spectral analysis. The Grand Magelan Telescope (GMT), scheduled for completion in 2018, could be used for this purpose.
The possibility that DM annihilation might have effects on today's stars was actually considered in the 1980s by various authors, such as Krauss et al (1985), Bouquet and Salati (1989) and Salati and Silk (1989). More recently, the effect on today's stars was re-examined under the assumption that DM is made of WIMPs (Moskalenko and Wai 2007, Scott et al 2007, Bertone and Fairbairn 2008, Scott et al 2008, or within the hypothesis of inelastic DM (Hooper et al 2010).

Boosted leptophilic annihilation motivated by PAMELA data
Recent measurements of cosmic ray positrons and electrons in the GeV-TeV range could have significant implications on our understanding of DM. The PAMELA collaboration (Adriani et al 2009b(Adriani et al , 2010 reported an e + flux excess in the cosmic energy spectrum from 10 to 100 GeV, reinforcing what was previously observed up to an energy of 50 GeV by the HEAT experiment (Barwick et al 1997). The FERMI-LAT collaboration (Abdo et al 2009a(Abdo et al , 2009b has found an excess in the e + + e − flux in the 100-1000 GeV energy range above the background given by a conventional diffusive model, albeit in conflict with a much larger excess in flux in the 500 GeV range previously reported by the ATIC collaboration 4 . It is worth mentioning that a simple power-law fit of the FERMI-LAT electron energy spectrum is possible, consistent with astrophysical sources or DM annihilation (e.g. Di Bernardo et al 2009.
If confirmed, there are several possible explanations for the positron excess. The signals could be generated by astrophysical sources, such as pulsars or supernova shocks (e.g. Boulares 1989, Profumo 2008, Hooper, Blasi and Serpico 2009, Ahlers et al 2009, Blasi 2009, Fujita et al 2009, Malyshev et al 2009, Shaviv et al 2009. Uncertainties in cosmic ray propagation in the galaxy leave open the possibility that standard cosmic ray physics could explain the signal (Delahaye et al 2009). Another possible interpretation of the data could be in terms of a signal from DM annihilation (Baltz et al 2002, Kane et al 2002, de Boer et al 2002, Hooper et al 2004, Hooper and Silk 2005, Nomura and Thaler 2009, Grajek et al 2009, Arkani-Hamed et al 2009, Barger, Gao, Keung, Marfatia and Shaughnessy 2009, Bergstrom et al 2009, Bai et al 2009, Bi et al 2009, Cholis, Dobler, Finkbeiner, Goodenough and Weiner 2009, Cholis, Goodenough, Hooper, Simet and Weiner 2009, Cirelli et al 2009, Fox and Poppitz 2009, Meade et al 2009, Zurek 2009, Meade et al 2010 or decay (Chen and Takahashi 2009, Chen et al 2009a, Ishiwata et al 2009, Yin et al 2009, Bajc et al 2010, Chen et al 2010. Others point out constraints against such an interpretation (Abdo et al 2010, Abazajian et al 2010. The prospect that DM has been detected, although indirectly, has stirred a great deal of excitement and a flurry of interest. There are several model-building challenges, however. If DM is to explain the positron flux excess in the cosmic-ray spectrum, in most models the products of decay or annihilation must be primarily leptonic since an excess in the anti-proton fluxes is not found (Adriani et al 2009a).

5
Moreover, in order to explain the observed signals, the annihilation rate needs to be of the order of 10-10 3 larger than for thermal relic. This could be explained by non-relativistic enhancements of the cross-section, such as the Sommerfeld enhancement (Sommerfeld 1931, Hisano et al 2004, Profumo 2005, March-Russell et al 2008, Cirelli et al 2009, Lattanzi and Silk 2009 or a Breit-Wigner enhancement (Feldman et al 2009, Guo and Wu 2009, Kadota et al 2010. Scenarios involving non-standard cosmologies have also been considered as a possible solution (e.g. Gelmini and Gondolo 2008, Catena et al 2009, Pallis 2010. Alternatively, the annihilation rate could be greater than the standard prediction if the dark matter is not smoothly distributed in the local halo of the Milky Way (Diemand et al 2008, Hooper, Stebbins and Zurek 2009, Kamionkowski et al 2010. For a review of the DM explanation for the cosmic ray e ± excess, see He (2009).
In light of these new DM models, this paper will study the possible modifications to the evolution of a DS for the case of leptophilic boosted DM. In Spolyar et al (2009), a comprehensive study of DSs for various WIMP masses was performed, but the annihilation cross-section was kept at its fiducial unboosted value. We return here to this problem, including the possibility of a boost factor for the cross-section. Our starting point will be Bergstrom et al (2009), who find regions in the (m χ , B) parameter space (mass of WIMP versus boost factor) that fit PAMELA/FERMI data based on three different classes of DM models. For simplicity, we consider only two of those models in the present study: from the leptophilic class, the µ channel case, where 100% direct annihilation to µ + µ − is assumed, and the AH4 model, a subclass of the Arkani-Hamed model (Arkani-Hamed et al 2009) that postulates that the Somerfeld enhancement is due to the exchange of a new type of light (sub-GeV) particle. In the AH4 case, the new force carrier is a scalar, φ, which subsequently decays 100% to µ + µ − . Somerfeld enhancement is now generated naturally via ladder diagrams for the χχ → φφ process, producing 4 muons per annihilation versus 2 muons in the direct annihilation case.

Effect of concentration parameter
A second focus of this paper will be to study the dependence of DS properties on the concentration parameter of the initial density profile of the halo within which the first stars form. To remind the reader, the concentration parameter (c) characterizes how centrally condensed the initial density profile is: a larger concentration parameter means that more of the mass is concentrated at the centre rather than at the outside of the DM halo (the precise definition is given below in section 2.3.2). Previous studies of DSs considered c = 10 and c = 2. Here we systematically examine a sensible range of concentration parameters. Very recent results from numerical simulations of structure formation seem to favor a 'floor' of c = 3.5 on the concentration parameter of early structures (Zhao et al 2003, Tinker et al 2010 (some halos can have a smaller concentration parameter if star formation begins before the halo is fully formed). At any rate, the value will differ from halo to halo even at the same redshift. Thus, to explore how the properties of the DS are affected by a change in the concentration parameter, we run the same simulation with three different values for it, c = (2, 3.5, 5), which characterizes the overall range for the concentration parameter DM halos hosting the first stars.

Canonical values
We assume a redshift of formation z = 20 for the first stars, a DM halo of 10 6 M and concentration parameters c = (2, 3.5, 5). We take the annihilation cross-section to be where the boost factor B varies between 100 and 5000 (depending on the particle physics model). The corresponding WIMP mass is taken from the models of Bergstrom et al (2009) that we described above and ranges from 100 GeV to 4 TeV.
A key question for the DS is the final mass, as the DS accretes more and more material. As long as there is a reservoir of DM to heat the DS, the star continues to grow. In the original work of Spolyar et al (2009), an assumption was made that the initial DM inside the DS annihilates away in ∼ 500 000 years for a spherical DM halo; here the DS grow to ∼1000M . In a later work by Freese et al (2010), this assumption was questioned due to the fact that DM haloes are instead triaxial, so that a variety of DM orbits can keep the central DM density higher for longer periods of time and the DS can grow supermassive >10 5 M . In reality, DSs will form in a variety of DM environments and will grow to a variety of masses. For the purpose of illustrating how DSs vary due to differences in the halo concentration parameter and also due to enhanced annihilation rates, we restrict ourselves to the first option for adiabatic contraction (AC), in which the DM originally in the star (due to AC) is the only DM available to the DS; that is, the DS can grow to ∼1000M .
In addition to this simple AC, we also consider the effect of captured DM on the first stars .
In this case, DM passing through the first stars can scatter off the baryons multiple times, lose energy and become bound to the star (direct detection experiments are based on the same physics: the scattering of DM particles off nuclei). Subsequently, the star builds up a reservoir made up of captured DM, which can power stars. In the 'minimal capture case', the DM heating from captured DM and fusion powers the star in equal measure once it reaches the main sequence.
In section 2, we describe the elements necessary for studying the stellar structure of the DS. In section 3, we present the results on the influence of varying the concentration parameter and boost factor on the formation and evolution of DSs and on their properties. We summarize in section 4.

Initial conditions and accretion rates
The standard picture of Population III star formation starts with a protostellar gas cloud that is collapsing and cooling via hydrogen cooling into a protostar (Omukai and Nishi 1998) at the centre of the halo. However, as was found by Spolyar et al (2008), DM heating could alter the evolution of the first stars significantly. As soon as the protostellar gas reaches a critical core density, DM heating dominates over all possible cooling mechanisms. The cloud condenses a bit more, but then stops collapsing and becomes a DS in equilibrium. At this point, DM annihilation can power the DS.
As the initial conditions for our simulations, we take a DS in which the baryons are fully ionized. This luminous object powered by DM annihilations has a mass of 3M , a radius of 7 1-10 AU and a central baryon number density of ∼ 10 17 cm −3 . We look for equilibrium solutions as described below. From this starting point, we follow the evolution of the DS as it accretes baryonic mass from its surroundings. We use the accretion rate of Tan and McKee (2004), which decreases from 1.5 × 10 −2 M years −1 at 3M to 1.5 × 10 −3 M years −1 at 1000M . At each stage in the accretion process, we again find solutions in hydrostatic and thermal equilibrium. Eventually the accretion is cut off by feedback effects; as the DM runs out due to annihilation, the DS heats up to the point where it emits ionizing photons that shut down further accretion. Feedback turns on once the surface temperature T eff reaches 50 000 K and is accounted for by introducing a linear reduction factor that shuts off accretion completely once the stars' surface temperature reaches 10 5 K.

Basic equations
We use the numerical code that was previously developed in a paper by two of the current authors . In this section, we review the ingredients of the code before demonstrating the modifications relevant to this paper.
One key requirement is the hydrostatic equilibrium of the star. This is imposed at each time step during the accretion process, dP dr Here, P denotes the pressure, ρ(r ) is the total density and M r is the mass enclosed in a spherical shell of radius r . We assume that the star can be described as a polytrope with where the 'constant' K is determined once we know the total mass and radius (Chandrasekhar 1939). The energy transport is initially convective with polytropic index n = 3/2, but as the star approaches the zero age main sequence (ZAMS) it becomes radiative with n = 3. The code interpolates between n = 3/2 and n = 3 to account for the shift in energy transport as the star grows in mass. We can determine the temperature at each point in the radial grid via the equation of state of a gas-radiation mixture, Here, k B is the Boltzmann constant, m u is the atomic mass unit and µ = (2X + 3/4Y ) −1 is the mean atomic weight. In all of the resulting models, T 10 4 K except near the surface, so we use the mean atomic weight for fully ionized H and He. We assume an H mass fraction of X = 0.76 and an He mass fraction of Y = 0.24, and that they will remain constant throughout the simulation.
At each point in the radial grid, T (r ) and ρ(r ) are used to determine the Rosseland mean opacity κ from a zero metallicity table from OPAL (Iglesias and Rogers 1996) supplemented at low temperatures by opacities from Lenzuni et al (1991) for T < 6000 K. The location of the photosphere is determined by the hydrostatic condition, 8 where g is the acceleration due to gravity at that particular location. This corresponds to a point with inward integrated optical depth τ ∼ 2/3; here the local temperature is set to T eff and the stellar radiated luminosity is therefore with R * being the photospheric radius. The thermal equilibrium condition is where L tot is the total luminosity output from all energy sources, as described below in section 2.4. Starting with a mass M and an estimate for the outer radius R * , the code integrates equations (4) and (6) outward from the centre.
The total luminosity output L tot is compared with the stellar radiated luminosity, as in equation (8), and the radius is adjusted until the condition of thermal equilibrium is met (a convergence of 1 in 10 4 is reached).

Dark matter (DM) densities
2.3.1. Initial profile and concentration parameter. The first stars form inside ∼10 6 M haloes. Simulations imply that DM halos have a naturally cuspy profile, but there is still some uncertainty about the exact inner slope of a DM halo: Diemand et al (2007), Springel et al (2008) and Klypin et al (2010). Luckily, a previous paper  showed that a DS results regardless of the details of the initial density profile, even for the extreme case of a cored Burkert profile (such a Burkert profile is completely unrealistic). In this paper, we use a Navarro, Frenk, and White (NFW) profile (Navarro et al 1996) for concreteness. We assume that initially both the baryons (15% of the mass) and the DM (85% of the mass) can be described with the same NFW profile, where ρ 0 is the 'central' density and r s is the scale radius. Clearly, at any point of the profile, baryons will make up only 15% of the mass. The density scale, ρ 0 , can be re-expressed in terms of the critical density of the universe at a given redshift, ρ c (z), via where c ≡ r vir /r s is the concentration parameter and r vir is the virial radius of the halo. We assume a flat CDM universe with current matter density m = 0.24 and dark energy density = 0.76. One of the main points of this paper is to study the dependence of DS properties on the value of the concentration parameter c. Hence we consider a variety of values, ranging from c = 2 to c = 5.

Adiabatic contraction (AC).
As the baryons start to collapse into a protostellar cloud at the centre of the DM halo, the DM responds to the changing gravitational potential well and becomes compressed. As described in our previous work, we use AC to calculate the effect of baryons on the DM profile. With an initial DM and gas profile and a final gas profile, 9 we can use the adiabatic invariants to solve for the final DM profile. We use the Blumenthal method ( Barnes and White 1984, Blumenthal et al 1986, Ryden and Gunn 1987 to calculate the adiabatic compression of the halo. In this case, the simplifying assumption of circular orbits is made. Angular momentum is the only nonzero invariant. Its conservation implies that In the case of circular orbits, M is the mass interior to the radius r of an orbit, and the indices f and i refer to the final and initial orbits, respectively. As mass grows inside the orbit, its radius must shrink and the DM profile steepens. The validity of this method in this context has been checked in work by Freese et al (2009), where a more precise algorithm, developed by Young (1980), has been used and a difference within at the most a factor of 2 was found. Whereas the Blumenthal method assumes circular orbits, Young's method only assumes spherical symmetry of the system. Therefore, only one of the three conserved actions is identically zero in this case. Namely the plane of each orbit does not change. The other two actions, angular momentum and the radial action, respectively, are nonzero and conserved. In view of Freese et al (2009), we are confident that the simple Blumenthal method is sufficiently accurate for our purpose.
In this paper, we assume that the entire DM moves on circular orbits. The DM will become exhausted once the entire DM on orbits interior to the DS has been depleted. The timescale for this to happen is of the order of a million years. This is probably an unduly cautious assumption. In a recent paper , we studied the case of triaxial haloes with large numbers of centrophilic orbits, which provided a much larger DM reservoir than we are considering, leading to supermassive DSs. In reality, DSs will form in a variety of DM environments and will grow to a variety of masses. In this paper, we could have grown the stars to supermassive sizes, but do not consider it necessary to do so for showing the dependence of DS properties on the two effects that we are studying.

DM capture.
Until now, we have only discussed the DM inside the DS due to AC. However, a DS's DM reservoir can be refueled by DM capture. This refuelling requires an additional piece of particle physics: the scattering of DM off the atomic nuclei inside the star. Some of the WIMPs from far out in the halo have orbits passing through the star. This DM can be captured and sink to the star's centre, where they can annihilate efficiently. The capture process is irrelevant during the early evolutionary stage of the DS, since the baryon density is not high enough to efficiently capture DM. However, at the later stages as the DS contracts towards the ZAMS, the baryon densities become great enough for substantial capture to be possible. This mechanism was first noticed simultaneously by  and .
The capture rate is sensitive to two uncertain quantities: the scattering cross-section of WIMP interactions with the nuclei and the background DM density. In terms of the relevant particle physics, we consider only spin-dependent (SD) scattering cross-sections with σ c = 10 −39 cm 2 .
We have previously presented in  details of our capture study and will not repeat them here. We wish to emphasize that we assume the same case of minimal capture that we studied previously in Spolyar et al (2009), in which case the heating from fusion and that from DM heating are taken to be comparable. A more extreme and interesting case of dominant capture, which could last as long as the DS continues to exist in a high DM density environment, was studied elsewhere  and can lead to supermassive >10 5 M DSs.

Energy sources
There are four possible contributions to the DS luminosity, from DM annihilation, gravitational contraction, nuclear fusion and captured DM, respectively.

DM annihilation.
The heating due to DM annihilation is given in equations (1) and (2) and dominates from the time of DS formation until the adiabatically contracted DM runs out. In order to compute the luminosity generated by DM heating, one needs to know what fraction of the total energy generated by WIMP annihilations is deposited in the star. This quantity, which we name f Q , will be different for various models of DM. In previous works (e.g. in Spolyar et al 2009), a fiducial value for f Q of 2/3 was used as appropriate for many typical WIMPS, under the following assumptions. Firstly, the final products of DM annihilation, after all unstable particles have decayed to the lightest states, are taken to be of three types: (i) stable charged particles (i.e. e ± ), (ii) photons and (iii) neutrinos. Secondly, the energy distribution was assumed to be roughly equal among the three final products listed above. The electrons and photons were taken to thermalize in a very short timescale with the star and deposit their energy, whereas the neutrinos escape; hence f Q 2/3. On the other hand, given a specific model of DM, one could compute the precise value of f Q using the energy distribution of all final annihilation products for the model under consideration. While this procedure could be performed for the specific models in this paper, still we use the standard f Q = 2/3 in order to make simple comparisons with our previous work 5 . Since the aim of this paper is to understand the effect of the boost factor on the DS properties, we fix the value of f Q to 2/3, the same as in Spolyar et al (2009). The differences between the boosted and unboosted cases will now be due only to the boost factor itself and the different masses of the WIMP in the two cases (rather than to the detailed values of f Q ).
Since f Q appears always multiplied by σ v , any ambiguity we have introduced by fixing the energy deposition factor to 2/3 can be traded for an ambiguity in the annihilation cross-section. To further justify our assumption of f Q ∼ 1, we use a result of Gondolo et al (2010), who have shown, using PYTHIA, that in the case of leptophilic DM models used for explaining PAMELA excess positrons as well as the excess in Fermi-LAT electrons, 'one gets f Q 0.56 almost constant in the range 200 GeV m 2 TeV . Both of the models we are investigating in this paper, namely the AH4 and the µ channel, belong to the leptophilic class, as the annihilation products in both cases are muons.

Gravitational energy.
Once the DS runs out of DM, it begins to contract; gravity turns on and powers the star. Although relatively short, this Kelvin-Helmholtz (KH) contraction phase has important consequences: it drives up the baryon density and increases the temperature, therefore leading to an environment where nuclear fusion can take place. For a polytrope of index n = 3, the gravitational contribution to luminosity was found in Freese et al (2009) using the virial theorem where E rad = dV aT 4 is the radiation energy.

Nuclear reactions.
Subsequently, nuclear reactions become important. We include the following nuclear reactions, which are typical of a zero metallicity star during the pre-main sequence evolution: (i) burning of primordial deuterium (assumed to have a mass fraction of 4 × 10 −5 ), which turns on rapidly once the stars' central temperature reaches ∼ 10 6 K, (ii) the equilibrium proton-proton cycle for hydrogen burning and (iii) the alpha-alpha reaction for helium burning. Since we track the evolution of the DS only until it reaches ZAMS, we do not need to consider the changes in stars' atomic abundances. To calculate the nuclear luminosity, we define We use the methods described by Clayton (1968) to obtain the energy generation rate, nuc . For the proton-proton reaction, we use the astrophysical cross-section factors of Bahcall (1989) and the He burning parameters of Kippenhahn and Weigert (1990).

Luminosity due to captured DM.
As we have seen in section 2.3.3, during the later stages of the pre-main sequence evolution, captured DM can become an important energy source. The luminosity due to DM capture is where and the factor of 2 in equation (16) reflects the fact that the energy per annihilation is twice the WIMP mass. In equation (17), ρ cap stands for the captured DM density profile. In all of the simulations, we consider the case of 'minimal capture', which corresponds to equal contributions to luminosity from capture and nuclear fusion when the star reaches the main sequence.

Results
On the whole, for all of the values of the boost factor and concentration parameter that we have considered in this paper, the results are roughly the same: the final DS is roughly ∼1000M , ∼10 7 L , and lives for ∼ 10 6 years. Thus if the e + excess seen in PAMELA is due to WIMP annihilation, the required leptophilic boosted cross-section is certainly compatible with the DS picture. However, there are interesting differences between the models, which we will discuss. Other than in the subsection immediately following this one, we consider four WIMP models. As motivated below, we focus on one boosted model denoted by AH4 with the following set of parameters: B = 1500, m χ = 2.35 TeV and c = 3.5. As our unboosted models, we take 100 GeV WIMPs with the canonical cross-section of 3 × 10 −26 cm 3 s −1 , and we consider three values of the concentration parameter, c = 2, 3.5 and 5. For comparison, the 'canonical case' considered by Spolyar et al (2009) was the unboosted 100 GeV case with c = 2. The relative boost factor between the AH4 model and the unboosted models is best described as follows. Since equation (1) tells us that DM heating scales as σ v /M χ , one can see that the AH4 is exactly equivalent to a 100 GeV WIMP mass with boost factor 150/2.35 64. In other words, the relative boost factor between the AH4 model and the unboosted cases is actually 64. It is important to stress that all results we will be presenting in this paper depend only on this relative boost factor, and not the actual boost B for the annihilation rate. Therefore it might be somewhat misleading to refer to the models we are considering in terms of a specific particle physics model, as any two models that have the same relative boost factor will produce identical results for the DS. However, we keep the terminology established in the literature, i.e. AH4 and µ channel models for conciseness.
Also, we assume throughout that the value of the boost factor is such that it gives the σ v required in order to explain the PAMELA/Fermi data. Namely σ v ds = σ v gal , where the two correspond to the value of the DM annihilation cross-sections in the DS and our galaxy, respectively. If we assume the Sommerfeld enhancement as the source of the boost factor, we might have σ v ds σ v gal , as the velocity inside the star, roughly ∼ 10 km s −1 , is much lower than v gal 300 km s −1 and the Sommerfeld enhancement is higher at lower velocities, up to a certain saturation limit. Therefore, the values we take for the boost factor should be seen as conservative. For clarity, we fix the boost factor to the values given by Bergstrom et al (2009).

'Final' luminosity for DSs with boosted DM
We investigate the 'final' luminosities and masses of DSs as they enter the main sequence. We study a variety of WIMP parameter ranges capable of explaining the PAMELA data; in particular, we follow the work of Bergstrom et al (2009) and take various boosted DM annihilation rates and WIMP masses from figure 1 or their paper, which gives 2σ contours in the enhancement factor-mass plane needed to fit PAMELA and Fermi data. In this section, we consider two of the particle physics models they consider: the µ-channel and AH4-type models.
In figure 1, we plot the 'final' luminosity at the end of the DS lifetime, when the DS has accreted to its maximum mass: the DM has run out and the star is about to enter the ZAMS, where it will be powered by fusion. We plot the luminosity as a function of the boost factor. The points represent individual simulation outputs. The solid (dashed) lines connect simulation outputs generated with the DM mass and boost factors from the 2σ contours in figure 1 of Bergstrom et al (2009) corresponding to the Fermi (PAMELA) confidence regions in Bergstrom et al (2009). We note that even if the boost factors range over more than one order of magnitude, from 100 to 5000 (the corresponding WIMP masses take values in the 1-4 TeV range), the final luminosities in all cases are relatively similar, ranging from ∼7 × 10 6 L to ∼9 × 10 6 L (the variation in luminosity for a fixed boost factor corresponds to the range of allowed m χ in Bergstrom et al (2009). For instance, for a boost factor of 1000, the upper bound of the Fermi contour in the left panel of figure 1 corresponds to m χ = 1.4 TeV; the lower bound corresponds to m χ = 1.9 TeV). As the boost factor increases, we note that the final luminosity gets slightly smaller. This is due to the fact that the ratio σ v /M DM typically also increases in most cases. This leads to faster depletion of the adiabatically contracted DM from the star, shortening the  Bergstrom et al (2009). The contour regions are generated using the range of WIMP masses for a fixed boost factor taken from the corresponding 2σ contours in figure 1 of Bergstrom et al (2009). The central star-shaped point in the right panel will be taken to be our designated boosted AH4 model and is consistent with both Fermi and PAMELA. time taken to reach the main sequence, and therefore reducing the amount of mass accreted and consequently the luminosity at that point.

Structure and evolution of DSs
In this section, we analyze the two effects of (i) boosting the annihilation cross-section and (ii) a variety of concentration parameters on the structure of a DS. As mentioned above, for the 'boosted' case, we choose one sample point from the AH4-type models, which corresponds to the large star-shaped point at the center of the right panel of figure 1. Henceforth, we denote by AH4 this point, which has B = 1500, m χ = 2.35 TeV and c = 3.5.

Luminosity evolution.
In figure 2, we compare the luminosity evolution for the four cases we have just described. In all cases, DM annihilation heating provides the dominant contribution to the DS luminosity until the DM runs out. At this point, Kelvin-Helmholtz contraction sets in and briefly provides the dominant heat source, until the star becomes hot and dense enough for fusion to begin. In the cases where we include capture, DM annihilation may again become important at the later stages leading to a new DS phase.
In both upper panels of figure 2, we take c = 3.5: the left panel is the boosted AH4 model, while the right panel is unboosted (the relative ratio of boost factors is effectively 64, as mentioned before). Due to the ambiguity in the value of the concentration parameter, we also study its implications on the evolution of the DS by running the unboosted 100 GeV case for c = 2, 3.5 and 5, respectively; the lower two panels are the unboosted case with c = 2 and 5. The lower left panel is the same as the canonical case studied in Spolyar et al (2009).

Effects of boost:
The DM heating is the most powerful in the boosted AH4 case, since equation (1) indicates that heating scales with cross-section. Thus at any given time during the DS phase, this model has the brightest luminosity, as can be seen in figure 2. Consequently, the DM is burned up more quickly in the boosted AH4 case, leading to the shortest DS lifetime.
In order to balance the higher DM heating, a larger radius is required, which leads to a lower central temperature and density (see figure 9). As the adiabatically contracted DM runs out, the evolution is similar in all cases, the final luminosity being of the order of 10 7 L .
Effects of concentration parameter: When we compare the three different concentration parameter cases, we note that the nuclear fusion sets in earlier for smaller values of c. The higher the concentration parameter, the more adiabatically contracted the DM available; therefore it takes longer to start the transition to the main sequence. This will lead to slightly higher final masses, which in turn translate to higher final luminosities. Nevertheless, as pointed out before, they are all very close to 10 7 L . In all three cases where σ v = 3 × 10 −26 cm 3 s −1 , we notice a flash in luminosity when both the gravitational energy and the energy due to annihilations of the adiabatically contracted DM are relevant. This happens at a time somewhere between 0.31 Myr (for c = 2 ) and 0.42 Myr (for c = 5). The same flash can be seen for the AH4 panel, but now the true maximum of the luminosity occurs after only 0.1 Myr. The maximum luminosity is ∼4 × 10 7 L . The main differences between the different WIMP models as regards the luminosity evolution in figure 2 are in the time during which the 'pure' dark star phase lasts. The higher the boost factor, the shorter this phase. Conversely, a larger concentration parameter c prolongs the DS phase since more DM is available.

Pre-main sequence evolution, HR diagram.
In figure 3, we plot the Hertzsprung-Russell (HR) diagrams for the four cases. One can see two distinct phases. First, the DS goes up the Hayashi track with a very steep increase in the luminosity yet relatively cool surface temperature, T eff 10 4 K. At the end of the Hayashi track, the star enters the Henyey track. This path corresponds to an almost constant luminosity, while the temperature increases fast, mostly due to the KH contraction phase. As a rule of thumb, once a star is on the Heyney track, its core should be fully radiative. The graphs end at a temperature of ∼ 10 5 K when the star reaches the main sequence.
Effects of boost: From the left panel of figure 3 one can see that the boosted AH4 case has the highest luminosity, due to the extremely efficient DM heating Q DM ∼ σ v /m χ forming a luminosity peak. However, as the AH4 case burns up its DM, its luminosity falls. The boosted and unboosted cases eventually cross over at a temperature of ∼ 10 4 K, and henceforth the unboosted case has a higher luminosity. Consequently, the boosted AH4 case has the lowest luminosity as the star moves onto the main sequence, as discussed above.
Effects of concentration parameter: In the right panel of figure 3, the trend is uniform: an increase in the concentration parameter leads to an increase in the luminosity. The difference is relatively small in the early stages of the evolution, at low temperatures. This is due to the fact that the adiabatically contracted DM density profile is not very sensitive to the concentration parameter; therefore about the same amount of DM heating will be generated in each case. However, for a lower value of the concentration parameter, the adiabatically contracted DM runs out faster as there is less DM available. This leads to a shorter 'pure' DS phase, as can also be seen from figure 2, and consequently to slightly lower final mass and luminosities.

Radius and effective temperature of the DS. Effects of boost on radii:
In the 'pure' DS phase, the star will be puffier for a higher σ v /m χ ratio, as can be seen from the left panel of figure 4. This is expected due to the much higher DM heating in that case: a larger radius is required in order to balance the DM energy production and the radiated luminosity, which scales as R 2 . For instance, in the boosted AH4 case, the maximum radius is about 2 × 10 14 cm, whereas for the 100 GeV non-boosted WIMP, the DS will have a maximum radius of 10 14 cm. However, as mentioned before, the KH contraction will set in earlier in the boosted case. This phase corresponds to the sharp decrease in radius seen in figure 4. The final radii as the DS enters the ZAMS are similar in both cases, around 6 × 10 11 cm.
Effects of concentration parameter on radii: From the right panel of figure 4, we notice a uniform increase in radius with the concentration parameter. Again, more efficient DM heating leads to a puffier DS. The maximum radii range from 8 × 10 13 cm to 10 14 cm for c = 2-5. After the KH contraction phase, the DS settles to a radius very close to 6 × 10 11 cm as it enters the ZAMS.
Effective temperatures: In figure 5, we plot the effective surface temperature as a function of time. At first, during the DS phase, T eff is relatively constant below 10 4 K. Once the DM starts to run out, the star contracts and heats up, leading to a sharp increase in T eff due to the onset of the KH contraction phase. Once nuclear fusion becomes the dominant energy supply and the star ceases to contract, the surface temperature reaches a plateau. The final value of T eff ∼ 10 5 K is always the same for all cases, regardless of the value of the boost factor or concentration parameter. The AH4 case (left panel of figure 5) starts as a cooler star, again typical of more efficient DM heating in that case due to the boost factor. In addition, in the AH4 case, the DM runs out more quickly, leading to an earlier increase in temperature. Regarding the various concentration parameters: during the DS phase, the surface temperature is roughly the same in all cases: ∼ 7 × 10 3 K. The DS phase lasts the longest for the highest value of c as this case has the most DM to begin with; thus the temperature starts to rise later for ever larger concentration parameters c.
3.2.4. Energy transport near the core. The DS starts with a fully convective structure, modelled by a fluid with a polytropic index n = 1.5. Then a radiative core starts to develop, which grows outwards, until most of the star is described by a polytrope of index n = 3. In figure 6, we plot the radiative gradient in the innermost zone at the centre of the DS. The dashed horizontal line illustrates the critical value for convection. Models above the line have a convective core, while models below the line have radiative cores. Models with more efficient DM heating-i.e. the models with higher values of c or the AH4 model-transition to radiative energy transport later; stars with more efficient DM heating require a larger radius. At a fixed luminosity with more efficient DM heating, the star must have a larger radius and lower densities to keep DM heating and the stellar luminosity in balance. With a lower density, the star has lower central temperatures. At lower temperatures, the number of bound states increases, which increases the number of bound-free transitions. Also, the number of free-free transitions increases. These two effects increase the opacity, which produces a larger radiative gradient and delays the transition from convective to radiative transport. Once the original DM in the star runs out and nuclear fusion begins, a convective core develops in all cases. The central gradient is high due to the fact that nuclear fusion takes place primarily in the core of the star. Similarly, once the star repopulates its DM due to capture, this new DM population is thermalized with the star and its density is also sharply peaked at the centre of the star. Thus both nuclear fusion and captured DM lead to a large radiative gradient in the core and therefore favour convection.
Again, the transition happens later for higher values of c. The time when the convective core develops corresponds almost precisely to the time when the captured DM heating becomes significant, as can be seen by comparing with figure 2: 0.31, 0.37 and 0.42 Myr for c = 2, 3.5 and 5, respectively. figure 7 for the four cases we considered.

Baryonic central density. The baryon central density is plotted in
The higher cross-section of the AH4 case at first leads to a puffier star (larger radius to keep the radiated luminosity at a level to balance the higher DM heating); thus it is not surprising that the AH4 case initially has a relatively lower central density. However, the initial DM in the DS runs out earlier in the AH4 case due to more efficient burning; hence in the left panel the two lines cross due to the earlier onset of the KH contraction (marked in the plot by the sharp increase of the densities) in the AH4 case.
The central density ρ b (0) scales inversely with the concentration parameter, as can be seen from the right panel of the same plot. Again, more DM heating (higher c) will lead to larger radii and therefore smaller densities. However, as opposed to the situation depicted in the left panel, in the right panel, the curves do not cross since the models with a larger concentration parameter have more DM, which delays the onset Helmholtz contraction. We come back to this in section 3.2.7. In all cases considered here, once the star goes onto the main sequence, the central density is close to 100 g cm −3 .
3.2.6. Mass as the DS enters the ZAMS. In figure 8, we have plotted the DS mass as a function of time. In all cases, the final mass when the DS enters the main sequence is ∼ 700-1000M ; however, there are slight differences for different models. The models with more effective DM heating-i.e. the AH4 model compared with the 100 GeV unboosted case-burn up their original adiabatically contracted DM the quickest and enter the KH contraction phase the soonest. This in turn leads to a smaller final mass, as feedback effects will shut off accretion sooner. Comparing the cases with different concentration parameters, we notice that the DS final mass is an increasing function of c. As previously explained, an increase in the concentration parameter leads to a longer DS phase and hence to more mass accreted. In the case of the unboosted 100 GeV WIMP with c = 3.5, the final mass is around 900M , whereas in the AH 4 case it is close to 700M . For c ranging from 2 to 5, the DS will have a mass in the range 800M -1000M when it reaches the main sequence.
3.2.7. Density profiles for DM and baryons inside the DS: amount of adiabatically contracted DM. In figures 9 and 10, we have plotted the density profiles of the adiabatically contracted DM and baryons, respectively. The AH4 model for the same stellar mass has a lower DM density than the canonical unboosted 100 GeV case by roughly an order of magnitude and also has a more extended profile. For instance, in the case of a DS of 300M , the values are 5 × 10 −11 g cm −3 (AH4) and ∼ 1 × 10 −9 g cm −3 (canonical), respectively. This effect can be attributed to the fact that at higher annihilation cross sections or lower particle masses 6 , a larger radius and a lower density are needed in order to balance the DM heating and the stellar luminosity (which scales as R 2 ).
During the early stages of the DS evolution, the dependence of the adiabatically contracted DM density on the concentration parameter is very small, at least for the range we have considered here. As can be seen from figure 7, prior to the onset of the KH contraction phase, the central baryon densities for models with different values of the concentration parameter have similar baryon density and DM density as well. Models with a larger concentration parameter have slightly more dark matter, have slightly lower central DM densities and are also more extended. Before the contraction phase, the central DM density of the c = 5 case's density is 10% lower than the c = 2 case. The radius is 20% larger. Models with different concentration parameters only begin to dramatically diverge once the star begins to contract and enters its KH contraction phase. At this point, the star begins to shrink, which cause the DM densities to increase dramatically. DS in halos with a larger concentration parameter have more DM and thus delay the onset of the KH phase. For c = 2, the contraction phase begins at t ∼ 0.28 Myr (see figure 4). At this time, the stellar mass has reached ∼700M . For c = 3.5 and 5, the contraction phase begins later. In the case of c = 5, the star has a mass of ∼850M . Thus the contraction begins once the star is more massive.
At a fixed stellar mass, the DM densities will differ dramatically between models that are in the contraction phase compared to those that are not contacting. For instance, let us consider a 750M DS. In the case of c = 2, the star has entered the contraction phase, while the cases with a larger concentration parameter (3.5 and 5) are yet to begin their contraction phase. Hence the stars with a larger concentration have a lower DM density and are more extended, which can be seen in figure 9.
Finally, in figure 11 we have plotted the amount of adiabatically contracted DM inside the DS as a function of time. One can also see that DM densities are many orders of magnitude lower than their baryonic counterparts at all times. Although the amount of DM never exceeds 0.4M , this is sufficient to power the DS all the way up to ∼1000M (where most of the mass is baryons) and 10 7 L . Indeed, DM heating is a very powerful energy source. Adiabatically contracted DM density profiles. Each line corresponds to a fixed value of the mass of the DS during its evolution. Note that certain lines that are mentioned in the legend do not appear plotted in all four panels. This is because at that stage the DS has exhausted all of the DM.

Summary and conclusions
In this paper, we have considered the effect on DSs of leptophilic boosted DM annihilation cross-sections, as would be typically required in order to explain PAMELA data. Secondly, we have varied the concentration parameter in a host of DS models. We have restricted our study to include the following two sources of DM: (i) the DM originally contained in the star due to AC and (ii) the minimal capture scenario. We have not included the additional DM due to extended AC or to maximal capture models discussed in Freese et al (2010). Nonetheless, the dependence of DS properties on the boost factor and concentration parameter can easily be seen in this paper. As our prototypical boosted case, we have focused on the AH4 model with the following parameters: B = 1500, m χ = 2.35 TeV and c = 3.5. In the unboosted cases, we have taken M χ = 100 GeV and three values of the concentration parameter, c = 2, 3.5 and 5.
We have found that if the positron excess observed by PAMELA is indeed due to leptophilic boosted DM, then there would be an early DS phase of stellar evolution powered by DM heating, lasting long enough to bring the star to substantially higher mass and luminosity than predicted for regular Population III zero-metallicity stars. Our basic results are that the final stellar properties, after the DS runs out of its original adiabatically contracted DM fuel, undergoes KH contraction and enters the main sequence, are always roughly the same: ∼1000M , ∼10 7 L , lifetime ∼ 10 6 years. Similarly, these same DS properties are also the basic result independent of the value of the concentration parameter in the range between c = 2 and c = 5.
We reiterate that these values are only for the case of simple AC and minimal capture; if we were to include the additional DM due to extended AC or maximal capture, then these values would be different by many orders of magnitude. However, the basic dependence on the parameters of interest in this paper would generalize. In particular, the result that the final mass, luminosity and lifetime are relatively independent of boosted cross-section or the value of c would still hold. In addition, the slight differences from model to model, discussed in the next paragraph, would also hold.
We have found that the lifetime, final mass and final luminosity of the DS, although roughly similar in all cases, do show some dependence on the boost factor and concentration parameter. We found that the DS lifetime is shorter in the boosted case, since the DM is exhausted sooner. On the other hand, the lifetime is longer for higher concentration parameter since there is more DM to begin with. Thus nuclear burning becomes effective earlier in the case of a boosted cross-section or a low concentration parameter. The DS accretes matter continuously while it remains powered by DM heating. Hence, the largest final stellar mass results are for the longest living DS, i.e. the unboosted case with the highest values of c. In all cases, the final mass is ∼1000M . We have shown the HR diagram for all cases, studying both the Hayashi track and the approach to fusion power. We have shown that the 'final' luminosity at the end of the DS lifetime is ∼10 7 L in all cases, with the luminosity decreasing slightly as a function of increasing boost factor or decreasing c, again because of the more rapidly depleted pool of DM. We have also examined the luminosity evolution of the DS as a function of time. During the DS phase itself, the luminosity is higher/lower for boosted/unboosted cross-sections. The reduced luminosity during the DS evolution in the unboosted case is a consequence of the reduced energy production, since DM heating is proportional to annihilation cross-section and the square of the DM density. Then at lower cross-section (unboosted), a smaller radius is needed to balance the DM energy production and the radiated luminosity. Consequently, the unboosted case has higher central densities (both for baryons and DM). Similarly, lower values of c have lower luminosity during the DS phase, smaller radii and higher central density. In all cases, the DM density is a minute fraction of the total density, with baryons dominating the gravitational potential; we have shown the density profiles of both components. Again in all cases, the amount of DM inside the star never amounts to more than 0.4M , a tiny fraction of a star that grows to ∼1000M ; yet this DM is sufficient to power the star. This is 'the power of darkness'.