Electroweak phase transition in the 2HDM: Collider and gravitational wave complementarity

The knowledge of the Higgs potential is crucial for understanding the origin of mass and the thermal history of our Universe. We show how collider measurements and observations of stochastic gravitational wave signals can complement each other to explore the multiform scalar potential in the two Higgs doublet model (2HDM). Accounting for theoretical and current experimental constraints, we analyze the key ingredients in the shape of the Higgs potential triggering the transmutation in phase transition, from the smooth crossover to the strong first-order phase transition ($\xi_c>1$), focusing on the barrier formation and the upliftment of the true vacuum. In particular, we observe that the $\xi_c>1$ regime is favored for lower scalar masses, rendering strong extra motivation for collider searches. We contrast the dominant collider signals at the HL-LHC (high-luminosity LHC) with observable gravitational wave signals at LISA. We obtain that the HL-LHC will be able to cover a vast range of the $\xi_c>1$ parameter space, with scalar decays to heavy fermions $(H,A,H^\pm\to tt, tb)$ being the most promising smoking gun signature of a strong first-order electroweak phase transition in the 2HDM.


I. INTRODUCTION
The structure of the Higgs potential is deeply connected with the origin of mass and thermal history of our Universe. When the Universe was cooling down, at a temperature of order 100 GeV, it went through a transition from a symmetric phase to an electroweak (EW) broken phase, where the Higgs field(s) acquired nonvanishing vacuum expectation values.
The evolution of this process strongly depends on the shape of the Higgs potential. Distinct profiles for the Higgs potential result in contrasting courses for the electroweak symmetry breaking (EWSB) in the early Universe, ranging from the smooth crossover transition in the Standard Model (SM), with the observed 125 GeV Higgs boson [1], to the strong first-order phase transition, with new physics.
The dynamics of the electroweak phase transition (EWPT) could have profound consequences for particle physics and cosmology. Most notably, it may be behind the matter and antimatter asymmetry puzzle. This asymmetry can be quantitatively featured by the baryon-to-photon ratio measurement n B /n γ ≈ 6 × 10 −10 [2], which is several orders of magnitude larger than the expected for the symmetric scenario, indicating an asymmetry in the early Universe between baryons and antibaryons. The bulk of the antibaryons have been annihilated in the thermal history, resulting in the large present density of baryons. The ingredients required to generate the baryon asymmetry of the Universe are theoretically well understood and summarized by the three Sakharov conditions [3]. They impose that our fundamental theory displays baryon number violation C along with CP violation, and departure from thermal equilibrium. Whereas the SM satisfies baryon number and C violation, the source of CP violation from the Cabibbo-Kobayashi-Maskawa (CKM) matrix is too small, and the observed Higgs mass m h = 125 GeV is too high to generate the out-ofequilibrium conditions from a strong first-order phase transition [1,4]. Thus, baryogenesis requires physics beyond the SM to grant these missing ingredients. In the present work, we focus on the latter problem, generating the out-of-equilibrium conditions at the electroweak scale [5][6][7][8].
The transmutation of the EWPT pattern from the smooth crossover to the strong firstorder phase transition usually requires new degrees of freedom around the EW scale, with sizable interactions with the Higgs boson [9]. Therefore, it generally accommodates beyondthe-SM scenarios with exciting phenomenological prospects, both at collider and gravita-tional wave (GW) experiments. At colliders, the Higgs pair production pp → hh usually plays a leading role in this discussion, as it grants a direct probe of the Higgs potential [10][11][12]. It provides access to the triple Higgs coupling via nonresonant Higgs pair production, as well as to the interactions of the SM Higgs boson with new heavy scalars through resonant di-Higgs searches. The high-luminosity LHC projections indicate that the triple Higgs coupling will be bounded to 0.1 < λ h 3 /λ SM h 3 < 2.3 at 95% confidence level (C.L.) [13]. Resonant searches are also a main focus, leading to significant limits [14]. Complementing the collider searches, the space-based GW experiments, such as LISA [15], will provide a new window to the Higgs potential. First-order phase transitions, that emerge from a scalar field tunneling through an energy barrier in the potential generate a significant source of gravitational waves. The correspondent signal spectrum displays a characteristic peak associated with the temperature at which the phase transition occurred. For phase transitions at the EW scale, this leads to a GW spectrum around the mHz frequency, after redshifting the signal to the present time [16,17]. This prompts exciting prospects to access the nature of EWPT at LISA, as this is precisely the frequency band that this experiment is sensitive to.
In this work, we study the EWPT pattern in the two Higgs doublet model (2HDM) [18][19][20][21][22][23][24][25][26][27][28], where the SM is augmented by an extra doublet. Instead of focusing on benchmarks or a particular parameter space region, a general scan is performed on the theoretically and experimentally allowed parameter space. We divide our analysis into two main stages. First, we scrutinize the new physics modifications to the shape of the Higgs potential that lead to a strong first-order phase transition. We devote particular attention to the barrier formation and to the true vacuum upliftment with respect to the SM case. The obtained results work as a guide for the phenomenological studies derived in the second part of the manuscript, where we perform the respective analysis at both the HL-LHC and LISA. In this last part, besides the commonly discussed channel A → ZH, other promising channels which can cover the EWPT parameter space are also investigated. These studies highlight the leading collider signatures for first-order EWPT at the LHC, as well as the complementarity of probes between collider and GW experiments.
The structure of this paper is as follows. In Section II, we briefly describe the 2HDM. The one-loop effective potential at finite temperature is discussed in Section III. It is followed by an introduction of EWPT and GW signals in Section IV. In Section V, we discuss how the shape of the potential will affect the EWPT, focusing on the barrier formation and vacuum upliftment. In Section VI, inspired by our shape analysis, we tailor the collider studies to the most promising channels. In addition, we derive the sensitivity to the correspondent GW signals generated in the early Universe. Finally, we summarize in Section VII. Some useful relations for the parameters in 2HDM are listed in Appendix A.

II. TWO HIGGS DOUBLET MODEL
The two Higgs doublet model displays one of the most minimalistic extensions of the SM that is compatible with the current experimental constraints [29]. In this work, we consider the CP -conserving 2HDM with a softly broken Z 2 symmetry. 1 The tree-level potential is given by where the mass terms m 2 11 , m 2 22 , and m 2 12 along with the couplings λ 1 ...λ 5 are real parameters from Hermiticity and CP -conservation. The required Z 2 symmetry, which is softly broken by m 2 12 , transformations Φ 1 → Φ 1 and Φ 2 → −Φ 2 guarantee the absence of the dangerous tree-level flavor changing neutral currents (FCNCs) [30,31]. After EWSB, the neutral components of the two SU (2) L doublets develop vacuum expectation values (VEVs). Expanding around the VEVsω i , the scalar doublets Φ i may be written as where the zero-temperature vacuum expectation values v i ≡ω i | T =0 are connected to the SM The CP -conserving 2HDM leads to five physical mass eigenstates in the scalar sector: two CP -even neutral scalars h and H, a neutral CP -odd scalar A, and a charged scalar pair H ± . The relation between the mass and gauge eigenstates is established by the rotation angle β for the charged and CP -odd sectors, where tan β ≡ v 2 /v 1 , and by the mixing angle 1 Baryogenesis requires physics beyond the SM to generate new sources of CP violation and out-ofequilibrium conditions. In the present work, we focus on the latter issue. α in the CP -even sector  The rotation matrix is defined as with s x ≡ sin x and c x ≡ cos x. G ± and G 0 represent the charged and neutral massless Goldstone bosons.
Instead of the eight parameters in the Higgs potential m 2 11 , m 2 22 , m 2 12 , λ 1 ...λ 5 , a more convenient choice of parameters is The conversion between these two sets of parameters can be found in Appendix A. The parameters t β ≡ tan β and c β−α ≡ cos(β − α) are of critical phenomenological importance.
They control the coupling strength of scalar particles to fermions and gauge bosons. Given the current experimental constraints, a particular relevant regime is the alignment limit c β−α = 0 [32], where the 125 GeV CP -even scalar Higgs boson couples to SM particles precisely as the SM Higgs boson.
In general, there are four types of Z 2 charge assignments in the Yukawa sector that avoid FCNC at tree level. In this work, we focus on the type-I and type-II scenarios. In the first case, all fermions couple only to Φ 2 , whereas in the latter, only the up quarks couple with Φ 2 , leaving the down quarks and charged leptons to couple with Φ 1 . For both types I and II, we perform a uniform scan over the parameter space region, tan β ∈ (0.8, 25) , m 2 12 ∈ (10 −3 , 10 5 ) GeV 2 , m H ∈ (150, 1500) GeV , The observed 125 GeV Higgs boson is identified with the h scalar. The parameter space scan is performed with ScannerS [33,34]. Using this framework, we impose the constraints from perturbative unitarity [35][36][37], boundedness from below [38], vacuum stability [39,40], electroweak precision, and flavor constraints. In addition, HiggsBounds and HiggsSignals are used to incorporate the searches for additional scalars as well as the constraints from the 125 GeV Higgs boson measurements [41,42].

III. ONE-LOOP EFFECTIVE POTENTIAL AT FINITE TEMPERATURE
To study the electroweak phase transition in the early Universe, we use the loopcorrected effective potential at finite temperature. In addition to the tree-level potential V 0 from Eq. (1), the effective potential displays one-loop corrections at zero temperature from the Coleman-Weinberg potential V CW and counterterms V CT . Finite-temperature corrections V T are also included. The effective potential reads The Coleman-Weinberg potential can be written in the Landau gauge as [43] where the index i sums over all particles in the thermal bath with field-dependent mass In general, the one-loop Coleman-Weinberg corrections shift the scalar masses and mixing angles with respect to the tree-level potential. To optimize our parameter scan, we adopt a renormalization prescription that enforces these parameters to match with their tree-level values [19,46]. In this setup, the counterterm part of the potential can be written as with the following on-shell renormalization conditions at zero temperature: The fields φ i (i = 1, ..., 8) denote scalar components from the Φ 1 and Φ 2 doublets, ω generically represents the ω i values, and ω tree generically stands for the minimum of the tree-level potential for the φ i fields. We followed the prescription of Ref. [46] to consistently calculate the first and second derivatives of V CW . The first renormalization condition Eq. (10) imposes that the zero-temperature minimum is not shifted with respect to the tree-level value. Similarly, the second condition Eq. (11) ensures that the zero-temperature masses and mixing angles remain the same as their tree-level assignment.
The last term in Eq. (7), the one-loop thermal corrections V T , can be expressed as [47] where the sum extends over fermions f and bosons, with the latter subdivided in transverse modes of gauge boson V T = W T , Z T and longitudinal modes of gauge bosons and scalars The resummation of the n = 0 Matsubara modes of V L result in thermal corrections to their masses [48,49]. The second line in Eq. (12) indicates the Daisy contributions, where m V L is the thermal Debye mass following the Arnold-Espinosa scheme [19,47]. Lastly, the thermal functions for fermions (J + ) and for bosons (J − ) read While the effective potential displays theoretical uncertainties arising from the choice of gauge parameter [44,[50][51][52][53][54], gauge-independent probes can be constructed exploiting the Nielsen identities [55]. These identities state that the gauge dependence vanishes at the extrema of the potential where ξ is the gauge fixing parameter. This motivates us to adopt two distinct methods in our manuscript for the phenomenological analyses. The first approach encompasses the calculation of the finite-temperature effective potential and the subsequent numerical scan.
The second one focuses on the calculation of the vacuum upliftment at T = 0. As we will highlight in Sec. V B, the upliftment of the true vacuum with respect to the symmetric one at zero temperature works as an effective probe to the strength of phase transition. Whereas the first method displays uncertainties rooted in the choice of gauge parameter, the latter approach is gauge invariant, as guaranteed by the Nielsen identities [22,50]. Notice that we introduce extra counterterms to preserve the position of the EW vacuum, as well as the masses, at one-loop order. The phenomenological agreement between our numerical scan with the profile derived from the vacuum upliftment will, in particular, evince that the numerical scan is well grounded, despite its uncertainties. 3

IV. ELECTROWEAK PHASE TRANSITION AND GRAVITATIONAL WAVES
The effective potential at finite temperature determines the dynamics of the phase transition. The two Higgs doublet model exhibits multiple phase transition processes. For successful baryogenesis, the sphaleron process inside the bubble should be heavily suppressed to prevent the net baryon number generated around the bubble wall from significant washout.
This condition requires that the EWPT be of strong first order [49] where v c ≡ ω 2 1 + ω 2 2 | Tc is the Higgs VEV at the critical temperature T c , which is defined when the would-be true vacuum and false vacuum are degenerate. The approximate inequality denotes the theoretical uncertainty in this condition [50].
The transition from false to true vacuum takes place via thermal tunneling. It results in the formation of bubbles of the broken phase that expand in the surrounding region of symmetric phase, converting the false vacuum into true vacuum. The tunneling probability can be written as [56,57] 3 Note that when the coupling in the scalar potential is large, as usually required by a strong first-order phase transition, one should check the reliability of the perturbative calculations. Lattice simulations for 2HDM are performed for two particular benchmark points in Ref. [24], where the authors made comparisons among different methods. It is shown that the perturbative estimation of the strength of the phase transition where S 3 is the three-dimensional Euclidean action corresponding to the critical bubble Here, the scalar field φ is the bubble profile of the critical bubble. It is obtained as a solution to the following differential equation We use the publicly available code CosmoTransitions to solve the differential equation and compute the Euclidean action S 3 [58].
The first-order phase transition completes around nucleation temperature T n , where one bubble nucleates per unit volume [59] ∞ Tn dT T This ensures that the bubbles percolate even in the inflating universe. For the electroweak phase transition, the above condition can be roughly approximated as [49] One of the important consequences of the strong first-order EWPT is the production of stochastic gravitational waves. The GW signals from phase transition have three main sources: collision of the vacuum bubbles, sound waves, and turbulence in the plasma. For each source, the GW spectrum can be expressed as numerical functions in terms of two important parameters determined from the phase transition dynamics [16,60]. The first parameter is α ≡ /ρ rad , the latent heat released in the phase transition ( ) to the radiation energy density (ρ rad ). The latent heat and ρ rad are obtained from where ∆ means the difference between the stable and metastable minima, and g is the number of relativistic degrees of freedom in the plasma. The second key parameter characterizing the spectrum of gravitational waves is the inverse time duration of the phase transition β/H n . This quantity is defined as where H n denotes the Hubble constant at the nucleation temperature T n . Strong GW signals are typically associated with large latent heat release (large α) and slow phase transition Finally, to estimate the sensitivity of GW experiments, we adopt the signal-to-noise ratio where Ω Sens is the sensitive curve of the considered GW experiment [15] and T is associated with the duration of the mission. In the present study, we focus on the LISA experiment as a benchmark, assuming T = 5 years and the threshold of detection as SNR = 10 [60].

V. THE SHAPE OF THE HIGGS POTENTIAL
The extension of the SM Higgs sector with another Higgs doublet in the 2HDM can promote the phase transition pattern from a smooth crossover to a strong first-order phase transition. 4 To study the key ingredients triggering this transmutation in the EWPT, we need to probe the shape and thermal evolution of the effective potential. In this section, we focus on the barrier formation and the upliftment of the true vacuum [22,61,62]. The parameter space analysis is organized based on the most relevant components of the effective potential generating these changes in the profile of the Higgs potential [63]. This characterization is used as a key ingredient to pin down the leading phenomenological parameters for strong first-order phase transition.
In our studies, we focus on the type-I and -II scenarios for the Yukawa couplings. The main difference between these two cases in the effective potential comes from the bottom Yukawa coupling in the Daisy terms. In terms of phase transition, these two scenarios result in negligible differences. In Fig. 1, we show the ratio of ξ c between type I and II as a function of t β using the same input parameters. We focus on the points that satisfy all the current constraints for both scenarios. We observe depleted differences between ξ I c and ξ II c with most points differing only in the subpercent level. While there is a small enhancement for the ratio ξ I c /ξ II c toward larger t β (as in type II, the Daisy contributions result in a slightly deeper potential at the true vacuum), the difference is phenomenologically insignificant. Since the considered scenarios display a similar phase transition profile, with only subleading differences, we will mostly focus on the type-I case in the present section.
Nevertheless, when we discuss the experimental sensitivities, we will show the results for both scenarios as they can present distinct collider phenomenology due to their different fermionic couplings.

A. Barrier formation
Moving forward, we scrutinize the 2HDM phase transition pattern, analyzing three classes of contributions to the potential barrier: tree-level (V 0 ), one-loop (V 1 ≡ V CW + V CT ), and thermal effects (V T ). Our main target is to identify which of these terms plays a crucial role in introducing the barrier between the broken and unbroken vacua, granting the possibility of strong first-order electroweak phase transition (SFOEWPT) in the 2HDM [63]. The correlations among these contributions to the potential barrier are presented in Fig. 2 illustration, we focus on the type-I 2HDM with ξ c > 1. The barrier of the potential is the position where the effective potential obtains the maximum value in the tunneling path obtained by solving Eq. (18). We defined the height of the barrier as the difference between the effective potential at the barrier and false vacuum. The position of the barrier at T c is approximated by the point where the potential attains maximum value in the line connecting the true and false vacua. We observe that the potential barrier in the 2HDM is dominantly generated by a coalition between the one-loop and thermal components. These terms display positive contributions to the potential barrier δV b 1 , δV b T > 0 for 99% of the parameter space points. In contrast, the tree-level term V 0 typically works against the barrier formation.

. For
In Fig. 3 (left panel), we present the fraction between the two leading terms to the potential barrier as a function of the order parameter ξ c . We focus on the region with positive contributions δV b 1 , δV b T > 0 enclosing the bulk of the parameter space points. Two comments are in order. First, in the strong first-order phase transition regime ξ c > 1, the phase transition is mostly one-loop driven; i.e., the effective potential barrier is dominantly generated by the one-loop term. In this case, the loop corrections can generate relevant nonpolynomial field dependencies, such as h 4 ln h 2 , that contribute to the barrier formation [63,64].
Second, if the fraction of the barrier height provided by the one-loop contribution is close to 100%, the tunneling from the false vacuum (metastable vacuum) to the true vacuum is more challenging. For this reason, the universe with ξ c 2.5 is trapped in the false vacuum and electroweak symmetry breaking does not occur. We should notice that this feature is associated with the dominant phase space regime δV b 1 , δV b T > 0. Conversely, the rarer tree-level or thermally driven setup can still generate stronger phase transition ξ c 2.5.

B. Vacuum upliftment
After looking at the general new physics contributions producing the barrier, we now focus on the effects on the potential at the vacua. It has been shown that the strength of the phase transition is correlated with the upliftment of the true vacuum compared to the symmetric one at zero temperature [22,61,62]. That is, if the Higgs potential is shallow at T = 0, the required thermal upliftment for SFOEWPT, making the true vacuum degenerate with the false one, is reduced. Following a similar notation to Ref. [22], we define a dimensionless parameter to measure the true vacuum upliftment where F 0 is the vacuum energy density of the 2HDM at T = 0 defined as   Fig. 3 (right panel), we note that the barrier height provided by the one-loop contribution is correlated with ∆F 0 /|F SM 0 |, which measures the vacuum upliftment at zero temperature. The larger the one-loop contribution, the higher the vacuum upliftment. This correlation is especially prominent for ξ c > 1. Since the one-loop effects are dominant with respect to thermal corrections for ξ c > 1, ∆F 0 /|F SM 0 | works as a good first approximation to study some general properties of the EWPT, even though it is a zero-temperature quantity. In particular, it is possible to define a typically necessary condition for a first-order phase transition with the minimum threshold ∆F 0 /|F SM 0 | > 0.34 [22]. Although this condition encapsulates most of the ξ c > 1 points, we should note from Fig. 3 that it does not work as a sufficient requirement for SFOEWPT.
While ∆F 0 /|F SM 0 | does not provide a one-to-one correlation with the strength of EWPT, it is helpful to explain critical features of the 2HDM parameter space. In the alignment limit, the important contributions to F 0 come from H, A, and H ± which can be written as Using these expressions, we can write the shift in vacuum energy density with respect to the SM value with Remarkably, the individual contributions from H, A, and H ± to F 0 are of the same form.
Thus, in this sense, there should be no difference in the preferred region in terms of m H , m A , and m H ± . On the other hand, m 2 12 /(s β c β ) also plays an important role. In Fig. 4, the individual contributions from A to F 0 are shown (H and H ± have the same form), from which it is easy to find that F 1,Φ 0 will be negative when m 12 / √ s β c β is larger than the scalar mass m Φ , where Φ = H, A, or H ± . The larger the difference is, the more negative it will be. Contrarily, when the scalar mass is larger than m 12 / √ s β c β , F 0 will tend to be positive, uplifting the true vacuum, and favoring the strong first-order phase transition. However, the vacuum upliftment ∆F 0 is limited from above and below by perturbative unitarity constraints [37].
For illustration, as a benchmark, we denote the perturbative unitarity constraints as the gray shaded region in Fig. 4 for t β = 1, c β−α = 0, and m A = m H ± = m H + 100 GeV.
As the allowed region (unshaded area) becomes narrower toward larger scalar masses and Before discussing the phase transition pattern presented in Fig. 5, we would like to point out several theoretical constraints that will be important for our analysis. Especially, we want to highlight that m 12 / √ s β c β cannot be too different from the scalar masses. First, we consider the perturbative constraints. We start by writing the λ 1 and λ 2 couplings in the alignment limit Because of their strong t β dependence, perturbativity limits for λ 1 and λ 2 demand m 2 H ≈ m 2 12 /(s β c β ) for t β significantly different from 1. Second, the boundedness from below limits [38]  • m H < m H ± ≈ m A (left branch): This mass configuration displays numerous first-order phase transition points. Since m 2 12 /(s β c β ) cannot be significantly larger than m 2 H , this leads to large and positive contributions from A and H ± to F 0 associated with typically small effects from H. These properties promote the m H < m H ± ≈ m A regime as one of the most likely configurations to achieve SFOEWPT in the 2HDM. Further, when the mass difference |∆m H | is large, the contributions to F 0 from H ± and A are also sizable. Hence, large |∆m H | is more likely to grant SFOEWPT. 5 Remarkably, this results in important phenomenological consequences for LHC searches. In particular, it favors new physics searches via the A → ZH channel that will be discussed in Sec. VI B.
• m H ≈ m H ± < m A (top branch): This regime also presents a sizable number of ξ c > 1 points. As m 2 12 /(s β c β ) cannot be much larger than m 2 H , the leading positive contribution to F 0 arises from the pseudoscalar A. In addition, when m 2 12 /(s β c β ) is smaller than m 2 H(H ± ) , H (H ± ) can also contribute positively, while at a subleading level when compared to A. Hence, the order parameter ξ c tends to be slightly suppressed, in comparison to the left branch. Finally, due to similar arguments as for the left branch, sizable |∆m A | is more likely to yield SFOEWPT.
• m H ≈ m H ± > m A (bottom branch): This region can provide SFOEWPT, as long as m 2 12 /(s β c β ) is much lower than m 2 H ≈ m 2 H ± . In this regime, the constraints from λ 1,2 imply t β → 1 to generate a large order-parameter, ξ c > 1.
• m H > m H ± ≈ m A (right branch): This mass configuration renders a suppressed number of first-order phase transition points. Similar to the bottom branch, m 2 12 /s β c β has to be much lower than m 2 H to achieve the ξ c > 1 regime. However, in this parameter space, only H can contribute significantly to F 0 , which leads to a lower chance to achieve SFOEWPT.  [67].

VI. COLLIDER AND GRAVITATIONAL WAVE SIGNALS
In this section, we focus on the complementarities between collider and gravitational wave experiments to probe the phase transition pattern in the early Universe. In Fig. 6 It includes a tabulated parameterization of the next-to-next-to-leading-order (NNLO) QCD gluon fusion and bb-associated Higgs production obtained from SusHi [69,70]. It also encompasses the next-to-leading-order (NLO) QCD top quark and charged Higgs boson associated production parametrized within HiggsBounds [41,[71][72][73][74][75][76].

A. Resonant and nonResonant di-Higgs searches
While in the alignment limit, the tree-level Higgs self-coupling matches the SM value λ h 3 = λ SM h 3 , the one-loop corrections in the 2HDM can significantly disrupt this equality. In Fig. 7, we show the new physics effects on the triple Higgs coupling λ h 3 /λ SM h 3 at one loop. We observe that the higher-order effects can produce extremely high deviations from the SM, as large as λ h 3 /λ SM h 3 ≈ 7, even in the alignment limit and in view of the theoretical and on the decoupling of the heavy states [78]. In view of the preference for relative light scalar modes m H 750 GeV, producing sizable order parameter ξ c > 1, the EFT does not provide a robust framework to study the preferred SFOEWPT parameter space regime for the 2HDM at the LHC [79]. Hence, as our limits strongly depend on the LHC results, we refrain from using the EFT approach in this study.
The Higgs pair production pp → hh provides a direct probe for the Higgs self-coupling at colliders [11,80,81]. However, the limited production rate, large destructive interference between the triangle and box diagrams, and sizable backgrounds make this analysis extremely challenging at the LHC. The ATLAS and CMS high-luminosity projections constrain the triple Higgs coupling at 95% C.L. to [13] 0.
This limited precision prompts the Higgs self-coupling as a key benchmark for future colliders. In particular, the rapid increase of the gluon luminosity at higher energies translates in a sizable pp → hh cross section at the 100 TeV Future Circular Collider. In such a setup, the large number of signal events would transform the di-Higgs production into a precision measurement, allowing for the full kinematic exploration, that is central for a better resolution of λ h 3 [82]. The improvement to the Higgs self-coupling measurement at higher energy colliders would allow for a more meaningful global fit analysis that is sensitive to a more complete set of new physics modifications to the Higgs potential [83]. Despite the limited experimental constraints, we observe in Fig. 7 that the HL-LHC will be sensitive to a large range of the 2HDM parameter space. In particular, it will be able to probe a substantial fraction of points with a large order-parameter that generally correlates with sizable triple Higgs coupling.
While the HL-LHC will be mostly sensitive to SFOEWPT with large order-parameter ξ c 2.5, LISA will be broadly sensitive to GW signals in the complementary regime ξ c 2.5, as shown by the color points of Fig. 7 (right two panels). As explained in Section V, extremely large ξ c leads to the configuration where the system is trapped in the false vacuum precluding successful nucleation. This is reflected in Fig. 7, where the points with large λ h 3 (and thus large ξ c ) in the left two panels do not appear in the right two panels.
Resonant di-Higgs searches provide another prominent probe for the phase transition pattern in the early Universe. As discussed in the last section, SFOEWPT is usually associated with light extra scalars m H 750 GeV resulting in a favored energy range for pp → H → hh production at the LHC. In Fig. 8, in addition to the current theoretical and experimental limits, we present the projected HL-LHC sensitivity at 95% C.L. to the resonant di-Higgs cross section (red dashed line). The results are obtained scaling the current sensitivity presented by ATLAS in Ref. [84], according to the luminosity, to the high-luminosity LHC with L = 3 ab −1 . This experimental study focuses on the leading 4b final state channel.
While this analysis explores the gluon fusion production mode, we note that the resonant production through weak boson fusion can provide additional relevant extra sensitivity [85].
We observe in Fig. 8 that the projected resonant Higgs pair production measurement will In the type-II scenario, the latter channel does not provide extra sensitivity at the HL-LHC to SFOEWPT. This is due to the stronger constraints on c β−α in the type-II 2HDM, pushing it further toward the alignment limit (see, e.g. Fig. 6 or Fig. 7). Remarkably, even exploring the dominant channel, where we have A → ZH and H → bb, the sensitivity is still somewhat weak, being limited mostly to the parameter space region m H 350 GeV. The main reason is that the bottom quark decay channel quickly becomes subdominant once the scalar mass is beyond the top-pair threshold. 7 Therefore, the fermionic channel with the top quark will be more promising in the high-mass region. We will explore the heavy scalar decays to top pairs in the next subsection.
The flipped channel H → ZA can a priori also provide strong limits [88]. It corresponds to the mass regime m H > m A , which is associated with the right and bottom branches of Fig. 5. In the right branch, the SFOEWPT is suppressed due to limited positive contributions from heavy scalars to F 0 . Conversely, the bottom branch could still provide ξ c > 1 points with t β → 1, as discussed in the previous section. However, this t β regime precludes   Ref. [91] according to the luminosity to the HL-LHC with 3 ab −1 . We observe that this channel is capable of covering the relevant region of the parameter space that can trigger strong first-order EWPT and produce detectable gravitational wave signals.

D. Combined results
Here, we compare the sensitivity to SFOEWPT for the three aforementioned search channel categories: is not yet performed by ATLAS and CMS. We leave the referent phenomenological study for future work [92].
In Fig. 13, we present the global HL-LHC and GW complementarities to probe the currently available SFOEWPT parameter space based on our uniformly random scan. We see that the HL-LHC searches will be able to cover ≈ 80% of the remaining ξ c > 1 parameter space for the type-I 2HDM and an impressive ≈ 96% for the type-II scenario. At the same time, LISA will be able to access a complementary parameter space region with a typically low production cross section at the HL-LHC for the considered processes. The requirement for small scalar masses to induce positive contributions for F 0 plays a crucial role in the sizable HL-LHC sensitivity to SFOEWPT. These fractions present only a lower bound. Adding other complementary 2HDM search channels at the HL-LHC, beyond the three considered classes, should push the quoted sensitivities to an even higher level.

VII. SUMMARY
Reconstructing the shape of the Higgs potential is crucial for understanding the origin of mass and the thermal history of electroweak symmetry breaking in our Universe. In this work, we explore the complementarity between collider and gravitational wave experiments to probe the scalar potential in the 2HDM. We scrutinize fundamental ingredients in the profile of the Higgs potential, namely the barrier formation and upliftment of the true vacuum that promote the transmutation of phase transition from the smooth crossover to the strong first-order phase transition. In addition, accounting for the theoretical and current experimental measurements, we study the prospects for the HL-LHC to probe the ξ c > 1 regime focusing on three prominent classes of searches [resonant and nonresonant di-Higgs, A(H) → ZH(A), and heavy scalar decays to fermions] and contrasted with the GW sensitivity at LISA. We summarize our novel results as follows: • When comparing the parameter space points that survive the theoretical and experimental constraints for type-I and type-II 2HDM, these scenarios result in an akin phase transition pattern.
• The barrier formation in the Higgs potential of the 2HDM is driven by the one-loop and thermal corrections, with the dominance of the one-loop terms for large orderparameter ξ c > 1.
• The strength of phase transition is correlated with the upliftment of the true vacuum with respect to the symmetric one at zero temperature [22,61,62]. This arises as a result of the dominance of the one-loop effects with respect to the thermal corrections for ξ c > 1. Based on this result, we shed light on the phase transition pattern analytically. In particular, we observe that larger vacuum upliftment is favored for lower scalar masses which is in accordance with the results from a generic discussion in [9] . This provides strong extra motivation for scalar searches at the LHC. Besides scalar masses below the TeV scale, the analytical structure of the new physics effects on the vacuum upliftment, leading to SFOEWPT, result in a peculiar hierarchy of masses among the new scalar modes. These findings work as a guide for collider and gravitational wave studies.
• We obtain that the scalar decays to heavy fermions (H, A, H ± → tt, tb) are the most promising smoking gun signature for SFOEWPT at the HL-LHC, followed by the di-Higgs searches. Based on the projections from the current ATLAS and CMS searches, the widely discussed channel A(H) → ZH(A) is still relevant, whereas to a smaller fraction of the parameter space. The main reason for such an observation is that the current experiments focus on the bb and W W decay channels [87]. These two decay modes only cover a small portion of the parameter space. We leave for future work a direct phenomenological comparison of the gluon fusion gg → H(A) and A(H) → ZH(A) channels, considering the promising tt heavy fermion final states [92].
• In contrast to the HL-LHC, LISA is going to be sensitive to a significantly smaller parameter space region, whereas it renders to complementary sensitivities where the correspondent LHC cross section is suppressed. Based on our parameter space scan, the combination of the LHC searches with gravitational wave studies presents exciting prospects to probe the vast majority of first-order phase transition points in the 2HDM.
Adding other complementary 2HDM search channels at the HL-LHC, beyond the three considered classes, should push new physics sensitivity to an even higher level.
In conclusion, the study of the thermal history of electroweak symmetry breaking is a crucial challenge for particle physics and cosmology. We demonstrate that the well-motivated 2HDM leads to a rich phase transition pattern favoring SFOEWPT below the TeV scale.
This renders exciting physics prospects at the HL-LHC and upcoming gravitational wave experiments, such as LISA.