Di-Higgs Production in the $4b$ Channel and Gravitational Wave Complementarity

We present a complementarity study of gravitational waves and double Higgs production in the $4b$ channel, exploring the gauge singlet scalar extension of the SM. This new physics extension serves as a simplified benchmark model that realizes a strongly first-order electroweak phase transition necessary to generate the observed baryon asymmetry in the universe. In calculating the signal-to-noise ratio of the gravitational waves, we incorporate the effect of the recently discovered significant suppression of the gravitational wave signals from sound waves for strong phase transitions, make sure that supercooled phase transitions do complete and adopt a bubble wall velocity that is consistent with a successful electroweak baryogenesis by solving the velocity profiles of the plasma. The high-luminosity LHC sensitivity to the singlet scalar extension of the SM is estimated using a shape-based analysis of the invariant $4b$ mass distribution. We find that while the region of parameter space giving detectable gravitational waves is shrunk due to the new gravitational wave simulations, the qualitative complementary role of gravitational waves and collider searches remain unchanged.


I. INTRODUCTION
The first direct detection of the gravitational waves (GW) by the LIGO and Virgo collaborations [1] has triggered a revived interest in using the stochastic GW from the electroweak phase transition (EWPT) to learn more about particle physics, in particular, to probe possible hints of physics beyond the standard model (BSM). A stochastic background of GW can be produced from a cosmological first-order phase transition. The EWPT is required to be first-order to provide a non-equilibrium environment for generating the observed baryon asymmetry in the universe, in the framework of electroweak baryogenesis (EWBG) (see [2] for a recent review). EWBG is one of the popular mechanisms for solving the long-standing baryon asymmetry problem. In this framework, BSM physics provides new sources of CP-violation [3][4][5] and facilitates strong first-order EWPT, both of which cannot be achieved within the SM, though it is possible to violate the baryon number in the SM through the weak Sphaleron process. Therefore, GW measurements can provide a new window to BSM physics.
A relevant theoretical benchmark construction is the so-called "xSM", which is a minimal extension of the SM by adding a new gauge singlet scalar. The xSM has been under extensive phenomenological studies due to its simplicity and also been used to study the GW signals due to its ease to accommodate a strong first-order EWPT [6][7][8][9][10][11][12]. In a previous paper, we have performed a full analysis of the xSM parameter space and identified the features of the parameter space that can give detectable gravitational waves [13]. One of the most important results in the GW literature since our previous study is a very recent numerical simulation of GW production, which found a significant deficit of the produced GW signals from sound waves [14]. The suppression arises presumably from the slowing down of expanding broken phase bubbles due to formation of reheated droplets of unbroken phase. The GW production can be reduced by a factor as small as 0.001 and it has profound implications as all previous studies might have overestimated the GW signal strengths. Hence, in this paper we checked the impact of the findings of Ref. [14] on a part of arXiv:1909.05268v1 [hep-ph] 11 Sep 2019 the xSM parameter space and explore how it affects the complementarity of the future space-based GW experiments and the collider searches at the LHC. As a consequence of the above suppression of GW production, a detectable GW signal generation will require highly supercooled EWPT. Such supercooled EWPT may lead to vacuum energy dominated universe resulting in the universe stuck in a false vacuum. We address this issue with detailed analysis for our GW benchmarks in this paper, in contrast to our previous studies on GW [13,15]. We also present a hydrodynamics analysis of the fluid velocity profiles to determine bubble wall velocities that are consistent with EWBG, similar to our previous papers. This is important for the EWBG to generate the observed baryon asymmetry successfully.
In light of above changes, in this work we perform a dedicated study of the resolved 4b decay channel of the di-Higgs production to explore the parameter space giving large detectable GW. In previous papers, we have studied the GW complementarity in the h 2 → h 1 h 1 →bbγγ [15] and h 2 → V V (V = W, Z) [13] channels. In contrast, in this work, we analyze another potentially major channel that can discover h 2 at the LHC complementing the collider studies of our previous two papers mentioned above. Our estimate of the LHC sensitivity to this channel benefits from a reliable background estimation from the recent ATLAS search for double Higgs production and decay to 4b-tagged jets [16]. There has been a similar 4b analysis in the literature in connection with EWPT [17] using several benchmarks, which, however, does not include GW analysis. We find that the High Luminosity LHC (HL-LHC) measurements in the 4b channel will be able to probe the xSM parameters space predicting not too small h 2 → h 1 h 1 branching ratio while LISA can typically detect strong GW even with tiny h 2 → h 1 h 1 branching ratio thus establishing the complementary roles of colliders and satellite experiments in constraining models predicting first order EWPT.
The paper is organized as follows. We first give a brief introduction to the xSM in the Sec. II, followed by a description of the GW calculations in Sec. III. We present the dedicated 4b collider analysis in Sec. IV and conclude in Sec. V.

II. THE MODEL
The xSM model is defined by adding a gauge singlet real scalar to the SM with the following scalar sector potential [6,7,9]: Here is the SM Higgs doublet and S = v s + s is the additional singlet scalar. The parameters in this potential are all real. Of these parameters, two (µ, b 2 ) can be replaced by v s and v EW through the minimization conditions of the scalar fields; another three parameters (λ, a 1 , a 2 ) can be replaced by the masses and mixing angle of the physical scalars (m h 1 , m h 2 , θ). The physical scalars are defined by where h 1 is identified as the SM Higgs while h 2 is a heavier scalar. With this setup, the potential is fully specified by the following five unknown parameters: The parameter space defined by the five parameters above can be subjected to broadly two categories of constraints. The first set of constraints comes directly from various theoretical requirements imposed on the scalar potential, including boundedness of the potential from below, the stability of the EW vacuum, and perturbative unitarity of 2 → 2 scattering processes. All the other constraints are phenomenological. Higgs signal strength measurement [18] constrains the mixing angle θ: | sin θ| < 0.33 at 95% CL [18]. Another set of constraints comes from EW precision measurements such as the oblique S, T, U parameters [19,20] and correction to the the W boson mass m W [21]. Both the EW precision measurements mentioned above constrain only (m h 2 , θ) at one-loop level, with the m W measurement providing a more stringent bound [21,22]. For the details of this model and impact of the constraints used on the parameter space, we refer the reader to our previous paper [13].

III. GRAVITATIONAL WAVES
During the EWPT, a stochastic background of GW can be generated. In contrast to GW from a binary system, the amplitude of the stochastic background is a random variable, which is unpolarized, isotropic, and follows a Gaussian distribution [23]. Therefore, it is characterized by the two-point correlation function and is proportional to the power spectral density, Ω GW (f ). Due to its stochastic origin, the detection method is also different. With one detector, this signal would behave as another source of noise, making its identification difficult. Thus, the detection of this kind of GW depends on cross-correlating the outputs from two or more detectors (for recent reviews on cosmological sources of stochastic GW, see Ref. [24] and for detection methods, see Ref. [25]).
Given a particle physics model, the starting point of calculating the GW is the finite temperature effective potential. From this potential, a set of portal parameters, which characterize the dynamics of the EWPT may be calculated. Here, T n is the nucleation temperature and quantifies the time epoch when the bubbles are nucleated with a large probability; α is the energy density released from the EWPT normalized by the total radiation energy density at T n ; β/H n describes, approximately, the inverse time duration of the EWPT and also serves as a length scale for GW spectra, such as the peak frequency; v w is the bubble wall velocity; κ v is the fraction of released energy transferred into the kinetic energy of the plasma; and κ turb is the fraction of energy going to the Magneto-Hydrodynamic-Turbulence (MHD). It is this set of portal parameters that determine the GW spectra. Since we have given a detailed description of the calculations of these parameters in the xSM in Ref. [13], here we present only a summary of the key physical highlights adopted in our calculations, and the new features of the current study.
• v w and Hydrodynamics -A stronger GW production usually requires a larger v w , while a successful EWBG needs a small, generally subsonic v w (for EWBG calculations, see e.g., [26][27][28][29][30]). Hence, a large v w , which can source detectable GW signals, is detrimental to the process of generating the observed baryon asymmetry of the universe. The solution to this tension adopted in this work follows Ref. [31] and hinges on the recognition that v w may not be the quantity that enters EWBG calculations, due to the non-trivial plasma velocity profile surrounding the bubble wall. This can be further understood by noting that there exist three fluid velocity profiles for a single bubble in the phase transition: deflagration, supersonic deflagration (also known as hybrid), and detonation (see Ref. [32] for a recent combined analysis). In the cases of deflagration and supersonic deflagration, the fluid has a non-zero velocity outside the bubble wall. From the perspective of the bubble wall frame, the fluid would head towards the wall with a velocity (≡ v + ) that is smaller than v w . While The gray region in the top left is theoretically inaccessible and the gray region in the bottom right gives a detonation fluid profile which fails EWBG. The green region gives deflagration mode and the brown region gives hybrid mode. All (v w , α) on the curve labelled v + = 0.05 give a plasma relative velocity 0.05 heading towards the wall from outside the bubble. All the scanned points fall on this curve, with the red (blue) points giving SNR larger (smaller) than 50, which are calculated without including the suppression factor (see the text for the details of the suppression factor). The dots denote the pairs of (v w , α) for which numerical simulations of GW production are performed in Ref. [14], with the color of each point showing the value of the suppression factor of the GW signal as can be read from the color legend. The number close to the top-most point in each of the three columns of points denotes the corresponding suppression factor for that point.
the justification of this argument still needs a combined analysis of both macroscopic bubble behavior and microscopic particle transport dynamics, we assume tentatively that it is true in our series of papers on this subject [13,15,33,34]. This implies that in the plane of (v w , α), as shown in Fig. 1, only certain regions are compatible with EWBG. In this plot, the green and brown regions denote the deflagration and supersonic deflagration (hybrid) modes of the plasma, while the gray regions are either theoretically not allowed or not compatible with EWBG. The curve denotes the location of all the scanned points in this study, corresponding to v + = 0.05, the usually adopted benchmark for EWBG calculations.
• Efficiency factors κ v and κ turb -The energy fraction transferred to the fluid kinetic energy κ v is an important parameter as it directly determines the strength of the GW signal coming from the sound waves. In the past, it has been calculated as a function of (v w , α) by first solving the fluid velocity profile for a single bubble and then calculating the kinetic energy from the energy-momentum tensor [32,35]. The resulting value agrees well with values inferred from numerical simulations for relatively weak phase transition, defined such that α 1. This agreement had motivated the generalization of the obtained GW formula for arbitrary values of α.
This naive generalization, however, has been proven to be wrong by a very recent numerical simulation result [14]. In this study, several sets of simulations have been performed for strong phase transitions, defined by the authors as α ∼ 1. Their study shows a significant deficit in the produced GW signals when the plasma is of the deflagration mode while for detonations, it is less affected. There is currently no simulation result available for supersonic deflagrations.
The suppression for the deflagration plasma is found to be due to the reduction of κ v . The physical reason is that the formation of reheated droplets of the unbroken phase slows down the bubble expansion. This disagreement reveals a missing piece in the understanding of the fluid dynamics and requires a more accurate theoretical modeling in the future, which, however, is beyond the scope of our work. Therefore, we use the results from the limited set of numerical simulations given in [14] and quantify their effects for generic choices of (v w , α) 1 .
We denote the reduction of GW amplitude as δ and show its values for the set of simulations performed in Ref. [14] in the (v w , α) plot in Fig. 1. The three columns of dots denote the values of δ for the corresponding values of (v w , α), with the color characterizing its value, which can be read from the legend to the right of the plot. The number associated with the point at the top of each column explicitly denotes the corresponding value of δ. We can see that as v w is increased, δ is less reduced. Since simulation results covering the whole plane of (v w , α) are currently unavailable, we choose δ = 0.01 as a conservative estimate of the GW signals 2 .
For the energy fraction going to the MHD, we note that the current simulations are not long enough for MHD to fully develop. It is found that κ turb ≈ (0.05 ∼ 0.1)κ v [36]. We choose κ turb = 0.1κ v in our work.
• Supercooled EWPT -Strong GW signals generally require more energy released from the EWPT, corresponding to more supercooled phase transitions. This is primarily so due to the previously mentioned reduction in GW production, which makes generating a detectable GW signal more difficult and requires even more supercooled EWPT.
The environment within which the bubbles are nucleated may then be vacuum energy dominated rather than dominated by radiation. In this case, the space outside the bubbles will inflate, and it may happen that the bubbles never meet each other. The universe would be trapped in a metastable vacuum. Therefore there exists a maximal phase transition strength, above which the EWPT is not feasible [37]. To guarantee that EWPT does, in fact, complete, one must ensure the stronger condition that the physical volume (as opposed to the comoving volume) of the unbroken space is shrinking.
The temperature at which the above condition is imposed is chosen to be the percolation temperature in Ref. [37]. Here the percolation temperature (≡ T p ) is defined such that approximately 30% of the spatial volume is in the symmetry broken phase. We note that T p is typically slightly less than T n . In the xSM, we found that for most of the points we used in Ref. [13], T n is less than T p by 2%, and the stronger percolation criterion can be satisfied. line. The right panels shows the ρ V ≡ ∆V normalized by the total radiation energy density ρ R for this benchmark. Here the red vertical line corresponds to T p and the slightly higher (magenta dashed line) is obtained by simply neglecting ∆V in the Hubble expansion rate. Here "R" and "V" denote the contribution of radiation and vacuum energy density respectively.
As an example, we show the details of a super-strong EWPT with α = 7.43 in Fig. 2. The left panel shows S 3 (T )/T as a function of T and the positions of T n (cyan dashed) and T p (red), which are barely distinguishable from each other. The right panel shows ∆V , the difference between the potential energy density at the false and true vacua, normalized by the radiation energy density ρ R . Here aside from T p which was solved using the full Hubble expansion rate, we also present a higher temperature obtained by neglecting ∆V in H, to show the impact on the obtained percolation temperature with and without the vacuum energy. So for this parameter choice, when α is large, the two temperatures do not differ substantially. We numerically checked that for this benchmark the physical volume of the unbroken phase does shrink at T p .
• One or Two-Step EWPT -The EWPT in this model can proceed in two different patterns. The first pattern is a direct transition into the EW broken phase, while the second involves two stages, with the first one giving the extra singlet a non-zero vev followed by a second transition to the EW broken phase [13,[38][39][40]. Most of the parameter space in this model involves a single step transition, while regions yielding two-stage transitions reside in a totally different part of parameter space [13]. We thus focus on points sampled only from the parameter space of one-step EWPT.
The strategy described above is followed in calculating the set of parameters in Eq. 4. The GW spectra can be obtained by plugging these parameters into a set of analytical formulae, obtained by fitting to results from numerical simulations of different GW production mechanisms. It has been realized in recent years [41] that the dominant contribution comes from sound waves, although significant advances have also been made in both analytical modeling and numerical simulations of pure bubble collision contributions [42][43][44][45]. By evolving the scalar-field and fluid model on a 3-dimensional lattice, the gravitational wave energy density spectrum can be extracted [36]: Here g * is the number of relativistic degrees of freedom and H * is the Hubble parameter at T * , both evaluated at a time when the phase transition has just completed. Moreover, f sw is the present day peak frequency of this spectrum: Hz.
We have inserted an ad hoc factor δ in this formula, to take into account the reduction of the gravitational waves, as discussed earlier. We choose a conservative value of δ = 0.01, to be consistent with results of the sets of simulations performed in Ref. [14].
In addition, the fully ionized plasma can result in the formation of magnetohydrodynamic (MHD) turbulence. This can be calculated analytically by assuming a proper power spectrum for the turbulence as well as the primordial magnetic field [46][47][48][49] or numerically by evolving the magnetic field governed by the MHD differential equations and also coupled to gravity [50,51]. We use the result presented in [47,52], with f turb being the peak frequency: Hz.
With the energy density spectrum obtained and used as a matched filter, the outputs from a pair of gravitational wave detectors can be cross correlated to search for the signal. The noise of each detector drops out of the ensemble average (which is equivalent to a time average) of this cross correlation, leaving the desired signals 3 . The detectability of the gravitational waves is then quantified by the signal-to-noise ratio (SNR) [53]: where T is the duration of the time series of the detector ouput in years and Ω exp (f ) is a similarly defined detector power spectral density for this pair of interferometers. Several such detectors have been proposed, including the Laser Interferometer Space Antenna (LISA) [54], the Big Bang Observer (BBO), the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [55], Taiji [56] and Tianqin [57]. In this study, we will focus on LISA as a benchmark detector and assume T = 5. It was suggested in Ref. [53] that for a four-link LISA configuration, the threshold for detection is SNR = 50 and for a six-link configuration, an SNR as low as 10 can be used.

IV. COLLIDER ANALYSIS
Having discussed the production of GW, in this section, we present our collider analysis. Our focus is on double Higgs production at the 13 TeV HL-LHC. We access this production mode via the pp → h 1 (bb)h 1 (bb) channel, demanding four resolved b-tagged jets. The main backgrounds for this search are the QCD multi-jet andtt.
In previous papers, we have studied the GW complementarity in the h 2 → h 1 h 1 →bbγγ [15] and h 2 → V V (V = W, Z) [13] channels. In Ref. [15], some representative xSM benchmark points with favorably large SNR were investigated and the potential of the LHC to discover a double Higgs signal inbbγγ using machine learning techniques were explored. On the other hand, the analysis of Ref. [13] established the complementary role of the LHC searches of h 2 in diboson channel by merely extrapolating recent ATLAS and CMS results to 3 ab −1 . The 4b channel investigated in this work complements the results of those work for two reasons. First, h 1 h 1 and W W, ZZ are the main decay modes of h 2 in the considered parameter space. Second, bbbb is the double Higgs decay mode with the largest branching ratio.
In order to constrain the xSM in the GW context, we performed a shape analysis of the 4b mass distribution [58,59]. The signal samples are generated with MadGraph5 [60]. Hadronization and underlying event effects are accounted for with Pythia8 [61] and detector effects are simulated with Delphes3 package [62]. Higher order corrections are included with a next-to-next-to-leading order QCD K-factor [63][64][65] for the non-resonant component of the signal. The resonant contribution to the the signal is not expected to receive a very different QCD correction [66,67] so we keep the same K-factor to all contributions.
In our analysis, we closely follow the ATLAS study, Ref. [16]. ATLAS models the background with a data-driven approach. This is a more reliable procedure as the multi-jet component displays very large QCD corrections that are challenging to realistically account for in a Monte Carlo simulation. Hence, we use the backgrounds from ATLAS in our analysis. We validate our signal simulations using the SM double Higgs production presented by ATLAS [16]. We observe excellent agreement for the SM double Higgs mass distribution.
We start our analysis by requiring four isolated b-tagged jets with p T j > 40 GeV and |η j | < 2.5. Jets are defined with a cone radius of 0.4 using the anti-k t jet algorithm implemented in Fast-Jet [68,69]. The b-tagging requirements use the working point with a b-tagging efficiency of 70% associated to a mistag rate of 15% for c-quarks and 0.2% for light flavours.
The four b-tagged jets reconstruct the two SM Higgs boson candidates. The pairings of jets into Higgs boson candidates are required to satisfy: where ∆R jj,lead (∆R jj,subl ) is the angular distance between the jets that reconstruct the leading (sub-leading) SM Higgs boson candidate. In this first step of the analysis, the leading Higgs boson candidate is chosen to have the highest scalar sum of jet transverse momentum. To further reject the multi-jet background, we impose a pseudorapidity difference between the two Higgs candidates of |∆η hh | < 1.5.
Mass-dependent selections on the Higgs boson candidates transverse momenta are imposed to further control the backgrounds is computed and the pairing with the smallest D hh is chosen. These two pairs are then associated to m lead 2j and m subl 2j which are used to impose invariant mass selections around the SM Higgs boson mass for the leading and sub-leading Higgs boson candidates according to To further suppress thett background, all possible combinations of three jets with one being b-tagged and a constituent of the Higgs boson candidate are considered. The two light jets are considered as forming a hadronically decaying W boson candidate. A measure of the compatibility with the top-quark candidate can be defined as where m t is the invariant mass of the three jet top candidate and m W is the two jet W boson candidate. Events with the smallest X W t < 1.5, from all possible three jet combinations, are The left and right panels differ in that the points in the right panels used a reduction factor of δ = 0.01 in calculating the GW spectra while those in the left panel is obtained without this factor, i.e., δ = 1. Theoretical requirements on the potential such as the perturbative unitarity, perturbativity, and vacuum stability at zero temperature generally impose an upper bound on the branching ratio, corresponding to the green-colored region.
vetoed. To improve the signal m 4j resolution, the four-momentum of each Higgs boson candidate is multiplied by the correction factor m h 1 /m 2j . This was found to improve the signal mass resolution by 30% with sub-leading impact on the background m 4j distribution [16]. In Fig. 3, we display the double Higgs invariant mass m hh distribution for the Gluon Fusion (GF) contributions: SM, resonant xSM, and full xSM (that accounts for the resonant and non-resonant contributions). The total background component obtained from ATLAS with a datadriven approach is also shown [16]. The xSM distribution is illustrated with the parameter choice m h 2 = 800 GeV and sin θ = 0.2. We observe that the full xSM invariant mass m hh distribution display a significant contribution from the non-resonant GF terms. Hence, instead of only accounting for the resonant signal contribution, we describe the signal component as a deviation of the full xSM GF distribution from the SM GF one.
To quantify the sensitivity of the LHC towards the xSM model, we perform a one dimensional binned log-likelihood ratio analysis, exploring the m hh distribution. In Fig. 4, we display 2σ exclusion bound for the BR(h 2 → h 1 h 1 ) as a function of the heavy Higgs mass m h 2 . We assume the Higgs mixing sin θ = 0.2. The high luminosity LHC with 3 ab −1 will be sensitive to BR(h 2 → h 1 h 1 ) < 0.5 in the full mass range 300 GeV< m h 2 < 1000 GeV, using the 4b channel but a branching ratio as small as ∼ 0.05 can be probed for heavy Higgs masses around 500 GeV.
We also overlay on this plot the region allowed by the theoretical potential requirements 4 , obtained by a scan with sin θ = 0.2 for each m h 2 . Also shown here are points that can give a detectable GW signal at LISA, where green points give 10 < SNR < 50 and red points gives SNR > 50. To see the impact of the GW suppression on the SNR, we use δ = 1 in the left panel and δ = 0.01 for the right panel. The reduction of GW production from sound waves leads to a significant shrinking of the parameter space capable of generating detectable GW. However, the overall behavior of the remaining parameter space in affecting the branching ratio remain qualitatively unaffected. Compared with the theoretically allowed green region, this parameter space leads to an overall reduction of the di-Higgs branching ratio, especially for heavy h 2 . We can thus see clearly the complementary role played by colliders and GW detectors in probing the xSM. For lighter h 2 , i.e., 300 GeV m h 2 600 GeV, colliders will be able to probe, in the 4b channel alone, almost the entire parameter space that is capable of generating a detectable GW signal. For heavier h 2 , it is difficult for the HL-LHC to explore this regime because of phase space suppression. However, as the GW signal spectra with m h 2 > 600 GeV present very small BR(h 2 → h 1 h 1 ), a dedicated study of the h 2 → W W, ZZ channel to determine the potential of the LHC to probe those points, beyond that performed in Ref. [15], might be interesting.

V. SUMMARY
Gravitational waves from the EWPT provide a new window for probing the physics beyond the standard model, complementing the current direct collider searches. We continue this complementarity study in this work by focusing on the di-Higgs production in the 4b channel, choosing the benchmark model xSM. In calculating the gravitational wave spectra, we carefully accounted for several subtle issues, such as the bubble wall velocity, the supercooled phase transitions, and especially the reduction in the gravitational wave production from sound waves outlined in recently conducted numerical simulations [14]. These constitute important ingredients towards a faithful characterization of the GW signals from the EWPT and its detection at future gravitational wave detectors. The most important advance in recent understandings of these problems is the reduction of the GW produced from the sound waves, which invalidates the previous naive generalization of these formulae to arbitrary values of v w and α. We incorporated this effect by applying a conservative reduction factor of 0.01.
In order to establish the complementary role of collider searches, we performed an analysis in the resolved h 2 → h 1 h 1 → 4b channel. We found that the 13 TeV HL-LHC is able to probe the xSM parameter space with BR(h 2 → h 1 h 1 ) < 0.5 for 300 GeV < m h 2 < 1000 GeV. It is clear from our analysis that due to the significant reduction of the gravitational wave signal strength, the xSM parameter space, which is capable of giving a detectable stochastic GW background, have shrunk. However, the qualitative complementarity role of future space-based GW detectors in assisting BSM physics searches at current and HL-LHC remains unchanged.