First shot of the smoking gun: probing the electroweak phase transition in the 2HDM with novel searches for A → ZH in ℓ + ℓ − t ¯ t and ννb ¯ b final states

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 ℓ + ℓ − t ¯ 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 ννb ¯ 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 ℓ + ℓ − t ¯ t final state shows a local 2 . 85 σ 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 ¯ 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.


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][9][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 ± .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][12][13][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. 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][20][21][22][23][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 Sect. 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 ≪ m A ≈ m H ± .As a consequence of this mass splitting between H and A the decay A → ZH is kinematically open [20], which (in contrast to the decay A → 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 → A → 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 → 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 ± or A and H ± [26]. 2 Accordingly, a possible detection of a A → ZH signal would provide information 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 β (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 2  12 , 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 BR(B → X s γ) = 3.49 ± 0.19 from the HFLAV collaboration [31], one finds a limit of m H ± ≳ 500 GeV [32,33]. 3The 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 → 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 bottomquark, tau-lepton or W -boson pairs were considered [35][36][37][38].Notably, until recently no searches for the A → ZH decay existed assuming the decay of H into top-quark pairs, which are favored for low tan β.Based on the combination of LHC searches and other experimental constraints, however, in the type II and type IV 2HDM (see Sect. 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 → t 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 → ZH → ℓ + ℓ − t 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. 4ecently, the experimental situation has changed with the first public results of searches for the signature A → ZH → ℓ + ℓ − t 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 ννb b final state.CMS has not yet released a result in the ℓ + ℓ − t 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 → ZH → ℓ + ℓ − t 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 σ (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 ± , then decaying as H ± → W ± H → ℓ ± νt 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 → 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 Sect. 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 Sect. 3 we present the numerical discussion of the impact of the new ATLAS A → ZH searches on the 2HDM parameter space.We divide our discussion of the new constraints into two parts, focusing on a low-tan β and a high-tan β regime in Sect.3.2.1 and Sect.3.2.2,respectively.Then, Sect.3.3 is devoted to an evaluation of the future prospects of the ℓ + ℓ − t t searches.The excess observed by ATLAS in the ℓ + ℓ − t t final state as a hint of a FOEWPT in the 2HDM is analysed in Sect.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 Sect. 4. App.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.App.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.

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.

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 Φ 1 → Φ 1 , Φ 2 → −Φ 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 2  12 .Given these constraints, the most general form of the Higgs potential is where the two doublets are conveniently parameterized as . ( u-type d-type leptons 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. Here v i is the vacuum expectation value (vev) acquired by the respective doublet.The eight degrees of freedom ϕ + i , ρ i and η i mix to give rise to the massive scalars h, H (CP-even), A (CP-odd), H ± and the three Goldstone bosons G 0 , G ± .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 α and β for the CP-even and the CP-odd/charged sector, respectively.
The coefficients m 2 11 , m 2 22 , m 2 12 , λ i (i = 1, ...5) in Eq. ( 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 β = v 2 /v 1 , and cos(β − α).In the limit cos(β − α) → 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, where two of these parameters are already fixed by experiment, namely m h ≈ 125 GeV and v = v 2 1 + v 2 2 ≈ 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 Tab. 1.It should be noted that due to the opposite Z 2 parities of the Higgs doublets each fermion can only be coupled to either Φ 1 or Φ 2 .In Tab. 1 we also give the tan β dependence of the coupling modifyers of the CP-odd Higgs boson A of the 2HDM to up-type and down-type quarks, ξ u,d A (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 modifyers ξ u,d A .Since ξ u A = cot β for all types, the dominant production mode for values of tan β ≈ 1 is gluon-fusion production involving a top-quark loop.For values of tan β ≳ 10 the production in association with bottom-quark pairs becomes important in type II and type IV, in which ξ d A = tan β, whereas in type I and type III this production mode is always substantially smaller than gluon-fusion production.

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 The temperature-independent part of the potential comprises the first three terms, where V tree is defined in Eq. ( 1), V CW is the one-loop Coleman Weinberg potential [48] incorporating the radiative corrections, and V 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 T , incorporates the one-loop thermal corrections in terms of the well-known J-functions (see e.g.Ref. [49]).The second term, V 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][52][53][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.

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. ( 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. 5We identify the 2HDM parameter space regions which, as the universe cools down, feature an EWbreaking 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]) 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 (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]).

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.
where v n is the EW vev at the nucleation temperature T n .The above condition (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 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).

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][64][65][66], as well as the turbulence generated after the collisions [67][68][69][70][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 . 6(ii) the phase transition strength α, 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, β/H.(iv) the bubble wall velocity in the rest frame of the fluid (and far from the bubble), v w .To compute α we follow Refs.[17,18], where ∆V (T * ) is the potential difference between the two vacua evaluated at the temperature of the phase transition, and ρ R is the radiation energy density of the universe.The inverse duration of the phase transition β/H can be generally calculated as where S(T ) is (as in Eq. ( 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 ≈ c s (c s being the speed of sound of the plasma)7 or v w → 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 from sound waves h 2 Ω sw and turbulence h 2 Ω turb , where sound waves are the dominant GW source for the FOEW-PTs 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 where T is the duration of the experiment, h 2 Ω Sens is the nominal sensitivity of the detector, computed according to the mission requirements [76], and h 2 Ω GW = h 2 Ω sw +h 2 Ω 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 T = 7 years, and consider a GW signal to be detectable if SNR > 1 (more stringent SNR detection criteria could also be considered [17]).

Numerical analysis
In this section we analyze in detail the impact of the recent ATLAS searches for the A → 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][79][80][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.

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]: 13 TeV including 139 fb −1 from Ref. [38] -gg → A → ZH → ℓ + ℓ − t t at 13 TeV including 140 fb −1 from Ref. [41] -gg → A → ZH → ννb b at 13 TeV including 140 fb −1 from Ref. [41] -b b → A → ZH → ννb b at 13 TeV including 140 fb −1 from Ref. [41] We note that the results from Ref. [38] update the previous ATLAS results based on 36.1 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 → 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][84][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.

New constraints on the 2HDM parameter space
In our analysis we target two separate 2HDM parameter regions: a low-tan β regime, corresponding to values of tan β ≤ 3, and a high-tan β regime with values of tan β ≥ 10.In the low tan β-regime, A can be produced with sizable cross sections via gluon-fusion, and we will show that searches in the ℓ + ℓ − t t final state exclude significant parts of previously unconstrained parameter space, whereas the ATLAS searches utilizing νν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 β regime, the state A can be produced with sizable rates in b b-associated processes for the 2HDM of Yukawa type II and IV.In this regime, the ννb b searches cover regions of parameter space that have already been probed via the ℓ + ℓ − b b final state.
For intermediate values of tan β, in between the low-and high-tan β 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 β for intermediate tan β values, such that the searches in the ℓ + ℓ − t t final states quickly loose sensitivity as tan β increases.At the same time the b b-associated production of A is enhanced in type II and type IV by factors of tan 2 β, such that the searches relying on this production mechanism become more sensitive for larger values of tan β.Yet, in our numerical analysis we find that one needs roughly an enhancement corresponding to tan 2 β ≈ 100 to achieve cross sections in the b b-associated production mode that are comparable to the exclusion limits resulting from the gg → A → ZH searches.Consequently, we chose the value tan β = 10 as a representative benchmark scenario in order to assess the sensitivity of the new searches in the high-tan β regime.The impact of the searches for even larger values of tan β can be extrapolated from the discussion of this scenario, as will be shown in detail in Sect.3.2.2.
As already discussed in Sect. 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 ± are assumed to be mass-degenerate, i.e. m A = m H ± , and the squared mass scale M 2 = m 2 12 /(s β c β ) is set equal to the mass of the heavy CP-even scalar H, i.e.M = m H . 8 In addition, the alignment limit cos(α − β) = 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 (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 ± as well as tan β.In the following, we discuss our results in the (m H , m A )-plane for different values of tan β within the two regimes discussed above.

Low tan β-region
In this section we present our results for the low-tan β regime, which focuses on the range 1 ≤ tan β ≤ 3. The lower bound on tan β was chosen because values of tan β 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 β depending on the 2HDM Yukawa type and the mass of the charged Higgs boson.As discussed above, we do not carry out a detailed investigation of the indirect limits from flavor physics in the following.
In Fig. 1 we show the impact of the new A → ZH searches from ATLAS [38,41] in the (m H , m A )-plane for tan β = 1 (upper left), tan β = 1.5 (upper right), tan β = 2 (lower left) and tan β = 3 (lower right).The upper left plot is valid independently of the chosen 2HDM Yukawa type.However, for tan β ̸ = 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 → A → ZH in the ℓ + ℓ − t t and the νν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 [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 Sect.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 Sect.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 → A → ZH → ℓ + ℓ − t 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 → H → t 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 Fig. 1 the parameter regions that exhibit a strong FOEWPT as defined in Sect.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 β = 1, all types
The results for tan β = 1 are shown in the upper-left plot of Fig. 1.One can see that the new A → ZH → ℓ + ℓ − t t ATLAS search (red) excludes the region 350 GeV ≲ m H ≲ 450 GeV and 650 GeV ≲ m A = m H ± ≲ 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 ± → tb decay [88,89], searches for neutral scalars decaying into top-quark pairs [90], and searches for the A → ZH decay in the ℓ + ℓ − b b final state [38], the mass range 300 GeV ≲ m H ≲ 450 GeV and 450 GeV ≲ m A = m H ± ≲ 700 GeV is now excluded.Fig. 1 also highlights that for tan β = 1 the parameter region with a strong FOEWPT to which the new ATLAS search is sensitive, assuming m A = m H ± , is already excluded by the charged Higgs-boson searches.Yet, we stress that if the condition m A = m H ± were relaxed by allowing for an additional mass gap between these states, i.e. m H ± > 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 β = 1.5, type II The results for tan β = 1.5 are shown in the upper right plot of Fig. 1.While for low tan β 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 β = 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 β = 1 (upper left plot), one can see that the region excluded by the searches for the charged scalars via pp → H ± tb → tb tb, where the cross section times branching ratio roughly scales with 1/tan 2 β in the low-tan β regime, is substantially smaller.This search loses even more sensitivity where the decay H ± → W ± H is kinematically allowed, giving rise to the slope of the corresponding excluded region for m H ≲ 500 GeV (which is more pronounced than for tan β = 1 because of the reduced H ± tb coupling).As a consequence, for tan β = 1.5 the H ± → 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 ≈ 2 m t by searches for H → τ + τ − [91,92] and by searches for the A → ZH

ATLAS obs.: gg
HiggsBounds comb.95% C.L. excl.decay using the ℓ + ℓ − b b final state [38].Above the di-top threshold, the decay H → t t very quickly dominates, and the new ATLAS search in the ℓ + ℓ − t t final state is the most sensitive one.In contrast to the tan β = 1 case, for tan β = 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 ≳ 400 GeV and m A = m H ± ≳ 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σ excess for m A ≈ 650 GeV and m H ≈ 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 Fig. 1, and the corresponding cross section are such that they fall into the strong FOEWPT region.In Sect.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 ± → W ± H → ℓ ± νt t, which so far has not been performed. 9In particular, producing the charged scalar via pp → tbH ± 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 β = 1.5 the new smoking gun search using the ννb 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 ≳ 200 GeV.However, in contrast to the searches using the ℓ + ℓ − t t final state indicated by the red shaded region, the blue shaded region indicating the new exclusion region resulting from the search using the ννb 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 → ZH with Z → ℓ + ℓ − and the decay of H into bottom-quark pairs [38].One should note, however, that the new A → ZH search in the ννb 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 ℓ + ℓ − b b and ℓ + ℓ − 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 Fig. 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 ννb 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 ννb b final state can potentially provide new constraints.
• tan β = 2, type IV We show the results for tan β = 2 in the lower left plot of Fig. 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 → b b is enhanced in type IV for tan β > 1, such that the searches in the ννb b final state become more important with increasing values of tan β.Moreover, in type IV the decay width for H → τ + τ − is suppressed approximately by 1/tan 2 β, whereas it is enhanced by about a factor of tan 2 β in type II.Hence, while in type II the parameter region below the di-top threshold, i.e. m H < 2m t , is entirely excluded by the searches for di-tau resonances, in type IV the ννb b search can potentially yield stronger constraints.One can see in the lower left plot of Fig. 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 ≲ tan β ≲ 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 → H ± tb → tb tb, suppressed by factors of about 1/tan 2 β in the low-tan β regime, cannot probe the parameter plane in this case.Below the di-top threshold, we find that the A → ZH searches in the ℓ + ℓ − b b (gray) and the ννb 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 → ZH using the decay Z → ℓ + ℓ − have been performed by ATLAS [38], which are more powerful than the new searches using the Z → νν decay (the corresponding CMS search using the Z → ℓ + ℓ − 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 → νν (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 → b 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 → t t is in fact the only channel that can probe the parameter plane.As indicated with the red shaded area, the searches in the ℓ + ℓ − t t final state are able to exclude masses smaller than m H ≈ 400 GeV and m A ≈ 750 GeV for the lighter and the heavier BSM resonance, respectively.As it is also visible in the plots for tan β = 1, tan β = 1.5 and tan β = 2 of Fig. 1, the difference between the expected (red dashed line) and the observed (red solid line) exclusion region resulting from the searches using the ℓ + ℓ − t 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 β = 1.5 and tan β = 2, where the observed limit is stronger than the expected one).
• tan β = 3, type IV As a final step of the discussion of the low-tan β regime we consider a value of tan β = 3.The results of our analysis are shown in the lower right plot of Fig. 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 ℓ + ℓ − t 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 Sect.3.3 the prospects for probing the benchmark plane for tan β = 3 in future runs of the LHC, in which roughly 20 times more integrated luminosity will be collected by both ATLAS and CMS. 10 At and below the di-top threshold m H ≈ 2m t the results are similar to the case of tan β = 2, where the smoking gun searches relying on the decay H → b 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 → t t are suppressed approximately by 1/tan 2 β, and the partial width for the decay H → b b is conversely enhanced by (approximately) tan 2 β.As a result, the gray exclusion region from the searches in the A → ZH → ℓ + ℓ − b b channel extends to slightly larger masses for tan β = 3 compared to tan β = 2 (lower left plot).

High tan β-region
In the discussion above, we investigated the low-tan β 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 β ≳ 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 b-associated production, which is enhanced by about tan 2 β in these types.Consequently, focusing on the high-tan β regime here, we expect the searches for the signature A → ZH assuming b 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 b-associated production utilizing the ννb b final state, whereas the smoking-gun search utilizing the ℓ + ℓ − t t final state was considered only assuming gluon-fusion production of the heavy BSM resonance.Thus, the only relevant searches for the A → ZH decay in the following discussion will be the previously reported searches utilizing the ℓ + ℓ − b b final state [35,38] and the new searches utilizing the ννb b final state [41].
In Fig. 2 we show our results for tan β = 10 as a representative benchmark scenario for the hightan β regime.The color coding of the exclusion regions and the scatter points is the same as in Fig. 1, except for the yellow dashed and solid lines indicating the expected and observed exclusion limits resulting from the recent ATLAS search for b b → A → ZH → ννb 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 b → A → ZH → ℓ + ℓ − b 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 β ≳ 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 ννb 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 ℓ + ℓ − b 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 Fig. 2), the new searches using the ννb b final state could potentially give rise to new constraints.
The two LHC searches relevant in Fig. 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 β.We can therefore extrapolate based on the results for tan β = 10 shown in Fig. 2 that also for larger values of tan β the searches making use of the Z → ℓ + ℓ − decay mode are more promising to probe the considered benchmark plane compared to the searches using the Z → νν decay mode.It should also be taken into account that for larger values of tan β other LHC searches become relevant in type IV. 11In particular, searches for new resonances produced in b 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 β ≳ 15.Moreover, for such values of tan β 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.

Future prospects for ℓ + ℓ − t t searches
In Sect.3.2.1 we have demonstrated that the new ATLAS smoking-gun searches targeting the ℓ + ℓ − t t final state exclude sizable parts of previously allowed parameter space of the 2HDM assuming values of tan β 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 ≲ tan β ≲ 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 → ZH decay in the ℓ + ℓ − t t final state during future runs of the LHC and the High-Luminosity LHC (HL-LHC).As input 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 App.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.
Here, σ exp.95% CL Run 2 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 σ exp.95% CL proj.
is the future projection of the expected cross-section limits depending additionally on the assumed integrated luminosity 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 Fig. 3 we show our projections in the 2HDM benchmark plane introduced in Sect.3.2 for the Yukawa type II with tan β = 1.5 in the left plot and tan β = 3 in 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 Fig. 1, and the red dashed lines indicate the expected exclusion regions for different values of the integrated luminosity, ranging from L = 300 fb −1 (end of LHC Run 3) to L = 3000 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 L = 140 fb −1 , and the magenta star indicates the masses for which ATLAS has observed the most pronounced local excess (see Sect. 3.4).As already discussed in Sect.3.2.1,currently the smoking-gun searches are not able to probe the benchmark plane for tan β = 3 (see the lower left plot of Fig. 1).Accordingly, no red shaded region is visible in the right plot of Fig. 3.
One can observe in the left plot of Fig. 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 β = 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 Sect.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 ≲ 550 GeV and m A ≲ 700 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 β 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 Sect.2.5 that for the considered values of tan β 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 β, where the gg → A → ZH → ℓ + ℓ − t t searches lose their sensitivity.A GW signal potentially detectable with LISA therefore cannot be fully probed with the searches in the ℓ + ℓ − t 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 ± 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][103][104] for related investigations in other extended scalar sectors).
In the right plot of Fig. 3, in which we show the projections for tan β = 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 ≈ 550 GeV and m A ≈ 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 → A → ZH → ℓ + ℓ − b b (both for type II and type IV), as was discussed in detail in Sect.3.2.1 (see the lower right plot of Fig. 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 → b b is strongly suppressed for m H > 2m t .
In summary, the fact that the smoking-gun search for gg → A → ZH → ℓ + ℓ − t t will be able to probe masses of m H > 2m t in the low tan β 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 1 st -order EW phase transition in the 2HDM?
We now turn to the analysis of the local 2.85 σ 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. 12To this end, we first determine the cross section that would be associated with the excess.Since 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 σ(gg corresponding to the excess observed by ATLAS at the above-mentioned mass values for the two types of Higgs bosons. 13n 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 ± = 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 ± 50 GeV which we use as a rough estimate for the potential signal region.

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 Sect.3.2 that values of tan β ≈ 1.5 are required for obtaining sufficiently large cross sections to describe the excess (see the upper right plot of Fig. 1).We therefore vary tan β in the interval 1 ≤ tan β ≤ 2, and in addition we consider deviations from the alignment limit in terms of the free parameter cos(β − α), see the discussion in Sect.2.1.Since we only consider tan β values close to tan β = 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(β − α), 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 Fig. 4 we show the parameter plane with cos(β − α) on the horizontal axis and tan β on the vertical axis and the additional 2HDM parameters set according to the discussion above, 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 χ 2 -fit to the various measurements.Specifically, we demand that the χ 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 σ (∆χ 2 125 < 6.18) away from the SM result (χ 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(β − α) = 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 β ≲ 1.2 (both plots) are excluded by the cross section limits from searches for charged Higgs bosons using the H ± → tb decay [88,89].The gray regions at |cos(β − α)|≳ 0.025 (both plots) are excluded by searches for resonant Higgs-boson pair production using the H → hh decay in the b bτ + τ − 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(β − α) = 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. σ(gg → A) × BR(A → ZH) × BR(H → t t).In the green shaded areas the predicted values are within the interval given in Eq. ( 11), corresponding to a description of the excess within 1 σ.
In both plots one can see that for values of 1.2 ≲ tan β ≲ 1.8 and values of cos(β − α) 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 → V V and H → hh, become relevant outside of the alignment limit.As a result, for a fixed value of tan β the predicted cross sections decrease with increasing values of |cos(β − α)|.This gives rise to the the shape of the green lines peaking at cos(β − α) ≈ 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 ≲ cos(β − α) ≲ 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(β − α)|≳ 0.025 are excluded by the cross section measurements of the SM-like Higgs boson and/or by the searches for H → hh → b bτ + τ − .In type II the exclusion regions of both these constraints are largely overlapping, as shown in the right plot of Fig. 4, where the condition ∆χ 2  125 < 6.18 yields an exclusion for values of cos(β − α) ≲ −0.02 and cos(β − α) ≳ 0.045, while the searches for H → hh exclude |cos(β − α)|≳ 0.025.The situation is different in type I, where for negative values of cos(β − α) the parameter region excluded by HiggsSignals features values of cos(β − α) ≲ −0.15, whereas the H → hh searches exclude already values of cos(β − α) ≲ −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 σ 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.

Gravitational wave detection
In the discussion in Sect.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 ± 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 Fig. 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 Sect.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 SNR ≳ 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 The mass values for which the excess observed by ATLAS in the ℓ + ℓ − t t final state was most pronounced (2.85 σ 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.
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 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 Fig. 6 the spectral shape of the GW backgrounds produced during a FOEWPT for several parameter points with masses of the heavy scalars specified in Tab. 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 Fig. 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 Tab. 2 the parameters for the point (m H , m A ) = (450, 650) GeV although we omit its GW spectrum in Fig. 6 because of the smallness of the SNR.The spectral shapes of the GW backgrounds are computed as discussed in Sect.2.5, where the solid curves depict the sound-wave contribution h 2 Ω sw only, whereas the dashed curves depict the sum of sound-wave and turbulence contributions, i.e. h 2 Ω sw + h 2 Ω 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 spacebased GW detectors.One can see that only for the smallest value of m H = 417.2GeV, 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 modeldependent) 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. 14

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 β is of interest, which is preferred by EW baryogenesis.Since the searches in the ℓ + ℓ − t 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 cross-section limits reported by ATLAS.We focused on 2HDM benchmark scenarios assuming the alignment limit, cos(β −α) = 0, mass degeneracy between the charged Higgs boson and the pseudoscalar state, m H ± = m A , and setting the decoupling scale M (defined by the relation m 2 12 = M 2 sin β cos β) 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 β regime with 1 ≤ tan β ≤ 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 ℓ + ℓ − t t has high sensitivity.We found that for tan β = 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 ± → tb, H → τ + τ − and H, A → t t, masses of 300 GeV ≲ m H ≲ 450 GeV are now entirely excluded in this scenario, while previously for tan β = 1 a wide parameter space region with 350 GeV ≲ m H ≲ 460 GeV and 650 GeV ≲ m A ≲ 800 GeV was left unconstrained.For increasing values of tan β, the charged scalar searches lose sensitivity, and in particular for tan β = 1.5 the ℓ + ℓ − t 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 > 2 m t .For tan β = 2, and irrespectively of whether the presence of a FOEWPT is required or not, the search for the smoking gun signature A → ZH with H → t 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 β = 3, the largest value considered in the low-tan β regime, the ℓ + ℓ − t 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 → ZH with H → b 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 ≈ 400 GeV.We have furthermore pointed out a novel, not yet performed LHC search that would be complementary to the smoking gun A → ZH search in the low tan β region as a probe of the 2HDM parameter region featuring a strong FOEWPT.The channel that we propose for experimental analysis consists of H ± production (e.g. via pp → H ± tb) followed by the decay H ± → W ± H → ℓ ± ν t t.
In addition to the searches in the ℓ + ℓ − t t final state, ATLAS also reported for the first time searches for the A → ZH decay making use of the decay of the Z boson into neutrino pairs and the H → b b decay mode.Here, both the gluon-fusion and the b b-associated production of A were considered.These searches in the ννb b final state may become important for large values of tan β.Investigating a representative 2HDM benchmark scenario with tan β = 10 (and the remaining 2HDM parameters as described above), we found that in the 2HDM the earlier searches in the ℓ + ℓ − b b final state give rise to stronger exclusions than the new ννb b ATLAS search.Nevertheless, it should be taken into account that the new searches utilizing the Z → νν 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 ννb 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 ℓ + ℓ − t t searches described above, we investigated in Sect.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 > 2 m t and tan β ≥ 3 will become accessible.We therefore anticipate that the smoking gun signature A → ZH with H → t 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 β 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 ℓ + ℓ − t 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 σ.We have demonstrated that the excess can be described at the 1 σ level in the 2HDM type I and II (as representative scenarios) in the approximate range of tan β values between 1.2 and 1.8 for the alignment limit (a 1 σ description of the excess is also possible for small departures from the alignment limit, with the viable tan β 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 ± → W ± H followed by H → t 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. 15Therefore 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 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 App. A.
Our analysis shows that the smoking gun signature A → ZH (with H → t 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.
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. ( 3)).During the installation of this package the following external codes will also be downloaded and installed: AnyHdecay [83][84][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][79][80][81]) and with existing measurements of the detected Higgs boson at about 125 GeV (HiggsSignals [86,87,109]).
THDM EWPOS [120][121][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.
Note that HiggsTools will only be installed if not already installed in the current python environment of the user.
A. 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: 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π to the eigenvalues of the scattering matrix.One can change the upper limit by giving an argument in the function ('Unitarity to 1:', pt.check_perturbative_unitarity(1.)) where the value of the argument multiplied by 16π is the applied upper limit (thus 0.5 is the default value).One can also check for NLO perturbative unitarity in the high enery limit according to the expressions derived in Ref. [127]: ('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 → sγ and B s → µ + µ − are checked.One can access the specified allowed region by typing 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 χ 2 fit to the experimental EWPO in terms of the oblique parameters S, T , and U by typing 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 σ confidence level.
Regarding the collider constraints, for compatibility with the 125 GeV Higgs-boson measurements one can specify the allowed range of ∆χ 2  125 to perform a test using the best χ 2 fit performed with HiggsSignals.The default is to allow parameter points with χ 2  125 -values that are not more than 2 σ away from the SM χ 2 125 -value, which corresponds to ∆χ 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 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.

print('sigma(gg -> h)', pt.XS['Hl']['gg'])
where Hl 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 b-associated production, vector-boson fusion, t tassociated 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 where the predictions for the branching ratios are computed via the interface to the AnyHdecay library.

B Comparison to previous CMS projections
In Ref. [24], we estimated the projected (HL-)LHC sensitivity for the process A → ZH in the Z t t final state for several integrated luminosities.We used the results for the expected sensitivity in this channel obtained in a Master thesis for the CMS Collaboration [42,43] and applied them to the (m H , m A ) parameter plane also investigated in this paper, with tan β = 3 and m H ± = m A .Here we aim to compare these prior sensitivity projections with those based on the ATLAS expected limits, as detailed in Section 3.3.In Fig. 7 we present the resulting expected 95% C.L. exclusion limits corresponding to integrated luminosities of 300, 600, 1000, 3000 fb −1 , projected for future (HL-) LHC runs.On the left-hand side, we display the exclusion regions derived from a straightforward rescaling of the CMS expected limits for different luminosity values, thus not accounting for changes in systematic uncertainties.On the right-hand side we show the exclusion regions resulting from a similar rescaling process, albeit based on ATLAS expected limits, again without accounting for possible changes in systematic uncertainties.The color code shows the phase transition strength ξ n for parameter points featuring a FOEWPT with ξ n > 1.The blue region indicates the area that features electroweak symmetry non-restoration at high temperatures (see Ref. [24] for details).The comparison demonstrates good agreement between both sets of projections, reinforcing the robustness of our conservative estimate of the future prospects.

2 12 = m 2 H
sin β cos β m A = m Z + m H p e r t .u n it a r it y v a c u u m s t a b .ATLAS obs.: gg → A → ZH → ℓ + ℓ − tt ATLAS exp.: gg → A → ZH → ℓ + ℓ − tt HiggsBounds comb.95 %C.L. excl.

2 12 = m 2 H
sin β cos β m A = m Z + m H p e r t .u n it a r it y v a c u u m s t a b .

2 12 = m 2 H 1 Figure 1 :
Figure 1: Impact of the new ATLAS searches for the A → ZH signature in the (m H , m A )-plane for tan β = 1 (upper left), tan β = 1.5 and type II (upper right), tan β = 2 and type IV (lower left), and tan β = 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 ℓ + ℓ − t t and ννb 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 σ local significance) observed by ATLAS in the ℓ + ℓ − t t search are indicated with a magenta star in the upper right plot.

2 12 = m 2 Hv n /T n 1 Figure 2 :
tan β = 10, cos(β − α) = 0, mA = m H ± , m sin β cos β m A = m Z + m H p e r t .u n i t a r i t y v a c u u m s t a b .ATLAS obs.: b b → A → ZH → ννb b ATLAS exp.: b b → A → ZH → ννb b HiggsBounds comb.95% C.L. excl.As in Fig. 1, but for tan β = 10 in type IV.Parameter space regions excluded by the new ννb b searches in the b b-associated production channel are indicated in yellow, while the yellow dashed line indicates the expected exclusion limit.

18 σ 2 12 = m 2 H 1 − 18 σ 2 12 = m 2 H 1 Figure 4 :
Figure4: 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σ 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 σ 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 → hh decay in b bτ + τ − final states.

1 SNRFigure 5 :
Figure5: The signal-to-noise ratio (SNR) of the scanned parameter points in the (m H , m A )-plane for tan β = 1.5.The mass values for which the excess observed by ATLAS in the ℓ + ℓ − t t final state was most pronounced (2.85 σ 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.

Figure 6 :
Figure6: Gravitational wave spectra for parameter points specified in Tab. 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.

( 1 Figure 7 :
Figure 7: Projected exclusion regions in the (m H m A plane with tan β = 3 and m ± = m A and for integrated luminosities of 300, 600, 1000, 3000 fb −1 , expected to be collected in future runs of the LHC.The displayed limits are derived from rescaled CMS (left) and ATLAS (right) expected limits for the ℓ + ℓ − t t final state.The color bar indicates the strength of the phase transition.The blue points indicate the parameter region that features electroweak symmetry non-restoration at high temperatures (see Ref. [24] for more details).

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 Fig.6.Dimensionful parameters are given in GeV.The SNR values evaluated for LISA include the turbulence contribution.
Installation requires a python3 environment and compilers for Fortran and C++.The package and all its dependencies are installed by executing:To import the package in an example notebook type from thdmTools.parampointimport Parampoint To define the parameters of the point as a dictionary in python 1 Installation thdmTools is publicly available at: https://gitlab.com/thdmtools