Gravitational Wave Echo of Relaxion Trapping

To solve the hierarchy problem, the relaxion must remain trapped in the correct minimum, even if the electroweak symmetry is restored after reheating. In this scenario, the relaxion starts rolling again until the backreaction potential, with its set of local minima, reappears. Depending on the time of barrier reappearance, Hubble friction alone may be insufficient to retrap the relaxion in a large portion of the parameter space. Thus, an additional source of friction is required, which might be provided by coupling to a dark photon.The dark photon experiences a tachyonic instability as the relaxion rolls, which slows down the relaxion by backreacting to its motion, and efficiently creates anisotropies in the dark photon energy-momentum tensor, sourcing gravitational waves. We calculate the spectrum of the resulting gravitational wave background from this new mechanism, and evaluate its observability by current and future experiments. We further investigate the possibility that the coherently oscillating relaxion constitutes dark matter and present the corresponding constraints from gravitational waves.


I. INTRODUCTION
The weakness of the gravitational force ensures that gravitational waves (GWs) can store information about the history of our Universe. On the other hand, though, it also implies that only dramatic and cosmological events may leave imprints that can be observed at present or in the near future. GW detectors are, for instance, sensitive to primordial phase transitions, inflation, and oscillatory motion of ultralight fields (see, e.g., Refs. [1,2] for a recent discussion of potential signals of cosmological origin).
In this work, we introduce a different mechanism that leads to the production of GWs in the early Universe. In the following, we demonstrate the idea using the relaxion framework [3]; however, the core scheme may be applicable to other forms of new physics. At early times, the relaxion is trapped in a local minimum. However, assuming that the reheating temperature is above the electroweak (EW) scale, 1 we expect the EW symmetry to be restored. Hence, the minimum in which the relaxion was originally trapped disappears, and the relaxion starts rolling down its potential. At some temperature below the EW scale, we expect the set of local minima to reappear. There is a limited region of the model's parameter space such that the relaxion would be trapped still within the same min- * abhishek.banerjee@weizmann.ac.il † eric.madge-pimentel@weizmann.ac.il ‡ gilad.perez@weizmann.ac.il § w.ratzinger@uni-mainz.de ¶ pedro.schwaller@uni-mainz.de 1 The requirement of a high reheating temperature, motivated by a large class of inflation models (see, e.g., Ref. [4] and references therein), is generic and also needed in most models that explain the observed baryon abundance (see ,e.g., Ref. [5] and references therein).
imum [6,7]. In this case, even the minimal realization of the relaxation mechanism could lead to a viable ultra light dark matter (DM) candidate [6], which may be tested in the future due to the presence of relaxion-Higgs mixing [8,9]. Despite the attractiveness of this minimal setup, there is a sizeable region of parameter space where the postreheating displacement of the relaxion would be too large and the relaxion would not be trapped but instead would descend in an uncontrolled way toward the global minimum. This runaway can be avoided if an additional source of friction is added to the model. Coupling a scalar field to a dark U (1) gauge field strength is known to lead to an efficient energy loss mechanism in the form of tachyonic dark photon field production [10,11] (see Refs. [12,13] for the case of a relaxion or other axionlike particles (ALPs) in the mass range O(MeV), where self-friction may be induced, or Ref. [14] for stopping via an instability in a modified relaxion potential, albeit outside the region of interest of this work). Dark photon production quickly reaches a quasi-steady-state where the friction balances the slope of the potential. At each time, the dominantly produced momentum-mode k = ξaH is about to exit the tachyonic band, where ξ is an O(10 − 100) parameter, H is the Hubble rate, and a is the scale factor.
The energy density stored in this mode is roughly constant, ρ X ∼ m 2 φ f 2 φ [6], where m φ and f φ are the relaxion mass and the decay constant, respectively, and is the source of GW production. The rolling of the relaxion and the dark photon production stop around the time when the potential barriers reappear, which traps the relaxion and ends this epoch. 2  The grey dashed line encloses the parameter space in which relaxation can be realized without dark photon friction, which is discussed in more details in Appendix A. The prospective GW sensitivity of µAres (green) as well as SKA after an observation period of 5 years (turquoise) and 20 years (blue) is indicated by the respective coloured regions. In the purple coloured region, a subrange of the viable reappearance temperatures can be excluded using current NANOGrav data from the 11 year data set. The regions bounded by the coloured dotted lines need super-Planckian decay constants to be probed by the respective experiment. In the lower panel, the grey shaded region can reproduce our best-fit spectrum (at Tra ∼ 20 MeV) to the potential stochastic GW background seen in the recent NANOGrav data. An animated version of the lower panel of the figure, explicitly showing the dependence on the reappearance temperature, can be found in the ancillary material of this paper.
Given the above mechanism, we can estimate the GW signal as follows. First, the peak frequency at present time is obtained by redshifting the dominant k-mode at the time of reappearance, cf. Eq. (35), where T ra is the temperature when the potential barriers reappear. Here and in the following, quantities indexed ance process of the backreaction potential and assume instantaneous reappearance. Various considerations regarding barrier reappearance can also generate gravitational waves, a scenario considered in Ref. [15].
"0" and "ra" are evaluated today and at the time when the potential barriers reappear, respectively. Second, the GW amplitude is roughly given by the square of the energy density stored in the dark photon just before T ra , cf. Eq. (28), where ρ 0 c = 3M 2 Pl H 2 0 is the critical energy density of the Universe today, f ra peak = ξH ra is the peak frequency without red-shifting, and c 2 eff 1/4 is the efficiency factor for converting dark photon energy into GWs obtained from the analytic calculation in Appendix C. Note that c eff ∼ O (1) as the relaxion potential energy at each time is predominantly deposited in a narrow range of exponentially growing dark photon momentum modes, resulting in large time-dependent inhomogeneities and an efficient generation of anisotropic stress sourcing GWs, in a similar manner as in the case of tachyonic preheating [16][17][18][19]. In the last line, we have reexpressed the signal strength in terms of the relaxion-Higgs mixing sin θ hφ and the cutoff scale Λ. This result shows that some of the parameter space of the model may lead to a visible signal in near-future GW experiments, allowing us to probe parameter regions that are currently unexplored by other experiments, as discussed below. In addition, we note that the relation between the physical parameters of the models and the GW amplitude is given by , showing a rather mild dependence on the actual amplitude.
For convenience, the range of relaxion masses, m φ , and mixing angles with the Higgs boson, sin θ hφ , which can be probed via current or future GW experiments as well as the corresponding constraints on the parameter space are summarized in Fig. 1. The green and blue/turquoise coloured regions can be accessed with µAres and the Square-Kilometre Array (SKA) observatory, respectively, depending on the temperature of barrier reappearance. In the purple region, the reappearance temperature is restricted by current data from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). In addition to that, in the grey shaded region, we present the parameter range in which our model can account for the potential GW signal recently observed in NANOGrav data. The grey dashed line encloses the region in which the relaxion can be trapped without dark photon friction. On the other hand, as the figure illustrates, there is a large fraction of parameter space where an additional source of friction is required for the viability of the relaxion mechanism and thus motivates us to add a relaxion-dark photon coupling to the model. A more detailed discussion of the figure is deferred to Sec. IV. This paper is organized as follows. We briefly review the relaxion and the dark photon dynamics in Sec. II. Sections II A and II B contain a brief discussion of the interplay of the relaxion-dark photon dynamics and the possibility of the relaxion being DM, respectively. In Sec. II C we review the constraints on our model and discuss the available parameter space. Subsequently, the production mechanism of the GW background is studied in Sec. III. We derive the GW spectrum in Sec. III A and briefly describe how the detectability of the signal is evaluated in Sec. III B. The results of this paper are then discussed in Sec. IV. Section V concludes the paper. A brief discussion of the minimal relaxion scenario (without a dark photon coupling) is deferred to Appendix A. Further details regarding the calculation of the dark photon and GW spectra are provided in Appendices B and C, respectively.

II. SETUP
In this work, we consider the relaxion φ coupled to a dark photon field X µ , with the potential of the relaxion field φ and Higgs doublet H given by where λ is the Higgs quartic coupling and Here, c is an O(1) number, g is a dimensionless parameter, Λ is the Higgs mass cutoff scale, Λ br is the backreaction scale, v H = |H| = 174 GeV is the Higgs vacuum expectation value, and f φ is the decay constant of the relaxion. During inflation, the relaxion rolls down the linear slope of its potential V roll . It thereby scans the Higgs mass parameter µ 2 (φ). Once µ 2 crosses zero, the Higgs acquires a nonvanishing vacuum expectation value, triggering the breaking of the EW gauge symmetry. The Higgs then backreacts creating wiggles in the relaxion potential via V br . Once the Higgs backreaction balances the rolling potential, the relaxion is trapped in the first minimum it encounters. Choosing cgΛ 3 f φ ∼ Λ 4 br , we end up with a weak-scale expectation value for the Higgs boson, solving the hierarchy problem. The relaxion mass and the relaxion-Higgs mixing angle can then be written as [7,9] in terms of the theory parameters. Here, m h = 125 GeV is the Higgs mass.

A. Relaxion and dark photon evolution
After reheating, the EW symmetry will be restored due to thermal corrections to the potential, provided that the reheating temperature is above the EW phase transition temperature. As a consequence, the relaxion will start rolling again, leading to exponential production of dark photon modes. To see the interplay, the coupled differential equations describing the evolution of the spatially homogeneous relaxion and the dark photon modes are given by where θ = φ/f φ and primes denote derivatives with respect to conformal time τ with a dτ = dt. We have written the dark photon in terms of its Fourier modes X λ (k, τ ) in Coulomb gauge, ∇ · X = 0, aŝ Here, λ denotes the dark photon helicity, and ε λ are the corresponding circular polarization vectors. As evident from Eq. (8), dark photon modes with 0 < k < λr X θ are tachyonic and will experience exponential growth compared to the vacuum fluctuations. The resulting dark photon spectrum then features anisotropies in its energy-momentum tensor which will act as a source for GW production, leading to a stochastic GW background [17][18][19][20][21][22]. Furthermore, since only modes of the helicity with the same sign as the relaxion velocity can become tachyonic, the rolling relaxion will produce a circular polarized dark photon background. In our case, as we choose θ > 0, only the positive-helicity modes are exponentially produced. For these modes, the solution to the equations of motion is given in the Wentzel-Kramers-Brillouin (WKB) approximation by where Ω 2 (k, τ ) = k r X |θ (τ )| − k 2 > 0 is the corresponding tachyonic frequency and g(k, τ ) = τ dτ Ω(k, τ ). The approximation holds for |Ω /Ω 2 | 1. At early times, just after reheating, the friction from dark photons can be neglected. Assuming radiation domination, the relaxion velocity can then be written as where we imposed θ (τ rh ) = 0. Subsequently, due to the exponential production, the dark photon friction becomes comparable to the other terms in the equation of motion, Eq. (7), and asymptotes to the slope of the rolling potential. The timescale τ pp at which particle production kicks in can be determined by wherek is the mode that dominates the XX term, given by the saddle point ∂g(k, τ )/∂k|k = 0. After τ pp , due to the balance between the potential slope and the backreaction from the dark photon, the relaxion field velocity becomes proportional to the Hubble rate and evolves as [6,10] with a small logarithmic correction ( 1). We defined here the parameter ξ = r X |θ | aH at τ pp . From Eq. (12) we then obtain ξ ∼ O(10 − 100) with a mere logarithmic dependence on the relaxion parameters [10]. The dominating k-mode at each epoch then becomes Once the Universe has cooled sufficiently, the EW phase transition occurs, and the wiggles of the backreaction potential reappear. The rolling of the relaxion between reheating and T ra leads to a displacement from the minimum in which it originally settled during inflation by For the relaxion to remain trapped in this minimum, we need to require that the displacement is less than the distance between the minimum and the next maximum, ∆θ 2δ, where δ = Λ 2 br /(Λv H ) [6,7]. This sets a lower bound on the coupling to the dark photons. For smaller couplings, the dark photon friction is insufficient to prevent the relaxion from rolling into one of the neighbouring minima. The relaxion then needs to traverse ∆θ ∼ O(n) to end up in the n-th minimum, where n = 1 denotes the minimum in which it stopped during inflation, extending the parameter space of the theory. However, going beyond the first minimum requires a careful adjustment of the initial conditions to let the relaxion stop exactly at the bottom of the n-th minimum at reappearance. Otherwise, the time required for the relaxion to reach the bottom would exceed the age of the Universe. We thus simply assume r X = ξ/(2δ) in the following. 3

B. Relaxion dark matter
After the reappearance of the Higgs backreaction potential, the displaced relaxion begins to oscillate around the minimum of its potential, providing a candidate for ultralight DM as discussed in Ref. [6]. Assuming the 10 . Such large couplings can be obtained in a technically natural way for example via the clockwork mechanism [23][24][25][26]. maximal displacement of ∆θ = 2δ, the relaxion relic abundance is given by where T osc ∼ min[T ra , m φ M Pl ] is the temperature when the relaxion starts to oscillate. Requiring that the relaxion provides all of DM, i.e., Ω φ h 2 = 0.12 [27], the relaxion decay constant can be expressed as As we are considering coherently oscillating relaxion DM here, its mass needs to be less than approximately 10 eV in order to be described by a classical field [28]. For this range of relaxion masses, the possible decay channels are into two photons and two dark photons, Γ φ = Γ γγ + Γ XX . The decay rate into a dark photon pair is given by while the decay width of φ → γγ is subdominant compared to that of φ → XX as it is suppressed by the square of the relaxion-Higgs mixing angle [6], which in turn is bounded from above by sin θ hφ v H /f φ [9]. The relaxion lifetime hence becomes where we have chosen r X = ξ/(2δ). Since the decay of DM into relativistic particles affects the spectrum of the cosmic microwave background (CMB) at low-multipoles, the lifetime is constrained as τ φ > 160 Gyr [29]. As shown in Refs. [22,[30][31][32], for r X (∆θ) sep ∼ O(10 2 ), an oscillating ALP may introduce a second phase of tachyonic dark photon production, which could suppress the DM abundance by up to two orders of magnitude. This condition is satisfied in parts of the parameter space where the field displacement is maximal, since there r X ∼ O(ξ/δ) and ξ ∼ O(10 − 10 2 ), and introduces some uncertainty in our estimate of the DM abundance in those regions.

C. Constraints
A successful cosmological relaxation of the Higgs mass requires the inflation sector to dominate the total energy density, 3H 2 I M 2 Pl Λ 4 , as well as that the classical motion of the relaxion dominates over quantum fluctuation during inflation, (∆φ) cl H I /2π. Here, H I is the Hubble scale during inflation. Combining these two constraints, we get an upper bound on the cutoff scale Λ, As we are considering a Higgs-dependent backreaction potential, we also require Λ br v H [3,9]. 4 The allowed range of the effective cutoff Λ of the theory is Also, for a generic backreaction potential which does not change the late-time cosmology, the range of reappear- For masses below the eV scale, the relaxion can further mediate long-range forces. Experiments looking for such interactions (fifth force experiments, inverse-square-law, and equivalence-principle tests) constrain the coupling of the relaxion to ordinary matter [33][34][35][36][37][38][39], which is induced by the relaxion-Higgs mixing angle given in Eq. (6). In a similar manner, for masses up to the keV range, the mixing is constrained from stellar cooling [40][41][42][43], as it induces relaxion-mediated contributions to the energy loss in stars. Slightly weaker limits on the mixing angle can furthermore be obtained from bounds on the solar relaxion flux as constrained by XENON1T and other darkmatter direct detection experiments [44].
Additional constraints arise when coupling the relaxion to a dark photon field (see Appendix A for a discussion of the minimal scenario without this coupling), with the coupling here chosen to saturate the trapping condition, i.e., r X = ξ/(2δ). For the dark photon induced friction to trap the relaxion, reappearance has to occur sufficiently late for the dark photon to be produced, i.e., H pp > H ra . This sets a lower bound on the relaxion mass, m φ √ 10 δH ra , for which the dark photon scenario can be applied.
If this condition is satisfied, throughout its evolution from reheating to reappearance, the relaxion continuously produces dark photons, depositing energy density into the latter. At the time of trapping, t ra , the dark photon energy-density can be estimated as For the Universe to remain radiation dominated, we require ρ X (T ra ) to be smaller than the radiation energy density at T ra , i.e. ρ X (T ra ) 3M 2 Pl H 2 ra ∼ T 4 ra ; otherwise, the relaxion dominates the total energy density of the Universe. After T ra , dark photon production stops, and the corresponding energy density redshifts as that of radiation.
The dark photon is further constrained by its contribution to the energy density of the Universe at low temperatures, typically parametrized in terms of the effective number of neutrino species N eff . For temperatures T m e , after electrons and positrons have annihilated, the dark photon energy density (or of any additional relativistic species for that matter) can be written in terms of its contribution ∆N eff to N eff as ρ X (T ) = 7 π 2 120 ∆N eff 4 11 where the 4/11 factor accounts for the difference between the neutrino and photon temperature. Using Eq. (22) and the current 95 % C.L. limit N eff = 3.27 ± 0.29 from CMB data [27] combined with local measurements of the Hubble rate [45], as well as N SM eff = 3.046 in the Standard Model, we obtain a lower limit on the reappearance temperature as a function of the relaxion mass and decay constant where g * s,ra is the number of entropic degrees of freedom at reappearance. If we assume that the relaxion accounts for the full DM abundance, then plugging Eq. (16) to Eq. (24), we get, Here, we also assume that the relaxion starts to oscillate at T ra , which is only true for a sufficiently heavy relaxion. As we require Λ 4πv H ∼ 1 TeV, this sets a lower bound of T ra 240 MeV on the reappearance temperature for the relaxion DM scenario to be realized. 5 Upon the same assumptions, Λ br v H further sets an upper bound on the reappearance temperature, T ra min v H , 80 GeV 10 6 GeV Λ 96 g * s,ra 1 3 . (26) In Fig. 2, we show the minimal and maximal allowed reappearance temperature for relaxion DM as a function of the cutoff Λ of the theory. Combining Λ br < v H (blue) and the N eff constraint (green), we see that the highest Λ for which the relaxion can be realized as coherently oscillating DM is Λ 10 7 GeV, which is in accordance with the constraints for T ra 6 GeV. Because of the rapid change in the radiative degrees of freedom around the time of the QCD phase transition, the ∆N eff limit on the reappearance temperature saturates at T ra ∼ T QCD for Λ 2 TeV. We also depict the weaker bound ρ X < ρ rad from the total energy density in orange.
In Fig. 3, we depict the constraints on the relaxion parameters as a function of the relaxion mass m φ and the mixing angle sin θ hφ , where we determined Λ and f φ from the DM abundance using Eqs. (6) and (16). The red and orange shaded regions are excluded by the indicated constraints, where the shape of the exclusions now partially differs from the corresponding ones in Fig. 1 as we require here Ω φ = Ω DM . Scanning over all allowed values of T ra , the full range of masses and mixing angles for which we can obtain coherently oscillating relaxion DM is indicated by the black-framed regions in Fig. 3. We obtain two separated islands of viable parameter space, one at low masses, 10 −11 eV m φ 10 −6 eV, with mixing angles around 10 −23 sin θ hφ 10 −18 , and another island at high masses, 10 −2 eV m φ 1 eV, with a narrow range of mixing angles around 10 −14 sin θ hφ 10 −11 . Note that, in the high-mass island, the coupling to dark photons is required to trap the relaxion, whereas in most of the low-mass island, relaxion DM can be realized without dark photon friction [6]. It shall, moreover, be emphasized that the low-mass DM regions in the minimal and dark photon scenario are separated in the reappearance temperature, since relaxion stopping via Hubble friction requires m φ √ 8 δH ra , whereas dark photon production only occurs for m φ √ 10 δH ra . See Appendix A for a detailed discussion on the minimal scenario.

III. GRAVITATIONAL WAVES
Let us now consider the stochastic GW background generated from anisotropies in the energy-momentum tensor of the dark photon produced during the postinflationary evolution of the relaxion. In particular, we focus here on the GWs sourced during the rolling of the relaxion between reheating and the EW phase transition. Further gravitational radiation can be generated when the relaxion oscillates around the minimum of its potential, as explored in Refs. [20][21][22][46][47][48][49] in the context of general ALP models, or in a confining phase transition that generates the backreaction potential [15]. In combination, the various sources may lead to an interesting and rich GW spectrum with multiple peaks. The spectrum of a stochastic background of GWs is characterized by its fractional energy density, i.e. the energy density normalized to the critical value, ρ 0 c = 3M 2 Pl H 2 0 , per logarithmic frequency interval, where the total GW energy density is given by Here, h ij denotes the GW metric perturbations, and the dot indicates the derivative with respect to cosmic time t. Switching to conformal time τ , dt = a dτ , where a is the scale factor of the Universe, the metric reads During radiation domination, the Einstein equations in the linear regime for the metric perturbations in momen-tum space using transverse-traceless gauge become where k = |k| is the comoving wave number. The anisotropic stress tensor Π ij relates to the energymomentum tensor T ij , via Π ij (k, τ ) = Λ ab ij (k)T ab (k, τ ), where Λ ab ij = P a i P b j − 1 2 P ij P ab with P ij = δ ij − k i k j /k 2 being the projector that extracts the transverse and traceless part [1]. The equations of motion are then solved by (neglecting the a term which vanishes in a radiation-dominated universe i.e. for a ∝ τ ) where G(k, τ, τ ) = sin[k(τ − τ )]/k is the causal Green's function. We denote here the operator form of any quantity Q byQ.

A. Gravitational wave production
The GW energy density per logarithmic interval in the comoving momentum k of a generic stochastic source at conformal time τ is given by [1] where τ i is the time at which the GW source starts operating and Π 2 (k, τ , τ ) is the unequal time correlator of the anisotropic stress, defined as 0|Π ab (k, τ )Π * ab (k , τ )|0 = (2π) 3 δ(k −k )Π 2 (k, τ, τ ). In our case, the GWs are generated between reheating and reappearance; hence τ i = τ rh and τ ≤ τ ra . As the GWs produced before the relaxion reaches its terminal velocity will, however, be subdominant, we can take τ i = τ pp , so to first approximation, the GW signature becomes independent of the temperature to which the Universe was reheated.
The dark photon anisotropic stress sourcing the GWs can be written in terms of the dark electric and magnetic fields aŝ Focusing on the dominant modes which have completed their phase of maximal tachyonic growth, q r X |θ |/2, we find that GWs with momentum k are dominantly produced at the time of reappearance. For frequencies below the peak, both dark photon modes in Eq. (33) have momenta close to the one that experiences maximal growth at reappearance, |q| ∼ |k − q| ∼ r X |θ ra |/2, whereas above the peak, one of the contributing modes must have a larger momentum and therefore have exited the tachyonic band before reappearance. Details of the calculation are deferred to Appendix C. The present-day GW spectrum can then be written as with the peak frequency where k ra is the mode that exits the tachyonic band at reappearance, and the peak amplitude We used Eq. (13) here as well as ξ 2δr X and Λ 4 br = m 2 φ f 2 φ /δ to recast Eqs. (35) and (36) into the estimations in Eqs. (1) and (2). The spectral shape S ξ is given by . (37) Note that, similar to the oscillating axion case [20][21][22]46], we obtain a GW spectrum with an unpolarized lowfrequency tail and a chiral spectrum above the peak.

B. GW detection
Having predicted the GW spectrum generated by the relaxion dynamics, we can now evaluate its detectability in current and future experiments. A stochastic GW background can be detected in a given experiment if its signal-to-noise ratio (SNR) exceeds a threshold value thr . The SNR is given by [50] where T obs is the period of observation, (f min , f max ) is the frequency range of the detector, and Ω n (f ) is the detector's noise spectrum converted to fractional energy density. For a cross-correlated measurement in a network of detectors, as, for instance, a pulsar timing array (PTA), the noise spectrum has to be replaced by the effective noise Ω eff of the network (see Ref. [50] for further details), and the SNR is given by Eq. (38) with an additional factor of 2.
We present here the projected sensitivities of the planned space-based Laser Interferometer Space Antenna (LISA) [51,52], the planned Square-Kilometre Array (SKA) [53] PTA observatory, 6 and the proposed microhertz experiment µAres [55]. In addition, we evaluate current limits from the NANOGrav 11 year dataset [56]. Other GW observatories such as ground-based interferometers or potential LISA-successor experiments in the decihertz regime do not cover the frequency range in which a signal is expected in our scenario.
Recently, NANOGrav has further reported strong evidence for a common-spectrum stochastic process across the pulsars in their 12.5 year dataset [57], which might be due to a GW background. However, a GW detection could not be established yet, due to the lack of conclusive evidence regarding the interpulsar correlations of this process. Nonetheless, we also present here an interpretation of this signal as the GW background generated in our model, fitting our signal to the NANOGrav data based on the procedure outlined in Ref. [58].
Furthermore, GWs contribute to the total amount of radiation in the Universe and are therefore subject to constraints on N eff . This sets an upper bound on the peak amplitude of the spectrum. Using Eq. (23), we obtain where I S (ξ) = d log x S ξ (x) 0.8 -1.5 is the integral over the spectral shape and T rec is the photon temperature at recombination.  Fig. 4 An even stronger bound is obtained from the contribution of the dark photons to ∆N eff . Since the GWs are sourced by dark photon anisotropies, this directly leads to an upper bound on the GW amplitude which is stronger than the one from the direct contribution of the GWs to ∆N eff . Plugging Eq. (24) into Eq. (36), we obtain h 2 Ω peak GW 10 −8 g * ,ra 106.75 where the superscript X indicates that this is the dark photon contribution to ∆N eff . A similar constraint on the GW amplitude was already observed in Ref. [58] in the audible axion scenario. There, it was noted that the parameter region that provides the best fit to the stochastic GW background that was potentially observed by the NANOGrav collaboration [57] is in tension with constraints on ∆N eff .

IV. DISCUSSION
In Fig. 4 we show some example GW spectra. The corresponding parameter values are listed in Table I. The coupling to dark photons is set to r X = ξ/(2δ). In addition, the projected sensitivities of LISA (red), µAres (green), and SKA with an observation period of 5 years (turquoise) and 20 years (blue), represented by the corresponding power-law integrated sensitivity curves [50], are depicted as shaded regions. The current 95 % upper limits from the 11 year NANOGrav data set [56] are shown in purple. The example spectra B -D are detectable with µAres, whereas the benchmark A is accessible via SKA and is excluded by NANOGrav. The benchmark labeled " * " shown in orange corresponds to our best fit to the NANOGrav 12.5 year dataset as discussed further below.
At low frequencies, the GW spectra behave approximately as f 3 , in accordance with the expectations based on causality arguments using that the anisotropic stress of a causal source cannot be correlated at scales above the horizon size at the time of production [1,59]. At high frequencies, the spectra fall approximately like f −4 , allowing a simple distinction from the much steeper falling GW background generated from oscillating [20][21][22]46] or rotating [60] axionlike fields. It should further be noted that, when the peak position is fixed, changing ξ barely affects the UV tail, while the IR tail goes as ξ 2 (cf. Fig. 9), potentially allowing to one disentangle the reappearance temperature and ξ in the peak frequency Eq. (35) and thereby facilitating the determination of the relaxion parameters from a hypothetical observed signal. Larger values of ξ further result in a flatter peak, although this may be an artifact of our analytic approximation; cf. Eq. (37).
The range of peak frequencies and amplitudes that can be attained for coherently oscillating relaxion DM is displayed in Fig. 5, considering the case when the relaxion starts oscillating at the time the wiggles reappear. The polygons are obtained from the allowed range for the cutoff Λ, for T ra between approximately 240 MeV and v H ; cf. Fig. 2. The solid and dashed lines assume ξ = 100 and ξ = 10, respectively. Peak positions inside the polygons can be realized. The sensitivities of µAres, LISA and PTAs are again indicated as shaded regions. Note that the shading and solid lines correspond to the sensitivity for ξ = 100. For ξ = 10, the detection reach is degraded to the correspondingly coloured dashed lines.
For a large part of the peak frequencies and amplitudes that can be realized with relaxion DM, an observable signal is obtained, although mostly requiring futuristic space-based interferometers such as µAres for observation. For low values of ξ ∼ 10, a detection with SKA is possible. The present-day sensitivity of NANOGrav and expected reach of LISA are, however, insufficient for a detection. While NANOGrav is able to exclude parts of the parameter space if the DM assumption is relaxed (cf. benchmark A in Fig. 4), the sensitivity of LISA will remain insufficient even in this more general case.
We also show the GW sensitivity for relaxion DM as coloured regions in Fig. 3, using the same colouring scheme as above. The (light) coloured regions enclosed by the solid green and blue lines here indicate the relaxion masses and mixing angles for which, at least in a subrange of the reappearance temperatures in accordance with the constraints, an observable signal in µAres and SKA can be obtained. While µAres covers the full low-mass island as well as the range m φ 0.1 eV in the high-mass island, the sensitivity of SKA is limited to slightly lower DM masses. Note that, as the GW spectrum strongly depends on the reappearance temperature, a nonobservation in these experiments would in most cases not rule out the coloured parameter space, as a detection can be evaded by adjusting the reappearance temperature to shift the signal outside the experiment's reach. In the dark-green coloured region bounded by the dotted lines, however, the generated stochastic GW background is observable in µAres for the full range of allowed reappearance temperatures, guaranteeing a detectable signal in this region. The temperature dependence of the relaxion constraints and GW limits can be seen explicitly in the animated version of the figure that can be found in the ancillary material.  Table I. The orange spectrum labeled " * " corresponds to our best fit to the NANOGrav 12.5 year dataset. The colored regions depict the projected power-law integrated sensitivity curves of LISA (red), µAres (green), and SKA (turquoise and blue, assuming 5 and 20 years of observation, respectively) and the current exclusion from the NANOGrav 11 year data set (purple).  5. Values of the peak frequency and amplitude of the GW spectrum which can be obtained in the relaxion DM scenario. The edges of the polygon correspond to the minimal and maximal amplitudes which can be obtained for ξ = 100 (solid lines) and ξ = 10 (dashed lines), considering the case when the relaxion starts to oscillate immediately after barrier reappearance.
Last but not least, let us now dismiss the assumption that the relaxion constitutes DM and return to Fig. 1, where we again indicate the parameter regions in which the GW background can be observed in µAres (green), NANOGrav (purple), or SKA (blue/turquoise) by the respective colouring. The coloured regions respect the f φ < M Pl constraint. If we allow for super-Planckian decay constant, the regions extend to the dotted lines.
Regarding the GW sensitivity, the reader needs to be aware that the figure shows the projection of a fourdimensional plot, as T ra and f φ (or Λ) are not fixed. While red and orange shading marks the values of m φ and sin θ hφ for which it is not possible to evade the respective constraints by adjusting the remaining parameters (i.e. these coloured regions are definitely excluded), the GW contours (blue, turquoise, purple and green) correspond to the maximal reach of the respective experiments. They are obtained by taking the maximal SNR that can be achieved in each experiment for the given values of m φ and sin θ hφ . Given the experimental sensitivities we assume here, a detection via GWs can be evaded in all of the viable parameter space. In particular, the purple colouring and lines do not indicate that the corresponding parameter points are excluded by NANOGrav data but that NANOGrav constrains the reappearance temperature (and the decay constant) in this region. The same also applies to the projections for SKA and µAres. Furthermore, although the sensitivities overlap in the plot, µAres and PTAs operate in different frequency regimes and are therefore typically sensitive to different ranges of T ra . Nonetheless, we can conclude that current and future GW experiments can potentially detect the stochastic GWs from relaxion-stopping via dark photon emission, and thereby constrain the parameter space. We again enclose an animated version of the figure in the ancillary material available online, where the dependence on T ra (but not on Λ) is shown.
In addition to the current and prospective exclusion range, we also work out the parameter range in which our model can account for the potential GW signal recently observed in NANOGrav [57]. Fitting our spectrum to the data using the same procedure as in Ref. [58], where we keep ξ fixed and only fit the peak frequency and amplitude, we obtain the best-fit points and the corresponding 1σ and 2σ contours shown in Fig. 6. We further indicate the minimal peak frequency dictated by the lower bound on the reappearance temperature, T ra T BBN , as well as the maximal peak amplitude consistent with the constraint on ∆N eff by dotted lines. While these bounds exclude an explanation of the observed stochastic process in terms of our model for large values of ξ ∼ 100 (green), we can reach into the 1σ and 2σ regions for intermediate ξ ∼ 25 (orange), and for ξ ∼ 10 (blue), we can account for the NANOGrav best-fit point.
In the ξ = 10 (lower) panel of Fig. 1, we indicate the values of m φ and sin θ hφ for which we can attain the bestfit point (f peak = 3.3 nHz, h 2 Ω peak GW = 1.2 × 10 −9 ) by the grey shaded region. Note that, as we fix ξ and f peak , this also fixes the reappearance temperature to T ra ∼ 20 MeV, while the peak amplitude then fixes f φ as a function of m φ . Hence, in the grey shaded region, Λ can be adjusted within the constraints to obtain the respective value of the mixing angle. The best-fit spectrum for ξ = 10 is also depicted in orange in Fig. 4.

V. CONCLUSION
In this work, we have considered the possibility of probing the relaxion, which was proposed to ameliorate the Higgs hierarchy problem, via gravitational waves. A coupling to dark photons tames the relaxion excursion after reheating and thus in turn opens up a large fraction of the parameter space which was excluded in the minimal scenario without dark photons. Furthermore, dark pho-ton production after reheating can act as a source for the generation of a stochastic gravitational wave background. The gravitational waves are sourced by the anisotropies that are induced by the tachyonic production of dark photons between the electroweak phase transition and big bang nucleosynthesis. Hence, instead of the inflationary dynamics responsible for solving the hierarchy problem, we are probing here the late-time dynamics of the relaxion.
We have shown that this stochastic gravitational wave background can be probed by various current (NANOGrav) and future (SKA and µAres) gravitational wave detectors. In addition, we also highlight the parameter range in which our gravitational wave signal can account for the common-spectrum process observed in the most recent NANOGrav data [57]. Alongside the existing theoretical constraints, we have presented the relaxion parameter space which can be detected or excluded by the gravitational wave observatories in Fig. 1.
We find that the spectral shape of the gravitational wave signal in our model falls as the fourth power of the frequency above the peak, unlike a steeper falling gravitational wave signals generated by other axionlike field dynamics [20][21][22]60], whereas it behaves like f 3 in the infrared, as expected based on causality arguments. An observation of the spectrum in the range around the peak should allow for a determination of the reappearance temperature, while the amplitude can be used to determine the product of the relaxion mass and decay constant.
Furthermore, we have shown that the relaxion can constitute dark matter in the present Universe in the mass range of 10 −11 eV m φ 10 −6 eV and 10 −2 eV m φ 1 eV. While this scenario cannot be constrained with current NANOGrav data, most of the dark matter parameter space will be accessible by SKA and/or µAres in the future. Hence, with the advent of gravitational wave astronomy, we are now facing promising prospects for probing the relaxation of the electroweak scale via the stochastic gravitation wave background generated when stopping the relaxion, independently of whether the relaxion constitutes dark matter or not. Here, we discuss the minimal relaxion scenario where a coupling to dark photons is absent. As discussed in Sec. II, due to the restoration of the EW symmetry, the relaxion starts rolling again if the reheating temperature is higher than the EW scale, v H . To trap the relaxion in the minimum where it stopped during inflation, we need to require that it is displaced by less than ∆θ 2δ, where 2δ is the separation between the minimum and the maximum [7]. In the absence of a dark photon coupling, the Hubble friction is sufficient to trap the relaxion if m φ √ 8 δH ra , where H ra is the Hubble scale at the time of barrier reappearance [6]. In combination with the constraints discussed in Eqs. (19) to (21), this trapping condition limits the mixing angle as, within which the relaxion can be trapped by the Hubble friction. Moreover, the above equations restrict the relaxion mass for which the minimal scenario can be realized to m φ 2 × 10 −5 eV. In addition to these constraints, for the relaxion not to overclose the Universe we need to require Ω φ Ω DM . This further constrains the range of the relaxion mass, and thus the minimal relaxion scenario can only be realized for m φ 5 × 10 −8 eV .
(A4) Figure 7 shows the available parameter space in the relaxion mass m φ vs. mixing angle sin θ hφ plane. We assume here the minimal model without coupling to dark photons and do not require the relaxion to account for DM. Again, the red and orange shaded regions are excluded by the indicated constraints. The condition that the relaxion remains trapped in the same minimum in which it settled before reheating, ∆θ 2δ, then additionally excludes the region shaded in blue, with only the white region remaining as viable parameter space. In the turquoise shaded region, fulfilling ∆θ 2δ without dark photons further leads to an overclosure of the Universe. The region in which the relaxion can constitute DM is enclosed by the dashed black contour. If the relaxion is coupled to a dark photon, we further open up the blue and turquoise shaded regions of the parameter space. With a dark photon, relaxion DM can be realized within the solid black contour (same as the black framed regions of Fig. 3). Note that the dark photon scenario can only be applied in the region below the grey dashed line. Above the line, the condition H pp > H ra required for dark photon production is not satisfied (cf. the discussion in Sec. II C).

Appendix B: Dark photon spectrum
As discussed in Sec. II A, since θ > 0 in our case, positive-helicity dark photon modes with k < k ex = r X θ are tachyonic and experience exponential growth, whereas modes outside the tachyonic band oscillate. In particular, the amplitude of the oscillating modes does not change in time. At any time, the peak of the dark photon spectrum is given by the mode that experiences maximal growth with k = k m (τ ) = |r X θ |/2. We therefore take the ansatz where we neglect the negative-helicity modes as well as all modes that did not experience maximal growth yet, since these are exponentially suppressed. We treat modes with k ex < k < k m like the oscillating modes in order to include the peak of the spectrum in our estimate. Neglecting the energy produced before particle production, the dark photon energy density is given by where k pp = k ex (τ pp ) is the mode that exits the tachyonic band at particle production. On the other hand, we can determine the energy density spectrum of the dark photons from the relaxion energy density as where we assumed that at each time energy is dominantly transferred into the maximally growing mode, which goes like k m ∝ 1/τ . Hence, we can rewrite A k = A X k −9/2 with A X = π 2 Λ 2 br ξ r X a 2 ra k 2 ra , where k ra is the mode that exits the tachyonic band at reappearance, k ra = ξ/τ ra . The dark photon spectrum Available parameter space in the minimal relaxion scenario when the reheating temperature is higher than the EW scale. The red and orange shaded regions are excluded by the indicated constraints or combinations thereof. The blue shaded region is further excluded if the relaxion is not coupled to the dark photon source term, leading to a field excursion of more than 2δ. The dark photon production cannot occur in the region above the grey dashed line. Relaxion DM can be realized within the solid and dashed black contours, respectively, depending on whether a coupling to dark photons is assumed or not. In the turquoise parameter region, the condition ∆θ 2δ further leads to an overproduction of relaxion DM in the minimal scenario.
To corroborate our estimation, we have simulated the dark photon and relaxion evolution after reheating, solving the equations of motion (7) and (8) numerically, using 5000 logarithmically spaced, discretized momenta. In Fig. 8, we present the result of the numeric simulation for the dark photon modes as a function of the momentum k. The simulation assumes m φ = 1 meV, f φ = 2.35 × 10 13 GeV, Λ = 100 TeV, and ξ = 77. The value of r X has been determined numerically from Eq. (12) [6]. We show the spectrum at a/a pp = 4 (red), 20 (green), 100 (orange), and 500 (blue), where the full and light coloured lines correspond to the positive and negative helicity, respectively. As expected, the positive helicity modes dominate over the negative helicity by far. Furthermore, the amplitude for positive helicities indeed follows a k −9/2 law [cf. Eq. (B5), i.e.
√ 2k|X + | goes as k −4 ] between the momentum k pp that exits the tachyonic band at particle production and the peak momentum k m (τ ) that experiences the largest growth rate at the respective time. The peak momenta and k pp are indicated by the coloured and black, dashed vertical lines, respectively.