Cosmological relaxation from dark fermion production

We consider the cosmological relaxation solution to the electroweak hierarchy problem using the fermion production as a dominant friction force. In our approach, neither super-Planckian field excursions nor a large number of e-folds arise, and scanning over thermal Higgs mass squared is avoided. The produced fermions from the relaxion source through the derivative coupling are SM-singlets, what we call dark fermions, and they can serve as the keV scale warm dark matter candidates.


Introduction
Despite the Higgs discovery at the LHC [1,2] and the continuous measurements of its properties, the smallness of the Higgs mass still remains mysterious within the quantum field theory of the Standard Model (SM). The traditional approach to the naturalness problem such as the supersymmetry [3][4][5][6] or Higgs compositeness [7][8][9][10][11][12] relied on the symmetry-based selection rules and the light new particles such as top partners in minimal scenarios [13,14] were expected to be observed at the LHC. However, the absence of the evidence for the New Physics near the electroweak scale so far only degrades the naturalness principle as a valid guiding principle that has been used for decades to extend the SM [15].
The relaxation meachanism is a recent new approach to the naturalness problem based on the cosmological dynamics of the Higgs mass squared [16]. In the relaxation approach, the smallness of the Higgs mass parameter is not associated with the symmetry, and in principle no new physics near the electroweak scale is required to show up. After the JHEP02(2020)135 first realization of the relaxation mechanism for the electroweak hierarchy problem in [16], a series of improvements attempting to resolve downsides in the original scenario have appeared . Among them, the one using the particle production [28] is of particular interest. We postpone the overview of the original idea and its variant using the particle production to section 2.
In this work, we newly investigate the cosmological relaxation solution to the electroweak hierarchy problem using the particle production of fermions. The application of the fermion production [47,48] in the Beyond the SM (BSM) scenarios has been somewhat limited mainly due to the Pauli-blocking (unlike the parametric resonance for scalars [49] or the particle production of the tachyonic gauge bosons [50]). A recent interesting application in cosmology is the axion inflation in [51,52] where the backreaction from the fermion production was shown to be more efficient than the dissipation via the Hubble friction in supporting the slow-roll of the inflaton and from where we adopted many technical results. The goal of this work is to investigate the plausibility of realizing the cosmological relaxation using the fermion production as a dominant way of dissipating relaxion energy while aiming to maintain the cutoff scale in a similar size to that of other variant models. Two benchmark BSM scenarios that we consider in this work are the non-QCD model in [16] and the double scanner mechanism proposed in [17]. While the non-QCD model does not look satisfactory (not conclusive though), we use it as a toy example for the simpler demonstration of the underlying physics, and the double scanner mechanism will be taken as our proof-of-concept example. We will show that our proof-of-concept example avoids downsides in the original relaxation scenario, but at the same time it brings new types of theoretical challenges. We will also demonstrate that the fermions have to be SM singlet for the mechanism to work and their masses are interestingly in the right ballpark to be the dark matter candidates.
The paper is organized as follows. In section 2, we review the original relaxation mechanism in [16] and some of its improvements. In section 3, we survey two BSM models, namely non-QCD model and double scanner model, focusing on whether the backreaction from the fermion production can support a slow-rolling relaxion while satisfying all theoretical constraints. In section 4, we discuss about the prospect for dark matter candidates and check their compatibility with the current phenomenological and astrophysical bounds. The concluding remarks are summarized in section 5. In appendix A, we provide all detailed derivations of our analytic expressions that we have used throughout our work.

Ralaxation overview
In this section, we briefly review the cosmological relaxation of the electroweak scale proposed by Graham, Kaplan and Rajendran [16] (which we refer to as GKR) as well as the issues in the GKR scenario. Relaxion is an axion-like particle (ALP) whose discrete shift symmetry is softly and explicitly broken by a small coupling. The relaxion potential in GRK scenario takes the form ∆V = −Λ 2 + gφ |h| 2 + gΛ 2 φ + · · · + Λ 4 c cos (φ/f ) , (2.1)

JHEP02(2020)135
where φ is the relaxion field, h is the Higgs doublet, Λ is the cutoff scale of the relaxion model, and g is a small dimensionful coupling. The ellipsis in eq. (2.1) refers to the higher order terms in gφ/Λ 2 . The small potential terms for φ are technically natural in the sense that the discrete shift symmetry φ → φ + 2πf is restored when the coupling g → 0. The periodic potential Λ 4 c cos (φ/f ) arises from model-dependent non-perturbative dynamics (either SM QCD or non-QCD strong gauge group), and the height of the potential barriers Λ 4 c ∝ M 4−n v n , where the integer n ∈ [1,4], 1 M is a parameter of mass dimension, and v is the Higgs vacuum expectation value (VEV). The QCD relaxion is problematic since it generates an O(1) shift in the θ-term, causing the strong CP problem. Either an additional mechanism (see [16,41] for example) must be introduced to remedy the problem or one needs to consider a new strong gauge group instead of QCD.
The relaxion field initially can be anywhere, while φ Λ 2 /g, such that the effective Higgs mass squared is positive, µ 2 ≡ −Λ 2 +gφ > 0, and electroweak symmetry is unbroken. Since the cosine potential is switched off in the unbroken phase, the evolution of the relaxion is driven by the linear potential. The mechanism is insensitive to the initial value of φ since the relaxion is slow-rolling due to the Hubble friction. The relaxion must scan O(1) fraction of the field space to naturally pass a critical point φ c = Λ 2 /g where µ 2 changes its sign, and the Higgs develops a nonvanishing VEV. The increasing Higgs VEV backreacts onto the potential by growing the potential barriers. Eventually it results in a compensation between the slope of the periodic potential Λ 4 c /f and the linear slope gΛ 2 , stabilizing the relaxion in a local minimum near a small value, namely, the electroweak scale v. Therefore, the smallness of the Higgs VEV is explained by the dynamical evolution of the relaxion field instead of fine-tuning or anthropics.
Certain conditions must be satisfied to make the relaxation mechanism work. The slow-roll of the relaxion must be long enough such that it can scan an O(1) fraction of the whole field range, and it sets a lower bound on the number of e-folds N e H 2 /g 2 : ∆φ ∼φ∆t ∼φ N e /H ∼ (gΛ 2 /H 2 ) N e Λ 2 /g where we have used the slow-roll condition 3Hφ + d∆V dφ ∼ 0. The vacuum energy should be greater than the typical relaxion energy density, H 2 M 2 P Λ 4 . The potential barrier must be formed within the Hubble sphere, H −1 > Λ −1 . In addition, it has to be the classical rolling instead of the quantum spreading that sets the vacuum into the correct one, leading to H 2 /φ < 1. After the relaxion stops rolling and the reheating of the universe occurs, a sufficiently high temperature may erase the barrier and cause the relaxion to roll again. Hence, either the reheating temperature must be low enough such that the barrier does not get melt or the traveling distance during the second rolling must be short enough not to overshoot the electroweak scale. Some downsides, however, exist for the original GKR scenario which relies on the Hubble expansion to dissipate the relaxion energy. The typical size of the dimensionful coupling g is tiny (e.g. g ∼ 10 −31 GeV for the QCD axion model), and it leads to the super-Planckian relaxion scanning during the inflation and the exponentially long e-folding. The fact that those downsides are linked to the inflation suggests looking for a more efficient JHEP02(2020)135 ψ ψ ϕ Figure 1. A cartoon of the relaxation with the dark fermion production.
way of dissipating the relaxion energy than via the Hubble friction. Such an example is found in the particle production sourced by the relaxion field. In ref. [28], Hook and Marques-Tavares (HMT) utilized the particle production of tachyonic gauge bosons through φFF as an efficient way of dissipating relaxion energy. The exponential production of the electroweak gauge boson naturally occurs to stop the relaxion at a value that generates the electroweak scale. Since the mechanism can work in a non-inflationary era, the issues mentioned above can be avoided. However, unlike the GKR scenario, the original HMT mechanism operates in the broken phase of the electroweak symmetry. Such a model requires a specific, usually non-trivial UV completion, e.g. from a left-right symmetric model as in the appendix of [28]. The cutoff scale of the HMT scenario is relatively low, up to 10 4∼5 GeV, and hence this scenario addresses the little hierarchy problem (see also [53,54] for more recent and updated analyses on the HMT mechanism).

Ralaxation from dark fermion production
In this section, we investigate the plausibility of the fermion production as an alternative to the Hubble friction as a dominant way of dissipating relaxion energy to achieve the slow-rolling of the relaxion. The cartoon picture for the situation is illustrated in figure 1. As in the original GKR scenario, we rely on the cosine potential being switched on only in the broken phase of the electroweak symmetry to stop the relaxion at the right place while a large field excursion of the relaxion occurs in the unbroken phase. Since the fermion production can not be exponential due to Pauli-blocking (unlike the case of the tachyonic gauge boson), it would be difficult to be implemented to work in a similar manner to [28] in the broken phase as a way of naturally selecting the electroweak scale.

Dark fermion production
The fermionic system that couples to the relaxion field φ through the derivative coupling in an expanding Universe can be written as

JHEP02(2020)135
where e µ a is a vierbein and D µ is the covariant derivative due to the spin connection from the scale factor a in the metric, The overall dependence on the scale factor a in the fermionic system can be removed via rescaling, ψ → a −3/2 ψ, after which the Lagrangian density for the fermions in the conformal time takes a simple form, Throughout this work, we will assume that the relaxion field φ is spatially homogeneous.
As was pointed out in [52] (see [55] for the related discussion), the fermion production in the basis with the derivative coupling as in eq. (3.3) is not accurately estimated. A prescription is to go to a new basis, where the fermion production is unambiguously estimated, via the rotation [52] ψ → e −i γ 5 φ/f ψ ψ . (3.4) In the new basis, the interaction between the relaxion and fermions takes the form, where m R = m ψ a cos(2φ/f ψ ) and m I = m ψ a sin(2φ/f ψ ). One clearly sees in eq. (3.5) that the fermion production vanishes in the massless limit since the fermions become free fields. The new basis is more suitable to define the fermion number mainly because the Hamiltonian for the fermion takes a simple quadratic form, where ψ is taken to be a quantum field in terms of the creation and annihilation operators. The fermion number density for a momentum k and helicity r is given by the vacuum expectation value of the number operator at a finite time, namely 0|a † r ( k)a r ( k)|0 . The detailed computation will be given in appendix A. Here, we simply take the final result to discuss about the main feature of the scenario. The fermion ψ in eq. (3.3) can not be a SM one since the SM fermion will be massless during the scanning era in the unbroken phase of the electroweak symmetry. We assume that ψ is a SM singlet massive Dirac fermion (what we call the dark fermion in the following) and that its mass is a Higgs VEV independent. 2 The non-SM fermion ψ could be a candidate for a dark matter as a bonus. We will investigate its plausibility in detail in section 4.

Slow-rolling from backreaction
In presence of the backreaction due to the fermion production, the equation of motion for the relaxion readsφ where V φ ≡ ∂V /∂φ and dot denotes the differentiation with respect to the cosmic time, for instance,φ = ∂ t φ. The backreaction B from the fermion production in eq. (3.7) is given by The exact evaluation of the backreaction B can be found in appendix A. In the limit µ 2 ≡ m 2 ψ /H 2 ξ and ξ 1, it is approximated to be where ξ is defined as the ratio of the velocity of φ to the scale f ψ in the Hubble unit, The parameter ξ plays a crucial role in controling the slow-rolling of the relaxion, or the strength of the fermion production or energy dissipation. A strong fermion production occurs when the adiabatic condition is strongly violated. This implies that the relaxion should maintain a sizable velocity to strongly depart from the adiabaticity. What controls the size of the velocity of φ would be the slope of the linear potential, namely g in our framework. Having said that, the fermion production assisted slow-rolling favors a bigger slope than that in the original GKR relaxation operating with the Hubble friction, whereas the slope is smaller compared to the case of the inflation through the fermion production in a steep axionic potential [52]. While the velocity of the slow rolling φ during the inflation in the absence of the backreaction is determined by equating the second and third terms on the left hand side of eq. (3.7), the velocity in our scenario is engineered to be determined by equating V φ (φ) with the backreaction term B in eq. (3.7): It implies that the size of the backreaction is solely determined by the size of the slope g as the left hand side of the slow roll equation in eq. (3.11) depends only on the slope for a given cutoff scale Λ.

Theoretical constraints for relaxation
In the following, we list various constraints for the successful relaxation from the dark fermion production, and we express them as a lower or upper bound on the dark fermion mass.
which is easily satisfied in the parameter space of interest.
In addition to the constraints listed above, when the relaxion is scanning over the effective Higgs mass squared, the temperature is required to be negligible compared to the electroweak scale, not to scan over the thermal Higgs mass squared. We help this issue by considering our relaxation mechanism during inflation driven by a separate inflaton sector (where the inflationary Hubble H is much higher than Λ 2 /M P ). The constraint can be avoided when the interaction rate between the relaxion and the SM sector can be made small enough to be diluted during the inflation. 3 The cosine potential for φ is turned on when φ passes the critical point φ c from which it enters into the broken phase of the electroweak symmetry. The relaxion is being trapped in one of minimum when the slope of the linear potential is balanced with the slope of the cosine potential, does not lead to the constraint that fits to the 3 We find that this is plausible only in the double scanner mechanism among two scenarios in sections 3.4 and 3.5. The dark fermions may thermally produce relaxions even during the inflation (it is a characteristic feature of the derivative coupling that induces strong backreaction at the higher energy or temperature) in both scenarios. However, the small non-derivative couplings between the relaxions and the SM sector in the double scanner mechanism lead to an inefficient interaction rate than H during the inflation. A similar suppression is not obvious in non-QCD model as the relaxion may thermally produce new non-QCD gauge bosons G above the confinement scale through the derivative coupling φG G , and then the new gauge boson will thermalize the SM sector. form in eq. (3.21). Using the inequality f > Λ in eq. (3.20) leads to g < Λ 4 c /Λ 3 , and it can be combined with the condition for sub-Planckian field excursion ∆φ ∼ Λ 2 /g < M p to obtain another upper bound on the cutoff,

JHEP02(2020)135
The excluded region in (Λ c , Λ) plane from eq. (3.21) for N e = 10 2 is illustrated in figure 2 where the strongest bound corresponds to the first one (red region) in eq. (3.21) followed by the third (blue region) as next strongest, and the second one is the weakest (not shown). We also show in figure 2 the case from eq. (3.22) (gray region) which roughly matches to eq. (3.21) for N e = 1. As is evident in figure 2, the higher barrier is more compatible with the higher cutoff scale.
In the next section, we will survey a few new physics models that predict different forms of Λ c , and we will determine the viable parameter space consistent with all the constraints in those models. While one could explore the parameter space with the maximally allowed cutoff scale within the ballpark for Λ c in a specific model, throughout our work, we will fix the cutoff scale to Λ ∼ 10 4∼5 (aiming to address only the little hierarchy problem) as our benchmark point. Note that another type of relaxation scenario with the particle production in [28] also has a similar low cutoff scale.

Non-QCD model
The model that we try first as an illustration is the non-QCD model [16] supplemented by dark fermions. New massive fermions L, N (and their conjugates L c , N c ) in the non-QCD model are charged under the new gauge group that gets strongly coupled in the low energy JHEP02(2020)135 scales, and they have Yukawa-type couplings with the Higgs field: (3.23) where L and N have the same electroweak quantum numbers as those of the lepton doublet and right handed neutrino. The L and L c fermions must be heavier than the electroweak scale to avoid the phenomenological constraints whereas N and N c can be made very light such that it can form a condensate below the confinement scale ∼ 4πf π . The hierarchy of m L f π m N will be assumed as in [16]. Assuming that the anomalous interaction (φ/f ) G µνG µν is allowed in the model, it can be traded for the phase of the mass term m N via the chiral rotation for N , where m N collectively refers to not only the bare mass in eq. (3.23) but also all kinds of corrections to the mass. An analogous term to eq. (3.24) will generate the periodic potential for φ of the type ∼ Λ 4 c cos(φ/f ) through the condensate N N c ∼ 4πf 3 π below the confinement scale. While m N in eq. (3.24) gets various contributions, the h-dependent contribution at tree level is estimated to be where f π is the chiral symmetry breaking scale of new confining gauge group. For the mechanism to work, all h-independent contributions to m N must be subleading [16]: 26) and it implies that the masses of L and L c are at order of a few hundred GeV (lighter masses will be constrained at the LHC). As was mentioned above, demanding that N and N c should be lighter than the confinement scale gives rise to Using the relation in eq. (3.20) with the expression of Λ c in eq. (3.25), the decay constant in the cosine potential can be expressed as and it needs to be bigger than the cutoff scale, f Λ, to be consistent with the EFT. The viable parameter space for the relaxation from the dark fermion production to work is illustrated in (g, m ψ ) plane in the left panel of figure 3, while fixing other parameters such as Λ, H, f ψ , and N e . A different choice of those parameters leads to a different allowed region in (g, m ψ ) plane. In the right panel of figure 3, we illustrate the allowed region of the parameters (from eqs. (3.27) and (3.28)), which are specific to the non-QCD model, Figure 3. The viable parameter space for non-QCD model. The light-red region in bottom-right corner of the first panel corresponds to the excluded region from the small energy density carried by the dark fermion. The constraint from the 'classing rolling beats the quantum spreading' is too weak to show up in the plot. We set Λ = 10 4 GeV, H = 5 × 10 −6 GeV, f ψ = 0.5 GeV, and N e = 100 in both panels (and g = 10 −6 as well for the right panel). The approximation µ 2 ξ is not valid in the region above the dashed line. in the (yỹ/m L , f π ) plane for the cutoff Λ = 10 4 GeV. Based on the result in figure 3, we present one benchmark point in table 1 for the illustration. The large field excursion of φ is definitely sub-Planckian, ∆φ ∼ Λ 2 /g < M p , for the range of g in figure 3 and cutoff scale Λ ∼ 10 4 GeV.

JHEP02(2020)135
If the energy released from the strong fermion production can be transferred to the visible sector with the new strong group, the barrier Λ c will disappear during the reheatingera and the relaxion will start rolling again, spoiling the mechanism. This unwanted property has led to non-trivial constraints on the model [16]. The situation is worse in our scenario with the dark fermion production as there could be a chance that the SM sector might be thermalized even during the inflation, spoiling the entire mechanism. One way to resolve these issues can be found in the so-called double scanner mechanism [17] where the new strong gauge group is engineered to get strongly coupled at the same scale as the cutoff value. The detailed analysis of the double scanner mechanism in the context of the strong dark fermion production will be the subject of the next section.

Double scanner mechanism
The double scanner mechanism [17] introduces an additional slow-rolling field σ whose main role is controlling the amplitude of the cosine potential, while scanning the Higgs mass parameter is still carried out by the original relaxion field φ. The relevant part of the potential is given by where we will assume the universal decay constants f ψ = f σ for φ and σ for simplicity. The amplitude A(φ, σ, h) in eq. (3.29) is given by where is a supurion that accounts for the shift symmetry breaking. Recall that the confinement scale in non-QCD model has to be much lower than the cutoff scale to suppress all h-independent contribution to the cosine potential individually (see eq. (3.26) for instance). In the double scanner mechanism, the confinement scale can be made as big as the cutoff scale, Λ con ∼ Λ, while keeping h-dependent contribution in eq. (3.30) as the dominant term. The newly introduced slow-rolling field σ cancels all h-independent terms of order ∼ Λ 4 in eq. (3.30) together during the cosmological evolution of interest with the appropriate choice of coefficients, β, c φ , and c σ of order one, and this cancellation is technically natural. The double scanner mechanism features multiple stages of cosmological evolution in (φ, σ) plane. In stage one, two fields φ and σ start evolving at some initial points while φ Λ/g and σ Λ/g σ such that the effective Higgs mass squared parameter is positive. The electroweak symmetry is unbroken as usual. Since the amplitude has a generic size of A ∼ Λ 4 , the potential for φ is dominated by the A cos (φ/f ) which causes φ to get stuck at some minimum. The σ field continuously rolls down while scanning the amplitude A. At some point, the cosine potential for φ becomes smaller than the linear potential for φ, and φ starts rolling down the linear potential (the second stage begins). The field φ continues rolling down along the trajectory, represented by φ * in [17], along which the evolution is dominantly driven by the linear potential. The second (third) stage ends (begins) when φ passes the critical point, φ c = αΛ/g, after which the sign of the effective Higgs mass squared flips, and the electroweak symmetry breaking is triggered. In stage three, the h 2 term in eq. (3.30) is switched on, and it dominates the amplitude. The amplitude scaling as ∼ h 2 keeps growing as φ keeps rolling down the potential since the negative Higgs mass squared term increases. The relaxion field stops rolling and it is trapped in one of the cosine potential wells when the steepness of the cosine potential is balanced with the slope of the linear potential: In the last stage, the field σ keeps moving down the potential 4 until it finds its minimum somewhere. Around the end of the last stage, we do not expect any cancellation in the amplitude A(φ, σ, h), and the typical size of the amplitude will be A ∼ Λ 4 which implies that the natural size of φ mass without fine-tuning is expected to be Whereas the mass of σ is given by For the successful scanning over φ tracking σ along the sliding trajectory φ * (where the cosine potential is smaller than the linear potential) before it reaches the critical point, φ c , the condition dφ(t)/dσ(t) = (g/g σ ) 1/2 > dφ * /dσ should be satisfied. 5 Otherwise, φ gets deviated from the trajectory before it reaches φ c , and it gets stuck at some minimum. To avoid this situation, we require Once it passes the critical point where the Higgs barrier is switched on, φ needs to exit the trajectory to continue its evolution along the path where the amplitude grows like [17] for the related discussion). Theoretical constraints for the double scanner mechanism to work are similar to those in non-QCD model in section 3.4. The constraints from eq. (3.13) to eq. (3.18) similarly apply to σ with g σ , and we take a stronger one in the numerical simulation to determine the viable parameter space. Since the cosine potential for φ contributes to the Higgs mass squared of order ∼ m 2 φ , we demand m 2 φ v 2 not to reintroduce the naturalness problem. Using the expression in eq. (3.32), we obtain the upper bound on g, While the double scanner mechanism relies on the field σ which cancels the amplitude in eq. (3.30), the cosine potential gets quantum corrections at quadratic order on cos(φ/f ) such as 2 Λ 4 cos 2 (φ/f ) etc [17]. The mechanism is spoiled unless those corrections to Higgs barrier (∼ Λ 2 v 2 ) remain subleading, and ignoring the loop factors, it requires Combing eq. (3.31) and eq. (3.36), we obtain another upper bound on g, where the last inequality is due to f Λ which also makes the above constraint stronger than eq. (3.35). The parameter space consistent with all the constraints mentioned above in (g, m ψ ) plane is illustrated in figure 4, and a benchmark point is presented in table 2. The cutoff scales were fixed to Λ = 10 4 or 10 5 GeV in figure 4 as our aim is to address the little hierarchy problem (see caption of figure 4 for other chosen parameters).

JHEP02(2020)135
As was explained in section 3.3, the inflation is driven by a separate inflaton sector to avoid scanning the thermal Higgs mass squared during the scanning era. This idea works for the double scanner mechanism due to the confinement scale as big as the cutoff scale unlike the case of the non-QCD model. However, the successful double scanner mechanism from the fermion production also has to be safe against being thermalized to the temperature above Λ during the reheating era, and it imposes a constraint on the inflaton sector. In this work, without getting into the details of the reheating, we assume that the energy stored in the inflaton sector is smaller or at most comparable to that of the relaxation sector (collectively denotes the entire sector of ψ, φ, SM, and the strong gauge group) when the relaxation sector gets in thermal equilibrium. 6 We also assume that the end of relaxation coincides 6 A possibility to make the inflaton sector energy density to be subdominant is to require that the inflaton decays into radiation after inflation and that this timescale is sufficiently short compared to the thermalization timescale of the relaxation sector (see [56] for a recent discussion about reheating in different sectors). In a situation that the reheating temperature in the inflaton sector is higher than the cutoff scale Λ and the relaxation sector is in thermal equilibrium with the inflaton sector during such a short time scale, thermal effects would erase the periodic potential barriers, leading to the stabilized relaxion to roll again. The inflaton sector and its interaction to the relaxation sector have to be constrained such that either the second scanning era itself does not happen or it does not overshoot the electroweak scale during the second scanning era (if it occurs). Another possibility is to consider a scenario where its energy density decays faster than the radiation [57,58]. with the end of inflation such that the fermion energy density ρ ψ would not be diluted away by inflation to allow the possibility of ψ as a dark matter candidate. This can be achieved by appropriately setting up the inflationary sector and choosing parameters therein.

Comments on f ψ Λ
Apparently, the scale f ψ in the derivative interaction in eq. (3.3) needs to be much smaller than the cutoff scale of the model, namely f ψ Λ, to make the fermion production efficient enough. Two hierarchical scales without a natural explanation can be considered to be either inconsistent from the EFT point of view or a strong coupling problem when f ψ is normalized to Λ. It has been a generic issue in applications which heavily rely on the fermion production as a main dissipation, and we are not an exception. While we do not have a solution for this issue, we will briefly comment on the possibility of the thermal decoupling between the dark fermion sector and the visible sector with the new strong gauge group.
The thermal decoupling is an attractive idea, if it can be realized, for two theoretical issues: the second scanning in the reheating era and two apparently inconsistent scales within an EFT. If the interaction rate between two sectors can be made negligible over the entire range of temperature with respect to the Hubble parameter, the dark fermion sector will have its own thermal history without interfering with the visible sector. If this is the situation, the barrier of the cosine potential will not be erased. Besides, it will be natural for two thermally disconnected sectors to have their own cutoff scales, and the EFT of one sector would not spoil the validity of the EFT of the other sector.
From our numerical simulation, we find that the interaction rate between two sectors via the intemediate φ-exchange such as diagrams in figure 5 can be made smaller than the Hubble parameter. It is basically because the diagrams are doubly suppressed: small m ψ parameter from the derivative coupling and h-φ mixing angle. There also could be processes that decay into the states of the new strong gauge group. Based on our naive dimensional analysis (NDA), we find that the interaction rate can be made smaller for the double scanner scenario (marginally smaller for the non-QCD model) than the Hubble parameter although it has to be confirmed by more exact numerical simulation.
The major threats to the idea of the thermal decoupling come from the processes such as those in figure 6. For instance, the dark fermion can first thermally produce the relaxion φ whose lifetime is not shorter than the typical interaction timescale. Then, those φ fields can thermally produce the SM particles. The typical diagrams such as those in figure 6 have only one suppression which is not enough to make the interaction rate smaller than JHEP02(2020)135 the Hubble parameter. 7 In this situation, the energy carried by the dark fermion will be efficiently transferred to the visible sector. On the other hand, when the temperature drops below m φ , the diagrams involving φ in figure 6 will be exponentially suppressed as φ becomes non-relativisitc particle, and thus two sectors will be thermally decoupled.
The above discussion suggests that one need to suppress those diagrams in figure 6 to achieve the thermal decoupling over the entire temperature range. One might try to dilute the energy density released by the dark fermion away below m φ such that the diagrams in figure 6 are exponentially suppressed. We leave it for the future work.

Prospect for dark matter
We investigate the compatibility of two models in sections 3.4 and 3.5 with the current phenomenological and astrophysical constraints.

Relic abundance and dark matter
Since the Higgs field dominantly mixes with φ, we diagonalize only the 2 × 2 mass matrix of (h, φ), and we treat σ separately.

JHEP02(2020)135
where the mixing angle is given by where we defined y ≡ 2 m hφ /(m 2 h − m 2 φ ), and m hφ = gv. Since m 2 h > m 2 φ is required not to reintroduce the fine-tuning, the mixing angle is roughly ∼ 2 m hφ /m 2 h ∼ 2gv/m 2 h . The decay rate of φ in non-diagonalized basis is given by where Γ φ→ψψ = 1 2π It is important to notice that the derivative coupling gives an effective coupling of m ψ /f ψ instead of E-growing effective coupling E/f ψ . One can see that the suppression by m ψ /E is originated from the equation of motion of fermions in the derivative coupling. The decay rate for σ is obtained by the replacement φ ↔ σ. While φ can always decay into dark fermions before the Big Bang Nucleosynthesis (BBN), the σ → ψψ channel might not be kinematically allowed in some allowed parameter space where m σ ∼ g σ < 2 m ψ . Once the σ → ψψ decay channel is kinematically opened, the σ field decays into dark fermions well before BBN. Otherwise, σ mainly decays into SM fields via the mixing with the Higgs, and its lifetime is longer than the age of the universe. Therefore, the candidates that could potentially serve as dark matters in our scenario are ψ and σ. The abundance of σ can get a contribution from the vacuum misalignment. The energy density ρ σ at the start of the oscillating regime is ∼ m 2 σ (∆σ) 2 where the misalignment of σ is given by ∆σ ∼ √ N e H. 8 The relic abundance for N e ∼ O(10 2 ) is expected to be negligible, where T osc = m σ M p is the temperature below which σ is in the oscillating regime. We conclude that the non-thermally produced σ has negligible contribution to the current dark matter abundance. The σ field can also be produced from the thermal bath via ψψ → σσ (via t and u-channels) process. The production of σ from the φ scattering is negligible. The order of magnitude estimate of the thermally averaged interaction rate, ignoring all numeric factors, in the limit T m ψ , m σ is roughly Figure 7. The thermally averaged interaction rate of ψψ → σσ, Γ ψψ→σσ (T ), (black) for two benchmark points in table 2 and the Hubble rate (red).
The thermally produced σ field is expected to decouple from the thermal bath below the temperature, which looks close to m ψ and m σ especially for the case of our benchmark point with Λ ∼ 10 4 GeV in table 2. The more exact numerical evaluation of the thermally averaged interaction rate is illustrated in figure 7 for two benchmark points in table 2, and the decoupling temperatures are found to be roughly two orders of magnitudes higher than our rough estimate in eq. (4.7). It implies that the thermally produced σ gets out of equilibrium while being relativistic. Among the annihilation channels of ψ into a pair of σ, φ, or SM particles, the dominant channels are ψψ → φφ, σσ (the channel to φφ will be shut off below m φ though). Similarly ψ decouples from the thermal bath while being relativistic. Both ψ and σ (or only ψ if σ decays before BBN) could be warm dark matter (WDM) candidates, and their abundances today are, assuming no entropy dilution with the constant relativistic degrees of freedom from its decoupling epoch, where we have used the fact that after decoupling the photon temperature T γ and the effective temperature of σ, ψ evolve identically due to the constant relativistic degrees of freedom, namely T γ = T σ, ψ . For our keV-scale σ and ψ, the relic abundance in eq. (4.8) is three orders of magnitude higher than what is needed to be consistent with the present relic abundance. This constraint is in a tension with those from the Lyman-α forest analysis on the free-streaming length as they put a strong lower bound on the thermal WDM mass m WDM 5keV [60][61][62]. If the warm dark matter is a partial component of the whole dark matter, the Lyman-α bound can be relaxed [63].

JHEP02(2020)135
The tension between the relic abundance (4.8) and the Lyman-α bound in thermal WDM can be ameliorated in our scenario as follows. For σ, the easiest solution will be to make it heavier than 2m ψ so that the channel σ → ψψ is kinematically opened and σ decays into dark fermions before the BBN. This option will be viable for the first benchmark point with the cutoff, Λ ∼ 10 4 GeV, in table 2 after a slight modification of parameters, but it will be difficult for the second benchmark point with the higher cutoff. For ψ (and σ if the decay to ψ is forbidden), for instance, one can modify the strength of ψψ → σσ to significantly increase the decoupling temperature such that the discrepancy of the effective degrees of freedom at the decoupling temperature T d and at T 0 leads to (up to) two-order of magnitude suppression. It can be achieved by adapting non-universal fermion couplings to φ and σ such as setting f σ = Λ while keeping f ψ as before. In this situation, the slow-rolling of σ will be maintained via the Hubble friction 9 while the fermion production is still responsible for the slow-rolling of φ. This choice makes the interaction rate Γ ψψ→σσ (T ) always smaller than H within the cutoff scale, and thus the decoupling temperature will be roughly m φ below which Γ ψψ→φφ is switched off. Using the total entropy conservation of the Universe, one can estimate g * S (T 0 )/g * S (T d ) ∼ O(10 −2 ) for T d ∼ m φ ∼ O(10 1−2 ) GeV. Consequently, the thermal WDM mass bound from the Lyman-α constraints gets slightly relaxed due to the change in the relativistic degrees of freedom at the decoupling temperature, Although predicted by our relaxation parameter space, the keV-scale dark fermion mass is better consistent with this bound. Similarly, the relic abundance of ψ is diluted as where we have introduced an extra suppression factor 1/S from the entropy dilution that might be originated from a short-period entropy injection right after the decoupling of ψ. It can help lowering the relic abundance to the acceptable level [64]. Since the interaction rates between σ and φ with f σ = Λ are negligible compared to H, the σ fields are likely non-thermal with a negligible abundance.
Before concluding this section, we briefly point out one significant difference between our fermion warm dark matter and the standard sterile neutrino warm dark matter. The dark fermion in our model is not only a SM singlet, but also a stable particle. The sterile neutrino, however, can decay. Apart from decaying into three left-handed neutrinos when the sterile neutrino mass is lower than twice of the electron mass, the sterile neutrino has the loop-level decay channel: N → ν + γ [65], which could be detectable from the X-ray observations [66,67].

Conclusion
We have investigated a scenario of the cosmological relaxation of the weak scale supported by the backreaction from the dark fermion production during inflation, with the cutoff scale 10 4∼5 GeV. Being a more efficient friction source than the Hubble expansion, the fermion production plays a significant role in removing downsides in the original GKR scenario. Hence, our models do not have any extremely small parameters, the number of e-folds is of appropriate size, and the relaxion field excursion is sub-Planckian.
A characteristic property of the models with the fermion production through the derivative coupling is the possible thermalization between the produced fermions and the axionic source field even during inflation. While this could be alarming in a relaxion model in which scanning over the thermal Higgs mass squared should be avoided, the double scanner scenario can survive (while the non-QCD model might not) by considering the relaxation during the inflation whose inflationary expansion rate is higher than Λ 2 /M p such that the thermal relaxion cannot thermalize the visible sector during inflation. Another generic feature, or unwanted downside, of those models lies in the appearance of a scale f ψ , associated with the axionic derivative coupling to the fermion, much smaller than the EFT cutoff scale. This may lead to strong coupling, non-perturbativity, or the EFT inconsistency problem, and should be solved separately by an extra mechanism. We leave it to the future work.
Our model allows the possibility of the observation-consistent keV scale warm dark matter. While the fermion mass parameter controls the strength of the interaction in various places in our model, its preferred value for our scenario to work is coincidently in the vicinity of the right ballpark for the keV scale warm dark matter.

Acknowledgments
We would like to thank Pedro Schwaller and Wan-il Park for useful discussions.

The model
We consider the theory of the relaxion coupled to the SM singlet Dirac fermion through the derivative coupling (see eq. (3.1)), on the FRW metric where a(t) is a scale factor of the Universe. We assume that the pseudo-scalar is spatially homogeneous. More explicit form of the action in eq. (A.2) on the FRW metric is The overall scale factor due to √ −g in the Lagrangian for fermions can be removed via rescaling, ψ → a −3/2 ψ. Under this rescaling, the covariant derivative due to the spin connection become partial derivative, and the resulting Lagrangian becomes The Lagrangian in eq. (A.5) is problematic for the estimation of the fermion occupation number since the fermionic part in the corresponding Hamiltonian is not entirely captured in the quadratic form in ψ: where m R = m ψ a cos 2φ f ψ and m I = m ψ a sin 2φ f ψ . The corresponding Hamiltonian is given by, Since the Hamiltonian takes a quadratic form in ψ, one can unambiguously estimate the fermion occupation number.

JHEP02(2020)135
A.3 Fermion production The first step to estimate the fermion occupation number is plugging the expression for the quantum field for ψ, in the Hamiltonian while keeping the relaxion as a classical source. The spinor function in the helicity basis can be parametrized as where two-component spinor χ r with the helicity r satisfies the eigenvalue equation, σ · k χ r = rk χ r , and u r and v r represent the relative amplitudes between two different chirality states. The Hamiltonian in terms of the creation and annihilation operators is where the expressions for A r , B r in terms of u r and v r can be found in [52], The Hamiltonian in eq. (A.11) can be diagonalized by the unitary transformation, where the mixing angles α r and β r are called Bogoliubov coefficients which are functions in terms of u r and v r . The Hamiltonian in eq. (A.11) has two energy eigenvalues, The occupation number is given by More details can be found in [52]. One notes that the fermion occupation number as a function of the cosmic time t is the same as the one in terms of τ . Alternatively, the occupation number can be derived purely based on the group theoretic property as recently shown in [55]. In the group theoretic approach in [55], the authors showed that the expression in eq. (A.15) collapses into the SO(3), 10 invariant form,

JHEP02(2020)135
where two SO(3) vectors ζ r and q are given by and they are subject to the equation of motion, which looks similar to, for instance, that of the classical precession motion. The equation of motion in eq. (A.18) reproduces the same result as the one from the traditional approach.

A.4 Solution of equation of motion
While the fermion production is unambiguously defined with the Hamiltonian from the Lagrangian in eq. (A.7), the analytic solution is more clearly obtained in the basis with the derivative coupling in eq. (A.5). The equation of motion from the Lagrangian in eq. (A.5) is where dot denotes the differentiation with respect to the cosmic time, for instance,φ = ∂ t φ.
Using the relation σ · kχ r = rkχ r , the equation of motion in eq. (A.19) reduces to those in terms ofũ r andṽ r , and they are To convert the equation of motion into more familar form where analytic solutions are manifest, we introduce a new set of variables, Then, the equation of motion in terms of a new set of variables becomes whereũ r andṽ r were used to refer to the basis with the derivative coupling while keeping u r and v r for the basis without the derivative coupling. We take two linear combinations, JHEP02(2020)135 s r = (ũ r +ṽ r )/ √ 2 andd r = (ũ r −ṽ r )/ √ 2, and iterate two coupled first-order equations to get two decoupled second-order differential equations, After making a rescaling ofs r = x −1/2 s r (andd r = x −1/2 d r ) and changing the variable, x = −z/(2i), the differential equations take the form of the Whittaker equation.
We choose the boundary conditions such that the fermion occupation number in eq. (A.15) vanishes in the limit x → ∞ (or equivalently τ → −∞ or t → −∞ in terms of the conformal time or the cosmic time). One notes that fermion occupation numbers in both basis (with and without the derivative couplings) become identical in the limit x → ∞ (see section 5 of [55] for the detail). Zero occupation number in the basis with the derivative coupling guarantees the same boundary condition in the basis without the derivative coupling. Therefore, we can safely take the solutions obtained in the basis with the derivative coupling, namely (A. 25) and use them in the new basis. On the other hand, the scanning process in the relaxation mechanism occurs during the finite time, roughly ∆t ∼ N e /H. While it implies that the initial condition for the zero occupation number should be imposed in principle at a finite time τ 0 ( = −∞), the approximate solutions in eq. (A.25) must be sufficient since the above time interval in the cosmic time corresponds to the exponentially separated time interval in the conformal time, Although we have checked numerically that the result assuming the finite τ 0 is similar to the asymptotic case with τ 0 → −∞, we provide a semi-analytic proof for the validity of our approximation in eq. (A.25).
Since the velocityφ is positive, the occupation number is dominated by fermions with the helicity r = −1 since the fermion production can actively occur whenever

JHEP02(2020)135
and it is easily satisfied for the helicity r = −1 around k ∼ −2ξ/τ . Away from k ∼ −2ξ/τ , the WKB approximation is valid, and the occupation number is roughly constant. Even in the absence of the interaction of the relaxion to fermions, the occupation number is non-vanishing since a free massive fermion system in the expanding Universe can cause the fermion production. The latter contribution is not captured by the relation in eq. (A.27) as it is a gravitational effect, and we find that it shuts off around k ∼ −µ/τ . Consequently, the occupation number for the helicity r = −1, assuming an initial condition at τ 0 → −∞, can be written approximately as where the momentum k is the one in the frame with the conformal time and µ 2 /ξ in the region of interest is typically much smaller than 1/2. The numerical validation of the approximated occupation number in eq. (A.28) is shown in figure 8 for two choices of ξ values, ξ = 10, 100, for the purpose of illustration (see also figure 1 of [52] or figure 4 of [55] for related discussions). While the numerical simulation with a larger ξ in our benchmark points is technically challenging due to the highly oscillatory behavior, we suspect the generic feature remains the same. where W (1) (W (2) ) is the Whittaker function of the first (second) kind. The coefficients A, B in eq. (A.29) are found to be functions of three dimensionless parameters, It can be shown that A, B converge into their asymptotics at x 0 → ∞ (as those in eq. (A.25)) as long as x 0 = −kτ 0 2ξ. Since τ /τ 0 1 and the parametrization in eq. (A.28), the contribution to the fermion production from the momentum interval k = (−2ξ/τ 0 , −2ξ/τ ) will agree with our approximation with eq. (A.25). It is the contribution from the interval k = (0, −2ξ/τ 0 ) that might potentially invalidate our approximation.
Defining an intermediate momentum scale, and taking n r (k < k int ) = 1 as a conservative choice, our approximation with the solutions in eq. (A.25) will be valid as long as the following inequality for the number densities JHEP02(2020)135 integrated over each momentum interval is satisfied, where the contributions are dominated by the case with r = −1 and we dropped the factor from the integration over the solid angle. Demanding inequality in eq. (A.32) gives rise to the constraint, The constraints in eq. (A.33) are satisfied in our benchmark points as long as the number of e-folding is bigger than ∼ O (10), namely N e O(10). Therefore, we conclude that our approximation using the solutions in eq. (A.25) is justified.

A.5 Backreaction
The equation of motion for the relaxion in presence of the coupling to the fermion ψ is obtained in the basis without the derivative coupling, and it is given bÿ where the term in the right-hand side is the backreaction, what we call B, due to the fermion production. Using the expression for the quantum field in eq. (A.9), it is given by in terms of u r and v r (also in terms ofũ r andṽ r , quantities in the basis with the derivative coupling), (A.35) The above expression differs from that in [52] by the overall factor of 2. We suspect that it is due to the missed normalization factor 1/ √ 2 in V r (k, t). Using the relations, u r = (s r + d r )/ √ 2x andṽ r = (s r − d r )/ √ 2x in section A.4, the backreaction can be expressed in terms of s r and d r whose analytic solutions were the Whittaker functions, The integration is non-trivial, but it can be done analytically using the integral representation of the Whittaker functions known as the Mellin-Barnes representation, where the contour C covers from −i∞ to i∞ with a deformation to separate the poles of We have checked the result in [52]. While we have a few disagreements with those in [52] in some detail (possibly due to typos), we have reproduced the same result as [52] up to a factor of 2 that was mentioned above. Since all necessary details can be found in [52], we do not repeat the computation here. Instead, we present only the final result for the backreaction. It is given by . (A.38)

JHEP02(2020)135
When µ 2 ξ and 1 ξ (note that this limit is consistent with the parametrization in eq. (A.28)), the backreaction in eq. (A.38) is approximately given by where the negative sign implies that the backreaction plays a role of drag force in the classical picture. Unlike what was discussed in [52], we find that the condition µ 1 is not necessary to derive the approximation in eq. (A.39). Since ξ 1, the condition µ 2 ξ that leads to the approximation in eq. (A.39) is always satisfied for µ 1. Although the approximation in eq. (A.39) apparently does not depend on the Hubble parameter H as µ, ξ ∝ 1/H, it should not be considered to be valid in the H → 0 limit. Since two expansion parameters in terms of H scale like 1/ξ ∝ H and µ 2 /ξ ∝ 1/H, the limit that leads to the approximation in eq. (A.39) is not consistent with the H → 0 limit.
For the purpose of illustration, we numerically evaluate the integrand of the backreaction in eq. (A.35), namely before the integration over the momentum, for the same set of ξ values as those in figure 8, ξ = 10, 100 (smaller than our benchmark values), and they are shown in the upper panel of figure 9 as a function of x = k/aH. While the x-dependence in d 3 k/a 3 is not included in the upper panel of figure 9, including it amounts to simulate the integrand of eq. (A.36). As is evident in the lower panel of figure 9, it is more pronounced in a larger x region though due to the multiplicative x. Similarly to the occupation number in figure 8, the backreaction shuts off around x ∼ 2ξ as is expected.

A.6 Fermion energy density
In this section, we estimate the fermion energy density using the previously obtained fermion occupation number. Unlike the case of the backreaction, we provide the full computation filling the gap in [52]. The total energy density summed over the fermion and its anti-particle is given by We can re-express ρ ψ (τ ) in terms of s r and d r that we know their solutions, (A. 41) In what follows, we will calculate each contribution in the above expression individually, providing the detail. The first term of the second line in eq. (A.41) is trivially evaluated, and it is given by  In the limit Λ µ, it is approximated to be The second term of the second line in eq. (A.41) is the real part of what we have already calculated in the backreaction (except the helicity r in the helicity sum):

JHEP02(2020)135
The third term of the second line in eq. (A.41) is the volume integration of |s r | 2 . Using the integral representation of the Whittaker function, we can the analytic expression for it.
where the contour C in the t integration separates the poles of Γ 1 2 + ia + t Γ 1 2 − ia + t from those of Γ − 1 2 − irb − t , and the contour C for the s integration separates the poles of Γ 1 2 − ia + s Γ 1 2 + ia + s from those of Γ − 1 2 + irb − s . Since the integrand goes to zero as [t] → ∞, we can take the clockwise contour C. Similarly for the contour C . While there exist infinitely many poles at t = n − 1/2 − irb (with n ≥ 0) enclosed by C and at s = n − 1/2 − irb (with n ≥ 0) enclosed by C , only few poles contribute as Λ → ∞ due to Λ 3−s−t . The poles in the contour C that contribute are at (A.47) We consider each of them individually and sum them up later.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.