Searching for High Frequency Gravitational Waves with Phonons

The gravitational wave (GW) spectrum at frequencies above a kHz is a largely unexplored frontier. We show that detectors with sensitivity to single-phonon excitations in crystal targets can search for GWs with frequencies, $\mathrm{THz} \lesssim f \lesssim 100 \, \mathrm{THz}$, corresponding to the range of optical phonon energies, $\mathrm{meV} \lesssim \omega \lesssim 100 \, \mathrm{meV}$. Such detectors are already being built to search for light dark matter (DM), and therefore sensitivity to high-frequency GWs will be achieved as a byproduct. We begin by deriving the absorption rate of a general GW signal into single phonons. We then focus on carefully defining the detector sensitivity to monochromatic and chirp signals, and compute the detector sensitivity for many proposed light DM detection targets. The detector sensitivity is then compared to the signal strength of candidate high-frequency GW sources, e.g., superradiant annihilation and black hole inspiral, as well as other recent detector proposals in the $\mathrm{MHz} \lesssim f \lesssim 100 \, \mathrm{THz}$ frequency range. With a judicious choice of target materials, a collection of detectors could optimistically achieve sensitivities to monochromatic signals with $h_0 \sim 10^{-23} - 10^{-25}$ over $\mathrm{THz} \lesssim f \lesssim 100 \, \mathrm{THz}$.

The prediction of gravitational waves (GWs) by Einstein in 1916 sparked a century long search for their existence.
The first indirect evidence came from measurements of the orbital decay of the Hulse-Taylor pulsar [1], which were found consistent with GW emission [2].Direct evidence of the existence of GWs would follow in 2016 when GWs from an inspiraling pair of black holes (BHs) were measured [3].Even more recently the NANOGrav collaboration reported evidence of a stochastic GW background [4][5][6].These initial direct detections mark the beginning of "GW astronomy" as a viable method to study the universe.
In this work we show that high-frequency GWs can be detected via conversion into single phonons in crystal targets.
Phonons are quanta of lattice vibrations, and most crystal targets have gapped phonon modes, with energies in the 1 meV < ∼ ω < ∼ 100 meV range, corresponding to frequencies 10 12 Hz < ∼ f < ∼ 10 14 Hz.Because these modes are gapped, an incoming GW can be kinematically matched to the phonon dispersion relation, and therefore can resonantly vibrate the ions in the lattice, creating a phonon.This is somewhat analogous to an incoming photon converting to a phonon in a "polar" crystal, i.e., crystals with charged ions in the unit cell.If an incoming photon is kinematically matched to a gapped phonon mode in a polar target, the photon will resonantly drive dipole oscillations, which is equivalent to exciting the gapped phonon mode. 1he scarcity of high-frequency GW detectors is in part due to the difficulty of generating sufficiently strong GWs with astrophysical or cosmological sources.The generic problem of creating abundant high-frequency GWs is that the rapid oscillations necessary to generate high frequencies can only be achieved with astrophysically small masses, thereby limiting the signal strength.For example, this is precisely the difficulty in generating high-frequency GWs from the inspiral of two BHs.As we discuss in Sec.IV, the maximum BH mass, M max , which can generate GWs at a frequency ω is M max ∼ 10 −8 M ⊙ (THz/ω).This is a much smaller mass, corresponding to a much smaller signal, than the O(10) M ⊙ BHs LIGO is sensitive to.However, it is not physically impossible to generate detectable high-frequency GWs, and studying the sensitivity of current technology, likely built for other reasons, e.g., DM direct detection, provides excellent motivation for further studies of high-frequency GW sources.Indeed, perhaps the most important feature of our proposed detection scheme is that the detectors are already being built for DM direct detection experiments.Direct detection of DM by singlephonon excitations has been shown to be a promising route to search for both scattering of low-mass (sub-MeV) DM [67][68][69][70][71][72][73][74][75][76][77][78][79] and absorption of sub-eV bosonic DM [68,69,77,[80][81][82].The TESSARACT experiment [83] plans to utilize both Al 2 O 3 and GaAs as targets.Therefore, similar to the repurposing of axion DM detectors as GW detectors [47,61], any future DM direct detection experiment utilizing single-phonon excitations can directly be used as a GW detector.This paper is organized as follows.In Sec.II we derive the GW absorption rate into single phonons in two steps.
First, in Sec.II A, we summarize the general-relativistic derivation for the forces on a lattice of point masses due to an incoming GW.Then in Sec.II B, using the interaction Hamiltonian generated by the GW force, we derive the absorption rate into single-phonon excitations in both polar and non-polar targets.In Sec.III we begin by carefully defining the detector sensitivity to deterministic signals, i.e., those for which the GW strain is non-stochastic.We then discuss the detector sensitivity for a wide variety of target materials previously considered in the context of DM direct detection [73].In Sec.IV we discuss the general difficulty with generating high-frequency GWs, and then consider specific examples of potential high-frequency sources, such as superradiant annihilation IV A and BH inspiral IV B, carefully understanding their signal strength in the context of counting experiments.We conclude in Sec.V. A brief discussion of constraints on stochastic sources comprises Appendix A.

II. THEORETICAL FOUNDATION
The interaction of GWs with solid objects can be described in multiple ways.If the wavelength of the incoming GW is much larger than the detector, the deformation of the detector can be described using the theory of elasticity [84][85][86].
However, when the frequency of the GW is much larger, comparable to the O(meV − 100 meV) gapped phonons, the low energy effective theory of elasticity is no longer appropriate.The GW wavelengths, O(10 µm − mm) are smaller than the size of bulk single crystals.While much smaller than the size of the detecting crystal, these wavelengths are still much larger than the interatomic spacing, O(10 −10 m).Therefore, similar to how O(meV) photons couple to the dipole moment of a unit cell in polar targets to generate phonons, the effect of an incoming GW can be intuitively understood as a coupling to the quadrupole mass moment of the unit cell, which then generates phonons.Since the GW interaction is more similar to the photon interaction, computing the number of phonons produced from GWs proceeds more naturally by inheriting methods from particle and condensed matter physics.
The purpose of this section is to derive the absorption rate of incoming GWs into single-phonon excitations.We begin in Sec.II A with a review of the derivation of the force on a point mass within general relativity, and then generalize to a lattice of point masses.Furthermore we detail the necessary assumptions needed to treat the effect of the incoming GW on each point mass as a classical force.The potential energy associated with this force is the starting point of Sec.II B, and defines the interaction Hamiltonian.For non-polar targets, the focus of Sec.II B 1, the absorption rate then follows simply from Fermi's Golden Rule, and we derive the necessary matrix elements.In polar targets, Sec.II B 2, the phonons mix with the photon, complicating the absorption rate derivation.To compute the absorption rate for these targets we utilize the formalism from Refs.[82,87,88], which has carefully accounted for these mixing effects in the context of light DM absorption.We will work in natural units, h = c = 1, and use a mostly positive, (−, +, +, +), metric signature.

A. General Relativistic Forces
Our goal is to understand how a weak GW interacts with a crystal lattice, and therefore the natural starting point is understanding how a single mass interacts with a GW.While this is a textbook discussion [86,89], the general coordinate invariance of GR renders this a delicate subject.For example, in transverse-traceless (TT) coordinates the coordinate position of a single mass is unaffected by a passing GW.However this is simply an artifact of TT coordinates, reminding us that coordinates have no inherent meaning in GR; physics is encoded in coordinate invariants, e.g., proper distances.Fermi-normal (FN) coordinates equip the coordinates with physical meaning by defining the coordinates as the proper distance to an observer. 2ollowing Ref. [90], the equation of motion for a test mass in FN coordinates, assuming a non-relativistic mass and observer at the origin, is given by, to linear order in h, where the GW perturbation, h µν , is defined as a perturbative addition to the metric, g µν ≈ η µν + h µν , and η µν is the Minkowski metric.Note that h ij in defined in TT coordinates and the right hand side of Eq. ( 1) is coordinate independent to linear order in h: TT and FN coordinates are related by O(h) [90], and therefore a coordinate transformation would transform the right hand side of Eq. (1) from O(h) → O(h) + O(h 2 ).From Eq. (1) we see that a passing GW acts like a fictitious tidal force, F i , This result is straightforwardly generalized to a crystal lattice, i.e., each mass experiences a force, F i I = m I ḧij x j I /2, where I indexes the mass.
In addition to gravitational forces, the crystal lattice dynamics are of course influenced by the electromagnetic forces binding the crystal.These determine the leading order, O(h 0 ), contributions to the position and are simply the equilibrium positions of the ions (or, equivalently, atoms in non-polar targets), x i I ≈ x i I,0 .Perturbations around this equilibrium, u i I (t), induced by the GW are then driven by the force in Eq. ( 2), and back to equilibrium by the crystal harmonic "spring" forces.The equation of motion, at O(h), is then given by, where V IJ is the spring constant connecting mass I to mass J.Note that the V ij IJ only need to be determined to O(h 0 ).The perturbations themselves are O(h), and therefore, analogous to the driving GW force, a coordinate transformation would only generate shifts from O(h) to O(h) + O(h 2 ).Therefore the effect of a GW on a lattice of masses is to introduce a fictitious force, Eq. ( 2), which can drive the phonon modes in a crystal, analogous to a passing electromagnetic wave.

B. Absorption Rate Calculation
Since phonons are the quanta of lattice vibrations, the natural formalism to discuss single-phonon absorption rates is quantum mechanics.In calculating the absorption rate, the incoming GW may be treated as either a classical GW background or an incoming graviton [91], analogous to equivalence of describing a photon classically or quantum mechanically when considering light-matter interactions.We will choose the latter option for consistency with the phonon description.Lattice vibrations are described as oscillations around an equilibrium position; the position of the j th ion in the ℓ th unit cell of the crystal is where x 0 ℓj is the equilibrium position, r 0 ℓ is the lattice vector of the ℓ th unit cell, x 0 j is the equilibrium position relative to the center of the unit cell, and u ℓj (t) is the displacement of the ion away from equilibrium.When quantizing the lattice vibrations [92], u ℓj (t) can be expanded in terms of phonon raising and lowering operators: where N is the number of unit cells in the lattice, m j is the mass of the j th ion, ν indexes the band number, k indexes the first Brillouin zone (1BZ) momentum vectors, ω νk is the phonon energy, a νk and a † νk are the raising and lowering operators, respectively, which satisfy the canonical commutation relations, [a νk , a † ν ′ k ′ ] = δ νν ′ δ kk ′ ; and ϵ νjk are the phonon polarization vectors and satisfy ϵ νj−k = ϵ * νjk .The interaction Hamiltonian between phonons and a GW is determined by the coupling of u ℓj to h ij , which is due to the force in Eq. (3), This interaction Hamiltonian is the basis for the absorption rate calculation.At first sight one might simply wish to apply Fermi's Golden Rule to compute the graviton induced single-phonon transition rate.However, in polar targets, where the ions have charge, Q j , the phonon mixes with the photon, and screening can occur.Therefore we split the absorption rate in to two calculations.In Sec.II B 1 we consider non-polar materials, where Q j = 0, and Fermi's Golden Rule may be straightforwardly applied.In Sec.II B 2 we consider general materials, and compute the absorption rate using the optical theorem, following the derivation in Ref. [82].While the approach in Sec.II B 2 is more general, and reproduces the results in Sec.II B 1 in the limit of a non-polar material, it is also more technically involved.
Before continuing it is worthwhile to introduce notation common to both sections.It will be useful to work with the canonically-normalized GW field as opposed to the perturbations in the metric directly.This involves a rescaling, which transforms the interaction Hamiltonian in Eq. ( 6) to where h ij now has mass dimension one, and is quantized in the TT frame as, where ω p = |p|, and the overall normalization is found by equating the energy density, ρ GW = ⟨ ḣij ḣij ⟩/2, to the quantum mechanical Hamiltonian density, H/V = (1/V ) λp ω p a † λp a λp [93] with polarization vectors normalized to e λ ij e ij λ ′ = δ λλ ′ . 3Additionally we take the incoming GW to have four momentum, Q µ = (ω, q) = (ω, ω n), where ω = 2πf is the GW frequency and q = ω n is the GW momentum.
The absorption rate per incoming graviton of polarization λ is given by Γ λ .The total absorption rate, per detector exposure, for any kind of GW with polarization λ, is given by, where ρ T is the target mass density, and n λ GW is the number density of the λ th polarization of GWs.The differential energy density is related to the number density by dρ λ GW = 2πf dn λ GW .Assuming the incoming GW is isotropic and independent of polarization, where ρ c = 3H 2 0 /8πG is the critical density, H 0 is the Hubble constant today and Ω GW (f ) ≡ (1/ρ c ) dρ GW /d log f .Substituting Eq. ( 10) in to Eq. ( 9), and averaging over the GW polarizations to find the total, averaged, GW absorption rate, R, gives where we have defined Γ(f ) as the polarization and angularly averaged absorption rate per GW.Therefore to compute R in Eq. ( 11), we need to compute Γ λ , which will be the focus of Secs.II B 1 and II B 2.
For both polar and non-polar targets, in order to compute Γ λ the phonon energies ω νk , polarization eigenvectors ϵ νjk , and equilibrium ion positions x 0 j are needed.We calculate these with first-principles methods similar to Refs.[70,73,74,76,80,82], which we briefly summarize here.Using VASP [94][95][96] the ions and electrons in the lattice are relaxed to their minimum energy configuration, which sets the equilibrium positions, and perturbing the ions away from these minima allows calculation of the "force constants", V ℓℓ ′ jj ′ , which define the harmonic phonon Hamiltonian, Then V ℓℓ ′ jj ′ is diagonalized with the help of phonopy [97,98] to find the phonon eigensystem, and the final absorption rate, R in Eq. ( 11), is computed with the help of PhonoDark-abs [82].

Non-Polar Targets
Since the atoms at each site in non-polar targets are electrically neutral, there is no phonon-photon mixing and the absorption rate can be computed straightforwardly with Fermi's Golden Rule.Specifically, we will compute the transition rate due to the interaction Hamiltonian δH in Eq. ( 7) from an initial state, |I⟩ = |λ, f, n⟩ ⊗ |0⟩, containing zero phonons and an incoming graviton, to a final state, |F ⟩ = |0⟩ ⊗ |ν, k⟩, containing zero gravitons and one phonon indexed by its band number ν and crystal momentum vector k.This rate is given by where δH is evaluated at t = 0, since the time dependence has been factored out to provide the energy conserving delta function as in ordinary time-dependent perturbation theory, and ⟨I|I⟩ = ⟨F |F ⟩ = 1.The phonon contribution to the matrix element is computed using Eq. ( 5), where u ℓj ≡ u ℓj (t = 0), |ν, k⟩ = a † νk |0⟩, we have defined the phonon transition matrix element, T jνk , and used the long-wavelength approximation |k • x 0 j | ≪ 1, which is true for the momentum transfers of interest here since |k| ∼ ω ∼ meV and |x 0 j | ∼ Å ∼ (keV) −1 .The graviton contribution to ⟨F |δH|I⟩ is then given by evaluating the matrix element of ḧij (x 0 ℓj ) ≡ ḧij (x 0 ℓj , 0) in Eq. ( 8), using the quantization of the graviton field from Eq. ( 8), where |λ, f, n⟩ = a † λq |0⟩ and in the last step we made the same long-wavelength approximation |q • x 0 j | ≪ 1 as in Eq. ( 15).Substituting Eqs. ( 15) and ( 16) into ⟨F |δH|I⟩ gives where Ω = V /N is the volume of the unit cell.Eq. ( 17) can be simplified further since the phonon transition matrix elements satisfy the "coupling to mass" [68,69,71,80] sum rule, j m j T jνk = 0, when ν runs over the gapped, optical modes.Therefore, when expanding x 0 ℓj the term proportional to r 0 ℓ vanishes, and using ℓ e i(q−k)•r ℓ = N δ q,k further simplifies Eq. ( 17) to 4 Lastly, substituting Eq. ( 18) in to Eq. ( 14), we obtain where we have approximated ω νq ≈ ω ν , as appropriate for momentum transfers well within the 1BZ. 4 Note that the coupling to mass effect also enforces independence of the final result from the absolute position of the unit cell center.To see this, imagine choosing a new center, shifted by some ∆x, such that x 0 j → x 0 j + ∆x.The contribution to ⟨F |δH|I⟩ in Eq. ( 18) from ∆x then vanishes due to the coupling to mass effect.

Polar Targets
If the ions on the lattice sites have net electric charge, then they will also couple to the photon in addition to any incoming GW.This mixing introduces screening effects which frequently arise when studying absorption or scattering of light DM [68,81,82,87,88,99,100]. Accounting for this screening will require introducing some formalism which has appeared in Refs.[82,87,88].While technically more complex, this formalism will also apply to non-polar targets and provide a rigorous method for smearing the delta function in Eq. ( 19), in addition to being able to account for mixing effects in polar targets.
Consider an effective Lagrangian containing the graviton and photon fields, which are mixed by self-energies containing states in the medium, e.g., phonons or electrons, where L 0 h , L 0 A are the free graviton and photon Lagrangians in vacuum, respectively, and the second term in Eq. ( 20) contains the self-energies, Π, which are 1PI diagrams induced by the medium [82,87,88].For example, the singlephonon contribution to Π hh is given by the diagram, where the internal double solid line indicates a phonon propagator, the vertex Feynman rule is determined by the interaction Hamiltonian in Eq. ( 7), and we have left indicies implicit for simplicity.The single-phonon contribution to the other self-energies, e.g., Π hA , Π AA , are analogous diagrams with the corresponding external h replaced with A.
The in-medium states are those which diagonalize the Lagrangian in Eq. ( 20), and can be separated into unmixed "graviton-like" ( ĥ) and "photon-like" ( Â) states, since the couplings between the graviton and photon are perturbatively induced via the medium.The absorption rate of an incoming vacuum graviton is then, approximately, given by the absorption rate of its graviton-like counterpart, which can be computed with the optical theorem, where Π ĥĥ (f, n) is from the graviton-like term in the diagonalized Lagrangian, L ⊃ − λ ĥλ Π λ ĥĥ ĥλ /2.To compute the absorption we must diagonalize Eq. ( 20) to find Π λ ĥĥ (f, n).Gauge freedom allows for immediate simplification in both the photon and graviton sectors, reducing the degrees of freedom from 10 and 4 to 6 and 3 for the graviton and photon field, respectively.We are focused on computing the absorption rate of an incoming vacuum graviton.Since any polarization mixing introduced by the medium is gravitationally suppressed, we can reduce the graviton system to just the two vacuum graviton modes.Projecting the graviton field into the reduced polarization basis as h ij = e λ ij h λ , the Lagrangian in Eq. ( 20) simplifies to where λ, λ ′ index the usual +, × graviton polarizations, and self-energies with polarization indices indicate projection onto the polarization vectors, e.g., Π λ hh ≡ e λ ij Π ij,i ′ j ′ hh e λ i ′ j ′ .Moving to the photon sector, DM absorption calculations are typically performed in the Lorenz gauge, ∂ µ A µ = 0.However, subtleties arise in computing the Feynman diagrams in Eq. ( 23) for graviton kinematics, Q 2 = 0, since the standard choice of the longitudinal polarization vector, e µ L = (ω, ω n)/ −Q 2 , is singular.Therefore we work in the Coloumb gauge, ∂ i A i = 0, where the polarization vectors are e µ L = (1, 0, 0, 0), e µ ± = (0, n± ), and n± are two vectors orthonormal to n. Projecting the photon field into these polarizations, A µ = e λ µ A λ , the Lagrangian in Eq. ( 25) becomes, where we have explicitly separated the longitudinal photon mode A L from the transverse modes A σ , where σ indexes transverse polarizations, ±, and polarization indices represent projections onto polarizations.Projection on to the longitudinal mode is equivalent to projection on to ni by the Ward identity, The Lagrangian in Eq. ( 26) can now be perturbatively diagonalized since the off-diagonal components mixing h and A are O( √ G), and therefore Π λ ĥĥ is given by, where α, α ′ ∈ {L, ±} index all the photon polarizations, and the inverse in Eq. ( 27) represents the matrix inverse.
Note that for an isotropic medium with Π σ AL = Π σ LA = 0, ∆ AA is diagonal, and the inverse in Eq. ( 27) is trivial.The absorption rate is therefore given by To compute the self-energies in Eq. ( 29) we follow Ref. [82].As mentioned earlier, these self-energies will generally receive contributions from both the phonons and electrons in the target.However, since our focus is on single-phonon excitations, some simplifications can be made.Notice that the dependence of the absorption rate, Γ λ , on Π hh is only through the imaginary part.Imaginary contributions to self-energies correspond to physical degrees of freedom going on-shell, and therefore below the electronic band gap, Im [Π hh ] will only receive a contribution from phonon excitations.Furthermore, we expect the electronic contribution to Π hA to be sub-dominant, since any GW-electron coupling will be suppressed by m e /m j ∼ 10 −4 relative to the GW-ion coupling.Therefore when computing both Π hh and Π hA we include only the phonon contribution.Calculation of Π AA is identical to that described in Ref. [82], and will include both electron and phonon contributions.
The single-phonon contribution to the self-energies can be written as [82] Π where Π ΦΦ ′ implicitly carries the Lorentz indices of the fields Φ and Φ ′ , D ν is the phonon propagator, and F Φ,j is a form factor encapsulating the coupling of the j th ion to the field Φ (and also carries an implicit Lorentz index).
The photon form factor is simply the charge of a given ion, and the graviton form factor can be read off from the interaction Hamiltonian in Eq. ( 7), where Q j = N p,j − N e,j , and N p(e),j are the number of protons (electrons) at site j.Substituting Eqs. ( 32) and (33) into Eq.(30) gives which can then be straightforwardly substituted into Eq.( 29) to compute the absorption rate.Note that if the mixing is removed by setting Π hA = Π Ah = 0 (which can be accomplished by setting Q j = 0 as for a non-polar target), and the phonon is assumed to be a perfect resonance, γ → 0, then we recover the absorption rate previously derived from Fermi's Golden Rule in Sec.II B 1. The main effect of the mixing contribution to Eq.( 29) is to slightly shift the resonance locations from their value in the limit of no mixing, i.e., ω ν .

III. SENSITIVITY TO GRAVITATIONAL WAVES
With the calculation of the GW absorption rate in Sec.II, we can now determine the experimental sensitivity.In the context of DM direct detection, the sensitivity is governed by the number of phonons produced, N ph , over the observation time, T .Requiring that N ph > 3 is then directly converted to a 95% C.L. limit on the DM coupling parameters, assuming negligible backgrounds.While N ph will still be the figure of merit here, the question of experimental sensitivity becomes more subtle for GW searches since the incoming GWs need not be coherent over T , nor have a monochromatic frequency spectrum.These two features are dramatically different than the signal due to absorption of non-relativistic DM, which occurs at frequencies equal to the DM mass and persists over the whole observation time.Due to the relativistic kinematics of the incoming GW, a closer analogy would be to a signal produced by a thermal cosmic axion background [101], or from DM produced in the Sun [102,103].Here, computing the number of phonons produced becomes intimately tied to the parameterization of the signal.Therefore for ease of comparison it is useful to create classes of signals with similar parameterizations.The sources of high-frequency GWs in this frequency range we focus on, for reasons discussed further in Sec.IV, may be referred to as deterministic signals, characterized by a known signal profile which may be parameterized as,5 Stochastic sources, characterized by their power spectral density, are considered in App. A. In Eq. ( 36), f s (t) is the instantaneous signal frequency at time t, and h 0 (t) is the signal amplitude at time t.Assuming that f s and h 0 are slowly-varying functions of time, the GW energy density for deterministic signals is where we have defined, ρ 0 GW (f ) ≡ (π/8)f 2 /G.The GW energy density parameter Ω GW for deterministic signals is then given by Note that while Ω GW (f ) is typically used in the context of stochastic GW signals, its use here can be understood as just a rewriting of the differential energy density, dρ GW /df , given by Eq. ( 10); ρ GW in Eq. ( 37) is related to Ω GW in Eq. ( 38) by ρ GW = (dρ GW /df ) df = ρ c Ω GW (f ) d ln f .The expected number of phonons produced during T is found by integrating the rate in Eq. ( 11), and multiplying by the detector mass, M , While Eq. ( 39) is the underlying quantity which determines the sensitivity, comparing the response to different deterministic signals, for example those with different time dependence, is subtle.Ideally, the signal strength and detector sensitivity would be completely decoupled, such that if the signal strength is greater than the detector sensitivity, then the signal may be seen.This decomposition is most clear in the context of monochromatic signals, for which f s is constant, and the amplitude is slowly varying.For these signals, we can define, where h s characterizes the signal strength and h det characterizes the detector sensitivity.This allows Eq. ( 39) to be written as Therefore when h 2 s > h 2 det (f s ) the number of phonons produced is greater than 3, and limits on the signal can be set at the 95% C.L.
There is redundancy in our definitions of h s and h det (f ), since only their relative magnitude is important, and therefore great care should be taken to compare the sensitivity of different experiments which may differ in definitions of h s and h det .The definitions chosen here are useful because they do not only apply to detection via single-phonon excitations.For any direct detection experiment where the primary observable is some number of excitations, and the interaction of a GW with that excitation is Γ(f ), Eq. ( 41) can be used to define the detector sensitivity.
Another type of deterministic signal is a chirp signal, where ḟs = df s /dt is known.For these signals the time dependence in h 0 (t) can be traded for frequency dependence by solving for f = f s (t).Defining the signal strength in this case as where τ f ≡ f s / ḟs is the signal coherence time, the number of phonons produced from this signal type is can then be calculated, using dt = df / ḟ , as where h det has the same definition as Eq. ( 41), and the frequency integral is performed over f s (0) ≤ f ≤ f s (T ).We see that if h s > ∼ h det over an e-fold in f , then the number of phonons generated will be larger than 3.With Eq. ( 41), the detector sensitivity, h det , can now be computed in a general target material.In this section we focus on nine targets that have been previously studied in the context of DM direct detection.GaAs and Al 2 O 3 are being utilized in the TESSARACT experiment [83], specifically designed to search for DM-induced single-phonon excitations, and Si is used in SuperCDMS CPD [104] which is sensitive to secondary athermal phonons produced from a DM-induced event.Si and Ge have been used in other experiments, e.g., CDEX [105], DAMIC [106][107][108][109][110], EDELWEISS [111][112][113], SENSEI [114][115][116], and SuperCDMS [117][118][119], searching for DM-induced electronic excitations and nuclear recoils.NaI (DAMA/LIBRA [120], KIMS [121], ANAIS [122], SABRE [123], DM-Ice [124]), CsI (KIMS [125]), and CaWO 4 (CRESST [126][127][128]) have also been used as target material in nuclear recoil experiments.
In addition to these we include two more speculative targets: SiO 2 , shown to have strong phonon responses to absorption and scattering of dark photon DM [73,74,76,80,82], and diamond, which has been independently proposed to search for electronic excitations and nuclear recoils [72,129].
While it is likely not all of these targets will be employed as single-phonon detectors, comparing the projected sensitivity for a variety of targets illustrates how they can complement each other.The number of gapped phonon modes, and their energy spectrum, varies across targets.For example, CsI has a single resonance at ω ∼ 10 meV, diamond has one at ω ∼ 175 meV, and CaWO 4 has many resonances between 8 meV and 110 meV.The total number of gapped modes is 3n − 3, where n is the number of atoms in the unit cell, and the energy of the lowest gapped mode correlates with the speed of sound, which varies across materials.The GW sensitivity will be peaked near these resonance frequencies, and therefore a judicious choice of targets is necessary to cover the broadest possible frequency range.See Ref. [73] for the phonon band structures for the targets considered here.Additionally, in contrast to the search for dark photon DM [68,69,77,82], or axion DM in a background magnetic field [80,81] the targets do not have to be polar, i.e., contain oppositely charged ions in the unit cell.This allows for targets like diamond, Si, and Ge to be used to search for GWs even though they are not useful in searching for specific benchmark sub-eV DM models, which instead must rely on multiphonon processes [77] at O(meV) frequencies.
Understanding and mitigating backgrounds is crucial to the success of any direct detection experiment utilizing single-phonon excitations.The dominant irreducible background is from coherent solar neutrino scattering; the cosmic neutrino background is negligible.Solar neutrinos are expected to produce a background of phonons at a rate of R ∼ 10 −2 /meV/kg • yr between 1 meV < ∼ ω < ∼ 100 meV [130,131].Other backgrounds are expected to be important, and in fact dominate current low threshold direct detection technology.Examples of such backgrounds are those induced by cosmic high energy particles [132], and from radiogenic sources in the detector or shielding [131], the latter generating a background of phonons at a rate of R ∼ 1 − 10/meV/kg • yr with current levels of radio-purity.
Improved shielding, active signal vetoing, and more radio-pure samples will be essential to reduce these backgrounds.
Vibrational noise is unlikely to be important since any source frequency is far from the phonon resonances of interest here, and thermal noise will be Boltzmann suppressed.
In Fig. 1 we compare the detector sensitivity, h det , of the nine target materials assuming a kg • yr exposure and backgrounds at the irreducible level, < ∼ 1/kg • yr, which are negligible.The phonon line widths are taken to be γ ν = 10 −2 ω ν , values which have been shown to reproduce the dielectric function reasonably well [82].In addition   46), which emit GWs with power PGW, and are coherent over the observation time, i.e., monochromatic sources with hs ≈ h0 in Eq. ( 40).The red line, labelled "Superradiant Annihilation", is a monochromatic signal, from the annihilation of two bosons in a superradiant produced population around low-mass (M < ∼ 10 −7 M⊙) BHs.The signal strength is computed with Eq. ( 53), assuming the BH mass is maximal at each frequency; see Sec.IV A for more details.The blue lines, labelled "Black Hole Inspiral", are chirp signals from the inspiral of two low-mass (M < ∼ 10 −8 M⊙) BHs.The light blue dashed lines correspond to hs in Eq. ( 55) for fixed BH mass, M .The solid line represents the high-frequency boundary of possible BH inspiral signal strengths; at each f , we choose the BH mass M in Eq. ( 55) whose fISCO(M ) = f .to comparing the sensitivity of different targets, we also compare h det to the signal strength, h s , for the sources discussed in Sec.IV, all assumed to be r = 1 AU away.Gray lines are shown as benchmark lines and correspond to idealized, monochromatic sources which emit GW radiation from a single source isotropically with power P GW over the observation time.The sensitivity could be further improved by increasing the detector volume, h det ∝ 1/ √ V , as well as finding targets with smaller phonon linewidths, γ ν , since at the phonon peak, h det is proportional to √ γ ν .
In Fig. 2 we compare the sensitivity of single-phonon detectors to other high-frequency GW detectors in the MHz < ∼ f < ∼ 100 THz frequency range.The curve labelled "Phonons" is simply an outline of the sensitivities in Fig. 1.The detectors are compared in their ability to detect monochromatic GWs, coherent over all their observation times, with (time-independent) amplitude h 0 .Constraints from recasting limits from current experimental data Comparing the sensitivity of high-frequency GW detectors in the MHz < ∼ f < ∼ 100 THz range.The comparison between detectors is only meaningful for monochromatic signals with constant amplitude, i.e., those of the form h(t) ∼ h0 e 2πif t , where h0 is a constant.For these signals, if the h0 produced from a source is greater than the detector sensitivity, then it may be detected.The gray lines correspond to example, benchmark, sources isotropically emitting GWs with power PGW, at frequency f , and at a distance of r = 1 AU.h0 is computed using Eq.(46).The solid red line labelled "Phonons" is an outline of the specific target sensitivities shown in Fig. 1, which assume T = 1 yr.The sensitivity of different microwave cavities from Ref. [61] (ADMX [133][134][135], CAPP [136], HAYSTAC [137], ORGAN [138], and SQMS [61]) is shown in teal.The sensitivity of the MAGO 2.0 proposal from Ref. [59] are shown in green.The sensitivities of DMRadio-GUT [139] and DMRadio-m 3 [140], assuming a figure-8 pickup loop, have been combined, labelled "DMRadio", from Ref. [62] and are shown in orange.The sensitivity of atomic clocks are from Ref. [60] and shown in purple.

IV. SOURCES
The universe naturally provides a variety of GW sources.The Standard Model predicts GW production over a range of frequencies, from the inspiral of supermassive black hole binaries [148][149][150][151][152] at f ∼ nHz, to populations of gravitons frozen in from the early universe at f ∼ THz [153][154][155][156][157]. Additionally, many theories beyond the Standard  53).The turnover occurs when the superradiant annihilation timescale, Eq. ( 51), is equal to the observation time, T , and the cutoff is from requiring GM ω < 0.8, corresponding to the point where numerical simulations indicate a breakdown of our analytic approximations.The BH inspiral signal is shown at the different points in its frequency evolution and computed with Eq. ( 55).The cutoff occurs for BH masses too large to emit GWs at ω, i.e., their ωISCO < ω.
The characteristic strain of these cosmological, stochastic GWs is given by [46,194], which, while not an exact "apples to apples" comparison (see App.A for more details), is many orders of magnitude smaller than the detector sensitivities shown in Fig. 2. To emphasize how un-detectable this is, note that the detector sensitivity in Eq. ( 41) scales as L −3/2 , where L is the length scale of the detector.A detector would need to be scaled to L ∼ 100 km, while still maintaining zero background, to reach h c ∼ 10 −33 from its sensitivity of 10 −24 at L ∼ 10 cm.Therefore it will be challenging for detectors based on single-phonon sensitivity to detect these stochastic, early-universe GWs.
GWs may also be produced on astrophysical scales, i.e., by processes occurring in the Milky Way.Independent of how the GWs are generated, the GW amplitude h 0 from any source emitting GWs isotropically with power P GW at frequency ω and a distance r away is given by, Achieving a detectable strain at galactic distances requires an extraordinary amount of power.Said another way, if in one year the entire mass-energy of the Milky Way (∼ 10 12 M ⊙ ) were converted to meV-frequency GWs at the galactic center (r = 8 kpc), the amplitude would be only slightly above the detector sensitivity.
Therefore it appears the only reasonable source of GWs must come from a source at sub-kpc scales.Normalizing Eq. ( 46) to AU distance scales, the required P GW is much smaller, This signal strength is shown as a gray line in Fig. 2, and is more emblematic of the benchmark signals that very high-frequency GW experiments may be sensitive to. 6  We now consider two specific realizations of these deterministic signals: superradiant annihilation of bosons surrounding a BH (Sec.IV A), and BH inspiral (Sec.IV B).Their signal strengths, h s , are shown in Fig. 3 as a function of the BH mass.

A. Superradiant Annihilation
Superradiance is the process where a rotating BH of mass M creates a high occupancy bosonic cloud of particles by losing mass and angular momentum [195][196][197][198].The bosons that are produced occupy bound states around the BH, which are analogous to bound states in a hydrogen atom since the gravitational potential is approximately 1/r far enough from the BH.The production rate of bosons in a given bound state depends on the number of bosons present, and is therefore exponentially enhanced until the BH is spinning too slowly for superradiance to occur.The maximum number of superradiantly produced bosons in the "n = 2, ℓ = m = 1" mode, to borrow the notation of atomic physics, is, where M pl = 1.22 × 10 19 GeV, and ∆a * is the difference in the BH spin parameter, a * = J (GM 2 ) −1 , where J is the BH angular momentum, before and after superradiance occurs.These superradiantly produced bosons may then annihilate in to monochromatic GWs at frequencies ω = 2µ, where µ is the mass of the boson. 7For superradiance to occur, the BH must be spinning fast enough, mΩ H > µ, where Ω H is the angular velocity of the BH at the horizon.This is known as the "superradiance condition" and, when written in terms of the BH mass, requires GM ω < 1 for the ℓ = 1 mode.Therefore there is a maximum BH mass which can emit GWs at frequencies ω by superradiant annihilation, M sr max = M 2 pl /ω ∼ 10 −7 M ⊙ (meV/ω). 6Note that if compact objects of mass M constitute all of the local DM then the typical separation is (ρ DM /M ) −1/3 ∼ 100 AU M/10 −12 M ⊙ 1/3 , for ρ DM = 0.4 GeV/cm 3 .While the distance scale is much closer to AU scale of interest, the required masses are too small to generate meaningful GW signals. 7Level transitions also produce GWs, although at a subdominant rate to annihilation, and therefore we focus on the annihilation process here.
Computing the annihilation rate of two bosons in any state is a difficult problem because the gravitational "fine structure constant", α ≡ GM µ, can be non-perturbatively large.Therefore while analytic solutions exist for α ≪ 1, one must resort to numerical relativity simulations [199][200][201][202] for α < ∼ 1/2, near the edge of the superradiance condition.These simulations indicate that the maximum annihilation signal occurs for 0.25 < ∼ α < ∼ 0.5 [200,201], depending on the spin of the boson and the initial a * of the BH.For simplicity we adopt a semi-analytic approach to model the annihilation rate similar to Ref. [196], using the scaling relations from the perturbative regime, with an overall coefficient set by the numerical simulations, and restrict α < 0.4.With this parameterization the annihilation rate of two bosons in the n = 2, ℓ = m = 1 state is given by and where the overall coefficient, C Γ ≈ 10 −10 , is found from the numerical simulations in Ref. [200].
where the N 2 sr factor is due to the fact that the GW is emitted via an annihilation process.This generates P GW = ωR sr, ann of power in outgoing GWs.The annihilation process continues until all the bosons have annihilated to GWs.
The timescale for this to happen is also set by Eq. ( 50), The amplitude of the GWs emitted by superradiant annihilation is given by [196] h sr,ann 0 which can then be substituted in to the definition of the signal strength, h s , for monochromatic signals given in Eq. (40).The signal strength has parametrically different behavior depending on whether the signal lasts over the entire observation time: Eq. ( 53) is shown in Fig. 3 ("Superradiant Annihilation") as a function of BH mass, for different choices of signal frequency, assuming T = 1 yr.The maximum BH mass is set by M sr max .The turnover for each curve occurs when τ sr, ann ∼ T .BH masses above the turnover have τ sr, ann ≪ T , and those below have τ sr, ann ≫ T .Lastly, we note that while it may seem beneficial to increase the annihilation rate, Γ sr, ann in Eq. ( 49), to increase h s and GW power output, the timescale over which these GWs are emitted decreases, which hinders the signal.In fact, these competing effects exactly compensate each other.If the Γ sr, ann is increased enough to drive τ sr, ann ≪ T , then h s only depends on the combination R sr, ann τ sr, ann = N sr and becomes independent of Γ sr, ann .Therefore the peak h s , which occurs at large M when τ sr, ann ≪ T , is independent of Γ sr, ann .

B. Black Hole Inspiral
In addition to superradiant annihilation discussed in Sec.IV A, high-frequency GWs may also be produced from the inspiral of compact objects such as BHs, boson and fermion stars [203][204][205][206], gravitino stars [207], gravistars [208], and DM blobs [209].For simplicity we will examine the most straightforward scenario, the GW emission from the inspiral of two BHs with equal mass, M , in a perfectly circular orbit [86].To lose energy to GWs the orbital radius must decrease, drawing the BHs closer and thus increasing the frequency of the GWs produced.The minimum radius for which we can describe the inspiral as an adiabatically-changing circular orbit is, R ISCO = 12GM , where ISCO stands for "Innermost Stable Circular Orbit" (ISCO).Before the BHs reach the ISCO the GW frequency is given by, ω ∼ 10 meV 10 −9 M ⊙ /M (R ISCO /R) 3/2 .Requiring that R > R ISCO bounds M ω to be small.Therefore to reach a given ISCO frequency the BH mass must be less than M bhi max ∼ 10 −8 M ⊙ (meV/ω ISCO ), where "bhi" is shorthand for Black Hole Inspiral.Similar to the superradiant annihilation signal discussed in Sec.IV A, generating high-frequency signals comes at the cost of smaller masses.
A crucial difference between the superradiant annihilation signal and BH inspiral is that while superradiant annihilation is monochromatic, the GW frequency from BH inspiral is changing.This signal is therefore a chirp signal according to the classification discussed previously.The timescale of frequency change is given by, , which, while relatively fast compared to the observation time, is still large enough compared to the phonon lifetime (assuming ω ν /γ ν ∼ 100) to justify treating the signal as deterministic.Moreover, since τ bhi f ≪ T the minimum frequency, ω min = ω(t = 0) will be much less than ω max = ω ISCO , allowing us to approximate the signal bandwidth as 0 < ∼ ω < ∼ ω ISCO .Within this bandwidth the GW amplitude is [86], which, while seemingly large, is penalized in its overall signal strength (h s in Eq. ( 43)) by its small frequency coherence time, τ f .The signal strength is given by and shown in Fig. 3 (labelled "Black Hole Inspiral") as a function of BH mass for different choices of ω assuming T = 1 yr.

V. CONCLUSION
Single-phonon excitations, with energies in the O(1−100) meV range, have been shown to be exceptionally sensitive to light DM [67-71, 73, 74, 76, 77, 80-82].Here we show that this sensitivity extends to high-frequency GWs in the 10 12 Hz < ∼ f < ∼ 10 14 Hz range.In Sec.II we derived the forces acting on a lattice of point masses due to an incoming GW, and then used this to derive the GW-single phonon interaction Hamiltonian given in Eq. (6).We then used this interaction to derive the absorption rate in non-polar, Sec.II B 1, and polar, Sec.II B 2, crystal targets.An important difference between signals generated by absorption of light DM and GWs is that the GW signal may not be monochromatic or coherent over the observation time.This introduces subtleties in characterizing the detector sensitivity relative to the signal strength.In Sec.III we define the detector sensitivity, h det , and signal strength, h s , for the deterministic (as opposed to stochastic) signals of interest here.Most of the discussion readily generalizes to other direct detection experiments, making the definitions of h det and h s useful outside the context of just single-phonon excitations.We then computed the GW absorption rate into single-phonon excitations from first principles in nine targets: GaAs, Al 2 O 3 , SiO 2 , Si, Ge, Diamond, NaI, CsI, CaWO 4 , which have been well studied as targets in DM direct detection experiments.In Fig. 1 we show the corresponding detector sensitivity for each individual target, and in Fig. 2 we compare the detector sensitivities to other high-frequency GW experiments in the 10 6 Hz < ∼ f < ∼ 10 14 Hz frequency range.The diversity of targets highlights how a combination of targets may be used to provide broadband coverage in GW frequency space.
In Sec.IV we studied potential sources of these high-frequency GWs.GWs from the early universe are far too weak to be detectable with single-phonon detection, and an extraordinary amount of power must be emitted in GWs at O(kpc) distances to be measurable.Therefore it seems any viable source must be within the local solar system, at O(AU) distances.As examples of sources, we studied superradiant annihilation, Sec.IV A, and BH inspiral, Sec.IV B, and illustrated their relative signal strength in Fig. 3.In addition to these specific examples, in Figs. 1 and 2 we show the signal strength corresponding to any source emitting monochromatic, coherent, GWs with power P GW .These serve as useful benchmarks since they provide both a physically motivated comparison between experiments, and a target for future studies of high-frequency GW sources.
While our focus has been on single-phonon excitations in crystal targets, there are other ways GWs may interact with phonons.For single-phonon excitations, the lowest accessible frequency is set by the energy of the lowest gapped mode.However gapless, acoustic phonon modes have even smaller energies.These can be utilized in multi-phonon processes in both crystal [210], and liquid [211,212] targets where the incoming GW is kinematically matched to two phonons.It may also be interesting to utilize photon read out strategies as studied in Ref. [213].Furthermore, future research studying the details of different GW signals are important for differentiation from a DM induced signal.
The introduction of external electromagnetic fields may also be advantageous.The same inverse Gertsenshtein effect which the axion DM experiments are utilizing may also be an avenue to excite phonons in crystal targets or optomechanical experiments [214,215].Indeed, if an incoming GW converts to an electromagnetic field, that electromagnetic field can then excite optical phonon modes in polar materials, analogous to axion absorption on phonons in magnetized media [80,81].However in the crystal targets of interest here, the rate seems to be parameterically suppressed.The effective coupling parameter is h 0 B 0 , where B 0 is the external magnetic field.When compared to the axion DM coupling, g aγγ √ ρ DM B 0 /m a , this corresponds to h 0 ∼ g aγγ √ ρ DM /m a ∼ 10 −22 g aγγ / 10 −12 GeV −1 at ω ≈ m a = 10 meV, leading to much worse detector sensitivity than shown here.
We have shown that DM direct detection experiments utilizing single-phonon excitations can be powerful probes of high-frequency GWs in the THz < ∼ f < ∼ 100 THz frequency range.More broadly this suggests a connection between DM direct detection and high-frequency GW detection beyond axion experiments.This is important since the understanding of high-frequency GW sources is still emerging, and therefore multi-purpose DM-GW experiments present the best opportunity to explore the high frequency frontier.Furthermore, demonstrating sensitivity with multi-purpose DM-GW experiments motivates further research of GW sources which populate the high frequency landscape.
e r r a d ia n t A n n ih il a t io n B la c k H o le I n s p ir a l

FIG. 1 .
FIG.1.Detector sensitivities, h det , Eq.(41), for experiments utilizing gapped, single-phonon excitations in a variety of crystal targets, assuming negligible backgrounds and a kg • yr exposure.The frequency, and number, of resonances are properties of the specific target.All line widths, γν , Eq. (31), are taken to be 10 −2 ων , where ων is the phonon frequency.The projected sensitivity is cut off away from the peak resonances for visual simplicity, and to emphasize that the simplistic modelling of γν discussed previously is unlikely to hold far from the resonances.Note that Si, Ge, and Diamond are non-polar, whereas the rest of the targets are polar.In addition to the detector sensitivity we show the signal strength, hs, for three sources of high-frequency GWs, all assuming sources r = 1 AU away and T = 1 yr.Gray lines represent idealized, benchmark sources, computed with Eq. (46), which emit GWs with power PGW, and are coherent over the observation time, i.e., monochromatic FIG.3.Comparison of the signal strength, hs, from the inspiral of BHs with equal mass M ("Black Hole Inspiral", blue, Sec.IV B), and superradiant annihilation of bosons surrounding a BH of mass M ("Superradiant Annihilation", red, Sec.IV A) at different GW frequencies, ω.All lines assume the source is r = 1 AU away and T = 1 yr.The superradiant annihilation signal is computed with Eq. (53).The turnover occurs when the superradiant annihilation timescale, Eq. (51), is equal to the The total annihilation rate is further enhanced by the number of bosons in the cloud, 15,