A fresh look at the gravitational-wave signal from cosmological phase transitions

Many models of physics beyond the Standard Model predict a strong first-order phase transition (SFOPT) in the early Universe that leads to observable gravitational waves (GWs). In this paper, we propose a novel method for presenting and comparing the GW signals that are predicted by different models. Our approach is based on the observation that the GW signal has an approximately model-independent spectral shape. This allows us to represent it solely in terms of a finite number of observables, that is, a set of peak amplitudes and peak frequencies. As an example, we consider the GW signal in the real-scalar-singlet extension of the Standard Model (xSM). We construct the signal region of the xSM in the space of observables and show how it will be probed by future space-borne interferometers. Our analysis results in sensitivity plots that are reminiscent of similar plots that are typically shown for dark-matter direct-detection experiments, but which are novel in the context of GWs from a SFOPT. These plots set the stage for a systematic model comparison, the exploration of underlying model-parameter dependencies, and the construction of distribution functions in the space of observables. In our plots, the experimental sensitivities of future searches for a stochastic GW signal are indicated by peak-integrated sensitivity curves. A detailed discussion of these curves, including fit functions, is contained in a companion paper [2002.04615]. The data and code that we used in our analysis can be downloaded from Zenodo [https://doi.org/10.5281/zenodo.3699415].


Introduction
One of the prime targets of gravitational-wave (GW) astronomy in the coming years is going to be the detection of relic GWs that were produced in the early Universe [3]. There is a wealth of possible sources of primordial GWs [4,5], among which cosmological phase transitions [6] represent a particularly well-motivated scenario. Indeed, many extensions of the Standard Model (SM) predict a strong first-order phase transition (SFOPT) in the early Universe that readily results in an observable GW signal [7]. This nourishes the hope that future GW observations will open up new avenues in the hunt for physics beyond the Standard Model (BSM) that are complementary to conventional experiments, such as particle colliders and experiments aiming at the direct or indirect detection of dark matter (DM).
In recent years, the GW phenomenology of a large number of BSM models has been investigated (see, e.g., the SFOPTs studied in Refs. ). This leads to the question in which way one may best present and compare the different predictions of these models. At the moment, there is no universal approach towards such a systematic model comparison. Most authors simply focus on the GW signal in a particular BSM model, while ignoring possible similarities or differences to the GW signal in other models. In view of the multitude of models studied in the literature, it is therefore difficult to pin down the characteristics of a given model and to precisely quantify how its GW phenomenology differs from the phenomenology of other models. In the present paper, we attempt to remedy this situation by proposing a novel method for presenting and comparing the GW signals in different models. Our approach is based on the observation that the GW signal from a cosmological phase transition can effectively be parametrized in terms of a set of three peak frequencies, {f b , f s , f t } and three peak amplitudes, {Ω peak b , Ω peak s , Ω peak t }, which respectively correspond to the three physical processes that source GWs during a SFOPT (see Sec. 2). Such a parametrization is feasible because the frequency dependence of the GW spectrum close to the peak frequencies is approximately model-independent. Based on this observation, one is thus able to construct scatter plots in the space of peak frequencies and peak amplitudes that provide a compact and comprehensive overview of the GW phenomenology of a particular model.
In this paper, we will illustrate our idea by constructing the signal region for a particularly simple BSM model, namely, the real-scalar-singlet extension of the SM (xSM) [46][47][48][49][50][51][52], which supplements the SM Higgs sector by an additional real gauge-singlet scalar, S, that also obtains a nonzero vacuum expectation value (VEV) during the electroweak phase transition [53][54][55][56][57][58][59]. The GW phenomenology of the xSM has been studied before [12,19,27,43,[60][61][62][63][64]; however, in our analysis, we are going to look at the GW signal in this model from a new perspective by projecting the results of our parameter scan into the space of peak frequencies and peak amplitudes. This results in sensitivity plots for future space-based GW experiments that are reminiscent of similar plots in the field of DM direct-detection experiments. There, it is standard practice to compare the sensitivity of ongoing and upcoming experiments with the signal regions of different models in the parameter space spanned by the DM mass and the (spin-dependent or spin-independent) DM-nucleon scattering cross-section. In this sense, the goal of this paper is to construct equivalent sensitivity plots in the case of GWs from a cosmological phase transition in the early Universe.
In our new sensitivity plots, we can no longer use the standard power-law-integrated sensitivity curves [65] that are typically shown in the literature. These curves are useful to illustrate experimental sensitivities directly in terms of the GW energy spectrum, Ω GW (f ), as a function of GW frequency, f . However, strictly speaking, they only apply to signals that are described by a pure power law. Moreover, our scatter plots do not show the full GW spectrum as a function of frequency by construction. For these reasons, the experimental sensitivities in our plots are indicated by what we will refer to as peak-integrated sensitivity curves (PISCs). We will define these curves and briefly discuss their properties in Sec. 2. A more detailed discussion of the concept of PISCs can be found in a companion paper [1].
Our PISC plots offer several advantages compared to conventional presentations of the GW signal from a cosmological phase transition. First of all, every model results in a distinct signal region in our plots. This sets the stage for a systematic model comparison based on the size, shape, parameter dependence, etc. of different signal regions. In addition, our approach allows one to project the dependence on underlying model parameters directly into the space of physical observables, i.e., into the space of peak frequencies and peak amplitudes. Our plots thus facilitate the exploration of the GW phenomenology in a given model, and in particular, the exploration of underlying parameter dependencies. The same is also true in the case of more sophisticated analysis, such as, e.g., global fits resulting in likelihood functions. Finally, one can use our plots to construct distribution functions (i.e., histograms) of peak frequencies and peak amplitudes by projecting the signal regions in our plots onto the x-and y-axes, respectively. These distribution functions also characterize the GW phenomenology of a given model and are helpful in the comparison and exploration of different models.
The rest of the paper is organized as follows. In Sec. 2, we will first summarize the different contributions to the GW signal from a SFOPT and introduce the concept of PISCs. In Sec. 3, we will then discuss in more detail the efficiencies of the different GW production mechanisms that are at work during a SFOPT. This will allow us to estimate more precisely the fractions of the latent heat that are respectively converted into gradient energy of the scalar field, kinetic energy of the thermal plasma, etc. during the phase transition. Next, we will review the xSM and outline all relevant theoretical and experimental constraints on its parameter space in Sec. 4, before finally presenting our main results in Sec. 5. In this section, we will show our new sensitivity plots, discuss the dependence of the GW signal on some of the underlying model parameters, and construct histograms of possible peak frequencies and peak amplitudes. In Sec. 6, we will compare our novel method with existing approaches of studying the GW signal from SFOPTs, thus summarizing the key features and advantages of our idea. Sec. 7 contains our conclusions and a brief outlook on possible future steps. In Appendix A, we collect the results of a partial-wave analysis that allow us to constrain the parameter space of the xSM based on the requirement of perturbative unitarity.

Peak-integrated sensitivity curves
SFOPTs give rise to three independent sources of stochastic GWs, namely, due to collisions of scalar-field bubbles (b), sound waves in the bulk plasma (s), and vortical motion in the latter, i.e., magnetohydrodynamic turbulence (t). The respective GW spectra can be approximately obtained from numerical and semianalytical calculations and are quantified in terms of a peak amplitude, Ω peak i , and a spectral shape, S i , which depends on the peak frequency, f i , (2.1) Here, i ∈ {b, s, t} and α, β/H * , T * , v w , κ i are quantities characterizing the phase transition that are closely related to the hydrodynamics of the process. Therefore, we will define them in Sec. 3 where κ b and κ v indicate the efficiencies of converting the latent heat released during the SFOPT into kinetic energy of the expanding scalar-field bubbles and the surrounding plasma, respectively. For the energy going into the surrounding plasma, ε determines the fraction going into turbulent bulk kinetic energy. The peak frequencies read [66] f b = 1.6 · 10 −5 Hz g * (T * ) 100 1 6 T * 100 GeV (2.3a) f s = 1.9 · 10 −5 Hz g * (T * ) 100 1 6 T * 100 GeV They enter the spectral shape functions, S i , which are given by [66] S 11 3 . (2.4c) The frequency h * corresponds to the wave number k * that equals the Hubble rate at the time of GW production redshifted by the expansion of the Universe up to the present time, These contributions constitute a signal, which needs to be compared to the noise spectrum of the experiment under consideration to obtain the signal-to-noise ratio (SNR) [3,68] Here, n det distinguishes between experiments that aim at detecting the stochastic GW background by means of an auto-correlation (n det = 1) or a cross-correlation (n det = 2) measurement. The n det = 1 case refers to situations where an experiment consists of just one detector, while the n det = 2 case refers to situations where it is possible to cross-correlate the two signals of a detector pair within a detector network. Further below, we will specifically consider three satellite-borne GW interferometers: the Laser Interferometer Space Antenna (LISA) [69,70], the Deci-Hertz Interferometer Gravitational-Wave Observatory (DECIGO) [71][72][73][74], and the Big-Bang Observer (BBO) [75][76][77]. Given the envisaged configuration of these experiments, we will set n det = 1 for LISA and n det = 2 for DECIGO and BBO. For a discussion of the noise spectra of these experiments, we refer the reader to the Appendix of Ref. [1].
Having specified an experiment and its noise spectrum, Eq. (2.6) can also be written as (2.7) Here, the integration over the frequency range has already been carried out implicitly, where we included a conventional factor of 2 for i = j that results from the square in Eq. (2.6) and i, j ∈ {b, s, t}. The mixed peak amplitudes are defined as the respective geometric means, The peak-integrated sensitivities, h 2 Ω i/j PIS , are functions of (at least one of) the peak frequencies. While Eq. (2.7) looks more complicated than Eq. (2.6) at first sight, it conveys an important message: For a given experiment and observation time, the SNR is uniquely determined by the peak energy densities and the corresponding peak frequencies, once the integrals in Eq. (2.8) have been carried out. These peak quantities only depend on the modelspecific SFOPT parameters, not on the GW frequency itself. Therefore, we can visualize the sensitivity, i.e., a certain SNR threshold, by drawing the corresponding contour lines in the f i -h 2 Ω peak i/j planes; these are the anticipated peak-integrated sensitivity curves (PISCs). For a parameter point of a given model, we can compute the quantities α, β/H * , T * , κ i , and hence also the peak quantities in Eqs. (2.2, 2.3). Thus, each parameter point of the model, corresponding to an entire GW spectrum, is reduced to a single point in the six f i -h 2 Ω peak i/j planes. This opens up possibilities for comprehensive parameter scans and intuitive model analysis or comparison. Any point above a single PISC will immediately yield a SNR above the chosen threshold, while points below it can only surpass the threshold if the sum of contributions is larger than the threshold. We will illustrate this procedure and its applications in more detail in Sec. 5, where we study the GW phenomenology of the xSM.

Efficiencies of the different GW production mechanisms
In order to obtain the amount of energy that is transformed into GWs during a SFOPT, the dynamics of the expanding scalar-field bubbles has to be analyzed. The velocity profile of these bubbles, v(ξ), satisfies the differential equation [27,67,78,79] 2 which is a function of ξ ≡ r/t due to the symmetry of the problem, with r being the distance from the center of the bubble and t the time since nucleation, and is the Lorentz boost factor from the bubble wall rest frame to the bubble center rest frame. In our analysis, we assume the ultrarelativistic value c 2 s = 1/3 for the speed of sound in the plasma. Depending on the boundary conditions, this equation has three qualitatively different types of solutions, namely, detonations, deflagrations, and hybrid solutions. For detonations, the bubble wall moves at a supersonic speed and hits the fluid, which is at rest in front of the wall. This type of transition leads to the strongest GW signals. In deflagrations, the wall is subsonic, and the fluid behind the wall is at rest. Deflagration-type transitions are important for scenarios of electroweak baryogenesis, where out-of-equilibrium dynamics in front of the wall are needed. The hybrid solutions are a combination of the two, i.e., supersonic deflagrations, where the fluid is at rest neither in front nor behind the wall.
On very general grounds, one can determine the energy stored in bubble walls, which is subsequently released into GWs via bubble collisions. In order to do so, one compares the driving force of the bubble dynamics, i.e., the difference in pressure across the wall due to the difference in potential energy, p 0 = ∆V , with the friction force that counteracts this acceleration, p fric = −∆P LO − γ ∆P NLO , which has a leading order (LO) and a next-toleading order (NLO) contribution [14]. Here, the appearance of the relativistic γ factor in the NLO contribution has important consequences, as it effectively limits the energy stored in the bubble wall (by limiting its speed and making a "runaway" scenario hard to realize). This limits / suppresses the strength of the GW signal from bubble wall collisions. For the efficiency factor, it follows that [67] In this equation, γ eq ≡ α−α∞ αeq is the relativistic factor reached when the LO and NLO frictional forces exerted on the bubble wall by the plasma equilibrate with the driving pressure, and the wall velocity reaches its final value. The quantities α ∞ and α eq are defined analogously to α as the amount of energy density relative to the radiation energy density (see Ref. [67]), where ∆m 2 ≡ i c i N i ∆m 2 i is the mass difference across the phase transition, weighted by internal number of degrees of freedom (DOF), N i , and a bosonic (c i = 1), or fermionic (c i = 1/2) factor. The NLO term is weighted by the gauge coupling, g 2 ∆m V ≡ i g 2 i N i ∆m i . For the SM and its singlet extension of interest to us, this yields [66,67,79] where φ * is the Higgs field value at the time of nucleation giving mass to the particles across the phase transition in our case. Conversely, γ * corresponds to the velocity that would be reached if only the LO friction term ∆P LO were present. Since this LO term is independent of γ, it can grow to larger values [67], is the initial radius at nucleation and R * = 3 √ 8π v w /β; see Ref. [14]. 2 For a more detailed discussion of these relations, we refer to Ref. [67].
With all formulas at hand, we can now quickly come back to the influence of the γ factor accompanying the NLO friction term mentioned before. Its appearance leads to the difference between γ eq and γ * and thereby to a suppression of the GW signal from bubble collisions by a factor γ eq /γ * ; see Eq. (3.2). The remaining energy released during the phase transition, i.e., a fraction 1−κ b , is converted into heat as well as kinetic energy of the plasma. While the thermal energy does not result in the production of GWs, the kinetic energy in the plasma will excite sound waves as well as turbulence, both of which act as GW sources. With the bubble velocity profile at our disposal, we are able to determine the amount of energy that is transformed into kinetic energy of the plasma [79] κ Figure 1: Distribution of efficiencies in the xSM for converting the energy released during the SFOPT into collisional energy of the bubbles (b, red), sound waves (s, yellow), and turbulence (t, blue). The normalization of the vertical axis is arbitrary. Note that, here and in the rest of our analysis, we have fixed the bubble wall velocity at v w = 0.9 in order to increase the strength of the GW signal. As a consequence, most of the parameter points in our scan describe phase transitions of the detonation type. For lower velocities, which typically result in sub-or supersonic deflagrations and which are relevant in the context of electroweak baryogenesis [43], the GW signal from sound waves can be significantly suppressed [80].
where w(ξ) denotes the enthalpy density and is given by For relativistic velocities, v w ∼ 1, the expression in Eq. (3.6) is well approximated by [79] κ In our numerical study, we will use the full approximations in Ref. [79] for all regimes of v w . Regarding the amount of bulk kinetic energy that is respectively transformed into turbulence and sound waves, there is currently no consensus in the literature. Many authors quote ε = 5 . . . 10 % of bulk kinetic energy in the form of turbulence (see, e.g., Refs. [27,66]), whereas others use an equal split between turbulence and sound waves (see, e.g., Ref. [81]). However, more recently, the authors of Ref. [67] found a significant increase in the turbulencesourced GW spectrum based on the analysis in Ref. [82]. We will follow Ref. [67] and compare the duration of the period of sound waves τ s with the inverse Hubble rate at nucleation, where the mean plasma velocity U f is given by [67,83] In conclusion, we can either use Eq. (3.9) in Eq. (2.2), or define where an appropriate power needs to be included. We emphasize that the estimate of Ref. [67] yields an upper limit for the turbulence contribution to the GW signal and could overestimate turbulence with respect to sound waves. However, further studies are needed to settle this matter, and here we use the upper limit to show that the turbulence contribution can be very important and even dominate the signal. In Fig. 1, we show the distribution of efficiencies for the xSM from our numerical analysis. This shows that, even for the high wall velocity close to the speed of light, v w = 0.9, collisions are negligible, while the sound waves and, most importantly, even the turbulence contributions can dominate the total GW signal.

Real-scalar-singlet extension of the Standard Model
Let us now consider the xSM as a concrete and simple example of a BSM model that results in a GW signal from a SFOPT. The tree-level scalar potential of the xSM reads where H is the SM Higgs doublet and S a real scalar singlet. Note that we also allow for odd powers of the field S in the scalar potential. That is, we do not impose an additional Z 2 symmetry on the scalar sector of the xSM, as it is sometimes done in the literature.
To make the parameters more easily accessible, we recast the Lagrangian parameters into low-energy observables. In a first step, we expand the fields around their electroweak VEVs by using Demanding that there is a minimum at (h 0 , s) = (0, 0) gives the conditions In a second step, we determine the mass eigenstates h 1 and h 2 by diagonalizing the mass matrix of the neutral scalars h 0 and s. We introduce the mixing angle θ and define h 1 to be the SM Higgs boson with m h = 125 GeV, such that h 2 has the mass m s . This results in where we made use of the shorthand notations c x ≡ cos (x) and s x ≡ sin (x). Since m h and v h are fixed by measurements, the free parameters in this parametrization are (v s , m s , θ, µ 3 , λ S ). However, also these parameters are (directly or implicitly) constrained by theoretical arguments and experimental measurements, namely: 1. Boundedness of the scalar potential from below : 2. Perturbative unitarity: The leading partial-wave amplitudes for a given 2 → 2 scattering process need to obey max {A i } ≤ 8π, where A i are the absolute values of the eigenvalues of the S matrix (which we spell out for completeness in Appendix A).
3. Vacuum stability: In order to judge whether the electroweak vacuum, which is a local minimum by construction, 3 is also the global minimum of the potential, one needs to study the remaining vacua either numerically or analytically (see Ref. [56] for details).

Measurement of Higgs coupling strengths:
Given that the neutral component h 0 of the scalar doublet is in general not identical to the light scalar h 1 , it is possible to derive constraints by noting that h 1 couples to SM particles with a reduced strength ∝ cos 2 (θ). It is thus found that |sin(θ)| 0.3 is excluded at the 95 % C. L. [48,84,85] 5. Electroweak precision tests: According to Ref. [27], measurements of the W -boson mass result in the strongest constraints on the xSM parameter space. As a rule of thumb, the W -boson mass is most constraining for scalar masses m s 300 GeV, where it translates into sin θ 0.2 for v s = 0.1 v h . However, new physics beyond the xSM can always be used to relax this constraint [85].
The dynamics of the phase transition are governed by its temperature-dependent effective potential. Since temperature effects appear at one-loop order, the corresponding zero-temperature corrections have also to be taken into account at the same order. The effective potential up to one-loop order therefore reads The tree-level potential is given by Eq. (4.1). The second term, V 0 1 , denotes the Coleman-Weinberg potential, where the upper (lower) signs correspond to bosons (fermions). In addition, there are two more corrections that we take into account. First, we work with the thermally enhanced or improved finite-temperature potential, which is obtained by adding to the field-dependent masses in Eqs. (4.9, 4.10) the leading thermal corrections (see Ref. [86] for a recent discussion), where the thermal masses c i T 2 can be found, e.g., in Ref. [56]. Second, we keep the scalar VEVs and masses at their tree-level values by including the following finite counter terms, where the coefficients are chosen so as to satisfy the following renormalization conditions, In our numerical analysis, we employ the code CosmoTransitions [87] to compute the nucleation temperature, tunneling action, and all other phase-transition-related quantities that are key to calculating α and β -the parameters necessary for determining the GW spectrum as described in Secs. 2 and 3. Furthermore, we do not take into account possible effects in small regions of parameter space that might arise in situations with very strong supercooling (potentially leading to an additional period of vacuum domination; see Ref. [14]).

New sensitivity plots, parameter dependencies, histograms
In order to appreciate the capabilities and advantages of the procedure outlined in Sec. 2, we will now demonstrate how to use Eq. (2.7) in practice. To this end, we will consider a characteristic set of parameter points sampled within the following ranges, We sample all parameters making use of a linear prior, except for m s , for which we use a logarithmic prior, and only accept points that fulfill the theoretical consistency constraints (perturbative unitarity and vacuum stability). For each point in parameter space thus obtained, we compute the phase transition parameters α, β/H * , T * , and κ i . In this way, we generate a data set consisting of roughly 6000 parameter points, all of which successfully result in a SFOPT. For illustrative purposes, we fix the wall velocity at v w = 0.9, the SNR  , and each point represents an entire GW spectrum of a particular physical origin. The colorful curves are the peak-integrated sensitivity curves that are at the heart of our approach, and the colorful bands indicate that some of the PISCs depend on more than just one frequency (in this case, these frequencies are varied according to the spread in the data; see insets). A point above any one of the PISCs / bands has ρ > 1 (black), while points in bands or below the PISCs must be checked individually. For dark gray points, the combined SNR is above the LISA threshold, while light gray points are not detectable by LISA. Dashed lines / lighter bands indicate the projected BBO and DECIGO sensitivities.
threshold at ρ thr = 1, and the observation time at t obs = 1 yr. It is straightforward to generalize our analysis to other values of ρ thr and t obs by rescaling all PISCs by (t obs /yr) 1/2 /ρ thr ; see Eq. (2.7). Considering each contribution in Eq. (2.7) separately, we draw each PISC as a function of the corresponding peak frequency (varying all other frequencies within the ranges spanned by our data set). Next, accompanying each PISC Ω i/j PIS , we compute the peak frequencies f i (see Eq. (2.3)) and peak amplitudes Ω peak i/j (see Eq. (2.2)) for each point in our data set and scatter these points in our six PISC plots. In this way, we obtain Fig. 2.
Let us now highlight the important features that can be extracted from Fig. 2. Each point in one of the six panels represents a choice of parameter values in the xSM. That is, conventionally, one would draw an entire GW spectrum for each such point, which one would then have to compare to the standard power-law-integrated sensitivity curves (see Ref. [1] for a more detailed discussion). In our approach, by contrast, we simply need to verify whether a given point is above any of the six PISCs (of one experiment). In that case, the SNR will automatically be larger than the predefined threshold (black points in Fig. 2), indicating that this parameter point will be probed by LISA (or BBO, or DECIGO). In the opposite case, the point may still surpass the SNR threshold, namely, if the sum of contributions is larger than the threshold. We indicate such points by a dark gray color. Finally, points that are not testable by LISA are shown in light gray. This procedure allows us to map the phenomenologically relevant parameter space into the space of GW observables. Note also that, for any contribution that depends on more than one peak frequency (i.e., the turbulence (t), sound wave / turbulence (s/t), and bubble collisions / turbulence (b/t) channels), we cannot draw a single PISC in a two-dimensional plot. As a consequence, we draw peak-integrated sensitivity bands in these plots, where the peak frequency that is not shown is varied in a range that can either be chosen freely or according to the respective spread found in the data. 5 If we had not fixed the wall velocity at v w = 0.9, we would also have to include a band in the bubble collisions / sound waves (b/s) channel; see Eq. (2.3).
As an illustration of the usefulness of the PISC approach, we next identify the most constraining PISC channel, which turns out to be the s/t channel in our case. In this channel, only 45 of the 390 points with ρ > 1 require support from the other channels to surpass the SNR threshold. For this particular channel, we visualize the distribution of certain parameters in the plane of GW observables in order to answer the question which regions of the model parameter space will be probed by future GW missions; see Fig. 3. To do so, we require all points in our data set to satisfy the theoretical constraints 1. -3. in Sec. 4; however, at the same time, we do not reject points in violation of the experimental constraints 4. -5. to highlight the complementarity of collider and GW probes. As a simple example, we show in Fig. 3a the s/t-PISC together with the xSM parameter points whose color indicate the numerically computed nucleation temperature, with their size being proportional to the value of the latent heat parameter, α. We thereby verify a couple of simple statements, namely, that lower nucleation temperatures (i.e., stronger supercooling) correlate with larger α values and thus a stronger phase transition, which in turn gives a louder GW signal. This is also illustrated by the accompanying histograms, which show the distributions of model points with T n > 100 GeV (green) and those with T n < 100 GeV (orange). This confirms our expectations and serves as a cross-check of our results. We find that most scenarios with extreme values of α can be detected (or ruled out) with LISA's design sensitivity.
Next, we can test the model for less obvious correlations as shown in Fig. 3b, where in the same plot the color now indicates the value of the trilinear portal coupling µ HS . While there is no information in the sign of µ HS , we do find that large trilinear couplings preferably occur for low frequencies and thereby induce strong GW signals making them accessible to the planned space-based GW missions. This impression is supported by the distribution of µ HS values according to the histograms, where blue indicates values of the trilinear coupling above 5 TeV, green 500 GeV < |µ HS | < 5 TeV, and orange |µ HS | < 500 GeV. These findings are consistent with the results in Ref. [52], where it has been found that larger trilinear couplings tend to lead to stronger phase transitions. The size of the points in Fig. 3b is determined by the value of the dimensionless portal coupling, λ HS , which shows only a mild preference for larger values in the case of strong GW signals.
To understand the influence of the physical mass of the second scalar after symmetry breaking, m s , we consult Fig. 3c, where this mass is represented by the color code. Again, we find a notable correlation; however, this time, lighter scalars typically induce stronger GW signals, as one can verify from the histograms, which show in red the distribution of scalar masses m s < 200 GeV and in blue the distribution of scalar masses m s > 200 GeV. This time, the size of the points is related to the scalar-singlet self-coupling, λ S , whose values Panel (a) shows the variation of the nucleation temperature T n (color) and the latent heat α (size); in the histograms, we distinguish T n > 100 GeV (green) and T n < 100 GeV (orange). Panel (b) shows the variation of the trilinear coupling µ HS (color) and the portal coupling λ HS (size); in the histograms, we distinguish |µ HS | > 5 TeV (blue), 500 GeV < |µ HS | < 5 TeV (green), and |µ HS | < 500 GeV (orange). Panel (c) shows the variation of the scalar mass m s (color) and the self-coupling λ S (size); in the histograms, we distinguish m s > 200 GeV (blue) and m s < 200 GeV (red). Panel (d) shows the variation of the mixing angle θ (color) and the scalar mass m s (size); in the histograms, we distinguish θ > 0 (blue) and θ < 0 (red). appear evenly distributed over the space of GW observables. Finally, Fig. 3d shows the distribution of mixing angles in the plane of GW observables. We observe that large mixing angles occur for all frequencies, however, mostly in a narrow band centered in the scatter plot. In Fig. 3d, the colored histograms show the distributions of points with θ > 0 (blue) and θ < 0 (red), which are evenly distributed across the GW observables. Recalling that |sin θ| > 0.3 is excluded, we see that collider experiments are, in fact, complementary to GW searches in the sense that they do not probe the same regions of parameter space.

Comparison with existing approaches
In the previous sections, we have introduced the novel concept of PISC plots, choosing the xSM as a concrete example to illustrate our idea. Based on the results in Fig. 2 and 3, we are therefore now able to compare our new idea to other approaches in the literature for presenting the sensitivity of future experiments to the GW signal from a SFOPT (for a more comprehensive discussion, see Ref. [1]). A commonly employed strategy, e.g., is to plot the full GW spectrum together with the standard power-law-integrated sensitivity curves, which were introduced by Romano and Thrane in Ref. [65]. This approach conveys a useful impression of a given experiment's sensitivity reach, but comes with a number of limitations. First of all, it requires one to draw an individual GW spectrum as a function of frequency for each parameter point of interest. This quickly becomes impracticable, resulting in very busy plots as soon as one intends to study O (10) or more points in the model parameter space. A possible way out of this problem would be to restrict oneself to representing each spectrum by merely a single point: a point indicating the peak amplitude at the peak frequency. In fact, such a strategy would result in plots similar to our PISC plots. However, the important difference in this case would be that a scatter plot of peak amplitudes and peak frequencies in combination with power-law-integrated sensitivity curves would no longer contain any information on the expected SNR. To see this, recall that power-law-integrated sensitivity curves only have a well-defined statistical meaning for GW signals that are described by a pure power law (hence the name). For a GW signal from a SFOPT, this assumption is maximally violated close to the most relevant part of the spectrum, namely, the peak in the spectrum, where the frequency dependence changes from a positive power law to a negative power law. A main motivation behind our PISC approach therefore is to remedy this shortcoming. Our PISC plots also feature observables such as GW frequencies and signal strengths on the axes; but in contrast to the standard power-law-integrated sensitivity curves, our PISCs are constructed such that they still retain the full information on the SNR. For a given point in a PISC plot, the (partial) SNR simply corresponds to the vertical separation between the point and the PISC of interest. We therefore argue that our PISCs are the better choice compared to the standard power-law-integrated sensitivity curves for this type of signal. As long as the shape of the signal is not precisely known, it is reasonable to stick to power-law-integrated sensitivity curves. However, as soon as more information on the signal shape is available, which is the case for the GW signal from a SFOPT, one should also make use of this extra information and account for it in the construction of the sensitivity curves.
A second approach often employed in the literature is to present plots of the SNR as a function of some of the underlying model parameters. We reproduce plots of this type in Fig. 4, where we show projections of our xSM data set onto the α -β/H * , α -T n , and β/H * -T n planes in combination with a color code for the expected SNR for LISA. Let us SNR ≥ 10 1 ≤ SNR < 10 SNR < 1 Figure 4: Projections of our full set of xSM parameter points onto the α -β/H * (top), α -T n (middle), and β/H * -T n (bottom) planes. The color code indicates the total SNR ρ for LISA. Note that, in none of the three plots, ρ is a smooth function of the respective parameters on the x-and y-axes. Fig. 4 in Ref. [27] shows a similar plot of the α -β/H * plane. now compare Fig. 4 to our PISC plots in Fig. 2 and 3. In doing so, we shall summarize the characteristic features of our PISC plots and point out the advantages of our new approach: 1. Our PISC plots retain the full information on the SNR and encode it on the y-axis.
A parameter point being separated from a PISC by factor ∆y along the y-axis simply corresponds to a partial SNR of ∆y. The total SNR for this point then follows from adding all partial SNRs in quadrature; see Eq. (2.7). This is particularly useful when one is interested in comparing different SNR thresholds, ρ thr , to each other. We also point out that additional color coding as in the three plots in Fig. 4 is not necessary to indicate the expected SNR. Instead, color coding can be used to include additional useful information; cf. Fig. 3. As an alternative to using a color code, it is also possible to present contour plots of the SNR as a function of α and β/H * , α and T n , etc. However, in this case, one is no longer able to work with projections onto twodimensional parameter subspaces. As evident from Fig. 4, such projections typically do not result in a smooth functional behavior of the SNR. Therefore, instead of projections, one has to restrict oneself to two-dimensional hypersurfaces (or slices) through the higher-dimensional parameter space, such that the SNR becomes a well-defined function with well-defined contour lines. The number of possible hypersurfaces that one may want to look at is arbitrarily large, especially, when one is interested in studying the sensitivity of an experiment without considering a particular particle physics model. Fig. 2 and 3 only depend on the experimental noise spectra and spectral shape functions in Eq. (2.4). In this sense, they represent truly experimental quantities that are insensitive to uncertainties on the theory side. This is not the case for SNR plots, where the information on the expected SNR is subject to all uncertainties entering the calculation, both on the experimental and theoretical side. For example, in Fig. 4, the distribution of the blue, green, and red points depends on how we compute the parameters α, β/H * , and T n , whereas the PISCs in Fig. 2 and 3 do not depend on this step in the analysis. Furthermore, our PISC plots -indicating the projected sensitivities of LISA, DECIGO, and BBO in terms of observables that are experimentally accessible, namely, peak frequencies and peak amplitudes -may be regarded more useful from an experimental perspective, as they are based on quantities that will likely play an important role in the experimental data analysis. SNR plots such as the plots in Fig. 4, on the other hand, are very useful from a model-builder's perspective, as they immediately illustrate a handful of important physical relations.

The PISCs in
3. It is straightforward to generalize our PISCs to other signal shapes. That is, for any signal that comes with a universal shape function S, one may simply repeat our procedure in Sec. 2 and construct an analogously defined shape-integrated sensitivity curve. Of course, this will only work up to some level of generality. The GW signal from inflation, e.g., can be described by a large range of different shapes, depending on the underlying model. In this case, it is not possible to construct a universally applicable sensitivity curve. The same is true for the GW spectrum from cosmic strings, which can not be fit by a universal shape function S. This being said, it should still be possible to construct meaningful sensitivity curves (for GWs from inflation, cosmic strings, etc.) if one is willing to restrict oneself to a more model-dependent analysis, such that the GW spectrum is well described by a specific template function after all. In addition, we point out that our approach is very flexible in the sense that it can be easily updated if our understanding of the shape functions S b , S s , and S t should improve in the future.

4.
A key idea behind our PISC approach is to decompose the total SNR into six partial SNRs, which respectively represent the three physical contributions to the GW signal as well as their three cross-correlations; see Eq. (2.7). Consequently, we end up with six different PISCs that we need to draw for each experiment. This is a helpful feature of our approach that allows for an easy comparison of the six different signal channels (s, b, t, s/b, s/t, b/t) that potentially contribute to the total signal. Not only do our plots illustrate the relative importance of these six channels, they also allow one to ask questions such as, e.g.: What happens if one completely ignores the contributions from scalar-field bubbles and turbulence to the signal? In the case of SNR plots, such a question would prompt one to redo the entire analysis, now focusing on the sound wave contribution to the signal only. In our PISC plots, on the other hand, the answer is trivial. All one would have to do would be to discard all but the s-PISC plot.

5.
To generate SNR plots such as those in Fig. 4 Based on these five points, we argue that our PISC plots are particularly well suited to illustrate the sensitivity of future experiments to the GW signal from a SFOPT. They especially allow for an easy comparison of the sensitivities of different experiments and provide a novel way of visualizing GW sensitivities that is reminiscent of plots that one often encounters in other fields of experimental physics, such as the standard sensitivity plots for DM direct-detection experiments. We believe that our PISCs have the potential to develop into a comparable standard with regards to the GW signal from cosmological phase transitions.

Conclusions and outlook
In this paper, we proposed a novel method for visualizing and exploring the GW phenomenology of BSM models that result in a SFOPT in the early Universe. Our approach is based on the observation that the spectral shape of the GW signal from a cosmological phase transition is approximately model-independent. Hence, it is feasible to encode the entire relevant and model-specific information in a set of characteristic observables, namely, three peak frequencies and three peak amplitudes. For a particular model of interest, one is thus able to construct scatter plots directly in the space of observables, which define the signal region of the respective model, in analogy to similar plots in other fields of experimental high-energy physics. In these scatter plots, each individual GW spectrum is represented by a set of points, which enables one to compare and explore the characteristics of nearly arbitrarily many individual spectra at the same time. This needs to be compared to the traditional approach, according to which one would simply plot all GW spectra as functions of the GW frequency and which is clearly limited in terms of the number of spectra that one may study at once.
We demonstrated our new procedure by means of a simple example, i.e., the GW signal from the electroweak phase transition in the real-scalar-singlet extension of the Standard Model (xSM). The main results of our analysis are shown in Figs. 2 and 3. These plots represent a projection of a sample of viable xSM parameter points into the space of peak frequencies and peak amplitudes. The experimental sensitivities of upcoming GW observatories are illustrated by what we call peak-integrated sensitivity curves (PISCs) in these plots; see Eq. (2.8). Specifically, we presented the PISCs of three future space-based GW interferometers: LISA, DECIGO, and BBO. A more detailed discussion of these sensitivity curves can be found in the companion paper [1]. The main message from our PISC plots in Figs. 2 and 3 is that they provide a bird's eye view of the GW phenomenology of the xSM. Not only do our PISC plots allow us to assess which parts of the xSM signal region will be probed by LISA, DECIGO, and BBO, respectively, they also set the stage for the study of underlying model-parameter dependencies as well as for the construction of histograms that reflect the relative rate of occurrence of particular peak frequencies and peak amplitudes.
In view of the large number of studies in the literature that discuss the GW signal from a SFOPT, we hope that our approach will open up new avenues for comparing the signals predicted by different models. In the next step, it would be crucial to repeat our analysis for as many BSM models as possible, so as to create a solid database for the systematic and quantitative comparison of different scenarios. This task is beyond the scope of the present paper, in which we merely aimed at outlining our basic idea; but we hope that the example analysis in this paper will stimulate further community efforts in this promising direction. Similarly, it will be worthwhile to continue the exploration of the GW signal in the xSM. In Ref. [88], e.g., we showed that the combination of the xSM with the type-I seesaw extension of the Standard Model has interesting implications for neutrino and Higgs physics, including the possibility to generate the baryon of the asymmetry of the Universe at energies in the TeV range. It would thus be interesting to use our PISC approach to study the complementarity of GW searches, collider searches, and the phenomenology of low-scale leptogenesis in the xSM. Alternatively, one could combine our analysis with a global fit of the xSM supplemented by a suitable (fermionic or bosonic) DM candidate. In this case, one could perform a profile likelihood analysis in the space of DM parameters, e.g., along the lines of Refs. [89,90], and project the resulting likelihood function into our PISC plots. Finally, it is important to note that it would be straightforward to generalize the PISC approach to other types of stochastic and cosmological GW signals whose spectral shape is also approximately model-independent. We, however, leave this and all other open tasks mentioned above for future work. Instead, we conclude by stressing that our approach bears the potential to develop into a useful new standard tool for studying GWs from cosmological phase transitions in the early Universe.
The amplitudes for charged initial and final states read and S 2 = −2λ H . Notice that the latter two simply imply λ H < 4π and |λ SH | < 8π.