A novel search for high-frequency gravitational waves with low-mass axion haloscopes

Gravitational waves (GWs) generate oscillating electromagnetic effects in the vicinity of external electric and magnetic fields. We discuss this phenomenon with a particular focus on reinterpreting the results of axion haloscopes based on lumped-element detectors, which probe GWs in the 100 kHz-100 MHz range. Measurements from ABRACADABRA and SHAFT already place bounds on GWs, although the present strain sensitivity is weak. However, we demonstrate that the sensitivity scaling with the volume of such instruments is significant - faster than for axions - and so rapid progress will be made in the future. With no modifications, DMRadio-m$^3$ will have a GW strain sensitivity of $h \sim 10^{-20}$ at 200 MHz. A simple modification of the pickup loop used to readout the induced magnetic flux can parametrically enhance the GW sensitivity, particularly at lower frequencies.

The present gravitational wave (GW) program is focussed on the nHz to kHz frequency range, motivated by the signals expected from the merging of known compact astrophysical objects. This focus leaves the ultra-high frequency (UHF) range, above a kHz, largely unexplored, despite its unique opportunity to probe the physics of the very early Universe. For a recent summary of the challenges and opportunities at high frequencies, see Ref. [1].
In this work we propose a novel strategy for the UHF range based on GW electrodynamics: the modified version of electromagnetism applicable in the spacetime metric of a GW. We show there exists a close analogy to axion electrodynamics -the appropriate formalism when working in the background of an ultralight axion -and exploit this connection to convert axion haloscopes into GW telescopes. In the vicinity of static electric and magnetic fields, a GW sources a small electromagnetic signal oscillating at the GW frequency. This opens up the possibility of using low-mass axion haloscopes with a frequency range of 100 kHz -100 MHz as detectors for UHF GWs. Our proposed search strategy will utilize the anticipated rapid progress being made in this field by lumpedelement axion detectors. We will show how the existing results of ABRACADABRA [3][4][5][6] and SHAFT [7] can already be recast as novel limits on GWs. Going forward, DMRadio [8][9][10][11] will dramatically extend both the axion and GW reach at these frequencies. As we will demonstrate, the reach of DMRadio can be enhanced with a simple modification to the signal readout, by using a semicircular "figure-8" loop to measure the magnetic flux. With this adjustment, the future reach of DMRadio, as shown in Fig. 1, will represent both a competitive and complementary approach to UHF GWs. Looking even further into the future, the scaling of the GW reach with the instrument volume is particularly advantageous, which we highlight with the larger DMR 8 -100. Figure 1 also depicts existing proposals for the UHF band, including optically levitated sensors [12,13], bulk acoustic wave (BAW) devices [14] (see also Ref. [15]), interferometers such as the holometer [16][17][18], current and future microwave cavity instruments, such as the axion FIG. 1. The UHF-GW experimental landscape, with the approach introduced in this work shown in color. DMRadio 8 shows the projected reach of the full suite of DMRadio instruments (50L, m 3 , and GUT) adopting our advocated figure-8 pickup loop geometry. Looking to the far future, we also show the reach of an upscaled DMRadio with a magnetic field volume of 100 m 3 , labelled DMR 8 -100. A subset of existing proposals in this frequency range are shown in grey, taken from Refs. [1,2], as well as an estimate for the required sensitivity to see one signal from primordial black hole (PBH) binaries per year. Additional specifics are provided in the text. haloscope approach introduced in Ref. [2] for ADMX [19][20][21] and SQMS (see also Ref. [22]), as well as cosmological probes of GWs [23] based on observations by the radio telescopes EDGES and ARCADE [24,25]. We emphasize that the different proposals we show are at varied levels of maturity, and therefore caution against overly quantitative comparisons. Instead, we refer to the specific references for details. Our proposal is related to the (inverse) Gertsenshtein effect [26,27], which describes the conversion of GWs into photons, a concept that has been proposed as method to search for GWs, in different frequency regimes, in the laboratory [27][28][29][30] and cosmology [31][32][33][34][35][36]. For electromagnetic GW detectors in a broader sense, see also pioneering work in Refs. [37][38][39][40][41][42][43].

arXiv:2202.00695v2 [hep-ph] 3 May 2022
We organize the discussion as follows. We begin with several general comments on the modification to electrodynamics induced by a passing GW, before specializing to the case of interest: the sensitivity of instruments that use a toroidal magnetic field. We then outline how we can exploit this to recast existing ABRA and SHAFT results, and future DMRadio searches. In the Supplementary Material (SM) we provide the full details of our calculations and a brief discussion of GW sources in the UHF band. Gravitational Wave Electrodynamics. We describe the spacetime in the presence of a GW by the linearized metric g µν = η µν + h µν , with |h µν | 1. The perturbation to the flat-space metric generates a correction to the kinetic term of electromagnetism, which can be written as the following effective current, Here F µν is the electromagnetic field-strength tensor, and as we demonstrate in the SM, the effective polarization and magnetization vectors are where h = h µ µ . Manifestly, near large external electric or magnetic fields, GWs will source oscillating fields, the detection of which is the focus of this work.
The above formalism facilitates a comparison with axion electrodynamics. In particular, the coupling between the axion, a, and electromagnetism is also described by Eq. (1), but with P = g aγγ a B and M = g aγγ a E, as follows from L ⊃ g aγγ a E · B [44][45][46], where g aγγ is the axion-photon coupling.
We can extend this analogy in order to estimate the expected GW sensitivity of axion haloscopes. For both the axion and GW, the magnitude of the induced fields is controlled by a dimensionless combination, either the strain h ∼ |h µν | or g aγγ a. The axion darkmatter program aims to probe the QCD axion, for which m a f a ∼ m π f π , where f a is the axion decay constant, with g aγγ = α/2πf a . Matching the strain sensitivity to the average g aγγ a ∼ g aγγ √ ρ DM /m a for the QCD axion, we This estimate sets the scale for the GWs that can be potentially detected (see Fig. 1). In this argument we introduced the electromagnetic fine structure constant α, the local dark-matter density ρ DM , as well as the pion mass and decay constant, m π and f π . Cosmological GW sources, which are typically isotropic and incoherent, are bounded by constraints on the total amount of radiation in the Universe to satisfy h 10 −29 (100 MHz/f ) ∆N 1/2 eff , well below our estimated reach. More promising search targets are rare exotic astrophysical GW sources. As a concrete example, if primordial black holes (PBH) with m PBH M contribute to the dark-matter density, then some fraction of these will exist in binaries and emit high frequency GWs through their inspiral and eventual merger phase. In order to estimate the size of this signal, we take the most up-to-date estimates of the fraction of PBHs in binaries and their expected merger rate (see Sec. 4.6 of Ref. [47] for a recent review, although we emphasize that uncertainties remain such as the impact of accretion on the merger rate). Combining the merger rate with the assumption that PBHs saturate the relevant microlensing constraints [48], and accounting for the local overdensity of binaries given by the Milky Way halo, we arrive at the estimated sensitivity shown in Fig. 1 required in order to see one event per year at that frequency (marginalizing over m PBH ). See the SM for further details of this calculation and a discussion of other putative signals. We thus primarily focus on localized, approximately coherent GW signals in the following.
In the transverse-traceless (TT) gauge, the nonvanishing components of a plane GW read where φ h and θ h are, respectively, the azimuthal and inclination angles of the GW, and throughout we employ the shorthand s α = sin α and c α = cos α. h + and h × are the strain amplitude associated with the plus and cross polarizations. The choice of V is a convention, any unit vector perpendicular to k can be adopted. We note that different choices will mix the definitions of h + and h × .
In the TT gauge, the form of Eq.
(3) appears to considerably simplify Eq. (2), although for our purposes this can be deceptive because the experimentally generated electric and magnetic fields in the equation are not naturally defined in the TT frame. Instead, throughout this work we will operate exclusively in the frame where all detector quantities are defined. This is the proper detector frame, in which following Ref. [2] (see also Refs. [49][50][51]), we find . See the SM for details. On very general grounds, Eq. (4) shows that the effective current and therefore the fields it induces are rapidly suppressed for frequencies below the inverse length scale of the instrument. As we show, there are further suppressions if the experiment is highly symmetric.
Application to a Toroidal Magnetic Field. We consider now an explicit experimental configuration to detect the oscillating fields sourced by a passing GW. In particular, we follow the original ABRA proposal The detector geometry we consider for the detection of GWs. The experimental setup follows ABRA [3]: a toroidal magnetic field B 0 , as given in Eq. (5), is generated inside a toroid of inner radius R, width a, and height H. In the presence of a GW or axion, a magnetic flux is generated in a pickup loop placed at the center of the toroid, where B 0 = 0. Axion detection makes use of a circular pickup loop of radius r, whereas for optimal GW detection we advocate the figure-8 configuration depicted, made of two oppositely oriented semicircles. establishing a DC toroidal magnetic field, and zero otherwise. The corresponding geometry is depicted in Fig. 2. The combined effect of the magnetic field and a GW is the effective current in Eq. (1) that, according to the Biot-Savart law, sources a magnetic field in the region ρ < R such that where j ρ = j eff ·ê ρ and j φ = j eff ·ê φ . As we review in the SM, corrections to the Biot-Savart law from the displacement current enter only at O(ω 4 ). To detect this magnetic field, a pickup loop is placed at the center of the toroid, which will be sensitive to a magnetic flux equal to Eq. (6) integrated over the area of the loop. We next consider two different loop geometries, beginning with the approach used in existing axion instruments.
Circular pickup loop: For a complete circle, the j ρ contribution vanishes -independent of the form of j µ effleaving A shift of φ → φ + φ removes the φ dependence in the integrand in all terms except j φ . Using the results in the previous section, we obtain Independent of the incident GW direction, the h + component decouples. Furthermore, in contrast to what would be naively expected from Eq. (4), the flux receives no contribution at O(ω 2 ) because the Bessel function of second order satisfies J 2 (x) = x 2 /8 + O(x 3 ). Instead, in the limit R H 1/ω the leading contribution to the flux is Nonetheless, a GW will generate a flux, and existing axion searches for such a flux can constrain UHF-GWs. In this regard, Eq. (9) may be compared against the equivalent quantity for axions [3], for which j ρ = 0, j φ = g aγγ (∂ t a)B max R/ρ, and expanding in R/H yields To isolate the fate of the expected ω 2 contributions, let us note that, at leading order in ω, where we take φ h = 0. (It can be restored by replacing φ → φ − φ h .) All of these terms will vanish for a circular pickup loop geometry: as noted above j ρ cannot contribute in this case in general, and from the explicit expression we see that at O(ω 2 ), 2π 0 dφ j φ = 0. Having identified this, however, we can see that with an alternative geometry, these leading terms will survive. To be explicit, using the above currents we can compute the leading contribution to the magnetic field in the plane at the vertical center of the toroid, from Eq. (6), The sinusoidal variation of both polarizations demonstrates that the maximum flux is achieved with a pickup loop with oppositely oriented semicircles for φ ∈ [0, π) and φ ∈ [π, 2π)-the "figure-8" shown in Fig. 2.
The "figure-8" configuration: Integrating Eq. (12) over the figure-8 configuration yields The result is now O(ω 2 ) and sensitive to both polarizations. While it will maximize the GW sensitivity, the figure-8 pickup loop is insensitive to the axion signal, as the B z generated by the latter is independent of φ . A single semicircle is sensitive to both, with fluxes given by half of Eqs. (10) and (13), respectively. Comparing the two expressions, we also see that to compensate the ω 2 factor, the GW flux scales with an extra power of r, indicating the improved volume scaling over the axion flux. Accounting for the minimal inductance of the pickup loop [3], we see Φ 8 ∝ V 7/6 whereas Φ a ∝ V 5/6 . Ultimately, the beneficial volume scaling can be traced back to the fact that our measurement is linear in the induced fields, which themselves are proportional to h for the GW or a for the axion.
From Eq. (13), we note that the response to the two polarizations differs and depends on the direction of the incoming GW. With two identical detectors angled appropriately, polarization measurements as well as sky localization become possible (see also Ref. [52,53]). For a sufficiently coherent source, one may hope to use the Earth's rotation or even a mechanical rotation of the experimental setup to break these degeneracies with a single detector.
Gravitational Wave Sensitivity. Low-mass axion haloscopes perform a search for anomalous magnetic flux, and in the absence of a significant signal above background, interpret the results as limits on g aγγ through the use of Eq. (10). With the same equation we can convert existing and projected limits on g aγγ into limits on Φ a , which we can recast as strain sensitivities when compared with our predictions for the GW flux. The procedure is not quite as simple as equating the latter to Φ a , however. The sensitivity also depends on the relative coherence time of the two signals-a longer coherence time corresponds to a narrower signal in the frequency domain, which in general can be more sensitively detected over the background (for an extended discussion, see Ref. [54]). The coherence time of a signal with mean frequencyω is τ ∼ 2πQ/ω, where Q is the quality of the signal, a dimensionless measure of the inverse width of the frequency distribution. The non-relativistic nature of axion darkmatter implies a highly coherent signal with Q a ∼ 10 6 , such that τ a ∼ (1 neV/m a ) µs. Considering experimental runtimes longer than τ , then the flux sensitivity will scale as Φ ∝ Q 1/4 [55], and so our actual limit on the GW flux is given by Φ a (Q a /Q h ) 1/4 . Beyond this we assume the signal persists during the relevant experimental runtime, and fix the incident GW direction ask =ê y .
What remains is to determine a value for Q h . This will

DMRadio GW Sensitivity
The GW strain sensitivity of low-mass axion haloscopes. We recast the existing limits obtained by ABRA [6] (green) and SHAFT [7] (purple). For DMRadio we use the projected future sensitivity of the three instruments that will make up that program: 50L (blue), m 3 (cyan), and GUT (pink) [10,11]. In each case, results are shown for two choices of the GW signal coherence, Q h = 1 (opaque) and Q h = 10 depend on the specific source. As mentioned already, the localized sources that are our focus can be coherent, but are not expected to be as extremely coherent as a darkmatter signal (see the SM for further discussion). As benchmarks, we therefore consider sensitivity to both a coherent Q h = 10 3 and incoherent Q h = 1 signal, to indicate the dependence upon this choice.
With these choices, in Fig. 3 we show the reach of existing and future haloscopes, taking the default circular pickup loop geometry. This roughly amounts to using Eqs. (9) and (10), except that in all figures we use the full flux expressions (rather than these leading order results), which are provided in the SM. For ABRA-10 cm [6] and SHAFT [7], we recompute the GW flux for the explicit geometries specified in those works, as well as accounting for the ferromagnetic core adopted by SHAFT. The future projections of DMRadio-50L, m 3 , and GUT are obtained by scaling up the toroidal geometry of the ABRA-10 cm instrument to the volumes of these future experiments specified in Refs. [10,11]. As the exact parameters of these instruments have yet to be specified, these results should be interpreted as representative of the parametric reach-the sensitivities will change by O(1) amounts once the geometry of the instruments is known. For instance, while the 50L instrument will adopt a toroidal geometry, m 3 will be solenoidal, and the geometry of DMRadio-GUT has not yet been specified. In all cases, there is a suppression in the sensitivity at lower frequency, which is a result of the ω 3 factor in Eq. (9).
To overcome this, in Fig. 1 we show the combined sensitivity if the DMRadio instruments adopted the figure-8 configuration (assuming Q h = 10 3 ). While such DMRadio results are still more than a decade away, in order to highlight the significant volume scaling of our approach, we also show the reach of a scaled up version of DMRadio-GUT. In particular, DMR 8 -100 is an instrument with the same toroidal geometry as above, but scaled up to a magnetic field volume of 100 m 3 , the largest instrument suggested in the original ABRA proposal [3]. (The magnetic field volume, V B , is also referred to as the effective volume [7] or denoted by GV , where G is the geometric coupling V is the volume of the toroid [4][5][6].) Conclusions. We provide a formulation of GW electrodynamics which demonstrates that low-mass axion haloscopes are also UHF-GW telescopes. Through the use of an optimized figure-8 pickup loop geometry, the DMRadio program may discover not only the dark matter of our Universe, but also exotic sources of GWs.
Going forward, there are several questions opened by our results that warrant further study. Determining the optimal experimental configuration that allows DMRadio to unlock the ω 2 GW scaling, while maintaining full axion sensitivity will be critical. To that end, including both a circular and figure-8 loop in the detector may be required, and we note that the mutual inductance between the two loops would vanish [56]. Further, lumpedelement instruments have at present only performed measurements for ω 1/R. Given the ω 2 scaling of the GW effective current, experimental results at higher frequencies would be particularly welcome. From a theoretical perspective, it would also be interesting to understand the limit where ω ∼ 1/R, where a calculation based on the Biot-Savart law alone is insufficient. Nevertheless, the exact behavior has yet to be quantified even for axions, and would be particularly interesting to consider for GWs given the ω 2 scaling of the effective current. Un-derstanding the extension to other geometries will be important, particularly the solenoidal configuration already used by ADMX SLIC [57] and that will be adopted by DMRadio-m 3 . Finally, the discussion we presented has been couched purely in terms of the raw GW strain sensitivity. However, as our understanding of the sources at UHF improves, it will be necessary to reinterpret these results in terms of the specific model parameters describing the source populations. At such a time, it would also be worthwhile to understand the additional information that multiple instruments could uncover, as we briefly discussed.
While these issues should be resolved, a larger question looms. Figure 1 highlights that a number of distinct experimental proposals have coalesced on a strain sensitivity of h ∼ 10 −22 for f ∼ MHz GWs, a level that is still many order of magnitude away from any signal of the early Universe. Whether we can hope to probe such strain sensitivities remains to be determined. and Jan Schuette-Engel for helpful discussions. Additionally, we thank Jonathan Ouellet and the authors of Ref. [2] for comments on the manuscript. C.G.C is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy EXC 2121 "Quantum Universe" -390833306 and by the Alexander von Humboldt Foundation.
Note added. For an extended discussion of the GW signal from PBH binaries including also the stochastic component see Ref. [58], which presents results obtained in parallel and independently of the results shown here. Our results are consistent with their findings.

SUPPLEMENTARY MATERIAL "A novel search for high-frequency gravitational waves with low-mass axion haloscopes"
Valerie Domcke, Camilo Garcia-Cely, and Nicholas L. Rodd In this Supplementary Material we provide detailed derivations of the results quoted in the main text, including somewhat lengthy but useful analytical expressions which can be readily used to analyze more general detector geometries. We begin in Sec. I by deriving the expression for the GW in the proper detector frame, Eq. (4) from the main text, before turning to the heart of our calculations, the effective current induced by a passing GW in Sec. II. Finally, in Sec. III we review possible GW sources in the UHF band.

I. GRAVITATIONAL WAVES IN THE PROPER DETECTOR FRAME
The two most commonly used prescriptions for fixing the gauge redundancies in linearized general relativity are the TT gauge and the proper detector frame. In the former, the coordinate system is defined by freely falling test masses. The GW tensor takes a particularly simple form obeying This simplicity comes at a cost. In the TT frame, the description of rigid bodies including the various aspects of the detector, become relatively unintuitive, as their coordinates are deformed by a passing GW due to the motion of the coordinate system. On the other hand, in the proper detector frame the coordinate system is defined by rigid rulers and closely matches the intuitive description of an Earth-based laboratory, with the GW acting as a Newtonian force. The challenge associated with this frame is that the form of h µν will be more involved, particularly when the wavelength of the GW is comparable to the size of the detector. Of the two possible complications, we find the latter preferable, and therefore work in the proper detector frame throughout. We refer the reader to section 1.1.3 of Ref. [59] for a comprehensive discussion of both frames, and to Ref. [2] for specific examples in the context of resonant cavity axion haloscopes. In this paper, we assume the entire experimental apparatus to be rigid. The proper detector frame can be implemented with the use of Fermi Normal coordinates [49][50][51]. Doing so, the GW tensor is given by n + 2 (n + 3)!R 0kil,m 1 ...m n r k r l r m 1 ...r m n , with R µνρσ denoting the Riemann tensor. The m i indices appearing after the comma indicate that a spatial derivative with respect to the direction m i must be taken. The notationR indicates that the Riemann tensor and its derivatives are to be evaluated at a specific reference point, which we take to be the origin of the coordinate system. To linear order in h, the Riemann tensor is gauge invariant and can thus be evaluated in the particularly simple TT gauge. Accordingly, we compute Riemann tensor for a GW using k µ = (ω, k) and the conventions we established in (3). This procedure leads to Furthermore, for a GW we have R µνρσ,m 1 ...m n r m 1 ...r m n = (ik · r) n R µνρσ . Using this and the auxiliary function we can resum the series in Eq. (S2). We find which leads to Eq. (4) in the main text. These expressions agree with those previously shown in Ref. [2] for the particular casek =ê z (noting the different sign convention for plane waves adopted in that work). Corresponding results for arbitrary incident angles were moreover used in the numerical evaluations performed in Ref. [2]. From the results above, it is clear that in the proper detector frame, at leading order h µν ∝ ω 2 . As the experimentally measurable flux is linear in h µν , this implies the best scaling we can hope to achieve with frequency is also Φ ∝ ω 2 , and indeed this is attained by the figure-8 pickup loop.

II. THE EFFECTIVE CURRENT
The experimentally measurable effect discussed in the main body originates from the effective current a GW generates in the presence of electromagnetic fields. In this section we first discuss this current in general, introducing the language of effective polarization and magnetization for GW electrodynamics and demonstrating the close analogy with axion electrodynamics. We then show how a GW generates an effective current, giving explicit analytical expressions for all relevant components. We conclude this section by turning to the integration of this current, deriving the induced magnetic field and the resulting flux through the pickup loop.
Before we begin, let us be explicit about the order to which we compute our results. In the main text, we showed analytic expressions to O(ω 3 ) and expanded assuming (R + a)/H 1, in terms of the geometry in Fig. 2. All sensitivity estimates shown, however, did not expand in the height of the toroid. Here, using the results of Eq. (S5) we will compute the effective current to all orders in ω, and with no assumption as to the size of H. However, as we will review below, that we compute the induced magnetic field through the Biot-Savart law implies that we will not be able to determine B z -and consequently Φ -to O(ω 4 ) or higher. As lumped-element axion haloscopes operate mostly in the regime where ω 1/R, this will be sufficient for our purposes; deriving the full result valid even for ω ∼ 1/R would be an interesting future direction. (We note that, in contrast, the results in Ref. [2] do extend to all orders in ω, as required for resonant cavity instruments.) Unless otherwise stated, all results outside the main text will not assume (R + a)/H 1, and it is those results we used to compute all our projected limits.

A. The Effective Magnetization-Polarization Tensor
As we will show below, in the absence of ordinary electromagnetic currents, the effect of axions or GWs on electromagnetic fields can be described by a conserved effective current The conservation of this current motivates the introduction of a skew tensor where M νµ eff is the effective magnetization-polarization tensor. Its six components can be expressed in terms of two vectors, These are the effective polarization and magnetization, respectively. In terms of these two vectors, the effective current then becomes as in Eq. (1). From this, we can identify an effective charge ρ eff = −∇·P and an effective 3-current j eff = ∇×M+∂ t P.
Maxwell's equations take the familiar form (S10) Both axions and GWs can be formally conceived as a continuous medium [44][45][46]60], as the expressions written above make clear. Before turning to GWs, let us briefly review the explicit case of axion electrodynamics. There, the effective current is given by where the derivative will only act non-trivially on the axion due to the Gauss-Faraday law, ∂ νF µν = 0. Comparing with Eq. (S7), we see M µν eff = g aγγ aF µν , so that [44][45][46] P = g aγγ aB, M = g aγγ aE. (S12) Substituting these into Eq. (S10) recovers the equations of axion electrodynamics. If the axion field is non-relativistic, such that |∂ t a| |∇a|, then the leading contribution arises from the effective 3-current j eff ∂ t P = g aγγ ∂ t (aB), which is the effect that underpins much of the axion dark-matter program. For relativistic axions, however, all contributions are relevant and the equations match those considered in Ref. [54].

B. Gravitational Waves as a Source of the Effective Current
We will show now that, in a GW background and in absence of ordinary charges, Maxwell's equations can be written in the form of Eq. (S10). In curved space-time, the homogenous Maxwell equations read [60] whose solution is F αβ = ∂ α A β − ∂ β A α , independent of the background metric. This demonstrates that the homogeneous Maxwell's equations are not affected by the presence of the GW. On the other hand, their inhomogeneous counterparts become where j µ describes conventional charges and currents. Using the properties of the divergence, this can be written as The smallness of |h µν | implies that all expressions can be expanded perturbatively. To first order in h, we have which leads to In regions where j µ = 0, the source term on the right-hand side can be written as an effective current induced by the GW Conceptualizing the effect of the GW as an effective current was discussed in Ref. [22], working in the TT gauge for a uniform field. The extension to the proper detector frame and the importance thereof was pointed out in Ref. [2]. Note that, while there will be O(h) corrections to the electromagnetic fields -indeed, these are the fields we aim to search for experimentally -these do not enter Eq. (S18), as such terms would be O(h 2 ). Instead, the fields that enter are those that were established by the experimental apparatus. To keep track of this, one can introduce a bookkeeping notation F µν = F µν 0 + F µν h + O(h 2 ), in which case only F µν 0 enters the effective current. (A similar notation is often employed in the axion literature, where the expansion parameter is g aγγ , see e.g. Refs. [46,54].) We will occasionally make use of this notation, but where we do not, the relevant order of terms can always be determined from context. Furthermore, in the presence of ordinary currents, care must be taken. As follows from Eq. (S17), unless h = 0, one cannot simply add j µ and the effective current in order to account for GW effects. In fact, while ∂ µ j µ eff = 0, for ordinary currents we have instead, ∂ µ ((1 + h/2)j µ ) = 0, as a consequence of charge conservation in an arbitrary space-time, i.e 0 = √ −g ∇ µ j µ = ∂ µ ( √ −g j µ ).
Comparing Eqs. (S7) and (S18), we see that the expression in parentheses is a skew 2-index tensor, and can therefore be interpreted as an effective magnetization-polarization tensor. We can then immediately draw on the results of the previous section, and find as stated in Eq. (2). As suggested in the main body, an avenue towards the detection of this effective current is repurposing axion haloscopes, where one has a large magnetic field, but E 0 0. In this case, For a general magnetic field configuration, we will have both an effective charge and 3-current, and accordingly corrections to both sourced Maxwell's equations as in Eq. (S10).

C. The Biot-Savart Law and Induced Magnetic Fields
Having discussed general aspects of GW electrodynamics in the previous section, we now turn to the concrete calculations that underpin the detection scheme proposed in the main text. We focus on an external toroidal magnetic field -as appropriate for ABRA, SHAFT, and DMRadio-50L -and compute the induced magnetic field generated by the effective current.
Let us begin by justifying that the appropriate starting point for our calculation is the Biot-Savart law in (6). To do so, we show that the effective charge is not required to compute the induced magnetic field to all orders in ω, and that we can neglect the displacement current to O(ω 4 ). At O(h), we can write Maxwell's equations as These equations couple B h and E h . In order to isolate B h , we take a curl of the Ampère-Maxwell equation, and then substitute in Faraday's law and the magnetic Gauss' law. Doing so results in The general solution includes an additional term that depends on the boundary conditions of the instrument. We note, however, that determining the effect of the latter applies for any effective current, both axions and GWs, and this issue has yet to be settled for the conceptually simpler axion problem. If, at that time, the issue is settled using the language of j eff , then the result can simply be extended to GWs [62]. When 1/ω is much larger than the characteristic scale of our instrument -usually a good approximation when using low-mass axion haloscopes -we can simplify the result further. In particular, as the effective current for a plane GW satisfies j eff ∝ ω 2 e −iωt , Eq. (S23) implies that Accordingly, if we work only to O(ω 3 ), it is consistent to neglect displacement currents and start our calculation with the Biot-Savart law. This is assumed in the main text and will be adopted for remainder of the SM. That we can neglect the displacement current is equivalent to the magnetoquasistatic approximation conventionally adopted in low-mass axion haloscope calculations, see for example Ref. [3]. For a pickup loop parallel to the x − y plane, the only relevant component for computing the magnetic flux is For a toroidal magnetic field as given in Eq. (5), the azimuthal and radial components of j eff can be determined by direct calculation, where κ = ik · r. The components of h TT ij are determined by, for example, h TT ρρ = (ê ρ ) i (ê ρ ) j h TT ij , and take the following explicit forms All expressions within square brackets in Eq. (S26) are O(ω 0 ) at leading order, and therefore each line, except the last one, contributes at ω 2 . To leading order and taking φ h = 0, these equations reduce to those presented in the main text in Eq. (11). Similar expressions can also be derived for more complicated forms of the magnetic field, as, for instance, used in the SHAFT experiment. We will not state the corresponding relations here, although we note that these were accounted for in computing the sensitivity shown in Fig. 3. Before moving on to compute the magnetic flux, we note in passing that in our discussion in the main text and in the next subsection, we placed the pickup loop at the vertical center of the apparatus for simplicity. From the expressions above, we see that this simplifies the discussion as the terms odd in z vanish once we perform the volume integral over the effective current. However, more general placements of the pickup loop may be of interest, for example to include multiple pickup loops with different geometries, which would allow one to distinguish different incident directions and to differentiate between an axion and a GW signal. For this purpose, we observe that the terms proportional to z in Eq. (11) typically oscillate more strongly with φ than those proportional to ρ, such that the additional signal component which can be recovered away from the vertical center is (mildly) suppressed. In the case of a single pickup loop, this leads to a mild preference of placing the loop at the vertical center. We emphasize, however, that the expressions for the effective current and the induced magnetic field given above are fully general and can be used compute the flux through any pickup loop geometry.

D. The Magnetic Flux
Having calculated the magnetic field, we now turn to the magnetic flux through the pickup loop, The integration over the toroid is fixed by the geometry of the magnetic field, and we will restrict our attention to the detector geometry depicted in Fig. 2. What remains, however, is to fix a configuration for the pickup loop. As in the main text, we will consider two cases. In order to highlight the distinction, in Fig. S1 we plot the leading contribution to B z (r ) for a GW incident along the x-direction, taking r = R = a = H/4. The sinusoidal variation emphasized from the simplified result in Eq. (12) is clearly visible. From this, it follows that a circular pickup loop will be insensitive to the leading flux-the contribution from opposite sides of the pickup loop are cancelled. Nevertheless, as this is the loop geometry used by axion experiments we will consider it below. We will then consider the case of a partial loop, from which the optimal readout of a figure-8 geometry, where the contributions from each side can be added coherently, can be inferred. Circular loop: This is a highly symmetric configuration, allowing significant simplification. Let us first notice that φ in Eq. (S28) only enters in the combination r − r: it does not appear in the effective current. Explicitly, we have In particular, we observe that φ only appears in the combination φ − φ. We can thus largely remove φ by performing a shift φ → φ + φ, leaving The shift will also impact the range over which φ is integrated. If, however, we integrate around a complete circle (as required for a circular pickup loop), then this shift has no impact: the range remains φ ∈ [0, 2π). There are then two simplifications that occur. Firstly, the term proportional to j ρ vanishes as can be confirmed by performing the φ integral (note, for example, that the integrand is odd in φ ). This is true independent of the form of j ρ , as mentioned in the main text. Secondly, the only φ dependence left in the integrand sits in j φ , so that we can rewrite the flux as Using the form of j φ in Eq. (S26) we can perform the integral over φ to obtain FIG. S1. The leading ω 2 contribution to the z component of the GW induced magnetic field inside the torus, as a function of the loop coordinates (ρ , φ ). We have explicitly fixed the incident direction of the GW to bek =ê x , as shown. The sinusoidal variation of the leading magnetic field can clearly be seen, which if integrated over a circular pickup loop will vanish. Oppositely oriented semicircles -the figure-8 configuration -will instead be optimally sensitive to such a field. For the figure we show only the contribution from the × polarization, assume an ABRA-10 cm like configuration with r = R = a = H/4, show the magnetic field at the vertical center of the torus, and evaluate the e −iωt factor at t = 0. In the limit H R + a the quantity we plot can be determined from Eq. (12).
pickup loop, and that the leading contribution is O(ω 3 ).
In the main text we presented results derived assuming R, a H 1/ω. In this limit, e iωz cos θ h 1, the terms in Eq. (S31) odd in z vanish for a pickup loop placed at the vertical center of the apparatus, and the only remaining integral over z to perform is Together with Eq. (S32), this leads to Eq. (9) in the main text.
Partial loop: As shown by the calculation leading to (S32), the decoupling of h + together with the absence of an O(ω 2 ) contribution is a consequence of the reflection symmetry of the circular pickup loop. Motivated by this, we consider a pickup loop forming a circular arc subtending the angles 0 < φ < φ max . As explained above, the integration along the z-axis may be formally performed by replacing dz/|r − r| 3 → 2/|r − r| 2 z=0 and setting z = 0 everywhere else. Under the this approximation, the flux in Eq. (S28) takes the form The flux for specific geometries can be immediately obtained from this result. This includes the figure-8 configuration discussed in the main text; in particular, Eq. (13) follows from this result.

III. GRAVITATIONAL WAVE SOURCES IN THE ULTRA-HIGH FREQUENCY BAND
Here we summarize the possible sources that the search proposed in the main text may detect, with a particular focus on the case of PBH binaries. The discussion is not intended to be exhaustive, and we refer to Ref. [1]. There are no known sizable astrophysical sources in this frequency range, implying that any GW detection at UHF is a smoking gun signal of new physics. Many models that posit a completion of the Standard Model of particle physics at high energies also predict additional dynamics in the early Universe (such as phase transitions, formation of topological defects, or non-perturbative (p)reheating dynamics after inflation) which source gravitational radiation that can, depending on model parameters, constitute a sizable fraction of the total radiation energy present in the early Universe. Taking into account the cosmological red-shift, these stochastic GWs are observable at frequencies f ∼ 100 MHz/ * (T * /10 15 GeV) where T * is the temperature of the Universe when the GWs are sourced and * < 1 indicates the GW wavelength in units of the Hubble horizon at the time of production. This result suggests that the UHF band could be ideal for searching for new dynamics present at energies well beyond what we can hope to probe directly. However, such cosmological GWs contribute to the radiation energy budget of the Universe and thus impact big bang nucleosynthesis and the decoupling of the cosmic microwave background. They are therefore constrained by bounds on the effective number of additional neutrinos, ∆N eff [63,64], in particular ρ GW /ρ c 10 −5 ∆N eff , where ρ c is it critical density. For a broadband spectrum, this constrains the characteristic strain to be at most h c,sto 10 −29 (100 MHz/f ) ∆N eff . This is several orders of magnitude below the most optimistic sensitivities we presented in Fig. 1, suggesting that even with the advantageous volume scaling of the detection strategy suggested in this work, such a signal is out of reach in the foreseeable future.
Such strong constraints do not apply to isolated GW sources in the late Universe. To highlight this, we discuss in more detail the possible signal from PBH binaries, before briefly commenting on other possible sources in the late Universe at the end of this appendix. For simplicity, let us consider a circular binary of two PBHs which have equal mass m PBH . Then, the amplitude of the GW signal emitted from a binary at a distance D along the symmetry axis of the circular orbit is [59] We emphasize that from the above it can be seen that frequencies in the MHz band and beyond correspond to light black holes, with m PBH M , which excludes a stellar-origin for the black hole sources in this frequency band. The probability of observing such a PBH binary is determined by the formation rate of the binaries and their local density. The light PBHs we are interested in are dominantly formed in the early Universe. Through the distribution of all PBHs, a subset will be close enough that their gravitational attraction leads the pair to decouple from the Hubble flow and form a gravitational bound object. (In this scenario, a third PBH is in fact required in the process to ensure the conservation of angular momentum.) Assuming a Gaussian distribution for the primordial density perturbations, normalized to give a PBH abundance which constitutes a fraction f PBH of the total dark-matter density, PBHs are formed in rare events following Poisson statistics and the binary formation rate is obtained as [65,66] where for simplicity we have taken the PBH mass distribution to be concentrated at m PBH . This formation rate receives corrections from different effects, see Sec. 4.6 of Ref. [47] for a detailed summary and relevant references. We include two effects which suppress the merger rate, denoted by S early and S late in Eq. (S37). The first of these is a FIG. S2. Strain sensitivity h th required to observe a signal from PBH binary systems (either from the inspiral or merger phase) at rates between ten events per year and one per ten-thousand years. The case for one per year was reproduced in Fig. 1. Details of how these sensitivities were computed are provided in the text, however, we caution that there remain uncertainties in this computation (for instance from the impact of accretion on the merger rate) that we have neglected.
backreaction from the surrounding matter in the early Universe, which suppresses the merger rate for f PBH < 0.01 [66].
Since microlensing bounds largely constrain f PBH to lie below this value [48], this effectively changes the exponent of f PBH in Eq. (S37) from 53/37 to ∼2. (For an alternative treatment of this effect, see Refs. [58,66], although we confirmed adopting this alternative form for S early does not significantly alter our results.) Interactions of the binary with matter in the late Universe can also disrupt the system, again suppressing the merger rate, as encoded in S late [66]. It is also expected that accretion will impact the merger rate, however the degree to which it does so is less certain, and so for the simple estimate we provide here, we neglect the effect of accretion. We note, however, that a recent study in Ref. [67] suggested that for the PBH masses we consider, significant accretion is unlikely, so that neglecting its effect may well be a reasonable assumption. From Eq. (S37), we obtain the number of expected merger events per year in a volume of 1 kpc 3 , for a given m PBH and f PBH . If we have a GW telescope with a strain sensitivity h th at a frequency f , then using Eq. (S35), we can determine what distance to which we can detect such a binary signal, for a given m PBH . Combining the two results, we can determine the expected rate at which the telescope will observe a signal from PBH binary events, There are two additional features we have added to this expression beyond the above discussion. Firstly, the quality factor or coherence of the signal, Q, is proportional to the time remaining until the merger given an emitted GW frequency f and PBH mass m PBH , yielding Q ∝ f 2 /ḟ ∝ m −5/3 PBH [59]. For a binary system very close to merger, f /f 2 ∼ 1 and hence Q ∼ 1. Consequently, we can normalize Q to the mass which would be at the point of merger for the frequency considered, using Eq. (S36). As discussed in the main text, coherent signals are more readily detected, and so we have included this factor as an enhancement of Q 1/4 . Note that the scaling Q 1/4 is specific to the detector concept proposed here. When comparing to generic GW detectors, one often instead uses the characteristic strain, h c = fh(f ) withh(f ) denoting the Fourier transform of the GW signal. For a PBH binary sufficiently far from merger, h c Q 1/2 h PBH +,× . We caution that this caveat also applies when comparing to some of the detector sensitivities shown in grey in Fig. 1. (More generally, the distinction between Q 1/2 and Q 1/4 scaling for axion haloscopes depends on whether the observation time is less than or greater than the signal coherence time, and we have assumed the latter throughout.) Secondly, R 0 is the rate averaged over the cosmological volume of the binaries, whereas locally we live in a significant matter overdensity, the Milky Way. This reality is included in the factor of δ(r). We assume that the PBH binaries simply track the overall DM abundance, resulting in a simple linear dependence on δ(r) in the integrand. It enters in particular with a different power than a primordial overdensity or f PBH enter, since only the latter impact the binary formation processes in the early Universe. As pointed out in Ref. [68], for light PBHs the local DM overdensity in the Milky Way halo can boost the observed merger rate significantly. We parametrize the Milky Way halo by a NFW mass profile [69,70] up to the virial radius r 200 207.19 kpc, with a local density ρ 0.31GeV/cm 3 at the location of the solar system, r 8.13 kpc [71], where the halo parameters are determined from Ref. [72]. Fig. S2 shows the detector sensitivity required to detect a fixed expected number of merger event per year, one curve of which was included in Fig. 1. A relatively large event rate, Γ = 10/year requires a very good detector sensitivity, whereas rare event may be detected with a more moderate detector sensitivity. Here we have marginalized over the PBH mass in Eq. (S38), i.e. for any given frequency we chosen m PBH such as to maximize the rate Γ. The fraction of PBH dark matter f PBH is taken to saturate the microlensing bound at that PBH mass [48]. Due to the mild scaling with Q, this rate is dominated by merger events rather than the early inspiral phase in most of the parameter space. For fixed values of m PBH and f PBH these results are consistent with Ref. [58] after recasting the dimensionless GW amplitude h in terms of the characteristic strain h c . (Note that such high-frequency GW bursts can also be constrained by the non-observation of corresponding GW memory signals in ground-based interferometers [73]; see also Ref. [74].) In regions where the density profile of the Milky Way halo becomes irrelevant because the effective detector volume set by the Heaviside function in Eq. (S38) is much larger (at small frequencies) than the Milky Way halo, the strain sensitivity scales as h ∝ m 79/111 PBH ∝ f −79/111 , where the last relation only holds for signals close to merger, see Eq. (S36), and we have neglected Q. The difference in overall normalization between these two regions is given by the local DM overdensity at the position of the solar system. When the effective detector volume becomes comparable to the size of the halo, the DM density profile leads to the feature visible at intermediate frequencies in Fig. S2.
We end with several brief comments on additional late time sources. PBH can furthermore source GWs when they scatter off each other without merging. The resulting spectrum is not monochromatic and peaks at a frequency determined by the specific hyperbolic orbit that the PBHs follow. Observable signals require significant clustering or very rare events such as close to parabolic encounters [75][76][77]. Another possible late time source is axion superradiance [78]. Clouds of axion-like particles could form around black holes when the axion Compton wavelength matches the Schwarzschild radius. This would lead to an emission of GWs with a wavelength set by this same scale, assuming that the axion cloud constitutes 0.1% of the black hole mass. The signal is expected to be highly monochromatic. These examples illustrate that different search strategies will need to be implemented to optimally search for different possible GW sources. Generally, however, the GW signal is expected to be less coherent than axion dark-matter, and so in general strategies searching for broader signals in the frequency domain will need to be devised. (In this sense there are strong similarities with the search for relativistic axions, for further discussion see Ref. [54].)