Challenging the Stability of Light Millicharged Dark Matter

We investigate the cosmological stability of light bosonic dark matter carrying a tiny electric charge. In the wave-like regime of high occupation numbers, annihilation into gauge bosons can be drastically enhanced by parametric resonance. The millicharged particle can either be minimally coupled to photons or its electromagnetic interaction can be mediated via kinetic mixing with a massless hidden photon. In the case of a direct coupling current observational constraints on the millicharge are stronger than those arising from parametric resonance. For the (theoretically preferred) case of kinetic mixing large regions of parameter space are affected by the parametric resonance leading at least to a fragmentation of the dark matter field if not its outright destruction.


I. INTRODUCTION
Very light (sub-eV) dark matter (DM) must consist of bosonic particles. Common examples are axions, axionlike particles or dark photons [1][2][3][4][5]. By virtue of their tiny couplings and very low masses, it is often taken for granted that they make good and cosmologically stable candidates for DM. However, the fact that very light bosonic DM is long lived is far from trivial. Due to their low mass and low velocity, DM made from light bosons features very high occupation numbers which can dramatically enhance interaction rates with other particles and lead to parametric resonance phenomena [6][7][8][9]. For example, in significant parts of the parameter space the stability of axion-like particles towards their decay into photons requires a non-trivial interplay of the expansion of the Universe as well as plasma effects [1,2,10] (cf. [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] for some situations where Bose enhancement from high occupation numbers may lead to interesting signatures for axion-like particles). Of course, these particles do not carry a conserved charge that would naturally render them stable towards decay. One may therefore wonder what happens if the light DM particles are charged and only annihilations with suitable antiparticles are possible.
In general, having a conserved electric charge, stability towards particle decay is ensured. However, DM should not carry any net electric charge 1 . Therefore, it should be composed of an equal number of particles 1 In a scenario where very light bosonic DM does carry an ap-and antiparticles, opening up the possibility of annihilations, for example into radiation 2 . Naively, this seems to be strongly suppressed by the tiny value of the relevant charge that is required by phenomenology (for a review see, e.g., [55]), independent of the question of cosmological stability. However, for low masses, enhancements due to high occupation numbers may set in. Indeed, in this work, we argue that the coherent nature of the very light DM particles can drastically enhance the interaction rates with gauge bosons. This, in turn, can cause an annihilation into photons, even for tiny electromagnetic charges.
The basic reason for such an annihilation is that the DM coupling acts as an oscillating mass term for the photons in the coherent DM background. This can drive the gauge bosons into a parametric resonance [6][7][8][9], such that certain momentum modes are excited quite rapidly. The enhancement of the momentum modes then corresponds to an explosive production of photons. This phenomenon can lead to an efficient depletion of the DM energy density, that seriously challenges the cosmological stability of the DM candidate 3 .
This work is structured as follows. In Section II we discuss the phenomenon of rapid photon production in a millicharged DM background via a parametric resonance. Furthermore, we carefully examine plasma effects in the early Universe that are able to stop the depletion of the proximately conserved global charge, however, today's net charge density may be non-vanishing [53]. 2 In principle, a spatial separation of positive and negative charges may be possible [54]. While we do not have a conclusive argument excluding this possibility, we strongly suspect that such a situation is not viable in the case where the particles in question carry a gauge charge and are supposed to be the dominant form of DM (e.g. the presence of long range gauge interactions between the regions may modify the equation of state). Moreover, we note that this is, by definition, connected to a very inhomogeneous situation. 3 The amplification of gauge fields in the expanding Universe by a parametric resonance from charged scalars has also been considered as a source of large-scale primordial magnetic fields [56].
DM energy density. In Section III we consider the theoretically preferred situation where the millicharge arises from a hidden (or dark) photon kinetically mixed with its electromagnetic counterpart. A brief summary and discussion can be found in Section IV.

II. RESONANT DEPLETION OF MILLICHARGED DARK MATTER
In a scenario where DM carries electromagnetic charge, it is subject to annihilating into visible particles and, in particular, photons. In this section, we discuss that this depletion of the DM energy density can be drastically enhanced by parametric resonance phenomena [6][7][8][9], even for tiny charges. This implies that it is necessary to reconsider the cosmological stability of a millicharged DM candidate.
In first approximation the observed DM energy density can be understood as coherent oscillations of a spatially homogeneous 5 complex scalar field φ in the expanding Universe. We can decompose the field, In these coordinates, charge neutrality is ensured by trivial dynamics for the angular degree of freedom 6 , χ = const, and without loss of generality we can take χ = 0. The radial mode then oscillates with decreasing amplitude (3) 4 Here, we absorbed the gauge coupling into the definition of the charge, as there is only a single field involved. 5 Indeed, in the case of the misalignment mechanism and if the field already exists during inflation, any fluctuation is stretched out by inflation, thereby providing for homogeneity. 6 During inflation any non-trivial dynamics of χ is diluted quickly, Here, a(t) is the scale factor and t 0 denotes the time today with scale factor a 0 . Moreover, is the (average) oscillation amplitude today where we use ρ DM = 1.3 keV/cm 3 [78] as a reference value. The energy density associated to the scalar field then dilutes as ρ φ ∼ a −3 , as appropriate for a cold DM particle.
A. Rapid photon production via parametric resonance Obviously, the requirement that the energy density of φ scales like that of pressureless matter is not entirely sufficient for making it the DM. In addition, a viable DM candidate also has to be cosmologically stable. Crucially, since in our scenario φ carries electromagnetic charge, annihilation channels to photons are open, eventually challenging its stability. In the simple theory (1), the main example of a depletion mechanism would be the pairwise annihilation of DM particles, φφ → AA. As we will now show, even for tiny electromagnetic charges, the interaction rates of this channel can be significantly enhanced by resonance effects that lead to an explosive production of photons. To see this, let us consider photon modes in the classical DM background during the evolution of the Universe. These satisfÿ where A denotes a polarization mode of momentum k and H is the Hubble parameter. A collectively describes the spatial components of the gauge potential A µ , while we fix the temporal components by the Lorentz gauge condition ∂ µ ( √ −gA µ ) = 0. Once the DM field has overcome the Hubble friction, H m, it oscillates according to Eq. (3) with amplitude Φ(t). In this scenario, the equation of motion (5) can be rewritten as a differential equation of Mathieu type [79], with x = mt and we have defined The interaction with the millicharged DM acts as an oscillating mass term for the photons. It is therefore possible that some mode functions are enhanced by resonance effects, known as parametric resonance [6][7][8][9]. As the mode functions determine the occupation number of the gauge bosons, this process corresponds to a resonant production of photons.
Crucially, the solution of (6) contains an exponential factor, A ∝ exp (µ k x), with Floquet exponent µ k . This exponent is in general a complex number which, importantly, can have a positive real part 7 . For the purpose of our work, we will exclusively focus on the case where µ k is purely real and positive. This corresponds to an exponential growth of the respective momentum modes, That is, the rate of photon production is governed by the Floquet exponent µ k , which, in general, is a function of A k and Q. Depending on the dynamics of the Mathieu equation, a parametric resonance can be considered in two different regimes. In a narrow resonance (Q 1) only very few momentum modes are enhanced, while the opposite is true in a broad resonance (Q 1). Both regimes feature resonant instabilities, which ultimately lead to an explosive production of photons. For simplicity, let us collect a few important properties of these resonance bands for our discussion. For details on the dynamics of the Mathieu equation in general and its instability bands in particular we refer the reader to [79].

Narrow vs broad resonance
Let us first understand the characteristic behaviour of the Mathieu equation, while neglecting the expansion of the Universe. This serves as the basis on top of which we can later include the consequences of expansion.
In the narrow resonance regime, where Q 1, the instability bands of the Mathieu equation feature a small width. In our case, this means that only a limited range of momentum modes is resonantly enhanced. The first instability band is the dominant one, as it contributes to the exponential growth with the largest Floquet exponent. The latter is given by [79] Furthermore, to good approximation, in momentum space the resonance bandwidth reads [79] ∆k ∼ mQ .
In combination with the requirement Q 1 this justifies the name narrow resonance. Up to corrections of order Q, the Floquet exponent is maximal for momenta of the order of the DM mass, k * m, such that µ k * = Q/2. This is essentially the same result that one would obtain from a perturbative approach to interaction rates in φφ → AA processes, if Bose enhancement is taken into account (see, e.g., [80]).
In the broad resonance regime, where Q 1, the situation is more complicated. Here, exact expressions for the Floquet exponents are not available. Therefore, we will use analytic approximations of µ k given in [81,82]. As their precise form is not very enlightening, we do not quote the full expressions here. Instead, we give some approximate numbers and behaviours to facilitate the discussion. That is, typical values of the exponent for a wide range of momenta are µ k ∼ 0.15 and it can obtain a maximum value of µ k * ∼ log 1 + √ 2 /π ≈ 0.28 [81,82]. As a rough approximation one can therefore estimate [81,82] within the resonance bandwidth. Importantly, this does not strongly depend on Q. The width of the instability band is then typically of the order [8] ∆k ∼ mQ which can be parametrically large for Q 1. Moreover, we note that in the broad resonance regime we typically have to take multiple instability bands into account, once the expansion of the Universe is considered. Luckily, the above rough expressions hold for all of them. In addition, we remark that the typical distance in Q between resonance bands of the exponent µ k (A k (Q), Q), for momenta of the order k ∆k, is The narrow and broad resonance regime can behave very differently with respect to the depletion of the DM energy density. Generally speaking, a broad resonance is typically more efficient, as more momentum modes are within the resonance band at the same time and the exponent is typically larger. That said, for our discussion the dependence of the rate of exponential growth, given by µ k , as well as the width of the resonance ∆k on the parameter Q is important. In the narrow resonance regime, the former is parametrically given by µ k ∼ Q while for a broad resonance it is typically of the order µ k ∼ 0.15, largely independent of Q. The bandwidth of both regimes is also considerably modified, i.e. ∆k ∼ mQ for a narrow and ∆k ∼ mQ 1/4 for a broad resonance, respectively. As we will see in the following section, both aspects are crucial for an explosive production of photons.

Including the expansion of the Universe
The previous discussion of the dynamics of the Mathieu equation only applies to a static situation. However, in the early Universe the expansion cannot be neglected. In this case, the parameters A k and Q explicitly depend on the scale factor. Strictly speaking, the concept of (static) resonance bands then ceases to be meaningful. Nevertheless, if the changes are sufficiently slow (compared to the time-scale of the DM oscillations) we can still get a reasonable picture by imagining the movement of a given momentum mode k along the trajectory (A k (t), Q(t)) through the instability chart. For instance, the system might start to evolve within a broad resonance, but as Q decreases with time, it eventually ends up in a narrow resonance regime before it terminates.
In a situation where Q may change significantly between consecutive oscillations of the driving DM field φ, one would instead have to move from a parametric resonance to a so called stochastic resonance [8]. Here, with a single oscillation of φ the phase of each photon mode is drastically altered, such that they are practically uncorrelated at any stage of photon production. Therefore, due to interferences, the number of photons produced typically increases but can also decrease with progressing DM oscillations, thereby slightly reducing the efficiency of the resonance. For a detailed discussion of the phenomenon of a stochastic resonance in an expanding Universe, we refer the reader to [8].
Such a thorough treatment of stochastic resonance is beyond the scope of this work. We therefore follow the more intuitive approximate approach outlined above, which has already been pursued in [10]. For our purpose, the most important difference between the static and the dynamical situation is that the photons experience a redshift as the Universe expands. Mathematically, the Mathieu parameter A k directly depends on the physical momentum of the mode, ∝ k/a, which changes with time. Each mode therefore only spends a finite amount of time in the resonant region. This can prevent the DM from efficiently annihilating into photons, if the latter are shifted out of a resonance quickly enough.
After a short amount of time δt, the momentum of a photon mode is shifted by [10] δk k Hδt , where δt is thought to be differential on cosmological scales, δt H −1 . That is, the momentum modes of the photons can only grow exponentially for a short amount of time, δt exp ∼ 1/(2µ k m), before they get shifted out of the resonance band, δt exp δt. While this is a universal feature that eventually terminates the resonant production of gauge bosons, its physical manifestation within the Mathieu dynamics has to be established carefully. In particular, there can be differences between a narrow and a broad resonance due to the significant modifications of the instability chart in these regions. Hence, we will discuss both scenarios separately.
a. Narrow resonance In a narrow resonance the Floquet instabilities occur at integer values of k/m (see also Fig. 3). In this regime, the lowest instability band, corresponding to the momentum k * m, is the most effective and we focus on this in our analysis. Hence, effectively, there is only a single resonance band that can induce an exponential growth of the photon modes. At the same time the expansion of the Universe can redshift the photons out of this instability, thereby preventing their resonant enhancement. That is, naively, there is a competition between the characteristic time of exponential growth determined by the Floquet exponent and the time the photon modes spend inside the resonance band. Reversing this argument, it can be written as a naive condition to avoid the rapid production of photons in the DM background, Obviously, this requirement is time dependent. Loosely speaking, the condition has to be satisfied at all times in the narrow resonance regime in order to avoid the complete fragmentation of the DM field and thereby to guarantee the cosmological stability of the DM candidate. The stability condition (16) so far only takes into account the growing exponential factor of the photon mode functions. However, as there can still be a small prefactor in front, a single short burst of rapid photon production may not be sufficient to trigger a complete annihilation of the DM field. Instead, the latter is only effective, if a significant amount of energy is transferred from the DM to the photons, i.e. if ρ A /ρ φ grows sufficiently that it becomes of order unity 8 . This, in turn, can be used to obtain a more precise condition for its cosmological stability. While the DM energy density schematically reads ρ φ m 2 Φ 2 /2, the energy density of the photons can be obtained by summing over all modes 9 where ω k denotes the energy of each momentum mode. Indeed, if the sub-exponential prefactor of n k in (9) is small, the resonance needs to be active for a considerable amount of time to transfer a significant amount of energy from the DM field into photons, In practice, the factor ζ depends on the initial conditions associated to the photon mode functions. Following [10], the sub-exponential correction can be estimated via a saddle-point approximation of (17). It schematically reads where n 0 is the initial occupation number of the photon modes, which can be determined by vacuum fluctuations or CMB photons, n 0 = 1/2 or n 0 2T CMB /m, respectively. The prefactor ζ can then be chosen conservatively, i.e. corresponding to the larger of both options.
Finally, requiring that the photons are redshifted out of the resonance quickly enough to avoid a complete fragmentation of the DM field yields a condition for the stability of the DM candidate, Using ∆k/k * = Q and µ k * = Q/2, this condition reads ζ (m/H)Q 2 (see also [8]). Therefore, in a narrow resonance, to avoid a fragmentation of the DM field its electromagnetic charge q has to satisfy In principle, to guarantee stability this inequality must hold at all times of the cosmic evolution. However, to ensure consistency, the system has to be in the narrow resonance regime, Q 1. In general, this is not necessarily the case. For example, as typically ζ ∼ 10 − 100, this condition is not consistent with the narrow resonance regime for m/H ∼ O(1) [8]. It is therefore worthwhile to also consider the broad resonance regime.
b. Broad resonance In a broad resonance the instability chart of the Mathieu equation is more complicated. In contrast to the case of a narrow resonance, the instability bands are not sharply localized around integer values of k/m. Instead, for a fixed Q, they can extend from a typical scale of k * ∼ mQ 1/4 down to possibly even vanishing momentum (see also Fig. 3). This means that, as the Universe expands, a given momentum within a certain instability band is not redshifted out of a single resonance, but, because Q decreases similarly as Q ∼ a −3 , it can cross multiple instability bands before it enters the regime of a narrow resonance. Therefore, the photon mode can experience multiple resonant enhancements on its trajectory through the instability chart. As an approximation we can model this by summing up all resonances that a given k-mode crosses, Since the Mathieu parameter Q is now a function of time, in a radiation-dominated Universe the exponent can be written as As pointed out above, within the instability bands for the range of momenta k ≤ ∆k ∼ mQ 1/4 , the Floquet exponent takes typical values of µ k ∼ 0.15−0.28, which are mostly independent of Q to good approximation. In this case, the bands can be assumed to be of width and distance of order √ Q in Q. Hence, the integration of µ k is dominated by small values of Q. As a rough approximation we can therefore choose the integration boundaries to be Q − = 1 (beginning of the broad resonance region) and Q + = ∞. Evaluating this numerically we find with monotonically decreasing values within those limiting cases. Similar to the case of a narrow resonance, the fragmentation of the DM field is efficient (ultimately leading to a breakdown of our description of the DM oscillations), if a significant fraction of energy is transferred from the DM to the photons, ρ A /ρ Φ ∼ 1. Again, this requirement can be translated into a stability condition for the DM candidate, where we take account of modes up to k = mQ 1/4 in the energy density. This is the broad resonance equivalent to the stability requirement (21), which is applicable in the narrow resonance regime. Note that in determining this expression we have made several rough approximations. As already stated above we have neglected the dependence on the upper integration boundary in (24). This ignores contributions suppressed by an inverse power of Q. We have also used that, due to the exponential growth in the modes, the biggest drain in energy occurs at late times and therefore evaluated the energy drain only at the end of the broad resonance regime. Moreover, using a fixed κ we have neglected that during the evolution the physical momentum of each mode decreases as k ∼ a −1 and therefore the exponent changes. Finally, in line with all these approximations we have simply dropped terms logarithmic in Q.
c. Discussion The stability conditions (21) and (25) put strong constraints on the value of the electromagnetic charge of the DM. Before explicitly evaluating them, let us first get some analytical understanding. In principle, in order to avoid the resonant depletion of the DM, both For comparison, the dotted blue line shows the narrow resonance condition evaluated close to the earliest possible time, 1000 t * , while plasma effects are neglected, mA = 0. The observational constraints are given by CMB observations [83], SN1987A [84] and stellar cooling (SC) [85], pulsar timing arrays [86] or by interactions with magnetic fields in galaxies and clusters [87,88]. In these limits, the dashing indicates regions where we have used very naive extrapolations from the high-mass regime.
conditions have to complement each other such that either one of them is satisfied at all times of the cosmic evolution. As we have pointed out before, we expect the fragmentation of the DM field to first be governed by a broad resonance. Then, as Q decreases with time in an expanding Universe, Q ∼ a −3 , the system will enter the narrow resonance regime before the fragmentation eventually terminates (given that the coherent field is not completely destroyed at that point). Therefore, in practice, one has to carefully establish which stability condition gives the correct, i.e. self-consistent, constraint on the millicharge at each time. This depends on the charge as well as on the time when the stability condition is evaluated.
In general, it is a priori not at all obvious, what time gives the strongest possible constraint on q. In fact, both regimes (21) and (25) behave very differently with the scale factor. To see this, we insert the evolution of the DM amplitude, Φ ∼ a −3/2 , and the behaviour of the Hubble constant during radiation domination, H ∼ a −2 , into the exponential factor (i.e. the right hand side) of both stability conditions. For a narrow resonance we obtain m H qΦ 2m while, in contrast, the broad resonance regime behaves as This suggests that in the narrow resonance regime the strongest constraint arises when evaluating at the earliest possible time, whereas the broad resonance case appears to be independent of the scale factor. Let us consider both scenarios. In case of a narrow resonance, Eq. (26) suggests that we should evaluate the stability condition when φ just starts to oscillate, i.e. at t * when H m. However, as noted before, fulfilling (21) with ζ ∼ 10 − 100 is inconsistent with the narrow resonance regime at t * . Hence, we either have to evaluate at a somewhat later time when H m or go into the regime of a broad resonance. For the strongest selfconsistent constraint on the millicharge in the narrow resonance regime we can evaluate at t ≈ 1000 t * .
In contrast, in the broad resonance regime, Eq. (27) suggests that the constraint on the millicharge is independent of the time when the stability condition is evaluated. This suggests that the limit does not strengthen much when approaching the broad resonance regime 10 . Indeed, we have checked that (25) roughly provides for the same limit as (21) evaluated at t ≈ 1000 t * .
A numerical evaluation of the stability condition (21) at t 1 ≡ 1000 t * is shown as the dotted blue line in Fig. 1. However, this estimate is probably too optimistic as we still need to include plasma effects, which we will do next.

B. Photons inside the early Universe plasma
So far, we have assumed that the photons are moving freely through the Universe. However, during the cosmological evolution, the early Universe is filled with a hot plasma that modifies their propagation. Indeed, for example in the case of axion-like particles, this is the dominant effect ensuring their stability [1,2,10]. Therefore, it is sensible to also consider this effect for the case of millicharged DM (a discussion of parametric resonance in charged cosmological scalars can also be found in [56] but they were focused on primordial magnetic fields rather than DM).
Naively, the photons interact with the charged particles of the medium such that they acquire a modified dispersion relation (and wavefunction renormalization), see e.g. [89]. Effectively they acquire a mass, 10 We note that the complete independence of the evaluation time is, of course, due to the simplistic approximations we employ. m A . In Eqs. (6) and (7) this leads to the replacement k 2 /a 2 → k 2 /a 2 + m 2 A . Therefore, A k is larger and, if the plasma mass is too high, the resonance becomes inefficient. In particular, in a narrow resonance regime the instabilities become ineffective if the plasma mass exceeds the mass of the DM candidate, m A m. For nearly all k the rate of exponential growth µ k becomes imaginary, corresponding to an oscillating rather than a growing mode function. In contrast, in a broad resonance, the production of photons with masses m A m is in principle possible 11 . However, this process requires comparatively large couplings in general. In particular, photon production in this regime can only be efficient for charges satisfying [8] Overall, we therefore expect that this possibility will result in a weaker constraint on the electromagnetic millicharge (see also the example below). The plasma mass of the photon depends on the temperature of the medium. That is, in an expanding Universe, it is time dependent 12 . As noted above, in practice, the condition m A m sets the earliest time at which the (narrow resonance) stability condition (21) can be evaluated and turns out to stabilize the scalar DM candidate in large parts of parameter space of the vanilla theory of millicharged DM. We show this as a solid blue line in Fig. 1.
In the broad resonance regime, the production of photons with masses above the DM mass is possible, but for this, larger charges are needed. In particular, as a minimal requirement, the broad resonance must be strong enough to overcome the mass threshold. This requires fulfilling the condition (28). This is usually already a weaker requirement than evaluating the narrow resonance stability condition at the point at which the photon mass is small enough for the narrow resonance to be active. As an example, this is demonstrated by a dashed and a dashdotted blue line in Fig. 1, where we choose times when the plasma mass is of the order of m A /m = 1 and m A /m = 100, respectively.
Looking at Fig. 1 we can see the drastic impact of the plasma effects. Comparing the naive estimate that completely neglects plasma effects (dotted blue line) with the constraint taking into account the plasma effects (solid blue line), the former turns out to be many orders of magnitude stronger. Indeed, we observe that in the regime of very low masses, the stability condition on millicharged DM is a weaker requirement than current observational 11 We thank Paola Arias, Ariel Arza and Diego Vargas for very useful discussions (triggered by the helpful comments of an anonymous referee for [90]) on this issue in a similar system. 12 We use the cosmological evolution of m A as given in [10] which is based on [89,[91][92][93][94].
constraints [83][84][85][86][87][88] 13 . While this is desirable from a physical point of view, the simple model we have considered in this section is disfavoured from a theoretical point of view, as quantization of electromagnetic charge would be hard to justify. Let us therefore turn to a more realistic and appealing theory involving a hidden photon that is kinetically mixed with the visible sector.

III. MILLICHARGED PARTICLES ARISING FROM A KINETIC MIXING OF A HIDDEN PHOTON WITH THE VISIBLE SECTOR
From a theoretical point of view, a fundamental millicharge is unappealing with respect to charge quantization. A well motivated alternative is provided by kinetic mixing [38]. A simple example of this scenario is an additional massless hidden photon X µ that is kinetically mixed with the electromagnetic photon A µ [38], For a hidden sector matter particle φ (in our case the DM candidate) carrying a (quantized) charge q under X µ , a small effective electromagnetic charge appears after diagonalizing the kinetic term. To see this, one can rotate the gauge fields by A µ → A µ and X µ → X µ + A µ . While this redefinition leads to canonically normalized kinetic terms of the gauge fields, it also appears in the gaugecovariant derivatives of the hidden sector matter field, where g is the hidden sector gauge coupling. As a consequence, the DM carries an effective electromagnetic charge [38] q eff = qg , where we have again absorbed the factor of the electromagnetic coupling e into the charge. In this way, a small (and possibly also g) can lead to a tiny electromagnetic charge, even if q is integer 14 .
In general, the amount of kinetic mixing is a free parameter and may even be of order one. However, if we consider the hidden photon to be part of a hidden sector we usually expect that the mixing is small. For instance, the hidden gauge group may be understood as a low-energy remnant of a UV theory with a unified gauge symmetry broken at some high scale [48]. After symmetry breaking, some degrees of freedom, for instance a heavy fermion, usually carry a charge both under electromagnetism as well as the hidden gauge group. Quantum mechanically, a kinetic mixing between both gauge fields is then induced by a fermion loop of the UV theory. At low energies, the kinetic mixing parameter is determined by the corresponding one-loop Feynman diagram and parametrically reads (see, e.g., [38,48]) where m ψ is the mass of the heavy fermion and µ is the regularization scale of the loop integral. This typically gives a small kinetic mixing, which is particularly tiny if also the hidden sector gauge coupling is small, g 1.
We can now apply the arguments of the previous section to this scenario. As we assume no plasma to be present in the hidden sector, we expect that annihilation into hidden photons remains possible and therefore guaranteeing stability may put stronger constraints on the millicharge.
We will examine this scenario in two separate steps. First we completely neglect the (small) kinetic mixing effects, i.e. we consider a secluded hidden sector without kinetic mixing, = 0. Then we argue that the main conclusions also hold in the phenomenologically more interesting case with a small but non-vanishing kinetic mixing parameter.

A. Secluded hidden sector
In the case = 0, our discussion in Section II completely carries over. In particular, the DM stability conditions (21) and (25) in the narrow and the broad resonance regime can be applied at all times of the cosmic evolution. Most importantly, as there is no effective mass of the hidden photon that could block the resonant enhancements, they can in principle be satisfied at the earliest possible time. As discussed in the previous Section II, a reasonable estimate is obtained by evaluating the stability condition from the narrow resonance regime at t 1 ≈ 1000 t * . This is shown as the solid blue line in Fig. 2, where we display the allowed value of g as a function of the DM mass m. (Note that here, we have normalized the field to unit charge, q = 1.) The requirement of avoiding a resonant depletion of the DM energy density into hidden photons puts severe constraints on the hidden gauge coupling for small masses. Allowed hidden gauge coupling g to the scalar DM candidate as a function of its mass m. The blue lines correspond to the stability condition due to the parametric resonance, evaluated close to the time where the field starts to oscillate, t1 ≈ 1000 t * , (solid) and at matter-radiation equality (dashed). The red line is given by a coherence condition, discussed in the main text. For comparison, the right axis shows the typical corresponding effective millicharge induced by a fermion-loop of a UV theory, q eff ∼ g ∼ eg 2 /(6π 2 ). Along the same lines, the light-shaded grey area in the upper region corresponds to observational constraints provided by interactions with a magnetized intergalactic stellar medium [87,88], see Fig. 1.

Backreaction effects
The above evaluation might be an overestimation of the constraint posed by the DM stability requirement. This is because, so far, we have neglected the backreaction of the parametric resonance on the DM field. Obviously, a first effect is the depletion of the DM field. This is what we have implicitly used to set our constraint, i.e. using energy conservation to determine the depletion from the produced gauge bosons. However, if the energy density in the hidden photons is comparable to that in the DM field, ρ A ∼ ρ φ , it is conceivable that energy starts to be transferred back to the DM field, slowing down the depletion. Although this is non-trivial due to the fact that most of the produced hidden photons have momenta k m, processes involving multiple hidden photons may be possible due to the high occupation numbers and the resonantly enhanced interaction rates. With the present analysis we cannot exclude this possibility. A thorough analysis of this effect would need to involve some careful numerical simulations, which is beyond the scope of this work. The overall allowed value of the hidden gauge coupling might therefore be higher.
That said, let us obtain a very conservative estimate of the point where the resonance should shut off. Allowing for the backreaction effect, we nevertheless expect that in such a situation a significant fraction of the total energy in the DM-hidden photon system, possibly ∼ 1/2, will be in hidden photons and therefore in the form of dark radiation. At matter-radiation equality such a large fraction of dark radiation is certainly excluded (cf., e.g. [97]). Therefore, we can evaluate the DM stability condition at matter-radiation equality. At this stage, at the latest, φ is required to behave as standard cold DM. At the same time, we expect the system to be in a narrow resonance regime, such that the stability condition (21) is valid. The resulting constraint on the hidden gauge coupling is shown as a dashed blue line in Fig. 2. It is considerably weaker than the original estimate, but still affects an appreciable region of parameter space, bearing in mind that this estimate is probably overly conservative.
This has been discussed in detail in [90] (see also [17,22,24,29,30] for discussions in the context of DM structures) from which we summarize the main implications. In order to preserve coherence of the hidden photons produced by the resonance, the width of the resonance in momentum space has to be larger than the momentum spread of the DM, ∆k ∆k φ . The latter can be estimated to be ∆k φ ∼ mv mr (a mr /a), where we require that φ should be non-relativistic at matter-radiation equality, v mr ∼ 10 −3 (cf., e.g. [98][99][100][101]). Therefore, the condition for preserving coherence can be written as ∆k mv mr a mr a .
As pointed out in Section II, the width of the resonance bands in momentum space, i.e. the left-hand side of this inequality, depends on the value of the Mathieu parameter Q. Evaluating (33) at matter-radiation equality we see that the required width is much smaller than the mass m. Therefore, we can use the narrow resonance regime where ∆k ∼ mQ. This can be immediately translated into a constraint on the hidden gauge coupling, which we similarly evaluate at matter-radiation equality. This is shown in red in Fig. 2.

Discussion
In general, our results, shown in Fig. 2, demonstrate that the stability requirement for a very light DM candidate in a secluded hidden sector affects sizeable regions of parameter space. The strongest constraint is posed by the parametric resonance stability condition evaluated close to the time, when the DM field starts to oscillate (solid blue). However, this neglects backreaction effects and therefore needs to be taken with caution. A more conservative estimate is given by evaluating the stability condition at matter-radiation equality (dashed blue). We expect that a more careful numerical analysis would most likely reveal a stability condition that lies in between those possibilities. Aside from backreaction effects, additionally, a non-trivial initial velocity distribution of the DM particles, possible in some models for their production, may allow to weaken the constraint as the resonance requires a sufficient amount of coherence. Taking into account that the DM velocity must be small enough to allow for successful structure formation, we find the most conservative constraint shown as the red line. Intriguingly, we are still able to probe large regions of parameter space. In fact, it seems challenging to motivate or construct models featuring gauge groups with such tiny gauge couplings. A famous example allowing for small gauge couplings is provided by the large volume scenario (LVS) of type IIB string compactifications [102,103]. Here, the gauge theory is supported on D-branes wrapping cycles of the internal Calabi-Yau manifold of ten-dimensional spacetime. The volume of these internal cycles, in turn, determines the gauge coupling, g ∼ V −1/3 [104]. Therefore, hyper-weakly coupled gauge theories can be engineered by choosing an appropriate Calabi-Yau geometry that supports large D-brane worldvolumes. In Fig. 2, we show typical values of the gauge coupling achieved in a generic (orange) and low string-scale, M S ∼ 1 TeV, (green) LVS [104]. Strikingly, even the extreme case of TeV-scale strings is rendered unstable for a wide range of masses even using the most conservative constraint. We note, however that the bound provided by the weak gravity conjecture [105], shown in grey, is not reached. It would therefore be interesting to see whether consistent models with such small gauge couplings can be constructed.
As we will argue in the following section, similar conclusions still hold, if the hidden sector is not completely secluded but has a small kinetic mixing parameter connecting it to the visible world.

B. Non-vanishing kinetic mixing
Along the lines discussed at the beginning of the section we focus on a situation with small kinetic mixing. In a homogeneous background ϕ, the equations of motion are linear in the photon A and hidden photon X and therefore couple distinct momentum modes of both fields. In fact, after having redefined the gauge fields by A µ → A µ and X µ → X µ + A µ , their mode functions satisfÿ Here, we have already included an effective mass for the photon, m A . Naively, the equations of motion imply that both the photon as well as the hidden photon modes can be enhanced by a parametric resonance induced by the oscillating DM background. However, a resonant enhancement of A is now parametrically weaker as compared to X, because its coupling contains an additional factor of the mixing parameter . This means that, for instance, modes of the hidden photons might be growing rapidly due to a broad resonance, while the photon modes already are in a very narrow resonance regime and not amplified efficiently. At the same time, this amplification might also act as an oscillating driving force on the right hand side of (34). Eventually, the growing modes of both fields will converge to the same resonance frequency after a certain period of time. Therefore, in general, the DM may largely annihilate into hidden photons and also, to a smaller fraction, into visible photons which follow shortly after. The above observations suggest that an efficient depletion of the DM energy density is possible in a theory featuring kinetic mixing. In practice, this is important, as the visible photon can obtain a non-negligible plasma mass, m A = 0, while the hidden photon is still massless. However, as the DM mainly annihilates into hidden photons, we are able to avoid plasma effects of the visible photons in the early Universe almost entirely. This is illustrated in Fig. 3 where we compare the instability chart of a completely secluded hidden photon (blue) with that for non-vanishing kinetic mixing, = 0.1 and m A /m = 100 (orange points denoting the boundary). We obtain these by numerically solving the coupled equations of motion for X and A for different momenta. As an example, we choose a DM mass of m = 10 −3 eV and both instabilities are evaluated when the DM field starts to oscillate, t * . Inside the instability bands, an exponential growth of the momentum modes of X, i.e. rapid production of hidden photons, is possible. We can see that the unstable regions are almost identical. The plasma mass in the visible sector does not prevent the resonant annihilation of DM into hidden photons. We expect this to be true everywhere in parameter space for kinetic mixing parameters smaller 16 than 0.1. Therefore, the allowed values of the hidden gauge coupling are almost identical to what is shown in Fig. 2. To give an impression of the constraints of the effective millicharge, we show on the right-hand axis of Fig. 2 indicative values of the millicharge obtained by combining Eqs. (31) and (32). The light grey region indicates the experimental and observational constraints on the effective charge as also shown in Fig. 1. There are large regions where even our most conservative estimate of the unstable region poses a stronger constraint on the effective millicharge than current observational bounds.

IV. CONCLUSIONS
The microscopic nature of dark matter (DM) that comprises large parts of the cosmic fabric remains elusive. As suggested by its name, so far, there is no experimental evidence of DM interacting with electromagnetism. While this naively rules out any sizable electric charge assigned to DM particles, it is still possible that their charge is tiny, thereby strongly suppressing interactions with photons. In this work, we have investigated the cosmological longevity of such DM particles in the sub-eV mass regime. In this mass range the DM particles must be bosonic and, for concreteness, we have chosen them to be scalar.
The millicharged particles are either minimally coupled to photons or their electromagnetic interaction is mediated via kinetic mixing with a massless hidden photon. In both cases, due to the large occupation numbers of the light DM field, even for tiny charges the DM may efficiently annihilate into gauge bosons via a parametric resonance [6][7][8][9].
We find that, in the case of a direct coupling to photons current observational constraints on the millicharge are stronger than those arising from parametric resonance, as shown in Fig. 1. This is mainly due to the plasma mass that the photons acquire in the hot medium of the early Universe, which essentially terminates the resonance if it is larger than the DM mass. In contrast, in the case of a theory featuring kinetic mixing, plasma effects are practically absent. Therefore, even employing conservative estimates large regions of parameter space are affected by the parametric resonance as illustrated in Fig. 2. In fact, in particular for very low DM masses, its electric millicharge has to be orders of magnitude below what has been typically obtained in UV models, e.g. in type IIB string compactifications in the large volume limit [102][103][104] (we note that, indeed, already the experimental and observational constraints rule out this region of parameter space for loop-induced kinetic mixing). That said, the limits do not yet reach the smallest possible values suggested by the weak gravity conjecture [105], which are many orders of magnitude below the smallest values found in the concrete realizations discussed above.
We conclude that it is far from trivial that very light bosonic DM carrying a tiny electric charge is long lived. Instead, its high occupation numbers can dramatically enhance interaction rates with gauge bosons leading to parametric resonance phenomena. Therefore, its cosmological stability cannot be taken for granted but has to be considered carefully.