# First shot of the smoking gun: probing the electroweak phase transition in the 2HDM with novel searches for

$A \to ZH$ in $\ell^+\ell^- t\bar{t}$ and $\nu\bar{\nu} b\bar{b}$ final states

Thomas Biekötter, Sven Heinemeyer, Jose Miguel No, Kateryna Radchenko, María Olalla Olea Romacho, Georg Weiglein

► To cite this version:

Thomas Biekötter, Sven Heinemeyer, Jose Miguel No, Kateryna Radchenko, María Olalla Olea Romacho, et al.. First shot of the smoking gun: probing the electroweak phase transition in the 2HDM with novel searches for $A \to ZH$ in $\ell^+\ell^- t\bar{t}$ and $\nu\bar{\nu} b\bar{b}$ final states. Journal of High Energy Physics, 2024, 2024 (01), pp.107. 10.1007/JHEP01(2024)107. hal-04246821

HAL Id: hal-04246821

https://hal.science/hal-04246821v1

Submitted on 9 Feb 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# First shot of the smoking gun: probing the electroweak phase transition in the 2HDM with novel searches for $A \to ZH$ in $\ell^+\ell^- t\bar{t}$ and $\nu\bar{\nu} b\bar{b}$ final states

Thomas Biekötter,ª Sven Heinemeyer,ᵇ Jose Miguel No,ᵇ,c Kateryna Radchenko,ᵈ
María Olalla Olea Romachoe and Georg Weigleinf

ªInstitute for Theoretical Physics, Karlsruhe Institute of Technology,
Wolfgang-Gaede-Str. 1, 76131 Karlsruhe, Germany

ᵇInstituto de Física Teórica UAM-CSIC,
Cantoblanco, 28049, Madrid, Spain

ᶜDepartamento de Física Teórica, Universidad Autónoma de Madrid (UAM),
Campus de Cantoblanco, 28049 Madrid, Spain

ᵈDeutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany

ᵉLaboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité,
F-75005 Paris, France

fII. Institut für Theoretische Physik, Universität Hamburg,
Luruper Chaussee 149, 22607 Hamburg, Germany

E-mail: thomas.biekoetter@kit.edu, sven.heinemeyer@cern.ch,
josemiguel.no@uam.es, kateryna.radchenko@desy.de,
mariaolalla.olearomacho@phys.ens.fr, georg.weiglein@desy.de

**ABSTRACT:** Recently the ATLAS collaboration has reported the first results of searches for heavy scalar resonances decaying into a $Z$ boson and a lighter new scalar resonance, where the $Z$ boson decays leptonically and the lighter scalar decays into a top-quark pair, giving rise to $\ell^+\ell^- t\bar{t}$ final states. This had previously been identified as a smoking-gun signature at the LHC for a first-order electroweak phase transition (FOEWPT) within the framework of two Higgs doublet models (2HDMs). In addition, ATLAS also presented new limits where the $Z$ boson decays into pairs of neutrinos and the lighter scalar resonance into bottom-quark pairs, giving rise to the $\nu\bar{\nu} b\bar{b}$ final state. We analyze the impact of these new searches on the 2HDM parameter space, with emphasis on their capability to probe currently allowed 2HDM regions featuring a strong FOEWPT. We also study the complementarity of these new searches with other LHC probes that could target the FOEWPT region of the 2HDM. Remarkably, the ATLAS search in the $\ell^+\ell^- t\bar{t}$ final state shows a local $2.85\sigma$ excess (for masses of about 650 GeV and 450 GeV for the heavy and light resonance) in the 2HDM parameter region that would yield a FOEWPT in the early universe, which could constitute the first experimental
hint of baryogenesis at the electroweak scale. We analyze the implications of this excess,
and discuss the detectability prospects for the associated gravitational wave signal from
the FOEWPT. Furthermore, we project the sensitivity reach of the $ℓ^{+}ℓ^{-}t\bar{t}$ signature for
the upcoming runs of the LHC. Finally, we introduce the python package thdmTools, a
state-of-art tool for the exploration of the 2HDM.

KEYWORDS: Multi-Higgs Models, Specific BSM Phenomenology

arXiv ePRINT: 2309.17431
# Contents

1 Introduction 1

2 Theoretical background 5

2.1 The Two Higgs Doublet Model (2HDM) 5

2.2 Thermal history analysis 6

2.3 Vacuum tunneling 7

2.4 Strong first-order electroweak phase transitions and baryogenesis 7

2.5 Gravitational waves 8

3 Numerical analysis 9

3.1 Implementation of new limits into HiggsBounds 9

3.2 New constraints on the 2HDM parameter space 10

3.2.1 Low tan $\beta$-region 11

3.2.2 High tan $\beta$-region 17

3.3 Future prospects for $\ell^+\ell^- t\bar{t}$ searches 18

3.4 A hint of a strong 1$^{st}$-order EW phase transition in the 2HDM? 21

3.4.1 Preferred parameter regions 22

3.4.2 Gravitational wave detection 24

4 Summary and conclusions 27

A The python package thdmTools 30

A.1 Installation 31

A.2 Code example for thdmTools 31

B Comparison to previous CMS projections 34

## 1 Introduction

The Standard Model (SM) of particle physics predicts the existence of a fundamental scalar particle as a consequence of the Higgs mechanism. In the year 2012 a scalar particle with a mass of about 125 GeV was discovered at the Large Hadron Collider (LHC) [1, 2]. So far, the properties of the detected scalar are compatible with the predictions of the SM. However, the experimental precision of coupling measurements of the detected Higgs boson is still at the level of 10% at best [3, 4], such that there is ample room for interpretations of the detected Higgs boson in various theories Beyond the Standard Model (BSM).

The SM predictions have been tested at the LHC at various energy scales, and so far the measurements are in remarkable agreement with the SM [5]. On the other hand, the SM fails to address various major existing observations in Nature. One of the most pressing open questions in this context concerns the origin of the matter-antimatter asymmetry of the Universe, which (according to the measured value of the mass of the Higgs boson) cannot
be explained in the SM [6]. BSM theories can address the shortcomings of the SM. In particular, models featuring extended Higgs sectors could allow for the generation of the baryon asymmetry of the Universe (BAU) via electroweak (EW) baryogenesis [7]. One of the simplest constructions to realize EW baryogenesis is based on a Higgs sector containing a second Higgs doublet [8–10], i.e. the Two Higgs Doublet Model (2HDM). This model contains two CP-even Higgs bosons, $h$ and $H$, one CP-odd Higgs boson, $A$ and a pair of charged Higgs bosons, $H^{\pm}$. In the 2HDM the phase transition giving rise to EW symmetry breaking can be rendered to be a sufficiently strong first-order transition, providing the out-of-equilibrium conditions required for EW baryogenesis [11–14]. Since the Higgs potential has to be significantly altered w.r.t. the one of the SM in order to yield a first-order EW phase transition (FOEWPT), this generally implies the existence of new physics at the EW scale that can be searched for at the LHC.<fnref n="1" /> Moreover, a cosmological first-order phase transition would lead to the generation of a stochastic primordial gravitational wave (GW) background that could be detectable with future space-based gravitational wave observatories, such as the Laser Interferometer Space Antenna (LISA) [17, 18].

The possibility of accommodating a strong FOEWPT in the 2HDM has been studied abundantly in the past (see e.g. [11, 12, 14, 19–24]). It was found that sizable quartic couplings between the scalar states are required to generate a radiatively and thermally induced potential barrier between the symmetry-conserving and the symmetry-breaking vacua of the 2HDM Higgs potential, thus facilitating the presence of a FOEWPT. In addition, assuming that the lightest Higgs boson $h$ is the one that is identified with the detected Higgs boson at 125 GeV, the strength of the transition is maximized (for other parameters fixed) [22] in the so-called *alignment limit* of the 2HDM, where the properties of the state $h$ resemble the ones of the Higgs boson as predicted by the SM (see section 2.1 for a more detailed discussion), and which is accordingly also favored by the measurements of the signal rates of the Higgs boson at 125 GeV. The sizable quartic couplings required to facilitate the FOEWPT imply, in combination with other experimental and theoretical restrictions, that the strongest phase transitions occur in scenarios with relatively large mass splittings between the BSM Higgs bosons [20, 23, 24], in particular characterized by a mass spectrum $m_H \ll m_A \approx m_{H^{\pm}}$. As a consequence of this mass splitting between $H$ and $A$ the decay $A \to ZH$ is kinematically open [20], which (in contrast to the decay $A \to Zh$) remains unsuppressed in the alignment limit of the 2HDM. Due to this feature, together with the result that a larger cross section for the process $pp \to A \to ZH$ is correlated with a stronger phase transition [24, 25], this process has been coined a *smoking gun signature* for a FOEWPT at the LHC [20, 24].

Independently of the nature of the EW phase transition, the fact that two BSM states are involved in the process $A \to ZH$ implies that searches for this decay are of particular phenomenological importance in the 2HDM. It should be noted that measurements of the EW precision observables constrain the amount of weak isospin breaking induced by the second Higgs doublet, enforcing approximate mass degeneracy between either $H$ and $H^{\pm}$ or $A$ and $H^{\pm}$ [26].<fnref n="2" /> Accordingly, a possible detection of a $A \to ZH$ signal would provide information

<footnote n="1">¹See refs. [15, 16] for counter examples.</footnote>

<footnote n="2">²Regarding the experimental value of the W-boson mass that is used in this context, we refer to the combined value excluding the CDF measurement, see ref. [27] for a detailed discussion. Conversely, using the</footnote>
about the entire mass spectrum of the 2HDM. Furthermore, the cross section corresponding to the hypothetical signal would rather precisely fix the otherwise free parameter $\tan \beta$ (defined as the ratio of the vacuum expectation values (vev) of the Higgs fields in the considered type of the 2HDM), which controls the couplings of the BSM scalars to fermions and gauge bosons. Assuming CP conservation in the Higgs sector, the only remaining free parameter of the 2HDM is then the soft mass parameter $m_{12}^2$, which would however be restricted to a specific interval based on theoretical constraints from vacuum stability and perturbativity [30]. Further constraints on the 2HDM parameter space can be obtained from flavor-physics observables. In particular, the measurements of transition rates of radiative $B$-meson decays give rise to a lower limit on the mass of the charged scalars in the type II and the type IV. Using the current experimental world average value $\text{BR}(B \to X_s \gamma) = 3.49 \pm 0.19$ from the HFLAV collaboration [31], one finds a limit of $m_{H^\pm} \gtrsim 500$ GeV [32, 33].$^3$ The application of the flavor-physics constraints relies on the assumption that there are no other BSM contributions to the flavor observables in addition to those of the new 2HDM scalars. As such, these constraints should be regarded as complementary to the constraints from direct searches at colliders, where the latter can be considered as more direct exclusions that are largely independent of other possibly existing BSM effects. Therefore, in this paper, we focus on the constraints from collider searches and do not carry out a detailed evaluation of the exclusion limits from flavor-physics observables and their complementarity with the limits from direct searches.

Searches for the decay $A \to ZH$ have been performed by both the ATLAS and the CMS collaboration utilizing the 8 TeV and 13 TeV data sets. These searches made use of the leptonic decay modes of the $Z$ boson, while for the lighter BSM resonance $H$ decays either into bottom-quark, tau-lepton or $W$-boson pairs were considered [35–38]. Notably, until recently no searches for the $A \to ZH$ decay existed assuming the decay of $H$ into top-quark pairs, which are favored for low $\tan \beta$. Based on the combination of LHC searches and other experimental constraints, however, in the type II and type IV 2HDM (see section 2.1) the allowed region for the scalar spectrum is such that typically $H$ is heavier than twice the top-quark mass, which implies that the decay $H \to t\bar{t}$ is dominant (except for parameter regions with a large enhancement of the couplings to bottom quarks). It was therefore pointed out [24, 25, 39] that the searches for $A \to ZH \to \ell^+\ell^-t\bar{t}$ hold great promise to probe so far unconstrained parameter space regions of the 2HDM, in particular, as discussed above, the regions suitable for a realization of a FOEWPT.$^4$

Recently, the experimental situation has changed with the first public results of searches for the signature $A \to ZH \to \ell^+\ell^-t\bar{t}$ by the ATLAS collaboration utilizing the full Run 2 dataset collected at 13 TeV [41]. In addition, ATLAS also presented new searches using the decay of the $Z$ boson into pairs of neutrinos and assuming the decay of the lighter scalar resonance into bottom-quark pairs, giving rise to the $\nu\bar{\nu}b\bar{b}$ final state. CMS has not yet

<footnote>combined value including the CDF measurement quoted in ref. [27] yields a preference for a non-vanishing mass splitting between $H^\pm$ and the neutral BSM states $H, A$, see e.g. refs. [28, 29]. However, we do not consider this possibility here.</footnote>

<footnote>$^3$We note that the 2022 HFLAV average value for $\text{BR}(B \to X_s \gamma)$ used here does not include the most recent Belle-II result [34], but which is in good agreement with the current HFLAV average value.</footnote>

<footnote>$^4$See also refs. [23, 40] for recent discussions of the smoking gun signature and its possible impact on the 2HDM based on expected sensitivities at the LHC.</footnote>
released a result in the $ℓ^+ℓ^-t\bar{t}$ final state but the experimental analysis is ongoing [42, 43]. In this work we will demonstrate that already the first released LHC experimental result on searches in the $A \to ZH \to ℓ^+ℓ^-t\bar{t}$ channel has a large impact on the parameter space of extended Higgs sectors and on possible scenarios giving rise to a strong FOEWPT, constraining sizable regions of so far allowed 2HDM parameter space. Besides, we discuss the possible phenomenological and cosmological implications of an excess of events with $2.85\sigma$ (local) significance at masses of 450 GeV and 650 GeV for the lighter and heavier BSM resonances, respectively, reported by ATLAS in this search and compatible with a strong FOEWPT in the 2HDM. Furthermore we analyze the strength of the GW signals in the parameter space that is compatible with the observed excess, and discuss the future detectability by LISA. In addition, we identify another potentially promising LHC search to target the region indicative for strong FOEWPT, namely the production of the charged Higgs $H^\pm$, then decaying as $H^\pm \to W^\pm H \to ℓ^\pm νt\bar{t}$, for which so far no results exist.

Performing an extrapolation of the 95% confidence-level expected cross section limits based on the LHC Run 2 presented by ATLAS in [41], we also investigate the future discovery reach of the smoking gun search. If this search were to lead to the detection of a signal, it will be rather straightforward to assess whether this signal can be interpreted within the context of the 2HDM. If so one can infer the allowed ranges of the 2HDM parameters according to the discussion above and make predictions about the associated phenomenology regarding other LHC searches and the nature of the EW phase transition. On the other hand, if a signal is observed that cannot be accommodated by the 2HDM (e.g. because the required mass splitting is too large), this would be a clear indication for physics beyond the 2HDM. Economical extensions would be, for instance, the complex (CP-violating) 2HDM or singlet-extended 2HDMs, where the cross sections for the smoking gun signature and the mass splittings between the scalars are less restricted as a consequence of additional parameters [25, 44]. In this manner, the detection of a signal in the $A \to ZH$ channel would provide invaluable information about the discrimination between different scenarios featuring extended scalar sectors and, therefore, about the underlying physics governing the dynamics of EW symmetry breaking.

The outline of this paper is as follows: in section 2 we introduce the 2HDM and specify our notation. Moreover, we briefly summarize the basis for the description of the EW phase transition in the early universe with regards to strong FOEWPTs, EW baryogenesis and GWs. In section 3 we present the numerical discussion of the impact of the new ATLAS $A \to ZH$ searches on the 2HDM parameter space. We divide our discussion of the new constraints into two parts, focusing on a low-tan $\beta$ and a high-tan $\beta$ regime in section 3.2.1 and section 3.2.2, respectively. Then, section 3.3 is devoted to an evaluation of the future prospects of the $ℓ^+ℓ^-t\bar{t}$ searches. The excess observed by ATLAS in the $ℓ^+ℓ^-t\bar{t}$ final state as a hint of a FOEWPT in the 2HDM is analysed in section 3.4, and its connection to the possible existence of a stochastic GW signal from the EW phase transition is discussed. We summarize our findings and conclude in section 4. Appendix A is devoted to presenting the python package `thdmTools`, a state-of-art tool for the exploration of the 2HDM which we have used for our analyses. Appendix B contains a comparison of the future projections obtained here based on the new ATLAS results to earlier projections based on expected cross-section limits from CMS.
## 2 Theoretical background

In order to specify our notation, we provide an overview of the 2HDM Higgs sector below. We will also briefly discuss the methods used to investigate the occurrence of a strong FOEWPT, and the related phenomenological consequences, such as the realization of EW baryogenesis or the generation of a primordial GW background.

### 2.1 The Two Higgs Doublet Model (2HDM)

As theoretical framework for our study we consider an extension of the SM Higgs sector by one additional complex SU(2) doublet. Models with two Higgs doublets have been widely studied in the literature [8, 45] and give rise to a rich phenomenology including the possibility of a strong FOEWPT. In the present analysis we assume a CP-conserving 2HDM with a softly broken $Z_2$ symmetry. This discrete symmetry between the doublets leaves the potential invariant (up to the soft breaking) under transformations of the type $\Phi_1 \to \Phi_1$, $\Phi_2 \to -\Phi_2$ in order to avoid the occurrence of large flavour changing neutral currents that would violate the existing bounds. We allow for a soft $Z_2$ -breaking via a mass parameter $m_{12}^2$. Given these constraints, the most general form of the Higgs potential is

$$V(\Phi_1, \Phi_2) = m_{11}^2(\Phi_1^\dagger\Phi_1) + m_{22}^2(\Phi_2^\dagger\Phi_2) - m_{12}^2(\Phi_1^\dagger\Phi_2 + \Phi_2^\dagger\Phi_1) + \frac{\lambda_1}{2}(\Phi_1^\dagger\Phi_1)^2 + \frac{\lambda_2}{2}(\Phi_2^\dagger\Phi_2)^2 \\ + \lambda_3(\Phi_1^\dagger\Phi_1)(\Phi_2^\dagger\Phi_2) + \lambda_4(\Phi_1^\dagger\Phi_2)(\Phi_2^\dagger\Phi_1) + \frac{\lambda_5}{2}[(\Phi_1^\dagger\Phi_2)^2 + (\Phi_2^\dagger\Phi_1)^2], \quad (2.1)$$

where the two doublets are conveniently parameterized as

$$\Phi_1 = \begin{pmatrix} \phi_1^+ \\ \frac{1}{\sqrt{2}}(v_1 + \rho_1 + i\eta_1) \end{pmatrix}, \quad \Phi_2 = \begin{pmatrix} \phi_2^+ \\ \frac{1}{\sqrt{2}}(v_2 + \rho_2 + i\eta_2) \end{pmatrix}. \quad (2.2)$$

Here $v_i$ is the vacuum expectation value (vev) acquired by the respective doublet. The eight degrees of freedom $\phi_i^+$, $\rho_i$ and $\eta_i$ mix to give rise to the massive scalars $h$, $H$ (CP-even), $A$ (CP-odd), $H^\pm$ and the three Goldstone bosons $G^0, G^\pm$. In this work we identify $h$ with the Higgs boson at 125 GeV detected at the LHC. The rotation from the gauge basis to the mass basis involves two mixing matrices with the mixing angles $\alpha$ and $\beta$ for the CP-even and the CP-odd/charged sector, respectively.

The coefficients $m_{11}^2, m_{22}^2, m_{12}^2, \lambda_i$ ($i = 1, \dots, 5$) in eq. (2.1) can be mapped to the physical basis of the masses and mixing angles. We employ the usual convention of the parametrization of the mixing angles as $\tan \beta = v_2/v_1$, and $\cos(\beta - \alpha)$. In the limit $\cos(\beta - \alpha) \to 0$ — the alignment limit — [46] the light Higgs boson $h$ in the 2HDM has the same tree-level couplings to the SM fermions and gauge bosons as the SM Higgs boson. In the numerical discussion, we use the following basis of free parameters of the model,

$$\cos(\beta - \alpha), \quad \tan \beta, \quad v, \quad m_h, \quad m_H, \quad m_A, \quad m_{H^\pm}, \quad m_{12}^2, \quad (2.3)$$

where two of these parameters are already fixed by experiment, namely $m_h \approx 125$ GeV and $v = \sqrt{v_1^2 + v_2^2} \approx 246$ GeV.

Extending the $Z_2$ symmetry to the Yukawa sector leads to four different 2HDM types depending on the $Z_2$ parities of the fermions, summarized in table 1. It should be noted
<table><caption>**Table 1.** Couplings of the two Higgs doublets to the SM fermions in the four types of the 2HDM. Right: tan β-dependence of the couplings of the up-type and down-type quarks to the CP-odd Higgs boson A.</caption><thead><tr><th></th><th>u-type</th><th>d-type</th><th>leptons</th><th>ξ<sup>u</sup><sub>A</sub></th><th>ξ<sup>d</sup><sub>A</sub></th></tr></thead><tbody><tr><td>type I</td><td>Φ<sub>2</sub></td><td>Φ<sub>2</sub></td><td>Φ<sub>2</sub></td><td>cot β</td><td>-cot β</td></tr><tr><td>type II</td><td>Φ<sub>2</sub></td><td>Φ<sub>1</sub></td><td>Φ<sub>1</sub></td><td>cot β</td><td>tan β</td></tr><tr><td>type III (lepton-specific)</td><td>Φ<sub>2</sub></td><td>Φ<sub>2</sub></td><td>Φ<sub>1</sub></td><td>cot β</td><td>-cot β</td></tr><tr><td>type IV (flipped)</td><td>Φ<sub>2</sub></td><td>Φ<sub>1</sub></td><td>Φ<sub>2</sub></td><td>cot β</td><td>tan β</td></tr></tbody></table>

that due to the opposite $Z_2$ parities of the Higgs doublets each fermion can only be coupled to either $\Phi_1$ or $\Phi_2$. In table 1 we also give the tan β dependence of the coupling modifiers of the CP-odd Higgs boson A of the 2HDM to up-type and down-type quarks, $\xi_A^{u,d}$ (see ref. [47] for a formal definition of these parameters). To a very good approximation the cross sections for the production of A at the LHC are proportional to the square of the coupling modifiers $\xi_A^{u,d}$. Since $\xi_A^u = \cot \beta$ for all types, the dominant production mode for values of $\tan \beta \approx 1$ is gluon-fusion production involving a top-quark loop. For values of $\tan \beta \gtrsim 10$ the production in association with bottom-quark pairs becomes important in type II and type IV, in which $\xi_A^d = \tan \beta$, whereas in type I and type III this production mode is always substantially smaller than gluon-fusion production.

## 2.2 Thermal history analysis

In order to study the physics of the EW phase transition, we will use the finite-temperature effective potential formalism. The one-loop, daisy resummed, finite-temperature 2HDM effective potential is given as

$$V_{\text{eff}} = V_{\text{tree}} + V_{\text{CW}} + V_{\text{CT}} + V_{\text{T}} + V_{\text{daisy}}. \qquad (2.4)$$

The temperature-independent part of the potential comprises the first three terms, where $V_{\text{tree}}$ is defined in eq. (2.1), $V_{\text{CW}}$ is the one-loop Coleman Weinberg potential [48] incorporating the radiative corrections, and $V_{\text{CT}}$ is a UV-finite counterterm potential introduced in order to keep the physical masses and the vevs of the Higgs fields at their tree-level values at zero temperature [21]. The thermal corrections to the scalar potential are split into two terms. The first one, $V_{\text{T}}$, incorporates the one-loop thermal corrections in terms of the well-known J-functions (see e.g. ref. [49]). The second term, $V_{\text{daisy}}$, is an additional piece accounting for the resummation of the so-called daisy-diagrams which signal the breakdown of fixed-order perturbation theory at finite temperature. As resummation prescription, we follow the Arnold-Espinosa method [50], which resums the infrared-divergent contributions from the bosonic Matsubara zero-modes. We emphasize that the computation of the finite-temperature effective potential, at the order performed in this work, is affected by sizable theoretical uncertainties, see [51–54] for a detailed discussion. As a consequence, the regions studied in the following should only be regarded as indicative for the presence of a strong FOEWPT.
## 2.3 Vacuum tunneling

In order to find regions in the parameter space of the 2HDM that feature a FOEWPT we track the co-existing vacua as a function of the temperature using the effective potential from eq. (2.4), by means of a modified version of the public code CosmoTransitions [55]. Typically, the universe evolves starting from an EW symmetric vacuum configuration at the origin of field space.<fnref n="5" /> We identify the 2HDM parameter space regions which, as the universe cools down, feature an EW-breaking minimum of the Higgs potential that is separated from the minimum in the origin by a potential barrier. The universe reaches the critical temperature $T_c$ when these two coexisting vacua are degenerate. At later times, when $T < T_c$, the minimum corresponding to the EW vacuum drops below the minimum in the origin, and thus becomes energetically more favorable. At this point, the onset of the first-order phase transition from the minimum at the origin to the EW vacuum depends on the transition rate per unit time and unit volume. The transition rate in turn depends on the temperature-dependent Euclidean bounce action $S(T)$ of the (multi-)scalar field configurations. The onset of the phase transition occurs when (see e.g. ref. [18])

$$S(T_n)/T_n \approx 140, \qquad (2.5)$$

which arises from the comparison of the transition rate and the expansion rate of the universe. $T_n$ is the nucleation temperature, which very accurately corresponds to the temperature at which the transition takes place. If the condition (2.5) is not fulfilled at any temperature $T < T_c$, the phase transition cannot complete successfully, and the universe remains trapped in the false vacuum at the origin [24] (see also refs. [25, 56]).

## 2.4 Strong first-order electroweak phase transitions and baryogenesis

The origin of the baryon asymmetry of the universe (BAU) is one of the main open questions of modern particle physics. In extensions of the SM such as the 2HDM it is possible to dynamically generate an excess of matter over antimatter in the universe by means of EW baryogenesis. According to the Sakharov conditions [57], a first-order EW phase transition is a necessary ingredient for the realization of EW baryogenesis as it provides the required conditions to bring the thermal plasma out of equilibrium. In order to avoid the washout of the asymmetry after the phase transition, the EW vev after the transition should be larger than the transition temperature, i.e.

$$\xi_n = \frac{v_n}{T_n} \gtrsim 1, \qquad (2.6)$$

where $v_n$ is the EW vev at the nucleation temperature $T_n$. The above condition (2.6), which defines a strong FOEWPT, yields the parameter region satisfying baryon-number preservation [58] after the transition, which is therefore suitable for EW baryogenesis. It should be noted that a successful realization of EW baryogenesis also requires BSM sources of CP violation. In the present paper we restrict ourselves to the condition for a strong FOEWPT

<footnote n="5">⁵EW symmetry non-restoration in the high-temperature regime $T \gg m_A, m_{H^\pm}, m_H, M$ (with $M^2 \equiv m_{12}^2/(s_\beta c_\beta)$) is possible in the 2HDM [24]. For the sake of simplicity, we omit this possibility in the discussion of FOEWPT here.</footnote>
and do not carry out a detailed investigation of the actual realization of EW baryogenesis (see ref. [59] for a recent study). Accordingly, we make the assumption that the additional sources of CP violation that are needed for EW baryogenesis do not have a significant impact on the properties of the FOEWPT (see e.g. ref. [14] for an example in the 2HDM).

## 2.5 Gravitational waves

It is well-known that a cosmological first-order phase transition gives rise to a stochastic gravitational wave signal [60, 61]. Since the EW phase transition would have happened at temperatures comparable to the EW scale, the GW signal spectrum would be largest around milli-Hz frequencies, thus in the best-sensitivity range of the planned LISA space-based GW interferometer [18, 62]. The GWs in a FOEWPT are sourced by the collision of bubbles and the surrounding plasma motions in the form of sound waves [63–66], as well as the turbulence generated after the collisions [67–71] (see ref. [17] for a review). In the case of the 2HDM, the GW contribution from bubble collisions themselves can be neglected, and the GW power spectrum may be modeled with the sound waves as dominant source [14]. There are four phase transition parameters that characterize the corresponding GW signal [17]: (i) the temperature $T_*$ at which the phase transition occurs, which we identify here with the nucleation temperature $T_n$.<fnref n="6" /> (ii) the phase transition strength $\alpha$, defined as the difference of the trace of the energy-momentum tensor between the two vacua involved in the transition, normalized to the radiation background energy density. (iii) the inverse duration of the phase transition in Hubble units, $\beta/H$. (iv) the bubble wall velocity in the rest frame of the fluid (and far from the bubble), $v_w$. To compute $\alpha$ we follow refs. [17, 18],

$$ \alpha = \frac{1}{\rho_R} \left( \Delta V(T_*) - \left( \frac{T}{4} \frac{\partial \Delta V(T)}{\partial T} \right) \bigg|_{T_*} \right), \qquad (2.7) $$

where $\Delta V(T_*)$ is the potential difference between the two vacua evaluated at the temperature of the phase transition, and $\rho_R$ is the radiation energy density of the universe. The inverse duration of the phase transition $\beta/H$ can be generally calculated as

$$ \frac{\beta}{H} = T_* \left( \frac{d S(T)}{dT} \frac{S(T)}{T} \right) \bigg|_{T_*}, \qquad (2.8) $$

where $S(T)$ is (as in eq. (2.5)) the temperature-dependent (3-dimensional) Euclidean bounce action. Finally, based on recent results indicating that phase transition bubbles preferentially expand with either $v_w \approx c_s$ ($c_s$ being the speed of sound of the plasma)<fnref n="7" /> or $v_w \to 1$ [73, 74] (see also ref. [75] for a further discussion of bubble wall velocity estimates in BSM theories) we choose to fix $v_w = 0.6$ as a representative case.

Based on the four quantities introduced above, the primordial stochastic GW background produced during a cosmological phase transition can be computed using numerical power-law fits to results of GW production obtained in hydrodynamical simulations of the thermal plasma. In our numerical analysis, we include the contributions to the GW power spectrum

<footnote n="6">⁶We could instead consider $T_*$ to be the percolation temperature [72], at which the phase transition completes from the percolation of bubbles, yet the numerical difference compared to $T_n$ is very small.</footnote>

<footnote n="7">⁷For a relativistic perfect fluid, $c_s = 1/\sqrt{3} \approx 0.577$.</footnote>
from sound waves $h^2\Omega_{sw}$ and turbulence $h^2\Omega_{turb}$, where sound waves are the dominant GW source for the FOEWPTs considered here. The specific formulas used in our analysis for the computation of the GW spectral shapes, their amplitudes and the peak frequencies can be found in ref. [24], which closely follows refs. [17, 71]. Whether a stochastic GW signal is detectable at a GW observatory depends on the signal-to-noise ratio (SNR), which can be computed for a specific parameter point and a specific GW experiment as

$$ \text{SNR} = \sqrt{\mathcal{T} \int_{-\infty}^{+\infty} df \left[ \frac{h^2 \Omega_{\text{GW}}(f)}{h^2 \Omega_{\text{Sens}}(f)} \right]^2}, \qquad (2.9) $$

where $\mathcal{T}$ is the duration of the experiment, $h^2\Omega_{\text{Sens}}$ is the nominal sensitivity of the detector, computed according to the mission requirements [76], and $h^2\Omega_{\text{GW}} = h^2\Omega_{\text{sw}} + h^2\Omega_{\text{turb}}$ is the spectral shape of the GW signal. For the present analysis, we focus on the GW detectability with LISA, for which we will assume an operation time of $\mathcal{T} = 7$ years, and consider a GW signal to be detectable if SNR > 1 (more stringent SNR detection criteria could also be considered [17]).

## 3 Numerical analysis

In this section we analyze in detail the impact of the recent ATLAS searches for the $A \to ZH$ signature [38, 41] on the 2HDM parameter space. Before we discuss the results of our analysis, we briefly describe the implementation of the new ATLAS limits into the HiggsTools package [77] (which contains HiggsBounds [78–81]), discuss the public numerical tools that were used in our analysis and introduce the software package thdmTools that we have developed as part of this work.

### 3.1 Implementation of new limits into HiggsBounds

In order to confront the 2HDM predictions for the different regions of the parameter space with the new cross section limits reported by the ATLAS Collaboration, we have implemented the 95% confidence level (C.L.) expected and observed cross section limits into the HiggsBounds dataset, corresponding to the following new ATLAS results that were not yet contained in the public HiggsBounds dataset [77]:

$$ - gg \to A \to ZH \to \ell^+ \ell^- b \bar{b} \text{ at } 13 \text{ TeV including } 139 \text{ fb}^{-1} \text{ from ref. [38]} $$

$$ - b\bar{b} \to A \to ZH \to \ell^+ \ell^- b\bar{b} \text{ at } 13 \text{ TeV including } 139 \text{ fb}^{-1} \text{ from ref. [38]} $$

$$ - gg \to A \to ZH \to \ell^+ \ell^- t\bar{t} \text{ at } 13 \text{ TeV including } 140 \text{ fb}^{-1} \text{ from ref. [41]} $$

$$ - gg \to A \to ZH \to \nu\nu b\bar{b} \text{ at } 13 \text{ TeV including } 140 \text{ fb}^{-1} \text{ from ref. [41]} $$

$$ - b\bar{b} \to A \to ZH \to \nu\nu b\bar{b} \text{ at } 13 \text{ TeV including } 140 \text{ fb}^{-1} \text{ from ref. [41]} $$

We note that the results from ref. [38] update the previous ATLAS results based on $36.1 \text{ fb}^{-1}$ collected during the first two years of Run 2 [36], which are contained in the public HiggsBounds dataset (and which will now be replaced by the full Run 2 results). The corresponding CMS results include first-year Run 2 data [35] and are also implemented in
HiggsBounds, but since they are based on less data the extracted limits are weaker than the limits from the new ATLAS analyses.

In our analysis below we use HiggsBounds to determine the parameter regions in the considered 2HDM scenarios that are excluded at the 95% C.L. by the existing limits from Higgs searches. In order to ensure the correct statistical interpretation of the excluded regions as a limit at the 95% C.L., HiggsBounds applies for each BSM scalar only the specific observed limit that has the highest sensitivity according to the expected limits of all considered searches. For illustration, we furthermore display the regions that would be excluded by different searches if each of these searches were applied individually. In this way the impact of the new $A \to ZH$ searches from ATLAS in the different final states becomes clearly visible, and one can assess to what extent these searches probe parameter regions that were so far unexplored. We note, however, that the requirement for a certain parameter point to be simultaneously in agreement with the 95% C.L. exclusion limits of all available searches would in general result in an excluded region that would be too aggressive in comparison to the parameter region corresponding to an overall exclusion at the 95% C.L. (which we determine using HiggsBounds as described above).

For the application of the HiggsBounds cross section limits, one has to provide the relevant model predictions for the cross sections and branching ratios of the scalar states. To this end, we utilized cross-section predictions from the corresponding part of the HiggsTools package [77] based on the effective-coupling input (see refs. [77, 82] and references therein for details), and branching ratios of the Higgs bosons as obtained using AnyHdecay [83-85].

In order to automate the interface to these numerical tools in 2HDM analyses, we have developed the public thdmTools package. Beyond the interfaces for HiggsTools (containing HiggsBounds and HiggsSignals [77, 86, 87]), where the latter tests the properties of the detected Higgs boson at about 125 GeV against the LHC Higgs boson rate measurements) and AnyHdecay, thdmTools also includes additional interfaces for assessing 2HDM parameter points against various experimental constraints, including those from electroweak precision observables and flavor physics (where the latter are not explicitly applied in the following analyses, see the discussion below). Moreover, thdmTools can be used to check against theoretical constraints from boundedness from below and perturbative unitarity. A brief description of the functionalities of thdmTools as well as instructions for download and installation can be found in appendix A.2.

## 3.2 New constraints on the 2HDM parameter space

In our analysis we target two separate 2HDM parameter regions: a low-tan $\beta$ regime, corresponding to values of $\tan\beta \le 3$, and a high-tan $\beta$ regime with values of $\tan\beta \ge 10$. In the low tan $\beta$-regime, $A$ can be produced with sizable cross sections via gluon-fusion, and we will show that searches in the $l^+l^-t\bar{t}$ final state exclude significant parts of previously unconstrained parameter space, whereas the ATLAS searches utilizing $\nu\bar{v}\bar{b}b$ final states only provide exclusions on the 2HDM parameter space that were already excluded by other LHC searches or theoretical requirements from perturbative unitarity. In the high-tan $\beta$ regime, the state $A$ can be produced with sizable rates in $b\bar{b}$-associated processes for the 2HDM of
Yukawa type II and IV. In this regime, the $\nu\bar{\nu}b\bar{b}$ searches cover regions of parameter space that have already been probed via the $\ell^+\ell^-b\bar{b}$ final state.

For intermediate values of $\tan\beta$, in between the low- and high-tan $\beta$ regimes considered here, the new ATLAS searches cannot probe substantial parts of the 2HDM parameter space. The reason is that the gluon-fusion production cross sections of $A$, dominantly generated via a top-loop diagram, are roughly proportional to $1/\tan^2\beta$ for intermediate $\tan\beta$ values, such that the searches in the $\ell^+\ell^-t\bar{t}$ final states quickly lose sensitivity as $\tan\beta$ increases. At the same time the $b\bar{b}$-associated production of $A$ is enhanced in type II and type IV by factors of $\tan^2\beta$, such that the searches relying on this production mechanism become more sensitive for larger values of $\tan\beta$. Yet, in our numerical analysis we find that one needs roughly an enhancement corresponding to $\tan^2\beta \approx 100$ to achieve cross sections in the $b\bar{b}$-associated production mode that are comparable to the exclusion limits resulting from the $gg \to A \to ZH$ searches. Consequently, we chose the value $\tan\beta = 10$ as a representative benchmark scenario in order to assess the sensitivity of the new searches in the high-tan $\beta$ regime. The impact of the searches for even larger values of $\tan\beta$ can be extrapolated from the discussion of this scenario, as will be shown in detail in section 3.2.2.

As already discussed in section 1, the 2HDM parameter region that we explore in this work is motivated by the possibility of realizing a strong FOEWPT giving rise to EW baryogenesis in the early universe. In this scenario the CP-odd scalar $A$ and the charged scalars $H^{\pm}$ are assumed to be mass-degenerate, i.e. $m_A = m_{H^{\pm}}$, and the squared mass scale $M^2 = m_{12}^2/(s_\beta c_\beta)$ is set equal to the mass of the heavy CP-even scalar $H$, i.e. $M = m_H$.<fnref n="8" /> In addition, the alignment limit $\cos(\alpha - \beta) = 0$ is assumed, in which the properties of the (in this case) lightest Higgs boson $h$ with mass $m_h = 125$ GeV are the same as for the SM Higgs boson at tree level. These conditions on the parameter space allow for sizable $m_A - m_H$ mass splittings, driven by the quartic couplings in the 2HDM scalar potential (2.1), facilitating the presence of a FOEWPT [20, 24] while being in agreement with the LHC measurements of the properties of the detected Higgs boson at 125 GeV as well as with the results for the EW precision observables and further theoretical constraints. After imposing the above-mentioned conditions, the only remaining free parameters are the masses $m_H$ and $m_A = m_{H^{\pm}}$ as well as $\tan\beta$. In the following, we discuss our results in the $(m_H, m_A)$-plane for different values of $\tan\beta$ within the two regimes discussed above.

### 3.2.1 Low $\tan\beta$-region

In this section we present our results for the low-$\tan\beta$ regime, which focuses on the range $1 \le \tan\beta \le 3$. The lower bound on $\tan\beta$ was chosen because values of $\tan\beta$ below 1 are in strong tension with constraints from flavor physics observables. The indirect limits from flavor physics also constrain the 2HDM parameter space for slightly larger values of $\tan\beta$ depending on the 2HDM Yukawa type and the mass of the charged Higgs boson. As

<footnote n="8">⁸We note that in the 2HDM interpretation presented by ATLAS almost the same benchmark scenario was considered [41]. However, therein the condition $M = m_A$ was used instead of $M = m_H$ as applied here. We apply the latter condition in order to have a theoretically consistent form of the Higgs potential for $m_A = m_{H^{\pm}} \gg m_H$, whereas the condition used by ATLAS gives rise to an unbounded Higgs potential, and thus an unstable EW vacuum, in the parameter space regions in which the $A \to ZH$ decay is kinematically allowed.</footnote>
discussed above, we do not carry out a detailed investigation of the indirect limits from flavor physics in the following.

In figure 1 we show the impact of the new $A \to ZH$ searches from ATLAS [38, 41] in the $(m_H, m_A)$-plane for $\tan\beta = 1$ (upper left), $\tan\beta = 1.5$ (upper right), $\tan\beta = 2$ (lower left) and $\tan\beta = 3$ (lower right). The upper left plot is valid independently of the chosen 2HDM Yukawa type. However, for $\tan\beta \neq 1$ the relevant cross sections and branching ratios depend on the Yukawa type, and the specific choice of type that is specified in the upper right and the lower plots will be further discussed below. In each plot we indicate the parts of the parameter space that are excluded by vacuum stability and perturbative unitarity with pink and cyan colors, respectively. The regions excluded by the new ATLAS search [41] for $gg \to A \to ZH$ in the $\ell^+\ell^-t\bar{t}$ and the $\nu\bar{v}\bar{b}b$ final states are indicated with red and blue shaded contours, respectively, whereas regions excluded from previous LHC searches (including the recent ATLAS $gg \to A \to ZH \to \ell^+\ell^-b\bar{b}$ search [38]) are indicated in gray. In each case the search channel giving rise to the exclusion (under the assumption that this search is applied individually, see the discussion in section 3.1) is stated in the plots. For the new ATLAS searches we show in addition the expected exclusion regions with dashed lines in the same colors. By comparing the gray shaded areas with the red and blue shaded areas, one can determine to what extent the new ATLAS searches probe previously unexplored parameter space regions.

While, as discussed in section 3.1, the excluded regions resulting from applying each of the searches individually are shown for illustration, we obtain the region that is overall excluded at the 95% C.L. from the application of HiggsBounds. This region, which is obtained by applying for each BSM scalar only the observed limit of the search that has the highest expected sensitivity, is indicated with the black dotted-dashed lines. If observed limits show a significant excess or a significant underfluctuation in comparison to the expected limits, the overall limit at the 95% C.L. obtained in this way can be weaker than the exclusion that would result from the requirement that each limit should be fulfilled individually. This feature is visible in the upper left plot for $m_A$ values slightly above 700 GeV. Here the $gg \to A \to ZH \to \ell^+\ell^-t\bar{t}$ channel has the highest sensitivity, but since the observed limit is weaker than the expected one, a parameter region that could nominally be excluded by the $gg \to H \to t\bar{t}$ search (indicated in gray) remains unexcluded because of the adopted procedure for obtaining a 95% C.L. exclusion.

Finally, we show in the plots in figure 1 the parameter regions that exhibit a strong FOEWPT as defined in section 2.2 (based on the one-loop thermal effective potential with daisy resummation in the so-called Arnold-Espinosa scheme). As discussed above, we note the sizable theoretical uncertainties in the predictions for a strong FOEWPT using this approach, and thus the regions shown should only be regarded as indicative for the presence of such strong transitions. The color coding of the points indicates the ratio between the vev $v_n$ in the broken phase at the onset of the transition and the nucleation temperature $T_n$.

### $\tan\beta = 1$, all types

The results for $\tan\beta = 1$ are shown in the upper-left plot of figure 1. One can see that the new $A \to ZH \to \ell^+\ell^-t\bar{t}$ ATLAS search (red) excludes the region 350 GeV $\lesssim$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Impact of the new ATLAS searches for the $A \to ZH$ signature in the $(m_H, m_A)$-plane for $\tan\beta = 1$ (upper left), $\tan\beta = 1.5$ and type II (upper right), $\tan\beta = 2$ and type IV (lower left), and $\tan\beta = 3$ and type IV (lower right). Parameter space regions excluded by vacuum stability or perturbative unitarity are indicated with pink and cyan colors, respectively. Regions excluded from previous LHC searches are indicated in gray, and regions excluded by the new $\ell^+\ell^-t\bar{t}$ and $\nu\bar{\nu}b\bar{b}$ searches are indicated in red and blue, respectively, where the dashed lines indicate the corresponding expected exclusion limits. Parameter space regions featuring a FOEWPT with $v_n/T_n > 1$ are indicated with the scatter points, where the color coding indicates the values of $v_n/T_n$. The mass values of the most significant excess (2.85 $\sigma$ local significance) observed by ATLAS in the $\ell^+\ell^-t\bar{t}$ search are indicated with a magenta star in the upper right plot.</figcaption></figure>
$m_H \lesssim 450$ GeV and $650$ GeV $\lesssim m_A = m_{H^\pm} \lesssim 800$ GeV, which was so far allowed. This demonstrates the exclusion power of such smoking-gun signature for masses above the di-top threshold. In addition, when combined with searches for the charged scalars using the $H^\pm \to tb$ decay [88, 89], searches for neutral scalars decaying into top-quark pairs [90], and searches for the $A \to ZH$ decay in the $\ell^+\ell^- b\bar{b}$ final state [38], the mass range $300$ GeV $\lesssim m_H \lesssim 450$ GeV and $450$ GeV $\lesssim m_A = m_{H^\pm} \lesssim 700$ GeV is now excluded. Figure 1 also highlights that for $\tan\beta = 1$ the parameter region with a strong FOEWPT to which the new ATLAS search is sensitive, assuming $m_A = m_{H^\pm}$, is already excluded by the charged Higgs-boson searches. Yet, we stress that if the condition $m_A = m_{H^\pm}$ were relaxed by allowing for an additional mass gap between these states, i.e. $m_{H^\pm} > m_A$ (which however would lead to a tension with the electroweak precision observables), the searches for the charged scalars would become less sensitive, such that the smoking gun search would have the highest sensitivity in an even larger region of parameter space.

### $\tan\beta = 1.5$, type II

The results for $\tan\beta = 1.5$ are shown in the upper right plot of figure 1. While for low $\tan\beta$ values the gluon-fusion production cross sections of $A$ are dominantly mediated by the top-quark loop, making the cross sections still very much independent of the type, the branching ratios of $A$ and $H$ differ depending on the chosen type. However, for $\tan\beta = 1.5$ the differences between the types are mild, and we focus on the Yukawa type II for definiteness. Comparing to the results for $\tan\beta = 1$ (upper left plot), one can see that the region excluded by the searches for the charged scalars via $pp \to H^\pm tb \to tb \bar{b} t$, where the cross section times branching ratio roughly scales with $1/\tan^2\beta$ in the low-$\tan\beta$ regime, is substantially smaller. This search loses even more sensitivity where the decay $H^\pm \to W^\pm H$ is kinematically allowed, giving rise to the slope of the corresponding excluded region for $m_H \lesssim 500$ GeV (which is more pronounced than for $\tan\beta = 1$ because of the reduced $H^\pm tb$ coupling). As a consequence, for $\tan\beta = 1.5$ the $H^\pm \to tb$ searches are not sensitive anymore to the parameter space region indicative of a strong FOEWPT. Instead, this region is excluded up to masses of $m_H \approx 2m_t$ by searches for $H \to \tau^+\tau^-$ [91, 92] and by searches for the $A \to ZH$ decay using the $\ell^+\ell^- b\bar{b}$ final state [38]. Above the di-top threshold, the decay $H \to t\bar{t}$ very quickly dominates, and the new ATLAS search in the $\ell^+\ell^- t\bar{t}$ final state is the most sensitive one. In contrast to the $\tan\beta = 1$ case, for $\tan\beta = 1.5$ the new search is able to exclude a significant parameter region featuring a strong FOEWPT that was previously allowed. The new search substantially pushes the lower limit on the Higgs boson masses to larger values of about $m_H \gtrsim 400$ GeV and $m_A = m_{H^\pm} \gtrsim 550$ GeV. We also stress that, based on the expected cross section limits, an even larger mass region would be excluded, as indicated with the dashed red line. However, ATLAS observed a local $2.85\sigma$ excess for $m_A \approx 650$ GeV and $m_H \approx 450$ GeV, giving rise to a weaker observed cross section limit. The masses corresponding to the excess, indicated with a magenta star in the upper right plot of figure 1, and the corresponding cross section are such that they fall into the strong FOEWPT region. In section 3.4 we will discuss in
greater detail the tantalizing possibility of such an excess to be the first experimental hint of a strong FOEWPT within the 2HDM. We will give a broad characterization of the FOEWPT predicted by this benchmark scenario, focusing on whether the scenario might be suitable for a realization of EW baryogenesis, and whether the associated GW signal might be detectable with LISA.

As an important outcome of the above discussion, a promising complementary LHC search to target the strong FOEWPT region consists of charged scalar production followed by the decay $H^{\pm} \to W^{\pm}H \to l^{\pm}vt\bar{t}$, which so far has not been performed.<fnref n="9" /> In particular, producing the charged scalar via $pp \to tbH^{\pm}$ would in this case lead to a 4-top-like (or 3-top-like, depending on the signal selection) signature, which has very recently been performed by CMS [94] and ATLAS [95] (but not interpreted in terms of the scenario discussed here), yielding a mild excess over the SM expectation.

Finally, it can be seen that for $tan\beta = 1.5$ the new smoking gun search using the $\nu\bar{\nu}b\bar{b}$ final state starts to probe the considered parameter plane. An exclusion region is visible below the di-top threshold regarding $m_H$ and for a minimum amount of mass splitting of $m_A - m_H \gtrsim 200$ GeV. However, in contrast to the searches using the $l^+l^-t\bar{t}$ final state indicated by the red shaded region, the blue shaded region indicating the new exclusion region resulting from the search using the $\nu\bar{\nu}b\bar{b}$ final state is already excluded by previous LHC searches, namely searches for $H$ decaying into tau-lepton pairs [91, 92] and searches for the smoking gun signature $A \to ZH$ with $Z \to l^+l^-$ and the decay of $H$ into bottom-quark pairs [38]. One should note, however, that the new $A \to ZH$ search in the $\nu\bar{\nu}b\bar{b}$ final state covers larger masses up to $m_H = 600$ GeV and $m_A = 1000$ GeV [41], extending the reach of previous ATLAS searches in $l^+l^-b\bar{b}$ and $l^+l^-W^+W^-$ final states [38] in the region with $m_H > 350$ GeV and $m_A > 800$ GeV. In the 2HDM constraints from perturbative unitarity (cyan area in figure 1) exclude large mass splittings between states from the same SU(2) doublet. As a consequence, the extended mass reach of the new searches in the $\nu\bar{\nu}b\bar{b}$ final state (not visible in the plot) does not give rise to new constraints on the 2HDM for $m_A > 800$ GeV. However, in other models allowing for larger mass splittings between the BSM states, the searches in the $\nu\bar{\nu}b\bar{b}$ final state can potentially provide new constraints.

### $tan\beta = 2$, type IV

We show the results for $tan\beta = 2$ in the lower left plot of figure 1. From here on, we focus our discussion on the Yukawa type IV, in which the new ATLAS searches have the highest potential for probing parameter regions that were unconstrained so far. In particular, compared to type I and III the decay width for $H \to b\bar{b}$ is enhanced in type IV for $tan\beta > 1$, such that the searches in the $\nu\bar{\nu}b\bar{b}$ final state become more important with increasing values of $tan\beta$. Moreover, in type IV the decay width for $H \to \tau^+\tau^-$ is suppressed approximately by $1/tan^2\beta$, whereas it is enhanced by about a factor of $tan^2\beta$ in type II. Hence, while in type II the parameter region below the di-top

<footnote n="9">⁹Searches targeting the $H^{\pm} \to W^{\pm}H$ decay have been performed by CMS assuming the decay $H \to \tau^+\tau^-$ and assuming a fixed mass of $m_H = 200$ GeV [93].</footnote>
threshold, i.e. $m_H < 2m_t$, is entirely excluded by the searches for di-tau resonances, in type IV the $\nu\bar{\nu}b\bar{b}$ search can potentially yield stronger constraints.

One can see in the lower left plot of figure 1 that in this case only three LHC searches give rise to excluded regions in the parameter plane. This is a manifestation of the fact that the so-called wedge-region of the 2HDM, with intermediate values of $2 \lesssim \tan\beta \lesssim 8$, is difficult to probe at the LHC [96]. As an example, we note that the searches for the charged scalars via the signature $pp \to H^\pm tb \to tb\bar{tb}$, suppressed by factors of about $1/\tan^2\beta$ in the low-tan $\beta$ regime, cannot probe the parameter plane in this case. Below the di-top threshold, we find that the $A \to ZH$ searches in the $\ell^+\ell^- b\bar{b}$ (gray) and the $\nu\nu b\bar{b}$ (blue) final states exclude the entire region allowed by the theoretical constraints. As discussed above, for $m_A < 800$ GeV searches for the decay $A \to ZH$ using the decay $Z \to \ell^+\ell^-$ have been performed by ATLAS [38], which are more powerful than the new searches using the $Z \to \nu\nu$ decay (the corresponding CMS search using the $Z \to \ell^+\ell^-$ decay covers masses up to $m_A = 1$ TeV, but is based on first-year Run 2 data only [37]). For $m_A > 800$ GeV ATLAS limits exist only from the new searches using the decay $Z \to \nu\nu$ (the resulting exclusion regions are not visible in our plots since in the 2HDM such large mass splittings are excluded by perturbative unitarity, indicated by the cyan area). Above the di-top threshold, the searches relying on the decay $H \to b\bar{b}$ quickly lose their sensitivity to the 2HDM parameter plane. Accordingly, for masses of $H$ substantially larger than twice the top-quark mass the new smoking gun search for the decay $H \to t\bar{t}$ is in fact the only channel that can probe the parameter plane. As indicated with the red shaded area, the searches in the $\ell^+\ell^- t\bar{t}$ final state are able to exclude masses smaller than $m_H \approx 400$ GeV and $m_A \approx 750$ GeV for the lighter and the heavier BSM resonance, respectively. As it is also visible in the plots for $\tan\beta = 1$, $\tan\beta = 1.5$ and $\tan\beta = 2$ of figure 1, the difference between the expected (red dashed line) and the observed (red solid line) exclusion region resulting from the searches using the $\ell^+\ell^- t\bar{t}$ final state arises from the excess observed in the ATLAS search (except for the upper right part of the red region in the plots for $\tan\beta = 1.5$ and $\tan\beta = 2$, where the observed limit is stronger than the expected one).

### $\tan\beta = 3$, type IV

As a final step of the discussion of the low-tan $\beta$ regime we consider a value of $\tan\beta = 3$. The results of our analysis are shown in the lower right plot of figure 1. Again, we focus on the Yukawa type IV (see the discussion above). One can see that in this case the smoking gun searches in the $\ell^+\ell^- t\bar{t}$ final state cannot probe the parameter space as a consequence of the suppression of the gluon-fusion production cross section of $A$. We will discuss in section 3.3 the prospects for probing the benchmark plane for $\tan\beta = 3$ in future runs of the LHC, in which roughly 20 times more integrated luminosity will be collected by both ATLAS and CMS.<fnref n="10" /> At and below the di-top threshold $m_H \approx 2m_t$ the results are similar to the case of $\tan\beta = 2$, where the smoking gun searches relying on the decay $H \to b\bar{b}$ essentially exclude the whole parameter region. One should note that in type IV (and type II) the partial widths for the decays $A, H \to t\bar{t}$ are

<footnote n="10">¹⁰See ref. [24] for an earlier projection based on expected cross section limits reported by CMS.</footnote>
suppressed approximately by $1/\tan^2\beta$, and the partial width for the decay $H \to b\bar{b}$ is conversely enhanced by (approximately) $\tan^2\beta$. As a result, the gray exclusion region from the searches in the $A \to ZH \to \ell^+\ell^-b\bar{b}$ channel extends to slightly larger masses for $\tan\beta = 3$ compared to $\tan\beta = 2$ (lower left plot).

### 3.2.2 High $\tan\beta$-region

In the discussion above, we investigated the low-$\tan\beta$ regime in which the CP-odd Higgs boson $A$ can be produced with a sizable cross section via gluon-fusion. On the other hand, for large values of $\tan\beta \gtrsim 10$ the gluon-fusion production mode is suppressed. In the Yukawa types II and IV of the 2HDM, $A$ is then produced more efficiently via $b\bar{b}$-associated production, which is enhanced by about $\tan^2\beta$ in these types. Consequently, focusing on the high-$\tan\beta$ regime here, we expect the searches for the signature $A \to ZH$ assuming $b\bar{b}$-associated production to become relevant in type II and type IV. Since in type IV the limits from searches for scalar resonances decaying into tau-lepton pairs are substantially weaker (see the discussion above), we investigate here the impact of the new ATLAS searches on the type IV 2HDM parameter space.

It should be noted that the new ATLAS searches reported in ref. [41] only considered the $b\bar{b}$-associated production utilizing the $\nu\bar{\nu}b\bar{b}$ final state, whereas the smoking-gun search utilizing the $\ell^+\ell^-t\bar{t}$ final state was considered only assuming gluon-fusion production of the heavy BSM resonance. Thus, the only relevant searches for the $A \to ZH$ decay in the following discussion will be the previously reported searches utilizing the $\ell^+\ell^-b\bar{b}$ final state [35, 38] and the new searches utilizing the $\nu\bar{\nu}b\bar{b}$ final state [41].

In figure 2 we show our results for $\tan\beta = 10$ as a representative benchmark scenario for the high-$\tan\beta$ regime. The color coding of the exclusion regions and the scatter points is the same as in figure 1, except for the yellow dashed and solid lines indicating the expected and observed exclusion limits resulting from the recent ATLAS search for $b\bar{b} \to A \to ZH \to \nu\bar{\nu}b\bar{b}$, respectively. One can see that the parameter space region excluded by this search (yellow shaded area) lies within the gray shaded area indicating the exclusion from the searches for $b\bar{b} \to A \to ZH \to \ell^+\ell^-b\bar{b}$ [38], which were published previously. Hence, although the new searches based on the decay of the Z-boson into neutrinos are able to probe the 2HDM parameter space for values of $\tan\beta \gtrsim 10$, these regions are already excluded by the searches making use of the decay of the Z-boson into charged leptons. We stress, however, that the new searches using the $\nu\bar{\nu}b\bar{b}$ final state cover a larger mass interval of up to 1.2 TeV for the heavy BSM resonance (not visible in the plot), whereas the corresponding upper limit in the ATLAS searches using the $\ell^+\ell^-b\bar{b}$ final state is about 800 GeV. Therefore, in other models in which larger mass splittings between the heavier and the lighter BSM resonance are possible compared to the 2HDM (where perturbative unitarity implies an upper limit on such mass splittings, see the cyan region in figure 2), the new searches using the $\nu\bar{\nu}b\bar{b}$ final state could potentially give rise to new constraints.

The two LHC searches relevant in figure 2 differ in the targeted decay mode of the Z boson, whose branching ratios are precisely measured. As a consequence, the relative importance of both searches is independent of the 2HDM parameters, especially of $\tan\beta$. We can therefore extrapolate based on the results for $\tan\beta = 10$ shown in figure 2 that also
<figure><img src="image_3.png" /><figcaption>**Figure 2.** As in figure 1, but for $\tan\beta = 10$ in type IV. Parameter space regions excluded by the new $\nu\bar{\nu}b\bar{b}$ searches in the $b\bar{b}$-associated production channel are indicated in yellow, while the yellow dashed line indicates the expected exclusion limit.</figcaption></figure>

for larger values of $\tan\beta$ the searches making use of the $Z \to \ell^+\ell^-$ decay mode are more promising to probe the considered benchmark plane compared to the searches using the $Z \to \nu\bar{\nu}$ decay mode. It should also be taken into account that for larger values of $\tan\beta$ other LHC searches become relevant in type IV.<fnref n="11" /> In particular, searches for new resonances produced in $b\bar{b}$-associated production with subsequent decay into bottom-quark pairs [98], giving rise to four b-jet final states, start to exclude sizable parts of the benchmark plane for $\tan\beta \gtrsim 15$. Moreover, for such values of $\tan\beta$ searches for new resonances produced in association with a photon and decaying into two jets [99] are able to exclude parameter regions especially in the mass-degenerate regime.

## 3.3 Future prospects for $\ell^+\ell^-t\bar{t}$ searches

In section 3.2.1 we have demonstrated that the new ATLAS smoking-gun searches targeting the $\ell^+\ell^-t\bar{t}$ final state exclude sizable parts of previously allowed parameter space of the 2HDM assuming values of $\tan\beta$ not much larger than one. In particular, we have shown that for BSM scalar masses above the di-top threshold and values of $1.5 \lesssim \tan\beta \lesssim 3$ the smoking-gun searches arguably are the most promising of all LHC searches for probing so far unexplored parameter space regions, with the potential to discover additional Higgs bosons that are consistent with a 2HDM interpretation. Due to their exceptional importance, we briefly discuss here the projected sensitivity of the searches for the $A \to ZH$ decay in the $\ell^+\ell^-t\bar{t}$ final state during future runs of the LHC and the High-Luminosity LHC (HL-LHC). As input

<footnote n="11">¹¹In type II, for $\tan\beta \gtrsim 10$ the whole investigated parameter plane is excluded for masses up to about 1 TeV by searches for scalar resonances decaying into tau-lepton pairs [91, 92, 97].</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. As in figure 1 for $\tan\beta = 1.5$ (left) and $\tan\beta = 3$ (right), shown here for type II, but the red dashed lines indicate projected expected exclusion regions assuming integrated luminosities of 300, 600, 1000, 3000 fb$^{-1}$ from future runs of the LHC.</figcaption></figure>

for our projections we use the expected limits from the ATLAS analysis for an integrated luminosity of 140 fb$^{-1}$. This improves upon the previous projections presented in ref. [24] that were obtained based on an estimate of the expected sensitivities from the CMS collaboration.

In appendix B we provide a comparison of the two projections, showing that they are in good agreement with each other in view of the systematic uncertainties of the analyses.

The projected exclusion limits discussed in the following were obtained by re-scaling the expected cross-section limits reported in ref. [41] with future values for the integrated luminosity that will be collected during future runs of the (HL-)LHC, i.e.

$$ \sigma_{\text{proj.}}^{\text{exp. 95\% CL}}(\mathcal{L}, m_H, m_A) = \sigma_{\text{Run 2}}^{\text{exp. 95\% CL}}(m_H, m_A) \sqrt{\frac{140 \text{ fb}^{-1}}{\mathcal{L}}}. \qquad (3.1) $$

Here, $\sigma_{\text{Run 2}}^{\text{exp. 95\% CL}}$ is the expected cross-section limit at 95% confidence level reported by ATLAS based on 140 fb$^{-1}$ collected during Run 2 as a function of the masses of the probed BSM resonances, and $\sigma_{\text{proj.}}^{\text{exp. 95\% CL}}$ is the future projection of the expected cross-section limits depending additionally on the assumed integrated luminosity $\mathcal{L}$. Accordingly, in the projections we only account for the reduction of statistical uncertainties, whereas no assumption is made on possible improvements of systematic theoretical or experimental uncertainties. Moreover, we do not account for the slight increase of the center-of-mass energy at future runs of the LHC and the HL-LHC, operating at 13.6 TeV and 14 TeV, respectively, compared to the Run 2 dataset collected at 13 TeV. Taking this into account, we consider our projections as fairly conservative estimates.

The projected expected cross section limits can be cast into projected exclusion regions in the 2HDM. In figure 3 we show our projections in the 2HDM benchmark plane introduced in section 3.2 for the Yukawa type II with $\tan\beta = 1.5$ in the left plot and $\tan\beta = 3$ in the right plot.
the right plot. In both plots, the color coding of the scatter points and the definition of the pink and cyan regions is as in figure 1, and the red dashed lines indicate the expected exclusion regions for different values of the integrated luminosity, ranging from $\mathcal{L} = 300 \text{ fb}^{-1}$ (end of LHC Run 3) to $\mathcal{L} = 3000 \text{ fb}^{-1}$ (end of the LHC high-luminosity phase). Moreover, in the left plot the red shaded area indicates the currently excluded region based on the observed cross section limits obtained for $\mathcal{L} = 140 \text{ fb}^{-1}$, and the magenta star indicates the masses for which ATLAS has observed the most pronounced local excess (see section 3.4). As already discussed in section 3.2.1, currently the smoking-gun searches are not able to probe the benchmark plane for $\tan \beta = 3$ (see the lower left plot of figure 1). Accordingly, no red shaded region is visible in the right plot of figure 3.

One can observe in the left plot of figure 3 that with the prospective improvements of the integrated luminosity it will be possible to increase very significantly the regions that can be probed in the considered benchmark plane for $\tan \beta = 1.5$. While currently in the upper right part of the red shaded region the smoking-gun searches are able to exclude masses up to values slightly below 500 GeV for the lighter and up to 850 GeV for the heavier BSM scalar, in the future the LHC will be able to probe via this search masses up to about 700 GeV and 1 TeV for the lighter and the heavier BSM scalar, respectively. This improvement in sensitivity has a very important impact on the parameter region that is suitable for the realization of a strong FOEWPT according to the thermal effective potential approach (as described in section 2.2). In the case of the absence of a signal the exclusion within the region that is indicative for a strong FOEWPT would extend up to $m_H \lesssim 550 \text{ GeV}$ and $m_A \lesssim 700 \text{ GeV}$. It should be noted in this context that the strength of the phase transition diminishes with increasing masses of the BSM scalars. As one can infer from the color coding of the displayed points, the projected exclusion regions cover the parameter region for which the strongest phase transitions can be accommodated. As a result, and since in the 2HDM the generation of a sufficiently large BAU may be possible only for small values of $\tan \beta$ not much larger than one [12], the searches for the smoking-gun signature will provide a stringent test of the possibility to explain the BAU by means of EW baryogenesis in the 2HDM.

In this context it is also important to note that in the 2HDM the primordial GW background generated during the phase transition is only potentially detectable with LISA for the largest possible values of $v_n/T_n$, which are only reached in a very restricted region of the 2HDM parameter space and have a very strong dependence on the details of the scalar spectrum [24]. We have verified using the approach discussed in section 2.5 that for the considered values of $\tan \beta$ all parameter points predicting a GW signal that is potentially detectable with LISA would be probed by the projected exclusion limits from the HL-LHC. Hence, in the 2HDM the HL-LHC results will have an enormous impact on the possibility for a detection of a GW background with LISA consistent with a FOEWPT. This exemplifies that the HL-LHC has the potential to probe large parts of the relevant parameter space before the LISA experiment will have started its operation. Here it should be noted, however, that the presence of a strong FOEWPT, without demanding a realization of EW baryogenesis, is also possible for larger values of $\tan \beta$, where the $gg \to A \to ZH \to \ell^+\ell^- t\bar{t}$ searches lose their sensitivity. A GW signal potentially detectable with LISA therefore cannot be fully probed with the searches in the $\ell^+\ell^- t\bar{t}$ final state.
Besides the analysis of the potential of future runs of the (HL-)LHC for probing the 2HDM parameter space in terms of projected exclusion limits, it is also of interest to investigate the possible interplay between the LHC and LISA for the case where a smoking-gun signal would be detected. We note in this context that the magenta star indicating the mass values corresponding to the excess observed by ATLAS lies well within the discovery reach of the LHC, quite possibly already after the end of Run 3. The detection of the smoking-gun signal would allow for the determination of $m_H$ and $m_A$, and possibly also of $m_{H^\pm}$ via the corresponding cross sections in combination with the application of other constraints. The experimentally determined values of the BSM scalar masses could then be used in dedicated analyses of the phase transition dynamics. For instance, the experimental information about the mass hierarchy of the scalar spectrum would allow an analysis of the thermal potential in an appropriately chosen dimensionally-reduced effective-field theory, in which the heavier scalars have been integrated out in a systematic way in order to facilitate the incorporation of relevant higher-order effects, as well as dedicated lattice simulations (see refs. [100, 101] for recent efforts towards these directions in the 2HDM, and refs. [53, 54, 102–104] for related investigations in other extended scalar sectors).

In the right plot of figure 3, in which we show the projections for $\tan \beta = 3$, one can see that with more integrated luminosity the (HL)-LHC also in this case is able to probe substantial parts of the otherwise unconstrained parameter space regions. Interestingly, the red dashed lines indicating the expected reach of the LHC stretch out to the largest values of $m_H$ within the parameter regions which might be suitable for a realization of a strong FOEWPT. Assuming an integrated luminosity of 3000 fb$^{-1}$ collected by both ATLAS and CMS by the end of the LHC high-luminosity phase, masses of up to $m_H \approx 550$ GeV and $m_A \approx 800$ GeV can be probed. Here it should be taken into account that the parameter space region with $m_H$ below the di-top threshold is already excluded by di-tau searches (only for type II) and by searches for $gg \to A \to ZH \to \ell^+\ell^- b\bar{b}$ (both for type II and type IV), as was discussed in detail in section 3.2.1 (see the lower right plot of figure 1). However, the sensitivity of these searches to the parameter space regions above the di-top threshold will not improve significantly with increasing data, because the branching ratio for the decay $H \to b\bar{b}$ is strongly suppressed for $m_H > 2m_t$.

In summary, the fact that the smoking-gun search for $gg \to A \to ZH \to \ell^+\ell^- t\bar{t}$ will be able to probe masses of $m_H > 2m_t$ in the low $\tan \beta$ regime in the future is crucial for testing the 2HDM parameter space regions suitable for an explanation of the BAU via EW baryogenesis.

## 3.4 A hint of a strong 1st-order EW phase transition in the 2HDM?

We now turn to the analysis of the local $2.85 \sigma$ excess observed by ATLAS. We investigate to what extent a scenario with the mass values of $m_A = 650$ GeV and $m_H = 450$ GeV can be accommodated in the 2HDM and how this scenario can be further tested in the future.<fnref n="12" /> To this end, we first determine the cross section that would be associated with the excess. Since

<footnote n="12">¹²While here we consider the possibility of new physics being the origin of the observed excess, we note that an excess in the $t\bar{t}Z$ final state falls within a class of experimental discrepancies with respect to the SM predictions observed at the LHC in multi-lepton $t\bar{t} + X$ final states, e.g. for $t\bar{t}W$ [105, 106], $t\bar{t}h$ [107, 108] and $t\bar{t}t\bar{t}$ [94, 95] production. The possibility of a mismodeling of the SM expectation in multi-lepton $t\bar{t} + X$ final states should also be investigated in this context.</footnote>
information on the likelihoods has not been made public by ATLAS, we settle here for an approximation based on the reported 95% C.L. cross-section limits. For the mass hypothesis stated above, ATLAS found an expected and observed limit of 0.299 pb and 0.762 pb, respectively. The best-fit signal cross section can be estimated in Gaussian approximation as the difference between observed and expected limit. Furthermore, again assuming a Gaussian distribution of the underlying likelihood, we can determine a symmetric uncertainty of the cross sections in such a way that the background-only hypothesis deviates by the observed local significance from the central value. In this way we obtain a cross section of

$$ \sigma(gg \to A) \times \mathrm{BR}(A \to ZH) \times \mathrm{BR}(H \to t\bar{t}) = 0.46 \pm 0.16 \text{ pb}, \quad (3.2) $$

corresponding to the excess observed by ATLAS at the above-mentioned mass values for the two types of Higgs bosons.<fnref n="13" />

In order to investigate an interpretation of the excess within the 2HDM, we utilize the same benchmark scenario as before, but now fixing $m_H = 450$ GeV and $m_A = m_{H^\pm} = 650$ GeV. While we adopt the mass values for which the excess observed by ATLAS is most pronounced, we note that the mass resolution of the ATLAS search is rather coarse. Thus, the excess would also be compatible with mass values in the vicinity of the specified values, and the overall conclusions regarding a description of the excess in the 2HDM would be unchanged in this case. On the other hand, it is known that in the 2HDM the GW signals produced in FOEWPTs are very sensitive to the details of the spectrum of the scalar masses [24]. In the discussion of the GW signals of this benchmark point and the analysis whether a potentially detectable GW signal would be predicted based on the excess, we will therefore vary the masses of A and H in a mass window of $\pm 50$ GeV which we use as a rough estimate for the potential signal region.

### 3.4.1 Preferred parameter regions

Accommodating the observed excess within the 2HDM implies also constraints on the other 2HDM parameters (besides the BSM scalar masses). It follows from the discussion in section 3.2 that values of $\tan\beta \approx 1.5$ are required for obtaining sufficiently large cross sections to describe the excess (see the upper right plot of figure 1). We therefore vary $\tan\beta$ in the interval $1 \le \tan\beta \le 2$, and in addition we consider deviations from the alignment limit in terms of the free parameter $\cos(\beta - \alpha)$, see the discussion in section 2.1. Since we only consider $\tan\beta$ values close to $\tan\beta = 1$, the theoretical predictions for the cross section depend only marginally on the Yukawa type. However, it should be taken into account that the allowed ranges of $\cos(\beta - \alpha)$, which is constrained (among others) by the cross-section measurements of the Higgs boson at 125 GeV, can be different in the different types. We restrict the discussion here to type I and II as representative examples.

In figure 4 we show the parameter plane with $\cos(\beta - \alpha)$ on the horizontal axis and $\tan\beta$ on the vertical axis and the additional 2HDM parameters set according to the discussion above,

<footnote n="13">¹³Assuming a Gaussian likelihood, another way of determining the uncertainty is by dividing the expected limit by 1.64, which would result in an uncertainty of 0.18 pb, which is in good agreement with the value given in eq. (3.2) considering the relatively low significance of the excess and the resulting size of the cross-section uncertainty.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** For a description of the excess observed in the ATLAS search within the 2HDM the green shaded regions are preferred at the level of $1\sigma$ in type I (left) and type II (right). The olive shaded regions are disfavored by the cross-section measurements of the 125 GeV Higgs boson by more than $2\sigma$ compared to the SM. The grey shaded regions are excluded by LHC cross-section limits from searches for charged Higgs bosons and from searches for resonant Higgs-boson pair-production using the $H \to hh$ decay in $b\bar{b}\tau^+\tau^-$ final states.</figcaption></figure>

for the 2HDM type I in the left and for type II in the right plot. The olive colored regions are disfavored based on the LHC cross-section measurements of the Higgs boson $h$ at 125 GeV, where we used HiggsSignals [86, 87, 109] (incorporated in HiggsTools [77]) to perform a $\chi^2$-fit to the various measurements. Specifically, we demand that the $\chi^2_{125}$ value of a given 2HDM parameter point arising from the measured properties of the detected Higgs boson at about 125 GeV has to be less than $2\sigma$ ($\Delta\chi^2_{125} < 6.18$) away from the SM result ($\chi^2_{125, SM} = 117.7$). Since up to now the LHC measurements regarding the properties of the Higgs boson at 125 GeV are in agreement with the SM predictions, the allowed 2HDM parameter region is located around $\cos(\beta - \alpha) = 0$. The gray regions in the two plots are excluded by the cross section limits from LHC searches for new Higgs bosons. Specifically, we find that the exclusion regions in the plot arise from two LHC searches. The gray regions at $\tan\beta \lesssim 1.2$ (both plots) are excluded by the cross section limits from searches for charged Higgs bosons using the $H^\pm \to tb$ decay [88, 89]. The gray regions at $|\cos(\beta - \alpha)| \gtrsim 0.025$ (both plots) are excluded by searches for resonant Higgs-boson pair production using the $H \to hh$ decay in the $b\bar{b}\tau^+\tau^-$ final state [110, 111]. It should be noted here that searches for resonant pair production of the Higgs boson at 125 GeV cannot probe the alignment limit of $\cos(\beta - \alpha) = 0$, since the triple scalar $hhH$ coupling vanishes at tree level in this limit. The green lines are contour lines indicating the predicted values of the cross section times branching ratios for the channel in which the excess was observed, i.e. $\sigma(gg \to A) \times \text{BR}(A \to ZH) \times \text{BR}(H \to t\bar{t})$. In the green shaded areas the predicted values are within the interval given in eq. (3.2), corresponding to a description of the excess within $1\sigma$.
In both plots one can see that for values of $1.2 \lesssim \tan\beta \lesssim 1.8$ and values of $\cos(\beta - \alpha)$ near the alignment limit a description of the excess is possible in accordance with the limits from searches for additional scalars and with the measurements of the properties of the detected Higgs boson at 125 GeV. Other decay channels for $H$, such as $H \to VV$ and $H \to hh$, become relevant outside of the alignment limit. As a result, for a fixed value of $\tan\beta$ the predicted cross sections decrease with increasing values of $|\cos(\beta - \alpha)|$. This gives rise to the shape of the green lines peaking at $\cos(\beta - \alpha) \approx 0$. Regarding a possible description of the excess in combination with deviations from the alignment limit, one can observe in the left plot that a sufficiently large cross section of more than about 0.30 pb can be achieved for $-0.15 \lesssim \cos(\beta - \alpha) \lesssim 0.10$, corresponding to modifications of the couplings of $h$ compared to the SM predictions that can exceed 10%. However, for both types values of $|\cos(\beta - \alpha)| \gtrsim 0.025$ are excluded by the cross section measurements of the SM-like Higgs boson and/or by the searches for $H \to hh \to b\bar{b}\tau^+\tau^-$. In type II the exclusion regions of both these constraints are largely overlapping, as shown in the right plot of figure 4, where the condition $\Delta\chi^2_{125} < 6.18$ yields an exclusion for values of $\cos(\beta - \alpha) \lesssim -0.02$ and $\cos(\beta - \alpha) \gtrsim 0.045$, while the searches for $H \to hh$ exclude $|\cos(\beta - \alpha)| \gtrsim 0.025$. The situation is different in type I, where for negative values of $\cos(\beta - \alpha)$ the parameter region excluded by **HiggsSignals** features values of $\cos(\beta - \alpha) \lesssim -0.15$, whereas the $H \to hh$ searches exclude already values of $\cos(\beta - \alpha) \lesssim -0.025$ (as in type II). Overall, we find that the observed excess can readily be accommodated by the 2HDM within the level of $1\sigma$ while being in agreement with all cross-section limits from BSM scalar searches and the measurements of properties of the detected Higgs boson at 125 GeV.

### 3.4.2 Gravitational wave detection

In the discussion in section 3.2.1 we already pointed out that the considered benchmark point predicts a FOEWPT according to the thermal effective potential approach. The fact that the excess can be readily accommodated in the 2HDM, as discussed above, motivates a closer look at the FOEWPT and related phenomenological consequences. More concretely, we analyze whether the stochastic GW background predicted by a parameter point compatible with the excess could potentially be observed at future experiments, in particular by LISA.

To this end, we performed a dedicated scan in $m_A$ and $m_H$ within a $\pm 50$ GeV mass window around the values of $m_A = 650$ GeV and $m_H = 450$ GeV, fixing the other 2HDM parameters as discussed above. The prediction for the GW signal was calculated for all the parameter points featuring a FOEWPT. In the 2HDM the GW signals produced in FOEWPTs are very sensitive to the precise values of the scalar masses [24]. Notably, a mere 50 GeV variation in the scalar masses leads to a SNR spanning many orders of magnitude, as can be seen in figure 5, where we show the parameter points featuring a FOEWPT of the dedicated scan in the $(m_H, m_A)$-plane. The colors of the points indicate the SNR of the GW signal at LISA. For the computation of the SNR, we assume a bubble wall velocity of $v_w = 0.6$, and the LISA operation time is assumed to be seven years, see section 2.5 for details. We only depict parameter points for which SNR $> 10^{-11}$. One can see that the predicted values of the SNRs vary over ten orders of magnitude within the relatively small mass window considered here. If we consider as detectable GW signals the ones with
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The signal-to-noise ratio (SNR) of the scanned parameter points in the $(m_H, m_A)$-plane for $\tan \beta = 1.5$. The mass values for which the excess observed by ATLAS in the $\ell^+\ell^- t\bar{t}$ final state was most pronounced (2.85 $\sigma$ local significance) are indicated with a magenta star. Shown are parameter points that feature a FOEWPT and a predicted GW signal with a SNR at LISA that is larger than $10^{-11}$. We regard the region within the magenta dashed contour as a rough estimate of the part of the parameter space that is compatible with the description of the excess in view of the mass resolution of the ATLAS search. The expected and observed exclusion limits are indicated by the red dashed and solid lines, respectively. The colors of the points display the SNRs of the GW signals at LISA assuming an operation time of seven years and a bubble-wall velocity of $v_w = 0.6$.</figcaption></figure>

SNR $\gtrsim 1$, the parameter space that can be probed with LISA is confined to the small region highlighted by the bright points. The parameter region featuring the strongest FOEWPTs, which is indicated by the yellow points, is cut off by the onset of the phenomenon of vacuum trapping, see the discussion in ref. [24]. On the other hand, in the lower right area of the magenta dashed contour no points are present because this parameter region gives rise to either a very weak FOEWPT yielding GW signals that are far out of reach of LISA or does not feature a first-order phase transition at all.

As a consequence of the limited mass resolution of the observed excess and the very sensitive dependence of the GW signals on the scalar masses, no definitive conclusion can be drawn about whether the 2HDM interpretation of the excess would be associated with a primordial GW signal detectable with LISA. We emphasize that this statement holds irrespective of the substantial theoretical uncertainties present in the computation of the phase transition parameters and the GW signals (see refs. [52, 112, 113] for detailed discussions of the theoretical uncertainties). Even if the theoretical uncertainties from unknown higher-order contributions in the calculation of the SNR were negligible, the parametric uncertainties arising from the experimental error in the determination of the masses of the BSM Higgs bosons would be a limiting factor for assessing the question whether such an excess would
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Gravitational wave spectra for parameter points specified in table 2 that are compatible with the excess observed in the ATLAS search. The solid (dashed) lines show the prediction without (including) the turbulence contribution, using $v_w = 0.6$. The colored regions show the prospective sensitivities of future experiments.</figcaption></figure>

<table><caption>**Table 2.** Results for parameters characterizing the phase transition for example points of the 2HDM that are compatible with the excess observed in the ATLAS search. The corresponding GW spectra are shown in figure 6. Dimensionful parameters are given in GeV. The SNR values evaluated for LISA include the turbulence contribution.</caption><thead><tr><th>m<sub>H</sub></th><th>m<sub>A</sub></th><th>m<sub>A</sub> - m<sub>H</sub></th><th>T<sub>n</sub></th><th>v<sub>n</sub>/T<sub>n</sub></th><th>α</th><th>β/H</th><th>SNR</th></tr></thead><tbody><tr><td>417.2</td><td>660.8</td><td>243.6</td><td>79.44</td><td>3.10</td><td>0.0308</td><td>77</td><td>13.7</td></tr><tr><td>432.8</td><td>673.4</td><td>240.6</td><td>86.23</td><td>2.85</td><td>0.0206</td><td>134</td><td>3.8</td></tr><tr><td>453.6</td><td>686.4</td><td>232.8</td><td>110.89</td><td>2.19</td><td>0.0073</td><td>468</td><td>0.022</td></tr><tr><td>445.7</td><td>677.1</td><td>231.4</td><td>116.48</td><td>2.06</td><td>0.0062</td><td>674</td><td>0.004</td></tr><tr><td>450.0</td><td>650.0</td><td>200.0</td><td>145.08</td><td>1.42</td><td>0.0029</td><td>5399</td><td>&lt; 0.001</td></tr></tbody></table>

lead to GW signals detectable with LISA or other future GW detectors. We regard this finding as a generic feature of the interplay between the LHC and GW detectors: while in the case of the absence of a signal of BSM physics at the LHC the resulting limits will place strong and definitive constraints on the possibility of a detection of a GW signal produced in a FOEWPT within the considered class of models, a possible observation of BSM scalars may not provide sufficient information to make a clear prediction on whether a GW detection at LISA can be expected.

To further illustrate the impact of the experimental mass resolution of BSM scalar searches at the LHC on the predicted GW signals, we show in figure 6 the spectral shape of the GW backgrounds produced during a FOEWPT for several parameter points with masses of the heavy scalars specified in table 2 together with the parameters that characterize the phase transition. The remaining 2HDM parameters are kept fixed according to the previous discussion. We chose the point with the largest SNR found in figure 5 and allow for up to 10% deviations in the values of the masses $m_H, m_A$, which translates into deviations of the
SNR of several orders of magnitude. In addition, we show in table 2 the parameters for the point $(m_H, m_A) = (450, 650)$ GeV although we omit its GW spectrum in figure 6 because of the smallness of the SNR. The spectral shapes of the GW backgrounds are computed as discussed in section 2.5, where the solid curves depict the sound-wave contribution $h^2\Omega_{sw}$ only, whereas the dashed curves depict the sum of sound-wave and turbulence contributions, i.e. $h^2\Omega_{sw} + h^2\Omega_{turb}$. We also show the sensitivity curves of LISA [18], AEDGE [114], DECIGO [115, 116] and BBO [117], where the latter three are planned, but not yet approved space-based GW detectors. One can see that only for the smallest value of $m_H = 417.2$ GeV, i.e. the largest mass splitting between $H$ and $A$, the GW signal might be detectable with LISA, according to the predicted SNR. For values of $m_H$ only a few percent larger, the peak amplitudes of the GW signals drastically decrease and quickly drop to values far below the experimental sensitivity of the proposed GW detectors. We emphasize again at this point that the detectability of the GW signal for a single parameter point cannot be determined definitively with the methods applied here due to the substantial theoretical uncertainties in the prediction of the GW signals. However, the fact that in the case of a possible detection of BSM scalars at the LHC a mass resolution at the percent level would be required in order to draw conclusions about the detectability of a GW signal poses a challenge independently of the status of the remaining theoretical uncertainties at that time.

Of course, one can also turn this argument around. An LHC discovery, e.g. a signal in the smoking-gun signature, in combination with a GW detection at LISA that is consistent with a FOEWPT as interpreted in a UV-complete model, could be used for a more precise (but model-dependent) determination of the parameters of the considered BSM Higgs sector. In this way space-based GW astronomy could become a complementary tool to sharpen the precision of particle physics.<fnref n="14" />

## 4 Summary and conclusions

Recently, ATLAS has reported for the first time, and based on the full Run 2 data set collected at 13 TeV, the results for searches for additional Higgs bosons where a heavier (pseudo-)scalar resonance $A$ is produced via gluon-fusion and subsequently decays into a $Z$ boson and a lighter scalar resonance $H$. The search made use of the leptonic decay of the $Z$ boson, while for the lighter scalar the search focused on the decay into a top-quark pair. This signature is exceptionally promising for probing the 2HDM parameter space for the case where the masses of the neutral BSM scalars are above the di-top threshold and have a splitting that is at least as large as the mass of the $Z$ boson. Consequently, within the 2HDM this signature has been identified as a smoking-gun signature for a FOEWPT, whose presence relies on sizable mass splittings in order to generate a potential barrier separating the symmetry-conserving and the symmetry-breaking vacua. In this context in particular the region of low $\tan\beta$ is of interest, which is preferred by EW baryogenesis. Since the searches in the $\ell^+\ell^- t\bar{t}$ final state are able to probe parameter space regions of the 2HDM that were unconstrained up to now, we have performed a comprehensive analysis of the impact of the

<footnote n="14">¹⁴This would be similar in spirit to the present situation regarding the sum of neutrino masses, constrained most stringently using astrophysical observations, e.g. the measurement of the spectrum of the cosmic microwave background [118].</footnote>
cross-section limits reported by ATLAS. We focused on 2HDM benchmark scenarios assuming the alignment limit, $\cos(\beta - \alpha) = 0$, mass degeneracy between the charged Higgs boson and the pseudoscalar state, $m_{H^{\pm}} = m_A$, and setting the decoupling scale $M$ (defined by the relation $m_{12}^2 = M^2 \sin\beta \cos\beta$) equal to the mass of the heavier CP-even scalar, $M = m_H$.

In the first part of our analysis, we determined the parameter regions that are excluded by this new search at the 95% C.L. in the $(m_H, m_A)$ plane. We started by considering a low-tan $\beta$ regime with $1 \le \tan\beta \le 3$, preferred in view of a possible realization of EW baryogenesis and associated with relatively large couplings of $A$ and $H$ to top quarks, implying that the new ATLAS $\ell^+\ell^- t\bar{t}$ has high sensitivity. We found that for $\tan\beta = 1$ the new search excludes a large region of parameter space that so far was not constrained by LHC searches. In combination with LHC limits from searches for $H^{\pm} \to tb, H \to \tau^+\tau^-$ and $H, A \to t\bar{t}$, masses of 300 GeV $\lesssim m_H \lesssim 450$ GeV are now entirely excluded in this scenario, while previously for $\tan\beta = 1$ a wide parameter space region with 350 GeV $\lesssim m_H \lesssim 460$ GeV and 650 GeV $\lesssim m_A \lesssim 800$ GeV was left unconstrained. For increasing values of $\tan\beta$, the charged scalar searches lose sensitivity, and in particular for $\tan\beta = 1.5$ the $\ell^+\ell^- t\bar{t}$ search is the only LHC search that can currently probe the parameter space regions featuring a FOEWPT above the di-top threshold, $m_H > 2m_t$. For $\tan\beta = 2$, and irrespectively of whether the presence of a FOEWPT is required or not, the search for the smoking gun signature $A \to ZH$ with $H \to t\bar{t}$ is currently the only LHC search that is able to exclude parameter space regions in our benchmark plane with $m_H$ above the di-top threshold. For $\tan\beta = 3$, the largest value considered in the low-tan $\beta$ regime, the $\ell^+\ell^- t\bar{t}$ searches are currently not yet able to probe the 2HDM parameter plane, because the gluon-fusion production cross section of $A$ is too small. Instead, we demonstrated that searches for $A \to ZH$ with $H \to b\bar{b}$, which had previously been carried out by both ATLAS and CMS including the full Run 2 datasets, are more promising in this scenario, giving rise to an exclusion region reaching masses of up to $m_H \approx 400$ GeV. We have furthermore pointed out a novel, not yet performed LHC search that would be complementary to the smoking gun $A \to ZH$ search in the low $\tan\beta$ region as a probe of the 2HDM parameter region featuring a strong FOEWPT. The channel that we propose for experimental analysis consists of $H^{\pm}$ production (e.g. via $pp \to H^{\pm tb}$) followed by the decay $H^{\pm} \to W^{\pm}H \to \ell^{\pm}\nu t\bar{t}$.

In addition to the searches in the $\ell^+\ell^- t\bar{t}$ final state, ATLAS also reported for the first time searches for the $A \to ZH$ decay making use of the decay of the $Z$ boson into neutrino pairs and the $H \to b\bar{b}$ decay mode. Here, both the gluon-fusion and the $b\bar{b}$-associated production of $A$ were considered. These searches in the $\nu\bar{\nu}b\bar{b}$ final state may become important for large values of $\tan\beta$. Investigating a representative 2HDM benchmark scenario with $\tan\beta = 10$ (and the remaining 2HDM parameters as described above), we found that in the 2HDM the earlier searches in the $\ell^+\ell^- b\bar{b}$ final state give rise to stronger exclusions than the new $\nu\bar{\nu}b\bar{b}$ ATLAS search. Nevertheless, it should be taken into account that the new searches utilizing the $Z \to \nu\bar{\nu}$ decay mode cover a wider mass interval and larger mass splittings between the two BSM resonances. In the 2HDM the maximum amount of mass splittings between the BSM scalars is limited by the perturbativity constraint: the scalars are contained in the same SU(2) doublet, such that their masses are confined to lie not too far away from the overall decoupling scale $M$. In other models, in which additional mass scales are present and larger
mass splittings can be realized between different BSM scalars, the new searches in the $\nu\bar{\nu}b\bar{b}$ final state could be able to probe so far unconstrained parameter space regions.

As a further part of our analysis, motivated by the strong impact of the $\ell^+\ell^-t\bar{t}$ searches described above, we investigated in section 3.3 the future prospects in terms of projected exclusion regions, which were obtained via a simple rescaling of the reported expected cross-section limits of the new ATLAS search with integrated luminosities anticipated to be collected in the future at the (HL-)LHC. We found that the reach of the smoking gun signature will significantly improve and parameter regions with $m_H > 2m_t$ and $\tan\beta \ge 3$ will become accessible. We therefore anticipate that the smoking gun signature $A \to ZH$ with $H \to t\bar{t}$ decay will be the main LHC search channel in the future to probe the allowed 2HDM parameter space regions above the di-top threshold that feature a strong FOEWPT. Besides, given that successful EW baryogenesis prefers small $\tan\beta$ values (not much larger than one) [119], such searches will have a large impact on the possibility of explaining the matter-antimatter asymmetry via EW baryogenesis in the 2HDM.

The new ATLAS search in the $\ell^+\ell^-t\bar{t}$ final state showed an excess which is most significant for masses of $m_A = 650$ GeV and $m_H = 450$ GeV, with a local significance of $2.85\sigma$. We have demonstrated that the excess can be described at the $1\sigma$ level in the 2HDM type I and II (as representative scenarios) in the approximate range of $\tan\beta$ values between 1.2 and 1.8 for the alignment limit (a $1\sigma$ description of the excess is also possible for small departures from the alignment limit, with the viable $\tan\beta$ range shrinking accordingly), while being in agreement with all existing bounds from BSM scalar searches at the LHC. Further probes of the nature of the excess could be performed by searching for the decay $H^\pm \to W^\pm H$ followed by $H \to t\bar{t}$, as discussed above. Notably, the masses corresponding to the excess lie within the parameter space region indicative of a strong FOEWPT based on the thermal effective potential approach applied here. This parameter region could thus be suitable for successful EW baryogenesis. We emphasize, however, that in order to investigate whether the BAU can be predicted in agreement with observations one would have to take into account additional sources of CP-violation, required for the generation of the BAU according to the Sakharov conditions. These sources of CP-violation might have an impact on the cross section for the process in which the excess was observed. We leave a more detailed investigation of the predicted BAU and of the possible impact of new sources of CP-violation on the description of the excess for future studies. We analyzed the primordial GW signal that would be generated during the FOEWPT in the 2HDM parameter space region compatible with the observed excess. We found that, since the predictions for the GW spectra are highly sensitive to the precise values of the BSM scalar masses, the signal-to-noise ratio expected at LISA varies by

several order of magnitude for points within the region compatible with the ATLAS excess.<fnref n="15" /> Therefore at this stage no definitive statement can be made about whether the GW backgrounds would be detectable at LISA (or other future space-based GW detectors). Nevertheless, should a stochastic GW signal be detected by LISA, in combination with an

<footnote n="15">¹⁵We stress again that this is purely a consequence of the parametric uncertainty stemming from the experimental mass resolution of BSM scalar searches at the LHC, which will pose a challenge independently of whether the theoretical uncertainties on the predictions for the GW power spectra arising from higher-order effects can significantly be reduced.</footnote>
LHC signal as hinted by the ATLAS excess this would within the context of the 2HDM provide new and very precise information on the allowed values of the parameters.

In order to produce the results presented in this paper, we developed the software package `thdmTools`, which can be used for the phenomenological investigation of the CP-conserving 2HDM with softly-broken $Z_2$ symmetry. Accompanying this paper, we make the code `thdmTools` available to the public. Installation instructions and a brief discussion about the functionalities of the code and its interfaces to other public codes, in particular to `HiggsTools` are given in appendix A.

Our analysis shows that the smoking gun signature $A \to ZH$ (with $H \to t\bar{t}$) has great potential to further probe the viable 2HDM parameter regions, in particular those that may feature a strong FOEWPT as required for EW baryogenesis. In fact, the excess that has been observed in this search could be the first hint of such a strong transition in the 2HDM.

## Acknowledgments

We thank Philipp Gadow, Daniel Hundhausen, Benoit Laurent, Matthias Schröder and Matthias Steinhauser for helpful discussions. The work of T.B. is supported by the German Bundesministerium für Bildung und Forschung (BMBF, Federal Ministry of Education and Research) — project 05H21VKCCA. The work of M.O.O.R is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 101002846, ERC CoG "CosmoChart". J.M.N. was supported by the Ramón y Cajal Fellowship contract RYC-2017-22986, and acknowledges partial financial support by the grant PID2021-124704NB-I00, and by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements No 860881-HIDDeN and 101086085-ASYMMETRY. S.H. and J.M.N. acknowledge partial financial support by the Spanish Research Agency (Agencia Estatal de Investigación) through the grant IFT Centro de ExcePCIencia Severo Ochoa No CEX2020-001007-S funded by MCIN/AEI/10.13039/501100011033. The work of S.H. was supported in part by the grant PID2019-110058GB-C21 funded by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe", and in part by the Grant PID2022-142545NB-C21 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. K.R. and G.W acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2121 "Quantum Universe" - 390833306. This work has been partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 491245950.

## A The python package thdmTools

`thdmTools` is a python package for exploring the Two Higgs Doublet Model with real parameters and a softly-broken $Z_2$ symmetry, which incorporates tests of the relevant theoretical and experimental constraints. It allows the user to specify a parameter point in terms of the free parameters of the model (see eq. (2.3)). During the installation of this package the following external codes will also be downloaded and installed:

*AnyHdecay* [83–85]: Computes the branching ratios and decay widths of all Higgs bosons contained in the model
**HiggsTools** [77]: Checks compatibility with the experimental constraints on BSM scalars from LEP and LHC searches (HiggsBounds [78–81]) and with existing measurements of the detected Higgs boson at about 125 GeV (HiggsSignals [86, 87, 109]).

**THDM_EWPOS** [120–122]:Computes the prediction for some Electroweak Precision Observables (EWPO), in particular $M_W$ and the effective weak mixing angle at the Z-boson resonance, in the 2HDM at the two-loop level.

**SuperIso** [123, 124]:Computes the predictions for various flavour-physics observables in the 2HDM.

Note that HiggsTools will only be installed if not already installed in the current python environment of the user.

## A.1 Installation

thdmTools is publicly available at: https://gitlab.com/thdmtools. Installation requires a python3 environment and compilers for Fortran and C++. The package and all its dependencies are installed by executing:

```
make all
```

## A.2 Code example for thdmTools

To import the package in an example notebook type

```
from thdmTools.parampoint import Parampoint
```

To define the parameters of the point as a dictionary in python

```python
dc = {
    'type': 1,
    'tb': 3,
    'alpha': -0.3217,
    'mHh': 500,
    'mHl': 125.09,
    'mA': 500,
    'mHp': 500,
    'm12sq': 75000
}
```

To run the code and check all the included theoretical and experimental constraints

```python
pt = Parampoint语言模型 dc)
print('Stability:', pt.check_vacuum_stability())
print('Unitarity:', pt.check_perturbative_unitarity())
print('Flavour:', pt.check_flavour_constraints())
print('EWPO:', pt.check_ewpo_constraints())
print('Collider (HB,HS):', pt.check_collider_constraints())
```
These commands will print a boolean statement of True/False if the point is allowed/disallowed by the corresponding constraint, see below for further details on the criteria applied in each case. Vacuum stability can be checked with the conditions on boundedness from below of the tree level potential, see e.g. ref. [125]. An additional condition that ensures that the EW minimum is the global minimum of the tree level potential [126] can be applied by setting the optional argument `global_min` of the function `check_vacuum_stability()` to `True`:

```python
print('Stability:', pt.check_vacuum_stability(global_min=True))
```

In order to test perturbative unitarity the default configuration checks for the unitarity of the scalar four-point scattering matrix in the high-energy limit, applying an upper limit of $8\pi$ to the eigenvalues of the scattering matrix. One can change the upper limit by giving an argument in the function

```python
('Unitarity to 1:', pt.check_perturbative_unitarity(1.))
```

where the value of the argument multiplied by $16\pi$ is the applied upper limit (thus 0.5 is the default value). One can also check for NLO perturbative unitarity in the high energy limit according to the expressions derived in ref. [127]:

```python
('NLO Unitarity:', pt.check_perturbative_unitarity_nlo())
('NLO Unitarity to 1:', pt.check_perturbative_unitarity_nlo(1.))
```

Regarding constraints from flavor physics the predictions for $B \to s\gamma$ and $B_s \to \mu^+\mu^-$ are checked. One can access the specified allowed region by typing

```python
print(pt.b2sgam_valid)
print(pt.bs2mumu_valid)
```

Many other predictions for flavor-physics observables are available via the interface to SuperIso, such that users can define their own functions to exclude or accept parameter points including additional observables.

The function called above to check against constraints from EWPO calls the interface to THDM_EWPOS. This function in particular verifies whether the predicted values for the W-boson mass, the total decay width of the Z boson and the effective weak mixing angle at the Z-boson resonance are in agreement with the experimental measurements (by default the prediction for $M_W$ is checked against the average value from the LHC-TeV $M_W$ Working Group [27], which does not include the CDF $M_W$ result) within two standard deviations. The experimental values and their uncertainties can be changed by the user via optional arguments to the function `check_ewpo_constraints()`. As a further possibility, the user can perform a $\chi^2$ fit to the experimental EWPO in terms of the oblique parameters $S$, $T$, and $U$ by typing

```python
print('STU chi^2 fit:' pt.check_ewpo_fit())
```

In contrast to the analysis of the EWPO as specified above, the $S$, $T$, $U$ parameters are evaluated only at the one-loop level according to ref. [128]. The experimental fit values
of the oblique parameters, their uncertainties and correlations are taken from ref. [129]. By default, a parameter point is considered to be viable according to this approach if the predicted values of the S, T, U parameters are in agreement with the experimental fit values within the $2\sigma$ confidence level.

Regarding the collider constraints, for compatibility with the 125 GeV Higgs-boson measurements one can specify the allowed range of $\Delta\chi^2_{125}$ to perform a test using the best $\chi^2$ fit performed with `HiggsSignals`. The default is to allow parameter points with $\chi^2_{125}$-values that are not more than $2\sigma$ away from the SM $\chi^2_{125}$-value, which corresponds to $\Delta\chi^2_{125} = 6.18$ in two-dimensional parameter representations. For compatibility with cross-section limits from searches for additional scalars, one can call the `HiggsBounds` interface and access the result by means of

```python
print('HiggsBounds result:', pt.reshb)
```

which will print a table with the observed/expected ratios of the most sensitive channels for each of the scalars in the 2HDM. The full information of the `HiggsBounds` analysis is stored in the object `pt.reshb` and can be accessed using the various functionalities of `HiggsTools` [77]. The cross-section predictions from `HiggsTools` for the LHC operating at 13 TeV (based on the effective coupling input) can also be accessed by typing, e.g.

```python
print('sigma(gg -> h)', pt.XS['H1'] ['gg'])
```

where `H1` refers to the lighter CP-even scalar and `gg` selects the gluon-fusion production cross section. The cross sections of the second CP-even scalar, the CP-odd scalar and the charged scalars can be chosen by using the keys `Hh`, `A` and `Hp`, respectively. For the neutral states, the other available production modes stored in `pt.XS` are $b\bar{b}$-associated production, vector-boson fusion, $t\bar{t}$-associated production and $t$-associated production, which can be accessed with the keys `bbH`, `vbfH`, `ttH` and `tH`, respectively. Finally, it is also possible to access all branching ratios and total decay widths of the particles by typing

```python
pt.calibrateolation_ratios()
print('BR_h', pt.b_Hl)
print('BR_H', pt.b_Hh)
print('BR_A', pt.b_A)
print('BR_Hp', pt.b_Hp)
print('Total decay width h:', pt.wTot['Hl'])
print('Total decay width H:', pt.wTot['Hh'])
print('Total decay width A:', pt.wTot['A'])
print('Total decay width Hp:', pt.wTot['Hp'])
```

where the predictions for the branching ratios are computed via the interface to the `AnyHdecay` library.
References

[1] ATLAS collaboration, Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716 (2012) 1 [arXiv:1207.7214] [INSPIRE].

[2] CMS collaboration, Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC, Phys. Lett. B 716 (2012) 30 [arXiv:1207.7235] [INSPIRE].

[3] CMS collaboration, A portrait of the Higgs boson by the CMS experiment ten years after the discovery, Nature 607 (2022) 60 [arXiv:2207.00043] [INSPIRE].

[4] ATLAS collaboration, A detailed map of Higgs boson interactions by the ATLAS experiment ten years after the discovery, Nature 607 (2022) 52 [Erratum ibid. 612 (2022) E24] [arXiv:2207.00092] [INSPIRE].

[5] ATLAS collaboration, Standard Model Summary Plots February 2022, ATL-PHYS-PUB-2022-009, CERN, Geneva (2022).

[6] K. Kajantie, M. Laine, K. Rummukainen and M.E. Shaposhnikov, Is there a hot electroweak phase transition at $m_H \gtrsim m_W?$, Phys. Rev. Lett. 77 (1996) 2887 [hep-ph/9605288] [INSPIRE].

[7] V.A. Kuzmin, V.A. Rubakov and M.E. Shaposhnikov, On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe, Phys. Lett. B 155 (1985) 36 [INSPIRE].

[8] T.D. Lee, A Theory of Spontaneous T Violation, Phys. Rev. D 8 (1973) 1226 [INSPIRE].

[9] J.E. Kim, Weak Interaction Singlet and Strong CP Invariance, Phys. Rev. Lett. 43 (1979) 103 [INSPIRE].

[10] F. Wilczek, Problem of Strong P and T Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279 [INSPIRE].

[11] J.M. Cline and P.-A. Lemieux, Electroweak phase transition in two Higgs doublet models, Phys. Rev. D 55 (1997) 3873 [hep-ph/9609240] [INSPIRE].

[12] L. Fromme, S.J. Huber and M. Seniuch, Baryogenesis in the two-Higgs doublet model, JHEP 11 (2006) 038 [hep-ph/0605242] [INSPIRE].

[13] J.M. Cline, K. Kainulainen and M. Trott, Electroweak Baryogenesis in Two Higgs Doublet Models and B meson anomalies, JHEP 11 (2011) 089 [arXiv:1107.3559] [INSPIRE].

[14] G.C. Dorsch, S.J. Huber, T. Konstandin and J.M. No, A Second Higgs Doublet in the Early Universe: baryogenesis and Gravitational Waves, JCAP 05 (2017) 052 [arXiv:1611.05874] [INSPIRE].

[15] A. Ashoorioon and T. Konstandin, Strong electroweak phase transitions without collider traces, JHEP 07 (2009) 086 [arXiv:0904.0353] [INSPIRE].

[16] D. Curtin, P. Meade and C.-T. Yu, Testing Electroweak Baryogenesis with Future Colliders, JHEP 11 (2014) 127 [arXiv:1409.0005] [INSPIRE].

[17] C. Caprini et al., Detecting gravitational waves from cosmological phase transitions with LISA: an update, JCAP 03 (2020) 024 [arXiv:1910.13125] [INSPIRE].

[18] LISA COSMOLOGY WORKING GROUP collaboration, Cosmology with the Laser Interferometer Space Antenna, Living Rev. Rel. 26 (2023) 5 [arXiv:2204.05434] [INSPIRE].

[19] G.C. Dorsch, S.J. Huber and J.M. No, A strong electroweak phase transition in the 2HDM after LHC8, JHEP 10 (2013) 029 [arXiv:1305.6610] [INSPIRE].
[20] G.C. Dorsch, S.J. Huber, K. Mimasu and J.M. No, *Echoes of the Electroweak Phase Transition: discovering a second Higgs doublet through $A_0 \to ZH_0$, Phys. Rev. Lett. 113 (2014) 211802* [arXiv:1405.5537] [INSPIRE].

[21] P. Basler et al., *Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited*, JHEP 02 (2017) 121 [arXiv:1612.04086] [INSPIRE].

[22] G.C. Dorsch, S.J. Huber, K. Mimasu and J.M. No, *The Higgs Vacuum Uplifted: Revisiting the Electroweak Phase Transition with a Second Higgs Doublet*, JHEP 12 (2017) 086* [arXiv:1705.09186] [INSPIRE].

[23] D. Gonçalves, A. Kaladharan and Y. Wu, *Electroweak phase transition in the 2HDM: collider and gravitational wave complementarity*, Phys. Rev. D 105 (2022) 095041 [arXiv:2108.05356] [INSPIRE].

[24] T. Biekötter et al., *The trap in the early Universe: impact on the interplay between gravitational waves and LHC physics in the 2HDM*, JCAP 03 (2023) 031 [arXiv:2208.14466] [INSPIRE].

[25] T. Biekötter et al., *Fate of electroweak symmetry in the early Universe: non-restoration and trapped vacua in the N2HDM*, JCAP 06 (2021) 018 [arXiv:2103.12707] [INSPIRE].

[26] J.F. Gunion, S. Dawson, H.E. Haber and G.L. Kane, *The Higgs hunter's guide*, SCIPP-89/13, UCD-89-4, BNL-41644, Brookhaven Nat. Lab., Upton, NY (1989).

[27] S. Amoroso et al., *Compatibility and combination of world W-boson mass measurements*, arXiv:2308.09417 [INSPIRE].

[28] C.-T. Lu, L. Wu, Y. Wu and B. Zhu, *Electroweak precision fit and new physics in light of the W boson mass*, Phys. Rev. D 106 (2022) 035034 [arXiv:2204.03796] [INSPIRE].

[29] H. Bahl, J. Braathen and G. Weiglein, *New physics effects on the W-boson mass from a doublet extension of the SM Higgs sector*, Phys. Lett. B 833 (2022) 137295 [arXiv:2204.05269] [INSPIRE].

[30] F. Arco, S. Heinemeyer and M.J. Herrero, *Sensitivity and constraints to the 2HDM soft-breaking $Z_2$ parameter $m_{12}$*, Phys. Lett. B 835 (2022) 137548 [arXiv:2207.13501] [INSPIRE].

[31] HFLAV collaboration, *Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2021*, Phys. Rev. D 107 (2023) 052008 [arXiv:2206.07501] [INSPIRE].

[32] M. Misiak, A. Rehman and M. Steinhauser, *Towards $\bar{B} \to X_s \gamma$ at the NNLO in QCD without interpolation in $m_c$*, JHEP 06 (2020) 175 [arXiv:2002.01548] [INSPIRE].

[33] M. Steinhauser, private communication.

[34] BELLE-II collaboration, *Measurement of the photon-energy spectrum in inclusive $B \to X_s \gamma$ decays identified using hadronic decays of the recoil B meson in 2019-2021 Belle II data*, arXiv:2210.10220 [INSPIRE].

[35] CMS collaboration, *Search for neutral resonances decaying into a Z boson and a pair of b jets or $\tau$ leptons*, Phys. Lett. B 759 (2016) 369 [arXiv:1603.02991] [INSPIRE].

[36] ATLAS collaboration, *Search for a heavy Higgs boson decaying into a Z boson and another heavy Higgs boson in the $l\bar{b}b$ final state in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector*, Phys. Lett. B 783 (2018) 392 [arXiv:1804.01126] [INSPIRE].

[37] CMS collaboration, *Search for new neutral Higgs bosons through the $H \to ZA \to l^+l^-b\bar{b}$ process in pp collisions at $\sqrt{s} = 13$ TeV*, JHEP 03 (2020) 055 [arXiv:1911.03781] [INSPIRE].
[38] ATLAS collaboration, *Search for a heavy Higgs boson decaying into a Z boson and another heavy Higgs boson in the $ℓℓb\bar{b}$ and $ℓℓWW$ final states in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector*, Eur. Phys. J. C 81 (2021) 396 [arXiv:2011.05639] [INSPIRE].

[39] D. Gonçalves, A. Kaladharan and Y. Wu, *Resonant top pair searches at the LHC: a window to the electroweak phase transition*, Phys. Rev. D 107 (2023) 075040 [arXiv:2206.08381] [INSPIRE].

[40] R. Mammen Abraham and D. Gonçalves, *Boosting new physics searches in $t\bar{t}Z$ and $tZj$ production with angular moments*, Eur. Phys. J. C 83 (2023) 965 [arXiv:2208.05986] [INSPIRE].

[41] ATLAS collaboration, *Search for a CP-odd Higgs boson decaying to a heavy CP-even Higgs boson and a Z boson in the $ℓℓt\bar{t}$ and $\nu\bar{\nu}b\bar{b}$ final states using 140 fb$^{-1}$ of data collected with the ATLAS detector*, ATLAS-CONF-2023-034, CERN, Geneva (2023).

[42] D. Hundhausen et al., *Search for heavy Higgs bosons in the $Z + t\bar{t}$ final state with CMS*, contribution to DPG Spring Meeting, https://www.dpg-verhandlungen.de/year/2022/conferen ce/heidelberg/part/t/session/62/contribution/7

[43] Y. Fischer, *Search for heavy Higgs bosons in the $Z + t\bar{t}$ final state with the CMS detector*, Master's thesis, Universität Hamburg (2021), https://www.physik.uni-hamburg.de/iexp/grupp e-haller/scientific-output/documents/masterarbeit-yannick-fischer.pdf.

[44] T. Biekötter and M.O. Olea-Romacho, *Reconciling Higgs physics and pseudo-Nambu-Goldstone dark matter in the S2HDM using a genetic algorithm*, JHEP 10 (2021) 215 [arXiv:2108.10864] [INSPIRE].

[45] G.C. Branco et al., *Theory and phenomenology of two-Higgs-doublet models*, Phys. Rept. 516 (2012) 1 [arXiv:1106.0034] [INSPIRE].

[46] J.F. Gunion and H.E. Haber, *The CP conserving two Higgs doublet model: the approach to the decoupling limit*, Phys. Rev. D 67 (2003) 075019 [hep-ph/0207010] [INSPIRE].

[47] F. Arco, S. Heinemeyer and M.J. Herrero, *Triple Higgs couplings in the 2HDM: the complete picture*, Eur. Phys. J. C 82 (2022) 536 [arXiv:2203.12684] [INSPIRE].

[48] S.R. Coleman and E.J. Weinberg, *Radiative Corrections as the Origin of Spontaneous Symmetry Breaking*, Phys. Rev. D 7 (1973) 1888 [INSPIRE].

[49] M. Quiros, *Finite temperature field theory and phase transitions*, in the proceedings of the ICTP Summer School in High-Energy Physics and Cosmology, Miramare, Italy, June 29 – July 17 (1998), p. 187–259 [hep-ph/9901312] [INSPIRE].

[50] P.B. Arnold and O. Espinosa, *The effective potential and first order phase transitions: beyond leading-order*, Phys. Rev. D 47 (1993) 3546 [Erratum ibid. 50 (1994) 6662] [hep-ph/9212235] [INSPIRE].

[51] D. Curtin, P. Meade and H. Ramani, *Thermal Resummation and Phase Transitions*, Eur. Phys. J. C 78 (2018) 787 [arXiv:1612.00466] [INSPIRE].

[52] D. Croon et al., *Theoretical uncertainties for cosmological first-order phase transitions*, JHEP 04 (2021) 055 [arXiv:2009.10080] [INSPIRE].

[53] O. Gould and T.V.I. Tenkanen, *On the perturbative expansion at high temperature and implications for cosmological phase transitions*, JHEP 06 (2021) 069 [arXiv:2104.04399] [INSPIRE].
[54] O. Gould and T.V.I. Tenkanen, *Perturbative effective field theory expansions for cosmological phase transitions*, arXiv:2309.01672 [INSPIRE].

[55] C.L. Wainwright, *CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields*, Comput. Phys. Commun. **183** (2012) 2006 [arXiv:1109.4189] [INSPIRE].

[56] S. Baum et al., *Nucleation is more than critical: a case study of the electroweak phase transition in the NMSSM*, JHEP **03** (2021) 055 [arXiv:2009.10743] [INSPIRE].

[57] A.D. Sakharov, *Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe*, Pisma Zh. Eksp. Teor. Fiz. **5** (1967) 32 [INSPIRE].

[58] S. Dimopoulos and L. Susskind, *On the Baryon Number of the Universe*, Phys. Rev. D **18** (1978) 4500 [INSPIRE].

[59] D. Gonçalves, A. Kaladharan and Y. Wu, *Gravitational waves, bubble profile, and baryon asymmetry in the complex 2HDM*, Phys. Rev. D **108** (2023) 075010 [arXiv:2307.03224] [INSPIRE].

[60] E. Witten, *Cosmic Separation of Phases*, Phys. Rev. D **30** (1984) 272 [INSPIRE].

[61] C.J. Hogan, *Gravitational radiation from cosmological phase transitions*, Mon. Not. Roy. Astron. Soc. **218** (1986) 629 [INSPIRE].

[62] LISA collaboration, *Laser Interferometer Space Antenna*, arXiv:1702.00786 [INSPIRE].

[63] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, *Gravitational waves from the sound of a first order phase transition*, Phys. Rev. Lett. **112** (2014) 041301 [arXiv:1304.2433] [INSPIRE].

[64] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, *Numerical simulations of acoustically generated gravitational waves at a first order phase transition*, Phys. Rev. D **92** (2015) 123009 [arXiv:1504.03291] [INSPIRE].

[65] M. Hindmarsh, *Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe*, Phys. Rev. Lett. **120** (2018) 071301 [arXiv:1608.04735] [INSPIRE].

[66] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, *Shape of the acoustic gravitational wave power spectrum from a first order phase transition*, Phys. Rev. D **96** (2017) 103520 [Erratum ibid. **101** (2020) 089902] [arXiv:1704.05871] [INSPIRE].

[67] C. Caprini and R. Durrer, *Gravitational waves from stochastic relativistic sources: primordial turbulence and magnetic fields*, Phys. Rev. D **74** (2006) 063521 [astro-ph/0603476] [INSPIRE].

[68] G. Gogoberidze, T. Kahniashvili and A. Kosowsky, *The Spectrum of Gravitational Radiation from Primordial Turbulence*, Phys. Rev. D **76** (2007) 083002 [arXiv:0705.1733] [INSPIRE].

[69] C. Caprini, R. Durrer and G. Servant, *The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition*, JCAP **12** (2009) 024 [arXiv:0909.0622] [INSPIRE].

[70] A. Roper Pol et al., *Numerical simulations of gravitational waves from early-universe turbulence*, Phys. Rev. D **102** (2020) 083512 [arXiv:1903.08585] [INSPIRE].

[71] P. Auclair et al., *Generation of gravitational waves from freely decaying turbulence*, JCAP **09** (2022) 029 [arXiv:2205.02588] [INSPIRE].
[72] J. Ellis, M. Lewicki and J.M. No, On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal, JCAP 04 (2019) 003 [arXiv:1809.0824] [INSPIRE].

[73] B. Laurent and J.M. Cline, Fluid equations for fast-moving electroweak bubble walls, Phys. Rev. D 102 (2020) 063516 [arXiv:2007.10935] [INSPIRE].

[74] W.-Y. Ai, B. Laurent and J. van de Vis, Model-independent bubble wall velocities in local thermal equilibrium, JCAP 07 (2023) 002 [arXiv:2303.10171] [INSPIRE].

[75] G.C. Dorsch, S.J. Huber and T. Konstandin, Bubble wall velocities in the Standard Model and beyond, JCAP 12 (2018) 034 [arXiv:1809.04907] [INSPIRE].

[76] LISA collaboration, Science Requirements Document, Tech. Rep. ESA-L3-EST-SCI-RS-001 (2018), https://www.cosmos.esa.int/web/lisa/lisa-documents.

[77] H. Bahl et al., HiggsTools: BSM scalar phenomenology with new versions of HiggsBounds and HiggsSignals, Comput. Phys. Commun. 291 (2023) 108803 [arXiv:2210.09332] [INSPIRE].

[78] P. Bechtle et al., HiggsBounds: Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron, Comput. Phys. Commun. 181 (2010) 138 [arXiv:0811.4169] [INSPIRE].

[79] P. Bechtle et al., HiggsBounds 2.0.0: Confronting Neutral and Charged Higgs Sector Predictions with Exclusion Bounds from LEP and the Tevatron, Comput. Phys. Commun. 182 (2011) 2605 [arXiv:1102.1898] [INSPIRE].

[80] P. Bechtle et al., HiggsBounds-4: Improved Tests of Extended Higgs Sectors against Exclusion Bounds from LEP, the Tevatron and the LHC, Eur. Phys. J. C 74 (2014) 2693 [arXiv:1311.0055] [INSPIRE].

[81] P. Bechtle et al., HiggsBounds-5: testing Higgs Sectors in the LHC 13 TeV Era, Eur. Phys. J. C 80 (2020) 1211 [arXiv:2006.06007] [INSPIRE].

[82] LHC HIGGS CROSS SECTION WORKING GROUP collaboration, Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector, arXiv:1610.07922 [DOI:10.23731/CYRM-2017-002] [INSPIRE].

[83] A. Djouadi, J. Kalinowski and M. Spira, HDECAY: a Program for Higgs boson decays in the standard model and its supersymmetric extension, Comput. Phys. Commun. 108 (1998) 56 [hep-ph/9704448] [INSPIRE].

[84] HDECAY collaboration, HDECAY: twenty++ years after, Comput. Phys. Commun. 238 (2019) 214 [arXiv:1801.09506] [INSPIRE].

[85] M. Muhlleitner, M.O.P. Sampaio, R. Santos and J. Wittbrodt, The N2HDM under Theoretical and Experimental Scrutiny, JHEP 03 (2017) 094 [arXiv:1612.01309] [INSPIRE].

[86] P. Bechtle et al., HiggsSignals: confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC, Eur. Phys. J. C 74 (2014) 2711 [arXiv:1305.1933] [INSPIRE].

[87] P. Bechtle et al., HiggsSignals-2: probing new physics with precision Higgs measurements in the LHC 13 TeV era, Eur. Phys. J. C 81 (2021) 145 [arXiv:2012.09197] [INSPIRE].

[88] CMS collaboration, Search for charged Higgs bosons decaying into a top and a bottom quark in the all-jet final state of pp collisions at $\sqrt{s} = 13$ TeV, JHEP 07 (2020) 126 [arXiv:2001.07763] [INSPIRE].

[89] ATLAS collaboration, Search for charged Higgs bosons decaying into a top quark and a bottom quark at $\sqrt{s} = 13$ TeV with the ATLAS detector, JHEP 06 (2021) 145 [arXiv:2102.10076] [INSPIRE].
[90] CMS collaboration, Search for heavy Higgs bosons decaying to a top quark pair in proton-proton collisions at $\sqrt{s} = 13$ TeV, JHEP 04 (2020) 171 [Erratum ibid. 03 (2022) 187] [arXiv:1908.01115] [INSPIRE].

[91] ATLAS collaboration, Search for heavy Higgs bosons decaying into two tau leptons with the ATLAS detector using pp collisions at $\sqrt{s} = 13$ TeV, Phys. Rev. Lett. 125 (2020) 051801 [arXiv:2002.12223] [INSPIRE].

[92] CMS collaboration, Searches for additional Higgs bosons and for vector leptoquarks in $\tau\tau$ final states in proton-proton collisions at $\sqrt{s} = 13$ TeV, JHEP 07 (2023) 073 [arXiv:2208.02717] [INSPIRE].

[93] CMS collaboration, Search for a charged Higgs boson decaying into a heavy neutral Higgs boson and a W boson in proton-proton collisions at $\sqrt{s} = 13$ TeV, JHEP 09 (2023) 032 [arXiv:2207.01046] [INSPIRE].

[94] CMS collaboration, Observation of four top quark production in proton-proton collisions at $\sqrt{s} = 13$ TeV, Phys. Lett. B 847 (2023) 138290 [arXiv:2305.13439] [INSPIRE].

[95] ATLAS collaboration, Observation of four-top-quark production in the multilepton final state with the ATLAS detector, Eur. Phys. J. C 83 (2023) 496 [arXiv:2303.15061] [INSPIRE].

[96] S. Gori, I.-W. Kim, N.R. Shah and K.M. Zurek, Closing the Wedge: Search Strategies for Extended Higgs Sectors with Heavy Flavor Final States, Phys. Rev. D 93 (2016) 075038 [arXiv:1602.02782] [INSPIRE].

[97] CMS collaboration, Search for dilepton resonances from decays of (pseudo)scalar bosons produced in association with a massive vector boson or top quark anti-top quark pair at $\sqrt{s} = 13$ TeV, CMS-PAS-EXO-21-018, CERN, Geneva (2022).

[98] ATLAS collaboration, Search for heavy neutral Higgs bosons produced in association with b-quarks and decaying into b-quarks at $\sqrt{s} = 13$ TeV with the ATLAS detector, Phys. Rev. D 102 (2020) 032004 [arXiv:1907.02749] [INSPIRE].

[99] ATLAS collaboration, Search for low-mass resonances decaying into two jets and produced in association with a photon using pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector, Phys. Lett. B 795 (2019) 56 [arXiv:1901.10917] [INSPIRE].

[100] K. Kainulainen et al., On the validity of perturbative studies of the electroweak phase transition in the Two Higgs Doublet model, JHEP 06 (2019) 075 [arXiv:1904.01329] [INSPIRE].

[101] A. Ekstedt, P. Schicho and T.V.I. Tenkanen, DRalgo: a package for effective field theory approach for thermal phase transitions, Comput. Phys. Commun. 288 (2023) 108725 [arXiv:2205.08815] [INSPIRE].

[102] L. Niemi, M.J. Ramsey-Musolf, T.V.I. Tenkanen and D.J. Weir, Thermodynamics of a Two-Step Electroweak Phase Transition, Phys. Rev. Lett. 126 (2021) 171802 [arXiv:2005.11332] [INSPIRE].

[103] P.M. Schicho, T.V.I. Tenkanen and J. Österman, Robust approach to thermal resummation: Standard Model meets a singlet, JHEP 06 (2021) 130 [arXiv:2102.11145] [INSPIRE].

[104] P. Schicho, T.V.I. Tenkanen and G. White, Combining thermal resummation and gauge invariance for electroweak phase transition, JHEP 11 (2022) 047 [arXiv:2203.04284] [INSPIRE].

[105] CMS collaboration, Measurement of the cross section of top quark-antiquark pair production in association with a W boson in proton-proton collisions at $\sqrt{s} = 13$ TeV, JHEP 07 (2023) 219 [arXiv:2208.06485] [INSPIRE].
[106] ATLAS collaboration, Measurement of the total and differential cross-sections of $t\bar{t}W$ production in pp collisions at 13 TeV with the ATLAS detector, ATLAS-CONF-2023-019, CERN, Geneva (2023).

[107] ATLAS collaboration, Measurement of Higgs boson decay into b-quarks in associated production with a top-quark pair in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector, JHEP 06 (2022) 097 [arXiv:2111.06712 [INSPIRE]].

[108] CMS collaboration, Measurement of the $t\bar{t}H$ and $tH$ production rates in the $H \to b\bar{b}$ decay channel with $138 \text{ fb}^{-1}$ of proton-proton collision data at $\sqrt{s} = 13$ TeV, CMS-PAS-HIG-19-011, CERN, Geneva (2023).

[109] P. Bechtle et al., Probing the Standard Model with Higgs signal rates from the Tevatron, the LHC and a future ILC, JHEP 11 (2014) 039 [arXiv:1403.1582 [INSPIRE]].

[110] CMS collaboration, Search for a heavy Higgs boson decaying into two lighter Higgs bosons in the $\tau\tau b\bar{b}$ final state at 13 TeV, JHEP 11 (2021) 057 [arXiv:2106.10361 [INSPIRE]].

[111] ATLAS collaboration, Search for resonant and non-resonant Higgs boson pair production in the $b\bar{b}\tau^+\tau^-$ decay channel using 13 TeV pp collision data from the ATLAS detector, JHEP 07 (2023) 040 [arXiv:2209.10910 [INSPIRE]].

[112] P. Athron et al., Cosmological phase transitions: from perturbative particle physics to gravitational waves, arXiv:2305.02357 [INSPIRE].

[113] P. Athron, L. Morris and Z. Xu, How robust are gravitational wave predictions from cosmological phase transitions?, arXiv:2309.05474 [INSPIRE].

[114] AEDGE collaboration, AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space, EPJ Quant. Technol. 7 (2020) 6 [arXiv:1908.00802 [INSPIRE]].

[115] S. Kawamura et al., The japanese space gravitational wave antenna: DECIGO, Class. Quant. Grav. 28 (2011) 094011 [INSPIRE].

[116] S. Kawamura et al., Current status of space gravitational wave antenna DECIGO and B-DECIGO, PTEP 2021 (2021) 05A105 [arXiv:2006.13545 [INSPIRE]].

[117] V. Corbin and N.J. Cornish, Detecting the cosmic gravitational wave background with the big bang observer, Class. Quant. Grav. 23 (2006) 2435 [gr-qc/0512039] [INSPIRE].

[118] E. Di Valentino, S. Gariazzo and O. Mena, Most constraining cosmological neutrino mass bounds, Phys. Rev. D 104 (2021) 083504 [arXiv:2106.15267] [INSPIRE].

[119] L. Fromme and S.J. Huber, Top transport in electroweak baryogenesis, JHEP 03 (2007) 049 [hep-ph/0604159] [INSPIRE].

[120] S. Hessenberger and W. Hollik, Two-loop corrections to the $\rho$ parameter in Two-Higgs-Doublet Models, Eur. Phys. J. C 77 (2017) 178 [arXiv:1607.04610] [INSPIRE].

[121] S.Hessenberger, Two-loop corrections to electroweak precision observables in two-higgs-doublet-models, Ph.D. thesis, Tech. U., Munich (2018).

[122] S. Hessenberger and W. Hollik, Two-loop improved predictions for $M_W$ and $\sin^2\theta_{\text{eff}}$ in Two-Higgs-Doublet models, Eur. Phys. J. C 82 (2022) 970 [arXiv:2207.03845] [INSPIRE].

[123] F. Mahmoudi, SuperIso: a program for calculating the isospin asymmetry of $B \to K^*\gamma\gamma$ in the MSSM, Comput. Phys. Commun. 178 (2008) 745 [arXiv:0710.2067] [INSPIRE].

[124] F. Mahmoudi, SuperIso v2.3: a program for calculating flavor physics observables in Supersymmetry, Comput. Phys. Commun. 180 (2009) 1579 [arXiv:0808.3144] [INSPIRE].
[125] G. Bhattacharyya and D. Das, Scalar sector of two-Higgs-doublet models: a minireview, *Pramana* **87** (2016) 40 [arXiv:1507.06424] [INSPIRE].

[126] A. Barroso, P.M. Ferreira, I.P. Ivanov and R. Santos, *Metastability bounds on the two Higgs doublet model*, JHEP **06** (2013) 045 [arXiv:1303.5098] [INSPIRE].

[127] V. Cacchio, D. Chowdhury, O. Eberhardt and C.W. Murphy, *Next-to-leading order unitarity fits in Two-Higgs-Doublet models with soft $Z_2$ breaking*, JHEP **11** (2016) 026 [arXiv:1609.01290] [INSPIRE].

[128] W. Grimus, L. Lavoura, O.M. Ogreid and P. Osland, *A precision constraint on multi-Higgs-doublet models*, J. Phys. G **35** (2008) 075001 [arXiv:0711.4022] [INSPIRE].

[129] J. Haller et al., *Update of the global electroweak fit and constraints on two-Higgs-doublet models*, Eur. Phys. J. C **78** (2018) 675 [arXiv:1803.01853] [INSPIRE].