Dark Radiation constraints on minicharged particles in models with a hidden photon

We compute the thermalization of a hidden sector consisting of minicharged fermions (MCPs) and massless hidden photons in the early Universe. The precise measurement of the anisotropies of the cosmic microwave background (CMB) by Planck and the relic abundance of light nuclei produced during big bang nucleosynthesis (BBN) constrain the amount of dark radiation of this hidden sector through the effective number of neutrino species, Neff. This study presents novel and accurate predictions of dark radiation in the strongly and weakly coupled regime for a wide range of model parameters. We give the value of Neff for MCP masses between 100 keV and 10 GeV and minicharges in the range 10^(-11)-1. Our results can be used to constrain MCPs with the current data and they are also a valuable indicator for future experimental searches, should the hint for dark radiation manifest itself in the next release of Planck's data.


Abstract.
We compute the thermalization of a hidden sector consisting of minicharged fermions (MCPs) and massless hidden photons in the early Universe. The precise measurement of the anisotropies of the cosmic microwave background (CMB) by Planck and the relic abundance of light nuclei produced during big bang nucleosynthesis (BBN) constrain the amount of dark radiation of this hidden sector through the effective number of neutrino species, N eff . This study presents novel and accurate predictions of dark radiation in the strongly and weakly coupled regime for a wide range of model parameters. We give the value of N eff for MCP masses between ∼ 100 keV and 10 GeV and minicharges in the range 10 −11 − 1. Our results can be used to constrain MCPs with the current data and they are also a valuable indicator for future experimental searches, should the hint for dark radiation manifest itself in the next release of Planck's data.

Introduction
Extensions of the Standard Model (SM) of particle physics with an additional hidden U(1) h gauge symmetry have recently gathered a wealth of attention. In the simplest realization [1,2], the only new particle included is a gauge boson which has received many names: paraphoton [3], hidden photon [4] or dark photon [5] to name but a few. Despite the fact that the particles of the SM are all singlets under the new U(1) h (hence it being hidden), the hidden photon (HP henceforth) can have kinetic mixing with the hypercharge boson. It is radiatively generated if there exist "mediator" fields (irrespective of how large their mass is) charged under both gauge groups [1,6]. The most natural value of the kinetic mixing parameter is thus χ ∼ g ×O(10 −3 ) with g the hidden gauge coupling. Much smaller values occur when g is very small, like in [7][8][9], or when cancellations among different mediators happen, for instance if any of the U(1)'s is embedded in a non-abelian group. Typical values predicted in the literature yield a range χ ∼ 10 −12 − 10 −3 [9][10][11][12][13][14][15][16][17][18][19][20][21].
If the HP is massless, it has no phenomenological consequences whatsoever because the probability of producing a quantum of HP is proportional to the HP's mass. At the Lagrangian level, the only difference from a pure SM is a small renormalization of the hypercharge gauge coupling, g = g(χ). Since in practice g has to be measured one cannot know whether it contains a hidden contribution or not. However, if the hypercharge U(1) unifies with SU(2) weak isospin in a grand unified model, the value of g can be calculated theoretically given the weak coupling and one can constrain a HP contribution [53].
The unbroken hidden U(1) h case becomes very interesting when we consider additional particles with hidden U(1) h charge. Because of the small χ mentioned above, these particles appear as if they had a small electric charge = g χ/e [6] and we call them minicharged particles (MCPs henceforth) 2 . The existence of this type of MCPs does not challenge the standard arguments of the existence of magnetic monopoles and the quantization of charge but makes them more subtle [56].
Since the pioneer works [1,6] many experiments and phenomenological arguments have been devised to put the existence of MCPs into test. Direct laboratory searches for MCPs have been performed in accelerators [57], a dedicated beam dump experiment at SLAC [58] and ortho-positronium decays [59]. For MCPs of low mass (m f < 30 keV) the most relevant constraints come from stellar evolution. The stellar energy loss due to the emission of MCP pairs by plasmon decay has a number of consequences that can be constrained [57,60]. It delays the helium flash in red-giants (brightening the tip of the red-giant branch), accelerates the helium-burning stage [57,61,62] and the cooling of white-dwarves [57,62] and would have reduced the neutrino pulse of SN1987A [61,63].
In this paper we focus on cosmological probes of minicharged particles in models with a massless hidden photon. MCPs created in the early Universe can behave as dark matter (DM) and/or dark radiation (DR) and the HPs contribute to DR. Current cosmological data severely constrains the amount of DM, its possible interactions with the baryon+photon fluid and with itself, and the amount of DR. In order to translate this into bounds on the MCPs' and HPs' parameters one needs to accurately calculate the production and decoupling of MCPs and HPs in the early Universe. In the pioneer work of Davidson, Campbell and Bailey [57] these calculations were done analytically in simple approximations. They presented two important cosmological bounds. First, they derived an overclosure bound from requiring the relic density of MCP DM to be smaller than the critical density today (Ω < 1). Second, they used the contemporary constraint on dark radiation [64] (traditionally expressed in terms of the effective number of neutrino species N eff ) from the helium-4 produced in big bang nucleosynthesis (BBN). In the early Uni-1 Dark matter HPs can have their origin in the misalignment mechanism, in which case their mass can have much broader values [51,52]. 2 MCPs appear in different constructions in extensions of the standard model, see also [54,55].
verse, MCPs and HPs created during reheating or by interactions with the SM thermal bath quickly come into thermal equilibrium with each other (we assume that g is not hyper weak) constituting a thermal "dark sector" (DS henceforth). If the kinetic mixing is large enough, the DS and the SM thermalize with each other as well. If these DS particles are relativistic during BBN, their contribution to N eff is 2 + 8/7 = 3.14 (for a Dirac fermion MCP), which was ruled out by data back then. MCPs with masses m f MeV can avoid this bound. When the Universe's temperature reaches the MCP mass, the thermal abundance of MCPs becomes exponentially suppressed and, from that moment on, HPs have no means of interacting with the SM and decouple. All SM particles that become non-relativistic afterwards give their entropy to the SM bath heating it with respect to the DS. If m f is sufficiently above MeV (the key temperature range for the BBN bound) there are enough particle species in the SM to dilute the HP density below the observational bound. This lead to the very strong bound: m f > 200 MeV. Alternatively the MCPs should have never been in thermal equilibrium with the SM bath, which was estimated to happen for < 10 −8 . However, in that work it was incorrectly assumed that MCPs give all their entropy to HPs when they decouple at T ∼ m f . Since the amount of HPs was overestimated, so was the lower limit on the mass. In a later paper [61], the limit was corrected to m f > MeV after observing that for > 10 −8 the MCPs annihilate while still being in thermal equilibrium with the SM. Their energy is thus split between HPs and SM particles reducing the value of N eff . The situation for couplings around 10 −8 was never considered in any detail. It realizes an intermediate case between the assumption of all MCP entropy going to HPs (assumed in [57]) and being distributed equally into the SM and HP thermal populations. In this range of parameters one expects the bounds to strengthen because the DS comes close to equilibrium with the SM but the coupling between the SM and DS can be weak enough to favor the MCP entropy flow into the HPs, enhancing N eff .
The purpose of this paper is to update the cosmological constraints from dark radiation on MCPs in models with a HP, treating production and decoupling in full glory. This is timely because of the interest raised on these particles and hidden sectors in general and the considerable amount of cosmological data made available in the last decade. The nuclear reaction network of BBN is now better understood and brand new data on primordial element abundances has been collected and analyzed (especially deuterium [65,66] and helium [67,68]). The upper limit on the helium-4 abundance Y p < 0.2631 (95% C.L.) [69] stands as a reliable figure, regardless of assumptions on stellar processing or the uncertainties on the primordial baryon density and the neutron lifetime.
Furthermore, nowadays we have complementary information on the amount of dark radiation provided by the temperature anisotropies of the cosmic microwave background (CMB). Only recently the WMAP mission achieved enough precision to assess the existence of a cosmic neutrino background [70], i.e. N eff > 0. Combining CMB data with other late cosmology data sets -large scale structure, baryon acoustic oscillations (BAO) and direct measurements of the Hubble constant by the Hubble space telescope (HST)improves the measurements of N eff . The obtained values tended to be larger than 3, raising the excitement of a possible hint of new physics, see e.g. [71][72][73]. The latest CMB results of WMAP, SPT [74] and ACT [75] combined with BAO and HST gave N eff = 3.84 ± 0.40 at 68% C.L. [76]. Thus, although somehow controversial [77], rather than a constraint there seemed to be a 2-σ preference for a non-negligible amount of unaccounted dark radiation [78]. The recent results of the Planck mission [79] have unfortunately not clar-ified the issue. Combined with WMAP polarization (WP) maps, SPT and ACT, BAO and HST the Planck teams gives N eff = 3.52 +0. 48 −0.45 at 95% C.L. [79]. Although the error in N eff has decreased according to the expectations [80], the central value has done so too in such a way that the 2-σ excess remains. The Planck analysis have brought more information, revealing an increased tension between the HST direct measurement of the Hubble constant H 0 = 73.8 ± 2.4 km/(Mpc s) [81] and the lower estimate H 0 = 67.3 ± 1.2 km/(Mpc s) using Planck and other CMB data alone [79]. The value of H 0 is positively correlated with N eff in cosmological fits (see [82] and figure 21 of [79]) so that using the HST prior tends to push N eff to higher values. In other words, a high N eff softens the tension between CMB and local probes of H 0 but a systematic bias of the local measurements towards high-H 0 could be artificially triggering the excess 3 . It is also worth mentioning that although N eff > 3 reduces the tension of H 0 , it worsens the agreement between CMB and local measurements of the age of the Universe [84]. This reduction is in any case not significant [79,84,85]. When the HST prior is excluded from the analysis, the Planck team finds N eff = 3.30 +0. 54 −0.52 at 95% C.L. [79] and this is the value that we will use in this work. Planck has the potential to improve its sensitivity to N eff down to the ±0.2 level [80] and future polarization measurements can decrease this figure to the 0.05 level [86]. Thus there is still hope for a significant detection of DR in the future. We shall then present our results in a flexible and detailed way to allow the future user to derive stronger constraints or identify MCP parameters that fit an excess.
With this target in mind we have computed in detail all the processes leading to the production and decoupling of MCPs and HPs in the early Universe to track the amount of dark radiation present during BBN and later on during the CMB epoch. We can track the evolution of the energy density in the DS even for parameters where its coupling to the SM is only mild and thermalization with the SM bath is never complete. Pertaining this, we acknowledge a very comprehensive study of dark radiation in general extensions of the SM [87] which appeared recently. This excellent work covers partially the scope of this paper, touching on the MCP+HP case (sec 3.3.1 [87]). In comparison, we focus exclusively on it so we can explore a wider parameter space of couplings and masses, and discuss the role of the different production and decoupling channels. However, our works are complementary because [87] considers a range of different initial temperatures of the dark sector while we set it to zero to obtain conservative constraints.
Our results are summarized in figure 1 where we also show the most relevant constraints on MCPs in models where the minicharge arises as a consequence of kinetic mixing (g = 0.1). Other interesting constraints which are not shown have been discussed in [37,[88][89][90][91][92][93]. The Planck N eff constraint disfavors MCPs with m f < GeV down to ∼ 10 −8 − 10 −7 but leave a small region around m f ∼ 5 MeV (to be discussed further on). The BBN constraints cover this gap. They are similar to previous ones except for the region m f ∼ 100 MeV and ∼ 10 −7 where we find the mentioned strengthening of the BBN constraint due to the weak coupling between the hidden and SM sectors. Since our constraints have some overlap with astrophysical bounds at the lowest masses we computed the high-mass boundary of the helium-burning (HB), red-giant (RG) and white dwarf (WD) bounds more accurately. We also included the recent update [94] on MCPs' acoustic oscillations during recombination [95].
The plan of the paper is as follows. In sections 2 and 3 we present our definitions of The results of this work are: the constraint on N eff during BBN (dark blue) and on N eff by Planck (light blue). We have also improved the bounds from white dwarves (WD), red giants (RG) and horizontal branch (HB) with respect to the originals by calculating the high mass behavior. The remaining bounds are taken from elsewhere: LHC [23], DM [57], COLL [57], SLAC [58], OPOS [59], TEX [96] and CMB [94,95] (see also appendix C). the MCP+HP Lagrangian extending the SM. In section 4 we describe the equations and reactions ruling the evolution of the energy density of the hidden sector. In section 5 We present the bounds coming from N eff at the CMB epoch and explain different examples of the thermal histories encountered in different regions of parameter space. In section 6 we focus on the constraints from BBN and in section 7 we present our conclusions. The revision of the astrophysical bounds at high masses and the update on MCPs' acoustic oscillations is done in the appendix.

The model
In this article we extend the SM gauge group by an additional unbroken local U(1) h under which all SM particles are singlets. We also add a massive hidden fermion charged under the new U(1) h only. The additional terms to the SM Lagrangian then read where F µν , F µν are the field strength tensor of the HP and SM photon, respectively, f denotes the hidden fermion, m f its mass and χ is the kinetic mixing parameter after electroweak symmetry breaking. The covariant derivative is where g is the gauge coupling of the U(1) h and A µ is the vector potential of the HP.
Our results are only negligibly affected by physics at high energy scales so that we do not include mixing with the Z 0 boson or corrections from possible UV completions (e.g. SUSY). Since the U(1) h is unbroken, the HP remains massless. The kinetic part of the Lagrangian can be brought into the canonical form by rotating away the kinetic mixing through the redefinition A µ → A µ − χA µ . Without the coupling to f , the resulting A would be completely decoupled from the SM and thus unobservable. Including the hidden fermion, however, the redefinition induces a term −g χf / Af , i.e. the fermion gets an electric charge = g χ/e (e is the electron charge). Since χ is typically small, so is the electric charge. Therefore, f is called a minicharged particle (MCP).
The gauge coupling g can be of order unity and in the following we assume that the coupling between HPs and MCPs is strong enough to keep them in local thermal equilibrium (LTE) with a temperature T DS at all times. We check the limits of this assumption in section 4.2.

Number of effective neutrinos
The HPs and the MCPs contribute to the radiation density of the Universe. Assuming HPs and MCPs in thermal equilibrium the energy density of the dark sector can be computed in terms of a common dark sector temperature, T DS , and an effective number of DS relativistic degrees of freedom, g * DS , as where z = m f /T DS and the integral is defined to be 1 for z → 0. The HP is massless and always contributes, while the MCP contribution is exponentially suppressed once their temperature falls below their mass. The spectrum of the CMB is sensitive to the amount of radiation in the Universe at the epoch of matter-radiation equality and decoupling (T γ ∼ O(eV)), and Planck [79] was able to measure the energy density of radiation with unprecedented precision. The energy density of radiation is usually parametrized by the effective number of neutrinos N eff . At the CMB epoch this is defined through where g γ = 2 are the photon's degrees of freedom. In standard cosmology with only the SM particle content, N eff = 3.046 [97]. In our model N eff includes the contribution of the HPs and MCPs as well. The total N eff then reads where T DS , T ν and T γ are the temperatures of the DS, neutrinos and photons, respectively. The first term is the neutrino contribution, which in the MCP scenario can significantly deviate from the standard value.

Equations for the SM-DS energy transfer
To obtain N eff at the CMB epoch, we have to track the temperature ratios T DS /T γ and T ν /T γ during the evolution of the Universe. For the temperature of the DS, we study the evolution of the SM and DS energy densities with time. We use the following set of coupled differential equationṡ where ρ SM (ρ DS ) is the energy density of the SM (DS) particles, P SM (P DS ) is the pressure of the SM (DS) particles,˙denotes a derivative with respect to time, and the source term W encodes the energy transported from the SM to the DS sector per unit time by particle reactions, described further in section 4.1. H is the Hubble parameter, where M p = 1.22 × 10 19 GeV is the Planck mass. The neutrino energy density is contained in ρ SM up to T γ ∼ 3 MeV. We instantly decouple neutrinos at that temperature, afterwards tracking their energy density separately.
If the SM and DS are very strongly coupled (high ) we should have of course T γ = T DS and the cooling of the Universe is ruled by quasi-adiabatic expansion, in which the comoving entropy is conserved. In this case, one can compute N eff explicitly as a function of the SM-DS decoupling temperature with the formulas developed in [98]. The decoupling temperature in this case is never far from m f , when the MCP population gets exponentially suppressed. If we want to know the temperature more precisely, we need to accurately compute the time when the energy transfer between the sectors becomes inefficient by using the above equations (4.1). We have cross-checked our results with the entropy conservation hypothesis to find good agreement for large using a decoupling temperature T ∼ m f /10. For this comparison we have extended the formulas of [98] to cover smoothly the cases considered there, see appendix A.

Source term
The source term W of equations (4.1) is the particle physics' input on how efficiently the DS and the SM exchange energy. The most relevant reactions are 2 to 2 processes so their contribution can be written as a sum of terms of the sort where | M| 2 a+b→c+d is the matrix element for the reaction a + b → c + d summed over initial and final polarizations and E trans is the transported energy per collision. f x is the phase-space density, which for particles in LTE is either a Bose-Einstein or Fermi-Dirac distribution. To reduce the number of integrals analytically as in [99] we approximate the distributions by Maxwell-Boltzman type. Note that we neglect blocking and stimulation factors whose effects are always small. The one-particle phase space differential volume is where g x denotes the internal degrees of freedom of particle x besides spin, which is included in the matrix element squared.
To leading order, the following channels contribute to the source term W where particles can be replaced by the corresponding antiparticles for scattering, and e + e − can be replaced by other electrically charged particle/antiparticle pairs including mesons like π + π − . We divide the processes into three classes: Those that are efficient when the population of DS particles is very small ("production channels"), those that are most efficient when the DS population is sizable ("decoupling channels"), and other channels, which give small corrections.

Production channels
SM particle pair-annihilation (4.5a) and plasmon decay (4.5b) into MCPs are efficient even in the absence of a DS thermal bath. These channels produce an abundance of DS particles and bring the DS and the SM sector closer to equilibrium. The energy transfer normalized to the equilibrium value (∼ T 4 γ ) goes as W/ρ SM T 5 γ /T 4 γ when both species (the SM annihilated and DS created) are relativistic and decreases exponentially when one of their masses becomes smaller than T γ . The time interval is dt = dT γ /HT γ ∝ dT γ /T 3 γ (radiation domination) and thus the integrated energy transferred, Wdt/T 4 γ ∝ 1/T 1 γ , is dominated by the smallest temperatures where both species are still relativistic, T 1 γ ∼ max{m f , m SM }. The contribution from a heavy SM particle is thus inversely proportional to its mass unless it is lighter than the MCP mass. Since we cannot probe MCP masses much above the GeV scale, we neglect contributions of W ± and tt to the annihilation (4.5a).
Before the QCD phase transition (Λ QCD ∼ 180 MeV [100,101]) we include all contributions from elementary particles with masses smaller than the W -bosons. Afterwards, we should replace the contributions from quarks by mesons. As a compromise between simplicity and accurateness at the lowest energies we have only considered the contribution from the charged pions. Mesons and baryons more massive than the pions have their abundances already exponentially suppressed at the QCD phase transition already 4 .
For the pair production process (4.5a), the matrix element in the center of mass frame is (4.6) where s is the center of mass energy squared, m l the SM particle mass and Q l its electric charge. In practice, the lighter of the two masses has a subdominant effect on W so we neglect its contribution. In π + , π − annihilation, we have a similar expression from [103,104] where m ρ = 775.26 MeV is the mass of the ρ(700) meson [102], and Γ ρππ ≈ 149.1 MeV its decay width into pions [102]. This is the simplest form of F π (s) based on the vector dominance model [103,104]. The contribution of plasmon decay (4.5b) to Γ, i.e. the energy density transferred per unit time to the DS, takes the form where the plasmon decay rate in the comoving frame is where m γ is the photon plasma mass defined by the dispersion relation ω 2 − k 2 = m 2 γ and Z = Z(ω, k) is the renormalization factor [60]. We have to sum over photon polarizations, the two transverse and the longitudinal, for which m γ and Z are different. We are interested in a plasma made of relativistic particles where transverse plasmons dominate the decay rate [105] so we neglect the longitudinal mode and use ω m γ . In this case the plasma mass at first order in α is where the sum goes over all charged particle species, g i controls the spin and color multiplicity of the particle and Q i is its electric charge. The renormalization factor is Z ∼ 1 unless ω ∼ m γ .

Decoupling channels
While SM particle annihilation (4.5a) and plasmon decay (4.5b) decrease approximately like exp(−2m f /T ) once T ∼ m f , Compton scattering (4.5e) and Coulomb scattering (4.5g) only decrease as exp(−m f /T ). These channels are therefore prone to dominate when T m f . However, they are only important if T DS ∼ T γ since they need a sizable abundance of MCPs to be effective.
For Compton scattering (4.5e) we use the following matrix element in the rest frame of the MCP [106] | M| 2 where ω is the angular frequency of the incoming photon and ω is the angular frequency of the outgoing HP. These two frequencies are related by Compton's formula . (4.13) The calculation of W for Coulomb scattering (4.5g) is quite cumbersome due to the forward Coulomb divergence. In vacuum we find [107] | where s, u, t are the Mandelstamm variables and m is the highest mass involved in the process. For pions we find where m π is again the charged pion mass. We include a form factor [108] (4.16) Note that in the Coulomb energy transfer integral we neglect the mass of the lighter particle, which slightly underestimates W when both masses are similar. Matrices (4.14) and (4.15) diverge at low t. The energy transfer, E trans , has to vanish at t = 0 so E trans ∝ t at low t and, after the t-integration, W is logarithmically sensitive to the cut-off. Since energy transfer via Coulomb scattering is most important when the MCP is decoupling and hence in the verge of being non-relativistic, we can take the cut-off to be the Debye screening momentum k D which in a relativistic plasma is just m γ . We implement this minimum plasma screening by the substitution t 2 → (t − m 2 γ ) 2 in (4.14). We have checked the validity of that prescription in a few cases of interest by comparing with the energy transfer of a massive fermion in a QED plasma calculated in thermal-field theory including dynamical screening [109]. Since all SM particles can scatter on the MCP, we include all the channels discussed in 4.1.1.

Other channels
Another channel we include is the vector boson fusion process (4.5c). It is neither important for T DS T γ nor for T DS < m f but gives corrections for T γ ∼ m f . We include this channel using the following matrix element [106] | M| 2 γγ →ff = 8 2 e 4 s 2 (cos 4 θ − 1) + 16m 4 (4.17) Processes (4.5d) and (4.5f) are of order O( 4 ) and can be neglected unless the MCPs are extremely light [91], a case already excluded by stellar evolution except in more involved models [110] which we do not considered here.

Initial conditions
To compute the temperature ratios with the system of equations (4.1), we have to specify the initial conditions. In the spirit of a hidden sector we assume that the DS is absent after reheating (T DS = 0) and is dynamically created by SM reactions. As soon as some MCPs are produced via (4.5a) and (4.5b), they can generate HPs via ff → γ γ until the distributions of both MCPs and HPs are thermal with a common temperature T DS . Let us explore when this reaction is effective with a simple order-of-magnitude analysis, leaving O(1) factors aside. At high T γ , SM fermion pair annihilation creates a population of MCPs of density n MCP ∼ α 2 2 T 4 γ /H and typical momentum O(T γ ). The reaction ff → γ γ starts to be effective when this population is enough to ensure that one MCP suffers one annihilation per Hubble time, i.e. when which happens at Once this temperature is reached, the typical momentum of the DS particles is degraded due to the thermalization to T DS < T γ and the cross section σ(ff → γ γ ) increases, boosting the process. Thus g is sizable, this thermalization happens so fast that assuming a common temperature T DS is justified. One can show that the higher T γ the less our results will depend on the initial conditions. We, therefore, start our calculations at T γ ∼ 10 7 GeV and compute the initial condition T DS for given values of and m f using eq. (4.20). For MCPs with 10 −4 in the mass range of interest the DS and SM thermalize so fast that setting T DS = T γ is equivalent to starting with T DS = 0. For such large kinetic mixings we will therefore assume thermalization of the SM sector with the DS.

Numerical evaluation
Now that the initial conditions and the source term are fixed, we can track the temperature ratios down to the CMB epoch for different MCP masses and minicharge using (4.1). In the strong coupled regime ( > 10 −4 ) we linearize the source term W around T γ = T DS to improve the stability of the solver. We can then calculate N eff with eq. (3.3). The error in ∆N eff introduced by the linearization appears to be ∼ 5% for the phenomenologically interesting region of m f > 100 MeV (see section 6). As a cross check of our simulations, we have compared to [62]. We reproduce their results when we reduce our reaction set to plasmon decay and pair annihilation only.

Results for N eff at the CMB epoch
The MCP model we consider here has three parameters: m f , χ and g . It is more convenient to use = g χ/e instead of χ because this is the parameter that controls the energy transfer between the SM and DS in the early thermalization of the DS and during decoupling. Of the m f , , g set, the hidden coupling g is perhaps the least relevant. It controls the thermalization of the DS by itself -but the requirements are not very restrictive-and affects the decoupling through Compton scattering, which is dominant only for m f m e unless g is large. Note that, of all the reactions listed in Eqs. (4.5a)-(4.5g) only Compton (4.5e) and vector fusion (4.5c) depend explicitly on g . Thus we have decided to scan the parameter space m f , for just three representative values of g = 1, 0.1, 0.01. Based on these cases, we can extend our conclusions to further values of g . In this region, the minicharge is so small that the DS never reaches equilibrium with the SM and thus N eff does not deviate from the standard value. We see that T DS /T γ (dot-dashed line) increases slowly until T γ ∼ m f (in this particular case m f = 100 MeV) when the pair creation of MCPs becomes exponentially suppressed due to the lack of energetic-enough SM charged particles. At this moment T DS is only ∼ T γ /10 ∼ m f , the thermal population of MCPs becomes soon exponentially suppressed and the DS completely decoupled (no Compton or Coulomb processes are efficient because the MCP population is tiny). At T γ ∼ 3 MeV neutrinos decouple and the subsequent e ± annihilation heats the photon bath with respect to the decoupled neutrinos and HPs . This makes the temperature ratios T ν /T γ (dashed line) and T DS /T γ drop between T γ ∼1 and 0.1 MeV. The final temperature of neutrinos is the standard value T ν /T γ = (4/11) 1/3 . We marked this value as a dotted horizontal line here and in all the plots. For substantially larger couplings, the DS reaches equilibrium with the SM sector. If m f > 1 GeV, the DS freezes out when the light quarks and gluons are still present in the SM bath. These degrees of freedom will eventually heat up the SM bath with respect to the HP. The temperature of neutrinos has its standard history. In the particular case depicted, the DS and SM sectors reach thermal equilibrium when T γ ∼ m f , i.e. very close to their decoupling. The coupling of the sectors is still weak so that most of the MCP entropy heats the HP bath with respect to the SM. Thus, for a little while, until the QCD phase-transition triggers the disappearance of colored degrees of freedom, the HPs have a higher temperature than the SM.

• Region C (figure 4 left)
For large and m f in the approximate range (10−1000) MeV, the DS decouples after the QCD phase transition but before neutrino decoupling. In this temperature range the only SM particles that can heat the SM bath with respect to DS after decoupling are electrons, muons and pions, whose number of degrees of freedom are comparable to the DS. Thus, T DS ends up being close enough to T γ to contribute sizably to N eff . Below m f ∼ 100 MeV only electrons are relevant. The case depicted in figure 4 shows such a case. The MCPs have decoupled at T γ ∼ m f = 100 MeV in thermal    equilibrium with the SM so the SM and HPs have the same temperature. Later, the e ± annihilation epoch heats photons, but not neutrinos nor HPs which share the same temperature (4/11) 1/3 T γ . In this range we thus have typically N eff ∼ 3 + 8 7 = 4.14.
• Region D (figure 4 right) For smaller m f ∈ (1 − 10) MeV, the DS decouples before e ± annihilation but after neutrino decoupling. When the MCPs become non-relativistic, they deposit a part of their energy in the SM plasma. Thus, T γ /T ν becomes larger than in the standard scenario, and consequently the neutrinos contribute less than 3 to N eff . In figure 4 this is evident from the fact that T ν /T γ becomes smaller than (4/11) 1/3 . HPs get some of the MCP energy and might even get their share also from e ± 's so they have a sizable contribution to N eff .
• Region E ( figure 5 left) For m f < m e and > 2 × 10 −9 , very high values for N eff are realized. The reason is that e ± annihilation pumps energy not only into photons but into the DS as well. If the DS is weakly coupled, most of the MCP energy goes into HPs at decoupling and therefore we have T ν /T γ > (4/11) 1/3 and a large HP contribution. An example of this is depicted in figure 5. If the coupling is strong, the MCP annihilation dumps exactly the amount of energy needed to restore the standard value for T ν /T γ into the photon bath (a consequence of dealing with a fermionic Dirac MCP, which has the same degrees of freedom as electrons and positrons).
• Region F (figure 5 right) In this region, the SM-DS interactions are strong enough to bring the DS energy to values close to thermalization, but not sufficient to ensure the fast transfer of energy between the sectors when the MCP becomes non-relativistic and decouples. Thus most of the energy of the MCPs ends up in HPs, which get hotter than photons at least for some period of time. In figure 5 we depict such a case. MCP decoupling happens before neutrino decoupling so that the neutrino density adopts its standard value. T DS /T γ is higher than 1 after decoupling, but is eventually suppressed by e ± annihilation to end up falling below 1. Nevertheless, the DS contribution to N eff can be still quite sizable. In this example only the HP contributes as 2 effective neutrinos.
5.2 Cases g = 1 and g = 0.01 Most of the reaction rates depend on the combination (e ) 2 ≡ (g χ) 2 so a change in g can always be compensated by varying χ accordingly. This is the case for e ± annihilation (4.5a), plasmon decay (4.5b) and Coulomb scattering (4.5g), which are responsible for production and decoupling in most of the parameter space. On the other hand, Compton scattering (4.5e) and vector fusion (4.5c) are proportional to g 2 (e ) 2 . Thus, increasing (decreasing) g increases (decreases) Compton scattering and MCP annihilation relative to Coulomb scattering, plasmon decay and e ± annihilation for fixed (m f , ).
Our results for N eff can be seen in figure 6 for g = 1 (left) and for g = 0.01 (right). The case g = 0.01 is virtually indistinguishable from the g = 0.1 case. This finding corroborates the fact that the Compton and MCP annihilation processes do not play a significant role in the value of N eff , at least for g ≤ 0.1. The only difference is a slight increase of N eff at low at the lowest masses m f m e . Low mass MCPs with almost thermal abundance can still mediate some energy into the DS by the Compton process after electrons have annihilated (when the Coulomb process is inefficient).
The g = 1 case is more interesting. All the isocontours of N eff move down in by a factor 2 ∼ 4, depending on the region. This indicates that the Compton and MCP annihilation processes have become the dominant energy-transfer reactions in the DS-SM decoupling.

Implications from Planck
The 95% C.L. Planck upper limit, N eff < 3.84 (Planck+WP+highL+BAO) [79], is marked with a red-dashed line in figures 2 (g = 0.1) and 6 (g = 1 and g = 0.01). The constraints are independent of g for g 0.1 and thus figure 2 is valid for g < 0.1 as well. Note that for low g one needs to check that the MCPs and HPs thermalize (the thermalization T γ of eq. (4.19) has to be larger than m f for our constraints to be consistent).
Let us also recall that the MCP relic abundance behaves as self-interacting dark matter, which is excluded by a number of arguments to be a dominant component of the observed cold dark matter [111,112]. For very small g the relic abundance can be significant and these bounds have to be taken into account.
All in all, the Planck analysis disfavors MCPs with masses between 14 MeV < m f < 390 MeV for a wide range of minicharges > 10 −7 . For larger minicharges, the bound improves so that m f < 1190 MeV is excluded for, e.g., = 10 −1 . Interestingly, a broad range of is favored in the ∼ 5 MeV mass range. In the next section, we show that this region is however ruled out by BBN.
We highlight the Planck result excluding the HST bias because it shows the potential of future N eff measurements to exclude robustly MCP masses up to GeV. The same spirit showed for instance in [87]. Including HST data implies N eff < 4 (Planck+WP+highL +BAO+HST) at 95% C.L. [79], which changes our results very little. We should however remark the extreme sensitivity of the MCP mass bound to N eff : relaxing the constraint to N eff < 4.2 would shift the constraint to m f > 1 MeV or so as emphasized previously [61]. A novel result of this paper is that this is only true above ∼ 10 −7 . The large-N eff peninsula around ∼ 10 −8 − 10 −7 would still be excluded. Let us once more remark that the discrepancy on the value of H 0 favored by CMB and the one implied by local measurements (HST) prevents to make strong claims about exclusion limits on N eff .

Constraints from big bang nucleosynthesis
The energy content of the Universe drives its expansion influencing the effectiveness of the nuclear reactions of big bang nucleosynthesis. The extra radiation predicted in the MCP+HP model described here increases the expansion with respect to the standard case. In such a Universe electroweak reactions freeze-out earlier, which implies more neutrons during BBN, and BBN itself happens earlier in time (so less neutrons decay). Since eventually all neutrons end up forming part of 4 He nuclei, our scenario implies a larger-than-standard 4 He yield.
Note that MCPs and HPs do not affect directly the relevant electroweak or nuclear reactions relevant in BBN, they do it only indirectly by affecting N eff (and thus H) and the baryon to photon ratio η (we have seen that MCPs and HPs can take and give entropy to the photon bath). There exist very accurate calculations of the relic abundance of primordial elements (helium, deuterium, lithium, ..) as a function of N eff and η, which are the only unknowns in the standard BBN scenario. However, when we include MCPs and HPs both N eff and η can evolve during the temperature ranges relevant for BBN (T ∼ 100 keV-2 MeV) and a simple rescaling of standard results is not always possible. Thus, we have adapted the BBN code employed in [113] to compute the primordial abundances of nuclei. As input we have the thermal histories that we computed in the previous section Dark green coloring denotes regions far away from the upper limit Y p < 0.263 (95% C.L.) [69]. The limit is given by the red dashed contour line. Orange and red regions are excluded by more than 2σ level. to evaluate N eff at the epoch of the CMB. We find that the 4 He abundance gives an additional interesting constraint on the parameter space of the MCPs. Isocoutours of the yield Y p = 4n He /n B (normalized to the total baryon density) are shown in figure 7. Using the constraint Y p < 0.263 [69], we can exclude MCPs with m f < 16 MeV for > 1.4 × 10 −8 . This eliminates the region 2 MeV < m f < 14 MeV still allowed by N eff at the CMB epoch. Note that this constraint is slightly more conservative than the recently suggested Y p = 0.254 ± 0.003 [68] (actually it corresponds to a 99% C.L. exclusion). Hence the combination of BBN and Planck data disfavors MCPs with s in the range 10 −7 −10 −8 for masses m f < 390 MeV.
We checked that deuterium does not give us any further constraint. We do not consider lithium in this study since already in standard BBN the amount of 7 Li differs from the SM prediction by more than 4σ [114,115].
For g = 1 (g = 10 −2 ), the BBN results can be found in figure 8. Again, the contour lines for g = 1 are shifted towards smaller minicharges compared to g = 0.1 and the results for 10 −2 are indistinguishable from g = 0.1.

Conclusions
In this paper we report our detailed calculation of the contribution of minicharged particles and hidden photons to the dark radiation of the Universe. Using our results we can conclude that the recent Planck data together with BBN constraints disfavors the existence of MCPs lighter than ∼ GeV unless their minicharge is very small 10 −9 − 10 −7 (depending on the mass). Our results extend to a broad range of hidden sector gauge couplings g 0.1 and we have also covered the case g = 1. The next generation of cosmological probes will be able to assess the existence of dark radiation with an estimated 1-σ error of 0.05. Thus, we offer predictions of N eff for a broad range of MCP masses and minicharges that we will allow in the future to strengthen the constraints or, more importantly to pinpoint the possible parameters of these elusive particles in case of a firm discovery of dark radiation.

Acknowledgements
We would like to thank Yegor Goncharov, Thomas Hahn and Georg Raffelt for useful discussions. J.R. acknowledges support by the Alexander von Humboldt Foundation and partial support by the European Union through the Initial Training Network Invisibles.

A N eff from MCP decoupling in LTE
If one assumes instantaneous decoupling, one can use entropy conservation to compute N eff for the case where a particle species decouples during the annihilation of another particle species. In our case, we compute N eff to be where (4/11) 1/3 is the standard neutrino/photon temperature ratio, T d ν (T d DS ) is the decoupling temperature of the neutrinos (DS), g * Se + e − (T d DS ) are the entropy degrees of freedom of the electrons/positrons evaluated at the temperature of DS decoupling.

B Astrophysical bounds at high masses
The astrophysical bounds from red-giant, helium burning and white dwarf stars [57,[60][61][62] are based on constraints on stellar energy loss. MCPs are produced by pairs in plasmon decay γ * → ff in stellar interiors and leave the star unimpeded contributing to the energy loss more efficiently than photons (only emitted from the surface). Transverse plasmons in such non-relativistic plasmas have a dispersion relation ω 2 − k 2 = ω 2 p . The relevant values for the plasma frequencies in the interior of helium-burning, Red Giant and white-dwarf stars are ω p ∼ 2, 18, 23 keV [60]. The energy loss per unit volume of transverse plasmon decay into massive MCPs is where T γ is the plasma temperature and the plasmon decay rate into MCPs is given by eq. 4.10. We have denoted by K 2 = ω 2 −k 2 the 4-momentum squared of the plasmon. The above equation can be straightforwardly extended into off-shell plasmons once we know their self energy Π(ω, k) in the medium because off-shell excitations are also thermally distributed [116]. Thus we have Q = In our case we can take Π T ω 2 p +iωΓ T . Γ T is the rate of Thomson scattering into free nonrelativistic ambient electrons Γ T = n e σ T with n e the electron density and σ T = 8πα 2 /3m 2 e the Thomson cross section. Therefore we find The decay of the plasmon into the MCP pair requires K 2 > (2m f ) 2 . When ω p > 2m f the new factor behaves like a delta function dω 2 δ(K 2 − ω 2 p ) enforcing the dispersion relation ω 2 − k 2 = ω 2 p because typically Γ T ω p . The pole contribution dominates the ω−integral and we recover the results of [57,[60][61][62] which considered only on-shell plasmon decay. Contrarily, when ω p < 2m f the pole does not contribute much and as we consider larger MCP masses soon becomes irrelevant. In this regime, our calculation reflects the process γ + e − → e − + ff where the MCP pair is emitted through an off-shell photon after a common Thomson scattering of a thermal photon. We have neglected the contribution of electron-nucleus Bremsstrahlung e + Z → e + Z + ff because it is subdominant at high MCP masses. In order to built the bounds shown in figure 1 we have computed the integral (B.3) and match the constraint at low MCP masses with the already existing bounds [57,[60][61][62]. We have colored as excluded all the regions where the millicharge is larger than the upper bound obtained. However, it is not completely clear what happens when we consider values of much above this boundary. For sure, the physics of stars will be very strongly modified but computing a self-consistent bound is extremely complicated. In principle there could be islands of parameter space where MCPs are trapped inside the star with a corresponding HP thermal bath and the energy loss is somehow quenched. We consider this unlikely, as a sizable amount of radiation would in any case be radiated for instance by our Sun and these particles should have produced some kind of signature on Earthly experiments. We thus conclude that all the colored region is most likely excluded up to the largest values of . Of course, these constraints are valid for models in which the minicharge arises by means other than the kinetic mixing.

C Minicharged particles during recombination
Following [94,95] we recomputed the bound on the MCP abundance during recombination. If MCPs couple strongly enough during recombination, they participate in the acoustic oscillations of the photon-baryon plasma. Comparing cosmologies with MCPs to the Planck data, [94] finds that the bound on the relic density of MCPs is in the strongly coupled regime. This value agrees very well with the analytical prediction by [57]. Figure 1 shows this result.