Phenomenology of relaxion-Higgs mixing

We show that the relaxion generically stops its rolling at a point that breaks CP leading to relaxion-Higgs mixing. This opens the door to a variety of observational probes since the possible relaxion mass spans a broad range from sub-eV to the GeV scale. We derive constraints from current experiments (fifth force, astrophysical and cosmological probes, beam dump, flavour, LEP and LHC) and present projections from future experiments such as NA62, SHiP and PIXIE. We find that a large region of the parameter space is already under the experimental scrutiny. All the experimental constraints we derive are equally applicable for general Higgs portal models. In addition, we show that simple multiaxion (clockwork) UV completions suffer from a mild fine tuning problem, which increases with the number of sites. These results favour a cut-off scale lower than the existing theoretical bounds.


Introduction
Relaxion models offer a new perspective on the hierarchy problem [1]. The weak scale is obtained in a dynamical way as the Higgs mass depends on a time-dependent vacuum expectation value (VEV) of a scalar field, φ . This scalar evolves and eventually halts at a value rendering the effective Higgs mass much smaller than the cutoff. This is achieved due to the fact that the potential of φ consists of a backreaction term that is switched on once the Higgs mass square gets negative and electroweak symmetry breaking (EWSB) is induced. When compared with conventional models of naturalness, this class of models leads to a completely different phenomenology as there is no analog of top or gauge partners that can be discovered at colliders. Instead, as we will discuss in detail in this paper, for experimental verification of the relaxion mechanism over a broad mass range we need to go to the low energy, high precision frontier.
Let us first present a very brief review of the relaxion mechanism. In relaxion models the value of µ 2 , the mass squared term in the Higgs potential, changes during the course of inflation as it varies with the classical value of φ, which slowly rolls because of a potential, In these equations g is a coupling 1 and Λ is the scale where the Higgs quadratic divergence gets cut off. Note that the operator in Eq. (1.3) can be radiatively generated by closing the Higgs loop in the term gΛφ(H † H) in Eq. (1.1) and thus technical naturalness demands r 1/16π 2 . In canonical models the field φ slowly rolls down (during inflation) from some initial large field value φ > Λ/g, such that µ 2 is positive and the electroweak symmetry unbroken. It stops rolling shortly after the point φ c Λ/g, where µ 2 becomes negative, electroweak symmetry is broken and the Higgs gets a vacuum expectation value, v 2 (φ) = −µ 2 (φ)/λ. A crucial ingredient of the relaxation proposal is the feedback mechanism that triggers a backreaction potential once the Higgs gets a VEV, where 1 ≤ j ≤ 4 is an integer andĥ = (v + h)/ √ 2. 2 As φ continues rolling, |µ 2 (φ)| becomes larger, resulting in a monotonically increasing Higgs VEV and thus increasing the backreaction's amplitude. Eventually the barriers become large enough and the relaxion stops rolling at an arbitrary O(1) value of the phase φ 0 /f , , (1.5) 1 Note that the coupling g defined here is dimensionless and is related to gGKR, the one in ref. [1], via gGKR = gΛ. 2 For an alternative proposal where the rolling is stopped due to particle production see ref. [2]. where Note that for values of Λ 4 br (v(φ 0 )) rgΛ 3 f , Eq. (1.5) is satisfied for | sin(φ 0 /f )| 1, but asM 4−j v j grows monotonically and the rolling starts at a random phase, φ/f , the relaxion stops well before it reaches this stage. Therefore, generically the phase φ 0 /f has an O(1) value. It is basically the result of a balance between the two terms controlling the derivative of the potential in Eq. (1.5). We shall return to this point when discussing CP-violation. For a small enough value of g, the cut-off can be raised much above the electroweak scale. Such a small value of g can be radiatively stable as in the g → 0 limit we recover the discrete symmetry, (1.7) Higher dimensional terms such as g 2 φ 2 Λ 2 , g 3 φ 3 Λ . . . contribute at the same order for φ φ c Λ g and thus do not affect the above analysis. For the same reason, in variants of the above model where only even powers of φ appear one can proceed along the same lines to obtain essentially the same results.
It is important to emphasize that in order to achieve higher values of the cut-off Λ, higher values of Λ br are required. Let us review the three main reasons for this. First of all, cosmological considerations during inflation (classical rolling must dominate over quantum fluctuations, see ref. [1]) put an upper bound on the cut-off which decreases if Λ br is smaller, where from now onward, for simplicity, by v,M and Λ br we will refer to the final values of these quantities at the relaxion minima φ = φ 0 . An even stronger bound can be derived if we demand that the relaxion does not have transplanckian excursions, which once again favours a large Λ br . As the requirement of subplanckian field excursions depends on quantum gravity assumptions and can be possibly evaded by UV model building, we will not take this as a strict bound and extend our analysis also to the transplanckian region. Finally, as argued in ref. [3], if the relaxion is a compact field, as any pseudo Nambu-Goldstone boson (PNGB) must be, the ratio of the distance the relaxion rolls to the periodicity of the backreaction, is, for a single axion sector, generically expected to be an O(1) number. On the other hand, this ratio must be large to raise the cut-off substantially above the weak scale, as Eq. (1.5) implies, Λ n r 1/4 Λ br . (1.11) Thus smaller values of Λ br require larger values of n for a given cut-off scale.
In this work we derive several new results relevant for relaxion phenomenology. We emphasise the importance of relaxion-Higgs mixing that is expected in a large class of relaxion models and focus on its experimental and observational implications. As the relaxion-Higgs mixing turns out to be proportional to Λ 4 br , these observational constraints put an upper bound on the backreaction scale, Λ br . We also derive a theoretical upper bound on the backreaction scale Λ br . We further consider multiaxion (clockwork) models where n ∼ e N , N being the number of sites in these models, and show that for a too large value of N , the clockwork construction becomes tuned unless further structure is assumed. By eqs. (1.8)-(1.11) we see that together these considerations favour lower values of the cut-off scale.
In the following section we derive the expressions for the relaxion-Higgs mixing and in section 3 review existing backreaction models and bounds on the backreaction scale. In section 4 we consider bounds on models with compact relaxions. We find a rich variety of experimental and observational probes for the relaxion in the mass range 0.1 µeV to 50 GeV described in detail in section 5 and section 6. All our bounds are equally applicable to general Higgs portal models. As the relaxion couplings to SM particles via the mixing are like that of a CP-even scalar, in the sub-eV range fifth force experiments can constrain large parts of the relaxion parameter space. In the keV-MeV range constraints on the relaxion parameter space arise from astrophysical star cooling bounds and cosmological probes of late decays, including constraints from entropy injection, BBN observables, CMB distortions and distortion of the extragalactic background light (EBL) spectrum. In the MeV-GeV region we find that the most important bounds arise from cosmological entropy injection and BBN bounds, cooling rate of the SN 1987 supernova, beam dump experiments and from constraints on rare B-and K-meson decays. Finally for GeV scale masses the bounds arise from LEP Higgs-strahlung data and LHC Higgs coupling bounds on the h → φφ channel. We also discuss how presently unconstrained parts of the relaxion parameter space would be probed by future data from experiments such as the PIXIE detector for CMB distortions, the NA62 experiment and especially the SHiP beam dump experiment. In section 7 we discuss the implications of our bounds on the theoretical parameter space of relaxion models. In section 8 we briefly discuss how the characteristic CP violation of relaxion models can be probed and finally we conclude in section 9. Useful relations are derived in the appendices.

Relaxion-Higgs mixing
Relaxion models contain two sources of breaking of the shift symmetry, the one that allows the Higgs mass to scan as the relaxion rolls and the backreaction term. In this section we will see how the presence of both terms can lead to spontaneous CP-violation in the backreaction sector and a measurable relaxion-Higgs mixing (see also ref. [4]). The full relaxion potential is given by combining the terms in Eq. To obtain the mixing terms we expand around the minima of the fields φ and H, In models with even j, v H = v = 246 GeV. On the other hand, as we will see in section 3, the backreaction sector breaks electroweak symmetry in models with odd j. In these models, where v is the electroweak symmetry breaking (EWSB) VEV in the backreaction sector. The minimisation conditions and explicit mass matrices, M 2 ij , for the j = 1 and 2 cases can be, respectively, found in appendix A and B. We find in both cases that the leading contribution to the mass matrix elements can be written entirely in terms of the parameters of the backreaction sector. In particular, for both j = 1 and j = 2. In addition, as discussed below, we expect Λ br v H which implies that the relaxion-Higgs mixing angles is naturally small, sin θ 1. We find that, to leading order the relaxion Higgs-mixing angle θ and the relaxion mass m φ are, As anticipated, the mixing angle is proportional to the spontaneous CP-violating spurion sin(φ 0 /f ) in the backreacting sector. In more complicated relaxion models there can be mechanisms to suppress this phase. An example is the model with the QCD axion where a small phase is necessary to be compatible with the non-observation of a strong CP-phase. In such models (see also ref. [5] and ref. [6]), relaxion-Higgs mixing is also suppressed.
Couplings: As the relaxion mixes with the Higgs boson, it inherits its couplings to SM particles suppressed by a the mixing angle sin θ as a universal factor -such as in Higgs portal models. 3 For g φψ , the coupling to pairs of matter fields ψ, and g φV , the coupling to pairs of V = W ± or Z, the couplings are given by At the loop level, the relaxion couples via quark loops to gluons and quark and W ± loops to photons, (2.8) 3 Note that in j = 1 models the Higgs couplings themselves might differ from their SM values because of the reduced Higgs VEV, vH = √ v 2 − v 2 ; in the following we will assume that v 100 GeV so that these are at most 10 % effects which we would ignore (see also section 3). where where τ x = m 2 h /4m 2 x . Let us finally comment on the pseudoscalar couplings of the relaxion to Standard Model particles, as these may have a significant impact on the experimental probes discussed in the following. However, these couplings are model-dependent and as the relaxion potential can be controlled by a sequestered sector [1,3] these couplings could be in principle suppressed relative to the "Higgs-portal" couplings discussed above (which are at the core of the relaxion construction). As we show in appendix C, this is the case in existing backreaction models (see section 3) where we find that these couplings are in magnitude generally smaller than or equal to the Higgs portal couplings. An exception is the pseudoscalar coupling to photons which in some backreaction models (see appendix C) can be larger than the one induced via Higgs mixing while in other models is of the same size as the Higgs-portal coupling. As the presence of a large pseudoscalar coupling to photons is thus modeldependent, we will comment on its implications only qualitatively.

Review of backreaction models and existing bounds on Λ br
As both the cutoff of the theory as well as the relaxion-Higgs mixing depends polynomially on the back reaction scale, it is important to examine what is its allowed range. In this section we thus describe the different backreaction models in the literature and discuss various bounds on the size of the scaleM or Λ br that appears in the backreaction potential [see Eq. (1.6)].
Note, first of all, that for odd j = 1 or 3, a non-zeroM in Eq. (1.4) must break electroweak symmetry which already suggestsM v, but let us analyze this case in more detail. The simplest relaxion model [1] where the backreacting sector is QCD and the relaxion couples to gluons like the axion, φ f G µνG µν is an example of a j = 1 model.
Non-perturbative effects generate a potential for the axion, where Λ 4 br = 4πf 3 π m u = 4πf 3 π y u v/ √ 2 is set by the pion decay constant f π and the up quark mass. As the relaxion stops at a generic value of the phase, φ/f , QCD relaxion models are generally in conflict with the non-observation of a large value of the strong CP phase. This problem can be solved in more complicated variants where there is a dynamical mechanism to make the above phase small.
An alternative approach would be to give-up the solution to the strong CP problem and to consider an additional strong sector. For instance, a new technicolor-like strong sector would lead to an EWSB condensate of techniquarks, Ū L U R +D L D R v 3 , where U L,R and D L,R are quarks with the same electroweak charges as the SM quarks u L,R and d L,R , but charged under the new strong group and not QCD. If the relaxion is coupled to the operator G µνG µν , involving the strong sector gauge bosons (G µν corresponds to the new strong sector field strength), a backreaction is generated with j = 1 and where y is the smaller of the U or D Yukawa coupling with the SM Higgs, and v H is the VEV of the Higgs doublet so that It is worth pointing out that such models would not be as strongly constrained as typical technicolor models because the condensate does not need to explain the large top mass and because the presence of an elementary Higgs somewhat alleviates the tension with electroweak precision observables [7]. In this work we have assumed v 100 GeV and ignored O(v 2 /v 2 H ) effects. A less constrained model with j = 2 was presented in ref. [1]. In this model, φ couples to G µνG µν , the gauge bosons of an EW symmetry preserving strong sector. The Higgs couples to two vector-like leptons charged under this strong group as follows, where (L, N ) have the same quantum numbers as the SM lepton doublet and right-handed neutrino, respectively, and (L c , N c ) are in the conjugate representations. If we take m N 4πf π m L , only the fermion N forms a condensate that is EW preserving. Upon integrating out L, L c , the Higgs contributes to the mass of N as follows, ∆m N = y 1 y 2ĥ 2 /m L so that the relaxion potential gets the backreaction A perturbative j = 2 model was presented in ref. [3]. In this model the relaxion is a familon, the PNGB of a spontaneously broken flavour symmetry. Let us consider the Lagrangian, where L and L c have the same quantum numbers as before, and N is a SM singlet fermion. The one-loop Coleman-Weinberg potential of the relaxion φ reads wherem is the larger of m L , m N . A theoretical challenge that any j = 2 model faces is that at the quantum level the backreaction term [see Eq. (1.4)] generates the termM 2 Λ 2 c 16π 2 cos φ f upon closing theĥ loop. This term is independent of the Higgs VEV, which implies the presence of an oscillatory potential even before the Higgs condenses [5]. Thus, the relaxion stops rolling prematurely, before EWSB, unless the scale Λ c at which the Higgs loop is cut-off satisfies Λ c 4πv . (3.8) In axion-like models this is automatically satisfied because the instanton contribution are highly suppressed at energy scales larger than the confinement scale, 4πf π , so that Eq. (3.8) implies f π v. In the model of Eq. (3.4) there is actually another contribution to the potential that exists even before EWSB, ∆V N 4πf 3 π m N cos φ f , where technical naturalness requires that m N must be larger than y 1 y 2 m L log(Λ/m L )/(16π 2 ) . Demanding the above EW preserving contribution to be smaller than the backreaction generated upon EWSB, 4πf 3 π ∆m N cos φ f , we obtain m L 4πv/ log Λ/m L . Together with these bounds, the requirement ∆m N 4πf π m L so that N forms a condensate and L does not, implies that in this model where we have assumed m L v due to experimental bounds for the first inequality. In the perturbative familon case, of Eq. (3.6), a simple extension can ensure that the constraint in Eq. (3.8) is satisfied as followed. The Majorana mass m N is actually induced via a mini see-saw mechanism. A new heavier fermion N c is added to the theory, After N c is integrated out, the Majorana mass of N is induced, m N = m 2 D /m N c . One can show in this case that two-loop corrections to the relaxion potential do not get contributions from energies above the scale m N c so that we get Λ c = m N c , and Eq. where we have assumed y 1,2 < 4π. Now let us discuss some model-independent bounds on the backreaction scale. First note that Eq. (2.5) implies that for a non-tachyonic φ, Finally notice that in the presence of the mixing the Higgs-like eigenvalue would satisfy, m 2 h > M 2 h h . For the j = 1 case this leads to a bound onM simply arising from the expression for the Higgs mass that is given by (see Eq. (A.8)) where the inequality becomes an equality in the limit of no relaxion-Higgs mixing. We must have λ > 0 to ensure that the potential does not have a runaway direction which implies the following bound, (3.14)

New bounds on compact relaxions
In this section we consider simple multiaxion (clockwork) models and then show that these suffer from stability issues when the number of sites (axions), N , becomes too large. The instability is, in fact, related to the very same issue of highly irrelevant operators that plagues the two-site construction in ref. [3]. First of all note that in realistic relaxion models the coupling g in Eq. (1.2) and Eq. (1.3) is obtained from a compact term (at least in QFT constructions where the relaxion is a pNGB), but with a larger periodicity F [3], which allows us to make the following identifications for κ < 1: One can now directly obtain Eq. (1.11) by demanding V (φ) = 0 using Eq. (4.1), As shown in ref. [8] and [9], the Choi-Kim-Yun (CKY) alignment mechanism [10] 4 (also known as the clockwork mechanism in the relaxion context) for multiple axions (or PNGBs) can provide a relaxion potential having two periodicities with a large ratio F/f ∼ e N , N being the number of axions. Let us first review these multiaxion models. We describe here the realization of ref. [9]. Consider N + 1 complex scalar fields Φ 1 to Φ N +1 with the potential The above potential respects a U(1) symmetry under which the fields Φ 1 , Φ 2 . . . Φ N have charges Q = 1, 1 3 . . . , 1 3 N . For simplicity the symmetry preserving cross terms such as Φ † i Φ i Φ † j Φ j have been ignored and an approximate permutation symmetry has been assumed (for → 0) so that the masses m 2 and quartic couplings λ are equal for all the fields. For λ the radial parts of the fields obtain a VEV, Φ i =f √ 2 e iφ i /f wheref 2 = 4m 2 /λ, such that at low energies only the angular degree of freedom remains. N superpositions of the angular fields obtain masses, but the direction is a flat direction that describes a Goldstone boson. The Goldstone mode has an O(1) overlap with the first site and is exponentially suppressed overlap with the last site, is the norm of the vector defined by Eq. (4.5). Let us now introduce some anomalous breaking of the global U(1) at the first and last sites, where f = Nf and we have used eqs. (4.6), (4.7) to rewrite the first line in terms of φ 0 . Non-peturbative effects now generate the desired relaxion potential in Eq. (4.1) with F = 3 N f so that Eq. (4.3) now becomes Thus we see that the CKY/clockwork mechanism can give us a cut-off that grows exponentially with the number of axions, N . Note that the above analysis holds only if so that the potential generated by the anomalous breaking of the U(1) symmetry in Eq. (4.1) is subdominant compared to the potential generated from Eq. (4.4), the linear combination in Eq. (4.5) remains a Goldstone mode, and all the heavier modes can be decoupled.
We now show that, for a too large N , these models become finely tuned if we relax the approximate permutation symmetry in Eq. (4.4) that was assumed only for convenience of calculation. If we relax this assumption some of the mass square terms might be positive. Let us assume, for instance, that k − 1 consecutive fields, Φ n 1 +1 . . . Φ n 1 +k−1 have positive mass square terms so that there are no corresponding PNGB modes φ n 1 +1 . . . φ n 1 +k−1 for these scalars. At first sight, this breaks the link in the axion chain because -instead of one Goldstone mode as in Eq. (4.5) -there are now two decoupled Goldstones fields [in the absence of the subdominant boundary terms of Eq. (4.8)], where n 2 = N − n 1 − k and N 1 and N 2 are again normalisation constants. None of the above modes can be identified with the relaxion as no single mode above is subject to both the backreaction at the first site and the rolling potential at the last site. However, the link between the two chains is not completely lost, as a process like Φ n 1 −1 → 3Φ n 1 → . . . → 3 m Φ n 1 +k generates a higher dimensional operator that weakly couples the two sectors, is an exponentially small number due to /λ 1. More precisely, as the N − 1 heavier pseudo-Goldstone modes that have masses m 2 φ ∼ 3 2 f 2 /2 (see ref. [9]), must be much lighter than the radial modes ρ with m 2 ρ ∼ λf 2 /2, one needs (4.14) With the above term in the potential we once again obtain that φ 0 from Eq. (4.5) is a Goldstone mode. With the addition of the explicit breaking terms on the first and last site in Eq. (4.8), the terms relevant for the potential of φ 01 and φ 02 are, where we have appropriately redefined φ 01 and φ 02 so that the phase appears only in the last term. The two lightest modes now are superpositions of Eq. (4.11) and Eq. (4.12). The mass matrix of φ 01 and φ 02 is given by which results in, up to normalization factors, the two mass eigenstates, where s α = sin α, c α = cos α and α is the mixing angle. Let us first show that in the limit that contribution of the term proportional to Λ 4 to the gradient of the φ m2 potential is subdominant, i.e. Λ 4 3 n 2 +k ε f 4 , we recover the usual relaxion potential. In this limit we obtain tan α = 3 −n 1 −k N 2 /N 1 , and the first eigenstate in Eq. (4.17) becomes identical to the relaxion mode in Eq. (4.5). To obtain the Lagrangian for the lightest mode we first use the condition to stabilise φ m2 , which in this limit reads Substituting the solution φ m2 = 0 in Eq. (4.15) and using tan α = 3 −n 1 −k N 2 /N 1 yields the potential in Eq. (4.1).
In the opposite limit, i.e. Λ 4 3 n 2 +k ε f 4 , the gradient of the φ m2 potential is dominated by the term proportional to Λ 4 , which drives φ m2 to the global minimum of the rolling potential giving the Higgs an O(Λ 2 ) mass. Therefore for the relaxion mechanism to work one needs Λ 4 3 n 2 +k ε f 4 , which implies We see that for k = 4, λ = 1, r = 1/16π 2 and N = 28, using f < M Pl we get Λ 2 TeV for any positive integer n 1 , so that the relaxion mechanism cannot even address the little hierarchy problem in this case.
How long must the relaxion chain be so that a sequence of k − 1 = 3 consecutive positive masses becomes highly probable? To compute this probability we need to find the number, N 3 (N ), of sequences of N '+' or '-' signs with at least one chain of 3 consecutive positive signs '+++ ' . It can be shown that N 3 (N ) obeys the following recursion relation, Here the first term comes from the fact that if we already have at least one '+++' chain in a sequence of N axions, by adding either a '+' or '-' sign at the N + 1 th position we obtain an arrangement of size N + 1 satisfying our criterion. This does not include arrangements of N axions with no chain of 3 consecutive positive '+' signs but having a '-++' at the end such that we get a required arrangement if at the N + 1 th position we add a '+' sign; this is taken care of by the second term in Eq. (4.22). Finally in the last term we subtract the double counting resulting in cases where the sequence captured by the second term already includes a '+++' in the remaining subchain.
To obtain a successful relaxion model, however, we are interested in an N -site sequence with no '+++' chain, N 3 (N ), which is given by (4.23) It turns out that N 3 (N ) satisfies the following familiar relation, which is nothing but the recurrence relation of the 3-step Fibonacci sequence. 5 By inspection, N 3 (3) = 7 , the 5th element of the 3-step Fibonacci sequence so that we must have N 3 (N ) = fib 3 (N + 2). Our arguments can be easily generalized to find the number of arrangements with no chains of k − 1 positive masses which turns out to be just the (k − 1)-step Fibonacci sequence. Hence the probability to randomly obtain a sequence with at least one chain of k − 1 positive masses is We find that for N ≥ 28 the probability of having at least k − 1 = 3 consecutive positive masses in a chain of N axions is P(k − 1, N ) 90%. Thus for for N ≥ 28, from Eq. (4.20) generically we have Λ 3 TeV as already discussed above. For N 28 axions there is the possibility of raising the cut-off to a value of where we have used Eq. (4.9) and the numerical value above is for Λ br √ m h v. As we will show in sections 5 and 6, experimental probes can constrain Λ br to even smaller values as a function of f (or alternatively the relaxion mass) and this in turn would imply an even lower cut-off in accordance with Eq. (4.26).

Laboratory probes of relaxion-Higgs mixing
In this and the next section we discuss in detail the bounds and the future probes for relaxion-Higgs mixing, distinguishing between laboratory experiments, discussed here, and cosmological and astrophysical probes considered in section 6. As we will show below, the relaxion mass can range from far below the eV-scale to almost the weak scale. Therefore a variety of experiments is needed to look for the relaxion. As the couplings to SM particles are proportional to sin θ, a convenient plane to present the constraints is the sin 2 θ-m φ The lifetime τ φ also depending on m φ and sin 2 θ with thresholds (vertical gray lines) and example values of τ φ (horizontal gray lines). The lifetime for any other sin 2 θ value can be obtained from the sin 2 θ = 1 line using τ φ ∼ 1/ sin 2 θ.
plane. Before going into the details of the various constraints to be presented in figures 2-5, let us first identify the regions of the sin 2 θ-m φ plane that are relevant for relaxion models. For the convenience of the reader we repeat the expressions for the mass and mixing angle of the relaxion in the small-mixing approximation from Eq. (2.4) and Eq. (2.5) substituting for definiteness sin where (Λ max br ) 2 = 2 1/4 m h v/j is the maximal allowed value of Λ br that follows from Eq. (3.12). 7 Other O(1) choices of φ 0 /f lead to slightly modified numerical values for 6 Note that we have assumed v 2 H v 2 which amounts to ignoring, at most, O(10%) effects (see section 3). 7 It was shown in ref. [12] that in j = 2 models it is possible to have Λ br m h v with smaller then O(1) values for sin(φ0/f ) v 2 /Λ 2 br , such that Eq. (3.12) is still satisfied. In this work we take O(1) values of sin(φ0/f ), and in accordance with Eq. (3.12), Λ br m h v thus not considering this region of the parameter space. As the backreaction scale, Λ br , is in any case constrained to be less than a few times the weak scale (see section 3), this is actually a narrow region of the parameter space where the constraints are expected to be similar to those we obtain for Λ 2 br ∼ m h v. m φ , sin θ and Λ max br . These can be inverted to obtain We use the equations above to make contours of constant Λ br and f in the sin 2 θ-m φ plane in figure 2, 3 and 5. Although we have made the contours for the j = 2 case, using Eq. (5.2) one can easily translate to the j = 1 case by substituting f → 2f, Λ br → √ 2Λ br . For relaxions heavier than 5 GeV, the mixing can be O(1) and Eq. (5.1) and Eq. (5.2) are no longer valid. Thus in section 5.2.1 and figure 4 where we consider relaxions in this mass range, we exactly diagonalize the mass matrices in appendix A and B to obtain the Λ br and f contours.
We see from Eq. (5.1) and that if Λ br is much smaller than Λ max br , for a given f , both m φ and sin 2 θ increase with Λ br and we get sin 2 θ ∼ m 4 φ . This implies that in this regime, a light relaxion has typically a suppressed mixing. However, if we take values of Λ br close to Λ max br , this tendency does not hold anymore. This behaviour can be seen in figure 1(a) where we plot the relaxion mass as a function of Λ br for different values of f taking j = 2. We see that, for all f , the relaxion mass is maximum for Λ br = Λ * br = 2 −1/4 m h v/j, and for larger values of Λ br the relaxion mass drops rapidly with Λ br as the term within the parenthesis in Eq. (5.1) becomes smaller. The relaxion mass can, in fact, be made arbitrarily small by choosing a Λ br that is sufficiently close to its maximal value Λ max br , while hardly changing sin 2 θ. In the sin 2 θ-m φ plane this can be seen from the shape of the f contours in figure 2 (and subsequent figures) for which two branches can be clearly identified. The region Λ br > Λ * br corresponds to the top left part of figure 2 where the f contours become nearly horizontal as sin 2 θ hardly changes but the mass can become arbitrarily small. The thick grey line in figure 2 is the contour Λ br = 0.99Λ max br . The whole region above the this line, which we refer to as the "tuned region", corresponds to the narrow region in the theory space 0.99Λ max br < Λ br < Λ max br , marked by the thick black line in figure 1(a). Therefore, in the following we will mostly discuss the "untuned region" Λ br < Λ max br , which translates to and implies that in most of the theoretical parameter space, if we make the relaxion lighter, it also becomes more weakly coupled. We would like to point out that this is a general feature of Higgs portal models. For instance, consider the potential [13], where φ and h are as defined in Eq. (2.2) whilem 2 φ ,m 2 h and x parametrise the couplings in a general Higgs-portal model. For small mixing angles we getm h = m h , the observed Higgs mass, and the condition x < x max = m h /v to ensure that the lighter eigenstate does not become tachyonic. The region above the grey line in this case corresponds to the small range 0.98 The second restriction on the sin 2 θ − m φ parameter space arises from the fact that we consider only the range, (5.5) The lower bound on f arises from the fact that in our analysis of relaxion-Higgs mixing we ignored any new states (for instance radial modes) that must exist below the scale Λ = 4πf to UV-complete the backreacting sector. Thus our analysis holds only if both the Higgs boson and the relaxion have a mass much smaller than the mass scale of these UV states, i.e. for f m h . We will call the region defined by Eq. (5.3) and Eq. (5.5) the 'relaxion parameter space', i.e the region in the sin 2 θ − m φ relevant for relaxion models.
Let us now discuss the mass range the relaxion can have given these restrictions. In the untuned region the relaxion can be made lighter either by decreasing Λ br or increasing f . In our analysis we do not consider f > M Pl = 2 · 10 18 GeV, but as there is no strict lower bound on Λ br , the relaxion can be made as light as we want by taking sufficiently small values of Λ br . As discussed in section 1, however, lower values of Λ br are theoretically disfavoured. For instance if we require relaxion field excursions to be subplanckian this puts a bound sin 2 θ 10 −27 as shown in figure 2. In the untuned region this can be translated to m φ > 0.001 eV. As the requirement of subplanckian field excursions depends on quantum gravity assumptions and can be possibly evaded by UV model building, we will not take this as a strict bound and extend our constraints also to the transplanckian region.
We now turn to the question of how heavy the relaxion can be. The largest relaxion mass is obtained for the minimal value, f = m h , and weak scale values of Λ br where the small-mixing approximation in Eq. (5.1) no longer holds. In section 5.2.1, by exactly diagonalising the mass matrices in appendix A and B, we find an upper bound m φ 60 GeV (see figure 4).
For readers interested in general Higgs portal models our analysis provides the complete constraints in the untuned region of their parameter space apart from the area outside the region defined by Eq. (5.5). Whereas for f > M Pl , the constraints in the untuned part of the region arise only from fifth force experiments and have been discussed elsewhere (see for instance [13,14]), the region corresponding to f m h can be potentially constrained only by some cosmological probes that we will mention in the next section but not fully derive.
Before going into the details of the different experimental probes, a comment is in order. In the following we are going to study the constraints on the relaxion parameter space driven by its mixing with the Higgs. As it is impossible to include the effects of the pseudoscalar couplings of the relaxion in a model-independent way we do not consider these. In any case, in existing explicit models, these couplings are generally not larger than the Higgs portal couplings as discussed in appendix C. An exception is the pseudoscalar coupling to photons which can in some backreaction models (see appendix C) be larger than the one induced via Higgs mixing. In section 7 we qualitatively comment on how our constraints would change if a large pseudoscalar coupling to photons is present.
In the following we describe the constraints on the relaxion in different mass ranges as the relaxion mass spans a wide range from sub-eV values to tens of GeV. While relaxions heavier than a MeV can be potentially probed by collider searches, the only laboratory probes for sub-MeV relaxions are fifth force experiments. We discuss these two categories separately starting with sub-MeV relaxions.

The sub-MeV mass range
In this mass range the relaxion has a very large decay length making it impossible for collider searches to probe visible decays of the relaxion. This can be seen from figure 1(b) where we plot, using the expressions in appendix E, the relaxion lifetime as a function of its mass for different choices of sin θ. Eq. (5.3) implies for the considered mass range sin θ 10 −9 , and figure 1(b) shows the corresponding enormous rest frame decay length of cτ 10 14 m. Therefore the only possible laboratory probes are either fifth force experiments, or experiments looking for invisible particles. This last class of experiments, at least at the moment, is not sensitive enough to provide constraints on the very small Higgs-relaxion mixing in this mass range [15]. Fifth force experiments denote experiments which can detect the existence of a new degree of freedom by the corresponding new Yukawa-like force induced between two electrically neutral test bodies. A relaxion induces a spin-independent Yukawa force between two test bodies A and B, defined by the potential where m A , m B are their respective masses and α A , α B parametrise the couplings of the relaxion to the two bodies. In Higgs portal models, the couplings are given by [13] where g hN N 10 −3 and m nuc = 1 GeV. The sensitivity of the various fifth force experiments depends on the interaction length λ which is related to the mediator mass m φ via Let us start discussing probes of new long-range forces going down from macroscopic length scales to the pm scale of MeV particles. We present the bounds arising from these probes in figure 2. For very low masses (below 3 · 10 −15 GeV), the strongest constraint comes from the Eöt-Wash experiments [16,17] that looked for deviations from Einstein's weak equivalence principle (labelled as EqP in figure 2) by precision measurements of the longrange force between a heavy attractor and two different test bodies in a torsion balance. Let us notice that this experiment is able to constrain the Higgs portal down to very small couplings, but for masses lighter than 10 −16 GeV the probed parameter space belongs to the tuned region (for other potentially relevant discussion of cosmological and/or low energy probes see for instance [18,19]). Therefore, in figure 2 we do not show relaxion masses of m φ < 10 −16 GeV although the EqP bound extends even further. On shorter length scales, the mass range 3 · 10 −15 -10 −11 GeV, the strongest bounds arise from constraints on violations of the inverse square law (labelled as InvSqL) that have been obtained by various experimental groups [20][21][22][23][24][25]. The excluded region shown in figure 2 is an envelope that contains bounds from all these experiments with the strongest one coming from the Irvine experiment in the mass range 3 · 10 −15 -5 · 10 −14 GeV [20,21], from the Eöt-Wash 2006 experiments in the mass range 5 · 10 −14 -2 · 10 −12 GeV [25] and from the Stanford experiment [22,24] in the mass range 2 · 10 −12 -5 · 10 −11 GeV. Finally we also show the constraints from tests of the Casimir force [26,27], the force induced by the zero point energy of the electromagnetic field when two conductors are brought very close to each other. While these bounds from the tests of the Casimir effect are weaker than the bounds of the torsion balance experiments below 10 −11 GeV, they are the strongest bounds above this mass as shown in figure 2. The shaded area below the horizontal light gray, dotted line (sin 2 θ ≤ 10 −27 ) shows the region where the relaxion has transplanckian excursions for any value of the cut-off scale Λ > 2 TeV (see Eq. (1.9)). For heavier particles, i.e. shorter-range forces, the sensitivity is even lower. The intermediate region, between 10 eV and 1 MeV, is the most challenging region to probe in laboratories. The most sensitive experiment in this mass region are neutron scattering experiments that test the existence of a new sub-MeV boson based on their influence on the neutron-nucleus interaction. These experiments set a very weak bound, s θ 0.1 [28,29], and are therefore incapable of probing a relevant region of the parameter space. In a subset of this mass range from a keV to an MeV (shown later in figure 5) the relaxion parameter space can be probed only by astrophysical and cosmololgical observations to be discussed in detail in the next section. The 10 eV-keV mass range, on the other hand, is largely unconstrained as shown in figure 2.
Let us conclude this subsection by commenting that fifth force experiments are a unique probe of light states like relaxions that couple to electrons and nucleons as CP-even scalars. Axions, for instance, do not give rise to spin-independent long range forces at leading order because of their pseudoscalar nature and are thus only weakly constrained by fifth force experiments. Therefore, different laboratory probes have been proposed to circumvent this problem. This is the case for light shining through the wall (LSW) experiments [30], which are also sensitive to Higgs portal models [31]. However, their reach is too limited to compete with fifth force experiments and therefore these do not appear in our plot.

Relaxion masses between the MeV-and the weak scale
Let us now study the region of parameter space where the relaxion mass is above the electron threshold and thus it can decay into SM fermions. Furthermore, as shown in figure 1(b), in this region the relaxion has a shorter lifetime and can be directly searched for in laboratory facilities. Let us further distinguish two sub-regions based on the different relevant probes. The bounds in the MeV-5 GeV mass range are presented in figure 3, including also astrophysical and cosmological constraints which will be discussed in the next section. Figure 4 presents the bounds in the GeV region.

The 1 MeV-5 GeV range
This region of the parameter space is well covered by rare K-and B-meson decays at proton beam dump and flavour experiments. Crucial for both kinds of experiments is the possibility of producing a relaxion in rare decays of K-and B-mesons. In flavour experiments that probe rare decays, constraints are put on the branching ratios [32] where p φ is found using two-body kinematics and F K is defined in [32]. Even in proton beam dump experiments, rare mesons decays are the main the production mode of the relaxion. The smallness of the branching ratio is overcome by the large luminosity. Electron beam dump experiments do not have any sensitivity to Higgs-relaxion mixing due to the suppressed electron Yukawa coupling.
Beam dump experiments: In proton fixed target experiments, relaxion beams are produced from meson decays and constraints are imposed by looking at its visible decays. The region to the left of the cτ = 2 m line in figure 3 is roughly the region where the relaxion decay length in the lab frame is greater than about 100 m (assuming a relativistic boost factor γβ ∼ 50 [32]). Thus this is the region relevant for beam dump experiments looking for long lived particles. We will discuss here the sensitivity of the CHARM experiment and future experiments such as SHIP [33] and SeaQuest [34]. NA62 is also planning a beam dump run as proposed in ref. [35].
The CHARM beam dump experiment performed a search for long-lived axion-like particles decaying to e + e − , µ + µ − or γγ in collisions of a 400 GeV proton beam on a copper target [36] with a 35 m long detector located 480 m from the target. This search can be also reinterpreted in the context of Higgs portal models [32,33,37,38], where the scalar φ is predominantly produced in rare decays of K-and B-mesons. In such an experiment around 10 17 kaons and 10 10 B-mesons are produced per year. Figure 3 shows that CHARM (dark red) is able to constrain only masses below the kaon threshold. The limited sensitivity is due to the lower B-meson luminosity and the large distance between target and detector. The large distance, on the other hand, is good to probe the low mass region where the relaxion has a longer decay length as shown in figure 1(b). In figure 3 we show also the projections (in lighter red) for several future proton beam dump experiments such as SHIP (dotted) and SeaQuest (dash-dotted). While the reach of NuCal exceeds CHARM for a scalar/pseudoscalar with couplings only to photons [35], the NuCal bound in the presence of Yukawa-like couplings to fermions is in ref. [39] found to be weaker than the CHARM limit and therefore we omit NuCal in figure 1(b). If in future beam dump experiments the detector is closer to the target than in the case of CHARM, good improvements over CHARM can be achieved for relaxion masses heavier than the muon threshold where the lifetime is shorter. As already noticed, lighter relaxions have a longer lifetime and therefore the CHARM bound in this region can be improved at proton fixed target experiments by looking for invisible new particles. The present sensitivity, however, is limited in the region of the parameter space relevant to our scenario [40][41][42][43].
Rare meson decays: Rare decays of K-, B-and Υ-mesons can be mediated by a light scalar particle φ, and therefore bounds on their branching ratios constrain the relaxion-Higgs mixing angle. In figure 3 the turquoise region corresponds to the bounds on B-decays and the blue region to K-decays. We do not show bounds coming from rare Υ decays since they are always weaker than other existing bounds.
Let us first discuss how B-decays constrain the relaxion-Higgs mixing. Both Belle and LHCb are sensitive to the decay process of B ± → K ± φ → K ± l + l − with l = µ at LHCb and l = µ, e at Belle [44,45]. In the experimental analyses, the regions of m ll in [2.95 GeV, 3.18 GeV] and in [3.59, 3.77 GeV] are vetoed in order to suppress the background from the J/ψ and the ψ resonances, respectively. Figure 3 shows the constraints on sin 2 θ derived in [38] using the upper bound on the branching ratio as a function of the dilepton invariant mass provided by LHCb and Belle (both in turquoise). The figure indicates an almost comparable sensitivity of both experiments in the mass region above ∼ 300 MeV. 8 In addition, LHCb constrains BR(B 0 → K 0 * φ) BR(φ → µ + µ − ) as a function of the mass and the lifetime of a new boson φ. Figure 3 includes this LHCb search as a bound on (m φ , sin 2 θ) for m φ ≤ 1 GeV as presented in ref. [46] for a scalar mixing with the Higgs from the model of ref. [47], which also applies to the relaxion case. It improves the previous LHCb bound by appropriately an order of magnitude. Using the full 2-dimensional information provided in ref. [46], an extension of the bound up to m φ ≤ 4.35 GeV would be possible.
In the range of 0.212 GeV ≤ m φ ≤ 0.3 GeV [48], Belle has performed a dedicated study of B 0 → K * 0 µ + µ − , searching for a peak in the dimuon spectrum in order to enhance the sensitivity; this bound is also shown in figure 3. For even lighter masses, the limit on the B → K + invisible from Belle and BaBar (also indicated in turquoise) constrains the relaxion parameter space in a region where the relaxion has a long decay length. In both of the above cases we have used the constraints derived by ref. [32].
Let us now discuss the constraints set by the searches for visible and invisible rare Kmeson decays, using again the results of ref. [32]. These are shown in dark blue in figure 3. A search for K L → π 0 l + l − has been performed at KTeV/E799 [49,50] and translated into bounds on a pseudoscalar [39] and scalar [33] mediator of this decay. The corresponding constraint (dark blue) in figure 3 is stronger above the muon threshold, where it surpasses the current constraints from B-decays, and much weaker when the only visible decay mode is into electrons (shown as a dark blue line). The branching ratio of K ± → π ± µ + µ − has been measured by the NA48/2 fixed target experiment [51] at the CERN SPS. Despite the good agreement of the branching ratio with the SM prediction, the resulting bound on sin 2 θ is weaker than those derived from the B-decays into visible final states due to the relative CKM suppression of V ts · V td compared to V tb · V ts in the W − t -loop of the penguin diagram.
On the other hand, strong constraints can be set looking at invisible K-decays and this is a promising search for light relaxions. Indeed for small enough couplings and/or light enough masses -more precisely the region to the left of the cτ = 2 m contour in figure 3 -the relaxion decays outside the detector (see also figure 1(b)). Searches for the invisible K-decay K ± → π ± + invisible have been performed by the E787 [52] and the E949 [53] experiments at BNL, also considering two-body decays and providing limits on BR(K → π + φ) BR(φ → invisible) as a function of m φ . The constraints on the Higgs portal model coming from these analyses were previously studied in [32], which, however, focussed only on the region above 100 MeV, while the search is sensitive also to lighter relaxions. Therefore, we extended the analysis to lower masses as shown in figure 3; the gap in the mass range 0.1 GeV ≤ m φ ≤ 0.21 GeV is due to the fact that the region around m π ± has been vetoed. We find that the E949 experiment gives stronger limits than the CHARM beam dump experiment for relaxions with m φ ≤ 10 MeV. Indeed, lighter relaxions have a larger decay length, thus they are most likely detected as invisible particles. Therefore, this is one of the most promising regions for rare K-decay measurements to probe new physics. For instance the CERN experiment NA62 will improve the present limit on invisible K- decays by almost an order of magnitude. They expect to see 90 SM signal events and 20 background events in two years [54]. Using only this information about the total rate and no information about the differential distribution of the SM and background events, we show a conservative estimate of the 95% CL excluded region in light blue in figure 3 where we have assumed a 10% theoretical error [55]. The gap in the excluded region is again due to the veto around the charged pion mass, 100 MeV m φ 160 MeV [54].
Finally, for GeV-scale masses we see from figure 3 that some regions of the parameter space are bounded by LEP and LHC searches that we describe in detail in the next section.

The m φ > 5 GeV mass range
Finally we consider the mass region m φ > 5 GeV where the mixing angle sin θ can become O(1) and the expressions in Eq. (5.1) do not apply anymore. To compute the mixing angle, sin θ, and the mass, m φ , as functions of Λ br and f , we therefore exactly diagonalise the mass matrix in appendices A and B for the j = 2 (j = 1) case. We fix the value of the unknown λ by demanding that we obtain the observed Higgs mass for the heavier eigenvalue. This is how we obtain Λ br and f contours in figure 4. It is in this region that we obtain lowest values of f close to m h . As discussed in the beginning of this section, for even smaller values of f < m h our analysis of relaxion-Higgs mixing does not hold anymore.
LEP constraints: In the high-mass range, LEP and the LHC provide useful constraints on the mass and coupling of the relaxion. At LEP, the Higgs-strahlung process of, e + e − → Z → Z * h with Z * and h each decaying to a pair of fermions, is sensitive to the Higgsrelaxion mixing. If m φ < 2m µ , φ escapes the detector. For visible decays above the dimuon threshold, L3 [56] sets the most stringent bounds in the range m φ < 11.5 GeV whereas for 12 GeV ≤ m φ ≤ 116 GeV the combination of the four experiments ALEPH, DELPHI, L3 and OPAL at LEP [57] constrains this process most strongly. The experiments provide a mass-dependent upper bound on the ratio of cross sections [57], where σ max is the largest cross section σ(e + e − → Z → Z * φ) compatible within the 95% CL with the combined data sets, and σ SM is the SM reference cross section σ(e + e − → Z → Z * H SM ). In Higgs portal models, the ratio of φZ-production to the SM Higgs production cross-section for the same mass is just sin 2 θ, so that S 95 can be directly interpreted as the 95% CL upper bound on sin 2 θ. We show the parameter space excluded by LEP in green, labelled by "LEP hZ", in figure 4.
Higgs coupling bound on h → φφ: Finally we discuss how Higgs coupling measurements at the LHC constrain the h → φφ process. The strongest constraint on the partial width to this non-standard decay channel arises from the potential dilution it can cause to the visible decay channels of the Higgs boson to SM particles. While such a dilution of the visible decay channels may be compensated by increased scaling factors of the couplings [58], this is not the case in Higgs portal models (like the relaxion case we are considering) where the Higgs boson couplings are universally suppressed by cos θ with respect to their SM values. This configuration with one universal coupling modifier and non-standard decay channels has been considered in ref. [58]. Therefore we apply their upper limit on the Higgs branching ratio to non-standard channels from a fit to the data of ATLAS and CMS at 8 TeV with HiggsSignals [59]: We compute the partial width of h into φφ, using the coupling g hφφ that has been derived exactly in Eq. (D.2) in appendix D for j = 2, taking h − φ mixing into account. The hφφ coupling is parametrically different in j = 1 models and has not been considered here (see appendix D) . This allows to set bounds on the relaxion parameter space via (5.14) where Γ SM h = 4.12 MeV [60].
Higgs decays to two relaxions at the LHC: In addition, the explicit searches at the LHC for non-standard decays of the Higgs boson (see e.g. ref. [61]) with a mass of m h = 125 GeV include the decay channel of the Higgs boson into two lighter scalars (or pseudoscalars) φ that each further decay into a pair of fermions f or photons γ: h → φφ → 4f /4γ at ATLAS [62,63] and CMS [64][65][66][67]. Their results can be interpreted as bounds on the decay of the Higgs boson into two relaxions that further decay into the analysed final states. So far, data is only available from LHC Run 1 at 8 TeV. In ref. [66], the CMS searches with the final states 4µ, 4τ, 2µ2τ and 2µ2b have been translated into upper bounds on at the 95% CL under the assumption of g φf ∝ m f , which holds also in the relaxion case, see section 2. Therefore, the prediction of R hφµ depending on m φ and sin 2 θ compared to the experimental limits provides bounds in the (m φ , sin 2 θ) plane for those values of m φ that are covered by the set of searches. The mass range of 0.25 GeV ≤ m φ ≤ 3.55 GeV is covered by the 4µ final state, and the current data is sufficient to exclude parts of the parameter space shown in figure 3 (blue), but not stronger than the flavour bounds in this region. At higher masses, the 2µ2b final state is particularly sensitive due to the enhanced branching ratio of φ → bb compared to φ → µµ. However, the ATLAS and CMS searches based on the data from Run 1 do not constrain the relaxion parameter space beyond the constraints derived from the Higgs couplings fit. Assuming an improvement of the experimental limits by a factor of 10, which we very roughly estimate (neglecting the change of systematic uncertainties and not combining channels and experiments) as the reach during Run 3, we also show projections of the bounds for Run 3 (dark blue, dotted). While the 4τ final state at the LHC will not set stronger bounds than Higgs-strahlung at LEP, the constraints coming from the 2µ2b channel might provide an improvement comparable to the projected Higgs couplings fits. To summarise, figure 4 visualises that the bounds from LEP and the LHC are complementary in the sense that LEP is more constraining on sin 2 θ for m φ <25 GeV whereas the indirect constraint from the bound on the decay width into NP final states at the LHC sets a stronger constraint for m φ >25 GeV. Again we show contours of constant Λ br and f which, as we already mentioned, have been obtained by exact diagonalisation of the mass matrices in appendices A and B. We show the contours for Λ br = 120 GeV for j = 2 (gray, dashed) and j = 1 (brown, dashed), f = m h and f = 1 TeV for both the j = 2 (black) and the j = 1 case (brown).

Cosmological and astrophysical probes of relaxion-Higgs mixing
As discussed in the previous section, laboratory measurements can probe a significant region of the relaxion parameter space. However, in the sub-MeV region, before the fifth force experiments start to gain sensitivity in the sub-eV region, a large portion of the parameter space is left unconstrained. In this section we show how astrophysical and cosmological probes can explore part of this region of the parameter space, as shown in figure 5, and also provide relevant bounds if the relaxion mass is in the MeV-GeV range (also shown in figure 3). In order to identify the part of the parameter space most relevant for relaxion models and to gain an understanding of the theory contours in figure 5, we refer the reader to the discussion at beginning of section 5.

Cosmological probes
Late relaxion decays can be constrained by a variety of cosmological probes such as light element abundances, CMB spectral distortions and distortions of the diffuse extragalactic background light (EBL) spectrum. In this section we first compute the relaxion abundance generated by misalignment and thermal production and then use this result to study how these bounds apply to our scenario.

Relaxion abundance
Misalignment production: During inflation the expectation value of the field φ, φ , satisfies the classical equation of motion. Quantum fluctuations lead to a spreading of the field around this classical value. The spreading is given by (see for instance ref. [5]) where ∆φ = φ− φ and ∆φ 2 = (φ− φ ) 2 , H I is the Hubble scale during inflation and N e is the number of e-folds. We see that the spreading stops when the r.h.s. above vanishes, that is for where V (φ) = ∂V /∂φ, and to obtain the inequality we have used the requirement that the dynamics of the relaxion is dominated by classical rolling and not quantum fluctuations, H I < (V (φ)) 1/3 (see ref. [1]). This gives us the misalignment of φ from its classical value just after inflation. After this, the Universe goes through a phase of radiation domination. If the temperature of the Universe is below the temperature T 0 with the relaxion field oscillates around the minimum. This leads to an energy density, ρ φm , and an effective non-relativistic number density, n φm , given by 5) and thus results in a comoving number density, where ∆φ max is that maximal value of ∆φ given by Eq. (6.2), the entropy density, s = 0.44 g S * (T i ) T 3 i and g S * (T i ) is the effective number of degrees of freedom in entropy at the temperature T i . If the reheating temperature is larger than T 0 , then T i = T 0 , otherwise T i is the reheating temperature.
Thermal production: Relaxions can be thermally produced by the process HH → φφ at temperatures above the Higgs mass, by the processes q(g)+g → q(g)+φ at temperatures below the electroweak critical temperature, T EW , by the pion-relaxion conversion process N +π → N +φ at temperatures below Λ QCD , and finally by inverse decays. Let us consider these processes one by one.
At temperatures above the electroweak critical temperature, T EW ∼ m h v/m t the Higgs portal mixing in Eq. (2.5) is absent and the relaxion interacts only with the Higgs doublet. The main production mode of the relaxion is then the process HH → φφ via the coupling Note that any contribution to the process from the backreaction potential is absent, because in both the non-perturbative axion and the peturbative familon model, the backreaction term dissolves at high temperatures. In the axion case, the potential becomes negligible at high temperatures because instanton effects become very weak as the non-abelian gauge coupling becomes perturbative. In the familon model the Coleman-Weinberg potential gets no contributions from momenta above m N c so that for T m N c the backreaction potential vanishes also in this case. The comoving number density for φ resulting from this process has been computed in ref. [5] to be where we have not considered any contribution above the electroweak critical temperature and g * (T EW ) is the effective number of relativistic degrees of freedom in energy density at the temperature T EW . As we shall see in the following, this is negligible compared to the production via the relaxion-Higgs mixing in the EW broken phase. Now let us consider relaxion production in the EW broken phase, that is, production at temperatures much below the critical temperature of the electroweak phase transition, T EW ∼ m h v/m t = 180 GeV. In order to ensure that any finite temperature effects are negligible, we take T < T 0 = 20 GeV so that we always have (T /T EW ) 2 1. At these temperatures t, h, Z, W ± are not relativistic and their densities are Boltzmann-suppressed. We thus ignore any contribution to thermal production of relaxions from processes involving these states for T 20 GeV and ignore any contribution at all from the temperature range 20 GeV < T < T EW where finite temperature effects become important. We also do not consider any possible contribution from the backreaction sector as this would be impossible to compute model-independently. Consequently, our final result for the relaxion abundance will be a conservative lower bound and the cosmological bounds we derive can possibly be even stronger. For T 20 GeV, the dominant production processes are the Primakoff process q(g) + φ → q(g) + φ, involving the φgg vertex and the Compton photoproduction process q + φ → q + φ which involves the φqq vertex. Using the production rate for the Primakoff process computed in ref. [68], we get, where we have considered only the top loop for computing the φgg coupling as the loop contribution of lighter quarks vanishes for temperatures above their masses. For the Compton process the thermally averaged rate is given by [69], Clearly, the dominant contribution is from bottom quarks and the contribution from lighter quarks is negligible. The interference between the Primakoff and Compton processes also scales as m 2 f T , but is suppressed by another power of α s with respect to Γ f C in Eq. (6.10) and thus we ignore this contribution. We also ignore any contribution from the electromagnetic counterpart of the above processes (that is replacing gluons by photons in the respective diagrams) which are expected to be suppressed by powers of (α em /α s ). Thus, we finally obtain for the total production rate, With the knowledge of Γ we can now compute the abundance of thermally produced relaxions by solving the Boltzmann equation, where x = 1/T and the Hubble scale H t = √ 4π 3 g * (T ) 45 Integrating the above, we get where Y pr = 0.278/g pr * and g pr * 86.25 is the number of relativistic degrees of freedom in energy density in the 1-20 GeV temperature range. In the sum over fermion species we include only the c and b quarks as the contribution due to the other quarks is negligible (see Eq. (6.10)). We have taken the final temperature T f = 1 GeV for Primakoff production to justify our use of perturbative QCD and T = m f for the Compton process because below this temperature the respective fermions become non-relativistic. For s θ 10 −6 the relaxions have an equilibrium density given by Y = Y eq = 0.003 whereas for s θ 10 −6 , the relaxions have a much smaller density, Y hφ = 2.9 × 10 9 s 2 θ . (6.14) Once the Universe cools down to a temperature below the quark/hadron transition, i.e. T 200 MeV, relaxions can be produced via the pion-relaxion conversion process, i.e.
N + π → N + φ, N being a nucleon. Using T we obtain the following parametric estimate for this process, (6.15) One can check that and hence we ignore this contribution. Finally, inverse decays may become significant at temperatures just a bit larger than the relaxion mass. The ratio Γ φ /H t , Γ φ , being the relaxion decay width, is maximal for T m φ /5 as below this temperature, the relaxions become non-relativistic and the rate is Boltzmann-suppressed while above these temperatures H t increases. We check numerically that and thus the contribution from inverse decays can also be safely ignored. We now show that the contribution to relaxion abundance from the q(g)+φ → q(g)+φ processes in Eq. (6.13) by far dominates over the contributions in Eq. (6.6) and Eq. (6.8). First note that we can rewrite Eq. (6.8) as using Eq. (1.5) and Eq. (2.5). As we will discuss in detail in the next subsection, cosmological probes are sensitive only if the relaxion decays after 1 s. As one can see from figure 3, in the region of parameter space which lies in the untuned area defined in Eq. (5.3), if the relaxion decay time is greater than 1 s (below the purple curve) we must have sin θ < 10 −4 . In this region Y H 2 is clearly always smaller than Y hφ in Eq. (6.13), even for a cut-off as low as 3 TeV. As far as the contribution from misalignment, Y m , is concerned we have checked numerically that Y m Y hφ except in a region of the parameter space where none of the cosmological constraints apply as Y hφ < Y m < 10 −20 are both extremely small. Thus we conclude that, under our assumptions, the abundance is well approximated by Eq. (6.13).

Cosmological bounds on late decays
In this subsection we study the bounds on late decays of the relaxion. The earliest the relaxion has to decay to have any effect on cosmology is after 1 s, that is at the neutrino decoupling time, which in the relaxion parameter space corresponds to m φ < 150 MeV as shown in figure 5. On the other hand for relaxion masses m φ < 0.1 keV, Eq. (5.3) implies that sin 2 θ 10 −17 and thus a lifetime, τ φ 10 26 s (see figure 1(b)) much greater than the age of the Universe (10 17 s). This means that for masses m φ < 0.1 keV an exponentially small number of relaxions have decayed by the present time and, as we will soon show more rigorously, there are consequently no bounds in this region. To compute the various constraints from late decays it is important first to know whether the relaxion decays relativistically or non-relativistically at a given point in the parameter space. If the relaxions are relativistic, their temperature can be computed from their number density, If T φ (τ φ ), the relaxion temperature at the time of its decay, is smaller than m φ /5, it can be safely considered to have become non-relativistic before decaying. If it becomes nonrelativistic, it can even dominate the energy density of the Universe before decaying (as the energy density of non-relativistic matter decreases more slowly compared to that of relativistic matter). As we will see in this section, such a scenario is highly constrained. In most of the parameter space where various bounds on late decays are relevant, the relaxion decays non-relativistically and thus its energy density before decaying is ρ φ = m φ Y φ s. Thus the various bounds on late decays generally put an upper bound on m φ Y φ as a function of the lifetime τ φ . Let us now discuss the various constraints on the relaxion decays.
Entropy injection: If the relaxions decay after the neutrinos have fully decoupled, i.e. for τ φ 1 s, they increase the entropy of the SM plasma by ∆S, and thus decrease both the baryon-to-photon ratio η B and the effective number of neutrino species, N eff . Let us now proceed to compute ∆S/S. For τ φ 1 s, relaxions decay nonrelativistically except in a small region of the parameter space with sin 2 θ 10 −4 and m φ < 1 MeV which is outside the region of interest defined in Eq. (5.3). In any case for relativistic decays, and, as we will see, entropy injection smaller than a few percent is unconstrained. To obtain the last inequality above we have used Eq. (6.19). In the rest of the parameter space where the relaxion decays non-relativistically, we must differentiate between the scenario where the relaxion energy density as a fraction of the energy density of radiation, i.e., is smaller than unity, δ 1, from the scenario, δ 1, where the relaxion dominates the energy density. The entropy injection is given by where x = 1.50 [70] for δ 1 whereas x = 1.83 [71] for δ 1.
Having obtained the expression for ∆S/S, let us proceed to derive the constraints from η B and N eff measurements. We first discuss the bound from N eff . Entropy injection anytime after neutrino decoupling and before recombination leads to the reduction in N eff , that is: S before S after 4/3 (6.24) with N SM eff = 3.046. Following ref. [72] we use the bound N eff > 2.6 and show in pink the region excluded by this constraint in figure 3 and figure 5.
We now discuss bounds arising from the decrease in the baryon-to-photon ratio, η B , caused by relaxion decays. Since the baryon-to-photon ratio is inversely proportional to S, η B is reduced as follows due to entropy injection, A change of η B between BBN and CMB epoch is not supported by observation since the measured value of η B during the CMB epoch agrees well with the value after the end of BBN. Therefore, entropy release between these two epochs must be suppressed. In particular, CMB and BBN data constrain ∆S/S to be smaller than 2% [73]. In figure 3 and 5 we show the regions of parameter space excluded by this bound in orange.
Big-bang nucleosynthesis: Big-bang nucleosynthesis (BBN), the formation of light elements in the early Universe, might be altered by late relaxion decays into SM particles. The effect depends strongly on the relaxion mass, particularly whether or not it is heavy enough to cause electromagnetic or hadronic cascades. In our region of interest (i.e. for f > m h ) relaxions above the pion threshold have a lifetime bigger than 1 s (see figure 5), so they do not affect cosmology. Decays of lighter relaxions give rise to electromagnetic showers as long as their mass is bigger than twice the minimum photo-disintegration energy of light nuclei (m φ 5 MeV). In the relaxion parameter space (see the beginning of section 5) m φ < 150 MeV for τ φ > 1s, so we obtain our bounds from relaxion decays into electrons. BBN bounds put constraints on ρ φ /s = m φ Y φ as a function of the lifetime τ φ . We consider here the bounds presented in ref. [73] for the decay of a 140 MeV scalar.
The region f < m h in figure 5, while not relevant for relaxion models, can be interesting in general Higgs portal models. This region can be constrained, for instance, by BBN bounds on decays to pions, hadronic showers etc which can be easily derived using our expression for the abundance in Eq. (6.13).
Distortion of the CMB spectrum: The energy spectrum of the cosmic microwave background (CMB) allows also to constrain energy release in the early Universe. Constraints from CMB distortions become effective for relaxion decays that take place after 10 6 s as at earlier times the thermalization process is very efficient. There are two types of distortions: µ-distortions and y-distortions which dominate at different times. At τ DC = 10 6 s (T γ ∼ 750 eV), the photon number changing double Compton scattering process (γ + e → γ + γ + e) freezes out. As a result, the photons can no longer be in a Planck distribution (where the number of particles is fixed by the total energy). On the other hand, the Compton process is active until τ C = 10 9 s, thus the photons can still maintain a Bose-Einstein (BE) distribution, but with a chemical potential µ, whereas the observed Planck spectrum corresponds to an almost vanishing chemical potential. Therefore, |µ| is constrained by the COBE/FIRAS data which give a bound of |µ| < 0.9 × 10 −4 at 95% CL.
The chemical potential generated by these late decays can be computed to be [74], In the above equation the factor involving exponentials accounts for the fact that only decays in the time period between τ DC and τ C contribute to µ-distortions. If the fractional energy δ 1, one can use ρ φ = m φ Y φ s, and ρ γ = π 2 15 T 4 γ to find the constraints whereas the region δ 1 is excluded as it will lead to an O(1) value for µ which is excluded. We find that a large portion of the parameter space is excluded by this constraint as shown in figure 5 in green.
If the relaxion decays later than τ C = 10 9 s (T ∼ 25 eV), even the Compton process freezes out and this leads to a deviation of the CMB spectrum from a BE distribution. The degree of thermalization that the photons can still achieve can be parametrized by y [74], The region with δ 1 is directly excluded whereas in the region δ 1 we use ρ φ = m φ Y φ s, and ρ γ = π 2 15 T 4 γ to compute the bound. In figure 5, we show the region excluded by the bounds from µ distortions in a darker shade of green than the one denoting y distortions. We also show by dashed lines the projection for the region PIXIE can exclude at 5-sigma level, given by |µ| < 1 × 10 −8 and |y| < 5 × 10 −8 [75].
EBL and reionization: After recombination (τ RC ∼ 10 13 s) the nuclei capture almost all the electrons to form neutral atoms so that the Universe becomes nearly transparent to radiation. The photons injected by relaxion decay can be in principle directly detected, unless their wavelength lies in the ultraviolet range (13.6 eV-300 eV) and they are absorbed in the photoionization process of atoms. In this ultraviolet mass range bounds from reionization can be set. Photons emitted from very late decays that do not lie in this range, can be observed today as a distortion of the diffuse extragalactic background light (EBL). The above constraints can be used to bound the quantity m φ Y φ /τ φ of as a function of m φ . Together these bounds cover the wavelength range between 0.1 and 1000 µm, that is roughly the mass range between 0.1 eV and 1 keV. We show in figure 5 the excluded region using the bounds derived in ref. [74] and [76], but appropriately rescaled to the different abundance in our case.
Dark matter: if the relaxion decays after ∼ 10 17 s it forms a very small component of the present dark matter density.

Astrophysical probes
SN1987a supernova: In the core of a supernova, a relaxion can be produced via its couplings to nucleons and thereby contribute to its energy loss. The relevant process is bremsstrahlung N + N → N + N + φ. Requiring that the energy loss into the new scalar must be smaller than the measured energy loss into neutrinos leads to bounds on the Higgsrelaxion mixing as long as the relaxion is lighter than 20 MeV. In figure 5 and 3 we show (in light blue) the bounds derived in [77], using the results of ref. [78]. This computation is exponentially sensitive to some uncertainties (see ref. [77]) and thus should be interpreted only as an order of magnitude estimate. At a more conceptual level, even the idea of energy loss via neutrinos has been questioned in the literature [79]. New laboratory constraints that are able to explore this region are therefore required.
CAST experiment: The CERN Axion Solar Telescope (CAST) looks via X-rays for axion-like particles coming from the sun. The present limit on the photon-ALP coupling is [82,83]: for m φ < 0.02 eV. The limit is slightly weaker than the GC limit and well outside the region of interest in Eq. (5.3), hence we omit it in figure 5. In contrast, IAXO [84], the new generation experiment, will be able to improve the limit. However, despite the future progress in this technology this class of experiments is not likely to be relevant for our scenario since it probes a region of the parameter space where fifth force experiments provide very strong bounds.

Implications for the relaxion theory space
In this section, we collect all bounds from laboratory experiments, colliders, astrophysics and cosmology that were shown in figures 2, 3, 4 and 5 for different mass regions and translate them, using Eq. (5.2), into the underlying theory parameters Λ br and f in figure 6. 9 As a connection between both parametrisations, Λ br and f were shown as a grid of contours in the previous plots, whereas in the (Λ br , f ) plane of figure 6 we show contours m φ . While the values we provide are for the j = 2 case, as mentioned below Eq. (5.2) one can obtain the values for the j = 1 by the simple translation Λ br → √ 2Λ br , f → 2f . We show how these bounds push the cut-off to smaller values by the upper horizontal axis, where we translate the Λ br scale in the lower axis to cut-off values using Eq.  The overview presented in figure 6 shows that large areas in the Λ br − f plane are already well covered by existing experimental and observational probes, for instance the high-f region up to M Pl is probed by the fifth force experiments, on the other hand the cosmological, astrophysical, beam dump and collider observables constrain lower values of f . We see that in the above f ranges, the region with electroweak scale Λ br is practically completely ruled out apart from small gaps that still remain. We also show in figure 6 how some of these gaps in parameter space might be covered soon by future experiments such as SHiP, NA62 and PIXIE. However, the region between f ∼ 10 10 GeV and 10 14 GeV which corresponds to relaxion masses between 0.1 eV and 1 keV, is currently barely constrained by data.
For any f (or m φ ) value, all the constraints can be evaded for sufficiently small Λ br values (there are no bounds for Λ br 0.3 GeV). Small Λ br values are however theoretically disfavoured for several reasons. First of all, as we see from the Λ cq contours in figure 6, the constraints derived here push the relaxion to a region with somewhat lower values for the upper bound on the cut-off derived from cosmological considerations during inflation.
If one takes seriously the requirement that the relaxion should not have transplanckian excursions, our bounds have a much stronger impact. This is because, as we see from figure 6, our bounds already cover a large part of the parameter space outside the shaded region where the relaxion travels transplanckian distances for any cut-off larger than 2 TeV. Coming to the issue of the very large global charges that arises due to the compact nature of the relaxion, we see from the upper horizontal axis that even in CKY/clockwork models the number of sites required can become uncomfortably large for very small backreaction scales. For N 30 (see section 4) our bounds can significantly constrain the cut-off. For instance for f = 1000 TeV we find Λ 100 TeV. As shown in section 4 the simplest clockwork models start getting tuned for N 30. As far as the proposal to solve the little hierarchy problem using modest n values is concerned [3], we see that such a proposal would be completely ruled out outside the f ∼ 10 10 GeV-10 14 GeV (m φ ∼ 0.1 eV -1 keV) region, as contrary to the philosophy of this approach, too large values of n > (v/Λ br ) 4 would be required.
One should be keep in mind while interpreting these bounds within the clockwork framework that in these models one must have f Λ from Eq. (4.10). Thus even from this point of view the unconstrained f ∼ 10 10 GeV-10 14 GeV (m φ ∼ 0.1 eV-1 keV) window is an interesting region as here the cut-off can be high in these models.
Finally let us discuss what impact the pseudoscalar couplings of the relaxion might have on the overall bounds. As explained in appendix C, in the electroweak preserving [1,3] models discussed in section 3, the relaxion does not have pseudoscalar couplings larger than the Higgs-portal ones, hence our experimental bounds would be qualitatively unchanged. Let us briefly comment on the possible change in our bounds if the pseudoscalar coupling to photons is larger than the one induced by Higgs mixing. As already mentioned, among the models discussed in section 3 this holds only for the pseudoscalar diphoton coupling in the non-QCD j = 1 model where the relaxion has a pseudoscalar coupling to photons suppressed only by 1/f and not by the backreaction scale (see Eq. (C.3)). In this case the astrophysical and cosmological bounds discussed in section 6 will be affected. An analysis of how the cosmological bounds change in the presence of a largeg φγ coupling is beyond the scope of this work. The enhanced coupling to photons will lead to a stronger bound from globular clusters, that is f 10 7 GeV, Eq. (6.29). However, this is valid only provided that the relaxion mass is lighter than 30 keV, so we immediately see from figure 6 that this is relevant only for Λ br v. Furthermore, the CAST experiment can put a bound on f of similar order on the coupling to photon Eq. (6.30) in the sub-eV region. For large Higgsrelaxion mixing fifth force experiments are sensitive, hence the CAST bound is irrelevant. However for Λ br v, when the sensitivity to fifth force experiments ceases, the CAST bound on the pseudo-scalar coupling can be important for sub-eV relaxions.
8 Testing for the CP violating nature of the relaxion In this section we investigate the feasibility of detecting a signal of spontaneous CPviolation together with a Higgs mixing signal. This would represent a smoking gun for our scenario since what we discussed so far about relaxion phenomenology applies to any scalar mixed with the Higgs. However, as already discussed in appendix C, the strength of the relaxion pseudo-scalar couplings depend on the details of the back-reaction sector. Couplings to fermions are typically very suppressed (compared to the one from Higgsrelaxion mixing), while the coupling to photong φγ is in many cases only as large as the scalar one, that isg φγ ∼ 10 −5 sin θ. In the electroweak breaking non-QCD model discussed in section 3, instead, the coupling to photons is in principle larger since it is not suppressed by the backreaction scale. Despite the model dependence, it is still an interesting question whether a CP-violating signal could be detected at the precision frontier. Let us investigate the relaxion contribution to the electric dipole moments (EDM). In our scenario the leading contribution to the electric dipole moment is generated through its couplings to fermions via Higgs mixing and with the pseudoscalar coupling to photons,g φγ .We will focus on the electron EDM, following [85], but similar results hold for the neutron EDM.
The first step is to understand in which relaxion mass range this probe can be effective. To this end let us estimate the strength ofg φγ ×g φe since the relaxion one-loop contribution to the electron EDM will be proportional to it. The current upper bound on the electron EDM is d e /e ∼ 8 × 10 −29 cm [86], which corresponds tog φγ g e ∼ 5 × 10 −14 GeV −1 [85], and improvements of one order of magnitude are expected in the coming years [87]. Let us then see how this compares to relaxion models. For the non-QCD electroweak breaking model we get:g where we used Eq. (2.6) and Eq. (C.3). The electroweak preserving models [1,3] have an additional Λ 4 br /v 4 suppression from the backreaction scale due to the suppression ing φγ in Eq. (C.2) as compared to Eq. (C.3). We see that in both cases, a relaxion with m φ 1 GeV yields a contribution to the d e that is below the current (and near future) sensitivity. The parameter space constrained is therefore in the few GeV region.

Conclusions
We study various phenomenological aspects of relaxion models. We focus on models where the rolling of the relaxion field stops due to the presence of a Higgs-relaxion backreaction term. We show that the relaxion generically stops its rolling at a point that breaks the CP symmetry, leading to relaxion-Higgs mixing. We investigate then the implications of this mixing, and analyse current and near future probes involving laboratory, cosmological and astrophysical measurements in terms of reach and sensitivity. In most parts of the parameter space, these observational constraints put the most stringent bound on the backreaction scale, Λ br . On the theoretical front, we show that simple multiaxion (clockwork) UV completions suffer from a fine tuning problem, which increases with the number of sites.
Let us describe in more detail our main results on the observational probes of relaxion-Higgs mixing. The constraints/discovery prospects derived by us are summarised in figures 2-5. In the sub-eV mass range the relaxion lifetime is much larger than the age of the Universe and thus cosmological or direct laboratory probes are not effective. Fifth force experiments, however, are sensitive in large regions of the parameter space in the sub-eV region because of the low mass of the relaxion and the CP-even nature of its couplings to SM particles via Higgs mixing (see figure 2). The eV-MeV region is practically unconstrained by laboratory probes, but a subset of this region (keV-MeV) can be constrained by astrophysical and cosmological probes as shown in figure 5. The cosmological probes are relevant here because this is the region of parameter space where the relaxion lifetime is between 1 s and 10 26 s and thus is tested by a variety of cosmological probes, such as entropy injection constraints from N eff and η B measurements, BBN observables, CMB spectral distortions and EBL distortions. Turning to the MeV-GeV region we find that in some parts of this mass range, the relaxion lifetime is just right for beam dump experiments (O(100 m) in the lab frame) such as the CHARM experiment and experiments probing invisible rare meson decays. We also find that future data from beam dump experiments like SeaQuest and especially SHiP and the currently running ultra-rare kaon decay experiment NA62 can probe new and interesting regions of the relaxion parameter space. In other parts of this MeV-GeV mass region visible rare meson decays also put significant bounds. Finally, for relaxion masses above 5 GeV the constraints arise from LEP bounds on the Higgs-strahlung process and LHC Higgs coupling bounds on the new channel, h → φφ, as shown in figure 4. In figure 6, we translate these bounds to the relaxion theory space and discuss the theoretical implications. We finally comment that, while the relaxion-Higgs mixing requires CP violation, most of the probes discussed above do not form a strong test of the CP nature of the relaxion. The pseudoscalar couplings of the relaxion tend to be more model-dependent. For instance, in the familon model that was constructed in [3] the relaxion does not couple to F µνF µν (with F being the QED field strength) at one loop but only to the orthogonal combination of the electroweak field strengths. We find that, in existing models, probes of CP violation are sensitive only for GeV scale relaxion masses.
suppression (∼ Λ 4 br /v 4 ) and at the same order in perturbation theory as the Higgs portal coupling (in the non-perutbative model the coupling can arise via mixing with the analog of the η and in the perturbative familon model via a 2-loop level diagram), In the non-QCD j = 1 model, however, it is possible to haveg φγ g φγ because this backreaction sector is just a scaled-up version of QCD. Thus, as is the case for QCD axions, the relaxion will get an anomaly-induced coupling of the same order via mixing with the η and pion analogs of the new strong sector. This generates which can be larger than the Higgs portal coupling, g φγ , for values of Λ br v. It is important to mention that while the pseudoscalar coupling of the relaxion to photons is smaller than the Higgs-portal one in the existing j = 2 models, an anomaly induced coupling of the size in Eq. (C.3) would exist in simple variants where the relaxion couples directly to the the electroweak doublet fermions.
One can proceed along the same lines to show that the pseudoscalar coupling of the relaxion to gluons is at least one loop suppressed with respect to the Higgs portal induced coupling to gluons because of the sequestering. We see, therefore, that apart from thẽ g φγ coupling in the j = 1 model, the pseudoscalar couplings of the relaxion are either suppressed or of the same order as the Higgs portal coupling in the models in section 3. Our results would thus be qualitatively unchanged by the presence of these couplings apart from the one exception above, on which we comment in the text.

D The hφφ coupling in j = 2 models
In this appendix we present the expression for the hφφ coupling in j = 2 models. To obtain this we expand the potential V (h, φ) in Eq. (2.1) around the minimum (v, φ 0 ) to obtain all cubic terms and then substituting the gauge eigenstates in terms of the mass eigenstates φ = −s θ h + c θφ , h = c θ h + s θφ . (D.1) In order to reduce the complexity the full expression while accounting for the leading mixing effects, we take sin(φ 0 /f ) = cos(φ 0 /f ) = 1/ √ 2 to finally obtain

E Expressions for relaxion partial widths and lifetime
In this appendix we provide the expressions for the relaxion partial widths for different channels. The dilepton (ll) and diphoton (γγ) partial widths are given by . (E.1) As far as colored states are concerned we use the perturbative description above m φ = 1 GeV. The partial width to quarks (qq) and gluons (gg) is given by For m φ < 1 GeV the only hadronic state we consider is the decay to pions. Different estimates of the partial width to pions vary over nearly two orders of magnitude [32]. Here we use the leading order calculation of ref. [88] which gives .

(E.3)
For m φ > 1 GeV one should use the partial width to kaons, η-mesons etc, but as no reliable estimate exists in this regime [32], our perturbative estimate is sufficient in this context. For a given mass, the total width, Γ φ , can now be obtained by summing over all the kinematically relevant decay modes. Analyzing the ratio Γ φ /M φ , we find that the relaxion is very narrow throughout the whole parameter space of our interest. For example, Γ φ m φ 2 · 10 −13 , 10 −5 · sin 2 θ for m φ = {0.1, 5} GeV .

(E.4)
For lighter masses, this ratio becomes even smaller. Hence, potential width effects do not arise. The ratio for intermediate masses between the two example values in Eq. (E.4) highly depend on the thresholds of those particles that the relaxion can decay into. The relaxion lifetime, τ φ = 1/Γ φ , is crucial in determining the applicability of various observational constraints. We show the lifetime as a function of m φ for different sin 2 θ values in figure 1(b).