Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies

The first direct measurement of gravitational waves by the LIGO and Virgo collaborations has opened up new avenues to explore our Universe. This white paper outlines the challenges and gains expected in gravitational-wave searches at frequencies above the LIGO/Virgo band, with a particular focus on Ultra High-Frequency Gravitational Waves (UHF-GWs), covering the MHz to GHz range. The absence of known astrophysical sources in this frequency range provides a unique opportunity to discover physics beyond the Standard Model operating both in the early and late Universe, and we highlight some of the most promising gravitational sources. We review several detector concepts that have been proposed to take up this challenge, and compare their expected sensitivity with the signal strength predicted in various models. This report is the summary of the workshop “Challenges and opportunities of high-frequency gravitational wave detection” held at ICTP Trieste, Italy in October 2019, that set up the stage for the recently launched Ultra-High-Frequency Gravitational Wave (UHF-GW) initiative.


Introduction
Gravity and electromagnetism are the only two long range interactions in nature, but over the centuries we have explored the Universe only through electromagnetic waves, covering more than 20 orders of magnitude in frequencies, from radio to gamma rays. The discovery of gravitational waves in 2015 has opened a totally new window to observe our Universe (Abbott et al. 2016b).
Judging by what happens with electromagnetic waves, there should be interesting physics to be discovered at every scale of gravitational wave frequencies. Current and planned projects such as pulsar timing arrays, as well as ground-and spacebased interferometers will explore gravitational waves in the well-motivated range of frequencies between the nHz and kHz range. However, both from the experimental and the theoretical point of view it is worth to consider the possibility to search for gravitational waves of much higher frequencies, covering regimes such as the MHz and GHz, see for instance Cruise (2012).
A strong motivation to explore higher frequencies from the theoretical perspective is that there are no known astrophysical objects which are small and dense enough to emit at frequencies beyond 10 kHz. Any discovery of gravitational waves at higher frequencies would thus indicate new physics beyond the Standard Model of particle physics, linked e.g., to exotic astrophysical objects (such as primordial black holes or boson stars) or to cosmological events in the early Universe such as phase transitions, preheating after inflation, oscillons, cosmic strings, thermal fluctuations after reheating, etc., see Caprini and Figueroa (2018) for a recent review.
For early Universe cosmology, gravitational waves may be the only way to observe various events. In particular for the time between the Big Bang and the emission of the cosmic microwave background radiation, electromagnetic waves cannot propagate freely, whereas, due to the weakness of gravity, gravitational waves decouple essentially immediately after being produced and travel undisturbed throughout the Universe forming a stochastic background that could eventually be detected. Even though it may not be easy to unambiguously determine the concrete cosmological source of a gravitational-wave signal, its cosmological nature of the spectrum may be identified, similar to what happened with the original discovery of the cosmic microwave background.
In this context, the existence of a stochastic spectrum in the range from kHz to GHz is well-motivated: causality restricts the gravitational wave wavelength to be smaller than the cosmological horizon size at the time of gravitational wave production. This roughly implies a gravitational-wave frequency above the frequency range of the existing laser interferometers Virgo (Acernese et al. 2015(Acernese et al. , 2019, LIGO (Aasi et al. 2015;Abbott et al. 2016a;Buikema et al. 2020;Tse et al. 2019) and KAGRA (Akutsu et al. 2019;Aso et al. 2013) for any gravitationalwave production mechanism that happens at temperatures larger than 10 10 GeV, 1 assuming radiation domination all the way to matter-radiation equality. In particular, GHz frequencies correspond to the horizon size at the highest energies conceivable in particle physics (such as the Grand Unification or string scale) and phenomena like phase transitions and preheating after inflation would naturally produce gravitational waves with frequencies around the GHz range.
Established gravitational-wave detector designs are limited to frequencies up to the kHz range. In particular, resonant mass detectors, going back to the original bar design of Weber (1967), focused on isolated high frequencies, often targeting known millisecond-pulsar frequencies. Similarly, the well-established interferometric gravitational-wave detectors LIGO, Virgo and KAGRA cover parts of the highfrequency band up to a few kHz. For the purposes of this white paper, we shall therefore use the expression high-frequency gravitational waves to refer to frequencies that are above the LIGO detection band, i.e., starting from around 10 kHz. In particular, taking inspiration from the electromagnetic spectrum, we denote the MHz to GHz range by Ultra High-Frequency Gravitational Waves (UHF-GWs). Several proposals have been made for pushing the high-frequency end of interferometric detectors into this region, however, detectors for the MHz, GHz and THz frequency bands require radically different experimental approaches.
Over the years, there have been isolated attempts to search for gravitational waves of very high frequencies and a few proposals have been put forward. These new concepts have largely been suggested in the form of theoretical papers with no serious discussion of the potential experimental noise sources that might limit their performance, or occasionally, bench tests of early prototypes. The current status of many of these ideas must be regarded as highly preliminary. The published concepts span a wide range of technologies with no real consensus yet as to where to concentrate the community effort. In addition to the selection of suitable technological pathways towards a serious attempt at a detection at high frequencies, there needs to be an identification of the most realistic sources and thereby the waveforms and spectra for which such detectors should be optimised. This process demands a close collaboration of theorists and experimentalists.
The goal of this report is to summarise and start a dialogue among the specialised community regarding the importance and feasibility to explore searches for highfrequency gravitational waves. We are aware that this may be a long term goal but are convinced that the physics motivation is strong enough to start a systematic study of the different sources of high-frequency gravitational waves and their potential detectability. It is the purpose of this white paper to put together the different ideas both from theory and experiment to explore the importance of searching for high-frequency gravitational waves. The origin of this initiative was a workshop organised at ICTP in October 2019 'Challenges and Opportunities of High-Frequency Gravitational Wave Detection' where members of the theoretical and experimental communities interested on high-frequency gravitational waves got together to explore the motivations and challenges towards this search. This workshop and the present white paper set the stage for the launch of the Ultra-High-Frequency Gravitational Wave (UHF-GW) initiative, 2 whose goals include supporting the testing phase of currently existing detector proposals and stimulating the technological developments necessary to come up with new schemes for gravitational-wave detectors at high frequencies.
The remainder of this report is organized as follows: Sect. 2 introduces some basic concepts and notation to discuss different types of gravitational-wave sources and to relate them to experimental sensitivities. An overview over gravitationalwave sources in the late and early Universe is given in Sect. 3, followed by a discussion of different detector concepts in Sect. 4. We conclude in Sect. 5. For a summary of the various detector concepts and the corresponding sensitivities see Sect. 4.3 and Table 1. For a summary of the various sources see Sect. 3.1, Figs. 1, 2, Appendix 1 and Tables 2, 3.
We collect here a few acronyms that will be used throughout the paper: Gravitational Wave (GW), Ultra High-Frequency Gravitational Waves (UHF-GWs), Cosmic Microwave Background (CMB), Black Hole (BH), Innermost Stable Circular Orbit (ISCO), Big Bang Nucleosynthesis (BBN).

Setting up the notation: comparing different GW sources and detectors
Depending on the source/detector, the strength of GWs, detector noise, and signalto-noise ratio are described using various different metrics (Maggiore 2007). In general, before using any given metric, it is important to make sure that it is appropriately defined for the scenario under consideration. In this section we summarize the relevant quantities and notation. We follow the definition in Allen and Romano (1999) for stochastic strain sources, and definitions in Moore et al. (2015) for time-dependent strain sources.
2.1 Gravitational-wave sources at high frequencies 1. For stochastic GWs, for example those coming from cosmological sources, a spectral density prescription is most suitable. The most common models assume that they are approximately isotropic, unpolarized, stationary, and have a Gaussian distribution with zero mean. They can thus be fully defined by the second moment (Allen and Romano 1999): Hereh A ðf ; XÞ is the Fourier transform 3 of the time-dependent strain in the GW polarization A, solid angle X, evaluated at a frequency f. S h ðf Þ denotes the onesided power spectral density. The energy-density q GW in GWs per logarithmic frequency interval is represented by X GW , conventionally normalized by the critical energy density q c ¼ 3H 2 0 =ð8pGÞ with G denoting Newton's constant and H 0 denoting the Hubble parameter today. We will denote the current value of X GW by X GW;0 . The power-spectral density can be directly related to the 00-component of the stress energy tensor, in turn yielding: Often, a dimensionless characteristic strain is assigned to the normalized energy density for stochastic GWs (see e.g., Thrane and Romano 2013;Romano and Cornish 2017) 2. For inspiral sources, such as BH mergers, a time-dependent strain h(t) can be obtained directly from Einstein's equations. Inspirals have an evolving frequency evolution, so usually the stationary phase approximation is used to obtain an analytical form for the Fourier transformhðf Þ (Maggiore 2007). The characteristic strain for such sources with inspiralling frequency can be defined so as to take the frequency evolution into account in the GW strength (Moore et al. 2015), Assuming that h 0 is the amplitude of the GW from the inspiral, i.e., the amplitude of the periodic function h(t), this results in the characteristic strain:

Detectors
Each detector has a different way of searching for GWs, with different antenna patterns, frequency bands, binning, etc. This should be taken into account when defining the appropriate noise and signal-to-noise metrics. For interferometers (such as LIGO) the impact of spatial antenna patterns is of the order of unity. For simplicity, the detector noise floor is usually specified assuming the noise is stationary and Gaussian (even though in reality it is usually neither). Similar to the discussion of stochastic GWs, this noise floor is specified by using a power spectral density, The angular brackets denote an average over multiple realizations of the system, which is obtained repeating the measure of the noise over several well separated time intervals of the same length, see Maggiore (2007). 4 In order to measure this noise, a fast Fourier transform of the detector noise in the absence of the signal is performed. This measured noise is compared to a numerical model, comprising of the sum of all the noises in the detector. An analysis showing each individual noise source (measured or modeled) summing up to the total measured noise is called the noise budget. Unless otherwise specified, if a detector noise is specified in terms of spectral density, it should be treated as S n ðf Þ (with . For a visual comparison of signal strengths of inspirals and stochastic signals against detector sensitivities, conventionally a dimensionless noise amplitude as been introduced, denoted by h c;n (Moore et al. 2015) h c;n ffiffiffiffiffiffiffiffiffiffiffi ffi fS n ðf Þ p : Some experiments looking for long-lived sources (e.g. monochromatic sources or stochastic sources) can choose to average the noise over a long time. This gives an additional boost in signal-to-noise ratio, which is sometimes reported as an enhanced sensitivity S n;int ¼ S n N avg : If the detector is operating at a frequency f center and it is integrating for a time T obs then N avg ¼ T obs =T FFT where T FFT is the time duration of each spectrum (assuming the segments can be coherently averaged over the duration T obs ).

Signal-to-noise ratio
Understanding whether a signal is detectable using a particular detector requires development of a metric for the signal-to-noise ratio q.
• The most efficient signal-to-noise ratio metric for broadband detection of transient sources uses matched-filtering (Maggiore 2007;Allen et al. 2012;Moore et al. 2015): If the frequency ranges of both the signal and the detector are sufficiently broad, d ln f ¼ Oð1Þ, then h c;insp $ h c;n roughly corresponds to q ¼ Oð1Þ. This explains why h c;insp and h c;n from Eq. (8) are useful in assessing the reach of a particular broadband instrument looking for an inspiralling source. • For a resonant detector with no sensitivity outside a small bandwidth, this signalto-noise ratio simply collapses to the single frequency band of detection, indicating that a correspondingly larger threshold value ofhðf Þ is required to yield a detectable signal at fixed h c;n . • For detecting approximately monochromatic sources, the signal-to-noise ratio similarly collapses to a single frequency. In this case, Df in Eq. (11) is given by the frequency resolution, i.e., either the width of the signal or the detector resolution, whatever is the relevant limiting factor. For searches of monochromatic GWs that last over long times, various astrophysical effects like the Earth's motion need to be taken into account. • Detecting stochastic sources usually requires utilizing cross-correlation between two or more GW experiments to distinguish the GW background from the experiment's noise (see also Sect. 4.4). Therefore, defining a meaningful signalto-noise ratio for detection of stochastic sources requires careful consideration of the noise, location, and alignment of each individual experiment. The signal-tonoise ratio can be increased by using more independent experiments, observing for longer times, and optimizing the size of frequency bins. Usually, the strength of the signal will be much less than the detector noise, so cross-correlation can provide a signal-to-noise ratio greater than 1. For a simple case of M [ 1 colocated detectors, measuring in a frequency band from f min to f max for a time T obs , the signal-to-noise ratio can be written as (Thrane and Romano 2013): We redirect the reader to Allen and Romano (1999); Thrane and Romano (2013); Romano and Cornish (2017) for a full analysis.

Comparison of signal strength and noise for narrowband detectors
Since most high-frequency detectors are narrowband, 5 here we provide some handy expressions to compare signal strength and detector sensitivity for narrowband detectors. The most natural way to express a detector's sensitivity is in power or amplitude spectral density, S n ðf Þ or ffiffiffiffiffiffiffiffiffiffi ffi S n ðf Þ p . On the other hand for signal strengths, the most natural units can depend on the type of source -dimensionless characteristic strain for inspirals, amplitude or power spectral density for stochastic sources, and wave amplitude for long-lived monochromatic sources. In order to compare the signal strength and detector sensitivity, we often strive to convey them in the same units. Here we will provide two ways to achieve this for narrowband detectors. The underlying principle for both methods is to first write down a reasonable signal-to-noise ratio metric, and use that to derive the appropriate comparable quantity. The signal-to-noise ratio for stochastic and long-lived monochromatic sources will be enhanced due to the integration over the observation time, in contrast with signal-to-noise ratio for transient sources, which will depend on just a single observation.
If we are interested in assessing the utility of a given detector to search for GWs from various types of sources, it would be natural to include the integration-time information in the signal depending on the source. On the other hand, if we wish to compare various detectors' suitability to a given source, it is convenient to include the time and bandwidth information to convert the detector sensitivity to the source units. Here we provide ways to do both.
1. Inspirals: For inspiralling sources passing through the band of a narrowband detector, the signal-to-noise ratio is shown in Eq. (11), and the Fourier transformhðf Þ is related to the characteristic strain using Eq. (5). We desire S h;res;insp or h c;n;insp such that where Df represents the frequency range in which GWs are measured by a given detector, i.e., f AE Df =2. 2. Stochastic sources: Following the same prescription as above, we write the narrowband version of the signal-to-noise ratio in Eq. (12): Using Eqs. (15) and (4a), we can obtain 3. Monochromatic sources: For long-lived monochromatic sources, the natural unit to specify the signal strength is just the amplitude of the sinusoidal wave, h 0 . The signal-to-noise ratio buildup is from the long integration times and can be written as Using this, noise-equivalent signal ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S n;res;mono p or signal-equivalent noise h 0;n;mono can be written: 6 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S n;res;mono 6 As is the case in LIGO and Virgo, a more complicated analysis is needed if the coherence duration of the detector is shorter than T obs . In that case, fast Fourier transforms taken over the coherence duration are incoherently averaged. This leads to a modified scaling q 2 res;mono $ jh0j 2 N 1=2 avg TFFT Sn (Astone et al. 2014).
Note that for monochromatic sources, there is no need to invent a characteristic strain, the most 'characteristic' strain is the amplitude h 0 itself.

Sources
This section reviews various production mechanisms for GW signals in the highfrequency regime, typically in the range (kHz-GHz), that fall into two broad classes. In Sect. 3.2 we discuss sources in our cosmological neighbourhood, which emit coherent transient and/or monochromatic GW signals. In Sect. 3.3 we turn to sources at cosmological distances which typically lead to a stochastic background of GWs. We emphasize that all proposed sources, with the notable exceptions of the neutron star mergers discussed in Sect. 3.2.1 (kHz range) and the cosmic gravitational microwave background discussed in Sect. 3.3.3, require new physics beyond the Standard Model of particle physics to produce an observable GW signal. Thus, while being admittedly somewhat speculative, these proposals provide unique opportunities to shed light on the fundamental laws of nature, even by 'only' setting an upper bound on the existence of GWs in the corresponding frequency range.

Overview
Figures 1 and 2 summarize a representative selection of the sources which are discussed in more detail in the following subsections. The regions bounded by the colored curves illustrate the region of parameter space which may be covered by the corresponding source for appropriate parameter choices as specified below. Except for the cases of inflation with broken spatial reparametrization symmetry and the cosmic gravitational microwave background they should not be mistaken for GW spectra obtained for a fixed model parameter choice.
In the same figures, we also indicate the demonstrated (filled boxes) or expected (empty boxes) sensitivity of the detector concepts discussed in Sect. 4.1. In some cases we report two sensitivities for a single detector, using two different intensities of the same color (see for instance the case of levitated sensors), if the sensitivity depends on the details of the future implementation of the detector or on some assumptions needed to place the constraint. In the case of the levitated sensors the two colors refer to two different versions of the same detector concept: a 1 meter and a 100 meter implementation, see Sect. 4.2.1, the latter giving a better sensitivity and in the case of the radiotelescopes EDGES and ARCADE, the two sensitivities refer to the weakest and strongest possible cosmic magnetic field, whose value is needed in order to place the constraint, see Sect. 4.2.2.
The comparison of signal strength and detector sensitivity in these figures should be taken with great caution, and serves as a very rough illustration only. In particular, the signals in Fig. 1 are coherent (and partially transient) signals whereas the signals depicted in Fig. 2 are stationary isotropic stochastic signals. A given detector concept will be more or less suitable for these different types of signals, which is not accounted for in this illustration. Further restriction may apply. For example, the quoted sensitivity for the radio telescopes ARCADE and EDGES assumes a cosmological distance between source and observer. Regarding the possible signals, we have aimed to make realistic estimates of the largest possible    Fig. 1 for the reference to the various detector concept sections signals in different models. This does however not factor in the likelihood of such a signal occurring in the detector lifespan. This is in particular true for the coherent sources, see the respective subsections for details. Figure 1 shows representative examples of coherent sources. For simplicity, we take the factor ffiffiffi ffi 2f _ f q converting between the amplitude and characteristic strain of a GW to be unity, which is a good approximation at the merging frequency of compact objects. Moreover, we use a reference value of 10 kpc for the distance to all sources.
The amplitude of the oscillating GW signal is then given by Eqs. (21) and (30), respectively. • For signals from axion superradiance we consider both the axion annihilation and axion decay channel (see Sect. 3.2.4). The frequency of the signal is determined by the axion mass, which is in turn linked to the BH mass by the superradiance condition in Eq. (31). Inserting this into Eqs. (33) and (35) and taking a=l ¼ 1=2, ¼ 10 À3 and M BH [ M yields the curves depicted. Figure 2 shows models producing stochastic GW signal: interestingly, most of them are concentrated in the UHF band. These are produced in the early Universe and are thus subject to the cosmological constraint on the number of effective degrees of freedom during BBN and at CMB decoupling, see Sect. 3.3.
• In certain models, inflation (Sect. 3.3.1) can yield a signal stretching over a broad frequency range (see Eq. (38)), with an amplitude determined by Eq. (39) and Eq. (41), respectively. Here in the case of inflation with extra-species we have taken the parameter n (defined in Eq. (39)) to be bounded by the perturbative limit, and in the case of inflation described by an effective field theory with broken spatial reparametrization symmetry we have chosen the speed of sound and the spectral tilt to be c T ¼ 1 and n T ¼ 0:2, respectively. Moreover, inflation models with strongly enhanced scalar fluctuations (P f .10 À2:5 Þ can source GWs with X GW; 0 .10 À9 at second order in cosmological perturbation theory. • For preheating (Sect. 3.3.2), we show typical values for models with parametric resonance in quadratic (right green box) and quartic (left green box) potentials as well as oscillons. In the latter case the frequency is set by the mass of the scalar field through Eq. (45), where here we have chosen the mass of the scalar field to be 10 10 GeV\m\10 13 GeV with X ¼ 100, while the amplitude is the typical value inferred from numerical simulations. • For the cosmic gravitational microwave background, we show the spectrum given by Eq. (46) with T max ¼ 10 16 GeV, which is the upper bound on the reheating temperature set by the constraints on the tensor-to-scalar-ratio (Akrami et al. 2020). • For phase transitions (Sect. 3.3.4), we assume a fixed latent heat, number of relativistic degrees of freedom and wall velocity. We also assume that sound waves do not last a Hubble time, such that the amplitude scales as the square of the inverse time scale of the transition. The peak frequency and amplitude are then given by Eqs. (47) and (48), where we consider transition temperatures T Ã \10 16 GeV. • As an example for topological defects (Sect. 3.3.5) cosmic strings lead to a broad spectrum with an amplitude given Eq. (49), where the string tension for stable cosmic strings is bounded by Gl\10 À11 whereas for metastable cosmic strings it can by as large as Gl ' 10 À4 above the LIGO frequency range. The spectrum of gauge textures is described by Eq. (51), where here we have chosen the symmetry breaking scale to be 10 12 GeV\v\10 19 GeV.

Late Universe
In this section we revise various sources that are relevant for high-frequency GW production and are active in the late Universe. A summary of these sources is given in Fig. 1 and Table 2 in Appendix 1.

Neutron star mergers
For not too high binary masses the merger of two neutron stars avoids the prompt collapse to a BH and leads to the formation of a massive rapidly rotating and oscillating neutron star remnant. The oscillations of this remnant are very characteristic of the incompletely known equation of state of high-density matter and generate GW emission in the kHz range (see Fig. 3). For instance, the dominant oscillation frequency of the post-merger phase (f peak in Fig. 3) scales tightly with the radii of non-rotating neutron stars (Bauswein et al. 2016). These radii are uniquely determined by the equation of state of neutron stars, and are therefore particularly valuable messengers of the underlying high-density matter physics (see e.g., Oertel et al. 2017 for a review). Simulation results show a tight correlation between the dominant GW frequency and neutron star radii; a fit to the data for fixed binary masses describes the relation with a maximum residual of only a few hundred meters, allowing for accurate radius measurements. Subdominant features in the GW spectrum (see Fig. 3) contain additional information about the equation of state and may also reveal the dynamics of the remnant, which is indispensable for a complete multi-messenger interpretation of neutron star mergers. The presence or absence of post-merger GW emission from a neutron star remnant on its own informs about the outcome of the merger (neutron star or BH). In combination with the measured binary masses, this information allows to constrain the threshold binary mass for prompt BH collapse, which is somewhere in the range (2.9-3.8) M , depending on the equation of state. This threshold depends sensitively on the maximum mass M max of non-rotating neutron stars. Obtaining the threshold mass for prompt BH formation through post-merger GW emission will yield a robust determination of the unknown maximum mass M max of non-rotating neutron stars, which is another important equation of state property that probes the very high-density regime (Bauswein et al. 2013). A robust measurement of M max is also relevant for stellar astrophysics since it, for instance, affects the outcome of core-collapse supernovae. Pulsar observations only yield accurate lower bounds on M max .
Generally, equation of state inference from the post-merger stage is complementary to other constraints, e.g., from the inspiral phase. The complementarity concerns the probed density regime, which is generally higher in the post-merger phase, and methodological aspects. Hence, the detection of post-merger GW emission is of highest importance to understand properties of high-density matter including the opportunity to probe the presence of a phase transition to deconfined quark matter (Most et al. 2019;Bauswein et al. 2019).
The different features of the post-merger GW emission have frequencies in the range (1-5) kHz, with the dominant peak between 2 and 4 kHz (see Fig. 3).

Fig. 3
Typical GW spectrum of the cross polarization of a 1.35-1.35 M merger along the polar direction at a distance of 20 Mpc. h eff ;Â ¼h Â ðf Þ Á f with the Fourier transform of the waveform h Â and frequency f. f peak , f spiral and f 2À0 are particular features of the post-merger phase, which can be associated with certain dynamical effects in the remnant. Since the simulation started only a few orbits before merging, i.e. at a relatively high orbital frequency, the power at lower frequencies (below $ 1 kHz) is massively under-represented in the shown spectrum, and the low-frequency part of the spectrum does not show the theoretically expected power-law decay. The thin solid lines display the spectra of the GW signal of the post-merger phase only revealing that the peaks are indeed generated in the post-merger phase. Dashed lines show the expected unity signal-to-noise sensitivity curves of Advanced LIGO (red) and of the Einstein Telescope (black). Image reproduced with permission from Bauswein et al. (2016), copyright by SIF/Springer Simulated injections show that at a distance of 40 Mpc (comparable to that of GW170817) a strain sensitivity of roughly ffiffiffiffi ffi S n p ' 3 Â 10 À24 Hz À1=2 is required for a detection of the main features (Torres-Rivas et al. 2019). Hence, measurements can be anticipated with a small sensitivity improvement either of Advanced LIGO/ Virgo/KAGRA or with a dedicated high-frequency instrument like NEMO (Ackley et al. 2020) (see Sect. 4.1.1).

Mergers of light primordial black holes
The low effective spins and progenitor masses of the BH mergers detected by LIGO/Virgo have revived the interest for primordial BHs in the range (1-100) M (Bird et al. 2016;Clesse and García-Bellido 2017;Sasaki et al. 2016), which could constitute (part of) the observed dark matter abundance. In this context, detecting a sub-solar mass BH would almost clearly point to a primordial origin. 7 The frequency associated to the ISCO when the inspiral GW emission is close to maximal 8 , is given by with m 1 and m 2 the masses of the two binary components and M denoting the solar mass. A good estimation of the GW strain produced at a given frequency f is provided by the Post-Newtonian approximation (Antelis et al. 2018) where M ðm 1 Â m 2 Þ 3=5 =ðm 1 þ m 2 Þ 1=5 is the binary chirp mass, D is the distance to the observer, G is Newton's constant and c is the speed of light. For an equalmass binary and an experiment of strain sensitivity h det , the corresponding astrophysical reach D max is given by High-frequency GW detectors could therefore detect or set new limits on the abundance of light, sub-solar mass primordial BHs, in particular if they have an extended mass distribution. Frequencies in the range ð10 4 À 10 12 Þ Hz correspond to a primordial BH mass range ð10 À9 À 10 À1 Þ M . In particular, the planetary-mass range, in which recent detections of star and quasar microlensing events (Niikura et al. 2019;Hawkins 2020;Bhatiani et al. 2019) suggest a dark matter fraction made of primordial BHs of f PBH $ 0:01, could be probed in a novel way.
There are two possible formation channels of primordial BH binaries, introduced hereafter: 1. Primordial binaries: they come from two primordial BHs that were formed sufficiently close to each other for their dynamics to decouple from Universe expansion before the time of matter-radiation equality (Nakamura et al. 1997;Sasaki et al. 2016). The gravitational influence of one or several primordial BHs nearby prevent the two BHs to merge directly, leading to the formation of a binary. In some cases, the binary is sufficiently stable and it takes a time of the order of the age of the Universe for the two BHs to merge. If the primordial BHs have a mass spectrum qðmÞ and are randomly distributed spatially, and that early forming primordial BH clusters do not impact the lifetime of those primordial binaries (a criterion satisfied for f PBH . where f PBH is the integrated dark matter fraction made of primordial BHs, m 1 and m 2 are the masses of the two binary BHs components and qðmÞ is the density distribution of primordial BHs normalized to one ( R qðmÞd ln m ¼ 1). 9 Assuming f PBH Â qðmÞ ' 0:01 at planetary masses in order to pass the microlensing limits, and considering only almost equal-mass mergers (m 1 $ m 2 $ m PBH ) that produce the highest strain, one obtains a merging rate of In turn, using Eq. (20), one obtains the required GW strain sensitivity to detect one of these merger events per year, which can be typically targeted by GW experiments operating at frequencies from kHz up to GHz, see e.g., Wang et al. (2019). 2. Capture in primordial BH haloes: the second binary formation channel is through dynamical capture in dense primordial BH halos. As any other dark matter candidate, primordial BHs are expected to form halos during the cosmic history. The clustering properties typically determine the overall merging rate.
For instance, for a monochromatic mass spectrum and a standard Press-Schechter halo mass function, one gets a rate (Bird et al. 2016) that is independent of the primordial BH mass. However, for realistic extended mass functions, the abundance, size and evolution of primordial BH clusters is impacted by several effects: poissonian noise, seeds from heavy primordial BHs, primordial power spectrum enhancement, dynamical heating, etc. Those can either boost or suppress the merging rate and make it a rather complex and model-dependent process, subject to large uncertainties (see e.g., Clesse and García-Bellido ( Trashorras et al. (2021) for recent studies of primordial BH clustering). As an alternative of using uncertain theoretical predictions, one can instead infer from LIGO/Virgo observations an upper bound on the primordial BH merging rate of s % 1:2 Â 10 4 yr À1 Gpc À3 at m PBH % 2:5M ) and of s % 50yr À1 Gpc À3 at masses between 10 M and 50 M (Abbott et al. 2019a). The boost in primordial BH formation at the time of the QCD transition will induce a peak at the solar mass scale in any primordial BH model with an extended mass function (Byrnes et al. 2018;Carr et al. 2021). If one normalises the merging rates at the peak with the LIGO/ Virgo rates in the solar mass range, then one obtains an upper bound on the rate distribution ds dðln m 1 Þdðln m 2 Þ % 4 Â 10 3 Â qðm 1 Þqðm 2 Þ ðm 1 þ m 2 Þ 10=7 while being agnostic about the total primordial BH abundance, f PBH . Then, like for primordial binaries, one can obtain an upper limit on the merging rate for equal-mass sub-solar binaries in halos. Assuming qðm PBH Þ % 0:01, as suggested by microlensing surveys, one gets independently of the primordial BH mass. This is of the same order than the rates obtained with the theoretical prescriptions leading to Eq. (25), for a monochromatic distribution. One thus expects that Eq. (26) is a good approximation for a broad variety of primordial BH scenarios, for both sharp and wide mass distributions. One then obtains the required experimental strain sensitivity to detect one event per year, This GW signal is therefore typically lower than for primordial BH binaries formed in the early Universe. However, it is still debated which of the binary formation channel is dominant, especially if f PBH J0:1, i.e. if primordial BHs explain a significant or even the totality of the dark matter in the Universe.
While in the literature it is often assumed for simplicity that primordial BHs have a monochromatic mass function, this is not expected in a realistic scenario. Even in the limiting case of BH formation due to a sharp peak in the primordial power spectrum, these primordial BHs would have a relatively broad distribution, due to effects related to the critical collapse (Musco et al. 2009;Musco and Miller 2013).
For this reason and to be generic, we have also estimated in the above paragraph the rate distribution for the two possible binary formation mechanisms, without specifying the primordial BH mass function.

Exotic compact objects
Beyond the very well-known astrophysical compact objects, namely BHs and neutron stars, there are several candidates for stable (or long-lived) exotic compact objects that are composed of beyond the Standard Model particles (Giudice et al. 2016). For instance, they can be composed of beyond the Standard Model fermions, such as the gravitino in supergravity theories, giving rise to gravitino stars (Narain et al. 2006). Exotic compact objects can also be composed of bosons, such as moduli in string compactifications and supersymmetric theories . Depending on the mechanism that makes the compact object stable (or longlived), scalar field exotic compact objects have specific names such as Q-balls, boson stars, oscillatons, oscillons. There are also more exotic possibilities, such as gravastars (Mazur and Mottola 2004). Exotic compact objects can form binaries and emit GWs in the same way as BH and neutron star binaries do. During the early inspiral phase, the frequency of the emitted GWs is twice the orbital frequency. At the ISCO, the frequency for a binary system of two exotic compact objects with mass M and radius R is given by (Giudice et al. 2016) where C ¼ GM=R is the compactness of the exotic compact object. This expression is only slightly modified for a boson star binary with two different values of the masses. Note that for a BH the radius is given by the Schwarzschild radius R S ¼ 2GM, therefore C ¼ 1=2 is the maximum attainable value for the compactness.
The GW strain for a boson star binary formed by equal mass objects M can be estimated as done in Sect. 3.2.2 where D is the distance between the source and the observer. The exact waveform produced by the merger of two exotic compact objects is in general different from that of BHs and neutron stars and depends on its microphysics details (Giudice et al. 2016;Palenzuela et al. 2017). 10 Hence, the detection of GWs from an exotic compact object merger would give further valuable information about beyond the Standard Model physics.

Black-hole superradiance
This section focuses on GW emission from clouds of axions or axion-like particles created by the gravitational superradiance of BHs (Ternov et al. 1978;Zouros and Eardley 1979;Arvanitaki et al. 2010;Arvanitaki and Dubovsky 2011;Arvanitaki and Geraci 2013;Aggarwal et al. 2020;Detweiler 1980;Yoshino and Kodama 2014;Arvanitaki et al. 2015;Brito et al. 2015a, b). Superradiance is an enhanced radiation process that is associated with bosonic fields around rotating objects with dissipation. The event horizon of a spinning BH is one such example that provides conditions particularly suitable for superradiance (Arvanitaki et al. 2015). When the axion Compton wavelength, determined by the axion mass m a , is about the size of the BH, the axions can accumulate outside the BH event horizon and outside the ergosphere efficiently. The BH forms a gravitationally bound 'atom' with the axions, with different atomic 'levels' occupied by exponentially large numbers of axions. The primary candidate for GWs at high frequencies is the axion annihilation process (a þ a ! h), with frequencies of around 100 kHz for M BH JM . The GW frequency emitted by this process is twice the Compton frequency of the axion, i.e.
f ¼ 2 m a 10 À9 eV 10 6 Hz : The expected GW strain for this process is roughly (Arvanitaki and Geraci 2013) where a ¼ GM BH m a , l is the orbital angular momentum number of the axions that decay, D is the distance from the observer and \10 À3 denotes the fraction the BH mass accumulated in the axion cloud. The superradiance condition constrains a=l\0:5 (Arvanitaki and Dubovsky 2011). See Brito et al. (2015a) and Aggarwal et al. (2020) for more recent calculations of GW strain from BH superradiance leading to axion annihilations. If lighter spinning BHs exist (see Sect. 3.2.2), then this process can produce GWs of even higher frequencies. Note that this GW signal is predicted to be monochromatic and coherent (Arvanitaki et al. 2015), rendering it quite distinct from any other astrophysical or cosmological sources discussed here.
For heavier BHs and lighter axions, this process can produce GWs in the LIGO/ Vigro band. Additionally, these bosons created in BH superradiance may also undergo transition from one level to another, emitting a graviton in the process. Both these processes would produce GWs of lower frequencies and could be detectable in LIGO-Vigro and LISA, e.g., Brito  Finally, recently it has been postulated that axions might also decay into gravitons (a ! h h) (Sun and Zhang 2020). In such a process, the GW frequency would be half of the axion Compton frequency, i.e.
f ¼ 1 2 m a 10 À9 eV 10 6 Hz : The corresponding strain of the coherent signal has been calculated in Sun and Zhang (2020) to be where \10 À3 denotes the fraction the BH mass accumulated in the axion cloud.

Early Universe
We now turn to cosmological sources emitting GWs at cosmological distances, i.e., in the early Universe. For a summary of these sources see Fig. 2 and Table 3 in Appendix 1. In this case, the source is associated to an event in our cosmological history, triggered e.g., by the decreasing temperature T of the thermal bath, and typically occurs everywhere in the Universe at (approximately) the same time. This results in a stochastic background of GWs which is a superposition of GWs with different wave vectors. The total energy density of a GW background q GW R d log k dq GW =d log k ð Þ , with characteristic wavelengths well inside the horizon, decays as relativistic degrees of freedom with the expansion of the Universe, i.e. as q GW / a À4 . This implies that a GW background acts as an additional radiation field contributing to the background expansion rate of the Universe. Observables that can probe the background evolution of the Universe at some particular moment of its history, can therefore be used to constrain q GW at such moments. In particular, two events in cosmic history yield a precise measurement of the expansion rate of the Universe: BBN and photon decoupling of the CMB. An upper bound on the total energy density of a GW background present at the time of BBN and CMB decoupling can be therefore derived from the constraint on the amount of radiation tolerated at those cosmic epochs, when the Universe had a temperature of T BBN $ 0:1 MeV and T CMB $ 0:3 eV, respectively.
A constraint on the presence of 'extra' radiation is usually expressed in terms of an effective number of neutrinos species N eff after electron-positron annihilation. After electron-positron annihilation, the total number of Standard Model relativistic degrees of freedom was g Ã ðT\T e þ e À Þ ¼ 2 þ 7 4 N eff 4 11 À Á 4=3 , with N eff ¼ 3:046. As the radiation energy density for thermal degrees of freedom in the Universe is given by À Á 4=3 DN eff , with q c denoting the energy density in photons. Writing the fraction of GW energy density today 12 as q c ðTÞ , we obtain a constraint on the redshifted GW energy density today, in terms of the number of extra neutrino species (Caprini and Figueroa 2018) where we have inserted X rad;0 h 2 We recall that the above bound applies only to the total GW energy density, integrated over wavelengths way inside the Hubble radius (for super-horizon wavelengths, tensor modes do not propagate as a wave, and hence they do not affect the expansion rate of the Universe). Except for GW spectra with a very narrow peak of width Df ( f , the bound can be interpreted as a bound on the amplitude of a GW spectrum, X GW; 0 ðf Þ h 2 H . 5:6 Â 10 À6 DN eff , over a wide frequency range. The bound obviously applies only to GW backgrounds that are present before the physical mechanism (BBN or CMB decoupling) considered to infer the constraint on N eff takes place.
Constraints on N eff can be placed by BBN alone, and/or in combinations with CMB data. In particular, Cyburt et al. (2016) finds DN eff \0:2 at 95% confidence level. Eq. (36) then gives for a stochastic GW background produced before BBN, with wavelengths inside the Hubble radius at the onset of BBN, corresponding to present-day frequencies f ! 1:5 Â 10 À12 Hz. A similar bound can also be obtained from constraints on the Hubble rate at CMB decoupling (Smith et al. 2006;Sendra and Smith 2012;Pagano et al. 2016;Clarke et al. 2020). This translates into an upper bound on the amount of GWs, which extends to a greater frequency range than the BBN bound, down to f J10 À15 Hz. Since high-frequency GWs carry a lot of energy, X GW / f 3 S h , these bounds pose severe constraints on possible cosmological sources of high-frequency GWs.

Inflation
Under the standard assumption of scale invariance, the amplitude of the GWs produced during inflation is too small (X GW; 0 .10 À16 ) to be observable with current technology. 13 11 This is independent of whether the extra radiation is in a thermal state or not, as this is only a parametrization of the total energy density of the extra component, independently of its spectrum. 12 We write the current value of the Hubble parameter as H 0 ¼ h H Â 100 km s À1 Mpc À1 . Early Universe and late time observations report slightly different values for h H , see Bernal et al. (2016) for a discussion. For all our purposes, we will assume h H ¼ 0:7 when needed.
Various inflationary mechanisms have been studied in the literature that can produce a significantly blue-tilted GW signal (i.e., increasing towards higher frequency), or a localized bump at some given (momentum) scale, with a potentially visible amplitude. A number of these mechanisms have been explored in Bartolo et al. (2016b) with a focus on the LISA experiment, and therefore on a GW signal in the mHz range. However, these mechanisms can be easily extended to higher frequencies.
Assuming an approximately constant Hubble H parameter during inflation, a GW generated N Hubble times (e-folds) before the end of inflation with frequency H is redshifted to frequency f today, with where N CMB is the number of e-folds at which the CMB modes exited the horizon. The numerical value of N CMB depends logarithmically on the energy scale of inflation, which is bounded from above by the upper bound on the tensor-to-scalar ratio (Akrami et al. 2020), H . 6 Â 10 13 GeV. Saturating this bounds implies N CMB ' 60, and a peak at f ¼ 1 MHz then corresponds to the N ¼ 4:7, while the LIGO frequency f LIGO ¼ O 10 2 Hz ð Þcorresponds to about N ¼ 14. Such late stages of inflation are not accessible by electromagnetic probes. Bartolo et al. (2016b) discusses three broad categories: the presence of extra fields that are amplified in the later stages of inflation (so to affect only scales much smaller than the CMB ones); GW production in the effective field theory framework of broken spatial reparametrizations and GWs sourced by (large) scalar perturbations. In the following we will briefly summarize these three cases.
Extra-species Several mechanisms of particle production during inflation, with a consequent GW amplification, have been considered in the recent literature. Here, for definiteness, we discuss a specific mechanism in which a pseudo-scalar inflation / produces gauge fields via an axionic coupling / 4f a FF, where F lm is the gauge field strength,F lm its dual, and f a is the axion decay constant. The motion of the inflaton results in a large amplification of one of the two gauge field helicities. The produced gauge quanta in turn generate inflaton perturbations and GW via 2 ! 1 processes (Barnaby and Peloso 2011;Sorbo 2011). In particular, the spectrum of the sourced GWs is (Barnaby and Peloso 2011) where H is the Hubble rate. In this relation, H and _ / are evaluated when a given mode exits the horizon, and therefore the spectrum in Eq. (39) is in general scaledependent. In particular, in the n ) 1 regime, the GW amplitude grows exponentially with the speed of the inflaton, which in turn typically increases over the course of inflation in single-field inflation models. As a consequence, the spectrum in Eq. (39) is naturally blue. The growth of n is limited by the backreaction of the gauge fields on the inflaton. Within the limits of a perturbative description, n . 4:7 , GW amplitudes of X GW;0 ' 10 À10 can be obtained. Domcke et al. (2016); García-Bellido et al. (2016) explored the resulting spectrum for several inflaton potentials. In particular hill-top potentials are characterized by a very small speed close to the top (that is mapped to the early stages of observable inflation), and by a sudden increase of the speed at the very end of inflation. Interestingly, hill-top type potentials are naturally present (Peloso and Unal 2015) in models of multiple axions such as aligned axion inflation (Kim et al. 2005).
Effective field theory spatial reparametrizations The modification of the theory of gravity which underlines the inflationary physics can give rise to an extra production of GWs with a large amplitude (and blue tilt) rendering it accessible to highfrequency GW experiments. From the theoretical point of view, the effective field theory approach (Cheung et al. 2008) represents a powerful tool to provide a clear description of the relevant degrees of freedom at the energy scale of interest exploiting the power of symmetries and gives an accurate prediction of observational quantities. In the standard single-field effective field theory of inflation (Cheung et al. 2008) only time-translation symmetry (t ! t þ n 0 ) is broken according to the cosmological background expansion during inflation. However when space-reparameterization symmetry (x i ! x i þ n i ) is also broken (Bartolo et al. 2016a;Graef and Brandenberger 2015), scalar and tensors (GWs) acquire interesting features. In particular tensors can acquire a non-trivial mass m h and sound speed c T , making them potential targets for high-frequency detectors since in this case the spectrum gets enhanced on small scales. At the quadratic level in perturbations, in an effective field theory approach, the action for graviton fluctuations h ij around a conformally flat Friedmann-Lemaître-Robertson-Walker background can be expressed as in Cannone et al. (2015); Bartolo et al. (2016a); Ricciardone and Tasinato (2017): The corresponding tensor power spectrum and its related spectral tilt are: Hence, if the quantity m h =H is sufficiently large, we can get a blue tensor spectrum with no need to violate the null energy condition in the early Universe. Consequently X GW;0 $ X rad;0 P T is enhanced at high frequencies, making it a potential target for high-frequency GW detectors. The upper bound at on the spectrum at high frequencies is set by the observational BBN and CMB bounds, see Eq. (37). This scenario shows how GW detectors at high frequency might be useful to test the modification of gravity at very high-energy scales.
Second-order GW production from primordial scalar fluctuations In homogeneous and isotropic backgrounds, scalar, vector and tensor fluctuation modes decouple from each other at first order in perturbation theory. These modes can however source each other non-linearly, starting from second order. In particular, density perturbations can produce 'induced' (or 'secondary') GWs through a f þ f ! h process Mollerach et al. (2004); Ananda et al. (2007); Baumann et al. (2007) (see also Kohri and Terada 2018;Espinosa et al. 2018;Braglia et al. 2020). This production, which simply involves only gravity, is mostly effective when the modes re-enter the horizon after inflation. (Second order GWs are also produced in an early matter era, Inomata et al. 2019a, b.) The amplitude of this signal is quadratic in the scalar perturbations, and scale-invariant O 10 À5 À Á perturbations, as measured on large scales by the CMB, result in unobservable GWs due to too small amplitude. On the other hand, if the spectrum of scalar perturbations produced during inflation has a localized bump at some given scale (significantly smaller than the scales of CMB and the large scale structure), as required e.g., to obtain a sizable primordial BH abundance of some specific given mass, the height of the bump could be sufficiently high to produce a noticeable amount of GWs (Inomata et al. 2017;García-Bellido et al. 2017;Bartolo et al. 2019). The non-detection of the stochastic GW background can also be used to constrain fluctuations Byrnes et al. (2019); Inomata and Nakama (2019). The induced GWs have a frequency f Ã parametrically equal to the momentum k Ã and can hence be related to the e-fold N of horizon exit of the scalar perturbation through Eq. (38).
The precise amount of produced GWs depends on the statistics of the scalar perturbations (Nakama et al. 2017;García-Bellido et al. 2017;Cai et al. 2019;Unal 2019). A reasonable estimate is however obtained by simply looking at the scalar two-point function, where P f is the power spectrum (two-point function) of the gauge invariant scalar density fluctuations such that hf k f k 0 i / dðkþk 0 Þ k 3 P f ðkÞ. From this relation, the present value of the induced stochastic GW background is given by At the largest scales of our observable Universe, the density fluctuations are measured as P f ' 2 Á 10 À9 , resulting in X GW;0 $ Oð10 À22 Þ. Primordial BH limits are compatible with P f as large as .10 À2:5 at some (momentum) scale k Ã , in which case X GW;0 $ Oð10 À9 Þ.

(P)reheating
Preheating is an out-of-equilibrium production of particles due to non-perturbative effects (Traschen and Brandenberger 1990;Kofman et al. 1994;Shtanov et al. 1995;Kaiser 1996;Khlebnikov and Tkachev 1996;Prokopec and Roos 1997;Kaiser 1997;Kofman et al. 1997;Greene et al. 1997;Kaiser 1998), which takes place after inflation in many models of particle physics (see Allahverdi et al. 2010;Amin et al. 2014; Lozanov 2019 for reviews). After inflation, the interactions between the different fields may generate non-adiabatic time-dependent terms in the field equations of motion, which can give rise to an exponential growth of the field modes within certain bands of momenta. The field gradients generated during this stage can be an important source of primordial GWs, with the specific features of the GW spectra depending strongly on the considered scenario, see e.g., Khlebnikov and Tkachev (1997); García-Bellido (1998); Easther and Lim (2006); Easther et al. (2007); García-Bellido and Figueroa (2007); García-Bellido et al. (2008); Dufaux et al. (2007Dufaux et al. ( , 2009 ;Figueroa et al. (2011). If instabilities are caused by the field's own self-interactions, we refer to it as self-resonance, a scenario which will be discussed in more detail below. Here we consider instead a multi-field preheating scenario, in which a significant fraction of energy is successfully transferred from the inflationary sector to other fields. For illustrative process, let us consider a two-field scenario, in which the postinflationary oscillations of the inflaton excite a secondary massless species. More specifically, let us consider an inflaton with power-law potential  (Turner 1983). Let us now include a quadratic interaction term g 2 / 2 v 2 between the inflaton and a secondary massless scalar field v, where g is a dimensionless coupling constant. In this case, the driving postinflationary particle production mechanism is parametric resonance (Kofman et al. 1994Greene et al. 1997). In particular, if the so-called resonance parameter q H g 2 / 2 H =x 2 H obeys q H J1, the secondary field gets excited through a process of broad resonance, and the amplitude of the field modes grows exponentially inside a Bose-sphere of radius k.k H $ q 1=4 Ã x H . The GW spectrum produced during this process has a peak at approximately the frequency and amplitude (Figueroa and Torrenti 2017), where q H is the energy density at time t ¼ t H , g and d are two parameters that account for non-linear effects, and C is a constant that characterizes the strength of the resonance. The factor H ða H =a RD Þ 1À3w parametrizes the period between the end of inflation and the onset of the radiation dominated stage with a transitory effective equation of state w. If non-linear effects are ignored, the frequency and amplitude scale as f $ q 1=4 H and X GW;0 $ q À1=2 H respectively.
GWs can also be strongly produced if the species of the fields involved is different, or when the resonant phenomena driving preheating is different than parametric resonance. For example, GWs can be produced during the out-of-equilibrium excitation of fermions after inflation, for both spin-1/2 (Enqvist et al. 2012;Figueroa and Meriniemi 2013;Figueroa 2014) and spin-3/2 Benakli et al. (2019) fields. Similarly, GWs can also be generated when the decay products are (Abelian and non-Abelian) gauge fields. For example, the gauge fields can be coupled to a complex scalar field via a covariant derivative like in Dufaux et al. (2010); Figueroa et al. (2016); Tranberg et al. (2018), or to a pseudo-scalar field via an axial coupling as in Adshead et al. (2018Adshead et al. ( , 2020aAdshead et al. ( , 2020b. Preheating can be remarkably efficient in the second case, and the GW amplitude can scale up to X GW $ O 10 À6 À 10 À7 À Á for certain coupling strengths, see Adshead et al. (2020aAdshead et al. ( , 2020b for more details. Production of GWs during preheating with non-minimal couplings to the scalar curvature has also been explored in Fu et al. (2018). Finally, the stochastic background of GWs from preheating may develop anisotropies if the inflaton is coupled to a secondary light scalar field, see Bethke et al. (2013Bethke et al. ( , 2014. Oscillon production Oscillons are long-lived compact objects (Gleiser 1994) that can be formed in the early Universe in a variety of post-inflationary scenarios which involve a preheating-like phase (Amin and Shirokoff 2010;Amin et al. , 2012Zhou et al. 2013;Amin 2013;Lozanov and Amin 2014;Antusch et al. 2017Antusch et al. , 2018aLozanov and Amin 2018;Amin et al. 2018;Antusch et al. 2019;Sang and Huang 2019;Lozanov and Amin 2019;Fodor 2019;Hiramatsu et al. 2021). Their dynamics is a possible source of GW production. Oscillons are pseudo-solitonic solutions of real scalar field theories: their existence is due to attractive self-interactions of the scalar field that balance the outward pressure. 14 The real scalar field self-interactions are attractive if the scalar potential is shallower than quadratic at least on one side with respect to the minimum. Oscillons can be thought off as bubbles in which the scalar field is undergoing large oscillations that probe the non-linear part of the potential, while outside the scalar field is oscillating with a very small amplitude around the minimum of the potential.
As discussed in the previous section, during preheating the quantum fluctuations of the scalar field are amplified due to a resonance process. The Universe ends up in a very inhomogeneous phase in which the inflaton (or any other scalar field that produces preheating) is fragmented and there are large fluctuations in the energy density. At this point, if the field is subject to attractive self-interactions, the inhomogeneities can clump and form oscillons. While clumping oscillons deviate significantly from being spherically symmetric, therefore their dynamics produce GWs. After many oscillations of the scalar field they tend to become spherically symmetric and GW production stops. However, during their entire lifetime oscillons can produce GWs also due to the interactions and collisions among each other (Helfer et al. 2019). Oscillons are very long-lived: their lifetime is model-dependent but typically J10 4 =m (Gleiser and Sicilia 2008;Amin and Shirokoff 2010;Amin et al. , 2012Salmi and Hindmarsh 2012;Saffin et al. 2014;Antusch et al. 2019;Gleiser and Krackow 2019;Zhang et al. 2020), where m is the mass of the scalar field. Oscillons eventually decay through classical (Segur and Kruskal 1987) or quantum radiation (Hertzberg 2010).
The peak of the GW spectrum at production is centered slightly below the value of the mass of the field, that typically correspond to a frequency today well above the LIGO range 15 (Zhou et al. 2013;Antusch et al. 2018a;Lozanov and Amin 2019). In a typical situation, an oscillating massive scalar field forming oscillons quickly comes to dominate the energy density of the Universe until the perturbative decay of the field itself. For the simplest case of a gravitationally coupled massive field that starts oscillating at H ' m and decays at H $ m 3 =M 2 p ) the frequency today can be estimated as f ' X m 10 12 GeV 5=6 10 6 Hz ; where the factor X which is typically in the range X ' ð10 À 10 3 Þ is due to the unknown precise time of GW production and can be obtained in concrete models through lattice simulations: the equality would hold if GWs were produced immediately when the scalar field starts oscillating. 16 On the other hand, the later GWs are produced, the less the frequency is red-shifted and the larger is X. The maximum value of today's amplitude for these processes, inferred from numerical simulations, is in the range X GW;0 ' ð10 À13 À 10 À10 Þ (Antusch et al. 2017(Antusch et al. , 2018aAmin et al. 2018), see Dufaux et al. (2007) for a discussion on how to compute the GW amplitude. Depending on the model, gravitational effects can become important and play a crucial role for the existence/stability of the solution (Seidel and Suen 1991). In particular the requirement that the potential must be shallower than quadratic is no longer necessary, as the attractive force is provided by gravity (Urena-Lopez et al. 2002). In this case oscillons are equivalent to oscillatons, see Sect  Kitajima et al. (2018) for models that lead to a GW peak at lower frequencies. 16 This rough estimate assumes that the field starts oscillating when H ' m. Since the potential contains self-interactions, assuming that the field starts at rest, the actual requirement for the start of the oscillations is V 00 ð/ in Þ $ H, where / in is the initial value of the field. Please note that if the field is the inflaton, the initial conditions are different from those assumed in Eq. (45) and therefore this estimate does not necessarily hold, see e.g., Antusch et al. (2017).

Cosmic gravitational microwave background
The hot thermal plasma of the early Universe acts as a source of GWs, which, similarly to the relic photons of the CMB, peak in the $ 100 GHz range today. The spectrum of this signal is determined by the particle content and the maximum temperature T max reached by the thermal plasma in the Universe history (Ghiglieri and Laine 2015;Ghiglieri et al. 2020;Ringwald et al. 2021). Ignoring the dependence on the number of relativistic degrees of freedom, the energy density in GWs per logarithmic frequency interval can then be written as where T 0 is the temperature of the CMB today, whileĝðT max ; 2pf =T 0 Þ encodes the sources of GW production in the thermal plasma: it is dominated by long range hydrodynamic fluctuations at 2pf \T 0 and by quasi-particle excitations in the plasma at 2pf $ T 0 , see Ghiglieri and Laine (2015); Ghiglieri et al. (2020); Ringwald et al. (2021) for more details. The peak frequency of X GW;0 ðf Þ in Eq. (46) is in the (1-100) GHz range today and depends on the number of entropic relativistic degrees of freedom g Ãs ðT ¼ T max Þ. The peak of X GW;0 ðf Þ approaches the BBN bound if T max $ M p . The CMB constraints on the tensor-to-scalar ratio however constrain the maximal reheating temperature to T max \10 16 GeV (Akrami et al. 2020) under the assumption of slow-roll inflation and instantaneous reheating. Therefore the detection of the cosmic gravitational microwave background corresponding to T max [ 10 16 GeV would rule out slow-roll inflation as a viable pre hot Big Bang scenario. Note that since at leading order X GW;0 ðf Þ scales linearly with T max and the peak frequency depends on g Ãs ðT max Þ, the detection of the peak of the cosmic gravitational microwave background would determine both T max and g Ãs ðT max Þ, see Ringwald et al. (2021) for more details.

Phase transitions
A first order phase transition in the early Universe proceeds by the nucleation of bubbles of the low-temperature phase as the Universe cools below the critical temperature (Steinhardt 1982;Hogan 1983). Due to the higher pressure inside, the bubbles expand and collide, and the stable phase takes over. The process disturbs the fluid, generating shear stresses and hence GWs (Witten 1984;Hogan 1986). As the perturbations are mostly compression waves, they can be described as sound waves, and their collisions are the main source of GWs (Hindmarsh et al. 2014(Hindmarsh et al. , 2015(Hindmarsh et al. , 2017a. The peak frequency of an acoustic contribution to a relic GW background from a strong first order transition is controlled by the temperature of the transition T Ã , and the mean bubble separation R Ã . 17 Numerical simulations show for wall speeds not too close to the speed of sound that (Hindmarsh et al. 2017a) where H Ã is the Hubble rate at nucleation. The theoretical expectation is that 1.ðH Ã R Ã Þ À1 .10 4 . The intensity depends on H Ã R Ã , on the fraction of the energy density of the Universe which is converted into kinetic energy K and on the lifetime of the source, which can last for up to a Hubble time. Denoting the lifetime of the velocity perturbations by s v , the peak GW amplitude can be estimated as Hindmarsh et al. (2015); Guo et al. (2021) whereX GW is a simulation factor and s is the life time of the sound waves. Numerical simulations indicateX GW ¼ O 10 À2 ð Þ: Hence, X GW; 0 . 10 À7 today, with the upper bound reached only if most of the energy available in the phase transition is turned into kinetic energy. This is only possible if there is significant supercooling.
The calculation of the kinetic energy fraction and the mean bubble separation requires a knowledge of the free energy density f ðT; /Þ, a function of the temperature and the scalar field (or fields) / whose expectation value determines the phase. If the underlying quantum theory is weakly coupled, and the scalar particle corresponding to / is light compared to the masses gained by gauge bosons in the phase transition, this is easily calculated, and shows that first order transitions are generic in gauge theories in this limit (Kirzhnits 1972;Kirzhnits and Linde 1976), meaning that there is a temperature range in which there are two minima of the free energy as a function of /. The critical temperature is defined as the temperature at which the two minima are degenerate, separated by a local maximum.
The key parameters to be extracted from the underlying theory, besides the critical temperature T c , are the nucleation rate b, the strength parameter a and the bubble wall speed v w . The nucleation rate parameter b ¼ d log p=dt, where p is the bubble nucleation rate per unit volume, is calculable from f ðT; /Þ through an application of homogeneous nucleation theory (Langer 1969) to high-temperature fields (Linde 1983). This calculation also gives T Ã as the temperature at which the volume-averaged bubble nucleation rate peaks. The strength parameter is roughly, but not precisely, one quarter of the latent heat divided by the thermal energy (see Hindmarsh and Hijazi 2019 for a more precise definition) at the nucleation temperature, and also follows from knowing f ðT; /Þ. The wall speed is a nonequilibrium quantity, which cannot be extracted from the free energy alone, and is rather difficult to calculate accurately (see Dorsch et al. 2018; Laurent and Cline 2020 and references therein). In terms of these parameters, it can be shown that Enqvist et al. (1992) R Ã $ v w =b. The kinetic energy fraction K, can be estimated from the self-similar hydrodynamic flow set up around an isolated expanding bubble, whose solution can be found as a function of the latent heat and bubble wall velocity by a simple one-dimensional integration (Turner et al. 1992;Espinosa et al. 2010;Hindmarsh and Hijazi 2019) 18 and typically ranges between K ¼ 1 À 10 À6 .
Current projected sensitivities for the Einstein Telescope and the Cosmic Explorer can probe a cosmological first order transition occurring at a temperature that is at most a few hundred TeV assuming a modest amount of supercooling (Abbott et al. 2017a;Punturo et al. 2010;Hild et al. 2011) (i.e., when T Ã $ T C and ðR Ã H Ã Þ À1 J100). Recently, there has been much interest in high scale transitions motivated by U(1) BÀL breaking for leptogenesis and the seesaw scenario (Jinno and Takimoto 2017;Marzo et al. 2019;Brdar et al. 2019b;Okada and Seto 2018;Hasegawa et al. 2019;Okada et al. 2021), as well as multi-step grand unification breaking patterns such as a Pati-Salam model (Croon et al. 2019;Greljo et al. 2020;Brdar et al. 2019a;Huang et al. 2020). However, in both cases it is more natural to motivate significantly higher scale transitions and to most naturally probe first order transitions at these scales one requires a detector sensitive to frequencies in the range f ' ð10 À3 À 10 3 Þ GHz. Finally, it has been shown that zero-temperature phase transitions are likely to occurr in string models that include warped throats (García García et al. 2018). These lead to a GW spectrum whose peak falls in the high-frequency range for small values of the compactification volume and large values of the gravitational warp factor associated with the throat in which the phase transition occurs.

Topological defects
Cosmic strings are one-dimensional topological defect solutions which may have formed after a phase transitions in the early Universe, if the first homotopy group of the vacuum manifold associated with the symmetry breaking is non-trivial (Kibble 1976;Jeannerot et al. 2003). They can also be fundamental strings from string theory, formed for instance at the end of brane inflation (Dvali and Vilenkin 2004;Copeland et al. 2004), and stretched to cosmological scales. The energy per unit length of a string is l $ g 2 , with g the characteristic energy scale (in the case of topological strings, it is the energy scale of the phase transition). Typically, the tension of the strings is characterized by the dimensionless combination Gl $ ðg=M p Þ 2 , e.g., the current upper bound from the CMB is Gl . 10 À7 , whereas GW searches in pulsar timing arrays constrain the tension to Gl . 10 À11 . Cosmic strings are energetic objects that move at relativistic speeds. The combination of these two factors immediately suggests that strings should be a powerful source of GWs.
Whenever cosmic strings are formed in the early Universe, their dynamics drive them rather rapidly into an attractor solution, characterized by their energy density maintaining always a fixed fraction of the background energy density of the Universe. This is known as the 'scaling' regime. During this regime, strings will collide, possibly exchanging 'partners' and reconnecting afterwards. This is known as 'intercommutation'. For topological strings the intercommutation probability is P ¼ 1, whereas P\1 is characteristic in cosmic superstrings networks. Closed string configurations -loops -are consequently formed when a string selfintersects, or two strings cross. Loops smaller than the horizon decouple from the string network and oscillate under their own tension, which results in the emission of gravitational radiation (eventually leading to the decay of the loop). The relativistic nature of strings typically leads to the formation of cusps, corresponding to points where the string momentarily moves at the speed of light (Turok 1984). Furthermore, the intersections of strings generates discontinuities on their tangent vector known as kinks. All loops are typically expected to contain cusps and kinks, both of which generate GW bursts Vilenkin 2000, 2001). Hence, a network of cosmic (super-)strings formed in the early Universe is expected to radiate GWs throughout the entire cosmological history, producing a stochastic background of GWs from the superposition of many uncorrelated bursts. 19 A network of cosmic strings contains therefore, at every moment of its evolution (once in scaling), sub-horizon loops, and long strings that stretch across a Hubble volume. The latter are either infinite strings or in the form of super-horizon loops, and are also expected to emit GWs. However, the dominant contribution is generically that produced by the superposition of radiation from many sub-horizon loops along each line of sight.
The power emitted into gravitational radiation by an isolated loop of length l can be calculated using the standard formulae in the weak gravity regime (Weinberg 1972). More explicitly, we can assume that, on average, the total power emitted by a loop is given by P 1Loop ¼ C Â ðGlÞ Â l, with C is a dimensionless constant (independent of the size and shape of the loops). Estimates from simple loops (Vachaspati and Vilenkin 1985;Burden 1985;Garfinkle and Vachaspati 1987), as well as results from Nambu-Goto simulations (Blanco-Pillado and Olum 2017), suggest that C ¼ 50. The GW radiation is only emitted at discrete frequencies by each loop, x n ¼ 2pn=T, where T ¼ l=2 is the period of the loop, and n ¼ 1; 2; 3; . . . We can write P 1Loop ¼ Gl 2 P n P n , with P n characterizing the power emitted at each frequency x n for a particular loop, depending on whether the loop contains cusps, kinks, and kink-kink collisions (Burden 1985;Allen and Shellard 1992). Namely, it can be shown that for large n, P n ¼ C fðqÞ n Àq , where fðqÞ is the Riemann zeta function. The latter is simply introduced as a normalization factor, to enforce the total power of the loop to be equal to C ¼ P n P n . The parameter q takes the values 4/3, 5/3, or 2 depending on whether the emission is dominated by cusps, kinks or kink-kink collision respectively, see e.g., Vachaspati and Vilenkin (1985); Binetruy et al. (2009);Auclair et al. (2020).
The resulting GW background from the stochastic background emitted by the loops chopped-off along the radiation domination period is characterized by a scaleinvariant energy density spectrum, spanning over many frequency decades. The high-frequency cut-off of this plateau is determined by the temperature of the thermal bath at formation of the string network, with T max . 10 16 GeV implying a turn-over frequency of f D . 10 9 GeV (Gouttenoire et al. 2020). The amplitude of this plateau is given by (Auclair et al. 2020) X plateau GW;0 ðf Þ % 8:04 X rad;0 In particular, this does not depend on the exact form of the loop's power spectrum, nor on whether the GW emission is dominated by cusps or kinks, but rather depends only on the total GW radiation emitted by the loops. The amplitude in Eq. (49) indicates that the stochastic GW background from cosmic strings can be rather large. 20 Moreover, if the phase transition responsible for cosmic string formation is embedded into a larger grand unified group, then, depending on the structure of that larger group, cosmic strings may be only metastable, decaying via the (exponentially suppressed) production of monopoles (Vilenkin 1982;Voloshin 2008, 2010;Leblond et al. 2009). In this case, the low-frequency end of the spectrum, corresponding to GW emission at later times, will be suppressed and the signal may only be detectable in the high-frequency range (Leblond et al. 2009;Dror et al. 2020;Buchmuller et al. 2020a, b). In this case, the string tension is only bounded by the BBN bound, Gl . 10 À4 and the scale-invariant part of the spectrum may extend from 10 3 Hz (LIGO constraint) up to 10 9 Hz (network formation).
Finally, let us recall that long strings (infinite and super-horizon loops) also radiate GWs. One contribution to this signal is given by the GWs emitted around the horizon scale at each moment of cosmic history, as the network energy-momentum tensor adapts itself to maintaining scaling (Krauss 1992;Jones-Smith et al. 2008;Fenu et al. 2009;Figueroa et al. , 2020. This background is actually expected to be emitted by any network of cosmic defects in scaling, independently of the topology and origin of the defects . It represents therefore an irreducible background generated by any type of defect network. In the case of cosmic string networks modelled by the Nambu-Goto approximation (where the thickness of the string is taken to be zero), this background represents a very subdominant signal compared to the GW background emitted from the loops. In the case of field theory strings (for which simulations to date indicate the absence of 'stable' loops), it is instead the only GW signal (and hence the dominant one) emitted by the network. 20 Important remark: as the characteristic width d $ 1=g of a cosmic string is generally much smaller than the horizon scale, it is commonly assumed that strings can be described by the Nambu-Goto (NG) action, which is the leading-order approximation when the curvature scale of the strings is much larger than their thickness. The above plateau amplitude Eq. (49) applies only for the case of NG strings. For NG strings to reach scaling, the GW emission of the loops is actually crucial, as it is the loss of loops from the network that guarantees the scaling. The loops need therefore to decay in some way (this is precisely what the GW emission takes care of), so that their energy is not accounted anymore as part of the string network. However, in field theory simulations of string networks (Vincent et al. 1998;Hindmarsh et al. 2009;Daverio et al. 2016;Hindmarsh et al. 2017b), the network of infinite strings reaches a scaling regime thanks to energy loss into classical radiation of the fields involved in the simulations. The simulations show the presence of extensive massive radiation being emitted, and that the loops formed decay within a Hubble time. This intriguing discrepancy has been under debate for the last $ 20 years, but the origin of this massive radiation in the lattice simulations is not understood. The GW energy density spectrum of this irreducible background from long strings is predicted to be exactly scale-invariant, for the modes emitted during radiation domination . The power spectrum of this background therefore mimics the spectral shape of the dominant signal from the loop decay, but with a smaller amplitude. The amplitude of the irreducible GW background from string networks depends ultimately on the fine details of the so-called unequal-timecorrelator of the network's energy-momentum tensor. This correlator, however, can be obtained only accurately from sufficiently large scale lattice simulations. In the case of global defects, the scale-invariant GW power spectrum has been obtained numerically with massively parallelized lattice field theory simulations for global strings 21 as ) Despite the numerical prefactor being much larger than unity, the quadratic scaling proportional to ðGlÞ 2 suppresses significantly the amplitude, see e.g., Buchmüller et al. (2013) for a comparison among GW signals emitted from the same string network. This amplitude is clearly subdominant when compared to the amplitude of the GW signal from the loops, which scales as ðGlÞ 1=2 , see Eq. (49). A proper assessment of the power spectrum of this stochastic background requires further results not available yet; namely, lattice simulations of cosmic networks with a larger dynamical range. Finally, we point out that since the irreducible GW emission described before is expected from any scaling defect network, global texture networks also emit a GW background due to their self-ordering during scaling (Jones- Smith et al. 2008;Fenu et al. 2009;Giblin et al. 2012;Figueroa et al. , 2020. Textures are formed when the second (or higher) homotopy group of the vacuum manifold is non-trivial (Vilenkin and Shellard 2000). One can achieve such condition in either case of the symmetry breaking of a global or a gauge group. In the case of a global group, the GW spectrum is scale invariant for radiation domination (Jones- Smith et al. 2008;Fenu et al. 2009;, and exhibits a peak at the horizon today for matter domination (Figueroa et al. 2020). For a gauge group the GW spectrum is suppressed at low-frequencies as the massive gauge boson can prevent self-ordering as the gauge field can cancel the gradient field at large scales. The peak frequency and amplitude of a gauge texture is therefore set by the symmetry breaking scale v (Dror et al. 2020 Given the frequency and amplitude both increase with v, it is only reasonable to consider high-frequency signals.

Evaporating primordial black holes
Section 3.2.2 discussed the GW signal emitted by primordial BHs that merge in the late Universe. Note that light primordial BHs (with mass smaller than 10 11 kg) that evaporate before BBN could produce a stochastic spectrum of GWs by merging and scattering (Dolgov and Ejlli 2011). The typical frequency for such a source is at f J GHz. Here we consider yet another source of GWs tied to primordial BHs, namely the emission of gravitons as part of the Hawking radiation of primordial BHs. This is particularly relevant for light primordial BHs which evaporated before BBN and for primordial BHs which have life time ranging from BBN until today (10 11 kg . m PBH . 10 14 kg). The graviton emission from a population of primordial BHs induces a stochastic background of GWs (Anantua et al. 2009;Dong et al. 2016) that peaks at very high frequencies, between f $ 10 13 Hz and 10 22 Hz. The shape and amplitude of the GW frequency spectrum depends on multiple factors, such as primordial BH abundance at formation, their mass spectrum, their eventual spin, the number of degrees of freedom in the particle theory. Roughly speaking, due to the redshift of GW amplitude and frequency, the observed GW spectrum is dominated by the latest stages of the primordial BH evolution and the frequency is hence set by the evaporation time (and hence the initial mass) of the primordial BH.
Taking into account the limits on the primordial BH abundance from BBN and extra-galactic background radiation, the maximum amplitude can be up to X GW;0 % 10 À7:5 for primordial BHs evaporating just before BBN, m PBH . 10 11 kg. For heavier ones that might have not yet totally evaporated today, 10 11 kg . m PBH . 10 14 kg, it can be up to X GW;0 % 10 À6:5 (Dong et al. 2016), with a spectrum peaked at frequencies between 10 18 Hz and 10 22 Hz.
The most interesting case are possibly much lighter primordial BHs which would have completely evaporated well before BBN, e.g., if they are produced at an energy close to the grand unification scale, E $ 10 15 GeV. Because the primordial BH density decreases like / 1=a 3 and the radiation density is / 1=a 4 , such early decaying primordial BHs can be very abundant in the early Universe, leading to an early matter dominated phase. GWs produced in their decay could then constitute a sizable fraction of the subsequent radiation dominated phase, limited only by the BBN and CMB constraints (see Eq. (37)). For primordial BHs from the rand unification scale, the GW frequency spectrum has a peak around 10 15 Hz and can reach an amplitude X GW;0 $ 10 À8 for number of degrees of freedom n dof $ 10 3 (Anantua et al. 2009). It might therefore be in the range of detection of future instruments.

Miscellaneous
In this section we summarize a few more ideas that have been shown to be relevant for high-frequency GW production but require more exotic setups to be realized and do not find place in the classification proposed above.
• Brane-world scenarios: the brane-world scenario (Rubakov and Shaposhnikov 1983) proposes that the very weak force of gravity in our ð3 þ 1Þ-dimensional Universe is only a part of the full strength gravity which is felt in the fifth dimension at a level commensurate with the other forces. This scenario suggests that two ð3 þ 1Þ-dimensional branes-one of which represents our 4-dimensional Universe-are separated in a fifth dimension by a small distance (Randall and Sundrum 1999;Maartens and Koyama 2010). If violent gravitational events-such as BH mergers-take place on the 'shadow' brane, which is in close proximity to our brane -they would excite oscillations not only in the shadow brane but also in the fifth dimensional brane separation, leading to GW production on our visible brane Clarkson and Seahra 2007). • Pre-Big Bang cosmology: the pre-Big Bang scenario provides an alternative to cosmological inflation to provide the initial conditions for the hot Big Bang theory. The scenario exploits the fundamental symmetries of string theory to build a model in which the Universe starts in a cold and empty state in the infinite past and moves towards a state of high curvature through an accelerated expansion Veneziano 2003, 2016). The state of high curvature corresponds to a region in the parameter space in which the theory is strongly coupled. It is then assumed that the strongly coupled theory is able to match this initial accelerated expansion to the usual hot Big Bang cosmology. Interestingly, this scenario predicts a blue spectrum of GWs, with a peak at high frequency (Brustein et al. 1996). • Quintessential inflation: if the inflationary epoch is followed by a phase in which the equation of state is stiffer than radiation, the stochastic spectrum of GWs features a growth at high frequency, followed by a sharp cutoff (Giovannini 1999). This kind of behaviour is expected in quintessential models of inflation, such as the one investigated in Peebles and Vilenkin (1999). The position of the peak depends very weakly on the number of minimally coupled scalar fields of the model, but it is independent of the final curvature at the end of inflation. Therefore, it is located at 100 GHz and cannot be moved around. The amplitude of the GW spectrum can become very large: in Giovannini (1999) the authors present a choice of the parameters such that X GW;0 ' 10 À6 at the peak. • Magnetars: magnetars are neutron stars that feature very large surface magnetic fields $ 10 9 À 10 11 T. Wen et al. (2017) suggests that gamma-ray bursts produced by the magnetar itself or by a companion body forming a binary system, and interacting with the surface magnetic field of the magnetar could be the source of high-frequency GW, with frequency around 10 20 Hz and energy density up to X GW;0 $ 10 À6 . • Reheating: Ema et al. (2020) shows that there exists a model-independent contribution to the stochastic GW background, due to the oscillations of the inflaton (or any other scalar field that oscillates around its minimum after inflation) around the minimum of its potential during preheating. These oscillations provide a driving force in the equation of motion for the tensor modes, leading to GW production at high frequency J10 5 Hz. The amplitude of this signal is bound to be quite small: in Ema et al. (2020) the authors present a choice of parameters such that X GW;0 .10 À21 . • Thermal gravitational noise of the Sun: the thermal motion of charged particles in the plasma of protons and electrons present at the core of stars leads to the production of a gravitational radiation noise (Weinberg 1972;Bisnovatyi-Kogan and Rudenko 2004). The frequency of this radiation is determined by the frequency of the collisions and, in the case of the Sun, it falls in the range ð10 12 À 10 18 Þ Hz. • Plasma instabilities: in Servin and Brodin (2003) the authors studied the interaction of electromagnetic waves and GWs in a magnetised plasma. In the high-frequency regime, a circularly polarized electromagnetic wave travelling parallel to the background magnetic field present in a plasma generates a GWs with the same frequency of the electromagnetic wave.
In summary, a broad range of theoretically well motivated extensions of the Standard Model of particle physics predict the existence of GW sources at different stages in the evolution of our Universe. The corresponding GW frequency range and amplitudes are summarized in Fig. 2 and in Fig. 1 as well as in Table 3 and Table 2 in Appendix 1.

Detection of gravitational waves at high frequencies
After the first detection of GWs at frequencies in the range (0.1-2.0) kHz (Abbott et al. 2019b), the expansion into other frequency bands is a natural next step-as it was in the 1950s when radio, X-ray and UV observations became possible with new technology. As detailed in the previous section, many exciting questions in astronomy, cosmology and fundamental physics are tied to GW signals with frequencies (far) above the capabilities of current detectors or their upgrades. Figures 1 and 2 give an impression of the range of GW amplitudes expected for various coherent and stochastic sources. Even GW upper limits with no known source targets at the time of publication of this white paper may be valuable in restricting physical theories.
In this section, we will investigate the experimental possibilities for the detection of high-frequency GWs. First, we will give an overview of current GW detectors and their limitations, followed by the introduction of several concepts for the detection of high-frequency signals. Depending on the detector concept and the targeted sources, we quote detector sensitivities in terms of strain ffiffiffiffi ffi S n p or in terms of the dimensionless quantities h c;n (for stochastic signals) and h 0;n;mono (for monochromatic signals), see Sect. 2 for details. Careful consideration of operation and bandwidth is needed to convert between these quantities.

Laser interferometers and resonant mass detectors and their limitations
The first GWs were detected by the Advanced LIGO (Abbott et al. 2016a) detectors in the US and the Advanced Virgo detector in Italy (Acernese et al. 2015). In early 2020, the Japanese KAGRA detector (Aso et al. 2013) joined LIGO's third observing run. These detectors are all based on the principle of a Michelson interferometer, using large suspended mirrors with several kilometers distance between them. Several other detectors are in the design phase. These detectors typically have their peak sensitivity at frequencies of a few hundred Hz.
However, some future detectors are designed to particularly expand the detection band towards either low or high frequencies. To expand the detection band of Earthbound interferometers to frequencies below 10 Hz, cryogenically cooled mirrors, large beam diameters and operation underground are considered (Abernathy et al. 2011;Adhikari et al. 2020). LISA, also based on laser interferometry, is a planned, satellite-based detector to increase the arm length beyond the possibilities on Earth and to reduce environmental noise sources such as seismic noise (Amaro-Seoane et al. 2017). LISA will have its peak-sensitivity in the mHz range. To increase interferometer sensitivity towards higher frequencies, options are an increase of laser power and/or resonant operation. The planned Australian NEMO detector will be targeting frequencies of up to several kHz, see Sect. 4.1.1 below. 22 While increasing the arm-length of an interferometer increases strain signal in some frequency band, longer arms are only really beneficial as long as the GW wavelength is longer than the interferometer arms. For significantly higher frequencies (MHz) interferometers with arm-lengths of meters are more suitable, but are of course at the same time limited by the smaller strain sensitivity achievable with shorter arms. This constitutes the main limitation of laser interferometers, used as direct strain meters, towards higher GW signal frequencies.
A concept to detect GWs which existed prior to the interferometers are resonant bar detectors, initially proposed and built by Joseph Weber in the 1960s. Their modern successors, resonant spheres, have peak sensitivities at several kHz. In Sect. 4.1.3, we will give a summary of these resonant spheres.

Laser interferometers: Neutron Star Extreme Matter Observatory (NEMO)
The first detection of a binary neutron star merger in 2017 (Abbott et al. 2017b) has increased the interest in the development of GW detectors with sensitivity in the few kHz regime which will be capable of detecting the merger and ringdown part of the waveform (Martynov et al. 2019). It is expected that such detectors will need to have strain sensitivities approaching ffiffiffiffi ffi S n p ' 10 À24 Hz À1=2 in the range (1-4) kHz for events that are likely to occur a few times per year. This sensitivity should be achieved by the third generation terrestrial GW detectors that are anticipated to come online in the later half of the 2030s (Abbott et al. 2017a;Punturo et al. 2010). The Australian GW community is currently exploring the feasibility of a new detector, 'NEMO', dedicated to detecting this merger phase and the following ringdown as well as testing third generation technology on a smaller scale (Ackley et al. 2020;Bailes et al. 2019;Adya et al. 2020). The planned sensitivity of this detector would reach ffiffiffiffi ffi S n p ' 10 À24 Hz À1=2 in the range (1-2.5) kHz (Ackley et al. 2020). This detector will work in collaboration with the existing second generation GW detector network that will provide sky localization for electromagnetic followup.
The dominant high-frequency noise source for interferometric GW detectors is quantum phase noise or shot noise as it is otherwise called. The magnitude of this noise source is inversely proportional to the square of the product of the circulating power incident on the test masses and the length of the arms of the detector. This generally necessitates extremely high powers in the arms of the interferometers (% 5 MW in the case of NEMO). Such high circulating powers lead to technical issues such as parametric and tilt instabilities and thermal induced distortions. These issues can be challenging to deal with, however a dedicated high-frequency detector promises to make this easier. This is because low-frequency sensitivity limits the actuation that can be applied to the test masses to correct instabilities and distortions. Further, relaxing the low-frequency sensitivity relaxes the requirements on seismic isolation and test mass suspension systems that can significantly reduce the cost.

Interferometers up to 100 MHz
As was first pointed out by Mizuno (1995), in laser interferometers the overall stored energy in the form of circulating laser power sets a limit on the achievable sensitivity and bandwidth, which is a consequence of the quantum Cramér-Rao bound. For a given laser power, higher bandwidth needs to be traded in for an increase in sensitivity. While opto-mechanical resonances can be introduced in the signal response of interferometers to shape the sensitivity curve for specific frequencies (Somiya et al. 2016;Korobko et al. 2018), it appears unlikely that the stored laser power can be further increased by several orders of magnitude. Therefore, broadband interferometric detectors reaching into the MHz detection range (while maintaining LIGO or Virgo-level strain sensitivity) seem not to be a viable option when taking also the arm-length argument from above into account.
Nevertheless there are three notable efforts (two existing and one under construction) of laser interferometers in the MHz range, which currently set the best experimental upper limits on GWs in their respective frequency bands.
One option is to build kHz-bandwidth interferometric detectors that are centered around much higher frequencies. Akutsu et al. (2008) have published upper limits from such a system working at 100 MHz. The detector used a synchronous recycling architecture based on a resonant recycling cavity of dimension 75 cm and a Nd:YAG laser with a power output of 0.5 W. The limit on stochastic GW signals was reported to be ffiffiffiffi ffi S n p $ 10 À16 Hz À1=2 , placing a bound of h c;sto . 7 Â 10 À14 . A study of the potential of this technique  showed that a sensitivity of 10 À20 Hz À1=2 is possible at 100 MHz with a bandwidth of 2 kHz, but the sensitivity decreases with increasing frequency and is not competitive above 1 GHz.
The sensitivity of a single instrument can be surpassed by correlating two colocated instruments in the case of searching for stochastic signals from GWs or other sources. The Holometer experiment at Fermilab consists of two co-located power recycled Michelson interferometers with 40-meter long arms. While their primary research target has been signatures of quantization of spacetime, they reach a sensitivity of 10 À21 Hz À1=2 approximately in the band (1-13) MHz (Chou et al. 2017). Using a 704-hr dataset from the Holometer experiment, the authors of Martinez and Kamai (2020) concluded that there are no identifiable harmonic sources such as cosmic string loops and eccentric BH binaries emitting in the frequency range (1-25) MHz.
The experimental GW group at Cardiff University is planning a set of two wideband table-top interferometers sensitive in the band (1-100) MHz (Vermeulen et al. 2021). These will be able to set new upper limits on a stochastic GW background in this frequency band.

Spherical resonant masses
The principle of a resonant mass detector is that its vibrational eigenmodes can get excited by a GW. These mechanical oscillations are transformed into electrical signals, using electromechanical transducers, and amplified by electrical amplifiers. These resonant detectors have a relatively small bandwidth, usually of less than 100 Hz. Thermal noise, Johnson-Nyquist noise, pump phase noise (if the transducer is parametric), back-action noise, and amplifier noise are the internal noises of this kind of detector. Therefore, the resonant mass antenna and transducers are made of high-quality factor materials in order to decrease thermal (mechanical) and Johnson-Nyquist (electrical) noises.
The idea of a spherical resonant mass antenna for GW detection has a long history and was first proposed by Forward (1971) followed by several decades of exploration and proposals (Wagoner and Paik 1977;Hamilton 1990;Johnson and Merkowitz 1993). In 1991, Aguiar proposed a large spherical antenna project in Brazil (Aguiar 2011). This detector, Mario Schenberg, in São Paulo, Brazil (Da Silva Costa and Aguiar 2014), was started to be built in 2000, around the same time as Mini-GRAIL, in Leiden, Netherlands. These two spherical detectors were active for about 15 years. At present, they are decommissioned, but Schenberg is planned to be reassembled at INPE, in São José dos Campos, about 100 km from its initial site at the University of São Paulo. 23 Such detectors have a bandwidth of 50-100 Hz with peak frequencies around 3 kHz for the quadrupole modes. To increase the frequency range, a xylophone configuration of several spheres has been proposed (Harry et al. 1996).
Spherical antennas provide more information, compared to the classical bar antennas, because of their quadrupole modes, while also being significantly more sensitive due to their favorable geometry of having a larger cross-section at identical mass. From the output of six transducers tuned to the quadrupole modes of the 23 These detectors had much smaller masses and diameters than originally proposed in the 1990s of up to 120 tons and 3 m, respectively, resonant around $ 700 Hz.
sphere, a single sphere can obtain complete information about the polarization and direction of the incoming wave.
In 2004, Mini-GRAIL reached a peak strain sensitivity of ffiffiffiffi ffi S n p ' 1:5 Â 10 À20 Hz À1=2 at a frequency of 2942.9 Hz at temperatures of 5 K. Over a bandwidth of 30 Hz, the strain sensitivity was about ffiffiffiffi ffi S n p ' 5 Â 10 À20 Hz À1=2 (Gottardi et al. 2007). Schenberg operating also at 5 K, reached strain sensitivities of ffiffiffiffi ffi S n p ' 1:1 Â 10 À19 Hz À1=2 for its quadrupolar modes ( $ 3:2 kHz) and ffiffiffiffi ffi S n p ' 1:2 Â 10 À20 Hz À1=2 for its monopolar mode ( $ 6:5 kHz), in 2015. Both antennas could reach sensitivities around ffiffiffiffi ffi S n p ' 10 À22 Hz À1=2 when operating at 15 mK. Schenberg, because it uses parametric transducers, can reach higher sensitivities if it implements squeezing of the signal. In this case, it would have similar sensitivities as the ultimate sensitivities of Advanced LIGO and Virgo around 3.2 kHz. 24 The conceptual difficulties in pushing this technology to higher frequencies are similar to the issues discussed for laser interferometers. Searching for GWs at higher frequencies requires smaller resonating spheres and consequently requires measuring smaller absolute displacements to achieve the same strain sensitivity. Note also that contrary to laser interferometers, resonant mass detectors have not yet reached the standard quantum limit yet. It thus seems unlikely that this technology can be pushed significantly beyond the kHz region.

Detection at frequencies beyond current detectors
In this section we will introduce several ideas and concepts for the detection of GWs at high frequencies beyond the capabilities of currently existing GW detectors.

Optically levitated sensors
Optically levitated dielectric sensors have been identified as a promising technique for resonant GW searches spanning a wide frequency band from a few $ kHz to $ 300 kHz (Arvanitaki and Geraci 2013;Aggarwal et al. 2020). A dielectric nanoparticle suspended appropriately at the anti-node of a laser standing wave within an optical cavity will experience a force when a passing GW causes a time-varying strain of the physical length of the cavity. The particle will be displaced from the location of the trapping light anti-node, resulting in a kick on the particle at the frequency of the GW space-time disturbance. The trapping frequency and mechanical resonance linewidth are widely tunable based on the laser intensity and laser cooling parameters chosen. For possible sources in this frequency band see e.g., Sects. 3.2.1, 3.2.2 and 3.2.3.
When detecting the resulting displacement of the particle at the trapping resonance frequency, the sensitivity is limited by Brownian thermal noise in the particle itself rather than the displacement detection of the particle. This results in improved sensitivity at higher frequency (unlike traditional interferometer style detectors which decrease sensitivity at high frequency due to laser shot noise) (Arvanitaki and Geraci 2013). The low-friction environment made possible by optical levitation in ultra-high vacuum enables extremely sensitive force detection (Ranjit et al. 2016), which becomes ultimately quantum-limited by photon-recoil heating from discrete scattering events of individual trap laser photons (Jain et al. 2016).
A 1-meter prototype Michelson-interferometer configuration detector called the 'Levitated Sensor Detector' is under construction at Northwestern University in the US, with a target sensitivity of better than ffiffiffiffi ffi S n p $ 10 À19 Hz À1=2 at 10 kHz and ffiffiffiffi ffi S n p $ 10 À21 Hz À1=2 at 100 kHz (Aggarwal et al. 2020). With Barker and coworkers at partner institution University College London, fiber-based approaches are being investigated to permit longer cavities without the need for expensive optics (Pontin et al. 2018). The ultimate strain sensitivity of a 10-meter room-temperature instrument is estimated to be better than approximately ffiffiffiffi ffi S n p $ 10 À20 Hz À1=2 at 10 kHz and ffiffiffiffi ffi S n p $ 10 À22 Hz À1=2 at 100 kHz. For a cryogenic 100-meter instrument this can be improved by more than an order of magnitude across much of the frequency band (Aggarwal et al. 2020). A detailed analysis of the search reach for GWs produced by axions via the BH superradiance process is provided in Aggarwal et al. (2020).

Inverse Gertsenshtein effect
The Gertsenshtein effect describes the conversion of photons to GWs in the presence of a magnetic field and was considered already decades ago as a source of GWs (Gertsenshtein 1962). While the coupling constant for this process is too small to be of interest for experiments in the near to medium future, the inverse Gertsenshtein effect, frequently referred to as magnetic conversion, can indeed be used to search for GWs (Boccaletti et al. 1970;Füzfa 2016Füzfa , 2017. While such dedicated instruments do not exist yet (apart from small prototypes), a first step in this direction has been done by using existing data from axion-search experiments. In these experiments, typically a strong static magnetic field of several Tesla is set up with field lines perpendicular to some interaction region, through which a beam line passes. In their nominal usage these experiments would search for axion-like particles, which can convert to photons in the presence of the magnetic field. These photons would be detected at the end of the beam line by electromagnetic detectors within the frequency band of interest (e.g. photodetectors for optical radiation). The very same experimental arrangements can also be used to search for GWs, by reinterpreting the acquired data, as has been pointed out and performed by the work in Ejlli et al. (2019). This work could set first upper limits for GWs at optical and Xray frequencies (i.e., around 500 THz and 10 6 THz respectively).
In the future this class of experiment is expected to continue with larger detectors of this sort being constructed. Given the motivation in the context of searches for high-frequency GWs, a dual usage of these detectors could be imagined with dedicated instruments and operational modes to search for GWs. For example, the planned IAXO detector (Armengaud et al. 2014) aims at searching for axions produced in the core of the sun. If it were fitted with different electromagnetic receivers, from radio to optical frequencies, GW searches could be facilitated in these bands. A particular advantage of searches with IAXO would be that the device can be pointed (within some limits) to different points in the sky. It could be of particular interest to point to patches in the sky where, for example, a binary BH merger is predicted to happen. The latter is a real prospect once the LISA space interferometer will be operational, which can detect inspiralling BHs long before their merger.
Note that, in principle, the inverse Gertsenshtein effect might be exploited at all frequencies, and that one further advantage of this concept is the tunability of the solid angle, which changes according to the direction of the magnetic field, see Sect. 4.4.5. In particular, the inverse Gertsenshtein effect has substantial room for development especially at GHz frequencies where many of the early Universe signals converge. Using magnets developed for particle accelerators a conversion path length of 100 metres with uniform magnetic field of 5.6 T is quite realistic, and being implemented for the ALPS experiment searching for axion-like particles (Bähre et al. 2013). Electromagnetic detectors having a thermal noise equivalent to 0.1 Kelvin would lead to a sensitivity around h c;n ' 10 À26 . This sensitivity could be further enhanced by the inclusion of a Fabry Perot cavity in the conversion volume and a factor of 100 improvement might be possible. Unfortunately, this improved sensitivity would be gained at the expense of the wide bandwidth of the technique and this would limit the applicability to stochastic signals. Ringwald et al. (2021) estimates the sensitivity that can be reached in the GHz region by using Single Photon Detectors (SPD) and Heterodyne radio receivers (HET). The corresponding limits are reported in Figs. 1 and 2.
The conversion of GWs to photons by the inverse Gertsenshtein effect cannot only be exploited in a laboratory setting but also by considering astrophysical or even cosmological 'detectors', see e.g., Pshirkov and Baskaran (2009);Chen (1995); Dolgov and Ejlli (2012); Cillis and Harari (1996); Domcke and Garcia-Cely (2021). In this case, the magnetic field is weaker and the background is much harder to control, but cosmic magnetic fields can extend coherently over kpc or even Mpc, implying an enormous 'detector volume'. The frequency range of MHz to GHz coincides with the Rayleigh-Jeans tail of the cosmic microwave background, a target of existing and upcoming radio telescopes. For example, the data of ARCADE 2 (Fixsen et al. 2011) and EDGES (Bowman et al. 2018) can be recast to respectively constrain h c;sto \10 À24 ð10 À14 Þ in the range 3 GHz .f . 30 GHz and h c;sto ðf % 78 GHzÞ\10 À12 ð10 À21 Þ for the strongest (weakest) cosmic magnetic fields in accordance with current astrophysical data (Domcke and Garcia-Cely 2021).

GW to electromagnetic-wave conversion in a static electric field
Lupanov (1966) considered the inverse Gertsenshtein effect but using a static electric field rather than a static magnetic field. The physics is essentially the same in the two cases but since the intensity of electric fields in laboratory settings is limited by the tendency to pull electrons from nearby conductors (or dielectrics) and thereby cause local short circuits, the available energy densities in electric fields are about one millionth of that created by magnetic fields in the several Tesla range. Hence the use of electric fields seems not to offer any advantages. Cruise (1983) showed that a GW could induce a rotation of the plane of polarisation in electromagnetic waves in certain geometries, some of which might be relevant astronomically. In 2000 the idea of resonant polarization rotation was extended (Cruise 2000) to a situation in which the electromagnetic wave was a circulating wave in a microwave waveguide ring. The original effect was amplified by the (potentially significant) quality factor of the waveguide ring. A proof of concept apparatus was constructed by Ingley (2005, 2006). Such a device would be narrowband with a sensitivity ffiffiffiffi ffi S n p $ 10 À14 Hz À1=2 at frequencies of 100 MHz. It is difficult to see the sensitivity of this scheme for GW detection increasing very far beyond the published value.

Heterodyne enhancement of magnetic conversion
Li and Yang (2004); Li et al. (2006); Baker et al. (2008); Li et al. (2009) have suggested enhancing the conversion efficiency of magnetic conversion detectors such as those discussed in Sect. 4.2.2, see also Ringwald et al. (2021) for a recent discussion. This proposal has been specifically aimed at the detection of cosmological (relic) signals of a stochastic nature and with dimensionless amplitudes in the range h c;sto $ ð10 À30 À 10 À26 Þ at 5 GHz, about the highest signal consistent with the BBN limit. The conversion from GW to electromagnetic wave is enhanced by seeding the conversion volume with a locally generated electromagnetic wave at the same frequency as that being searched for. In conditions in which the gaussian local oscillator beam is parallel to the incoming signal and at right angles to the static magnetic field, an additional beam of electromagnetic waves is generated by the conversion process, travelling at right angles to the incoming beam and the locally generated beam. The technical challenge is then to distinguish this perpendicular beam of, say, 800 photons/second from the locally generated beam at the same frequency and carrying 10 24 photons/second at frequencies of several GHz (Woods 2012). This demands a geometric purity in the Gaussian beam of better than 10 À21 , far beyond the current state of the art. The authors have proposed this interesting idea over many years but the lack of laboratory results on the performance of the necessary subsystems leaves the feasibility of this concept an open question.

Bulk acoustic-wave devices
Bulk acoustic-wave devices are one of the pillars of frequency control and frequency metrology . In its simplest form, a piece of piezoelectric material is sandwiched between two electrodes, converting the acoustic waves inside the material into electrical signals. With its relatively compact size and robustness, this technology gives one of the best levels of frequency stability near one second of integration time. More recently, it was demonstrated that quartz bulk acoustic-wave devices exhibit extremely high-quality factors (up to 8 Â 10 9 ) at cryogenic temperatures for various overtones of the longitudinal mode covering the frequency range (5-700) MHz Goryachev et al. 2013). For this reason it was proposed to use the technology for various tests of fundamental physics ) such as Lorentz invariance tests (Lo et al. 2016), quantum gravity research (Bushev et al. 2019) and search for high-frequency GWs . For the latter purpose, a bulk acoustic-wave device represents a resonant mass detector whose vibration could be read through the piezoelectric effect and a Superconducting Quantum Interference Devices (SQUIDs). The approach has the following advantages: highest quality factor (high-sensitivity), internal (piezoelectric) coupling to SQUIDs , allows parametric detection methods, large number of sensitive modes ([100) in a single device, modes scattered over wide frequency range (1-700) MHz, well-established and relatively inexpensive technology (mass production), high-precision (insensitive to external influences such as seismic vibration and temperature fluctuations). On the other hand, it is shown that at low temperatures identical devices demonstrate significant dispersion in mode frequencies, thus, showing low accuracy. The level of sensitivity of bulk acoustic-wave detectors is estimated at the level of ffiffiffiffi ffi S n p ' 2 Â 10 À22 Hz À1=2 subject to the mode geometry . With additional investment into research and development, this level can be improved and the frequency range extended down to hundreds of kHz range.
A search for high-frequency GWs with a single bulk acoustic-wave devices and two modes at 4 K has been running in the University of Western Australia since November 2018.

Superconducting rings
The quantum properties of vortices in superfluids may interact with spin components of the GWs. In addition, an extension of the electromagnetic impedance at a boundary in the case of GWs in a superconducting fluid suggests that the impedance mismatch well-known in classical bar detector theory is much reduced for a GW arriving at a boundary in a superconductor, essentially creating a very efficient mirror that could be used as a building block for an interferometer or Sagnac ring. Anandan and Chiao (1982); Chiao (2002) proposed a new detector format which utilises these putative principles, reaching a sensitivity of h 0;n;mono $ 10 À31 . Resonant operation will restrict the bandwidth in the GHz range. A good review of the issues surrounding the interaction of mesoscopic quantum systems with gravity was prepared on an European Space Agency contract by Kiefer and Weber (2005). This review casts doubt on some of the assumptions made by Anandan and Chiao. Caves (1979) published a theoretical study of a microwave cavity with a high mechanical quality factor. Mechanical deformation of the cavity by a GW coupled two of the cavity's resonant microwave modes and transferred the electromagnetic excitation to a previously unexcited cavity mode. Reece et al. (1984) built a similar system with a higher resonant frequency of 1 MHz and one operating at 10 GHz (Reece et al. 1982), while Pegoraro et al. (1978) designed a system with a sharp resonance at about 1 GHz. These schemes certainly offer some sensitivity in the frequency range above 1 GHz but that is limited to around h $ 10 À21 by the thermal noise in the microwave sensors. Even at cryogenic temperatures the sensitivity will be many orders of magnitude away from the level required for detecting cosmological sources. The bandwidth of this scheme is nominally limited by the bandwidth of the cavity, up to effects like detection or electronics gains and noises.

Graviton-magnon resonance
As pointed out in Ito et al. (2020), a GW passing through a ferromagnetic insulator can resonantly excite magnons (collective excitations of electron spins), similar to the excitations of phonons in resonant bar detectors. The readout is achieved by placing the magnetic sample inside a microwave cavity, coupling the magnon to a photon mode. This idea builds on the technique of ferromagnetic haloscopes proposed to for axion searches (Crescini et al. 2018;Flower et al. 2019). The sensitivity of such detector reaches strain of ffiffiffiffi ffi S n p ' 7:6 Â 10 À22 Hz À1=2 at 14 GHz and ffiffiffiffi ffi S n p ' 1:2 Â 10 À20 Hz À1=2 at 8.2 GHz. The sensitivity of this approach can be greatly improved by incorporating single frequency counters that are already available. A few orders of magnitude in sensitivity have been shown for axion detection (Lamoreaux et al. 2013).

Summary of detector sensitivities
In Table 1 we summarize the existing and proposed technologies for high-frequency GW detection, reporting the corresponding sensitivities. For all experiments that quote their sensitivity in terms of a power spectral density noise S n ðf Þ, a conversion to dimensionless strain h c;n;sto has been performed using Eq. (16b) with 1 year as the observation time and the specified detector bandwidth. For the other detectors, we specify the dimensionless strain variable on a case by case basis: we use h c;n;sto (see Eq. (16b)) for detectors that look for a stochastic signal, while we use h 0;n;mono (see Eq. (18b)) for detectors that look for a monochromatic GW. In the case of microwave cavities (see Sect. 4.2.8) we denote the sensitivity by h as it can refer either to bursts or to long duration signals, and we refer the reader to the original papers for the details, while we report the best sensitivity estimates. We also specify whether each experiment has already been built, is under construction, is only devised or only the physical mechanism has been identified (theory). The sensitivity values labeled by an asterisk refer to the planned sensitivities that will be achieved ð10 À22 À 10 À20 Þ Hz À 1 2 by the proposed future improvements of the currently built setups. Note that a square bracket in the frequency column refers to the bandwidth of the detector, while a round bracket refers to the range of frequencies 25 that can be covered by the detector itself. We report the bandwidths used in Table 1 to obtain dimensionless strain from PSD: Mini-GRAIL had a bandwidth Df ' 30 Hz (Gottardi et al. 2007); the Schenberg antenna had a bandwidth Df ' 50 Hz (Aguiar 2011); the 0.75-m interferometer has bandwidth Df ' 2 kHz (Akutsu et al. 2008) 26 ; the optically levitated sensors have a bandwidth Df ¼ f =10 (Arvanitaki and Geraci 2013;Aggarwal et al. 2020); Df ' ð10 À 50Þ kHz for Cruise's and Ingley's detector Ingley 2005, 2006); for enhanced magnetic conversion the bandwidth is Df ' 1 Hz (Li et al. 2009); bulk acoustic wave resonators have bandwidth where f is in the ranges reported in Table 1  ); the bandwidth for Pegoraro's detector is Df ' 1 Hz. For references that specify their detector sensitivity in dimensionless strain without specifying the exact form of the dimensionless strain, the sensitivities are labeled as just h.

Cross-correlation detectors
For a coalescing binary the information about the waveform of GWs is available.
The best approach to the detection problem is to use this information by projecting the observed data over the set of expected signals. Usually this set can be parameterized by a small number of parameters, varying in some allowed range. A 'scalar product' between data and the set of 'templates' is evaluated, using its maximum as a detection statistics. It should be noted that the set of expected signals does not have a linear space structure, in the sense that the linear combination of two possible signals is not generally speaking a possible signal. For this reason the computational cost of a search over a 'template bank' grows very fast with the number of free parameters. When the number of parameters is large, or when a parameterization of the waveform of the expected signal is not possible at all, other detection methods must be used. In the present context this is the case for several cosmological processes, which are expected to produce a GW signal that can be described as an overlap of a very large number of contributions. There is not a waveform here, the expected signal is a stochastic process and the best approach for the detection is the cross correlation one.

Relic gravitational radiation
The detection of the CMB in 1965 by Penzias and Wilson gave the first experimental insight into the properties of 'relic' radiation, the remnants currently observable of the Big Bang. The CMB is a stationary, stochastic radiation field, basically isotropic down to levels as low as 10 À5 with a Gaussian distribution.
It is natural to suppose that a useful starting point in planning GW observations of a relic radiation would be to assume that cosmologically-sourced GWs can be modeled by a stochastic field h ij ðx; tÞ with relatively simple properties. In the spacetime volume of a given experiment, we can describe our background as a superposition of plane waves Here P labels the polarization degrees of freedom of the field, whose number can depend on the considered theory of gravitation. The 'parameters' of the signal are the amplitudes h P ðkÞ introduced in Sect. 3, one for each mode of the field. When these must be considered stochastic variables, the relic gravitational radiation field is completely described by their (joint) probability distribution, and is called GW stochastic background. For a Gaussian stochastic background this joint probability distribution is Gaussian, and is completely determined by the second order expectation value Further assumptions lead to the simplification of this general expression. For example, stationarity in a given reference frame requires the expectation value hh ij ðx; tÞh ij ðy; t 0 Þi to be a function of t À t 0 only, and as a consequence the correlation between h A ðpÞ and h B ðqÞ can be non-zero only when x p ¼ x q . Stationarity in every reference frame implies homogeneity, and only correlations with p ¼ q are allowed. Each stochastic model can have its peculiar signature (see also Sect. 3): in the simplest stationary, isotropic and Gaussian model a parameterization of the model can be given in term of an array of functions which are connected to contributions which is a generalization of Eq. (1) and allows for a non trivial polarization structure.
If the radiation is stationary, a temporal signature which could be exploited by a single-instrument detection procedure is not available. Moreover, isotropy and homogeneity do not allow for a signal modulation which could be obtained in principle by changing the orientation or the position of the detector. 27 It could still be possible to detect the stochastic background as an 'excess noise' in the apparatus. However in order to do that the amplitude of the signal must be large enough to make it evident given a theoretical estimate of the noise budget, which is always uncertain. This means that the strategy for detection will necessarily be different from the strategy for discrete source detection.
The most obvious approach is the use of spatial correlations. If a detector is to be developed for the detection of GW relic radiation then a decision to operate it as a correlation detector will have a far reaching influence on many aspects of its design.
An excellent review of GW relic radiation and appropriate methods of detection, in the context of HF band, has been published by Allen (1997) and important properties of correlation detectors have been explored by Michelson (1987). The basic principle is to compare the signal from two detectors. This is comparing a random signal with another stationary, stochastic, isotropic, gaussian signal from the same source. Similar to template matching as a means of detecting discrete sources, in this case the template itself is random and that affects the statistical gain from performing a cross correlation between two detectors.

Properties of correlation detectors
At the simplest level the outputs s i ðtÞ of two detectors can be written as s 1 ðtÞ ¼ h 1 ðtÞ þ n 1 ðtÞ ; where h 1;2 are the stochastic signals and n 1;2 the noises. We suppose both the signals and the noises to have typical amplitudes h i , n i and typical bandwidth Df . If we evaluate the simple correlation over a time T ) Df À1 between the two outputŝ we find an estimator which is distributed as a Gaussian variable for large enough T. Its mean is given by where the temporal cross-correlation C h 1 h 2 is defined as and we supposed that the two noises are uncorrelated. The variance ofŶ is given by 27 Note however that the motion of the detector with respect to the cosmic rest frame of a cosmological stochastic background breaks isotropy. This yields a (weak) signal modulation and can be exploited to extract polarization information from the stochastic background of GWs (Seto 2006(Seto , 2007Domcke et al. 2020).
Z tþT t dt 00 C n 1 n 1 ðt 00 À t 0 ÞC n 2 n 2 ðt 00 À t 0 Þ ' ' Z 1 À1 dsC n 1 n 1 ðsÞC n 2 n 2 ðsÞ ' ' TC n 1 n 1 ð0ÞC n 2 n 2 ð0Þ 1 Df ; and we see that the signal-to-noise ratio hŶi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hŶ 2 i À hŶi 2 q / C h 1 h 2 ð0Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C n 1 n 1 ð0ÞC n 2 n 2 ð0Þ p ffiffiffiffiffiffiffiffi ffi increases as the square root of the measurement time. The minimum detectable energy density X GW is and decreases with the square root of T. Of course the minimum detectable signal amplitude is proportional to ffiffiffiffiffiffiffiffiffi ffi X min GW q , so it decreases much slowly.

The overlap function
In the previous, very simple analysis we neglected all the effects that can be adsorbed in a proportionality factor. For example when we evaluated the signal, we wrote implicitly assuming that K ¼ Oð1Þ. As a matter of fact the correlation between the signals coupled to two different detectors can be reduced by several effects. The quantity of interest can always been obtained by the correlation between the gravitational strain at two different points after a contraction on appropriate tensors which describe the detectors. Without entering in the details of the calculation, we see that the value of the previous correlation is influenced by two effects. First of all, the detectors will not be in the same position. In this case the phase factor in the integral will oscillate, and the correlation will be reduced. It is clear that this reduction will start when the separation between the two detectors will be larger than the wavelength of the field at the frequency of interest. In the very highfrequency regime, this will happen at very small separations, d '1=ð2pf Þ.
A further reduction of the correlation can be generated by two detectors which are coupled differently to the modes, for example because they are oriented differently.
The reduction of correlation is quantified by the so-called 'overlap function' cðf Þ, a frequency dependent factor with modulus always less than 1, which is simply the coherence between the two signal of interest. In a quantitative analysis cðf Þ must be included so as to diminish the signal correlation and increase the minimum detectable amplitude by a factor 1 cðf Þ . Michelson (1987) has worked through the derivation of cðf Þ and Allen (1997) has outlined the process of optimising the detection efficiency by optimal filtering in the time domain for two detectors with arbitrary separation and orientation.

Exploiting c(f Þ at very high frequencies
One possibility that arises at high frequencies due to the small wavelength of the radiation being studied is to construct laboratory scale detectors which in principle could be moved relative to one another, hence changing the value of the overlap function c. In this way the correlation of the signal can be modulated, and a detection of this modulation pattern could provide credible evidence of a relic radiation detection. It is interesting to note that the overlap reduction function depends on the polarization structure of the field. This means that by looking at the modulation pattern it is in principle possible to disentangle non standard polarizations coming from extended gravitational theories.

Signal switching
The prospect of a new detection technology for high-frequency GW detection opens a window for a slightly novel form of correlation detector. Instead of correlating the output signal with that of a similar detector, it may be possible to turn off the sensitivity of a single detector to GWs without affecting its other performance properties. If the temporal pattern of switching on and off can be seen in the signal output at a statistically significant level then a credible claim for a detection might be made.
An opportunity of this kind presents itself for correlation detection in the case of magnetic conversion detectors (see below), in which case signal switching could also be achieved by modulating the amplitude of the field and its direction. The devices are arranged to enable the conversion of the GW to an electromagnetic one within an electromagnetic cavity crossed by a strong magnetic field. The electromagnetic signal produced is proportional to the length of the cavity in GW wavelengths and the cavity dimensions transverse to its length are critical in determining the efficiency of the process. This is because the generated electromagnetic wave finds itself travelling in an electromagnetic waveguide, with a phase velocity greater than that of the GW. To make the conversion efficient these two phase velocities must be sufficiently close that, in the whole length of the detector, only a very small phase difference develops between the two waves. If this is not the case then the conversion process is essentially rendered inoperative. Current thoughts on slowing the electromagnetic wave in the cavity include the introduction of a gas into the cavity and adjusting the refractive index at the electromagnetic frequency by means of the pressure to equalise the phase velocities. Because the conversion is so inefficient if this is not done, controlling the gas pressure can essentially switch the sensitivity to GWs off and on with very minimal disturbance to its other operating parameters and hence to the output signal. The detector output can then be correlated with the detector operating pattern. Statistically, this is a more effective correlation process than described above between two similar detectors because in this case the signal is compared with an a priori determined template and not a random template. The minimum detectable signal in this case is / T À 1 2 allowing a faster gain in sensitivity with time.

Issues related to data acquisition and long term storage
To detect correlated periodic events in nearby detectors at frequencies around 1 MHz at an signal-to-noise ratio of 8, systematic errors related to timing should be of the order of 20 ns. This necessitates the need for a careful understanding of the various factors that contribute to errors arising from timing calibration. From the hardware side, low noise amplifiers, anti-aliasing electronics, etc. add delay to the data acquisition system. Quantization errors arising from the analog-to-digital converters add further delay and to minimize the effect, their sampling period would have to be made lower than the desired timing resolution. At such sampling rates, making use of super-conducting oversampling ADCs which achieve high dynamic ranges over narrow frequency ranges by pushing quantization noise outside the band-of-interest could turn out to be viable option. Both stochastic and continuous wave analysis rely on cross-correlation based search strategies and make use of year long data-sets for their analysis. Moving from current audio-band frequencies to megahertz regime can easily scale up the storage requirements to few petabytes of data. However it must be noted that the increase of storage is not proportional to the frequency of interest, but on the bandwidth, as the typical observation frequency can always scaled down with an appropriate heterodyne technique.
This calls for real-time analysis as proposed in SKA where the raw data is discarded after the low latency retrieval of relevant information. Advantages of using cloud and web 2.0 technologies for the collaborative development of the various data analysis pipelines used for signal detection and parameter estimation also needs to be investigated. Making use of data folding techniques based on inherent symmetries, such as Earth's siderial day rotation has been shown to decrease data volume in stochastic background searches at audio-band frequencies (Ain et al. 2015). Stacking years long data into a single day while preserving all the statistical properties would enable us to carry out the final analysis on personal computing devices.

Discussion and conclusions
The search for high-frequency gravitational waves is a promising and challenging search for new physics. It provides a unique opportunity to test many theories beyond the Standard Model that could not be tested otherwise.
Various models proposed to address open questions in particle physics and cosmology predict a gravitational-wave signal in the frequency range f ' ð10 3 À 10 10 Þ Hz, which could be a coherent signal -e.g., from mergers of compact objects or from axion superradiance around black holes -or stochastice.g., originating from certain models of cosmic inflation, from a phase transition in the very early Universe or from oscillons, evaporating primordial black holes, etc. See Figs. 2 and 1 as well as Tables 2 and 3 for an overview. Many of these models can lead to relatively large signals, corresponding to an Oð1Þ fraction of energy density in the early Universe converted to gravitational waves. This energy is redshifted in the expanding Universe, rendering even these strong signals challenging to detect today. Moreover, in many cases the amplitude of the signal depends sensitively on the model parameters and may be significantly lower in large parts of the model parameter space.
The high-frequency band comes with particular challenges and opportunities. High-frequency gravitational waves carry a high-energy density, implying that bounds provided by cosmology on the fraction of energy contained in gravitational waves translate to stringent bounds on the characteristic gravitational-wave strain. This poses a severe challenge for detection, since the magnitude of any observable effects is typically governed by the strain and not by the energy density. This renders the detection of cosmological sources of high-frequency gravitational waves much more challenging than comparable searches at lower frequencies. On the other hand, the lack of known astrophysical gravitational-wave sources in this frequency range poses a unique opportunity for foreground free searches of new physics.
At the moment, there is no general consensus on the most promising detection strategy in this frequency band, though many proposals have been put forward in the past decades. The proposals that we are aware of are summarized in Table 1, together with their frequency and sensitivity range. We emphasize that the same sensitivity (in terms of characteristic strain) at a higher frequency typically implies a reduced sensitivity to the viable parameter space of a given cosmological source, as discussed above. In this sense, detectors based on magnetic conversion or on the deformation of microwave cavities seem to be particularly promising avenues, though a more careful study of noise levels and of the margin on improvement with foreseeable technology development is needed in many cases.
We hope that this document will stimulate the necessary discussion and we strongly encourage feedback regarding further proposals or critical assessments which we may have missed. We have the ambition to consider this document as a first step towards a coherent international collaboration to seriously consider the search for high-frequency gravitational waves.
None of the proposals listed in this report currently reach the sensitivity needed to probe the new physics outlined above. At best, achievable proposals are still at least six orders of magnitude beyond the required sensitivity. However, we recall that, one hundred years ago, the technological gap in both the LIGO and LISA frequency ranges was of about 16 orders of magnitude (Chen et al. 2017). Also, less than 50 years ago, Misner, Thorne and Wheeler, declared that 'such detectors have so low sensitivity that they are of little experimental interest' (Misner et al. 1973), referring to laser interferometers. The first laser interferometer gravitational-wave detector, built at Hughes Research Laboratories in the 1970s (Forward 1978) had a sensitivity which was eight orders of magnitudes below the design sensitivity of the currently operating LIGO/Virgo/KAGRA detectors. There are currently clear development paths leading to detectors operating at sensitivities of h c;n ' 10 À26 using, e.g., magnetic conversion (see Sect. 4.2.2) and several new ideas such as magnon devices and superconducting systems which have received little detailed design development so far.
We therefore take the past history of laser interferometry as an encouraging lesson for the development of gravitational-wave detectors in the high-frequency band. The challenges are strong but the opportunities are unique. As a rule of thumb, probing very early epochs of our cosmic history and consequently particle physics at very high-energy scales requires searching for gravitational waves at high frequencies with correspondingly small experimental devices. As can be seen from Figs. 1 and 2, the Ultra High-Frequency band, ranging from the MHz to the GHz is an exciting window to explore fundamental physics up to the grand unification or string theory scales of order ð10 16 À 10 17 Þ GeV, thereby probing our cosmological history right up to the Big Bang. A detection of a gravitational-wave signal in this frequency range would be smoking gun signal for new physics, since no known astrophysical processes can generate sizable gravitational-wave signals at these frequencies. It would be remarkable if the experimental test of fundamental physics at the highest energies and earliest times in the history of the Universe could eventually be achieved not with huge particle accelerators nor satellite interferometry, but with small size table-top experiments.
This white paper set up the stage for the launch of the Ultra-High-Frequency Gravitational Wave (UHF-GW) initiative 28 , that supports the creation of a network of researchers for the development of gravitational-wave science in the highfrequency range. One of the goals of the initiative is to stimulate the technological development that is needed to build successful gravitational-wave detectors at high frequency. Table 2 Summary of coherent sources. The characteristic strain is given in Eqs. (5) and (6)