New sensitivity curves for gravitational-wave signals from cosmological phase transitions

Gravitational waves (GWs) from strong first-order phase transitions (SFOPTs) in the early Universe are a prime target for upcoming GW experiments. In this paper, I construct novel peak-integrated sensitivity curves (PISCs) for these experiments, which faithfully represent their projected sensitivities to the GW signal from a cosmological SFOPT by explicitly taking into account the expected shape of the signal. Designed to be a handy tool for phenomenologists and model builders, PISCs allow for a quick and systematic comparison of theoretical predictions with experimental sensitivities, as I illustrate by a large range of examples. PISCs also offer several advantages over the conventional power-law-integrated sensitivity curves (PLISCs); in particular, they directly encode information on the expected signal-to-noise ratio for the GW signal from a SFOPT. I provide semianalytical fit functions for the exact numerical PISCs of LISA, DECIGO, and BBO. In an appendix, I moreover present a detailed review of the strain noise power spectra of a large number of GW experiments. The numerical results for all PISCs, PLISCs, and strain noise power spectra presented in this paper can be downloaded from the Zenodo online repository [1]. In a companion paper [2], the concept of PISCs is used to perform an in-depth study of the GW signal from the cosmological phase transition in the real-scalar-singlet extension of the standard model. The PISCs presented in this paper will need to be updated whenever new theoretical results on the expected shape of the signal become available. The PISC approach is therefore suited to be used as a bookkeeping tool to keep track of the theoretical progress in the field.


Introduction
Since the celebrated first direct detection of gravitational waves (GWs) in September 2015 [3], the Advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) [4,5] and the Advanced Virgo (aVirgo) experiment [6] have observed a multitude of GW events [7,8]. One of the signals recorded during the first two aLIGO / aVirgo observing runs [9] originated from the coalescence of a binary neutron star [10], while all other signals were due to mergers of binary black holes [3,7,[11][12][13][14]. First results from the third observing run are reported on in ref. [8]. After the breakthrough discovery of these transient astrophysical sources, one of the next key objectives in GW astronomy is going to be the detection of a stochastic GW background (SGWB) from cosmological sources [15,16]. A wide range of violent phenomena in the early Universe can give rise to a primordial SGWB [17,18], among which cosmological phase transitions [19,20] are a preeminent example. Strong first-order phase transitions (SFOPTs) occur in numerous extensions of the standard model (SM) of particle physics and can potentially lead to a strong GW signal [21]. Many scenarios especially predict a GW signal that is peaked at frequencies in the milli-Hertz range. This makes GWs from a SFOPT a prime target for future satellite-borne JHEP01(2021)097 interferometers in space, such as the Laser Interferometer Space Antenna (LISA) [22,23], which is scheduled for launch in the 2030's. Consequently, many authors have recently studied the prospects of probing the dynamics of a SFOPT in scenarios beyond the standard model (BSM) with LISA [24,25] and other upcoming GW experiments (see, e.g., the incomplete list of BSM models in refs. ).
The GW signal from a SFOPT is conventionally expressed in terms of a GW energy density spectrum Ω signal (f ) as a function of GW frequency f , while the instantaneous sensitivity of a GW experiment is quantified in terms of a noise spectrum Ω noise (f ) [see eq. (A. 24) for the precise definition of these two quantities]. With these two spectra at one's disposal, one is able to assess the chances that the predicted signal is going to be experimentally detected. In practice, this is typically done by adopting one or both of the following two strategies: • Strategy #1: compute the associated signal-to-noise ratio (SNR) by integrating over the experiment's total observing time t obs and accessible frequency range where n det distinguishes between experiments that aim at detecting the SGWB by means of an auto-correlation (n det = 1) or a cross-correlation (n det = 2) measurement.
Then, if turns out to be larger than some threshold value, > thr , one concludes that the GW experiment under consideration will be able to detect the predicted GW signal.
• Strategy #2: construct the power-law-integrated sensitivity curve (PLISC) Ω PLIS (f ) [85] based on the noise spectrum Ω noise (f ) (and some thr ) and compare it to the signal spectrum Ω signal (f ). Then, if the signal and the PLISC intersect, such that Ω signal (f ) > Ω PLIS (f ) for some f , one typically also concludes that the experiment will be able to detect the signal.
Both strategies have several advantages and disadvantages. The SNR approach, e.g., has a clearly defined statistical interpretation and is applicable for an arbitrarily shaped signal Ω signal . On the other hand, the SNR no longer contains any spectral information because of the integral over the frequency range [f min , f max ] in eq. (1.1). This is somewhat unfortunate. Ideally, one would like to indicate whether a certain signal is going to be experimentally detected or not directly in a plot of the signal spectrum Ω signal . In fact, this has been one of the main motivations behind the idea of constructing PLISCs. However, as is no longer a function of f , there is no canonical way of including information on in a plot of Ω signal . As a consequence, many authors resort to graphical representations of as a function of the underlying model parameters {p i } that determine the shape of the signal, = ({p i }). This typically results in plots of (a subspace of) the model parameter space that contain information on in the form of contour lines or a color code. While such information may be very useful for a number of reasons, it can easily JHEP01(2021)097 happen that one ends up with plots that show only as a function of quantities that are not directly accessible by experiments. This may be regarded as a disadvantage of the SNR approach, as one may prefer to indicate an experiment's sensitivity reach in terms of physical observables rather than in terms of auxiliary model parameters. Often times, one also looses information because one is forced to restrict oneself to lower-dimensional hypersurfaces (typically two-dimensional slices) in the higher-dimensional model parameter space. Finally, one may argue that another disadvantage consists in the fact that the SNR approach requires a larger computational effort. While the PLISC approach is essentially a graphical one, computing the SNR always involves the extra step of carrying out the frequency integration in eq. (1.1). This is part of the reason why many authors in the literature actually content themselves with a graphical analysis in terms of PLISCs and refrain from performing a proper SNR analysis.
The PLISC approach, by contrast, manages to convey a useful and graphical impression of an experiment's sensitivity directly in terms of plots of the GW spectrum. However, it is important to note that, by construction, PLISCs do not encode information on the expected SNR as soon as the spectrum deviates from a pure power law. Therefore, in realistic situations where the signal is expected to have a richer structure than just a simple power law, PLISCs should rather be regarded as a qualitative visualization than a quantitative statistical tool.
The aim of the present paper is to remedy the shortcomings of strategies #1 and #2 for the particular case of GWs from a SFOPT in the early Universe. The main observation in this case is that the shape of the signal is model-independent to first approximation, such that it is actually not necessary to perform the frequency integral in eq. (1.1) over and over again. Instead, it is possible to compute it once and for all, whereupon the numerical result may be used for the signal in any BSM model that one is interested in. In realistic scenarios, the GW signal from a SFOPT receives contributions from several different physical sources (see section 2). However, for illustration, let us suppose for a moment that there is just one physical source of GWs. In this case, the signal spectrum can be schematically written as Ω signal (f ) = Ω peak signal ({p i }) S (f, f peak ) , (1.2) where Ω peak signal ({p i }) denotes the peak strength of the signal at the peak frequency f peak , and where the shape function S describes the frequency dependence of the signal in the vicinity of the peak frequency. Here, the shape function S is assumed to be model-independent, whereas the peak amplitude Ω peak signal captures the detailed dependence on the underlying model parameters {p i }. For a signal of this form, it is then possible to rewrite eq. (1.1) as follows, where Ω PIS (f peak ) denotes what we shall refer to as the peak-integrated sensitivity (PIS),

JHEP01(2021)097
For any given experiment, the integral in eq. (1.4) only needs to be computed once. As soon as this has been done, one can construct a peak-integrated sensitivity curve (PISC) by plotting Ω PIS as a function of the frequency f peak . In this plot, the original (one-dimensional) spectrum Ω signal then reduces to a single (zero-dimensional) point f peak , Ω peak signal . Thus far, the standard procedure for phenomenologists and model builders interested in studying the GW signal from a SFOPT typically involved three steps. First, one had to consult the latest theoretical literature on the functional forms of the peak amplitude Ω peak signal and the shape function S. Then, in a second step, one had to consult the experimental literature on the shape of the noise spectrum Ω noise for a given experiment. And finally, one had to tie both ends together and perform the last step, namely, the frequency integration in eq. (1.1), by hand. The PISC approach now closes the gap between the input from the theory side and the input from the experimental side by rendering the final frequency integration obsolete. In fact, we argue that it unifies the advantages of the two standard strategies #1 and #2 that we outlined above, while at the same time avoiding their disadvantages: 1. As soon as the PISC has been constructed for a particular experiment, it is no longer necessary to perform a frequency integration on a parameter-point-by-parameterpoint basis. Instead, one can simply work out a fit function for the exact numerical PISC (see section 3.7), which enables one to write down a quasianalytical expression for the SNR in eq. (1.3). This property turns PISCs into a handy and ready-to-use tool for the phenomenological exploration of BSM models that predict GWs from a SFOPT.
2. The PISC approach works for an arbitrary signal shape S. Unlike the PLISC approach, it is not restricted to signals that are (reasonably well) described by a pure power law.
3. At the same time, it results in a useful visual representation of an experiment's sensitivity in terms of physical observables (as opposed to auxiliary model parameters), namely, the peak frequency f peak and the peak amplitude Ω peak signal of the GW spectrum.
4. According to the PISC approach, signal spectra are projected onto individual points in the f peak -Ω peak signal plane. This can be used to generate scatter plots that allow one to identify what one may call a model's signal region. Plotting this signal region in combination with a PISC then indicates to what extent the model is going to be probed experimentally. This aspect is elaborated on in more detail in the companion paper [2]. 5. If new theoretical results should require a revision of the spectral shape S, it suffices to update the experimental PISCs, while the model-specific signal regions remain unaffected. Similarly, new insights into the dependence of f peak and Ω peak signal on the underlying SFOPT parameters only cause the signal regions to shift, but leave the experimental PISCs unchanged. This facilitates the update of sensitivity plots com-

JHEP01(2021)097
paring the predictions of many models at the same time (see also section 3.6 and ref. [86]). 6. The fact that PISC plots correspond to projections onto the f peak -Ω peak signal plane distinguishes them from the usual plots of the SNR on two-dimensional slices through the higher-dimensional parameter space that one often encounters. In the PISC approach, it is not necessary to keep a subset of model parameters fixed at specific values.
7. Most importantly, the PISC approach retains the full information on the SNR. For a data point f peak , Ω peak signal , the expected SNR simply corresponds to the vertical separation between this point and the PISC of interest. This facilitates the implementation and comparison of different SNR thresholds thr . All the relevant information is encoded on the y-axis of our plots; no additional color code or contour lines are needed. In addition, our plots allow for an easy comparison of the PISCs and hence expected SNRs for different experiments. Such a comparison would not be feasible in plots using a color code for the SNR and significantly more complicated in SNR contour plots.
Before we turn to a more detailed presentation of the PISC approach, it is worth pointing out some similarities to alternative treatments in the literature. Refs. [56,66], e.g., also present scatter plots of possible peak frequencies and peak amplitudes in specific BSM models. These predictions are, however, combined with the usual sensitivity curves, such that the final plots do not contain any information on the expected SNR. Moreover, ref. [87] studies the sensitivity of future GW satellite experiments based on a Fisher matrix analysis that quantifies the precision with which one will be able to reconstruct observables such as f peak and Ω peak signal from real data. Studies of this kind will be an important part of the data analysis once a SGWB signal has been detected. We, however, note that, compared to ref. [87], the aim of the present paper is a slightly different one. Instead of performing a δχ 2 analysis, we actually go one step back and focus on the maximal SNR at which GWs from cosmological phase transitions can be detected in upcoming experiments. In this sense, we hope that the novel concept of PISCs will first and foremost prove to be a helpful tool for the systematic exploration and comparison of the GW phenomenology in different BSM models.
The rest of the paper is organized as follows. In section 2, we will first review the generation of GWs during a SFOPT in the early Universe and collect all expressions that are necessary to compute the signal spectrum Ω signal . In section 3, we will then introduce the concept of PISCs (see section 3.1) and highlight a few possible applications. This will include the detailed discussion of a benchmark point (BP) in a particular SM extension (see section 3.2), the comparison of the GW signals predicted by different BSM models (see section 3.3), a few remarks on how to relax the assumptions underlying the construction of PISCs (see section 3.4), a comment on runaway phase transitions in vacuum (see section 3.5), a discussion of how new theoretical results can be used to update our PISC plots (see section 3.6), and finally the presentation of semianalytical fit functions for JHEP01(2021)097 the exact numerical PISCs of three specific experiments: LISA, the Deci-Hertz Interferometer Gravitational-Wave Observatory (DECIGO) [88][89][90][91], and the Big-Bang Observer (BBO) [92][93][94] (see section 3.7). Section 4 contains our conclusions as well as an outlook on how to extend and generalize our approach in the future.

Gravitational-wave signal from a cosmological phase transition
In quantum field theory, a first-order phase transition is characterized by a scalar field φ (or a set of scalar fields φ i ) experiencing a discontinuous change in its vacuum expectation value, φ false → φ true . In the context of early-Universe cosmology, such a transition manifests itself in the nucleation of bubbles filled by the true vacuum configuration φ true in the ambient plasma, where the scalar field still resides in its false vacuum configuration φ false . The expanding and colliding scalar-field bubbles as well as their interaction with the thermal plasma then result in the production of a primordial SGWB via three different mechanisms: • (b) collisions of expanding bubble walls [123][124][125][126][127][128], • (s) compressional modes (i.e., sound waves) in the bulk plasma [129][130][131][132], and • (t) vortical motion (i.e., magnetohydrodynamic turbulence) in the bulk plasma [133][134][135][136][137][138].
In general, the total GW signal from a cosmological phase transition approximately follows from the linear superposition of the signals stemming from these three individual sources, In this section, we shall give a brief overview of these three signal contributions, following the review reports by the LISA Cosmology Working Group in refs. [24,25], which represent standard references on this subject. For more recent work on the dynamics of cosmological phase transitions and the shape of the resulting GW signal, we refer to refs. [139][140][141][142][143][144][145][146][147][148][149][150][151][152][153][154][155].

JHEP01(2021)097
We begin by pointing out that, among refs. [24,25], only ref. [24] presents semianalytical expressions for all three GW sources, Ω b , Ω s , and Ω t . The discussion in ref. [25] is more conservative in the sense that it solely accounts for the signal from sound waves, which is in many cases much stronger than the signal from bubble collisions and in general much better understood from a theoretical perspective than the signal from turbulence. Ref. [25] also distinguishes between two different expressions for Ω s , depending on whether the time of shock formation in the bulk plasma after the phase transition is longer or shorter than the Hubble time. Meanwhile, it assumes the same spectral shape function for the signal from sound waves as the analysis in ref. [24]. In view of this situation, we decide to consistently base our analysis on the expressions for Ω b , Ω s , and Ω t in ref. [24]. There are two main reasons for this decision: First of all, we intend to demonstrate how to apply our new PISC method in case there is more than just one contribution to the total signal. The point is that we expect to see significant progress on the theory side in the coming years, which will eventually prompt one to go again beyond the conservative approach of ref. [25]. Our analysis thus sets the stage for this moment when the understanding of all three sources has improved and more reliable expressions for Ω b , Ω s , and Ω t have been attained. A second reason is that the focus of our analysis is primarily on the construction of new experimental sensitivity curves; we do not have anything new to say on the theoretical aspects of the expected GW signal. For this reason, we refrain from participating in the ongoing debate on the correct treatment of shock formation and the corresponding energy transfer from sound waves to turbulence. An attractive feature of our new sensitivity curves is that their construction is for the most part anyway independent of these open questions on the theory side. As anticipated in eq. (1.4), our PISCs will only require knowledge of the experimental noise spectra and spectral shape functions, but will be independent of the exact theoretical predictions for the peak amplitudes entering the GW spectrum. Therefore, as ref. [24] and ref. [25] use the same spectral shape function for the signal from sound waves, our sensitivity curves are actually not affected by the different treatment of sound waves in ref. [25]. For our purposes, the only noticeable consequence consists in the fact that the analysis in ref. [25] causes some of the benchmark points to slightly shift in our PISC plots compared to the analysis in ref. [24]. We will comment on these shifts in more detail in section 3.6, where we illustrate how an improved theoretical understanding of the peak amplitudes can be used to update our PISC plots without the need to revise any of the experimental sensitivity curves. At the same time, we stress that, otherwise, any future update of the spectral shape functions will require an update of the experimental PISCs. In fact, we expect that regular updates of our PISC plots would provide a useful means to track the theoretical progress in the field.
The three contributions Ω b , Ω s , and Ω t can be parametrized in a model-independent way in terms of a set of characteristic SFOPT parameters, α, β/H * , T * , v w , κ b , κ s , and κ t ,

JHEP01(2021)097
Let us now go through the different quantities in this equation one by one. The dimensionless energy density parameters Ω i = ρ i /ρ c on the left-hand side (l.h.s.) of eq. (2.2) measure the fractions of the total (critical) energy density ρ c = 3 H 2 0 M 2 Pl that are contained in GWs of a particular physical origin, i ∈ {b, s, t}. Here, H 0 is the Hubble parameter in the present Universe, and M Pl 2.44 × 10 18 GeV denotes the reduced Planck mass. In the following, we will typically multiply all energy density parameters Ω by the square of the dimensionless Hubble parameter h, which is defined via the relation H 0 = 100 h km/s/Mpc. In this way, we make sure that quantities of the form h 2 Ω are not affected by the experimental uncertainty in the Hubble parameter H 0 . The SFOPT parameter α in eq. (2.2) is proportional to the change in the trace of the energy-momentum tensor, ∆T µ µ , across the phase transition [25] Here, ρ v and ∂ρ v /∂T are the temperature-dependent effective potential (i.e., free-energy density) in the scalar sector and its derivative with respect to temperature T , respectively, while ρ r denotes the energy density of relativistic radiation. The subscripts "false" and "true" indicate that these quantities are evaluated in the false and the true vacuum configuration in field space, respectively. The temperature T * is the characteristic temperature at the time of GW production, which roughly equals the temperature at the time of bubble nucleation, T * T n , unless the phase transition occurs in a strongly supercooled state. In this case, the nucleation temperature can be significantly suppressed compared to the reheating temperature after the completion of the phase transition, T n T rh T * . Thus, for strongly supercooled phase transitions, all quantities in eq. (2.3) need to be evaluated at T n rather than T * . The SFOPT parameter β/H * is the inverse of the duration of the phase transition in units of the Hubble time H −1 * at the time of GW production. Formally, it is defined in terms of the derivative of the bounce action S that controls the rate of bubble nucleation, Again, this definition only holds as long as there is no large hierarchy between the temperatures T * and T n . In the case of strong supercooling (i.e., for T n T * ), one has to work instead with β/H * = H n /H * T n dS/dT | T =Tn . The parameter v w in eq. (2.2) represents the velocity of the bubble wall in the plasma rest frame. κ b , κ s , and κ t are three efficiency factors that characterize the fractions of the released vacuum energy that are converted into the energy of scalar-field gradients, sound waves, and turbulence, respectively. It is customary to express κ s and κ t in terms of an efficiency factor κ kin that characterizes the JHEP01(2021)097 energy fraction that is converted into bulk kinetic energy and an additional parameter , i.e., κ s = κ kin and κ t = κ kin . The precise numerical value of is the subject of an ongoing debate in the literature. While some authors estimate to be quite small, 0.05 · · · 0.10 (see, e.g., refs. [24,47,64]), based on the results presented in ref. [132], other authors use values as large as = 1 (see, e.g. ref. [159]). In our analysis, we will set = 0.10. For a more detailed discussion on the efficiency factors κ s and κ t , see also the recent analysis in ref. [148].
In order to estimate the efficiency factors κ b and κ kin , one has to distinguish between three different types of phase transitions: (i) nonrunaway phase transitions in a plasma (NP), (ii) runaway phase transitions in a plasma (RP), and (iii) runaway phase transitions in vacuum (RV). In the first case, the bubble wall velocity saturates at a subluminal value, and the contribution to the total GW signal from bubble collisions is negligibly small, κ b 0. Meanwhile, the efficiency factor κ kin is well approximated by the following fit function [156], Here, we will use in practice the large-v w expression for κ kin for all velocities and similarly, the small-v w expression for κ kin for all velocities v w ≤ v α w . In our numerical analysis in section 3, we will moreover fix the bubble wall velocity at v w = 0.95 for all NP phase transitions. This is in accord with the analysis in ref. [24], where the choice v w = 0.95 simply serves the purpose to increase the strength of the GW signal.
For RP phase transitions, the picture is a slightly different one, as in this case, the energy deposited into the scalar field is no longer negligible. RP phase transitions occur for α values α > α ∞ , where α ∞ marks the threshold value at which the walls of the scalarfield bubbles begin to "run away" at the speed of light, v w = 1. The threshold value α ∞ is model-dependent and follows from the shifts in the bosonic and fermionic masses squared that are induced by the changing scalar-field background across the phase transition [156], Here, g * , g b , and g f represent the effective number of relativistic degrees of freedom (DOFs) in the unbroken phase at the time of bubble nucleation, the internal DOFs of boson species b, and the internal DOFs of fermion species f , respectively. For a strong first-order electroweak phase transition (SFOEWPT) in a model with SM-like particle content, one finds

JHEP01(2021)097
where φ * denotes the scalar field value in the broken phase at the time of bubble nucleation. For simplicity, we shall estimate α ∞ based on eq. (2.8) in the case of all RP phase transitions that we are going to be interested in. In SFOEWPT scenarios where the SM Higgs field is only weakly coupled to the new-physics sector, this is typically a reasonable approximation. The ratio of α ∞ and the actual α value of a RP phase transition determines the efficiency of energy transfer from the vacuum to the scalar field and to the bulk plasma, respectively, RP: where κ therm accounts for the fraction of vacuum energy that is converted to thermal energy (which does not source the production of GWs) and where the parameter κ ∞ is given by (2.10) The expression for κ b in eq. (2.9) has recently been updated in ref. [160] based on the new all-orders calculation in ref. [161]. We will come back to this point in section 3.6. The case of a RV phase transition, finally, corresponds to the limit of very large α, such that the energy transferred to the plasma becomes negligible. In this case, the only remaining free parameters are β/H * and T * . All other parameters are fixed at characteristic values, Having introduced the parameters α, β/H * , T * , v w , κ b , κ s , and κ t , we are now able to spell out how the different contributions to the GW signal in eq. (2.2) can be parametrized in terms of these quantities. We begin by writing down the three peak amplitudes [24], (2.12) Here, we emphasize that the effective number of relativistic DOFs, g * , needs to be evaluated as a function of T * . In our numerical analysis in section 3, we will approximate g * by its SM value, making use of the numerical data tabulated in ref. [162] in order to correctly describe its temperature dependence. The spectral shape functions S b , S s , and S t are given as [24] S

JHEP01(2021)097
Note that S b and S s are normalized to unity at the respective peak frequencies, whereas the value of S t at f = f t depends on what we shall refer to as the Hubble frequency h * , . (2.14) This normalization is merely a matter of convention. Alternatively, one could simply rescale both h 2 Ω peak t and S t by a factor 2 11/3 (1 + 8π f t /h * ), such that S t (f = f t ) = 1. The Hubble frequency h * corresponds to the particular wavenumber k * that equals the Hubble rate H * at the time of GW production. Its redshifted, present-day value is solely controlled by T * , Finally, we state the three peak frequencies [24], which completes our discussion of eq. (2. 2),

Definition
Our idea of peak-integrated sensitivity curves (PISCs) is based on the following observation: For a GW signal from a SFOPT and within the approximations outlined in section 2, the frequency integration in eq. (1.1) becomes independent of any model details. Therefore, by combining the expressions in eqs. (1.1), (2.1), and (2.2), the SNR can always be written as with the partial SNRs i/j on the right-hand side (r.h.s.) of this relation being defined as Here, the numerator is defined as the geometric mean of the corresponding peak amplitudes, while the denominator represents what we will refer to as the peak-integrated sensitivity (PIS),

JHEP01(2021)097
We normalize Ω i/j PIS to an observing time of one average year in the Gregorian calendar, 1 yr = 3.1556952 × 10 7 Hz −1 , (3.5) such that the actual observing time t obs appears as a simple rescaling factor in eq. (3.1). 2 A remarkable property of the six sensitivities Ω i/j PIS is that they can be explicitly computed without ever referring to a particular BSM model. In particular, they can be fully parametrized in terms of the peak frequencies f b , f s , f t and the Hubble frequency h * without ever specifying the values of the SFOPT parameters α, β/H * , T * , v w , κ b , κ s , and κ t , This renders them a handy and model-independent tool for discussing the sensitivity of future searches for GWs from a SFOPT. The sensitivities Ω i/j PIS can especially be used to construct peak-integrated sensitivity curves (PISCs) and peak-integrated sensitivity bands (PISBs) by plotting Ω b PIS as a function of f b , Ω s PIS as a function of f s , etc. Then, once these curves and bands are known, one can fit the exact numerical results by semianalytical fit functions, which allows one to write down quasianalytic expressions for the total SNR as functions of the SFOPT parameters α, β/H * , T * , v w , κ b , κ s , and κ t . In this way, the concept of PISCs and PISBs closes the gap between the numerical modeling of SFOPTs on the theory side and the instrumental properties of future GW searches on the experimental side.
In the remainder of this paper, we will now discuss the idea of PISCs and PISBs in more detail and highlight a few possible applications. In doing so, we will focus on the sensitivities of three proposed space-borne GW interferometers: LISA, DECIGO, and BBO. Among these three experiments, LISA is the most mature one, which was approved by the European Space Agency as its third large-class (L3) mission in 2017. According to the L3 mission concept, LISA will consist of three identical spacecraft in an equilateral triangular formation separated by 2.5 million km and connected by six active laser links [22]. In this configuration, LISA will be able to search for a SGWB signal by performing an auto-correlation measurement [163], which means that we have to set n det = 1 in eq. (3.4). The design concepts of DECIGO and BBO envision, by contrast, a hexagonal configuration of two triangular detectors (i.e., a "Star-of-David"-like configuration). Each of these two experiments will hence effectively represent a two-detector network, which will enable DECIGO and BBO to search for a SGWB signal by performing a cross-correlation measurement. This implies n det = 2 in eq. (3.4).

Benchmark point in the scalar-singlet extension
First, let us illustrate the philosophy behind our PISC method by means of a single benchmark point in a particular SFOPT scenario. To this end, we shall consider the simplest JHEP01(2021)097 example of a BSM model giving rise to GWs from a SFOEWPT, namely, the real-scalarsinglet extension of the standard model (xSM). This model, also known as the (scalar) Higgs portal scenario, features an extra gauge singlet in the scalar sector, which couples to the SM Higgs boson and which may or may not be charged under a Z 2 symmetry. More details on the xSM as well as a more complete list of references are contained in the companion paper [2], where we apply our PISC method to investigate the GW phenomenology of this model. In ref. [2], we include all renormalizable operators in the scalar potential that are allowed by gauge invariance, i.e., we do not require the scalar singlet to be charged under a Z 2 symmetry. However, for the purposes of the following discussion, it will suffice to restrict ourselves to the Z 2 -symmetric formulation of the xSM, which can also result in a SFOEWPT [164]. The xSM benchmark point that we are going to be interested in corresponds to benchmark point B in section 4.2.2 of ref. [24]. It is characterized by the following SFOPT parameter values, and describes a RP phase transition, such that all three GW sources during the phase transition (bubble collisions, sound waves, and turbulence) contribute to the total signal. 3 In the next section, where we will compare different SFOPT scenarios to each other, we will refer to this point as benchmark point #14 (see table 1). The numerical values in eq. (3.7) are all we need to evaluate the expressions for Ω b , Ω s , and Ω t in section 2 and plot the total GW signal as well as its individual contributions as functions of frequency f (see figure 1). In view of figure 1, several comments are in order. In the upper panel of figure 1, we show the strain noise spectra Ω noise of all current and future GW experiments that we consider in this paper. These strain noise spectra, which we review in more detail in appendix A, are the starting point for constructing both power-law-and peak-integrated sensitivity curves. As for the second-generation ground-based interferometers (aLIGO, aVirgo, and KAGRA), we consider three different detector networks that can be formed by these experiments: • Hanford-Livingston (HL): aLIGO Hanford + Livingston Observatories (aLHO + aLLO); • Hanford-Livingston-Virgo (HLV): aLHO + aLLO + aVirgo; • Hanford-Livingston-Virgo-KAGRA (HLVK): aLHO + aLLO + aVirgo + KAGRA. 3 According to the new results in refs. [160,161], the GW signal from bubble collisions is actually strongly suppressed in the xSM (see also refs. [141,148]). We nevertheless stick to the benchmark point in eq. (3.7) and its interpretation as a RP phase transition. On the one hand, this facilitates the direct comparison between our results and the sensitivity plots in ref. [24]. On the other hand, it is expected that similar values of the SFOPT parameters T * , α, β/H * , and φ * /T * can be easily obtained for a RP phase transition in a hidden scalar sector that does not couple to the SM (see, e.g., ref. [45]). In this sense, "xSM" is understood to refer to such a hidden-sector equivalent of the actual xSM in the following (see also our discussion in section 3.6).

JHEP01(2021)097
For the HL and HLVK networks, we indicate the respective design noise spectra [165,166], while for the HLV network, we also indicate, in addition to the design noise spectrum, a noise spectrum representative for observing run 2 (O2) [167][168][169]. This is also reflected in our use of color in figure 1. Projected noise spectra based on sensitivity estimates are represented by simple lines, whereas noise spectra based on existing data are highlighted by a color shading.
In the lower panel of figure 1, we show the power-law-integrated sensitivities Ω PLIS that can be constructed from the strain noise spectra Ω noise according to the algorithm outlined in appendix A [see eq. (A. 29) for the precise definition of Ω PLIS ]. All PLISCs in figure 1 are normalized to an SNR threshold thr = 1; for all future interferometers, we set the observing time to t obs = 1 yr; and for all future PTA experiments, we use t obs = 20 yr. These values are not necessarily realistic (see, e.g., the discussion and references in appendix B of ref. [45]). Our main motivation for setting thr and t obs to these values rather is to guarantee an equal normalization of our power-law-and peak-integrated sensitivity curves. For the purposes of the present paper, this is a reasonable strategy, which allows for a more direct comparison of our PLISC and PISC plots. In addition to the experimental PLISCs, the lower panel of figure 1 also displays the GW signal Ω SFOPT for the xSM benchmark point as well as its three individual contributions, Ω b , Ω s , and Ω t [see eqs. (3.8) The PLISCs in figure 1 convey a useful impression of the different sensitivities of ongoing and planned GW experiments. The PLISCs for LISA, DECIGO, and BBO leave in particular no doubt that these three experiments would have (very) good chances to detect the GW signal from the SFOPT in our xSM benchmark scenario. This is the important, qualitative message of the PLISC plot in figure 1. Its quantitative information content, on the other hand, is somewhat limited. Just by inspecting the GW signal and experimental PLISCs in figure 1, it is, e.g., impossible to precisely infer the expected SNRs in eq. (3.8). To see this, recall that the numerical results in eq. (3.8) follow from the frequency integral over the strain noise spectra in the upper panel of figure 1 [see eq. (1.1)]. It is therefore notoriously difficult to include information on the expected SNR in the lower panel of figure 1. Strictly speaking, the only type of signal for which a PLISC plot lends itself to a statistical interpretation is a pure power law; hence the name. In this case, the factor by which the signal curve needs to be rescaled (i.e., the amount by which it needs to be vertically shifted) in order to align it with an equally sloped tangent of a PLISC can be interpreted as the corresponding SNR (see appendix A). 4 However, for a GW signal from a SFOPT, the assumption of a pure power law is maximally violated in the most relevant part of the spectrum, i.e., close to the dominating peak(s) in the spectrum. This

JHEP01(2021)097
is also nicely illustrated by the signal in the xSM benchmark scenario, which features a double-peak structure at f ∼ f b , f s and thus clearly deviates from a pure power law in the frequency range where the signal strength is the largest.
In general, one should therefore take PLISC plots such as the one in figure 1 with a grain of salt. They typically represent a helpful qualitative visualization; however, for signals that notably deviate from a power law, they no longer contain useful information on the SNR. This observation is one of the main motivations behind our new PISC method. We argue that, as soon as more information on the expected shape of the signal is available, this information should also be made use of in the construction of sensitivity curves. This is exactly the situation in which we find ourselves in the context of GWs from a SFOPT, where the spectral shape of the signal is controlled by S b , S s , and S t in eq. (2.13). Given the large and increasing interest in this type of signal, we therefore deem it justified to construct new sensitivity curves -PISCs -that incorporate knowledge of the expected signal shape.
In figure 2, we present the PISCs for LISA, DECIGO, and BBO [see eq. (3.4)] in combination with the predictions for the peak frequencies and peak amplitudes in the xSM benchmark scenario. For each PISC that depends on more than just one frequency, we fix the remaining frequencies at their respective benchmark values. In each of the six plots in figure 2, we also indicate the respective partial SNRs in the format log 10 ]. These partial SNRs illustrate a characteristic feature of our PISC plots: They retain the full information on the SNR, encoding it on the y-axis. That is, in each of the six plots, the vertical separations between the benchmark point and the PISCs directly correspond to the respective partial SNRs. The total SNRs in eq. (3.8) then follow from adding these partial SNRs in quadrature [see eq. (3.1)]. The fact that our method splits the total SNR into six different contributions also enables one to easily combine and compare these contributions. From figure 2, we can, e.g., read off that LISA will be most sensitive to the b/s-channel, i.e., the overlap of the two signals from bubble collisions and sound waves, while DECIGO and BBO will be most sensitive to the b-channel, i.e., the signal from bubble collisions only, without any inference from other sources. At the same time, all three experiments will be least sensitive to the t-channel, i.e., the signal from turbulence. Similar qualitative conclusions can also be drawn from figure 1, just by looking at the relation of the signal curve and the three sensitivity curves for LISA, DECIGO, and BBO; however, the advantage of our PISC plots in figure 2 is that they make these conclusions quantitatively more precise. Thanks to the partial SNRs in figure 2, it now possible to assess precisely by how much a signal is dominated by a certain contribution or to what extent an experiment will be sensitive to the six individual signal channels. Likewise, it is straightforward to restrict oneself to a particular channel in our approach, while neglecting all others. Suppose, e.g., we were only interested in the signal from sound waves because we wanted to compare our analysis with the one in ref. [25]. In this case, the full information on the projected sensitivities of LISA, DECIGO, and BBO to the GW signal from a phase transition would be encoded in the middle left panel of figure 2, i.e., the plot of the three Ω s PIS curves as functions of f s . We argue that essentially all sensitivity plots presented in ref. [25] can be mapped onto this single PISC plot.

JHEP01(2021)097
It is also interesting to compare our PISC plots to the usual scatter or contour plots of the total SNR as a function of a subset of model parameters, = ({p i }), that one frequently encounters in the literature (see also section 6 in ref. [2]). 5 Such SNR plots typically show the dependence of the SNR on auxiliary SFOPT parameters, such as α, β/H * , T * , etc., on two-dimensional hypersurfaces in parameter space. In principle, one may produce arbitrarily many SNR plots because of the arbitrarily many possibilities to slice the higher-dimensional parameter space, and still, one would fail to really capture the entire available information. In our case, the model parameter space is, by contrast, projected onto only six frequency -amplitude planes without any loss of information. In consequence, we do not need to keep any SFOPT parameters fixed at specific values. In refs. [24,25], most sensitivity plots assume, e.g., a specific phase transition temperature, such as T * = 50 GeV or T * = 100 GeV. In the case of our PISC plots, this is not necessary. All benchmark points are simply projected onto the same six planes, irrespective of their associated T * value. The sensitivity curves in these planes (i.e., our PISCs) are moreover completely independent of the SFOPT parameters. By construction, they only depend on the experimental noise spectra and spectral shape functions in eq. (2.13). In this sense, they represent truly experimental quantities that can be studied without worrying much about questions related to theory and model building. 6 This property goes hand in hand with the fact that our PISCs are formulated in terms of physical observables that are experimentally accessible and that will likely play an important role in the experimental data analysis. SNR scatter and contour plots, on the other hand, do not disentangle experimental from theoretical uncertainties, as they are subject to all uncertainties entering the computation of the SNR. Still they provide useful information from a model builder's perspective. In the end, SNR and PISC plots are complementary to each other, with their combination being the most powerful approach.

Comparison of different BSM models
Another advantage of our PISC plots is that they set the stage for a systematic comparison of the GW phenomenology in different BSM models. To illustrate this point, we shall revisit the full set of BSM models and benchmark points investigated in ref. [24] in this section, recasting their respective predictions in terms of our new PISC method. An equivalent analysis for all models and benchmark points studied in ref. [25] can be found in ref. [86]. The present work and ref. [86] therefore cover together the full set of results in refs. [24,25], illustrating how the sensitivity plots in these two papers can be translated into our PISC language. Comparing the present paper with ref. [24] and ref. [86] with ref. [25] thus highlights how improvements in modeling the GW signal are reflected in the different types of sensitivity plots, i.e., PISC and SNR scatter plots, respectively (see also JHEP01(2021)097  section 3.6). In addition to the xSM, ref. [24] also considers SFOPTs in: (i) the two-Higgs-doublet model (2HDM) [170,171], (ii) the next-to-minimal supersymmetric standard model (NMSSM) [172], (iii) the standard model effective field theory (SMEFT) [173], 7 (iv) a strongly coupled hidden sector giving rise to composite dark matter (DM) [174,175], and an approximately conformal dilaton model in the context of warped extra dimensions [176,177]. A detailed discussion of these models and SFOPT scenarios can be found in refs. [24,25], which we shall not reproduce here. For our purposes, the entire relevant data describing these models is contained in table 1, where we list the values of T * , α, β/H * , and φ * /T * for all benchmark points that we shall include in our analysis. For each phase transition, we also indicate whether we expect it to be of NP, RP, or RV type (see section 2). Here, note that, for the NMSSM and SMEFT, we consider both NP and RP phase transitions, which explains why the corresponding benchmark points are duplicated in table 1. As for the DM and dilaton models, more quantitative analyses of the phase transition dynamics are still pending. The corresponding benchmark points in table 1 therefore represent educated guesses rather than precise numerical predictions.
Based on the data in table 1, we are able to repeat our analysis in the previous section for all 33 benchmark points that we are interested in. First of all, we plot again the total GW signal and its three individual contributions as functions of frequency for all JHEP01(2021)097 33 benchmark points. This results in the busy plot in figure 3, which highlights another limitation of the standard PLISC approach. PLISC plots such as those in figures 1 and 3 are not well suited for comparing a large number of spectra to each other. This is also part of the reason why most PLISC plots in the literature only show a handful of spectra. As soon as one intends to study O (10) or more spectra at the same time, the PLISC approach becomes highly impractical. One may argue that a possible way out of this problem might consist in restricting oneself to just plotting points of the form f i , Ω peak i . In this case, our busy collection of signal curves would reduce to a simple scatter plot that could be compared more easily to the various PLISCs in figure 3. Indeed, this is a strategy that one sometimes encounters in the literature. We, however, argue that such an approach corresponds to comparing apples and oranges. PLISCs do represent useful sensitivity curves; but there is no reason to believe that they are also automatically the best choice for indicating experimental sensitivities to an ensemble of peak frequencies and peak amplitudes, a purpose they are not specifically designed for. In the case of GWs from a cosmological phase transition, PLISC plots simply no longer encode information on the expected SNR (see our discussion in section 3.2), irrespective of whether one decides to combine them with a busy collection of signal curves or a f i , Ω peak i scatter plot. Our PISC method amounts, by contrast, to comparing apples and apples. That is, our PISCs are constructed in exactly such a way that they represent the optimal sensitivity curves to be used in scatter plots of peak frequencies and peak amplitudes. This is illustrated in figure 4, where we now combine our PISCs for LISA, DECIGO, and BBO with the predictions of the benchmark points listed in table 1. In contrast to figure 3, the plots in figure 4 are much easier to read, while at the same time, they still encode the full SNR information on the y-axis. A slight difference between figure 2 and figure 4, though, is that we are now no longer able to draw sensitivity curves for all six signal channels; for all channels involving the signal from turbulence, we have to draw sensitivity bands. The reason for this is the dependence of the shape function S t on the Hubble frequency h * [see eq. (2.13)], whose value is independent of the parameter combination β/H * /v w , unlike the values of f b , f s , and f t [see eqs. (2.15) and (2.16)]. In figure 4, the different parameter dependence of h * and the three peak frequencies is accounted for by the width of the PISBs, which reflects the variation of the ratio h * /f t in our data set. We also point out that, as a consequence of this finite width, it is not clear whether benchmark points inside PISBs do or do not have a partial SNR larger than one. To remedy this shortcoming, we distinguish between "empty" and "filled" points in figure 4, which respectively correspond to partial SNRs for LISA, LISA i/j , smaller or larger than one. As in the previous section, let us also compare our results to the standard approach of SNR scatter and contour plots. In contrast to these standard plots, our PISC plots in figure 4 do not require contour lines or a color code to indicate the expected SNR. This provides us with the freedom to use a color code for distinguishing between the predictions of different models. Similarly, the fact that we do not rely on SNR contour lines in figure 4 allows us to plot and compare the PISCs and PISBs of three future experiments at the same time. We argue that such a simultaneous comparison of different models and different experiments would be significantly more difficult if one were to work with SNR plots only.

JHEP01(2021)097
Finally, we stress once more that it would be trivial to restrict oneself to individual signal channels in our analysis. For instance, if one were interested in the signal from sound waves only, the entire relevant information would be readily contained in the middle left panel of figure 4.

Relaxing some of the underlying assumptions
The sensitivity curves and bands in figures 2 and 4 are constructed in such a way that they have a natural interpretation in terms of the expected SNR. Consider, e.g., a benchmark point that is separated from the PISC in the i/j-channel by a multiplicative factor ∆y along the y-axis. Thanks to the definitions and conventions adopted in section 3.1, this benchmark point predicts a partial SNR of exactly i/j = ∆y. Benchmark points that directly lie on a PISC correspondingly predict i/j = 1. Each of our PISCs is therefore canonically normalized to a unit threshold value for the corresponding partial SNR. Similarly, all of our PISCs and PISBs are normalized to an observing time of t obs = 1 yr, which we simply choose for convenience. Our PISBs rely in addition on the explicit expressions for the Hubble and peak frequencies in eqs. (2.15) and (2.16). In this section, we shall now demonstrate how our plots can be generalized if one is interested in relaxing some of these assumptions, i.e., if one is interested in a different normalization or a different relation among the Hubble and peak frequencies.
In a first step, let us generalize eqs. (3.1), (3.2), and (3.4) to arbitrary thr and t obs , , where thr now denotes the partial SNR threshold for a single channel. The generalized peak-integrated sensitivity in the i/j-channel can thus be interpreted as the minimal peak amplitude that is necessary to reach, after an observing time t obs , a partial SNR of thr in this channel. A hypothetical scenario predicting peak amplitudes that just reach the respective threshold sensitivities in all six channels would therefore predict a total SNR of = √ 6 thr . As evident from eq. (3.9), the generalized sensitivities scale as follows with thr and t obs , simply shifting these base curves up or down, i.e., by rescaling them by a factor thr . This needs to be compared to the situation in the case of SNR scatter and contour plots, which typically do not exhibit a universal relation between shifts in the expected SNR, ∆ , and shifts in one of the coordinate directions, ∆x or ∆y. The second message of figure 5 is that our PISC plots are reminiscent of plots that one often encounters in other fields of experimental physics, such as, e.g., the standard sensitivity plots for DM direct-detection experiments or the Brazil-band plots that were often shown by the experiments at the Large Hadron Collider prior to the discovery of the SM Higgs boson. What we mean by this is that our PISC plots show, in a manner intuitive for particle physicists, how future GW experiments will approach from above and cut into the signal regions of specific SFOPT scenarios. As in the case of the DM and Higgs plots, a better sensitivity reach of an experiment corresponds to a PISC extending to lower values along the y-axis in our plots. Thanks to the scaling with the observing time t obs in eq. (3.10), the expected experimental progress over the years can in particular be pictured as pushing our PISCs further and further down in the vertical direction. Again, the situation in the case of SNR scatter and contour plots is quite different, as these plots typically do not exhibit any simple scaling relation with respect to the observing time t obs .
Next, let us briefly discuss an alternative presentation of the sensitivity bands in figure 4. To this end, suppose that future progress will lead to a better theoretical understanding of the peak frequencies f b , f s , f t , in particular, to the realization that f t is in fact not simply proportional to h * times one power of the parameter combination β/H * /v w . In this case, it will be helpful to have a proper understanding of the full dependence of Ω t PIS on both frequencies, f t and h * [see eq. (3.6)]. Therefore, instead of constructing sensitivity bands as in figure 4, one may as well work with contour plots of Ω t PIS in the f th * plane (see figure 6). Similar contour plots can also be drawn for the sensitivities in the b/t-and s/t-channels. The advantage of this alternative presentation is that the Hubble frequency h * now no longer corresponds to a hidden parameter. In figure 4, it is unclear by JHEP01(2021)097  Hubble frequency h* [Hz] construction with which curves in a PISB one should respectively compare the individual benchmark points. The plots in figure 6 offer a trivial solution to this problem, as they explicitly feature h * on the y-axis. The disadvantage of this method, however, is that it requires contour lines to indicate the experimental sensitivity. In this sense, the partial SNR is now encoded in the z-direction, which means that one looses some of the attractive properties of our PISC plots. In figure 6, we use, e.g., empty and filled symbols to indicate which benchmark points predict a partial SNR smaller or larger than one. However, beyond that, the plots in figure 6 contain no further information on the expected SNR.
In particular, it is unclear how far the individual points are separated from the threshold sensitivity at their respective locations in the f th * plane.

Runaway phase transitions in vacuum
The only free parameters in the case of a RV phase transition are β/H * and T * (see section 2). All other SFOPT parameters are fixed at the values listed in eq. (2.11). At the same time, the only relevant source of GWs during a RV phase transition are bubble collisions. For this type of phase transition, there is hence a one-to-one correspondence between β/H * and T * on the one hand and the peak frequency and peak amplitude in the b-channel on the other hand. This relation is shown in figure 7, where we overlay the PISC plot for the signal from bubble collisions with contours of constant β/H * and T * . As can be read off from this plot, LISA will be sensitive to β/H * values as large as β/H * ∼ 10 4 if T * is close to 10 GeV, while DECIGO and BBO will be able to probe RV phase transitions up to β/H * ∼ 10 5 for temperatures T * in the same range. Our results for LISA in figure 7 are consistent with figure 6 in ref. [24], which shows LISA's sensitivity to the signal from bubble collisions in the β/H * -T * plane. The main difference between figure 7 and figure 6 in ref. [24], however, is that our PISC plot encodes the full SNR information on the y-axis. Despite the fact that both plots are related by nothing but a simple parameter transformation, this is a distinct advantage of our figure 7. The plot in figure 7 represents a simple example that demonstrates how information on the underlying SFOPT parameters can be included in PISC plots. Further, more sophisticated examples can be found in the companion paper [2], where we use our PISC method for a comprehensive analysis of the GW signal in the xSM. In ref. cially show how PISC plots can be used to investigate the dependence of the signal on the SFOPT parameters as well as on model parameters such as particle masses or coupling constants. In addition, we use our PISC plots as a starting point for constructing distribution functions (i.e., histograms) in the space of peak frequencies and peak amplitudes. These distribution functions provide a useful tool to characterize the GW phenomenology of the xSM, in particular, when combined with data on the underlying model parameters. In future work, it would be interesting to repeat our analysis in ref. [2] for a broad class of BSM models. Such a global analysis, including PISC plots and histograms such as those in ref. [2], would allow for a comprehensive model comparison at both the qualitative and quantitative level.

New signal predictions after theory updates
In this paper, we apply our PISC approach to all models and benchmark points in ref. [24], which allows us to demonstrate how to use our method if the total GW signal receives three individual contributions, Ω b , Ω s , and Ω t (see the discussion in section 2). An equivalent analysis for the models and benchmark points in ref. [25], which only considers the GW signal from sound waves, can be found in ref. [86]. In this section, we shall now demonstrate how the different modeling of the signal in refs. [24,25] is reflected in our PISC plots. As we shall see, this highlights another advantage of our method: Any update in the theoretical description of the signal that does not affect the spectral shape functions in eq. (2.13) leaves our sensitivity curves invariant. An improved theoretical understanding of the peak amplitudes and peak frequencies merely shifts the individual benchmark points in our plots. 8 This needs to be compared to the SNR scatter plots in refs. [24,25] as well as to the plots generated by the online tool PTPlot, which crucially depend on the precise JHEP01(2021)097 modeling of the GW source. At the same time, we emphasize once more that any future update of the spectral shape functions in eq. (2.13) will require an update of the PISCs presented in this paper.
An important difference between refs. [24,25] is that ref. [25] discards the GW signal from bubble collisions (as well as the GW signal from turbulence, which still requires a better theoretical understanding). The reason for this is the realization [141,148,160,161] that the efficiency factor κ b can be significantly suppressed compared to earlier estimates [see eq. (2.9)]. As pointed out in ref. [141], soft particle emission by particles passing through the bubble walls (so-called transition radiation or splitting processes) leads to an additional pressure acting on the bubble walls in proportion to the wall Lorentz factor ∆P split ∝ γ n . While the next-to-leading-order calculation in ref. [141] found a scaling exponent n = 1, the all-orders resummation in ref. [161] recently demonstrated that, in reality, n = 2. 9 Based on this result, ref. [160] obtained the following updated expression for the efficiency factor κ b , where γ eq denotes the Lorentz factor that is reached when all forces acting on the bubble walls have equilibrated, so that the walls no longer accelerate, andγ * is the Lorentz factor that the bubble walls would reach by the time of collision if there was no friction whatsoever, Here, ∆V is the change in the effective potential across the phase transition, ∆P stands for the leading-order friction term, and R 0 and R * denote the average bubble radius at the time of nucleation and collision, respectively. In many models, one finds that -unless the phase transition is strongly supercooled or transition radiation significantly suppressed (e.g., because the symmetry-breaking field does not couple to gauge bosons)γ eq is typically much smaller thanγ * . This means that many phase transitions that were originally believed to be of the RP type are actually NP phase transitions, where the bubble walls reach a terminal velocity, γ * = min {γ * , γ eq } = γ eq , and the efficiency factor κ b is strongly suppressed, κ b 1. The improved understanding of κ b affects benchmark points #08 to #15, #18, and #19. In principle, it would be desirable to reconsider all of these points, making use of the updated expression in eq. (3.11). However, this is complicated by the fact thatγ * depends on the initial bubble radius R 0 , which requires knowledge of the initial profile and Euclidean action of the symmetry-breaking scalar field at the time of nucleation [148]. In the companion paper [2], we perform such an analysis for the xSM, which confirms that κ b 1 across the entire parameter space. Therefore, instead of explicitly computing R 0 and reevaluating κ b for points #08 to #15, #18, and #19, we simply conclude that all   Figure 8. Update of the benchmark points in ref. [24] according to refs. [34,160,161]. See text.

RP phase transitions in table 1 should either be ignored or replaced by an equivalent NP phase transition.
A second difference between refs. [24,25] is that ref. [25] accounts for the formation of shocks at some time τ sh after the phase transition. If this happens within a Hubble time, H * τ sh < 1, the GW signal from sound waves picks up an extra suppression factor [34,180] (3.13) where H * τ sh can be computed in terms of the enthalpy-weighted root-mean-square of the plasma velocity,Ū f , which follows from the kinetic-energy fraction of the bulk plasma, K, (3.14) In figure 8, we summarize how the different treatment of κ b and Ω peak s in ref. [25] compared to ref. [24] affects our PISC plot in the s-channel. In this figure, we no longer show the predictions of the RP phase transitions in the NMSSM and SMEFT; we replace the RP phase transitions in the xSM by equivalent NP phase transitions; and we rescale all peak amplitudes by the suppression factor in eq. (3.13). These three steps remove and shift some of our benchmark points. However, for the purposes of this paper, the main message of figure 8 is that the three PISCs simply remain the same as in the middle left panel of figure 2, despite the comprehensive theory update. Of course, an update of the spectral shape functions would require a revision of the sensitivity curves. A more extensive version of figure 8, including almost 4000 benchmark points for ten different particle physics models, can be found in ref. [86]. Both plots will continue to serve as useful resources in the future when new theoretical predictions for f s and Ω peak s should become available (see, e.g., refs. [157,158]).

Semianalytical fit functions
In the previous sections, we studied 18 different sensitivities: We considered three different experiments (LISA, DECIGO, and BBO), and for each experiment, we constructed PISCs in six different channels (b, s, t, b/s, b/t, s/t). We shall now conclude our analysis by providing semianalytical fit functions for all of these 18 sensitivities. In doing so, let us assume that all sensitivities can be reasonably well approximated by power series of the following form, for an appropriate set of 4-tuples (a, b, c, d) ∈ R 4 . Based on this ansatz, we are then able to fit our numerical data and determine the coefficients c (a,b,c,d) for each of our 18 sensitivities. Below, we present our results for LISA, DECIGO, and BBO using the following notation, (3.16) Our fit functions are constructed such that they reproduce our numerical results to high precision across the entire range of relevant frequencies. This is, e.g., illustrated in figure 9, where we compare our numerical results and fit functions in the b-, s-, and t-channels to each other. In the other three channels, our fit functions are of an equally high quality.

LISA:
The above fit functions allow us to write down a quasianalytic expression for the SNR, for LISA, DECIGO, and BBO, respectively. In this sense, the concept of PISCs in combination with the fit functions presented in this section amount to a quasianalytic solution to the problem of computing the SNR for the GW signal from a cosmological phase transition. As evident from eq. (3.35), this analytic solution only depends on the SFOPT parameters α, β/H * , T * , v w , κ b , κ s , and κ t ; the frequency dependence of the signal as well as the experimental noise spectra are already take care of by our fit functions. Our result in eq. (3.35) therefore renders any further numerical step (i.e., integration) in the computation of obsolete. With the results in this section at hand, the SNR can be evaluated analytically. This is in particular also true if one is only interested in the signal from a single source. Suppose, e.g., we were only interested in LISA's sensitivity to the signal from sound waves. In this case, the full information on the expected SNR will be contained in the following expression, LISA s = t obs 1 yr where the numerical values of the coefficients c −4 , c −3 , etc. can be read off from eq. (3.18). This is an important result of our analysis. We stress that it is independent of the explicit form of Ω peak s , such that it can be compared to the results in both ref. [24] and ref. [25].

JHEP01(2021)097 4 Conclusions and outlook
For a DM direct-detection experiment, it is typically straightforward to answer the following two questions: (i) "What is the experiment's sensitivity to the DM cross section σ DM as a function of the DM mass m DM ?" (ii) "To what extend will the experiment explore the parameter regions preferred by specific DM scenarios?" This situation in the field of DM experiments needs to be compared to searches for GWs from a SFOPT in the early Universe. In this case, one can likewise ask: (i) "What is an experiment's sensitivity to the GW signal Ω SFOPT as a function of the GW frequency f ?" (ii) "To what extend will it explore the parameter regions preferred by specific SFOPT scenarios?" In contrast to DM searches, there are at present no commonly accepted answers to these questions in the GW community. Existing approaches, such as PLISC and SNR plots, certainly do convey a useful impression of the sensitivities of current and future experiments, but still suffer from a number of shortcomings. Graphical analyses based on PLISCs no longer contain information on the expected SNR and are only meaningful as long as the expected signal does not deviate too much from a pure power law. SNR plots, on the other hand, present the reach of GW experiments in terms of auxiliary parameters instead of observable quantities, are often times restricted to two-dimensional slices through the higher-dimensional parameter space, and require additional elements such as a color code or contour lines to indicate the expected SNR. In this sense, neither of these approaches is truly on par with the sensitivity curves for DM experiments. This observation is the basis for our PISC proposal. Suppose, e.g., one asked: "What is LISA's projected sensitivity to GWs from sound waves?" We argue that the best possible answer to this question would be LISA's The PISCs constructed in this paper exhibit twelve characteristic features (see also section 6 of ref. [2]). They (i) retain the full information on the SNR, encoding it on the y-axis of the PISC plots; (ii) do not require extra graphical elements such as a color code or contour lines to indicate the expected SNR; (iii) do not require one to slice the parameter space into hypersurfaces; (iv) only depend on the experimental noise spectra and spectral shape functions in eq. (2.13), which renders them insensitive to the theoretical uncertainties in calculating the peak amplitudes in eq. (2.12); (v) can be generalized to arbitrary shape functions S; (vi) indicate sensitivities in terms of observables that will play an important role in the experimental data analysis rather than auxiliary quantities; (vii) directly illustrate how GW experiments will approach from above and cut into the signal regions of specific models; (viii) set the stage for the systematic analysis of underlying model-parameter dependences (see the example study in ref. [2]); (ix) allow for an easy comparison of the projected sensitivities of different experiments; (x) allow one to combine and compare different signal contributions at one's convenience; and (xi) close the gap between experiment and theory by eliminating the model-independent and redundant frequency integration in eq. (1.1). It is also easy to (xii) approximate our PISCs by fit functions, such that the SNR can be written as a function of the SFOPT parameters only. In this sense, our PISC method provides a quasianalytical solution to the problem of computing the SNR for the GW signal from a SFOPT. At the same time, we stress the JHEP01(2021)097 important caveat that the PISCs presented in this paper are always only as good as our knowledge of the spectral shape functions in eq. (2.13). Each update of these functions will require an update of our sensitivity curves. This opens up the possibility to use PISC plots as a bookkeeping tool to keep track of the theoretical progress in the field. Similarly, the sensitivity curves constructed in this paper can be generalized to sensitivity bands indicating the theoretical uncertainty in the expected shape of the GW spectrum [86].
There are several natural directions in which our analysis in this paper could be extended. An obvious extension would, e.g., consist in applying our method to further experiments. To facilitate such an analysis, we review the strain noise spectra of a number of interferometer and PTA experiments in appendix A. Beyond that, our method could also be extended to experiments that we do not consider in appendix A, such as, e.g., AEDGE [182], AIGSO [183,184], AION [185], AMIGO [186], Taiji [187], TianGO [188], TianQin [189,190], etc. Furthermore, it would be desirable to repeat our xSM analysis in the companion paper [2] for as many BSM models as possible. The ultimate goal of such an effort would be a comprehensive database providing the necessary means for a systematic and quantitative model comparison. Similarly as in the case of DM experiments, such a database would allow one to construct the signal regions for various models and to illustrate, by means of our PISC plots, how these signal regions are going to be probed by future experiments. The plots in figure 4 provide a glimpse of how such a simultaneous comparison of different models and different experiments could eventually look like. However, to be able to make stronger and more quantitative statements, it will be necessary to consider a significantly larger number of benchmark points for each model. In ref. [2], e.g., we study roughly 6000 points, which is necessary to fully chart and trace out the signal region of the xSM in our PISC plots. Finally, we point out that our PISC method could also be extended to any SGWB signal whose spectral shape is described by a clearly defined shape function S. In this case, one would be able to construct shape-integrated sensitivity curves (SISCs) in analogy to our PISCs. We leave a more detailed discussion of this possibility for future work. Instead, we conclude by stressing that the novel concept of peak-integrated sensitivity curves bears the potential to develop into a new useful standard tool for model builders, phenomenologists, and experimentalists that are interested in the GW signal from a phase transition in the early Universe.

A Review: noise spectra of interferometer and pulsar timing experiments
The construction of both power-law-and peak-integrated sensitivity curves requires knowledge of the experimental strain noise power spectra. In this appendix, we shall therefore review the strain noise spectra of LISA, DECIGO, and BBO, and in addition, several other interferometer and PTA experiments. Specifically, we are going to consider: aLIGO, aVirgo, KAGRA, CE, ET, LISA, DECIGO, BBO, NANOGrav, PPTA, EPTA, IPTA, and SKA. As for the second-generation ground-based interferometers (aLIGO, aVirgo, and KA-GRA), we will derive the effective strain noise spectra of three different detector networks: HL, HLV, and HLVK (see section 3.2). Our analysis in this appendix is supposed to facilitate the generalization of our PISC method to experiments beyond LISA, DECIGO, and BBO. In addition, we hope that it will serve as a useful resource for a broader range of applications.
In section A.1, we will first introduce the necessary formalism and fix our conventions. In section A.2, we will then present all relevant transfer functions (i.e., signal response and overlap reduction functions), before turning to the individual detector noise spectra in section A.3. In section A.4, we will finally put everything together and construct the strain noise spectra. For an overview of all quantities playing a role in the following discussion, see table 2.

A.1 Formalism
We are interested in experimental searches for a stochastic, Gaussian, stationary, isotropic, and unpolarized GW background. A detailed review of the formalism to describe such searches can be found in ref. [191]. In the following, we will adopt the conventions of ref. [191], but restrict ourselves to a briefer exposition and customize our notation. Most SGWB searches aim at measuring a nonzero cross-correlation signal in the outputs of two detectors whose intrinsic noise spectra are uncorrelated [15,83,84]. Let us now outline the main quantities entering the description of such a measurement and derive the expected SNR.
The raw data d I of a single detector I amounts to a time series output of the form which receives a signal contribution s I and a noise contribution n I . Here, the signal contribution represents the detector response to the incoming GWs, which depends on both the properties of the GWs and the geometry of the detector. In the following, large parts of our discussion will refer to the frequency domain, where d I , s I , and n I are replaced by their Fourier transformsd I ,s I , andñ I . In our convention, these Fourier modes are defined via Here, the factor 1/2 reflects the fact that, in our convention, all power spectra are defined to be single-sided. The variance of the detector noise, σ 2 I , can therefore be written as follows, with D I noise being integrated over the frequency range [0, +∞) rather than (−∞, +∞).

JHEP01(2021)097
Next, let us consider a network of detectors, I, J = 1, 2, · · · , and derive a similar integral representation for its response to the signal. In this case, we now have to compute the covariance matrix C IJ = s I s J − s I s J = s I s J instead of just a single variance. The two-point correlation function s I s J accounts again for both the properties of the SGWB and the geometry of the detector network. The incoming GWs are described by tensor perturbations of the spacetime metric g µν . In transverse-traceless gauge, we can write where the tensor perturbations h ij can be decomposed into plane waves as follows, Here, h p n (f ) denotes the amplitude of a sinusoidal plane GW with frequency f , polarization p, and propagation direction n, while (e p n ) ij is the corresponding polarization tensor, Analogous to eq. (A.3), the quadratic expectation value of the Fourier modes h p n reads where the Kronecker delta and Dirac delta functions account for the fact that we assume the SGWB to be unpolarized, isotropic, and stationary. The nontrivial information on the r.h.s. of eq. (A.8) is encoded in the GW strain power spectrum S signal , which describes the total strain power of the SGWB summed over both polarization states, p = +, ×, and integrated over the entire sky as a function of frequency. Analogous to eq. (A.4), the strain power spectrum can be used to write down an integral representation of the strain variance, 11 The signal response s I of detector I follows from convoluting the tensor perturbation h ij with the detector's impulse response R ij I . In the frequency domain, this results iñ The factor of 2 on the r.h.s. of this relation follows from the normalization of the polarization tensors in eq. (A.7) and is therefore nothing but a matter of convention. It would be straightforward to avoid this factor by performing the rescaling: 1/ √ 2 (e p n )ij → (e p n )ij, √ 2 h p n → h p n , and 2 S signal → S signal . In passing, we also mention that this factor of 2 is sometimes ascribed to the fact that the metric perturbations hij receive contributions from two different polarizations (see, e.g., refs. [15,17]). However, this statement is slightly misleading since the strain power spectrum is already summed over both polarization states [191]. That is, if we wrote S signal = S + signal + S × signal , the contributions for both polarizations would still feature a factor of 2.

JHEP01(2021)097
where the response function R p n,I (f ) represents the signal response of detector I to (i.e., the convolution of its impulse response with) a sinusoidal plane GW with frequency f , polarization p, and propagation direction n. The graph of R p n,I as a function of n describes the antenna pattern of the detector for GWs with frequency f and polarization p. More details on the impulse response and the response function can be found in ref. [191]. For our purposes, the important message from eq. (A.10) is that it allows us to write the quadratic expectation value of the signal modess I in a similar way as the expectation value in eq. (A.3), where we introduced the so-called overlap reduction function Γ IJ of the detector pair IJ, As evident from eq. (A.11), the overlap reduction function acts as the transfer function between the GW strain power spectrum and the signal response cross power spectrum of the IJ detector pair, C IJ = Γ IJ S signal . In line with our assumption of an isotropic and unpolarized SGWB, Γ IJ is defined as the sky-and polarization-averaged product of the response functions for the detectors I and J. In the special case of a single detector, J = I, it reduces in particular to the sky-and polarization-averaged square of the antenna pattern, which is also known as the detector transfer function or simply signal response function R I . The function R I relates the GW strain power spectrum S signal to the signal response auto power spectrum D I signal and can be used to define the strain noise auto power spectrum S I noise , (A.14) Often times, one also works with the normalized overlap reduction function which is normalized to γ IJ (f = 0) = 1 for a pair of identical, co-located, and co-aligned interferometers with an opening angle δ between the two arms. Below, we will be interested in the cases δ = π/2 (aLIGO, aVirgo, KAGRA, CE) and δ = π/3 (ET, LISA, DECIGO, BBO), for which the conversion factor between γ IJ and Γ IJ amounts to 1/5 and 3/20, respectively. The relation in eq. (A.11) finally allows us to write down the covariance matrix, Again, we only integrate over positive frequencies because all power spectra are single-sided.

JHEP01(2021)097
We are now all set to derive the expected SNR for a cross-correlation measurement of the SGWB. The basic idea is to apply a filter function Q IJ to the cross-correlation signal S IJ that can be constructed from the data streams of the two detectors I and J, (A. 17) and to choose (i.e., match) this filter function so as to maximize the corresponding SNR, The solution of this optimization problem has been worked out in refs. [83,84]; here, we will only state the final result. In Fourier space, the optimal filter function Q IJ turns out to be where N is an irrelevant normalization constant that cancels in the expression for the SNR in eq. (A.18). Based on this result, the optimal SNR ends up acquiring the following form, where n det = 2 counts the number of detectors involved in the cross-correlation measurement and the frequency interval [f min , f max ] defines the bandwidth of the IJ detector pair. This result is valid in the weak-signal regime, which assumes that the integrand of the frequency integral in eq. (A.20) is smaller than unity for all frequencies, Γ 2 IJ S 2 signal D I noise D J noise . In view of eq. (A.20), several comments are in order. First of all, note that both the optimal filter and the optimal SNR depend on the strain power spectrum of the signal, S signal . In principle, one would therefore need to know the exact shape of the signal that one intends to measure if one really wanted to identify a cross-correlation signal S IJ whose SNR matches the one in eq. (A.20). This is of course impossible, which is why, in practice, one has to resort to a library of template spectra. Among these template spectra, the best approximation of the true signal will then result in the SNR value closest to the optimal SNR. A second comment is that eq. (A.20) remains in fact valid if we consider an idealized auto-correlation measurement in a single detector rather than a cross-correlation measurement using a pair of detectors. In the case of LISA, one will, e.g., be able to monitor the detector noise in real time. For an auto-correlation measurement, the optimal SNR after perfect noise subtraction is then again given by eq. (A.20), however, with n det set to n det = 1 (see ref. [85] and references therein). A third comment finally is that it is straightforward forward to generalize eq. (A.20) to an entire network of detectors, I, J = 1, 2, · · · . In this case, one simply has to compute the partial SNRs for all possible pairs of detectors and add them in quadrature, Given this expression for , it is convenient to define an effective strain noise power spectrum which generalizes the idea of the strain noise power spectrum S I noise to a detector network. With this definition, the SNR can now be written as follows, Finally, both the strain power spectrum of the signal, S signal , and the noise power spectrum of the detector network, S eff noise , can be expressed in terms of GW energy density spectra, Here, the GW energy density spectrum Ω signal is defined as the energy density contained in GWs per logarithmic frequency interval and normalized to the critical energy density ρ c , Note that f ref cancels in this definition, which renders Ω PLIS insensitive to the exact choice of this auxiliary quantity. The interpretation of the power-law-integrated sensitivity (PLIS) is as follows: Any power-law signal that intersects the PLISC, such that Ω signal (f ) > Ω PLIS (f ) for at least some frequency f , results in an SNR above threshold; all curves tangential to the PLISC result in an SNR of exactly thr ; and all curves that always stay below the PLISC have subthreshold SNR. The factor by which the signal curve needs to be rescaled in order to align it with a tangent of the PLISC can thus be interpreted as the expected SNR.
It is also interesting to compare the power-law-integrated sensitivity in eq. (A.28) to the peak-integrated sensitivity in eq. (3.9). Obviously, we can reproduce the expression in eq. (A.28) by setting i = j and S i = S j = (f /f ref ) p in eq. (3.9), i.e., if we assume a spectral shape function that is described by a pure power law. In the case of GWs from a cosmological phase transition, for which better estimates of the spectral shape exist, this is certainly not the best choice. Among other things, this is an important motivation for our PISC method.

A.2 Transfer functions
Signal response function for a single detector: let us now turn to the transfer functions (i.e., signal response and overlap reduction functions) of specific experiments. 12 We begin by computing the signal response function R I for a single equal-arm Michelson interferometer. A closed analytic expression for R I does unfortunately not exist; however, an explicit integral representation can be found in ref. [195], For recent work on transfer functions for GW experiments, see refs. [192][193][194]. These papers also discuss the response to vector and scalar GW polarization states, which only occur in models of modified gravity.

JHEP01(2021)097
where δ denotes again the interferometer's opening angle and u θ , u θ , and θ are defined as u θ = u cos θ , u θ = u cos θ , cos θ = cos δ cos θ + cos sin δ sin θ . Instead of the FSR frequency f fsr , one sometimes also encounters f * = f fsr /π, which is referred to as the transfer frequency [196]. Below, we will be interested in CE's and LISA's signal response functions. Given their current design concepts, we find the following transfer frequencies for these two interferometers, The integral function I (u, δ) in eq. (A.30) does not admit a closed analytic form and needs to be evaluated numerically. In figure 10, we present our numerical results for both δ = π/2 and δ = π/3. As expected, the signal response function approaches a constant value in the small-frequency limit for both opening angles, R I → 1/5 sin 2 δ [see the discussion below eq. (A.15)]. At larger frequencies, R I is subject to sinusoidal modulations with period f fsr = πf * , while its overall amplitude drops off like 1/f 2 . Ignoring the oscillations at high frequencies, R I can be well approximated by rational fit functions for both CE and LISA, Here, we multiplied R LISA by an extra factor of 2 to account for the fact that LISA's six active laser links will allow one to construct two independent data streams at low frequencies that can be used for an auto-correlation measurement [196]. This factor is sometimes missed in the literature. CE, on the other hand, envisions an L-shaped interferometer similar to aLIGO, in which case only one data channel will be available for an auto-correlation measurement.
Overlap reduction functions for detector pairs in a detector network: a crosscorrelation measurement using a detector network requires knowledge of the overlap reduction functions Γ IJ for all detector pairs that can be formed within the network. For Lshaped ground-based interferometers, it is possible to derive an analytic expression for Γ IJ ,  Table 3. Geometrical data describing the angular separation and relative orientation of the six detector pairs that can be formed within the HLVK detector network (see refs. [197,198] for details).
or alternatively, the normalized overlap reduction function γ IJ [see eq. (A.15)] [197,198], where j n (n = 0, 2, 4) represents the spherical Bessel function of the first kind of order n. The frequency dependence in eq. (A.34) is encoded inū = f /f * , wheref * is now defined as In this expression, d ⊕ 12 742 km denotes the mean diameter of the earth, and β is the angle between the detectors I and J in geocentric coordinates. Similarly, ∆ and δ quantify the orientation of the two detectors relative to the great circle on which they both lie. For more details on the definitions of β, ∆, and δ, see ref. [197], in particular, figure 5 in this work. Eq. (A.34) can be used to compute the overlap reduction functions for all six detector pairs that can be formed in the four-detector network consisting of aLHO, aLLO, aVirgo, and KAGRA: aLHO-aLLO (HL), aLHO-aVirgo (HV), aLLO-aVirgo (LV), aLHO-KAGRA (HK), aLLO-KAGRA (LK), and aVirgo-KAGRA (VK). The values of the geometrical angles β, ∆, and δ for these six detector pairs are listed in refs. [197,198] (see  table 3). The resulting normalized overlap reduction functions are shown in the left panel of figure 11.
For all other detector pairs that we are interested in (i.e., for experiments that are not based on the cross-correlation of L-shaped ground-based interferometers), we compile explicit numerical results for the respective overlap reduction functions from the literature: • ET will consist of three V-shaped ground-based Michelson interferometers with opening angle δ = π/3 in a triangular configuration, i.e., rotated with respect to each other by an angle ω = 2π/3. Consequently, ET will allow one to perform crosscorrelation measurements in a network of three identical and co-located (but not co-aligned) detectors. The overlap reduction function for a pair of ET detectors can be found in ref. [199]. We extract the function graph from figure 8 of ref. [199] and rescale it by a factor sin 2 β = 3/4, in order to adjust its overall normalization, γ IJ → cos (2ω) = −1/2 in the small-frequency limit. 13 The normalized overlap reduction function thus obtained is shown in the right panel of figure 11. • DECIGO is envisioned to consist of two satellite-borne triangular Fabry-Pérot (FP) interferometers with opening angle δ = π/3 in a hexagonal configuration (i.e., ω = π). The overlap reduction function for this pair of detectors has been computed in ref. [200] (see also appendix B of ref. [201]). In our analysis, we shall use the numerical results obtained in ref. [200], which were kindly provided to us by Sachiko Kuroyanagi (see the right panel of figure 11). As expected, the normalized overlap reduction approaches unity in the small-frequency limit. Similarly as in the case of LISA, we assume that DECIGO will allow one to construct two independent data streams. At small frequencies, DECIGO's overlap reduction function therefore approaches the same value as LISA's signal response function, Γ IJ → 2/5 sin 2 β = 3/10. This value is a good approximation up to frequencies of O (100) Hz, which is due to the fact that DECIGO's relatively short arm length, L DECIGO arm = 1000 km, results in a large characteristic frequency, f DECIGO * = c light / 2πL DECIGO arm 47.71 Hz. The laser travel distance in DECIGO's FP cavity is by contrast much larger than L DECIGO arm , which renders DECIGO most sensitive in the deci-Hertz range; hence DECIGO's name. In many cases, it thus suffices to simply assume a constant overlap reduction function, Γ IJ 3/10.
• BBO is planned to have a similar geometry as DECIGO (i.e., two satellite-borne triangular interferometers with δ = π/3 and ω = π). It is, however, supposed to have a larger arm length, L BBO arm = 50 000 km, and utilize Michelson instead of FP interferometers. The overlap reduction function of the two BBO units has been calculated in ref. [85], and the corresponding numerical results are available from ref. [202]. We plot the normalized overlap reduction function in the right panel of figure 11. Following ref. [85], we assume a single data channel, such that γ IJ → 1 and Γ IJ → 1/5 sin 2 β = 3/20 in the small-frequency limit.

JHEP01(2021)097
Overlap reduction functions for pulsar pairs in a pulsar timing array. The overlap reduction function for two pulsars I and J in a PTA is of the form [203,204], (A. 36) with R I being the response function for a single pulsar (i.e., J = I) and ζ IJ denoting the Hellings-Downs factor [205] for two pulsars that are separated by an angle ψ in the sky, Here, we introduced c ψ = (1 − cos ψ) /2 and made sure that ζ IJ is properly normalized, i.e., ζ II = 1 for a single pulsar. Note that Γ IJ and R I are now no longer dimensionless as in the case of the interferometer experiments. Both quantities are expressed in units of Hz −2 , which is line with the fact that the timing noise of a pulsar has dimension Hz −3 (see section A.3).
For future experiments, the pulsar distribution in the sky is still unknown. In order to estimate their sensitivity, it is therefore reasonable to replace ζ IJ by an expectation value that reflects the expected distribution of pulsars in the PTA. A common choice, e.g., is to replace ζ IJ by its root-mean-square (RMS) expectation value assuming a flat prior on cos ψ, ζ rms = 1 2 This is the value that we shall employ in our analysis. In order to estimate the sensitivities of IPTA and SKA, we will thus work with the following overlap reduction function, Γ IJ (f ) = ζ rms 12π 2 f 2 . (A.40)

A.3 Detector noise power spectra
In addition to the transfer functions discussed in the previous section, we also require the intrinsic noise spectra D I noise for each detector. The noise spectra of aLIGO, aVirgo, KAGRA, and CE can be downloaded from LIGO's Document Control Center (DCC) [dcc.ligo.org]. In table 4, we list the DCC documents that give access to the respective data files. Note that, for aLIGO and aVirgo, we also consider sensitivities representative for O2 in addition to their envisioned design sensitivities. The noise spectrum of ET can be downloaded from ref. [206] (see also DCC document P1600143 [207]). In our analysis, we shall employ the ET-D-sum noise curve. Moreover, as regards the noise spectra of the future spacebased interferometers and PTA experiments, we will make use of the following analytical estimates:
LISA: following ref. [196], we write LISA's noise spectrum as a sum of two contributions, In our analysis, we shall assume the following values for DECIGO's characteristic experimental variables [90]: arm length L DECIGO arm = 1000 km, laser output power P = 10 W, laser wavelength λ = 532 nm, mirror mass M = 100 kg, and mirror radius R = 0.5 m. These values determine the finesse F of the FP cavity as well as the effective laser output power P eff , F = π (r E r F ) 1/2 1 − r E r F 10.18 , P eff = r E t 2 where the quantities r E , r F , and t F are given as follows (see ref. [200] for details), Here, 1/T is the cadence of the timing observations and σ t denotes the RMS error of the timing residuals. We shall assume a future IPTA data set based on N = 20, T = 2 weeks, and σ t = 100 ns, which goes beyond the timing precision achieved in the first IPTA data release [118]. As for SKA, we are even more optimistic, assuming an ambitious PTA with N = 50, T = 1 week, and σ t = 30 ns. These values are inspired by the assumptions made in refs. [121,208]. We thus obtain the following timing noise spectra for IPTA and SKA,

JHEP01(2021)097
In figure 12, we plot the noise spectra of all present and future interferometer experiments that we consider this paper. The timing noise of PTA experiments, which has dimension Hz −3 instead of Hz −1 and which is frequency-independent in the idealized case, is not shown.

A.4 Strain noise power spectra
Based on the transfer functions and detector noise spectra introduced in sections A.2 and A.3, we are now able to compute the (effective) strain noise power spectra of all autoand cross-correlation searches for a SGWB signal that we are interested in [see eqs. (A.14) and (A. 22)]. Here, when considering the effective strain noise of this network at design sensitivity, we can again set D aLHO noise = D aLLO noise = D aLIGO noise . However, when evaluating S eff noise for the representative O2 noise spectra (see , which we will only evaluate at design sensitivity, such that D aLHO noise = D aLLO noise = D aLIGO noise . For more information on the HLV and HLVK networks, see also refs. [209,210].

JHEP01(2021)097
In section 3.2, we use these effective strain noise spectra to draw PLISCs for NANOGrav, PPTA, and EPTA (see figure 1), assuming the following total observing times [108,112,114], Refs. [108,112,114] also directly present upper limits on the strength of a power-law SGWB signal as a function of the spectral index p. We explicitly check that these limits are consistent with the amplitudes Ω p in eq. (A.28), which we require for the construction of the PLISCs.
IPTA, SKA: we estimate the effective strain noise spectra of the future PTA experiments IPTA and SKA based on the transfer function in eq. (A.40) and the timing noise in eq. (A.49), where D PTA noise (f ) = 2 T σ 2 t and R PTA (f ) = 1/ 12π 2 f 2 and where we assumed the same Hellings-Downs factor for all pulsar pairs, ζ IJ = ζ rms [see eq. (A. 39)]. In section 3.2, we use the estimate in eq. (A.61) to draw the PLISCs for IPTA and SKA, assuming an effective observing time of t eff obs = 20 yr. Here, t eff obs accounts for the fact that the scaling behavior of the SNR with t obs changes for large values of t obs [216]. We define t eff obs such that the SNR in the weak-signal regime after an effective observing time t eff obs [see eq. (1.1)] equals the true SNR in the intermediate-or strong-signal regime after an actual observing time t obs . A more detailed discussion of the respective t obs values for IPTA and SKA is left for future work.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.