Relaxion Dark Matter

We highlight a new connection between the Standard Model hierarchy problem and the dark matter sector. The key piece is the relaxion field, which besides scanning the Higgs mass and setting the electroweak scale, also constitutes the observed dark matter abundance of the universe. The relaxation mechanism is realized during inflation, and the necessary friction is provided by particle production. Using this framework we show that the relaxion is a phenomenologically viable dark matter candidate in the keV mass range.

duced in Ref. [2]. In this case, the barriers Λ 4 b cos(φ/f ) are constant, and the field's kinetic energy is large enough to overcome them. By assumption, the EW symmetry is broken early on, and all the SM particles are initially very heavy. When the relaxion approaches the critical point at which the Higgs VEV is zero, all SM particles become light. At this point, the relaxion stops due to the production of SM gauge bosons, due to a coupling that we will discuss below in Sec. II. The EW scale can be related to the parameters of the model, with v EW Λ in a technically natural way.
The cosmological consequences of the relaxion scenario have not been fully explored yet. This is of the maximum importance, as it would help both in constraining the properties of the relaxion field and in pointing out possible observable signatures of its existence. In this note, we address the question of whether the relic population of relaxion particles can constitute the current DM density. In previous literature, the answer was negative. 2 The constructions discussed in Ref. [1,3,41] use a Higgs-dependent barrier to stop the field evolution, which happens during inflation. In these models, the relaxion abundance is negligible [41]. A second field, which scans the barriers' amplitude in Ref. [3], can instead have a sizable misalignment yield. Oppositely, if relaxation proceeds after inflation, and the source of friction is the tachyonic production of gauge bosons, the relaxion is overproduced [20], and one has to impose that its lifetime is short enough to dilute this abundance before nucleosynthesis.
Here we assume that relaxation happens during inflation and the relaxion stopping mechanism is provided by particle production [2]. This construction does not require new physics close to the TeV scale, and, in a portion of the parameter space, it can be realized without a large number of e-folds or super-Planckian field excursions. As we will detail below, in this scenario the relaxion particles, produced after reheating by scatterings in the SM plasma, can account for the observed DM density.
The paper is structured as follows. In Sec. II we introduce the model and discuss in details the mechanism to generate a small EW scale. In Sec. III we discuss the conditions that need to be applied on the parameters of the model. Section IV discusses the production of the relaxion DM population in the early universe, whose properties are discussed in Sec. V. Finally, we draw our conclusions in Sec. VI.

II. RELAXATION WITH PARTICLE PRODUCTION
We now introduce the relaxion model that we will consider throughout this paper, which was first introduced in Ref. [2]. We will assume that relaxation takes place during inflation, and refer the reader to Ref. [20] for an analysis of the case in which relaxation happens after inflation. The Lagrangian is [2]: where φ is an axion-like field with decay constant f , h is the Higgs field, Λ is the cutoff of the theory, the dimensionless parameters g and g are assumed to be spurions that explicitly break the axion shift symmetry and λ is the Higgs quartic coupling. The scale Λ b is related to the confinement scale Λ c of some non-abelian gauge group, at which the φ cosine potential is generated as Λ 4 b ∼ Λ 4 c , where in general 1. For example, if the relaxion is the QCD axion, ≈ m q /Λ QCD , where m q is the mass scale of the up-and down-quark. In the following we will not specify any further the mechanism responsible for the generation of these barriers. The effective scale F controls the interaction of the relaxion with the SM gauge bosons. We assume that the relaxation dynamics takes place in the broken phase so that the Higgs mass parameter, , is large and negative when the scanning process starts, µ 2 h (φ ini ) ∼ −Λ 2 . The potential of Eq. (1) also induces a mixing of the relaxion with the Higgs [41]. After relaxation ends and the relaxion stops in one of the minima of its potential, the mixing angle is The last term in Eq. (1) is responsible to slow down the relaxion once the particle production is triggered. B and W are the SM gauge bosons with g 1 and g 2 being the corresponding U(1) and SU(2) gauge couplings. When expanded in the mass eigenstates this term reads In what follows we will only consider the tachyonic instability from the Z Z term and absorb the gauge coupling in the definition of the corresponding field such that 1/f = (g 2 2 −g 2 1 )/F. The contribution from the W W term is expected to be suppressed due to the self-interactions of the W , which induce an effective mass making particle creation inefficient. Tachyonic production of photons, as we will see, is suppressed in this model, thus we expect the ZA term in Eq. (3) to be subdominant with respect to the ZZ one.
As it is well-known in the context of axion inflation (see, e.g., Ref. [46] and references therein for a recent account), a coupling of the form φF F withφ = 0 leads to the exponential production of gauge bosons. In the case at hands, this production is suppressed by the mass of the Z, which is initially large. As we will see below, particle production starts only when m Z Λ. At this point it should be clear that a coupling to photons φF F must be suppressed, because it would lead to exponential photon production during the relaxion's evolution, thus slowing down the field independently of the Higgs VEV and spoiling the mechanism. The non-generic coupling structure of Eq. (3) is designed for this purpose. This can, for example, descend from a left-right symmetric UV completion [47]. Also, it appears in non-minimal composite Higgs models with coset SO(6)/SO(5) [48][49][50][51][52], even an attempt of embedding the relaxion mechanism in such a construction is still missing in the literature.
The interaction of Eq. (3) generates a coupling to SM fermions (at one loop) and to photons (at one and two loops) [47,53]: where N F c is the color factor, Q F is the electric charge of the fermion F with mass m F , and x i ≡ 4m 2 i /m 2 φ . The functions B 1,2 are: When the relaxion is light, which will turn out to be the case of interest for our DM scenario, these functions scale This implies, for instance, that when the relaxion is lighter than the electron mass, the induced coupling to photons originated from the coupling in Eq. (1) is suppressed.
Let us now describe how the relaxion mechanism works in this model. The φ field rolls down its potential until it reaches the critical point when there is an exponential production of gauge bosons, which makes the relaxion slow down due to the transfer of its energy to the gauge bosons. This back-reaction mechanism becomes apparent once we examine the equations of motion for φ and Z: where m(h) = g 2 1 + g 2 2 h/2, Z Z is the expectation value of the quantum operator and Z ± refers to the two transverse polarizations of Z µ . 3 In terms of the mode functions Z ± , Z Z can be written as Assumingφ > 0, Z + has a tachyonic growing mode when ω 2 k,+ ≡ k 2 + (m(h)) 2 − kφ f < 0. The first mode that becomes tachyonic is the one for which ω k,+ is minimum, k c =φ/(2f ), thus Z + grows exponentially foṙ The growth of the mode function Z + continues until Z Z becomes the dominant term in the equation of motion of φ, Eq. (9). After this point,φ ≈ − Z Z /(4f ) < 0 and the relaxion velocity decreases. After φ has slowed down, the constant cosine potential acting as a barrier can make the relaxion evolution stop. An example of such evolution was computed numerically in Ref. [20]. The Z mass depends on the Higgs field, implying that condition (12) is not satisfied when the Higgs field value is large and particle production is ineffective. To obtain the correct value of the electroweak scale, the back-reaction should be triggered when m(h) = m Z ≈ 90 GeV, i.e., foṙ In the following we will use this equation to write f as a function of the model's parameters. The back-reaction must turn on when φ is close to the critical value Λ/g that cancels the Higgs mass term in Eq. (1), generating a parametric hierarchy between the cutoff Λ and the electroweak scale. This happens if the classical Higgs field h follows closely the minimum of its potential, as we will detail in the next section. After the tachyonic growth starts, the Higgs and gauge fields undergo a complicate dynamics. First of all, we expect the system to quickly thermalize as soon as the energy density of the gauge fields becomes larger than the EW scale [2] (see also Ref. [20]). The interaction rate can be estimated as with ρ ∼φ 2 > v EW , thus Γ int > v EW , which should be compared with the tachyonic growth rate which is O(v EW ). Notice that, in the following, we will impose that the timescale for particle production is much shorter than the Hubble time.
Finite density effects affect both the Higgs and the gauge fields' evolution. A positive mass term for the Higgs is generated in the thermal bath, temporarily restoring the electroweak symmetry. The field h (and thus the mass of the gauge bosons) rolls to zero, making the tachyonic growth faster.
On the other hand, the presence of the thermal plasma affects the dispersion relation of the gauge bosons, which is modified into where, in a hard thermal loop (i.e., high temperature) limit [54], Here m 2 D = g 2 EW T 2 /6 is the Debye mass of the bosons in the plasma. The factor g EW is obtained by taking into account all the SM fermions and their respective hypercharges, as g EW ≈ (32/9) sin 2 θ W g 2 1 ≈ 0.2, where g 1 is the SM hypercharge coupling and the Weinberg angle sin 2 θ W projects the Z onto its abelian component.
The function Π[ω, k] is positive for imaginary frequency ω = iΩ, thus damping the tachyonic instability. Expanding Eq. (16) for Ω/k → 0 we obtain Ω is thus maximized for k = 2φ/(3f ), with From Eq. (18) we can estimate the typical timescale for the exponential growth as The temperature T is obtained by assuming that the relaxion kinetic energy is transferred to radiation: To summarize, we expect that particle production leads to the production of a thermal bath of SM particles, which temporarily restores EW symmetry and, at the same time, reduces the particle production rate. After the relaxion field stops, the temperature is rapidly erased by cosmic expansion, and the Higgs relaxes to its zero temperature VEV which is now fixed by the value of the relaxion field. Equation (13) ensures that the final Higgs VEV is the measured one.

III. PARAMETER SPACE
The parameter space is characterized by six parameters: the cut-off Λ, the couplings g and g , the barriers' height Λ b , the decay constant f , and the Hubble constant during inflation H I . The scale f can be fixed in terms of the other parameters using Eq. (13) f =φ/(2m Z ) and the value of the slow-roll velocitẏ φ = gΛ 3 /(3H I ), yielding The couplings g, g must satisfy g > g /(4π) 2 , otherwise a linear term g Λ 3 φ/(4π) 2 , generated through a Higgs loop, would dominate over the term gΛ 3 φ in the relaxion potential. Since the two spurions may be generated in a similar manner in the UV theory, it would be reasonable to assume g = g . Still, as we will show below, it is convenient to relax this assumption, and therefore in the following we will consider the benchmarks g/g = 1, 10 3 , 10 6 . Independently of the relaxion being the DM, there are a number of conditions that the model must satisfy to actually solve the hierarchy problem. First, the relaxion should not affect the inflationary dynamics, implying that the relaxion potential is subdominant compared to the inflaton one. This gives a lower bound on the inflation scale H I : In addition, the assumption that φ evolves classically is valid only if the classical evolution dominates over the quantum fluctuations during inflation. Therefore we impose that, over a Hubble time, (δφ) class (δφ) quant with (δφ) class ∼ V φ /(3H 2 I ) and (δφ) quant ∼ H I /(2π). This gives us an upper bound on the inflation scale, where we used that V φ ∼ gΛ 3 . Furthermore, inflation should last long enough such that the relaxion has time to scan the Higgs mass parameter. The minimal number of e-folds which is required to scan a field range ∆φ ∼ Λ/g , is given by where in the last step we used that the slow roll velocity isφ ∼ 2m Z f . We also need to make sure that the Higgs field is efficiently tracking the minimum of its potential during the scanning process. This ensures that the back-reaction from the exponential production of gauge bosons is triggered when the VEV is at the electroweak scale. Hence, where v = (Λ 2 − g Λφ) 1/2 / √ λ is the minimum of the Higgs potential given in Eq. (1). The need for this conditions can be understood as follows. The mass of the Z boson, and hence the time at which its tachyonic production starts, depends on the value of the Higgs field h. After relaxation is over, h will relax to the minimum v, so that v controls the current value of the EW scale. Therefore, it is important that the evolution of h and v match, otherwise the relaxion field would stop as soon as h satisfies Eq. (12), while having a value of v different from the current one. Eq. (25) needs to be satisfied until the Higgs field value has reached the electroweak scale.
Another necessary condition is that the average slowroll velocity during the scanning has to be large enough to overcome the barriers generated by the cosine potential in Eq.
being a contribution due to the cosine potential. At the same time, the average slow roll velocity should not exceed the cut-off, for the consistency of the effective theory: V φ /(3H I ) Λ 2 . The sharp cut in the green region of Fig.(1) at Λ 10 4 GeV descends from this condition, after fixing H I to get the correct relic abundance (see below for more details).
Additionally, once the back-reaction has turned on, the barriers must be high enough to stop the relaxion evolution, requiring that We should also ensure that once relaxion is slowing down, the Higgs mass does not change by an amount larger than the correct value, i.e., We impose that the kinetic energy lost by φ due to particle production is larger than the one gained by rolling down the potential, We estimate the two terms as ∆K pp ∼φ 2 /2 and ∆K rolling ∼ dK dt ∆t pp , where dK/dt = −dV /dt ∼ gΛ 3φ . To be conservative [see Eq. (26)], we takeφ 2 /2 ∼ Λ 4 b . On top of that, one should guarantee that the particle production is faster than the expansion rate, so that the energy dissipation efficently slows down the relaxion field. Furthermore, the scanning must have enough precision to resolve the electroweak scale. The mass parameter µ 2 h cannot vary more than the Higgs mass over one period of the cosine potential, In addition, it is crucial that the induced coupling to photons in Eq. (4) is suppressed enough, otherwise the dissipation from particle production would be relevant independently of the value of the Higgs mass. Then we have to impose that the produced photons are efficiently diluted by the cosmic expansion where ∆t γ ∼ T 2 f 3 γ /φ 3 with f γ given in Eq. (6). The relaxion induced coupling to photons through the Higgs mixing is very suppressed and the dilution requirement in Eq. (32) for this contribution is trivially satisfied.
A last condition concerns the restoration of the shift symmetry. After the relaxion has been trapped into one of the wiggles, the temperature cannot be larger than the confinement scale, T < Λ b , where, to be conservative, we assumed Λ b ≈ Λ c [see discussion below Eq. (1)]. This condition is only relevant if the sector which generates the cosine potential gets in thermal equilibrium with the SM model. This can be estimated as follows. We assume that the barriers are generated by some confining gauge group, which is coupled to the relaxion via a term (φ/f )G G . Then, we naively estimate the rate for g g ↔ ZZ interactions as Γ ∼ T 5 /(f 2 f 2 ), which must be larger than the Hubble rate H I to achieve thermalization.
All in all, the conditions that apply to the parameters of our model are the following: Higgs tracking the minimum (34) overcome the wiggles (35) small Higgs mass variation (37) The coloured region in Fig. 1 shows the values of Λ, g for which the relaxion mechanism can be realized successfully, for a fixed ratio g/g . To each point it corresponds a range in the other three free parameters. The colors correspond to different conditions which are imposed in order to make the relaxion a viable DM candidate, as we will detail in the next section.

IV. RELAXION AS DARK MATTER
The relaxion can be produced via vacuum misalignment and through thermal scattering. In the first case, after the relaxion gets stuck in one of the barriers, it will eventually start to oscillate freely when particle production becomes inefficient, leading to a energy density which red-shifts as non-relativistic matter. Since in our scenario the relaxation dynamics happens during inflation, the energy density stored in the field is diluted away and the misalignment contribution to the relaxion abundance is negligible [2] (see also Refs. [3,41]). The only possibility to produce a significant relaxion abundance is then via scattering.
A population of relaxion particles is produced through a + b ↔ φ + c interactions, where the species a, b, and c belong to the SM and are in thermal equilibrium. The relaxion abundance is controlled by the Boltzmann equa-tion [55] where Y φ = n φ /s and x = m φ /T with n φ being the relaxion number density and s the entropy density. The equilibrium number density of φ is Y eq φ = n eq φ /s ≈ 0.278/g * where g * is the number of relativistic degrees of freedom, and H is the Hubble rate. The quantity Γ is given by the sum over the interaction rates, Γ ≡ i Γ i with Γ i = n ci σv i , where the sum includes gluon scattering [56], Primakoff scattering (via φγγ, φZγ, and φBB, see, e.g., Refs. [41,57,58]), Compton scattering of leptons and quarks (via φll and φqq, see, e.g., Refs. [41,57]), and Primakoff and Compton processes through the mixing with the Higgs (see, e.g., Ref. [41]). The pion-conversion processes π 0 N → φ N and π 0 π 0 → φ π 0 can be neglected, as they are only active for a short time around the QCD phase transition, and their integrated rate is thus negligible. Assuming that the initial φ density is negligible, where T 0 = m φ /x 0 can be identified with the reheating temperature (we will comment more on this below). If the integrand is large, then Y φ (x) ≈ Y eq φ and the correct DM abundance can only be met for a very light (thus hot) DM component. We must therefore be in the opposite situation, in which the integrand is smaller than one, and we can approximate Y φ by All the processes listed above are suppressed by the same physical scale f (or by the mixing angle θ for the ones mediated by Higgs-relaxion mixing). Therefore, the integral in Eq. (48) is dominated by the process that is active for the longest time, which is Compton scattering off electrons γ + e ↔ φ + e. Using the expressions for the interaction rates Γ i in Refs. [41,[56][57][58], we checked explicitly that all the other processes give a subdominant contribution to the relic abundance. We can also neglect all the interference terms, since at each temperature the subdominant processes are highly suppressed compared to the main one. Compton scattering is active until electrons become non-relativistic. This is why we need a rather low T 0 , which guarantees that the interactions are out-of-equilibrium (Γ i /H < 1) and that the relaxion never enters in thermal equilibrium with the SM bath. The interaction rate is given by where f e is given in Eq. (5). The relaxion decays through the loop-induced couplings to photons and SM fermions as in Eq. (4), by the leading interaction with the electroweak gauge bosons in Eq. (1), and via the mixing with the Higgs [for which we used the results in Ref. [59] multiplied by θ 2 of Eq. (2)]. We only consider 2-body decays in our analysis.
As we shall see, relaxion dark matter is in the keV range. In this mass ballpark the relaxion can only decay into photons and neutrinos. The decay into photons proceeds through the mixing with the Higgs and through the loop-induced coupling of Eq. (4). For simplicity, we assume that neutrinos are Majorana fermions, in which case the decay in this channel is suppressed compared to the one into photons as it proceeds via higher dimensional operators (see, e.g., Ref. [53]). If neutrinos are Dirac fermions, this can be the dominant decay channel. Nevertheless, the bounds from indirect detection on the DM decaying into photons (see next section) imply stronger constraints on the relaxion lifetime.

V. RESULTS AND DISCUSSION
We performed a scan looking for points in {Λ, g , Λ b , f } which can satisfy the DM hypothesis. For each point, the value of H I is fixed by plugging Eq. (21) into Eq. (48), and then requiring to match the observed DM abundance. This value has to be compared with the range allowed from the conditions on particle production. The result is shown in Fig. 1. The green region is the one where the relaxion is stable, all the bounds on a successful relaxation with particle production are simultaneously satisfied and the relaxion abundance matches the observed DM one (for a given range of Λ b and f ). The light green part is the one in which, additionally, the constraints from indirect detection are satisfied. In the yellow region the relaxion can be stable, but it is overproduced. Finally, in the orange region the relaxion's lifetime is shorter than the age of the universe. Table I shows the allowed parameter space for three benchmarks.
On top of the above constraints, we applied a lower bound on the DM mass from structure formation. The free-streaming length λ NR is constrained by observations from Lyman-α forest [60,61], which are in tension with a thermal relic below a few keV. We estimate the free-streaming length as [55]: where a(t) is the scale factor as a function of time within a 0 being the scale factor today, t FI refers to the time when DM abundance freezes-in, after which it freestreams relativistically until t NR t FI when it becomes non-relativistic. The relaxion then free-streams nonrelativistically up to the matter-radiation equality at t EQ . The time t NR can be easily obtained if the DM velocity distribution is thermal. However, in our scenario the DM distribution function can depart from a thermal distribution as the relaxion is produced out-of-equilibrium, which may weaken these bounds (see Refs. [62,63] and references therein). It would be important to further explore such feature which we leave for future work. Here we simply impose that the relaxion should be heavier than 2 keV. Figure 2 shows how the allowed region depends on the value of the reheating temperature. A drawback of our scenario is that the reheating temperature is rather low, T 0 30 MeV, 200 MeV, 30 GeV for g/g = 1, 10 3 , 10 6 , to avoid overabundance. Ultimately, this is due to the UV sensitivity of the production mechanism. If thermal equilibrium is reached, in this mass range the relaxion would be overabundant by a factor of 10 − 100. Oppositely, the correct relic abundance could only be obtained for a correspondingly lighter particle, which would then be too light to comply with the warm DM mass lower bound.
Measurements of the abundance of light elements, large scale structure data, and anisotropies of cosmic microwave background temperature constrain late-time entropy production, which then restricts the reheating temperature to be larger than 1 − 4 MeV [64][65][66][67][68][69]. Below this temperature, the universe behaves like radiation, and only very small entropy injections are possible.
A low reheating temperature can be achieved even if the temperature of the SM plasma at the end of inflation rises to much higher values, during a phase of entropy injection [70][71][72]. This is indeed expected to happen when the inflaton decays perturbatively into SM particles. The temperature first rises to a maximal value T max ∼ T 1/2 0 (H I M Pl ) 1/4 , then it decreases with the typical dependence on the scale factor T ∝ a −3/8 . This behaviour proceeds until the decay of the inflaton ceases at a time of order the inverse decay width of the inflaton, at which radiation dominance begins with the standard T ∝ a −1 behaviour. During the reheating phase, entropy is continuously created, and the Hubble rate scales as H ∼ [g * (T )/g 1/2 * (T 0 )]T 4 /(T 2 0 M Pl ) (expansion is faster for lower reheating temperature). In such a scenario, the abundance of relic particles are altered compared to the standard radiation dominance calculation. On the one hand, particles with mass larger than the reheating temperature can be copiously produced [73]. On the other hand, which is the case relevant here, particles with a freeze-out temperature larger than T 0 are diluted by entropy injection, and their abundance is smaller than in the standard freeze-out computation. As an example, in Refs. [72,74] was argued that SM model neutrinos in the keV mass range (hence now excluded) could have the right relic abundance to be a warm DM candidate. The relic abundance of long-lived particles in a low reheating scenario is studied in the literature for many kind of DM candidates, such as sterile neutrinos [75][76][77], supersymmetric particles and more generic WIMPs [72,[78][79][80][81][82], heavy particles [73], axions [72,[83][84][85][86][87][88][89]. Here we just assume that at T 0 the relaxion abundance is negligible, and that its relic abundance is built up during radiation dominance.
Finally, let us mention that baryogenesis mechanisms that require a large temperature are also viable in this scenario. As an example, electroweak baryogenesis is viable in this case for T 0 as small as 1 GeV [72,90], thus favoring large values of the ratio g/g . While this is an interesting option, a concrete realization which connects the relaxion to DM and to baryogenesis is beyond the scope of the present work.
Strong constraints on the model come from the observations of the galactic and extra-galactic diffuse X-ray and γ-ray background. We consider the constraints on decaying DM from Ref. [96] which uses the diffuse photon spectra data from different satellites. For our parameter space, which comprises masses around the keV range, the relevant bounds are given by the satellites HEAO-1 [97] and INTEGRAL [98]. In Fig. 1 we show in light green the region in agreement with the bounds on the lifetime of a scalar DM decaying into two photons, τ φ 10 26−28 s for m φ > 5 keV [96]. This constrains the relaxion mass to be m φ {5 keV, 5 keV, 17 keV}, respectively for the three benchmarks in Tab. I. Extrapolating the bound from Ref. [96] to lower masses further constrains the parameter space, but the results are qualitatively similar. This places relaxion DM in a knife-edge position: on the one hand, new results from indirect searches in the keV mass range could rule out this scenario; on the other hand, a numerical solution of the Boltzmann equation could weaken the lower bound on the relaxion mass, thus opening the parameter space for lighter DM.
FIG. 2. Allowed dark matter region as a function of the reheating temperature. The region shrinks for higher temperature T0.
such bounds. On the other hand, the parameter space for successful relaxation with particle production is also subject to some variation. The requirements for a successful particle production mechanism could be relaxed by considering the relaxion velocity. This has two natural reference values: the first is the slow-roll velocitẏ φ roll ∼ V φ /(3H I ), while the second is the minimal velocity to overcome the wigglesφ b ∼ Λ 2 b . The two are related byφ roll >φ b . In deriving the relations, we always chose the value that lead to the most conservative bound. Some of the conditions on particle production could therefore be weakened by choosing a different value for the velocity, but we do not pursue this possibility further.

VI. CONCLUSIONS
In this work, we showed that the relaxion mechanism can naturally provide a phenomenologically viable warm DM candidate in the keV mass range. We identified the relevant parameter space in the scenario in which relaxation happens during inflation, using particle production as a source of friction. We discussed astrophysical and indirect detection constraints on the model.
Recently, there has been an increasing interest in DM direct detection experiments that can probe the sub-MeV mass range (see, e.g., Refs. [100,101]). The relaxion would be a well motivated DM candidate in the keV range, which encourages new studies in this mass ballpark.
It would be interesting to further explore the consequences of such a model on structure formation, and perform a dedicated analysis of the indirect detection bounds. We leave these studies for future work.

ACKNOWLEDGMENTS
We are happy to thank Géraldine Servant for encouragement, discussions and helpful comments. We are also grateful to Filippo Sala for valuable comments. NF thanks the organizers and participants of the CERN TH Institute on Physics at the LHC and Beyond for interesting discussions while part of this work was completed.