Fingerprinting models of first-order phase transitions by the synergy between collider and gravitational-wave experiments

We investigate the sensitivity of future space-based interferometers such as LISA and DECIGO to the parameters of new particle physics models which drive a first-order phase transition in the early Universe. We first perform a Fisher matrix analysis on the quantities characterizing the gravitational wave spectrum resulting from the phase transition, such as the peak frequency and amplitude. We next perform a Fisher analysis for the quantities which determine the properties of the phase transition, such as the latent heat and the time dependence of the bubble nucleation rate. Since these quantities are determined by the model parameters of the new physics, we can estimate the expected sensitivities to such parameters. We illustrate this point by taking three new physics models for example: (1) models with additional isospin singlet scalars (2) a model with an extra real Higgs singlet, and (3) a classically conformal $B-L$ model. We find that future gravitational wave observations play complementary roles to future collider experiments in pinning down the parameters of new physics models driving a first-order phase transition.


I. INTRODUCTION
The discovery of the Higgs boson h at the CERN Large Hadron Collider (LHC) is one of the most prominent scientific developments in the past decades [1,2], establishing the spontaneous symmetry breaking and mass generation mechanism experimentally. Nevertheless, the whole picture of the Higgs sector remains unclear. Namely, the type and the number of Higgs multiplets, the shape of the Higgs potential, and the dynamics of the electroweak phase transition are all unknown. Understanding the nature of the Higgs sector is important not only for establishing the mechanism for the origin of mass but also for unraveling its connection to physics beyond the standard model (SM), such as neutrino oscillations, the existence of dark matter, baryon asymmetry of the Universe, and cosmic inflation. For example, electroweak baryogenesis in the early Universe [3] is an excellent physics case in which the Higgs sector leads us to new physics.
The conventional way to explore new physics models is to discover new particles and/or measure deviations from the SM predictions at collider experiments. So far, no new particle other than the Higgs boson has been found at the LHC. As for the deviations in various Higgs boson couplings, the expected accuracy is of a few percent level at the High-Luminosity LHC (HL-LHC) and is improved to a permille level at future electron-positron colliders such as the International Linear Collider (ILC) [4][5][6][7][8], the Compact Linear Collider (CLIC) [9,10], the Future Circular Collider of electrons and positrons (FCC-ee) [11], and the Circular Electron Positron Collider (CEPC) [12,13]. With such precision, we may be able to detect deviations in various coupling constants of the Higgs boson with a distinct pattern, by which we can fingerprint new physics models indirectly.
The shape of the Higgs potential can be directly reconstructed by measuring the triple Higgs boson coupling (the hhh coupling), which is expected to be determined with an order of one deviation at the HL-LHC. If the ILC with the center-of-mass energy of 1 TeV is realized, the error for the hhh coupling can be reduced to 10% [14][15][16], which is sufficient to test the scenario of electroweak baryogenesis. However, it has recently been discussed that the collision energy of the ILC is reduced to 250 GeV with the integrated luminosity to be 2 ab −1 [17] to make it a Higgs factory, where the Higgs boson decays can be measured very precisely while the measurement of the hhh coupling and the top Yukawa coupling are left for the far future. If this is the case, there may be little hope for the precise determination of the Higgs potential for a long time.
Fortunately, observation of gravitational waves (GWs) provides us with an exciting possibility of probing the early Universe well before the big bang nucleosynthesis. The detection of GWs from black hole binaries [18][19][20] and from neutron star mergers [21] has already signaled a new era of GW astronomy, and, in the future, space interferometers such as Laser Interferometer Space Antenna (LISA) [22] and Decihertz Interferometer Gravitational Wave Observatory (DECIGO) [23] will open up an era of GW cosmology. 1 Especially, the LISA project has already been approved and will start its operation in 2034, making it possible to test various extensions of the SM that predict stochastic GWs. The first-order phase transition is one of the best-motivated GW sources not only because it is a crucial element for successful electroweak baryogenesis, but also because the resulting GW spectrum is typically peaked around the interferometer frequency band: milli-to decihertz. Particle physics models which generate detectable GWs have been vigorously discussed by many authors , and the resulting GW spectrum has been studied in great detail from both analytic and numerical viewpoints . Therefore, by using this accumulated knowledge, we may be able to explore the Higgs potential through the observation of GWs at future space-based interferometers. Importantly, around the time of the LISA project, precision measurements of the Higgs boson couplings can be made at future collider experiments such as the ILC250 [17], and hence we expect a great synergy between GW observations and collider experiments.
Although many papers have investigated the possibility of detecting GWs from a phase transition at future experiments, most of them perform a relatively simple analysis in which it is discussed whether the predicted GW spectrum comes above or below the sensitivity curves. This type of analysis gives a rough estimate on what kind of models or which parameter space generates a detectable number of GWs. However, it cannot quantify to what extent the model parameters can be measured once GWs are detected, or what constraints can be derived when future experiments actually give us the data. In view of the recent growing interest in GWs, it is of great importance to study the attainable precision of the future GW experiments in exploring the Higgs sector and their complementarity to collider experiments.
In light of these considerations, in this paper, we adopt the method of Fisher matrix analysis and study expected constraints in future GW experiments such as LISA and DECIGO. We also consider an experiment like Big-Bang Observer (BBO) [124]. We investigate possible future constraints on parameters characterizing the spectral shape and those characterizing the properties of the transition. Since these quantities are determined by fundamental parameters in the underlying particle physics model, we can also estimate the expected sensitivities to such parameters. Then we compare or add them with possible future constraints from collider experiments to investigate the synergy between GW and collider experiments.
The organization of the paper is as follows. In Sec. II, we summarize our setup for the Fisher matrix analysis and explain how we constrain model parameters by assuming the specifications of future GW experiments such as LISA, DECIGO, and BBO. In Sec. III, we perform a Fisher analysis on a general peaky spectrum, taking the peak frequency, its amplitude, and spectral slopes as free parameters. In Sec. IV, we perform a Fisher analysis on transition parameters, i.e., α, β=H Ã , T Ã and so on (which we define later), using the GW spectral shapes in the literature. In Sec. V, we adopt specific particle physics models to illustrate that their model parameters can indeed be constrained by future GW experiments and discuss their complementarity to collider experiments. We finally conclude in Sec. VI. Some results based on different model setups are also presented in the Appendixes.

II. SETUP
In this section, we summarize the formalism adopted in our analysis. The GW spectrum from first-order phase transitions is also briefly discussed.

A. Gravitational-wave spectrum
Gravitational waves h ij are given as the transversetraceless part of the metric: In the following, we consider quantities such as the GW spectrum at the present time t ¼ t 0 and take aðt 0 Þ ¼ 1. We expand h ij as with λ being the label for GW polarization, and we impose the normalization condition ϵ λ ij ðnÞϵ λ 0 ij ðnÞ ¼ 2δ λλ 0 and the reality condition ϵ λÃ ij ðnÞ ¼ ϵ λ ij ðnÞ on the polarization tensor. Then GWs h λ satisfy h Ã λ ðf;nÞ ¼ h λ ð−f;nÞ from the reality of h ij . Now we define the power spectrum S h by hh λ ðf;nÞh Ã λ 0 ðf 0 ;n 0 Þi ¼ 1 16π δðf − f 0 Þδðn −n 0 Þδ λλ 0 S h ðfÞ: ð2:3Þ Here hÁ Á Ái denotes the ensemble average, and we assume that the two polarizations of GWs are uncorrelated and have the same amplitude. This power spectrum satisfies The intensity of GWs is also expressed by the ratio of their energy density to the critical energy density of the Universe. The former is given by (see, e.g., Ref. [125]) where M P ¼ ð8πGÞ −1=2 is the reduced Planck mass and hÁ Á Ái osc means taking both the ensemble average and the oscillation average. Note that the lhs does not depend on ⃗ x. Also, t is implicitly taken to be around the present cosmic age and omitted in the following. We decompose the total energy density into the contributions from each frequency as We define Ω GW to be the ratio of the GW energy density to the critical energy density ρ c of the present Universe:

B. Statistical analysis
In this subsection, we summarize the formalism we use for the statistical analysis for GW experiments. We use the Fisher matrix analysis, which is essentially a Gaussian approximation of the likelihood function. As we see below, the Fisher information matrix F ab is given by the curvature of the logarithm of this Gaussian-approximated likelihood around the fiducial parameter point. The inverse of this Fisher matrix gives the covariance matrix, which characterizes the uncertainties in the parameters.
In Secs. III-V, we assume LISA, DECIGO, and BBOlike (which we denote simply as BBO in the following) experiments. For cross-correlated detectors such as DECIGO and BBO (here we assume the cross-correlated DECIGO detector), the signal-to-noise ratio and δχ 2 , the latter of which is given by the logarithm of the likelihood function L, are calculated, respectively, as [126,127] S N Here T obs is the observation period and S h ðf; fpgÞ denotes the GW spectrum realized with a set of fundamental parameters fpg. In Secs. III-V, we take different parameter sets for fpg. Throughout this paper, fpg denotes fiducial values for fpg. Also, I and I 0 run over different interferometer channels. In addition, Γ II 0 is the overlap reduction function, which accounts for the insensitivity to the GW signal due to the geometry of detectors I and I 0 . 2 In Eqs. (2.9) and (2.10), σ II 0 in the denominator is given by In this expression, we included the effect beyond the weaksignal limit [135]. Also, σ ðnullÞ II 0 is given by taking S h → 0 limit in Eq. (2.11). The Fisher information matrix F ab , or the inverse of the covariance matrix hΔp a Δp b i, can be obtained from the expression (2.10) as (see, e.g., Ref. [127]) For the calculation of the overlap reduction function, see, e.g., Refs. [128][129][130]. For the Fisher analysis including the overlap reduction function with two units of triangular configuration, see, e.g., Refs. [127,[131][132][133][134].
Here ∂ p a denotes the derivative with respect to parameter p a . As a result, δχ 2 is approximated as This approximation is justified as long as Γ II , Γ I 0 I 0 , and Γ II 0 are of the same order. Now we discuss the case of LISA. LISA is a single detector, and therefore the above expression for crosscorrelated detectors may not be applied directly. As briefly discussed in Ref. [136], in an ideal case of autocorrelation, we may use an expression for the signal-to-noise ratio which is similar to cross-correlated cases. In this paper, we assume that this is indeed the case. The signal-to-noise ratio in such cases reduces to Here the labels I and I 0 drop, and also the factor of 2 drops compared to Eq. (2.9), because LISA has only one detector instead of two [136]. The corresponding expression for the likelihood becomes The resulting approximate expression for δχ 2 reduces to Eq. (2.13). In Secs. III-V, we sometimes show expected constraints in two-dimensional planes. When the number of fundamental parameters is more than two, the results are obtained after marginalizing over the parameters other than those shown in the figures by following the procedure below. Denoting the marginalized parameters fp ⊥ g collectively, we first construct marginalized likelihoodL by integrating out fp ⊥ g: It is understood that fpg in the lhs does not contain fp ⊥ g. Then the marginalized δχ 2 is given by the likelihood ratio as δχ 2 ðfpg; fpgÞ ¼ −2 lnL ðfpg; fpgÞ Lðfpg; fpgÞ : ð2:24Þ

C. Effective sensitivity and foregrounds
In this subsection, we clarify our assumptions on the effective sensitivity S eff in Eqs. (2.14) and (2.20) and also explain astrophysical foregrounds which enter S h in addition to the signal we would like to observe. 3 This common definition does not take into account a relatively large factor T obs R df ∼ T obs × f typ (with f typ being the typical peak frequency of the GW spectrum) which appears in Eq. (2.15). To take this into account, one may instead use the power-law sensitivity curve; see Ref. [136].

Effective sensitivity
We use the fitting formulas in Ref. [137] for LISA 4 and the ones in Ref. [141] for DECIGO and BBO-like experiments.

Foregrounds
It is known that GWs from astrophysical sources form unresolvable foregrounds. In this paper, we incorporate their effects by including the following power spectrum to S h in addition to the signal from first-order phase transitions. 5 In the analysis in Secs. III-V, we assume that the spectral form of these foregrounds is already known from other studies and do not consider their uncertainties.
One such foreground comes from compact white dwarf binaries in our Galaxy in the millihertz regime. The noise spectrum adopted in Ref. [137] is In our analysis, we use the following smoothened noise spectrum: Also note that S 0 WD above corresponds to the foreground to 4 More up-to-date sensitivity curves after the LISA Pathfinder results [138] are available in, e.g., Refs. [139,140]. We checked in several examples in the following analysis that our results do not change drastically even if we use the sensitivity curve in Ref. [138]. 5 It should be noted that these astrophysical foregrounds are correlated among the detectors and their treatment might be modified in a more realistic situation. the N2A5 configuration of LISA in Ref. [137] and therefore might not be applicable to DECIGO and BBO in a strict sense. However, we adopt this expression also for these detectors as a reference value.
Another source of foreground is binary neutron stars and binary black holes. As recently discussed in Ref. [142], the merger rate of neutron stars and black holes inferred from the detections of GWs by the LIGO and Virgo Collaboration might lead to a significant amount of foreground to stochastic GWs. However, since there are still large uncertainties in this foreground, we do not take this into account in the results presented in Secs. III-V. For completeness, in Appendix A we show the results including this foreground by adopting the following function given in Ref. [142]: ð2:33Þ for 1 Hz < f < 10 3 Hz. Though in Ref. [142] the foreground is shown only for 10 Hz < f < 10 3 Hz, we have slightly extrapolated it down to 1 Hz to make conservative estimates.

D. GW spectrum from first-order phase transitions
In this subsection, we summarize the spectral form of GW signal from first-order phase transitions. In cosmological first-order phase transitions, bubbles of true vacuum first nucleate at some temperature, and then they expand due to the pressure difference between the true and false vacua. They eventually collide and merge with each other, and during this phase GWs are sourced by the energymomentum tensor of the system. The dynamics is mainly determined by the following parameters: Here T Ã is the temperature of the Universe just after the phase transition, η is the (symbolically denoted) coupling of the scalar field to the surrounding plasma, and α ≡ ρ 0 =ρ rad is the ratio between the released latent heat ρ 0 and the background plasma energy density ρ rad at the time of transition. Also, β=H Ã ¼ dðS 3 =TÞ=d ln Tj T¼T N is the logarithmic temperature derivative of the three-dimensional bounce action with H Ã being the Hubble parameter at the time of the transition. This quantity determines the bubble nucleation rate. The η dependence can be translated to the dependence on bubble wall velocity v w through the relation in Ref. [143], 6 because the wall velocity is determined by the balance between the released energy and the friction on the walls. Therefore, instead of the parameter set (2.34), we consider the following one in the analysis below: As a result of the scalar and plasma dynamics mentioned above, three types of GW sources arise (e.g., Refs. [46,145]): (i) bubble collisions, (ii) sound waves, and (iii) turbulence. The first one comes from the collision of walls, i.e., scalar field configurations. This contribution is well approximated by the envelope of the configurations with infinitely thin shells [100][101][102][103]. More recently, the resulting GW spectrum has been calculated by many-bubble simulations [106,121,122,146] and also by an analytic approach [119]. 7 This scalar field contribution becomes significant when the bubble walls run away [147], which occurs when the friction from the thermal plasma on the walls cannot stop the acceleration of the walls. However, it has recently been pointed out that such runaway bubbles are unlikely after taking into account particle splitting processes around the walls [148]. Therefore, this scalar contribution now is not considered to be a dominant source of GWs.
The second contribution arises from the dynamics of the fluid, in contrast to bubble collisions. During bubble expansion, a significant fraction of the released energy is converted to the bulk motion of plasma surrounding the walls. This plasma motion is launched into free propagation after bubbles collide with each other, and it propagates as sound waves at the level of linear approximation. These sound waves have been found to continuously source GWs with wave numbers corresponding to the thickness of the bulk fluid [113,115,117], and it has been proposed to model this GW production by sound shells [149]. The resulting GW spectrum is [117] ð2:37Þ 6 In the deflagration case, the released energy heats up the plasma in front of the bubble walls. This heating backreacts on the walls and decreases the pressure exerted on them, and as a result the wall velocity can change as the transition proceeds (see, e.g., Ref. [144]). In this paper, we do not consider such effects to make our analysis simple. 7 This approach also gives a rough estimate on the dependence of the GW spectrum on the nucleation rate, which can be used to distinguish particle physics models once we observe GWs from first-order phase transitions [120]. :

ð2:38Þ
Here g Ã is the number of relativistic degrees of freedom, which we take to be 106.75 throughout the paper. Also, κ sw is the fraction of the released latent heat which goes into the plasma bulk motion and contributes to sound-wave formation. The peak frequency f sw comes from the aforementioned thickness of the sound shell. The last contribution, turbulence, arises when the sound waves develop into the nonlinear regime at late times. In this paper, we adopt the spectral form given in Refs. [110,150], based on the Kolmogorov-type turbulence proposed in Ref. [151]: ; ð2:41Þ : ð2:42Þ Here κ turb ¼ ϵκ sw is the fraction of the released latent heat which goes into the turbulent motion of the plasma. The value of ϵ takes at least ϵ ¼ ð0.05-0.1Þ as confirmed within the numerical simulation time [46,115] but is still quite unknown. In our analysis we mainly fix ϵ ¼ 0.1 to give conservative estimates, though we also compare between ϵ ¼ 0.1 and 0.5 in Sec. IV using Fig. 10. We note in passing that the estimation of GW spectrum resulting from a first-order phase transition is an ongoing hot topic (e.g., Refs. [117,118,[120][121][122][123]149,152,153]), and therefore the above spectra might not be exact. For example, the long-lastingness of sound waves, which is assumed in deriving Eqs. (2.36)-(2.38), is somewhat questioned in a recent study [153]. Also, Eqs. (2.36)-(2.38) have only a single frequency scale coming from the sound shell thickness, which may miss possible infrared structures [118]. In addition, the turbulence fraction ϵ still has huge uncertainty as mentioned above. However, an important point is that the phase transition dynamics is, in principle, determined by a few parameters, and the resulting GW spectrum is also determined by them accordingly. Therefore, it would be interesting to ask what kind of information we can obtain by observing GWs if we know the exact form of the GW spectrum, which depends on the parameters related to the transition. In this paper, we illustrate this point by taking Eqs. (2.36) and (2.39) as examples.

III. FISHER ANALYSIS ON GENERAL SPECTRUM
In this section, we first perform a Fisher analysis on a general peaky GW spectrum, taking the peak amplitude, peak frequency, and spectral slopes as free parameters. 8 We assume that the signal takes the following form: ðf=f peak Þ n R ðf > f peak Þ: ð3:1Þ We also assume n L > 0 and n R < 0.
We first show the result of a Fisher analysis using δχ 2 given in Eqs. (2.15) and (2.21) with the effective sensitivities (2.25)-(2.30). We take several fiducial values for the parameters ðf peak ; Ω GW;peak ; n L ; n R Þ as examples. The sample points we consider are (i) point 1: ðf peak ; Ω peak Þ ¼ ð10 −2 Hz; 10 −7 Þ; (ii) point 2: ðf peak ; Ω peak Þ ¼ ð10 −1 Hz; 10 −10 Þ; (iii) point 3: ðf peak ; Ω peak Þ ¼ ð10 Hz; 10 −10 Þ; and (iv) point 4: ðf peak ; Ω peak Þ ¼ ð10 −1 Hz; 10 −14 Þ. For the spectral slope, we consider two cases with ðn L ; n R Þ ¼ ð3; −4Þ and ð1; −3Þ. The former corresponds . In these figures, we marginalize the two spectral indices n L and n R following the procedure in Sec. II B and show contours of f peak =f peak − 1 and Ω GW;peak =Ω GW;peak − 1 for δχ 2 ¼ 2.3, which corresponds to 1σ in the two-dimensional plane. The three contours in each panel correspond to T obs ¼ 1, 3, and 10 yr, respectively. Also, the panel is enlarged when the size of outermost ellipse (which corresponds to T obs ¼ 1 yr) far exceeds unity. For point 1, we expect parameter determination with a good precision with all three detectors. For points 2-4, parameter determination by LISA is challenging, but we still expect a good sensitivity for DECIGO and BBO. It is seen that for point 4 DECIGO can perform well, even though the signal is below the sensitivity curve. This is understood through Eq. (2.15): even if the signal S h is below the sensitivity of the detector S eff , we have an additional factor ðT obs × f typ Þ ∼ ðOð1Þ yr × 0.1 HzÞ with f typ being the typical signal frequency for the fiducial values off peak andΩ GW;peak .
In Fig. 4, we show contours for 1σ fractional error Δf peak =f peak and ΔΩ GW;peak =Ω GW;peak (where Δf peak and ΔΩ GW;peak correspond to δχ 2 ¼ 1 for one degree of freedom) after marginalizing the other parameters. Also, Fig. 5 shows the same contours with the foreground chosen to be S WD þ S NSBH , as mentioned in Sec. II C. It is seen that the precision gets worse for the latter for the frequency range 1 Hz ≲ f ≲ 10 3 Hz, but there are still possibilities for parameter determination for Ω GW;peak ≳ 10 −12 .

IV. FISHER ANALYSIS ON TRANSITION PARAMETERS
In this section, we perform a Fisher analysis on the transition parameters (2.35) with the spectrum provided in  In Fig. 6, we show the signal with these parameter points as well as the sensitivity curves and the foreground from white dwarfs.
Before showing the results, it should be mentioned that, if we take all four parameters in Eq. (2.35) completely free, it is generically difficult to determine their values at the same time. This is because of the following reason. Suppose that the detector sees only the sound-wave peak, Eq. (2.36). (Notice that the sound-wave peak amplitude is typically much larger than that of turbulence as long as we adopt the spectral shapes in Sec. II.) For the spectral shape given by Eq. (2.36), the information the detectors can obtain is the position (i.e., frequency and amplitude) of the peak, which is not enough to determine all four parameters. Therefore, in the analysis below, we limit the number of free parameters to two (α and β=H Ã ) or to three (α, β=H Ã and T Ã ). When we show two-parameter planes in a threeparameter analysis, we marginalize over T Ã following the procedure in Sec. II B. Figures 7-9 are the results of a Fisher analysis for the three fiducial points above. In these figures, the left and right columns correspond to two-and three-parameter analyses, respectively, and the three contours in each panel correspond to T obs ¼ 1, 3, and 10 yr, respectively. Also, the panel is enlarged when the size of the outermost ellipse (which corresponds to T obs ¼ 1 yr) far exceeds unity. We observe the following.
(i) First, for the two-parameter analysis (left columns in Figs. [7][8][9], it is seen that the parameters are well determined except for point C for LISA. Even for point C for LISA, one combination of α and β=H Ã is determined well. One may wonder why the parameters can be fixed even when the signal barely touches the sensitivity curves. This is because of the same reason as Sec. III: we have an additional factor ðT obs × f typ Þ ∼ ðOð1Þ yr × 0. , the shape of the GW spectrum (i.e., structure of the lowand high-frequency tails) is almost independent of the transition parameters (except for the weak dependence of the turbulence spectrum on h Ã ). Therefore, the only information which the detectors can obtain is the peak position (i.e., peak frequency and amplitude) of each component. This means that the detectors can determine maximally two transition parameters if they see only one component. On the other hand, when the detectors see both sound wave and turbulence components, they can infer the third parameter from the different dependences of the two components on the transition parameters. To illustrate this point, in Appendix B, we perform a three-parameter analysis, changing the strength of the turbulence component by hand. We indeed see that the contours start to close and degeneracies disappear as we increase the turbulence fraction ϵ. Therefore, we reach one important conclusion: in light of parameter determination, observing subdominant components is essentially important.
Given this, let us gauge the uncertainty from the turbulence fraction ϵ by using the two-parameter analysis. In Fig. 10 by orders of magnitude, their shapes change in some cases. General tendencies can be read off together with Fig. 6: the shapes of the contours are determined by the sound wave signal alone when the turbulence signal is far below the sensitivity curves (LISA for point A, LISA and DECIGO for point B), while the latter signal contributes to narrow down the contours once it starts to come around or above the sensitivity curves (DECIGO and BBO for point A, BBO for point B). From this observation, we expect that the contour shapes do not change for point C for ϵ ¼ 0.1 and 0.5, which we indeed confirmed. Figures 11 and 12 show 1σ fractional error for α and β=H Ã for a two-and three-parameter analysis, respectively. In these figures, we fixed T Ã ¼ 100 GeV. It is seen that a high-sensitivity spot appears in some of the panels. This is because, if one fixes α, there is a typical value of β=H Ã which makes the signal peak close to the frequency at which the detector is most sensitive. (Note that for too small β=H Ã the signal from sound waves starts to overlap with the foreground from white dwarfs.) Also, it is seen that the sensitivity on α becomes worse as α increases for fixed β=H Ã . This is because the spectral shape, Eqs. (2.36) and (2.39), becomes almost independent of α for α ≫ 1. Physically, this means that the transition dynamics looks almost the same when the released latent heat dominates the radiation energy density (i.e., α ≫ 1).

V. FISHER ANALYSIS ON MODEL PARAMETERS
In this section, we take some specific examples of particle physics models which give rise to a first-order phase transition in the early Universe and illustrate how the detection of GWs contributes to narrow down the fundamental model parameters. Below, we consider (i) models with additional isospin singlet scalar fields with and without the classical conformal invariance, (ii) a model with an extra Higgs singlet field, and (iii) a classically conformal B − L model, respectively.

A. OðNÞ singlet extensions of the SM
We first consider extensions of the SM in which N additional isospin singlet scalars ⃗ S ¼ ðS 1 ; …; S N Þ T with a global OðNÞ symmetry are added to the SM particle content. There are two classes in such extensions: with or without classical conformal invariance (CCI).

Model
The tree-level potential of the OðNÞ models with CCI is given by where Φ is the isospin doublet Higgs field. In the CCI models, the electroweak symmetry breaking occurs by the Coleman-Weinberg mechanism [155]. Predictions of this class of models have been studied in, e.g., Refs. [156,157], and a distinctive phenomenological feature has been found: the deviation in the hhh coupling is universally about 70% [157]. Gravitational-wave production in this class of models has been studied in, e.g., Refs. [42,51]. In the following analysis, we take the free parameters to be N and λ S . This is reasonable, because two of the original four parameters λ Φ , λ S , λ ΦS , and N are fixed by the observed Higgs mass m h and its vacuum expectation value v.
On the other hand, the tree-level potential for the OðNÞ models without CCI is given by ð5:2Þ Compared to the above models, we have two additional parameters μ 2 and μ 2 S . This makes the number of the free parameters four instead of two. In the following analysis, we take the free parameters to be N, λ S , m S , and μ 2 S , where m S is the singlet mass after the transition. In this class of models, m S can be translated into the deviation in the triple Higgs coupling Δλ hhh =λ SM hhh ≡ λ hhh =λ SM hhh − 1.

Analysis
We take the following benchmark points for the models with and without CCI, respectively: (i) ðN; λ S Þ ¼ ð2; 0.1Þ.-with CCI; and (ii) ðN; λ S ; m S ½GeV; μ 2 S ½GeV 2 Þ ¼ ð8; 0.1; 385; 0Þ and (12,0.1,385,0).-without CCI. Several comments are in order before moving on to the results. For the former model, we perform a two-parameter analysis with N and λ S . In this analysis, we regard N as a continuous parameter to make Eq. (2.22) directly applicable. For the latter model, we fix μ S to the fiducial value and perform a three-parameter analysis with N, λ S , and m S . When showing the final figures, we marginalize over N and translate m S into the triple Higgs coupling Δλ hhh =λ SM hhh as mentioned above. Finally, for both models, we assume fixed values for v w , since it is generically hard to calculate the wall velocity in a given model. Therefore, the GW spectrum used in our analysis reflects the model parameters only through α, β=H Ã , and T Ã . The dependence through v w should be taken into account in more realistic analyses.
In the following, we would like to illustrate the ability of GW detectors to pin down the model parameters. However, it is known that smaller values of v w (which are favored for producing the baryon asymmetry of the Universe [3]) are generically in tension with a sufficient amount of GW production for observations (however, see, e.g., Ref. [35]). Therefore, we assume different wall velocities depending on the observation to stress the following: LISA is already powerful in observing GWs from first-order phase transitions, especially when the wall velocity is somewhat larger. On the other hand, if we keep the generation of the baryon asymmetry of the Universe in mind, DECIGO and BBO dramatically enlarge the detection possibility even when the wall velocity is smaller. Figure 13 is the GW spectrum realized in each model. The corresponding electroweak phase transition parameters are (i) ðα; β=H Ã ; T Ã ½GeVÞ ≃ ð0.080; 1000; 82Þ.-with CCI; and (ii) ðα; β=H Ã ; T Ã ½GeVÞ ≃ ð0.10; 1700; 83Þ and (0.14, 1600,77).-without CCI. In this figure, we assumed v w ¼ 0.95 (top panels) and v w ¼ 0.1 (bottom panels). We use the former value of v w for LISA and the latter for DECIGO and BBO, respectively.
The results of a Fisher analysis are shown in Figs. 14 and 15, where Fig. 14 is for LISA while Fig. 15 is for DECIGO and BBO, respectively. As seen in the left panel in Fig. 14, LISA has the potential to contribute to narrow down the parameters for the model with CCI, even if the spectrum is somewhat below the sensitivity curve. This is because of the same reason as Sec. III: we have a factor of ðT obs × f typ Þ ∼ ðOð1Þ yr × 0.1 HzÞ in Eq. (2.15) with f typ being the typical frequency of the signal. However, in passing, it should be again noted that we have assumed the ideal case discussed in Ref. [136]. Also, for the model without CCI, it is somewhat challenging to constrain the model parameters, as seen in the right panel in the same figure. On the other hand, DECIGO and BBO perform excellently to pin down the model parameters as shown in Fig. 15. For the model with CCI, both N and λ S can be determined with a good precision even for v w as low as 0.1. Note that, after restricting N to be an integer, the uncertainty in λ S becomes significantly small. For the model without CCI as well, though degeneracy appears because of the relatively large number of model parameters, both detectors can contribute to determine a certain combination of the parameters.   9 Even in parameter regions where only a small number of GWs are produced from the transition dynamics, deformations in the primordial GW spectrum might also help pin down the model parameters: see, e.g., Refs. [158,159].
Finally, we show in Fig. 16 the result of a Fisher analysis in the α-β=H Ã plane in order to discuss the complementarity between collider and GW experiments. In this figure, we perform a Fisher analysis on α and β=H Ã with T Ã ≃ 93 GeV fixed at the fiducial value and also v w fixed to be 0.95. (In Appendix A, we show the results after marginalizing T Ã .) The resulting 1σ contours for LISA are shown as the red and blue lines for the OðNÞ models with and without CCI, respectively. The three contours correspond to T obs ¼ 1, 3, and 10 yr, respectively. Both the left and right panels use the fiducial point N ¼ 2 for the OðNÞ models with CCI, while they use N ¼ 8 (left) and N ¼ 12 (right) with ffiffiffiffiffi μ 2 S p ¼ 0 GeV for the OðNÞ models without CCI. In the latter models, we choose m S so that the triple Higgs coupling has the same value as the former: Δλ hhh =λ SM hhh ¼ 66.7%. As seen from the left panel, LISA may be able to distinguish N ¼ 2 with CCI from N ¼ 8 without CCI for T obs ¼ 10 yr, even when collider experiments cannot distinguish the two classes from the triple Higgs coupling. On the other hand, as seen from the right panel, LISA may differentiate N ¼ 2 with CCI and N ¼ 12 without CCI even in shorter observation periods.

B. Real Higgs singlet extension of the SM
We next consider an extension of the SM with a real singlet scalar field S which takes a nonzero expectation value at low temperatures. FIG. 16. LISA 1σ contours for the OðNÞ singlet models with and without CCI in the α-β=H Ã plane. The red points correspond to N ¼ 1, 2, 4, 12, and 60 for the OðNÞ models with CCI from left to right, while the gray points correspond to N ¼ 1, 2, 4, 12, and 60 with μ 2 S ¼ 0 GeV 2 for the OðNÞ models without CCI. The fiducial values for the red contours correspond to the OðNÞ models with CCI with N ¼ 2, while the ones for the blue contours correspond to the OðNÞ models without CCI with N ¼ 8 (N ¼ 12) and μ 2 S ¼ 0 GeV 2 in the left (right) panel. Also, the three contours correspond to T obs ¼ 1, 3 and 10 yr, respectively.

Model
The tree-level potential of this model is given by One of the eight parameters in the model can be removed by the redefinition of the singlet scalar field. In the following analysis, we take μ S to be 0 by the field redefinition of S. The electroweak phase transition in this model involves not only tree-level mixing effects between the scalar fields but also thermal loop effects. If the transition is strongly first order, the latter effects are imprinted in the resulting shape of the GWs. Therefore, we can test the model both by precision measurements of various Higgs boson couplings and by GW observations [56].
In the following analysis, we take the free parameters of the model to be the mass eigenvalue m H of the additional singlet scalar eigenstate H, the deviation κ in the Higgs couplings to the gauge bosons and fermions, and the vacuum expectation value v S of the singlet scalar field, μ ΦS and μ 0 S .

Analysis
As explained above, there are five parameters in this model: m H , κ, v S , μ ΦS , and μ 0 S . However, in the following, we fix v S and μ 0 S at their fiducial values and include only the parameters related to the H field (m H , κ, and μ ΦS ) in the analysis. This is because unremovable degeneracies appear if we take all the parameters as free. In other words, observing the GW spectrum may not be enough to pin down the model parameters. However, the detection of GWs indeed contributes to narrowing down the allowed parameter space, as we see below. Also, for the wall velocity v w , we fix its value as in Sec. VA. We take the following benchmark point for all of LISA, DECIGO, and BBO: (i) ðm H ½GeV; κ; μ ΦS ½GeV; v S ½GeV; μ 0 S ½GeVÞ ¼ ð166; 0.96; −80; 90; −30Þ. The electroweak phase transition parameters for this fiducial point become (i) ðα; β=H Ã ; T Ã ½GeVÞ ≃ ð0.085; 420; 93Þ. For LISA, we fix the wall velocity to be v w ¼ 0.95. The GW spectrum realized at this parameter point is shown in the left panel in Fig. 17. The result of a Fisher analysis is shown in the right panel in the same figure. The narrow contours are for fixed μ ΦS , while the wide contours are obtained after marginalizing over μ ΦS . It is seen that LISA can contribute to constraining m H and κ with a good accuracy.
For DECIGO and BBO, we fix the wall velocity to be v w ¼ 0.1. The GW spectrum realized at this parameter point is shown in the top panel in Fig. 18. The results of a Fisher analysis is shown in the bottom panels in the same figure. The narrow contours are the result for fixed μ ΦS , while the wide contours are the ones after marginalizing over μ ΦS . It is seen that, though the GW amplitude is much smaller than the previous parameter point, DECIGO and BBO can perform excellently in constraining the model parameters.
Finally, we see the synergy between collider and GW experiments in Fig. 19. The condition for a strongly firstorder electroweak phase transition is satisfied in the greenshaded region, and the same result as the right panel in Fig. 17 is shown in blue. The yellow region is the 1σ expected sensitivity of ILC with ffiffi ffi s p ¼ 250 GeV and L ¼ 2 ab −1 [17], while the right-bottom shaded region is excluded by the current experimental data for the direct search for a heavy Higgs [160]. It is seen that we can test the model by both future GW experiments and collider experiments.

C. Classically conformal B − L model
We next consider the classically conformal B − L model proposed in Refs. [161,162] based on the argument on classical conformal theories [163]. It is known that, in nearly conformal models, a large amount of GWs can be produced due to huge supercooling and slow change of the nucleation rate (see, e.g., Refs. [52,72,85,88,164,165]). Gravitationalwave production in the classically conformal B − L model has been studied in Refs. [52,72] for relatively large and small B − L gauge coupling, respectively. In the following analysis, we consider the former parameter region.

Model
The relevant part of the model is the scalar sector, whose tree-level potential is given by where only four-point couplings appear due to the assumption of the classical conformal symmetry. Here Φ is the SM Higgs doublet, and X is the B − L breaking scalar with B − L charge þ2. 10 The B − L scalar field X develops the vacuum expectation value M ≡ ffiffi ffi 2 p hXi due to the running of the coupling λ X . The mixing term   [160]. 10 In this model, there are also right-handed neutrinos. In this paper, we neglect their effects, assuming that their Yukawa couplings are small enough. λ ΦX jΦj 2 jXj 2 generates the negative mass term for the SM Higgs, and electroweak symmetry breaking is realized at zero temperature. We consider the parameter space where M is relatively larger than the electroweak scale. In such cases, the mixing coupling λ ΦX becomes negligible, and the potential for the X field is mainly determined by the B − L gauge interaction.
In this scenario, the phase transition in the X direction occurs in the early Universe and produces a large number of GWs. To understand this, first let us consider the situation where the scalar field X is trapped at the origin. Then the finite-temperature effective potential roughly consists of the thermal mass term and the temperature-dependent quartic coupling: Re½X parametrizing the transition direction. Also, both the B − L gauge coupling g B−L and the quartic coupling λ X are understood as dependent on the temperature T of the Universe. For a low enough temperature, the effective quartic coupling becomes negative and the origin X ¼ 0 becomes the false vacuum. Then, the resulting tunneling rate can be written just by the combination of the couplings, because there is no scale other than the temperature T: Note that it is only logarithmically dependent on the temperature. As a result, the parameter β=H Ã ¼dðS 3 =TÞ=dlnTj T¼T N becomes relatively small, and we expect a large amount of GW production. After the transition, the χ field settles down to the true vacuum generated by the Coleman-Weinberg mechanism [155]. (ii) green: excluded by Z 0 search (see Refs. [166]); and (iii) yellow: phase transition does not complete in sufficiently large regions (see Refs. [52,72]). It is seen that a significant supercooling occurs in this model (α ≫ 1). In such cases, the combustion mode of the  walls is likely to be a very strong detonation, where the wall velocity approaches almost unity. 11 (Note that even in this case most of the released energy is still carried by the fluid motion [148].) Therefore, in the following analysis, we fix v w ¼ 1.
Below, we see that, for particle models with such a small number of free parameters, the detection of GWs significantly contributes to pin down the model parameters.

Analysis
We first take two fiducial points:   Indeed, the walls reach a terminal velocity and the transition proceeds with detonation when is satisfied [118]. Here g typ , m typ , and χ typ are the typical gauge coupling, mass scale of the potential, and the distance between the false and true vacuum, respectively. While the first three factors in the right-hand side take ∼ð0.01-1Þ values, the last factor gives a huge number ∼10 10-16 .
The GW spectra realized for these parameter points are shown in Fig. 22. 12 It is seen that the resulting GW amplitude is extremely large due to the behavior of α and β=H Ã shown in Fig. 21. We show the result of a Fisher analysis for these parameter points in Fig. 23. It is seen that the model parameters are precisely determined except for point 2 with LISA, in which case the peak frequency of the GW spectrum becomes relatively high. However, even in such a case, GW detection still contributes to constraining the parameters as we see in the top-right panel in Fig. 23. Next, we show contour plots for ΔM=M and Δα B−L =α B−L for different fiducial values forM and α B−L in Fig. 24. The three rows show LISA, DECIGO, and BBO from top to bottom, respectively. It is seen that LISA can pin down the model parameters in a wide range of the parameter space, while such a parameter space becomes much wider for DECIGO and BBO. Also note that Z 0 searches can corner the parameter space from lower values of M, which is favored from the viewpoint of naturalness.

VI. DISCUSSION AND CONCLUSIONS
In this paper, we investigated to what extent future space-based GW detectors such as LISA, DECIGO, and BBO(-like one) can contribute to pin down new physics beyond the standard model through the detection of GWs from a first-order phase transition. In order to go beyond the naïve comparison between the GW signal and the sensitivity curves and quantify the attainable precisions, we adopted the method of Fisher analysis in this paper.
First, in Sec. III (and Appendix A 1), we studied the sensitivity of the detectors to the parameters which characterize a general peaky spectrum. We parameterized the spectrum with the peak frequency, peak amplitude, and spectral indices. We performed a Fisher analysis to see the attainable uncertainties, and it was found that not only ultimately sensitive detectors such as DECIGO and BBO but also LISA, a relatively near-future detector, can significantly contribute to study GW spectral shapes.
Next, in Sec. IV (and Appendix A 2), we performed a Fisher analysis on the parameters which characterize the phase transition such as the latent heat fraction α, time dependence of the bubble nucleation rate β=H Ã , and the transition temperature T Ã . We adopted a classification of GW sources in a first-order phase transition in the literature (i.e., bubble collisions, sound waves, and turbulence) and used the spectral shapes provided there. Though the classification and determination of the spectral shapes realized in a first-order phase transition is a still ongoing hot topic (e.g., Refs. [117,118,[120][121][122][123]149,152]), we illustrated in this paper the procedure to determine the transition parameters by detecting one or several spectral shapes which have different parameter dependences by adopting expressions in the literature. As a result, it was found that, though the detection of single spectral shape is indeed helpful, the degeneracies in the parameters are resolved and their precise determination is possible if more than one spectral shape is detected.
Finally, in Sec. V (and Appendix A 3), we studied how the detection of GWs contribute to the determination of fundamental model parameters. This is possible because the transition parameters above are determined by the parameters of the particle physics model which drive a first-order phase transition. We illustrated this point by taking three examples: (i) models with additional isospin singlet scalars, (ii) a model with an extra real Higgs singlet, and (iii) a classically conformal B − L model. We found that the detection of the GW spectrum is indeed extremely powerful in pinning down the model parameters. However, the exploration of new physics becomes truly interesting when GW searches are combined with collider experiments. We also illustrated this point by taking the above three examples. For the first models, the determination of the triple Higgs coupling helps to identify the existence of the classical scale invariance. However, even if its value takes similar values for the cases both with and without the classical scale invariance, GW detection can distinguish the two (Fig. 16). In this sense, we can narrow down the model candidates and finally identify one by using two different experimental methods. For the second model as well, colliders such as ILC give different constraints than GW observations (Fig. 19). For the last model, GW observations can corner the model with the help of Z 0 searches (Fig. 24).
In summary, we found that future gravitational-wave observations can play complementary roles to future collider experiments. Fortunately, the LISA project and precision measurements of the Higgs boson couplings come around the same time in the future: a great synergy between GW observations and collider experiments is awaiting us. In this Appendix, we show numerical results for different parameter sets or different noise assumptions from the ones in the main text.

Fisher analysis on general spectrum
This subsection supplements the results in Sec. III. We see how they change depending on the spectral indices ðn L ; n R Þ and the foregrounds.
We change the spectral indices in Eq. (3.1) to ðn L ; n R Þ ¼ ð1; −3Þ and see how the result changes. In fact, these spectral slopes are suggested in the study of GW production from thin bubbles [118,121]. First we show the sensitivity curves and GW signals in Fig. 25. The signals now have broader peaks compared to Fig. 1.
We show the results of a Fisher analysis in Figs. 26 and 27. Figure 26 corresponds to the case with S h ⊃ S WD only, while Fig. 27 corresponds to the one with S h ⊃ S WD þ S NSBH .

Fisher analysis on transition parameters
This subsection supplements the results in Sec. IV. We see how they change depending on the wall velocity v w and the foregrounds.

a. v w = 1
We first show how the result change if we include S NSBH , i.e., S h ⊃ S WD þ S NSBH , keeping v w ¼ 1 as in Sec. IV. Figures 28 and 29 are two-and three-parameter analyses, respectively. These figures are practically the same as Figs. 11 and 12. This is because, for T Ã fixed at 100 GeV, the GW spectrum tend to have its peak below 1 Hz, and therefore the foreground S NSBH in Eq. (2.33) becomes almost irrelevant.
b. v w = 0.3 We next discuss how the result changes if we take v w ¼ 0.3. Figures 30 and 31 are the results of two-and three-parameter analyses, respectively, with the foreground from white dwarfs only. On the other hand, Figs. 32 and 33 are the results of two-and three-parameter analyses, respectively, with the foreground from neutron stars and black holes also included. It is seen that the parameter region shifts towards lower β=H Ã compared to the v w ¼ 1 case. This is because lower v w makes the peak frequency higher, while lower β=H Ã compensates that by shifting the peak to a lower frequency. It is also seen that Figs. 30 and 31 and Figs. 32 and 33 are almost the same. This is because of the same reason as the previous subsection: for T Ã fixed around 100 GeV, there is essentially no effect of S NSBH .

Fisher analysis on model parameters
This subsection supplements the results in Sec. V.
a. OðNÞ singlet extension of the SM In Figs. 34 and 35, we show the results of a Fisher analysis after marginalizing over the temperature just after the transition T Ã . The bands stretch in the β=H Ã direction, but we still have the possibility to distinguish the models, as seen in the right panel in Fig. 35.

b. Classically conformal B − L model
In Fig. 36, we show the result of a Fisher analysis corresponding to the result in Sec. V C after including S NSBH . It is seen that the result changes only slightly compared to Fig. 24. This is because the number of GWs produced in this model is so large that they dominate the foreground S NSBH . The foreground is taken to be S WD þ S NSBH .

APPENDIX B: RESOLVING THE PARAMETER DEGENERACY
In this Appendix, we illustrate what we mentioned in Sec. IV, that is, observing the turbulence component is essentially important in the full determination of the transition parameters. For this purpose, we take points A and B for LISA and perform the three-parameter analysis, changing the turbulence fraction ϵ by hand.
The left and right panels in Fig. 37 are for points A and B in Sec. IV, respectively, with the dotted lines being the turbulence signal (2.39) with ϵ ¼ 0.1 (blue), 1 (red), and 10 (green). The green solid, red solid, and blue dashed lines are the same as in Fig. 6: they represent the LISA sensitivity curve, white dwarf foreground, and the signal from sound waves. The corresponding results for the threeparameter analysis are shown in Fig. 38. As we mentioned in Sec. IV, while there exists a tight degeneracy when ϵ is small, the contours start to close as we increase the turbulence fraction ϵ. Therefore, we conclude that observing different components in the GW spectrum can be essentially important in determining the transition parameters.
FIG. 37. LISA sensitivity curve (green solid), white dwarf foreground (red solid), and signals from sound waves (blue dashed) and from turbulence (blue, red, and green dotted). The left and right panels are for point A and B in Sec. IV, respectively. For the turbulence signal, we change ϵ ¼ 0.1, 1, and 10 from bottom to top.