A Time-Varying Fine Structure Constant from Naturally Ultralight Dark Matter

We present a class of models in which the coupling of the photon to an ultralight scalar field that has a time-dependent vacuum expectation value causes the fine structure constant to oscillate in time. The scalar field is assumed to constitute all or part of the observed dark matter. Its mass is protected against radiative corrections by a discrete $\mathbb{Z}_N$ exchange symmetry that relates the Standard Model to several copies to itself. The abundance of dark matter is set by the misalignment mechanism. We show that the oscillations in the fine structure constant are large enough to be observed in current and near-future experiments.


I. INTRODUCTION
With the discovery of the Higgs boson [1,2], the Standard Model (SM) of particle physics is now complete. However, cosmological observations tell us that visible matter constitutes only about 20% of the matter in the universe, the remainder being composed of some form of non-luminous dark matter [3]. Despite contributing five times more to the energy budget of the universe than visible matter, the basic nature of dark matter has thus far eluded us.
Moduli are among the most well-motivated and compelling dark matter candidates [4][5][6][7][8]. They often arise in ultravioletcomplete theories such as string theory, where their vacuum expectation values play a role in determining the values of fundamental parameters such as the fine structure constant. The properties of moduli make them ideal candidates to play the role of dark matter. During inflation, quantum fluctuations are stretched to super horizon length scales, with the result that all scalars lighter than the Hubble scale acquire large random vacuum expectation values. After inflation ends, a light scalar such as a modulus remains frozen away from its minimum by Hubble friction. At late times, the modulus begins to oscillate about its minimum, contributing to the energy density of the universe as a component of dark matter. This framework, termed the misalignment mechanism [9], is a natural way of generating an abundance of extremely cold bosonic dark matter from light fields such as moduli.
There has recently been a surge of interest in finding new ways to search for very light dark matter, see e.g. . One of the the most exciting approaches to finding modulus dark matter is applicable when it is so light that its Compton wavelength is macroscopic, so that the period of its oscillations, set by the inverse of its mass, can be directly observed. As moduli set the value of fundamental constants, they can be searched for by looking for time dependence ("chronovariance") of these parameters. As a wave of modulus dark matter washes over us, it oscillates with a character-istic cos ωt ≈ cos mt time dependence, where m is the mass of the modulus field. Therefore a natural way to search for ultralight modulus dark matter is to look for time dependence of the fundamental constants. For early work on time variation of fundamental parameters, see for example [33][34][35][36][37].
For dark matter to be observable in this way, it must be strongly coupled enough to be detected over noise while oscillating at a low enough frequency that its effects are not averaged away. The condition that the frequency be low translates into the requirement that the mass be small. Since the mass of the modulus receives radiative corrections that depend on the size of its couplings, there is some tension between the condition that the mass remain small and the requirement that the couplings are large enough to be observable, resulting in a naturalness problem.
To understand the naturalness problem, we parametrize the couplings of the modulus φ to the electromagnetic field in the form, where κ is defined as √ 4πG N , where G N is Newton's constant and d e is a dimensionless coupling constant that parametrizes the strength of the interaction. For a given mass of the modulus φ there is a bound on the amplitude of its oscillations, and therefore on the magnitude of the variation of α, from the condition that the modulus contribute no more to the energy budget of the universe than the observed dark matter contribution. Current experiments are most sensitive to a modulus mass in the range 10 −22 − 10 −5 eV. The corresponding bound on d e stands at d e ≲ 10 −8 − 10 −1 , with the exact value depending on the modulus mass. This coupling in Eq. (1) gives rise to a radiative contribution to the mass of the modulus at two loops 1 , where Λ represents an ultraviolet cutoff. As an example, consider a modulus mass of 10 −18 eV for which the bound stands at d e ≲ 10 −3 . For this value of d e there is a contribution to its mass coming from loops involving the top quark of order 10 −11 eV. Thus, there must be a cancellation in the mass squared of this scalar between the top quark contribution and the bare mass to at least one part in 10 14 . More generally, for any value of the modulus mass there is an upper bound on d e above which the conflict comes to the fore. Furthermore, any understanding of the abundance of the modulus requires a consistent ultraviolet completion, since the behavior of a modulus in the early universe is very complicated and is extremely sensitive to the presence of other fields.
In this paper, we present a framework in which the modulus that controls the fine structure constant can naturally remain ultralight while constituting all of the observed dark matter in the universe. In order to make the modulus naturally light, we employ the mechanism proposed in Ref. [38]. Accordingly, we introduce N copies of the SM, where N is a number of order a few. If the same modulus controls the fundamental parameters of all the N copies of the SM while nonlinearly realizing the Z N symmetry, then the sum of the contributions to the potential of the modulus from each of the copies of the SM cancels to a very high degree of accuracy. The end result is then a modulus with a parametrically smaller mass than naively expected.
A very attractive feature of this framework is that, for certain ranges of the modulus mass and couplings, the misalignment mechanism naturally allows the modulus to constitute all of the observed dark matter. For a given mass and couplings of the modulus, we can determine the thermal corrections to its potential in the early universe. The simplest region of parameter space is that in which the Hubble friction holds the modulus in place and oscillations only begin after thermal effects have become subdominant to the zero temperature potential. Even in this most simple of scenarios, there exists a vast region of parameter space in which the modulus can constitute all of dark matter. In general, finite temperature effects can drastically alter the behavior of the modulus at early times. They can act either to reduce or increase the abundance of dark matter by relaxing the modulus to the minimum or maximum of its potential at high temperatures. Including these effects expands the region of parameter space in which the modulus can play the role of dark matter.
Prospects for probing the time variation of fundamental constants due to ultralight scalars are very promising [15,24,39,40]. In fact, existing data from experiments [41][42][43][44][45][46][47][48] is already sufficient to place constraints on the scenario we propose in the mass range 10 −22 eV < m φ < 10 −5 eV. These limits can be significantly improved in the future by increasing the integration times in the optical-optical clocks and ultimately by taking advantage of the anticipated 229 Th nuclear-optical clock [39,49]. Further advancements can be achieved with the recently proposed earth and space based atomic gravitational wave detectors, which will rely on atom interferometry [40,[50][51][52]. The MAGIS and AION experiments [50,51] will be earth based interferometers that are planned to gradually increase in size to finally reach a length of 1 km. MAGIS is currently building a 100 m interferometer [50], which should be able to probe the proposed model in the 10 −16 eV < m φ < 10 −14 eV mass range. The 1 km stage should increase this range to 10 −16 eV ≲ m φ ≲ 10 −12 eV. Once built, the AION experiment will be able to improve the bounds set by MAGIS at both the 100 m and 1 km scale for the same range of masses m φ . The AEDGE experiment [52] is a proposed continuation of the AION experiment that will take advantage of satellites in order to increase the scale of the detector to thousands of kilometers. While it is a very distant prospect, it will have greatly improved sensitivity for masses in the range, 10 −19 eV < m φ < 10 −13 eV.
The outline of this paper is as follows. In Sec. II, we discuss the framework and explain the mechanism that protects the mass of the ultralight modulus. In Sec. III, we study the cosmological history of this class of models and show that the misalignment mechanism can allow the modulus to constitute all of the observed cold dark matter. In Sec. IV we present our results. We determine the current bounds on this class of models and outline the region of parameter space that will be explored by current and future experiments. We conclude in Sec. V

II. THE FRAMEWORK
In this section we construct a class of models in which the fine structure constant oscillates in time, and which are free of naturalness problems. We consider a complex scalar Φ that is charged under an approximate global U(1) symmetry. This field is assumed to acquire a vacuum expectation value, spontaneously breaking the global symmetry at a scale denoted by f . It is convenient to employ an exponential parametrization of Φ, Here ρ represents the radial mode in the potential for the scalar field after symmetry breaking, while ϕ denotes the pseudo-Nambu-Goldstone boson. We write ϕ = ϕ 0 + φ, where ϕ 0 is the vacuum expectation value (VEV) of ϕ while φ represents the fluctuation about this expectation value. The field φ is assumed to couple to the electromagnetic field strength as shown in Eq. (1). Although this interaction is nonrenormalizable, it can be generated by coupling Φ to some heavy charged fermions that are integrated out. As a consequence of this coupling, changes in the background value of ϕ will cause variations in the fine structure constant. However, this coupling violates the U(1) global symmetry explicitly. We therefore expect that, in general, it will generate a potential for ϕ, leading to a quadratically divergent mass. Since the parameter d e controls both the amplitude of the modulation as well as the magnitude of the potential, we require some mechanism that protects the potential against large radiative corrections while still admitting an observable signal. As we now explain, this can be done by employing a discrete Z N symmetry under which the pseudo-Goldstone ϕ transforms nonlinearly [38]. We introduce N copies of the SM, each with its own matter content and gauge groups. We label each of these different copies by an index i, where i runs from 1 to N . The N copies of the SM are related by a discrete Z N symmetry under which i → (i + 1). To each copy of the SM we add a heavy Dirac fermion Ψ i that carries unit charge under the corresponding U (1) Y hypercharge gauge symmetry, but not under any of the other gauge groups. Each of the N fermions Ψ i is assumed to have a Yukawa coupling to the scalar Φ. Under the Z N symmetry these fields transform as, The discrete symmetry forces the Yukawa couplings, gauge couplings and masses to be same across all the N sectors. The Lagrangian for the Ψ i is restricted to have the form, Here the Yukawa coupling y can be chosen to be real without loss of generality as any phase of complex y can be absorbed into the definition of the complex scalar Φ.
After symmetry breaking, the relevant part of the Lagrangian becomes where the effective mass of each fermion depends on the VEV of the pseudo-Nambu Goldstone boson ϕ Here the parameter is defined as ≡ √ 2yf M . After integrating out the fermions Ψ k we get a contribution at one loop to the effective potential for the pseudo-Nambu Goldstone boson, where Λ is an ultraviolet cutoff. This integral can be evaluated exactly.
After dropping terms that do not depend on ϕ and terms that vanish in the limit M Λ → 0, we obtain where j ≡ cos(ϕ f + 2πj N ). The first two terms in Eq. (9) are quadratically divergent and logarithmically divergent. However, they can be seen to independent of ϕ as long as N > 4 by virtue of the following identity [38] N j=1 where C is a constant. The finite contribution is, however, ϕdependent. It can be extracted by expanding the last term in Eq. (9) in the small parameter . To all orders in , we get contributions from each sector that only differ by the phase of the cosine. Due to the identity in Eq. (10), all the terms up to O( N −1 ) cancel leaving the leading ϕ−dependence to appear at O( N ). These cancellations are a direct consequence of the unbroken Z N symmetry. After performing the integral in (8), the leading ϕ-dependent piece arises at order N where F (N ) is given by This potential has N minima at the locations ϕ m f = (2m − 1)π N where m runs from 1 to N . Integrating out the fermions Ψ i also leads to an effective coupling of ϕ to the gauge kinetic terms of each of the N hypercharge gauge fields at low energies, One of these N copies corresponds to the hypercharge gauge boson of the SM. Since the effective potential Eq. (11) is invariant under ϕ → ϕ + 2πf N , we have freedom to shift ϕ so that the coupling to SM hypercharge gauge boson simplifies to Since in our normalization B µ = A µ −W 3 µ , the coupling to SM photons after spontaneous symmetry breaking takes the same form, When ϕ performs small oscillations about one of the minima ϕ m of (11), it gives rise to a linear term in the coupling to the photon, where ϕ m is the m-th minimum of the potential, Eq. (11), and ϕ = ϕ m + φ. Comparing to Eq. (1), we see that in terms of the parameters of our model, d e is given by, We are now at a stage where we can demonstrate the improvement in naturalness that our model affords. As discussed in the Introduction, the natural expectation is that the mass of φ scales as which is just Eq. (2) with the divergence cut off by the mass of the new charged fermion, in this case M . To see how our model compares to this, we expand Eq. (11) about the minimum to obtain The term in brackets represents the improvement with respect to the naive estimate. In order for this expression to be valid, we had assumed that N > 4. Because this correction factor is proportional to N −2 , we see that as long as is small the mass is parametrically smaller than expected from naive considerations. This demonstrates that our construction can indeed solve the naturalness problem discussed in the Introduction. We note that the mass term for the modulus is generated at one-loop order rather than being given by the naive 2-loop estimate in Eq. (2). This is because both the mass of the modulus and its coupling to the hypercharge gauge boson arise at the same one-loop order when the heavy fermions Ψ i are integrated out, rather than the mass being generated radiatively from the coupling.
At this stage, the new fermions Ψ k are electrically charged stable fermions. If the reheat temperature is larger than their mass, then these particles obtain a thermal abundance and will tend to overclose the universe if their masses lie above the TeV scale. In order to avoid this, we allow each of the N Ψ k to decay by introducing a small mixing with the right-handed τ lepton of the corresponding sector through the interaction For M = 10 TeV, we require m ≳ 100 eV in order to have the Ψ particles decay prior to Big Bang nucleosynthesis (BBN) and not pose a cosmological problem.

III. COSMOLOGICAL HISTORY
In this section, we will consider the cosmological history of this class of models. We work under the assumption that of the N sectors, ours is the only one that is reheated after inflation 2 . We further assume that ϕ is homogenized as a result of inflation. The scalar field evolution is governed by the equationφ where the potential is given by Here V 0 is the zero temperature potential given in Eq. (11) while V T (ϕ, T ) is the contribution to the potential from finite temperature effects. An explicit expression for V T may be found in Appendix A. At temperatures T ≫ M , the finite temperature contribution to the potential can be approximated as, Once the temperature falls below its mass, the vector-like fermion Ψ begins to exit the bath. Accordingly the form of the thermal contribution to the potential undergoes a change.
In the temperature range M ≳ T ≳ 100 GeV it is well approx-imated by the expression, (24) For T ≲ M 20, the exponential suppression of the second term in the bracket means that it can be neglected.
The evolution of the scalar field ϕ at early times is governed by the extent to which its oscillations are damped by Hubble friction. To keep track of whether the system is underdamped or overdamped at temperature T , we define the parameter where m(T ) represents the contribution to the mass of the modulus from finite temperature effects, For η(T ) < 1 the system is overdamped at temperature T , while for η(T ) > 1 it is underdamped. At early times η(T ) increases as the universe cools down. This can easily be seen from Eq. (23) by noting that when T ≫ M , the mass scales linearly with temperature, m(T ) ∝ T , while the Hubble parameter decreases faster, H ∝ T 2 . Eventually, at temperatures T of order M , the growth of η(T ) slows down, reaching its maximal value η p at a temperature T p = 2M 5 . The value of η p is given by Here M pl is the reduced Planck mass and g * = 106.75 is the effective number of degrees of freedom at T ≲ M , consisting of just the SM fields. For η p ≳ 1, the field oscillates about the minimum of the thermal potential at temperatures T of order M . The oscillations begin at a temperature T osc such that m(T osc ) = H(T osc ), given by Here g Ψ * = 110.25 is the effective number of degrees of freedom at temperatures T ≫ M , which consists of the SM fields and the associated fermion Ψ.
Below the temperature T p , η(T ) decreases rapidly as the exponential suppression of the second term in Eq. (24) takes effect. Eventually at temperatures T ≲ M 20, the contribution from the first term in Eq. (24) becomes dominant. In this temperature regime, the mass of the modulus and the Hubble parameter scale identically, m(T ) ∝ H(T ) ∝ T 2 . As a result η(T ) reaches a terminal value η, given by The parameters η p and η play an important role in the description of the system. The value of η p indicates whether the system undergoes oscillations when T ≳ M , while η determines the behavior of the system at temperatures T ≲ M 20. These parameters are related as From this relation we can see that in the regime which is always overdamped at early times, η ≪ 10 −4 , the modulus field is effectively frozen until the temperature-independent part of the potential becomes dominant. In contrast, for η > 1 the oscillations of the field begin at the temperature T osc ≳ M and continue till the present time. In the intermediate regime 10 −4 ≲ η ≲ 1, the field oscillates at temperatures T of order M . However, these oscillations have ceased by the time the temperature falls below M 20 and only resume once the temperature-independent contribution to the potential begins to dominate.
In what follows below we obtain analytic expressions for the contribution of the modulus φ to the energy density of the universe. We focus on the limiting cases of η ≪ 10 −4 and η ≳ 1, deferring the intermediate range of η to our numerical study.
For η ≪ 10 −4 the early time evolution is especially simple as Hubble friction freezes the field in place at high temperatures, so that its evolution is independent of V T . It follows that for sufficiently small η the field is effectively fixed until the mass approaches its zero temperature value, m = m φ . Oscillations begin when 3H ∼ m φ , corresponding to a temperature At this point, we make the standard approximation that ϕ transitions instantly from an overdamped harmonic oscillator to an underdamped one. At this stage the potential Eq. (22) is dominated by the zero temperature term, The mass m φ is constant so the field will oscillate as an underdamped harmonic oscillator around a minimum ϕ m , with the amplitude decaying as Since the field is frozen, the value of ϕ when 3H ∼ m φ can lie anywhere in the range ϕ s ∈ (0, 2πf ). Therefore ϕ can fall into any of the N minima ϕ m of the zero temperature potential, Eq. (32), with equal probability. Consequently the initial misalignment value of the excitation around that minimum ϕ m can lie anywhere in the range φ s ∈ (−πf N, πf N ). Due to this randomness in the initial condition and the homogeneity of φ, it is not possible to determine the exact value of the field today. However, what can be done instead is to average over all possible initial misalignment values, φ s ∈ (−πf N, πf N ), given that all are equally likely. The expectation value of the energy densityρ φ can be related to an effective initial amplitude of the field φ s,eff , where φ s,eff = ⟨φ 2 s ⟩. The value of ⟨φ 2 s ⟩ can be obtained by averaging over all values of the initial misalignment, leading to φ s,eff = πf √ 3N . As a result the expected amplitude of the field today is given by where T 0 represents the current temperature of the Universe. This leads to a final result for the energy density in φ today, In this subsection, we consider the range of parameter space in which the behaviour of ϕ in the early universe at temperatures T ≲ M 20 is like that of a harmonic oscillator that is either underdamped or close to critically damped. When discussing finite temperature effects, it is necessary to distinguish between the cases when N is even and N is odd. The reason for this can be seen from Fig. 2. In the case of even N , the thermal potential pushes ϕ towards a maximum of the zero temperature potential. At late times, this results in an increase in the amplitude of oscillations, leading to enhanced abundance of dark matter and a larger signal. For odd N , the thermal potential pushes ϕ towards a minimum of the zero temperature potential, decreasing the amplitude of oscillations at late times resulting in a suppressed abundance. In this case the situation is even worse because the minimum that ϕ is pushed towards is the one in which there is no linear coupling to the photon, so that the signal under consideration is greatly suppressed. Rather than a linear coupling, this minimum leads to a quadratic coupling, which gives rise to different phenomenology [53,54] that we do not consider here.
a. N odd In the early universe ϕ is driven towards the minimum of the finite temperature potential, ϕ m = πf . From Fig. 2 we see that the minimum of the finite temperature potential is also a minimum of the zero temperature potential. Consequently the range of values of ϕ at the time when oscillations begin, H = 3m φ , is significantly smaller than initial (0, 2πf ).
It follows from this that for each m φ , there is a minimal value of η above which this new range covers only the central minimum of the zero temperature potential, ϕ s ∈ ((π − π N )f, (π + π N )f ). Above this critical value of η, which we denote by η 0 , ϕ always ends up in the central minimum. Furthermore, from Eq. (16) we see that at this central minimum, ϕ m = πf , the linear coupling of φ to photons vanishes, resulting in a greatly suppressed signal. The value of η 0 can be well approximated by considering evolution only in the regime T ≲ M 20. The evolution at higher temperatures is less efficient at focusing the range of values of ϕ due to the sudden drop of the modulus mass m(T ) at T ∼ M , which has the effect of regenerating the amplitude of the field. The subsequent evolution of ϕ in the regime when T ≲ M 20 and η < 1 scales with temperature as It follows that η 0 satisfies the condition, From this we obtain an expression for η 0 , We see that in the case of odd N , for η ≥ η 0 , the field is pushed to the minimum in which φ does not have the coupling shown in Eq. (1), resulting in a greatly suppressed signal. The abundance of φ is also suppressed limiting its contribution to dark matter.
b. N even For even values of N the picture is completely different. The minimum of the finite temperature contribution to the potential, which is at ϕ m = πf , now coincides with the maximum of the zero temperature potential, as can be seen in Fig. 2. Therefore, in the regime η ≳ 1 the initial value of the modulus ϕ at the time when oscillations begin will be close to π and the initial amplitude will be very close to This is very different from the case when the system is overdamped.
The amplitude of the oscillations of φ today depends on how close the field is to the minimum of the finite temperature potential at the time that the oscillations begin. Given that the mass of the scalar changes as a function of temperature, the simplest way to determine the behavior of ϕ(T ) at these early times is to employ the conservation of the number density of ϕ. The comoving number density of ϕ is approximately conserved in the limit that ϕ is close enough to a minimum that its potential is quadratic and as long as the WKB approximation holds, dm dt ≪ m 2 . The conserved comoving number density is given by Based on these considerations the effective initial value of the modulus φ s,ef f when the zero temperature potential begins to dominate can be obtained as (41) We continue to track the evolution of the field after it rolls down to the new minimum, which happens at temperatures T of order T s .
In contrast to the case when η ≪ 10 −4 , the position of ϕ when the zero temperature contribution to the potential begins to dominate is very close to ϕ = πf . Since this point corresponds to an extremum of both the finite temperature and zero temperature contributions to the potential, the gradient of the potential at ϕ = πf vanishes independent of the temperature. The fact that the potential in the neighborhood of this point is very flat leads to a delayed onset of oscillations [55,56] and interesting phenomenological signatures [57,58]. As shown in Appendix B, for the parameters considered in this paper this delay is long enough to ensure that when the oscillations about the new minimum begin, the contributions to the potential from finite temperature effects are already negligible. Therefore, to a good approximation we can assume that the potential has the zero temperature form in Eq. (11) from the time that the oscillations begin.
It is tempting to assume that oscillations about the true minimum occur in a harmonic potential, starting from the temperature T s , defined in Eq. (31), with an initial amplitude φ s,eff = πf N and continuing till today. This leads to the following expression for the contribution of the oscillating field to the energy density today, However, this expression does not give the correct result for the energy density as it fails to account for the the delay in the time at which the oscillations begin. To include this effect we employ the following empirical approximation [56], and β(y) = ln y + ln 2 1 4 √ 3 π 1 2 Γ(5 4)N .
The coefficient C(y) defined in Eq. (44) corresponds to the correction arising from the delayed onset of oscillations. The difference between the energy density in Eq. (43) and the naive estimate in Eq. (42) turns out to be quite significant. For M = 10 TeV the correction factor isρ u φ ρ naive φ ≈ 9 for N = 6 andρ u φ ρ naive φ ≈ 10 for N = 10. Since the field is pushed very close to the central maximum of V 0 (ϕ) potential, one could imagine that instead of the field homogeneously rolling down to a unique minimum, quantum fluctuations would push the field in some regions of space to the minimum on the other side of the hill, resulting in domain walls. However, just after inflation, quantum fluctuations are given by [55], The ratio of fluctuations to the average value of the field is constrained by CMB observations, Because this quantity changes at most logarithmically between the end of inflation and now [59] and the emergence of domain walls requires δϕ s ∼ ϕ s,ef f − πf , this scenario does not take place in our model.

IV. SIGNAL
In this section, we determine the region of parameter space populated by our model and explore the implications for experiment. Following the conventions employed in [10,11,39], we parametrize the changes to α in terms of the variable d e defined in Eq. (1), Here ∆α represents the amplitude of oscillations of the fine structure constant while φ local denotes the amplitude of oscillations of the scalar field at the location of the earth. We introduce a parameter r that represents the fractional contribution of φ to the total energy density in dark matter, Experiments place a bound on ∆α α for a given frequency of oscillation m φ . Then, for a given dark matter fraction r, this can be translated into a bound on d e using Eq. (48). As can be seen from Eq. (17), the value of d e is a function of the parameters , f , M and N . From this equation we further see that the value of d e also depends on the minimum ϕ m that the field settles in, which is in turn determined by the initial conditions. In our analysis, we take sin (ϕ m f ) to be the average of available values, e.g. for N = 6 we have ⟨ sin (ϕ m f ) ⟩ = 2 3 in the regime η ≪ 10 −4 and ⟨ sin (ϕ m f ) ⟩ = 1 2 in the regime η > 1. Additionally for odd values of N we ignore the contribution of the central minimum (sin (ϕ m f ) = 0) to the average. This allows us to determine the expected value of d e consistent with a given set of parameters.
In order to make contact with experiment it is convenient to eliminate the parameter in favor of m φ . The expression for m φ in terms of the other four variables is given in Eq. (32). The parameter space of our model can then be described in terms of the four variables f , m φ , M and N . Using the results of Section III, the expectation value of the dark matter abundance can also be expressed in terms of these parameters, after averaging over the initial conditions for the scalar field ϕ. The requirement that, after this averaging, the field φ constitutes all of the observed dark matter places a restriction on the allowed parameter space and allows us to fix the value of f in terms of the three other variables. Then d e is determined in terms of m φ , M and N . Taking advantage of the results of Section III, we can obtain semi-analytical expressions for d e as a function of these three remaining parameters in various regimes. These will prove helpful in illuminating the main features of the detailed numerical results, which we will present later.  A. Analytic Results for η ≪ 10 −4 In the overdamped limit, we can use Eqs. (19), (35) and (17) to obtain an expression for d e in terms of m φ , M and N , Here and For convenience, a few values of A 1 (N ) are given in Table I. This analytic approximation reproduces our detailed numerical results up to an accuracy of around 10%. From Eq. (50) we can see that d e decreases as we raise the mass M of the fermions Ψ. Then, by setting this mass to the lowest value allowed by experiment, we can place an upper bound on d e as a function of m φ for any given value of N .
B. Analytic Results for η ≳ 1 As η increases, the scalar field eventually enters the regime η ≳ 1. In this regime, the case of odd N is very different from that of even N . For odd values of N , the region where η ≳ 1 does not give rise to an observable signal because the field is trapped in the wrong vacuum as discussed in Section III. The boundary of this region is marked by η 0 , defined in Eq. (38). The condition η < η 0 translates into a restriction on the allowed range of m φ , M and N . From Eqs. (19), (29) and (17) we can determine the value of d e at the boundary of this region η = η 0 to be  where The numerical values of A 2 (N ) for a few sample points were given in Table I. The allowed parameter space is restricted to values of d e lower than Eq. (52). The line obtained from Eq. (52) is found to reproduce our numerical results up to a factor of 2.
We now turn our attention to the case of even N . The expressions for the energy density in the case η ≳ 1, Eq. (43), and the case η ≲ 10 −4 , Eq. (35), differ only by their proportionality constant. It follows from this that d e scales with m φ and M in the same way in both regimes, Here m 1 , obtained from Eqs. (35), (19) and (29), corresponds to the mass of the modulus for which η = 1 and ρ φ = ρ DM .
The numerical values of A 3 (N ) are given in Table I for a few reference points. Equation (54) reproduces the detailed numerical solution up to an accuracy of about 20%.

C. Numerical Results
The analytic solutions found in the subsection above are valid in the limiting cases when the scalar field is either highly overdamped or highly underdamped. In order to determine the solution in the region of parameter space 10 −4 ≲ η ≲ 1 where the system transitions between these two regimes, we find it necessary to solve Eq. (21) numerically. We parametrize the model in terms of the four parameters f, m φ , N and M . As explained earlier, lighter fermion masses are associated with larger values of d e . In our study we therefore consider two different values of the fermion mass, M = 10 TeV and M = 1 TeV, which are close to the current lower bound from collider experiments. In addition, we consider four different values of N , the odd values N = 5 and 11 and the even values N = 6 and 10. We then scan over f for different values of m φ . For each (f, m φ ) pair, the field was evolved from T i = f to T = T 0 starting from 1000 random initial conditions. In the case of odd N , we discard any (f, m φ ) pair such that, for more than 90% of initial conditions, the theory ends up in the vacuum with vanishing signal. For each point we obtain the contribution of φ to the dark matter density. We also determine the value of ⟨ sin (ϕ m f ) ⟩, which is obtained by averaging over solutions with sin (ϕ m f ) ≠ 0, and use this to find the value of d e at each point from Eqs. (17) and (19).
The results of our numerical study are shown in Figs. 3 and 4, where we have plotted d e as a function of m φ for these theories, along with the current limits and the projected reach of future experiments. Our goal is to identify the region of parameter space that is naturally populated by these models. To this end, for each m φ we have singled out the value of d e such that, for any d e larger than this, more than 90% of points will lead to less than the observed abundance of dark matter, ρ φ < ρ DM . Separately, for each m φ we have singled out the value of d e such that, for any d e smaller than this, more than 90% of points will lead to more than the observed abundance of dark matter, ρ φ > ρ DM . These values have been plotted in Figs. 3 and 4 as the two light blue (M = 1 TeV) lines. The region between these lines, which corresponds to the natural parameter space for M = 1 TeV, has been shaded in. The natural paramater space for M = 10 TeV has also been shown, shaded in dark blue. We explicitly show the values of the parameters for a few reference points in Table II. In order to simplify the numerical computation, the number of effective degrees of freedom in the thermal bath has been held fixed at g * = 106.75. The coefficient of the thermal contribution to the potential from the first term in Eq. (24), which also decreases as heavier species go out of the bath, has also been held fixed. These simplifications have an effect on the temperature at which the final oscillations begin, resulting in an overestimate of d e in the case of the lowest values of m φ by a factor close to two. The error is smaller for larger values of m φ since the oscillations begin earlier, and so there are a greater number of degrees of freedom in the bath at the corre-

V. CONCLUSIONS
Light moduli are one of the most attractive dark matter candidates. Not only are they ubiquitous in ultraviolet-complete models such as string theory, but the misalignment mechanism provides a simple explanation for their abundance. Despite these appealing features, detecting modulus dark matter remains a challenge. The problem is that the large coupling required to find them tends to be at odds with their very light mass, resulting in a hierarchy problem. In this article, we have shown that a nonlinearly realized Z N symmetry can naturally protect the mass of a modulus against radiative corrections, provided the modulus itself transforms nonlinearly under the Z N . It remains an intriguing open question whether a model that exhibits these features can be constructed within a string theory framework.
Much like the QCD axion, these models have a region of parameter space in which they naturally reproduce the observed dark matter abundance via the misalignment mechanism. The regions of parameter space populated by these models were shown in Figs. 3 and 4 for a few values of N . We see from these figures that future experiments searching for chronovariance of the fine structure constant are expected to probe a large part of the preferred parameter space. This class of models is both simple and testable and provides additional theory motivation for exciting new quantum limited experiments. (A5) The leading finite temperature contributions to the free energy of the universe that involve the hypercharge gauge coupling arise at two loops. The problem therefore reduces to computing these two loop diagrams. There are three types of diagrams labelled by a, b and c that contribute at this order, where The diagrams that give rise to F i a , F b and F c are shown in Fig. 6. Here the index i runs over all quark and lepton flavors. Using the methods of thermal field theory [60,61] the first diagram, which represents the correction to the free energy from diagrams involving the SM fermions, can be evaluated as where we follow the conventions of [60]. In these conventions theγ µ represent the Euclidean Dirac matrices, The momenta in the loop integrals are in Euclidean space and have components P = (ω i n , ⃗ p), where i = b, f . The symbol ⨋ refers to integration over the spatial momenta along with summation over the Matsubara frequencies, which are defined as for fermions and bosons respectively. If the momenta below the symbol ⨋ appear between curly brackets {...} it corresponds to summation and integration over fermionic boundary conditions while the absence of the brackets indicate bosonic boundary conditions. After the usual Dirac algebra we can ex-press Eq. (A8) in terms of known standard integrals [62,63], In the T m → ∞ limit the last integral vanishes [63]. The remaining integrals are well known, Putting these together we arrive at the result in the high temperature/low mass limit, where we have expressed g ′ in terms of the fine structure constant α and the Weinberg angle θ W . Since the q i represent the hypercharges of the SM fermions and the heavy fermion Ψ from our sector, after including all the contributing particles we obtain a factor of ∑ i q 2 i = 6. The second and the third diagrams represent contributions to the free energy from the Higgs doublet. We focus on the high temperature/low mass limit prior to electroweak symmetry breaking. The contribution to the amplitude from the second diagram can be written as the product of the amplitude of two independent one loop diagrams, In the T m → ∞ limit we obtain, Finally, the third diagram has a form very similar to the first (A11), except that now all the integrals are bosonic, (A16) Now, using Eqs. (A3) and (A6) we put all the pieces together and see that for T ≫ M , the temperature dependent contribution simplifies to, When the temperature falls below T ∼ M , the contribu-tion given by Eq. (A3) becomes exponentially suppressed while the correction to the gauge coupling from Eq. (A5) becomes independent of temperature. Additionally, the vectorlike fermion Ψ of our sector no longer contributes to the Eq. (A13), resulting in ∑ i q 2 i = 5. Therefore, when the temperature enters the range M ≳ T ≳ 100 GeV, the potential takes the form (A18) After electroweak symmetry breaking, we can neglect the contribution from the Z boson as it is suppressed almost immediately. To get the contribution from photons we multiply the temperature dependent part in Eq. (A6) by cos 4 θ W , which effectively replaces the hypercharge gauge coupling g ′ by e Here q 2 ef f (T ) is a coefficient that changes with the number of species that contribute to the temperature dependent potential. As the universe cools down the number of species that contribute to the potential drops as heavier fields go out of the bath. This in turn, increases the relative strength of Hubble friction. However, this effect is partially compensated for by the rise in the temperature when heavy fields exit the thermal bath.
Therefore, the thermal potential does not change qualitatively until the temperature drops below the mass of the lightest charged particle, the electron, so that T ≲ m e . However, in order for this effect to play any role in the evolution of the modulus ϕ, the thermal part V T (ϕ, T ) has to be significant at this point in time, which requires where η 0 was defined in Eq. (38). On the top of that, V T (ϕ, T ) has to dominate over V 0 (ϕ) at T ∼ m e . Hence, we have to satisfy If we additionally require that these conditions overlap with the parameter space not excluded by experiments, we conclude that for M ≳ 1 TeV no solutions exist that satisfy the above criteria simultaneously. It therefore suffices to consider the high temperature limit of the temperature dependent potential.
wheret is a time when T =T . Using the fact that this transition always occurs during the radiation dominated epoch, when T ∝ t −1 2 , we can write . (B9) The ratio T 4 ef f T 4 tr determines the shape of the potential in Eq. (B1) at T = T ef f . Eq. (B9) relates this ratio to the parameter η. Using these equations we can determine the ratio of the height difference between the central maximum and the closest minimum at T = T ef f relative to the same height difference at zero temperature, Since the lowest value of log 1 x for both N = 6 and N = 10 is log 1 x ≈ 5, in order to satisfy ∆h h > 0.9 we need η ≪ 70. This condition is easily fulfilled for N = 6 as η ≲ 40. However, for N = 10 it is satisfied only for smaller values of the mass of the modulus, m φ ≲ 10 −12 eV. Nevertheless, as evident from the numerical simulations presented in Fig. 4, Eq. (54) continues to give a good approximation to the actual result even when this condition is not satisfied.

Appendix C: Parametric Resonance
One may worry that decays of φ into photons may deplete the abundance of dark matter, particularly if the decay rate is enhanced by parametric resonance. In this appendix we consider this question. We find that in the region of parameter space of interest, the condition for parametric resonance is not satisfied. Therefore these decays are slow on cosmological time scales and the abundance of dark matter is not affected.
The Lagrangian for the coupling of φ to photons can be written as The resulting equation of motion for the photon is given by Working in the gauge where A 0 , ∂ µ A µ = 0 this becomes, This leads to the modified dispersion relation, Because the plasma mass of the photon is large at the time of the phase transition, we are interested in parametric resonance at times when φ f ≪ 1. Perturbatively expanding the time-dependent dispersion relation, we obtain When the φ condensate decays, the resulting photons have the frequency ω = m φ 2. In phase space, the photons are emitted into a thin shell centered around m φ 2 with width 2δk. The occupation number of photonsñ k γ with momentum k can be related to the number density of photons n γ as The integrated Boltzmann equation for the production of photons from the decay of the φ condensate in the Bose enhanced regime (ñ k γ ≫ 1) is given by, d dt (a 3 n γ ) = 2a 3 Γ(1 + 2ñ k γ )n φ ≈ 4a 3 Γñ k γ n φ = a 3 Γn φ 2πk 2 δk n γ (C7) where Γ denotes the decay rate of the condensate into photons at zero background photon density. This decay rate Γ can be estimated as For all of the benchmark points in Table II, the lifetime is greater than 10 75 years. This is much greater than the age of the universe. Therefore, in the absence of parametric resonance, decays of φ can safely be neglected.
The parametrically-resonant decay rate can be easily read off from Eq. (C7) as where we have used n φ = m φ φ 2 andφ φ = m φ . In an expanding universe the proper momentum is redshifting as k = k comoving a. This means that the radius of the thin momentum shell is shrinking at the rate of kH. Since the momentum shell has a finite width 2δk, it takes time t k = 2δk (kH) for an emitted photon to redshift outside the momentum shell. In order to avoid Bose enhancement of decays in the ex-panding universe, this time scale t k has to be parametrically smaller than the decay rate Γ res at the time of decay, At early times prior to photon decoupling, the photon has a plasma mass which is larger than m φ . Consequently, decays of φ into photons are kinematically forbidden until after decoupling. This means that the condition in Eq. (C10) only has to be satisfied at times after T decoupling ≈ 0.3 eV to avoid parametric resonance. Since decoupling happens during matter domination, the Hubble expansion at these times can be related to the field φ as where r is the fraction of dark matter constituted by φ. Then, the condition to avoid parametric resonance becomes, 2 √ r 2 e 4 3π 5 M pl φ f 2 ≪ 1.
This condition is easily satisfied because the amplitude of the field φ has already undergone significant damping by the time of recombination.