Confronting the Broken Phase of the N2HDM with Higgs Data

Abstract

The broken phase of the next-to-two-Higgs-doublet model (N2HDM) constitutes an archetype of extended Higgs sectors. In the presence of a softly broken ${\mathrm{Z}}_2$ symmetry throughout the scalar and Yukawa sectors, as the additional gauge singlet field does not interact with fermions, the model admits four variants of Yukawa interactions between the doublets and Standard Model fermions. We confront each type with experimental Higgs data, especially those from CMS and ATLAS detectors at the LHC. Interfacing the models with the state-of-the-art package $\mathtt{HiggsTools}$, we perform a statistical ${\chi }^2$ analysis to determine the best-fit points and exclusion limits at the $95\%$ and $68\%$ C.L.’s and identify SM-like Higgs measurements that affect each type the most. We further analyze the exclusion bounds on the additional Higgs bosons at the $95\%$ C.L., paying special attention to searches for hypothetical non-SM Higgs resonances decaying into a pair of bosons or fermions. We show regions where the additional Higgs bosons do not satisfy the narrow-width approximation utilized in most experimental searches.

1. Introduction

The next-to-two-Higgs-doublet model (N2HDM) is a well-motivated Beyond the Standard Model (BSM) extension, providing a plethora of new possibilities for direct and indirect experimental searches for new physics (see [1,2] and references therein). This is due to the generic structure of its Higgs sector, which contains an additional real singlet compared to the 2HDM. Generally, there are different aspects that motivate the study of doublet and singlet extensions of the Standard Model (SM) [3]. Among these are the continuing efforts to understand the nature of electroweak symmetry breaking (EWSB) and searches for deviations from the SM due to the possibility of the presence of additional Higgs doublets and/or singlets. Given that the Large Hadron Collider (LHC) is consistently examining the properties of the discovered scalar and searching for additional scalars predicted by such extensions [4,5], it is vital to understand to what extent those deviations are allowed and to confront BSM extensions with experimental Higgs data.

In fact, several aspects of the extended SM with singlets and doublets have been considered in the literature, including collider phenomenology, dark matter, and cosmology [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. The N2HDM can be seen as a baseline model that captures a range of phenomenological effects that arise from adding new doublets and singlets to the Higgs potential. It is an extension of the CP-conserving 2HDM with an additional real singlet (the CP-violating variant was investigated in [22]), where two ${\mathrm{Z}}_2$ symmetries are imposed to eliminate flavor-changing neutral currents (FCNCs) at the tree level, and possibly provides a stable dark matter (DM) candidate. Its vacuum structure is more involved than CP-conserving 2HDMs or singlet extensions, since a vacuum expectation value (VEV) can be acquired by two Higgs doublets and the singlet (the broken phase), only two Higgs doublets (the dark singlet phase), only one Higgs doublet and the singlet (the dark doublet phase), or only one Higgs doublet (the fully dark phase). In the first case, no DM candidate is present, while the other cases can provide a singlet DM, a doublet DM, or two DM candidates (see [23] for details and references to earlier work along this direction).

As the two doublets can interact with SM fermions in the presence of a softly broken ${\mathrm{Z}}_2$ symmetry preventing FCNCs, the Yukawa sector inherits the four different types associated with the 2HDM [24,25,26,27,28,29,30,31]: Type 1, Type 2, Type X (Lepton-Specific), and Type Y (Flipped). Other more generic scenarios for Yukawa interactions, without ${\mathrm{Z}}_2$ symmetry, have been considered in the 2HDM [32,33,34]. Moreover, in the N2HDM, both the CP-odd state A and the charged Higgs pair $H^{\pm }$ have the same structure as in the 2HDM. However, the constraints specific to the N2HDM may have some indirect effects on such states.

From the theoretical side, the renormalization of the full model was carried out in [35], where it was shown that the effects of corrections can be sizable. The impact of electroweak corrections was thoroughly analyzed in [36] and implemented in [37]. Cosmological aspects of the Type 2 model were investigated in [38], demonstrating cases where electroweak symmetry is not restored. A comprehensive analysis of vacuum instabilities was provided in [39]. Additionally, the naturalness of the model was considered in [40].

From the phenomenological and experimental sides, an analysis of Types 1 and 2 in an approximated version of the model was performed in [1] with $H_2$ considered SM-like, while a systematic analysis of the same types, focusing on wrong-sign regions and the singlet component of the SM-like Higgs boson (which could be any of the CP-even Higgs bosons), was conducted in [2]. The CMS collaboration searched for resonant pair production of Higgs bosons in the $b\overline{b}ZZ$ final state and provided an interpretation of the results specific to the N2HDM scenario [41]. A dedicated investigation of di-Higgs production in the four-photon final state was carried out in [42], assuming Type 1, while a comprehensive analysis of the limits on di-Higgs production was presented in [43] for Types 1 and 2, with careful consideration of the distinction between resonant and non-resonant regions. The accommodation of an additional 96 GeV Higgs boson was considered in [44,45,46]. The prospects of the model for the future electron–positron collider were presented in [47]. The model is implemented in the public tool $\mathtt{ScannerS}$ [48], which facilitates sophisticated phenomenological studies.

Given the continuing interest and interchanging efforts from both the phenomenology and experimental sides, especially with the advent of several new results from LHC SM-like Higgs measurements and searches for additional Higgs bosons, our aim in this paper is to expand on previous works and confront all types of N2HDMs with the latest Higgs data available in the public code $\mathtt{HiggsTools}$ (HT). Specifically, we interface the model with HT, carry out an up-to-date statistical ${\chi }^2$ analysis, and obtain exclusion limits on the SM-like Higgs boson with $95\%$ and $68\%$ Confidence Levels (C.L.’s), as well as $95\%$ C.L. exclusion limits on the additional Higgs bosons, paying special attention to the pair production of bosons through a heavy scalar resonance.

This paper is organized as follows. In Section 2, we describe the theoretical aspects of the Higgs sector. In Section 3, we present the analysis scheme and the constraints taken into account. In Section 4, we provide the results and discuss them in the context of recent LHC searches and measurements. Finally, we conclude in Section 5.

2. The Higgs Sector

In terms of the two $SU{\left(2\right)}_L$ Higgs doublets ${\Phi }_1$ and ${\Phi }_2$ and the real singlet field ${\Phi }_S$, the Higgs sector of the CP-conserving N2HDM is described by the following scalar potential:

$$\begin{array}{cc}\hfill V_{\mathrm{N}2\mathrm{HDM}}& =V_{2\mathrm{HDM}}+V_{\sin \mathrm{glet}},\hfill \end{array}$$

(1)

where

$$\begin{array}{cc}\hfill V_{2\mathrm{HDM}}& =m_{11}^2|{\Phi }_1{|}^2+m_{22}^2{\left|{\Phi }_2\right|}^2-m_{12}^2({\Phi }_1^{†}{\Phi }_2+\mathrm{h}.\mathrm{c}.)\hfill \\ \hfill & +{\displaystyle \frac{{\lambda }_1}{2}}{\left({\Phi }_1^{†}{\Phi }_1\right)}^2+{\displaystyle \frac{{\lambda }_2}{2}}{\left({\Phi }_2^{†}{\Phi }_2\right)}^2\hfill \\ \hfill & +{\lambda }_3\left({\Phi }_1^{†}{\Phi }_1\right)\left({\Phi }_2^{†}{\Phi }_2\right)+{\lambda }_4\left({\Phi }_1^{†}{\Phi }_2\right)\left({\Phi }_2^{†}{\Phi }_1\right)\hfill \\ \hfill & +{\displaystyle \frac{{\lambda }_5}{2}}[{\left({\Phi }_1^{†}{\Phi }_2\right)}^2+\mathrm{h}.\mathrm{c}.],\hfill \end{array}$$

(2)

and

$$\begin{array}{cc}\hfill V_{\mathrm{singlet}}& ={\displaystyle \frac{1}{2}}{m}_S^2{\Phi }_S^2+{\displaystyle \frac{{\lambda }_6}{8}}{\Phi }_S^4+{\displaystyle \frac{{\lambda }_7}{2}}\left({\Phi }_1^{†}{\Phi }_1\right){\Phi }_S^2+{\displaystyle \frac{{\lambda }_8}{2}}\left({\Phi }_2^{†}{\Phi }_2\right){\Phi }_S^2.\hfill \end{array}$$

(3)

All parameters are assumed to be real:

• $m_{11}^2$, $m_{22}^2$, and $m_S^2$: mass-squared parameters for ${\Phi }_1$, ${\Phi }_2$, and ${\Phi }_S$.
• $m_{12}^2$: soft-breaking mass-squared parameter.
• ${\lambda }_1$–${\lambda }_8$: quartic couplings.

The structure of the potential is dictated not only by SM symmetries but also by two additional discrete symmetries: ${\mathrm{Z}}_2$ (${\Phi }_{1,S}$: even, ${\Phi }_2$: odd) and ${\mathrm{Z}}_2^{\prime }$ (${\Phi }_{1,2}$: even, ${\Phi }_S$: odd). The first symmetry is softly broken by the $m_{12}^2$ term, similar to the 2HDM case, while the other one is spontaneously broken once the singlet field acquires a VEV. This structure forbids cubic terms in the potential (see [2] for more details).

The broken phase of the N2HDM is defined as the case where all fields obtain VEVs,

$$\begin{array}{c}\hfill 〈{\Phi }_1〉=\left(\begin{array}{c}0\\ {\displaystyle \frac{v_1}{\sqrt{2}}}\end{array}\right),〈{\Phi }_2〉=\left(\begin{array}{c}0\\ {\displaystyle \frac{v_2}{\sqrt{2}}}\end{array}\right),〈{\Phi }_S〉=v_s,\end{array}$$

(4)

where the factor ${\displaystyle \frac{1}{\sqrt{2}}}$ sets the convention $v=\sqrt{v_1+v_2}=246.22$ GeV for the electroweak VEV. At its minimum, the potential takes the form

$$\begin{array}{cc}\hfill V& ={\displaystyle \frac{m_{11}^2{v}_1^2}{2}}+{\displaystyle \frac{m_{22}^2{v}_2^2}{2}}-m_{12}^2{v}_1{v}_2+{\displaystyle \frac{{\lambda }_1{v}_1^4}{8}}+{\displaystyle \frac{{\lambda }_2{v}_2^4}{8}}+{\displaystyle \frac{{\lambda }_{345}{v}_1^2{v}_2^2}{4}}\hfill \\ \hfill & +{\displaystyle \frac{m_S^2{v}_S^2}{2}}+{\displaystyle \frac{{\lambda }_6{v}_S^4}{8}}+{\displaystyle \frac{{\lambda }_7{v}_1^2{v}_S^2}{4}}+{\displaystyle \frac{{\lambda }_8{v}_2^2{v}_S^2}{4}},\hfill \end{array}$$

(5)

which is minimized by

$$\left({\displaystyle \frac{\partial V}{\partial v_1}}{|}_{\min },{\displaystyle \frac{\partial V}{\partial v_2}}{|}_{\min },{\displaystyle \frac{\partial V}{\partial v_S}}{|}_{\min }\right)=0,$$

(6)

resulting in three conditions:

$$\begin{array}{cc}\hfill m_{11}^2& =m_{12}^2{\displaystyle \frac{v_2}{v_1}}-{\displaystyle \frac{1}{2}}(v_1^2{\lambda }_1+v_2^2{\lambda }_{345}+v_S^2{\lambda }_7)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill m_{22}^2& =m_{12}^2{\displaystyle \frac{v_1}{v_2}}-{\displaystyle \frac{1}{2}}(v_1^2{\lambda }_{345}+v_2^2{\lambda }_2+v_S^2{\lambda }_8)\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill m_S^2& =-{\displaystyle \frac{1}{2}}(v_1^2{\lambda }_7+v_2^2{\lambda }_8+v_S^2{\lambda }_6),\hfill \end{array}$$

(9)

where ${\lambda }_{345}={\lambda }_3+{\lambda }_4+{\lambda }_5$.

The fields in Equations (2) and (3) can be parametrized by expanding their neutral components around the VEVs:

$${\varphi }_i^0={\gamma }_i(v_i+{\rho }_i+i{\sigma }_i),$$

(10)

where $i=1,2,s$ with ${\gamma }_{1,2}={\displaystyle \frac{1}{\sqrt{2}}}$, ${\gamma }_s=1$, and ${\sigma }_s=0$. Focusing on the CP-even neutral components, we have

$$\begin{array}{cc}\hfill \mathcal{R}\left[{\varphi }_1^0\right]& ={\displaystyle \frac{v_1+{\rho }_1}{\sqrt{2}}}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \mathcal{R}\left[{\varphi }_2^0\right]& ={\displaystyle \frac{v_2+{\rho }_2}{\sqrt{2}}}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\varphi }_S& =v_S+{\rho }_S\hfill \end{array}$$

(13)

In the $\{{\rho }_1,{\rho }_2,{\rho }_s\}$ basis, the elements of the mass-squared matrix ${\mathcal{M}}_{\rho }^2$ can be derived by

$$M_{ij}^2={\left.{\displaystyle \frac{{\partial }^{2}V}{\partial {\rho }_i\partial {\rho }_j}}\right|}_{{\rho }_k=0}$$

(14)

Substituting $v_1=vcos\beta$ and $v_2=vsin\beta$, where $tan\beta ={\displaystyle \frac{v_2}{v_1}}$, and using the minimization conditions, the matrix elements of ${\mathcal{M}}_{\rho }^2$ are

$$\begin{array}{cc}\hfill M_{11}^2& ={\lambda }_1{v}^2{cos}^2\beta +m_{12}^{2}tan\beta \hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill M_{12}^2& ={\lambda }_{345}{v}^{2}cos\beta sin\beta -m_{12}^2\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill M_{13}^2& ={\lambda }_{7}v{v}_{S}cos\beta \hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill M_{22}^2& ={\lambda }_2{v}^2{sin}^2\beta +{\displaystyle \frac{m_{12}^2}{tan\beta }}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill M_{23}^2& ={\lambda }_{8}v{v}_{S}sin\beta \hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill M_{33}^2& ={\lambda }_6{v}_S^2\hfill \end{array}$$

(20)

One can exploit the properties of this self-adjoint matrix to set an upper bound on the lightest eigenvalue (${\Lambda }_{\min }$) of its diagonal form ${\mathcal{M}}_H^2$ (corresponding to mass eigenstates squared: $m_{H_i}^2$). In particular, rotating the top-left $2\times 2$ part of ${\mathcal{M}}_{\rho }^2$ by a unitary matrix defined in terms of the angle $\beta$, the upper bound is

$${\Lambda }_{min}\le min(A,B),$$

(21)

where A and B are the diagonal elements of the rotated $2\times 2$ submatrix:

$$\begin{array}{cc}\hfill A& =v^2\left({\lambda }_1{cos}^4\beta +2{\lambda }_{345}{cos}^2\beta {sin}^2\beta +{\lambda }_2{sin}^4\beta \right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill B& ={\displaystyle \frac{1}{2}}{m}_{12}^2\left(3-cos2\beta \right)cot\beta +{sin}^2\beta \left(v^2\left({\lambda }_1+{\lambda }_2-2{\lambda }_{345}\right){cos}^2\beta +m_{12}^{2}tan\beta \right)\hfill \end{array}$$

(23)

This sets an upper limit on the lightest CP-even Higgs boson in the model (in this paper, we consider this state as the SM-like Higgs boson). Formally, ${\mathcal{M}}_{\rho }^2$ can be diagonalized using its eigenvalues and eigenvectors, taking us from gauge eigenstates to mass eigenstates ($H_1,H_2,H_3$). Given the excessively long expressions, it is more convenient to define an orthogonal rotation matrix $R({\alpha }_1,{\alpha }_2,{\alpha }_3)$ whose columns are related to the eigenvectors of ${\mathcal{M}}_{\rho }^2$. Different parameterizations are possible, and the one used in [2] is

$$R=\left(\begin{array}{ccc}{c}_{{\alpha }_1}{c}_{{\alpha }_2}& s_{{\alpha }_1}{c}_{{\alpha }_2}& s_{{\alpha }_2}\\ -c_{{\alpha }_1}{s}_{{\alpha }_2}{s}_{{\alpha }_3}-s_{{\alpha }_1}{c}_{{\alpha }_3}& c_{{\alpha }_1}{c}_{{\alpha }_3}-s_{{\alpha }_1}{s}_{{\alpha }_2}{s}_{{\alpha }_3}& c_{{\alpha }_2}{s}_{{\alpha }_3}\\ -c_{{\alpha }_1}{s}_{{\alpha }_2}{c}_{{\alpha }_3}+s_{{\alpha }_1}{s}_{{\alpha }_3}& -c_{{\alpha }_1}{s}_{{\alpha }_3}-s_{{\alpha }_1}{s}_{{\alpha }_2}{c}_{{\alpha }_3}& c_{{\alpha }_2}{c}_{{\alpha }_3}\end{array}\right),$$

(24)

with the shorthand notation $c_{\theta }\equiv cos\left(\theta \right)$ and $s_{\theta }\equiv sin\left(\theta \right)$, and the mixing angles reside between $(-{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}})$. Using R, we obtain

$${\mathcal{M}}_H^2=R{\mathcal{M}}_{\rho }^2{R}^T,$$

(25)

where ${\mathcal{M}}_H^2$ is diagonal and can be arranged to have the ordering $m_{H_1}^2<m_{H_2}^2<m_{H_3}^2$. Any of these states can be SM-like, as was considered in [23]; however, in this paper, we consider $H_1$ to be SM-like.

In the N2HDM, the CP-even Higgs mass eigenstates $H_i$ are mixtures of the gauge eigenstates:

$$H_i=\sum _{j=1}^3{R}_{ij}{\rho }_j,$$

(26)

where $i,j=1,2,3$. This introduces a singlet component $|R_{i3}{|}^2$ that features this extension.

Furthermore, the coupling of $H_i$ to the gauge bosons $V=W,Z$ can be derived from the kinetic part of the Lagrangian:

$${\mathcal{L}}_{\mathrm{kin}}={\left(D_{\mu }{\Phi }_1\right)}^{†}\left(D^{\mu }{\Phi }_1\right)+{\left(D_{\mu }{\Phi }_2\right)}^{†}\left(D^{\mu }{\Phi }_2\right)$$

(27)

where $D_{\mu }$ is the covariant derivative. The gauge singlet does not couple directly to gauge bosons; hence, its kinetic term is not included. After electroweak symmetry breaking, field expansions, and rotations from gauge eigenstates to CP-even mass eigenstates, using

$${\rho }_j=\sum _{i=1}^3{R}_{ij}{H}_i,$$

(28)

where the sum is over i and only $j=1,2$ contribute to this coupling based on Equation (27), one can directly derive the result:

$${\mathcal{L}}_{H_{i}VV}^{\mathrm{N}2\mathrm{HDM}}=c\left(H_{i}VV\right){\mathcal{L}}_{hVV}^{\mathrm{SM}},$$

(29)

where the effective coupling is

$$c\left(H_{i}VV\right)={\displaystyle \frac{v_1}{v}}{R}_{i1}+{\displaystyle \frac{v_2}{v}}{R}_{i2}=cos\beta R_{i1}+sin\beta R_{i2}.$$

(30)

For $H_1$, substituting for the elements of R, the coupling becomes

$$c\left(H_{1}VV\right)=cos{\alpha }_{2}cos(\beta -{\alpha }_1).$$

(31)

This shows that the deviation from the SM is driven by the mixing angles. In the limit where ${\alpha }_2\to 0$ and ${\alpha }_1\to \alpha +{\displaystyle \frac{\pi }{2}}$, we recover the 2HDM coupling, while in the limit where ${\alpha }_2\to 0$ and ${\alpha }_1\to \beta$, we recover the SM coupling (the alignment limit). More details on the alignment limit in 2HDMs can be found in [49,50,51].

The Yukawa Lagrangian of the N2HDM, before electroweak symmetry breaking, can be written similarly to the 2HDM, with the singlet field being inert with respect to fermions. Imposing a softly broken ${\mathrm{Z}}_2$ symmetry on the scalar and Yukawa sectors ensures the absence of flavor-changing neutral currents and leads to four possible types of Yukawa assignments:

• Type 1 (T1): All fermions are exclusively coupled to ${\Phi }_2$.
• Type 2 (T2): Up-type quarks couple to ${\Phi }_2$, while down-type quarks and leptons couple to ${\Phi }_1$.
• Type X (TX): Quarks couple to ${\Phi }_2$, and leptons couple to ${\Phi }_1$.
• Type Y (TY): Up-type quarks and leptons couple to ${\Phi }_2$, and down-type quarks couple to ${\Phi }_1$.

As an example, consider the top Yukawa coupling in Type 1. After electroweak symmetry breaking and rotating to mass eigenstates $H_i$, the effective coupling of the lightest CP-even Higgs boson $H_1$ to the top quark relative to the SM can be expressed as

$$C\left(H_{1}t\overline{t}\right)={\displaystyle \frac{R_{12}}{sin\beta }}.$$

(32)

In the SM limit, we have $R_{12}\to sin\beta$, such that $C\left(H_{1}t\overline{t}\right)\to 1$. The full list of effective couplings is presented in Appendix A.

Finally, one can express the quartic couplings in the scalar potential in terms of physical masses, mixing angles, and VEVs:

$$\begin{array}{ccc}\hfill {\lambda }_1& =& {\displaystyle \frac{1}{v^2{c}_{\beta }^2}}\left(-{\displaystyle \frac{m_{12}^2{s}_{\beta }}{c_{\beta }}}+m_{H_1}^2{R}_{11}^2+m_{H_2}^2{R}_{21}^2+m_{H_3}^2{R}_{31}^2\right)\hfill \\ \hfill {\lambda }_2& =& {\displaystyle \frac{1}{v^2{s}_{\beta }^2}}\left(-{\displaystyle \frac{m_{12}^2{c}_{\beta }}{s_{\beta }}}+m_{H_1}^2{R}_{12}^2+m_{H_2}^2{R}_{22}^2+m_{H_3}^2{R}_{32}^2\right)\hfill \\ \hfill {\lambda }_3& =& {\displaystyle \frac{1}{v^2}}\left(-{\displaystyle \frac{m_{12}^2}{c_{\beta }{s}_{\beta }}}+{\displaystyle \frac{1}{s_{\beta }{c}_{\beta }}}(m_{H_1}^2{R}_{11}{R}_{12}+m_{H_2}^2{R}_{21}{R}_{22}+m_{H_3}^2{R}_{31}{R}_{32})+2m_{H^{\pm }}^2\right)\hfill \\ \hfill {\lambda }_4& =& {\displaystyle \frac{1}{v^2}}\left({\displaystyle \frac{m_{12}^2}{c_{\beta }{s}_{\beta }}}+m_A^2-2m_{H^{\pm }}^2\right)\hfill \\ \hfill {\lambda }_5& =& {\displaystyle \frac{1}{v^2}}\left({\displaystyle \frac{m_{12}^2}{c_{\beta }{s}_{\beta }}}-m_A^2\right)\hfill \\ \hfill {\lambda }_6& =& {\displaystyle \frac{1}{v_S^2}}(m_{H_1}^2{R}_{13}^2+m_{H_2}^2{R}_{23}^2+m_{H_3}^2{R}_{33}^2)\hfill \\ \hfill {\lambda }_7& =& {\displaystyle \frac{1}{v{v}_S{c}_{\beta }}}(m_{H_1}^2{R}_{11}{R}_{13}+m_{H_2}^2{R}_{21}{R}_{23}+m_{H_3}^2{R}_{31}{R}_{33})\hfill \\ \hfill {\lambda }_8& =& {\displaystyle \frac{1}{v{v}_S{s}_{\beta }}}(m_{H_1}^2{R}_{12}{R}_{13}+m_{H_2}^2{R}_{22}{R}_{23}+m_{H_3}^2{R}_{32}{R}_{33})\hfill \end{array}$$

(33)

3. Parameter Spaces: Scans and Constraints

Exploring the parameter space of each of the four types was carried out using ScannerS v.2 [48]. We modified the package to enable Latin Hypercube Sampling (LHS) [52,53], which we used along with random sampling. This hybrid scanning technique ensures good coverage of the parameter spaces. In particular, LHS divides each range into N sections, where N is the desired number of samples, and guarantees that combinations from different sections of each parameter are systematically sampled. With that in mind, we collected around 120K samples for the four types, around a third of which were sampled using LHS. The scanned parameters are shown in

Table 1. Parameter ranges for the four types, where mass and VEV variables are in GeV.

Parameter	Min Value	Max Value
$m_{H_a}$	125.09	125.09
$m_{H_b}$	30	1500
$m_{H_c}$	30	1500
$m_A$	30	1500
$m_{H^{\pm }}$	580	1500
$tan\beta$	0.8	20
$c_{H_{a}VV}^2$	0.9	1
$c_{H_{a}tt}^2$	0.8	1.2
$\mathrm{sign}\left(R_{a3}\right)$	−1	1
$R_{b3}$	−1	1
$m_{12}^2$	${10}^{-3}$	$5\times {10}^5$
$v_s$	1	3000

.

Within $\mathtt{ScannerS}$, each point in the parameter space is subjected to a series of validity checks and required to pass the following constraints ($\mathtt{ScannerS}$ interfaces with the obsolete versions HiggsBounds 5 and HiggsSignals 2; hence, we did not utilize them):

• Boundedness: To ensure that the scalar potential is bounded from below as the fields approach infinity, the following conditions need to be met [2]:

  $$\begin{array}{ccc}\hfill \begin{array}{cc}\hfill & {\lambda }_1>0,{\lambda }_2>0,{\lambda }_6>0,\hfill \\ \hfill & \sqrt{{\lambda }_1{\lambda }_6}+{\lambda }_7>0,\hfill \\ \hfill & \sqrt{{\lambda }_2{\lambda }_6}+{\lambda }_8>0,\hfill \\ \hfill & \sqrt{{\lambda }_1{\lambda }_2}+{\lambda }_3+D>0,\hfill \\ \hfill & {\lambda }_7+\sqrt{{\displaystyle \frac{{\lambda }_1}{{\lambda }_2}}}{\lambda }_8\ge 0,\hfill \end{array}& \hfill \mathrm{or}& \hfill \begin{array}{cc}\hfill & {\lambda }_1>0,{\lambda }_2>0,{\lambda }_6>0,\hfill \\ \hfill & \sqrt{{\lambda }_1{\lambda }_6}>-{\lambda }_7\ge \sqrt{{\displaystyle \frac{{\lambda }_1}{{\lambda }_2}}}{\lambda }_8,\hfill \\ \hfill & \sqrt{{\lambda }_2{\lambda }_6}\ge {\lambda }_8>-\sqrt{{\lambda }_2{\lambda }_6},\hfill \\ \hfill & \sqrt{\left({\lambda }_7^2-{\lambda }_1{\lambda }_6\right)\left({\lambda }_8^2-{\lambda }_2{\lambda }_6\right)}>{\lambda }_7{\lambda }_8-(D+{\lambda }_3){\lambda }_6,\hfill \end{array}\end{array}$$

  (34)

  where $D={\lambda }_4-{\lambda }_5$ if ${\lambda }_4>{\lambda }_5$ and zero otherwise.
• Perturbative unitarity: Ensure that the largest eigenvalue corresponding to $2\times 2$ scattering matrices is below the upper limit $8\pi$.

  $$\begin{array}{cc}\hfill & |{\lambda }_3-{\lambda }_4|<8\pi ,\hfill \\ \hfill & |{\lambda }_3+2{\lambda }_4\pm 3{\lambda }_5|<8\pi ,\hfill \\ \hfill & \left|{\displaystyle \frac{1}{2}}\left({\lambda }_1+{\lambda }_2+\sqrt{{({\lambda }_1-{\lambda }_2)}^2+4{\lambda }_4^2}\right)\right|<8\pi ,\hfill \\ \hfill & \left|{\displaystyle \frac{1}{2}}\left({\lambda }_1+{\lambda }_2+\sqrt{{({\lambda }_1-{\lambda }_2)}^2+4{\lambda }_5^2}\right)\right|<8\pi ,\hfill \\ \hfill & |{\lambda }_7|<8\pi ,|{\lambda }_8|<8\pi ,\hfill \\ \hfill & {\displaystyle \frac{1}{2}}\left|a_{1,2,3}\right|<8\pi .\hfill \end{array}$$

  (35)

  where $a_{1,2,3}$ are roots of the cubic equation given in [2] (Equation 3.37).
• Vacuum stability (by $\mathtt{EVADE}$ [54]): Ensure that the EW vacuum is stable or at least metastable and, in that case, long-lived.
• B Physics: Stringent constraints arise from the following processes (see Table 2 and Figure 9 in [55]):

  $$\mathcal{B}(B\to X_s\gamma )=\left(3.32\pm 0.{15}_{\mathrm{stat}+\mathrm{syst}}\right)\times {10}^{-4}\pm 7\%\left(\mathrm{theo}\right).$$

  (36)

  $$\begin{array}{cc}\hfill \mathcal{B}{(B_s\to {\mu }^{+}{\mu }^{-})}_{\mathrm{LHCb}}& =(3.0_{-0.5}^{+0.6})\times {10}^{-9},\hfill \end{array}$$

  (37)

  $$\begin{array}{cc}\hfill \mathcal{B}{(B_d\to {\mu }^{+}{\mu }^{-})}_{\mathrm{LHCb}}& =(1.5_{-1.0}^{+1.2})\times {10}^{-10},\hfill \end{array}$$

  (38)
• Electroweak precision measurements: Restrict the oblique parameters S, T, and U, where, for the simplified scenario $U=0$, we have [55]

  $$S{|}_{U=0}=0.04\pm 0.08,T{|}_{U=0}=0.08\pm 0.07,$$

  (39)

  with a correlation coefficient of $+0.92$.

Finally, after imposing the previous constraints on the generated samples, we interfaced each type with $\mathtt{HiggsTools}$ (HT) [56] via python (we provide the interfacing and analysis codes upon request through our Github page [57]) to confront the model with the latest Higgs data via the HT subpackages $\mathtt{HiggsSignals}$$\mathtt{v}.\mathtt{3}$ (HS) with the HS repository $\mathtt{v}.\mathtt{1.1}$ and $\mathtt{HiggsBounds}$$\mathtt{v}.\mathtt{6}$ (HB) with the HB repository $\mathtt{v}.\mathtt{1.6}$.

4. Results and Discussion

In this section, we present the results of the best-fit analysis, as well as the analysis of the most relevant processes for constraining non-SM Higgs bosons.

4.1. SM-like Higgs Signals

One of the major tasks of the CMS and ATLAS detectors at the LHC is to precisely measure the couplings of the SM-like Higgs to fermions and gauge bosons, which is realized by measuring the production and decay channels. This can quantify any deviations from the S and sets stringent limits on BSM models with a scalar boson that resembles the SM-like Higgs boson to some extent. Indeed, having the same mass is not sufficient to claim SM-like Higgs bosons in a given BSM extension. The predictions of branching ratios and production cross-sections have to be within observed measurements. A practical way to set such limits is through the $\mu$-framework. For a specific production channel i followed by a specific decay channel f,

$${\mu }_{if}={\displaystyle \frac{{({\sigma }_i\times {\mathcal{B}}_f)}^{\mathrm{obs}}}{{({\sigma }_i\times {\mathcal{B}}_f)}^{\mathrm{SM}}}}.$$

(40)

CMS and ATLAS provide data on each measured ${\mu }_i$ and ${\mu }_f$, as well as a combined $\mu$ for all measurements. The latest combined $\mu$ presented by the Particle Data Group (PDG) is [58]

$$\mu =1.03\pm 0.04.$$

(41)

$\mathtt{HS}$ computes ${\chi }^2$ from the signal rates of a given model, normalized by the SM, as (Equation (6) in [59]).

$${\chi }_{\mu }^2={(\widehat{\overrightarrow{\mu }}-\overrightarrow{\mu })}^T{\mathbf{C}}_{\overrightarrow{\mu }}^{-1}(\widehat{\overrightarrow{\mu }}-\overrightarrow{\mu }),$$

(42)

where the $\overrightarrow{\mu }$ vectors contain individual signal strengths as predicted by the models, while $\widehat{\overrightarrow{\mu }}$ represents the corresponding measurement, and ${\mathbf{C}}_{\overrightarrow{\mu }}^{-1}$ is a covariance matrix encoding uncertainties. $\mathtt{HS}$ also computes ${\chi }_m^2$ for the mass of SM-like Higgs. Given that we set the mass of the lightest CP-even Higgs to the observed SM-like Higgs, only ${\chi }_{\mu }^2$ contributes; hence, we drop the subscript in the subsequent analysis. As a reference value, the SM with $m_h=125.09$ GeV gives ${\chi }^2=152.54$, as calculated by $\mathtt{HS}$ with 159 observables. For a given point (p) in the N2HDM, we define

$$\Delta {\chi }^2={\chi }_p^2-{\chi }_{\min }^2,$$

(43)

where ${\chi }_{\min }^2$ is the minimum value in the parameter space, representing the best-fit point.

To facilitate comparison with the 2HDM, we present the results by shifting ${\alpha }_1$ by $-{\displaystyle \frac{\pi }{2}}$ so that the shifted angle is equivalent to the 2HDM convention. In this case,

$$c_{H_{1}VV}=cos(\beta -{\alpha }_1)cos\left({\alpha }_2\right)\stackrel{{\alpha }_1-\frac{\pi }{2}}{\to }\tilde{c}\equiv sin({\alpha }_1-\beta )cos\left({\alpha }_2\right).$$

(44)

We note that, in the allowed data for all types, $cos{\alpha }_2\sim \mathcal{O}\left(1\right)$. However, we emphasize that we are not necessarily in the 2HDM limit of the N2HDM since ${\alpha }_3$ is not always close to zero in the valid parameter space.

Figure 1 shows the $\mathtt{HS}$ results and the ${\chi }^2$ analysis for the four types of N2HDMs. Points that are colored red are ruled out by $\mathtt{HS}$ since $\Delta {\chi }^2>5.99$. The yellow points represent the $95\%$ C.L. for which $\Delta {\chi }^2\le 5.99$ (see Table 4 in [59]), while the green points are for cases where $\Delta {\chi }^2\le 2.3$, which fall in the $68\%$ C.L. The black stars represent the best-fit points where $\Delta {\chi }^2$ is minimum. For completeness, the gray points represent cases that are ruled out by $\mathtt{HB}$.

For T1, we observe in the top-left panel of Figure 1 that the best-fit point corresponds to a small value of $tan\beta =0.96$, and that it is slightly shifted to the negative side $\tilde{c}=-0.057$, which coincides with the results of the 2HDM presented in [60] and the case study in [56]. Most of the green points reside within $-0.1<\tilde{c}<0$ with $0<tan\beta <3$, while the yellow points expand to the positive side where $\tilde{c}>0.05$, especially as $tan\beta$ becomes larger than 2. As $tan\beta$ becomes smaller, the positive side of $\tilde{c}$ is disallowed by HS, as indicated by the red points.

Next, in T2, which is shown in the top-right panel, the best-fit point is located at $\tilde{c}=0.013$ and corresponds to $tan\beta =4.95$, which is larger than T1. Points that fall within the $68\%$ C.L. region are located on the positive side, where $\tilde{c}\sim 0.05$ for $tan\beta <2$, and approach zero as $tan\beta$ becomes larger. On the other hand, the yellow region spans $0.02<\tilde{c}<0.07$ for small values of $tan\beta$ and converges near $\tilde{c}\sim 0$ as $tan\beta$ increases.

As for TX, in the bottom-left panel, we can see that the green points shifted slightly to negative for smaller values of $tan\beta$ where $\tilde{c}\sim -0.04$, and as $tan\beta$ increases, this region becomes more symmetric around $\tilde{c}=0$. The yellow points, on the other hand, extend between $-0.06<\tilde{c}<0.6$ and $tan\beta <2$ and shirk as it increases. The best-fit point lies at $tan\beta =4.37$, for which $\tilde{c}=0.008$.

In TY, as shown in the bottom-right panel, the allowed region is centered at zero. Specifically, the green points span $-0.01<\tilde{c}<0.04$ for $0.8<tan\beta <3$, while the yellow points extend slightly further and up to $\tilde{c}\sim 0.05$. Both regions become closer to $\tilde{c}=0$ as $tan\beta$ increases. The best-fit point has $tan\beta =4.66$ and $\tilde{c}=0.009$.

To better understand the specific measurements affecting the regions surrounding the best-fit points, we computed, using $\mathtt{HS}$, individual ${\chi }^2$ for two points at $tan{\beta }_{{\chi }_{\min }^2}\pm \delta$, where $\delta \le tan{\beta }_{{\chi }_{\min }^2}\times 2\%$, and $\tilde{c}={\tilde{c}}_{{\chi }_{\min }^2}\pm \gamma$, where $\gamma \le 0.05$, depending on the availability of the data. Then, we computed $\Delta {\chi }^2$ between the two selected points to find measurements that lead to significant deviations from the value of ${\chi }_{\min }^2$. We find that these regions are mostly affected by searches for

• $pp\to h\to VV\to 4l$ [61,62,63,64],
• $pp\to ht\overline{t}\to (h\to b\overline{b})(t/\overline{t}\to \mathit{semileptonic})$ [65],
• $pp\to Vh\to (h\to b\overline{b})(V\to ll/l\nu /\nu \nu )$ [66,67],
• $pp\to h\to {\tau }^{+}{\tau }^{-}$ [68,69],
• $pp\to h\to \gamma \gamma$ [70,71]

Figure 1 also indicates that deviations from the SM are allowed in the four types; however, this is more pronounced in T1, for which the allowed region does not converge sharply to the alignment limit with the increase in $tan\beta$, as is the case in the other three types. This distinct feature is common in 2HDMs, as indicated previously. However, it should be mentioned that the properties of SM-like Higgs in the N2HDM can differ from those in the 2HDM due to the presence of an additional singlet component, $|R_{13}{|}^2$. In particular, in terms of the model parameters, we observe that, in T1 and TX, the singlet component of SM-like Higgs is restricted to values below $10\%$, while it is below $15\%$ in T2 and below $18\%$ in TY. Furthermore, one can set upper and lower limits on the angles ${\alpha }_1$ and ${\alpha }_2$ based on the allowed results, as shown in

Table 2. Allowed ranges of the mixing angles ${\alpha }_1$ and ${\alpha }_2$.

	${\alpha }_1$ (Min, Max)	${\alpha }_2$ (Min, Max)
Type 1	$-1.556$, $1.563$	$-0.3139$, $0.2787$
Type 2	$0.6649$, $1.502$	$-0.4026$, $0.3853$
Type X	$0.6908$, $1.466$	$-0.3115$, $0.2923$
Type Y	$0.7366$, $1.469$	$-0.4235$, $0.4339$

, while ${\alpha }_3$ is found to be allowed to take the full range specified in the previous section.

4.2. Bounds on the Additional CP-Even Higgs Bosons

ATLAS and CMS systematically search for non-SM Higgs bosons decaying into lighter bosons or fermions in several final states. This case is relevant for this work, since we are considering the ordering $m_{H_1}^{\mathrm{SM}}<m_{H_2}<m_{H_3}$. However, experimental results usually assume the NWA; hence, they can set limits on the additional Higgs bosons in the N2HDM provided that ${\displaystyle \frac{{\Gamma }_H}{m_H}}\ll 1$, where $\Gamma$ is the total width of the additional Higgs. Nevertheless, certain searches presented results valid beyond the NWA, and these are taken into account, as will be discussed later. For each additional Higgs, $\mathtt{HT}$’s subpackage $\mathtt{HiggsPredictions}$ ($\mathtt{HP}$) computes the production cross-sections in the effective coupling approximation, including QCD corrections. We take all main production channels into account, especially the fact that Vector Boson Fusion (VBF) and production in association with a vector boson (HV) can be larger than gluon fusion (ggH) in regions where the coupling of the additional Higgs to up-type quarks is very small compared to its coupling to gauge bosons. Similarly, production in association with bottom and anti-bottom quarks (bbH) can be larger than ggH in certain regions. All these effects have been taken into account in our analysis, within the precision provided by $\mathtt{HT}$ and its subpackages.

Furthermore, for each instance of input parameters, $\mathtt{HB}$ determines the most sensitive measurement. This is implemented by computing the model’s prediction for each relative $\sigma \times {\mathcal{B}}_i$ divided by the corresponding expected limits. The limit that maximizes this ratio is considered the most sensitive. Next, it computes the observed ratio (the model’s prediction for a specific decay channel divided by the corresponding observed limits). If the ratio is larger than 1, then it determines that the point is excluded at the $95\%$ C.L. [72,73].

For each type, we discuss the most sensitive channels based on $\mathtt{HB}$, paying special attention to heavy Higgs resonances decaying into a pair of bosons or fermions. All points already pass the constraints from $\mathtt{HS}$ and the set of constraints described in the previous section. Since the four types are generally affected by a common set of measurements, we present the results with respect to each class of relevant measurements. Moreover, in the presented figures, we use the HEP Inspire biographical code.

The most stringent bounds on the additional CP-even Higgs bosons come from the classes of searches listed below, where we define $H_1^{\mathrm{SM}}\equiv h$, and H is an additional non-SM Higgs boson.

4.2.1. Class: $pp\to H\to hh$

An important class of LHC searches is that for an additional Higgs decaying into two SM-like Higgs bosons. In the analyzed parameter spaces of the four types, the most sensitive ones are

• $pp\to H\to hh\to {\tau }^{-}{\tau }^{+}b\overline{b}$ [74]
• $pp\to H\to hh\to b\overline{b}\gamma \gamma$ [75]
• $pp\to H\to hh\to \mathit{fermions}/\mathit{bosons}$ [76]

where the results corresponding to this class appear in the left panels of Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. In these panels, green points are allowed, and gray points are ruled out by a different class (i.e., $pp\to H\to VV$).

We observe that the limit from ATLAS (139 $f{b}^{-1}$) on the narrow resonance production of a pair of SM Higgs bosons in the ${\tau }^{-}{\tau }^{+}b\overline{b}$ final state [74] ($\mathtt{ATLAS}:\mathtt{2022}\mathtt{xzm}$) imposes strong constraints, as it contributes the most to ruling out points in the parameter spaces. This constraint is relevant for a mass range of $H_2/H_3$ between 251 and 1600 GeV. It is very sensitive to $hh$ searches, given the relatively low background accompanied by a branching ratio of $\mathcal{B}(hh\to {\tau }^{-}{\tau }^{+}b\overline{b})\sim 0.073$. The upper limits observed in $\sigma \times \mathcal{B}$ range from 0.9 pb to 0.021 pb. As can be seen in the left panels of Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, all types are affected by this measurement, and the corresponding ruled-out points appear in blue. In particular, we observe that this search is deemed to be the most sensitive for values of mass above 400 GeV. Regions below that are affected by other searches that will be discussed shortly. Furthermore, we note that in T1, for $H_2$ with a mass of around 900 GeV, Figure 2 (left panel) shows a few black points where $\mathtt{HB}$ selects this search to be the most sensitive; however, these points reside beyond the NWA.

Moreover, ATLAS (139 $f{b}^{-1}$) provided limits on new hypothetical heavy scalars in the mass range of 251 GeV to 1000 GeV, which are set through searches in the final states to $b\overline{b}$ quarks and a pair of photons [75] ($\mathtt{ATLAS}:\mathtt{2021}\mathtt{ifb}$). The observed upper limits range from 0.64 pb to 0.044 pb. This search is responsible for ruling out parameter points (orange) with large $\sigma \times \mathcal{B}$ and $m_{H_2}<400$ GeV, as seen in the left panels of Figure 2 (T1), Figure 4 (T2), and Figure 6 (TX).

Finally, for this class, limits on a new heavy scalar boson in the mass range between 270 and 3000 GeV are obtained from the results of CMS (35.9 fb⁻¹) [76] ($\mathtt{CMS}:\mathtt{2018}\mathtt{ipl}$), which searched for a new scalar decaying into two SM-like Higgs bosons, one of which subsequently decays into $b\overline{b}$, while the other one could decay into $b\overline{b}/{\tau }^{+}{\tau }^{-}/VV/\gamma \gamma$. The observed upper limits range from 0.68 pb to 0.002 pb. The effect of this measurement is visible in the left panel of Figure 3 for $H_3$ of T1 and appears in yellow, where the mass is below 400 GeV.

4.2.2. Class: $pp\to H\to VV$

A general search conducted at the LHC is that for a heavy resonance decaying into a pair of gauge bosons, which would subsequently decay into fermions. While some searches considered different fermionic final states, other searches focused on semileptonic or leptonic final states. The majority of constraints come from this class of measurements, and we find that the most sensitive ones are

• $pp\to H\to VV\to \mathit{fermions}$ [77,78,79]
• $pp\to H\to VV/Vh\to \mathit{semileptonic}$ [80,81]

where the results corresponding to this class appear in the right panels of Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. In these panels, green points are allowed, and gray points are ruled out by a different class (i.e., $pp\to H\to hh$).

ATLAS (36 $f{b}^{-1}$) performed a general search for a heavy scalar resonance producing two bosons, which subsequently decay into fermions [77] ($\mathtt{ATLAS}:\mathtt{2018}\mathtt{sbw}$); this search is relevant for a mass range between 300 GeV and 3000 GeV. For ggH, the observed upper limits range from 0.38 (300 GeV) to 0.0013 (3 TeV), while for VBF, the observed upper limits range from 0.13 (500 GeV) to 0.0033 (3 TeV). The effects on the four types are shown in red, specifically in the right panels of each plot corresponding to $H\to VV$ in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, where $m_H\ge 500$ GeV. Some black points where $m_H>700$ GeV, for which this measurement was deemed the most sensitive, do not satisfy the NWA and are shown in the right panels of Figure 2 and Figure 3 (T1), Figure 4 and Figure 5 (T2), Figure 6 (TX), and Figure 8 and Figure 9 (TY).

Moreover, ATLAS (139 $f{b}^{-1}$) presented results concerning the production of a resonant non-SM Higgs ($H_2/H_3$) decaying into two gauge bosons, which subsequently decay into leptons [78] ($\mathtt{ATLAS}:\mathtt{2020}\mathtt{tlo}$). In models with two doublets, the upper limits for ggH and in the mass range between 200 GeV and 400 GeV are given in the NWA, where the range is from 0.11 pb to 0.047 pb. For larger mass values, limits beyond the NWA are included, up to ${\displaystyle \frac{\Gamma }{m_H}}=0.15$. For VBF, the upper limits range from 0.031 pb (210 GeV) to 0.0017 pb (2 TeV). All points affected by this measurement are shown in purple, where we can see that the effects are severe in the mass regions below 600 GeV in the right panels of Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

Moreover, the search by CMS [79] ($\mathtt{CMS}:\mathtt{2018}\mathtt{amk}$) at 35.9 $f{b}^{-1}$ for new resonances decaying into ZZ bosons, which subsequently decay into $4l$, $2l2q$, or $2l2\nu$, is valid for mass regions between 130 GeV and 3 TeV. It covers a wide range of widths: $0\le {\displaystyle \frac{\Gamma }{m_H}}\le 0.3$. For ggH, in the NWA, the upper limits range from 0.235 pb (130 GeV) to 0.0012 pb (3 TeV), while for ${\displaystyle \frac{\Gamma }{m}}$ = 0.3, they range from 0.104005 pb (130 GeV) to 0.006 pb (3 TeV). For VBF, they range from 0.166 pb (130) to 0.0011 pb (3 TeV). As for larger widths, where ${\displaystyle \frac{\Gamma }{m}}=0.3$, the range is between 0.0243 pb (130) and 0.0018 pb (3 TeV). We can see in the right panels of Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 (pink points) that this measurement rules out some vertical points corresponding to a mass range below 200 GeV.

Additionally, CMS (139 fb⁻¹) searched for di-boson resonances [80] ($\mathtt{CMS}:\mathtt{2021}\mathtt{klu}$), where the resonance has a mass range between 1 TeV and 4.5 TeV. The final states are those containing leptons and hadrons. In the models, we observe that all ruled instances are associated with production via VBF for which the upper limits range from 0.0086 pb (1 TeV) to 0.00016 pb (4.5 TeV), as can be seen in Figure 5 (T2), Figure 7 (TX), and Figure 9 (TY). Some points in each type, especially for $H_3$, might evade this constraint due to the breakdown of the NWA, and these are shown in black.

Finally, for this class, the ATLAS detector (with 139 fb⁻¹) reported limits on heavy resonances decaying into a pair of gauge bosons in the semileptonic final states (i.e., one V decaying into leptons, and the other into hadrons) [81] ($\mathtt{ATLAS}:\mathtt{2020}\mathtt{fry}$). The relevant mass range for the heavy scalar is between 300 and 5000 GeV. The search considered different production topologies; the one relevant here is that of VBF, as selected by the $\mathtt{HB}$ criteria. The upper limit ranges from 0.604 pb to 0.00024. The effects of this search are seen as yellow points in the right panels of Figure 5 (T2) and Figure 8 and Figure 9 (TY) where $m_H>600$ GeV. Some points, indicated by black points, may pass this constraint due to having large widths, especially for TY ($H_2,H_3$) and T2 ($H_2,H_3$).

4.2.3. Class: $pp\to H\to f\overline{f}$

ATLAS and CMS searched for additional Higgs bosons decaying into fermions. CMS carried out a dedicated search (138 fb⁻¹) in the $\tau \overline{\tau }$ final states [82] ($\mathtt{CMS}:\mathtt{2022}\mathtt{goy}$). This search affects a mass range from 60 to 3500 GeV and considers production via ggH and bbH. The upper limits range from $\mathcal{O}\left(10\right)$ pb (60 GeV) to 0.003 pb (3.5 TeV). In our analysis, we find that this search affects the low-mass ranges of $H_2$ and $H_3$, especially for TX.

4.3. Constraints on Effective Couplings

To understand how the combined constraints affect the parameters of the four types, we delve into their effects on the effective couplings, focusing on $H_2$ and $H_3$. Not only do they determine the size of the branching ratios to different SM particles, but they also determine which production channel is dominant. In particular, sizable couplings to up-type quarks accompanied by small couplings to vector bosons mean that ggH is larger than VBF, and vice versa. Also, since, in T1, the effective couplings to fermions are equal, ggH is always larger than bbH. This is not the case for T2, for which bbH can be sizable, especially in the smaller mass range.

Restricting ourselves to the case where all additional Higgs bosons are narrow (including A and $H^{\pm }$), Figure 10 and Figure 11 show the effective couplings of $H_2$ and $H_3$ where the color code represents the mass of the additional Higgs boson, sorted in descending order such that the smaller values are plotted on top of heavier ones. We can see in Figure 10 that the smaller range of mass is associated with small effective couplings, in general. We can further see that the effective couplings to bottom quarks (T2, TY) and $\tau$ leptons (T2, TX) are closer to zero for the small mass range. As the mass increases, these effective couplings open up and can be significantly larger than the effective couplings to gauge bosons. This constraint is not seen for the $H_3$ depicted in Figure 11, where the small mass range can still have effective couplings to fermions much larger than those to gauge bosons.

The mass distribution shows distinctive patterns across the parameter space, with lower masses (400–600 GeV) concentrating in regions of smaller couplings for $H_2$, while $H_3$ exhibits a more spread distribution, forming characteristic triangular patterns in T1 that become increasingly asymmetric in other types.

In each figure, the subplots for T1, T2, TX, and TY are arranged to highlight how different types cluster in different coupling regions. For example, in T2, one observes that $c_{H_{2}bb}$ can become as negative as around $-5$, indicating a much larger deviation from SM-like behavior compared to T1. Similarly, in TY, the range of $c_{H_{2}bb}$ reaches around $-8.414$, and in TX, the range of $c_{H_2\tau \tau }$ can reach approximately $-8.811$. Meanwhile, the couplings to vector bosons (e.g., $c_{H_{2}VV}$) remain in a tighter range, rarely below $-0.376$ or above $0.305$ in all four types. The same pattern is even more pronounced for $H_3$, where, in T2, the effective coupling to the bottom quark can drop to $-13.499$ or rise above 7, and in TY, it can reach about $10.148$. These wide spreads for fermion couplings contrast with the relatively narrower window for $c_{H_{3}VV}$, spanning around $[-0.324,0.331]$.

These coupling patterns have direct implications for the dominant production mechanisms, with the large variations in fermion couplings particularly affecting the interplay between gluon fusion and bottom-quark-associated production across different mass ranges.

Finally,

Table 3. Ranges of allowed effective couplings for $H_2$ and $H_3$.

Particle	Type	${\mathit{c}}_{{\mathit{H}}_{\mathit{i}}\mathit{VV}}$	${\mathit{c}}_{{\mathit{H}}_{\mathit{i}}\mathit{uu}}$	${\mathit{c}}_{{\mathit{H}}_{\mathit{i}}\mathit{dd}}$	${\mathit{c}}_{{\mathit{H}}_{\mathit{i}}\mathit{ll}}$
$H_2$	T1	$[-0.270,0.267]$	$[-0.248,1.046]$	$c_{H_{2}dd}=c_{H_{2}uu}$	$c_{H_{2}ll}=c_{H_{2}uu}$
	T2	$[-0.376,0.300]$	$[-0.309,1.042]$	$[-5.359,0.260]$	$c_{H_{2}ll}=c_{H_{2}dd}$
	TX	$[-0.291,0.259]$	$[-0.285,1.039]$	$c_{H_{2}dd}=c_{H_{2}uu}$	$[-8.811,0.255]$
	TY	$[-0.357,0.305]$	$[-0.282,1.045]$	$[-8.414,0.258]$	$c_{H_{2}ll}=c_{H_{2}uu}$
$H_3$	T1	$[-0.261,0.284]$	$[-1.037,1.053]$	$c_{H_{3}dd}=c_{H_{3}uu}$	$c_{H_{3}ll}=c_{H_{3}uu}$
	T2	$[-0.318,0.329]$	$[-1.042,1.019]$	$[-13.499,7.429]$	$c_{H_{3}ll}=c_{H_{3}dd}$
	TX	$[-0.278,0.298]$	$[-1.047,1.047]$	$c_{H_{3}dd}=c_{H_{3}uu}$	$[-5.124,6.119]$
	TY	$[-0.324,0.331]$	$[-1.010,1.036]$	$[-7.138,10.148]$	$c_{H_{3}ll}=c_{H_{3}uu}$

presents the ranges allowed for effective couplings. One can see that the effective couplings to bottom quarks and tau leptons can significantly deviate to larger values compared with SM-like Higgs with the same mass. On the other hand, the coupling to top quarks can be smaller or close to SM-like Higgs of the same mass, while the effective couplings to vector bosons are always smaller than SM-like Higgs of the same mass. As indicated earlier, T2 and TY can reach substantial negative values for $c_{H_{i}bb}$, down to $-5.359$ and $-8.414$ for $H_2$, and even lower for $H_3$. The effective coupling with tau has sizable variations: $[-8.811,0.255]$ in TX for $H_2$, and $[-5.124,6.119]$ in TX for $H_3$. In contrast, $c_{H_{i}VV}$ remains within comparatively narrow bounds: $[-0.376,0.305]$ for $H_2$ in T2, and $[-0.324,0.331]$ for $H_3$ in TY.

4.4. Constraints on Pseudoscalar and Charged Higgs Bosons

As mentioned in the Introduction, the pseudoscalar A and the charged Higgs pair $H^{\pm }$ have the same structure as in the 2HDM. However, these can be indirectly affected by the overall constraints on the N2HDM. Since the effective couplings of A and $H^{\pm }$ depend on $tan\beta$, Figure 12 and Figure 13 show the allowed parameter spaces in the $m_{A/H^{\pm }}$-$tan\beta$ planes for the four types.

The distributions in the $m_A$-$tan\beta$ plane exhibit distinctive patterns across all types. T1 shows the widest mass range, allowing pseudoscalar masses as low as 97.8 GeV and extending up to about 1460 GeV, with $tan\beta$ values reaching approximately 12. T2 shows a more constrained parameter space, with $m_A$ starting at around 446 GeV, and most allowed points are concentrated at lower values of $tan\beta$. TX shows an intermediate-mass range starting from 350 GeV, while TY has a higher mass threshold of around 470 GeV. The charged Higgs masses follow similar patterns, with T1 showing the broadest range starting from 179.7 GeV, while T2 and TY have higher thresholds of around 600 GeV. TX allows for relatively lighter charged Higgs, with masses starting from 268 GeV. In all types, the parameter space becomes increasingly sparse at higher $tan\beta$ values, with the densest populations observed below $tan\beta \approx 4$.

In T1, the pseudoscalar A with a mass range between 225 and 1000 GeV is mainly affected by the CMS search for A decaying into Z ($Z\to ll/\nu \nu$) and an SM-like Higgs boson ($h\to bb$) [83]. As for T2, we observe that some points are ruled out because A is inconsistent with measurements from ATLAS [84], where it decays into a heavy Higgs (decaying into top pairs) and a Z boson. This is relevant for mass ranges for A between 450 and 1200 GeV and heavy Higgs between 350 and 800 GeV. Additionally, for TX, we observe that A is affected by the likelihood analysis presented by the CMS search for additional Higgs bosons decaying into a pair of $\tau$s [82] for a mass between 160 GeV and 3500 GeV. We also find that, for TY, A is mostly affected by ATLAS searches for A decaying into heavy Higgs and Z bosons [84,85]. The former is for $m_A>800$ GeV and $m_H>300$ GeV, while the latter is for $230\le m_A\left(\mathrm{GeV}\right)\le 800$ and $130\le m_H\left(\mathrm{GeV}\right)\le 700$.

On the other hand, for all types, the charged Higgs boson is affected by the results from ATLAS [86] searching for $pp\to tb{H}^{+}\to tbtb$, where $200\le m_{H^{\pm }}\left(\mathrm{GeV}\right)\le 2000$ ($\mathtt{HT}$ extends the range to start from $m_{H^{\pm }}\ge 145$ GeV), in final states consisting of jets and one electron or muon. For a mass range smaller than 145 GeV, which is only relevant for T1 and TX, the constraint on the charged Higgs is mainly set, based on $\mathtt{HB}$ selection, by the ATLAS search for charged Higgs decaying into $\tau {\nu }_{\tau }$ [87], where the charged Higgs is produced in decays of the top quark.

4.5. Prospects for Future Searches

The LHC is expected to be upgraded to the High-Luminosity LHC (HL-LHC) [88] by 2030, allowing for the allocation of 3000 fb⁻¹ of data during its operation. This will enable more precise measurements of the properties of the SM-like Higgs boson and expand the discovery reach for additional scalars. In particular, the projected precision for measuring the main Higgs production channels ranges from 1.6% (ggH) to 5.7% (WH). Meanwhile, dominant decay modes are expected to be probed with precisions of 2.6% ($\gamma \gamma$), about 2.9% (ZZ, $W^{+}{W}^{-}$, ${\tau }^{+}{\tau }^{-}$), and 4.4% ($b\overline{b}$). Rare decays to ${\mu }^{+}{\mu }^{-}$ and $Z\gamma$ are expected to be observed, but with larger uncertainties. In terms of coupling modifiers, the projected uncertainties are also at the a-few-percent level [89,90]. These measurements will have important implications for the parameter spaces of BSM extensions such as the N2HDM.

Furthermore, searches for BSM Higgs bosons will considerably improve, extending the reach in probing mass ranges and couplings by up to 50% [89,90,91]. For instance, limits on the decay of a heavy scalar resonance into a pair of Z bosons are anticipated to improve by a factor of ten. Processes such as $A\to ZH$ and $H\to ZA$ will become important probes, especially in regions with substantial mass splitting, and decays to a pair of $\tau$ leptons will provide complementary channels. Parameter regions farther from the alignment limit can be tested via $pp\to A\to Zh\to \ell \ell b\overline{b}$, which is expected to gain sensitivity in further upgrades at $\sqrt{s}$ between 14 and 27 TeV (the High-Energy LHC). The overall implications for the different types of N2HDMs will be significant in terms of restricting the allowed SM-like Higgs couplings, the mixing angles, and the singlet component, as well as the allowed mass range and properties of the additional Higgs boson. A detailed analysis of such effects will be the subject of future work.

5. Conclusions

All in all, we have considered the broken phase of the N2HDM where the two Higgs doublets (${\Phi }_1$ and ${\Phi }_2$) and the singlet (${\Phi }_S$) acquire VEVs. The model admits two discrete symmetries, one of which is spontaneously broken by the singlet VEV, while the other is softly broken by the $m_{12}^2$ parameter and extends to the Yukawa sector. This brings about four types that encode the different possibilities of Yukawa couplings between the Higgs doublets and SM fermions (i.e., Type 1, Type 2, Type X, and Type Y). The parameter spaces have been subjected to limits from theory and observations. The model was interfaced with the latest Higgs data repositories of the public code $\mathtt{HiggsTools}$ and its subpackages.

By performing a statistical ${\chi }^2$ analysis using $\mathtt{HS}$, we identified the best-fit point for each type. We found that T1 is slightly shifted from the alignment limit, with $tan\beta$ being as small as 0.96, while the best-fit points of the other types reside within the alignment limit with moderate values of $tan\beta \sim 5$. We analyzed the specific CMS and ATLAS measurements affecting the neighboring regions of the best-fit points, leading to deviations from ${\chi }_{\min }^2$. Furthermore, we determined the effects of the constraints on the parameters of the model, namely, the singlet component of the SM-like Higgs and the mixing angles.

We have also analyzed the bounds on the additional Higgs bosons using $\mathtt{HB}$, particularly from searches for the resonance production of a pair of bosons or fermions via an additional Higgs. Additionally, we have shown the effects of the constraints on the effective couplings of the additional Higgs bosons and their masses, including the pseudoscalar and charged Higgs pair.

In particular, we observe that the SM-like Higgs signals are mainly constrained by searches involving $h\to VV,\gamma \gamma ,b\overline{b},{\tau }^{+}{\tau }^{-}$. The mixing angles ${\alpha }_1$ and ${\alpha }_2$ are subject to the ranges summarized in Table 2, while the singlet component of SM-like Higgs is below $10\%$ in T1 and TX, below $15\%$ in T2, and up to $18\%$ in TY. Concerning the additional Higgs bosons, the most restrictive classes of measurements are those searching for heavy resonances decaying into $hh$ and $VV$.

We find that the constraints on the effective couplings to vector bosons, $c_{H_{i}VV}$, are quite stringent in all types, requiring them to lie in comparatively narrow ranges, such as $-0.376\le c_{H_{2}VV}\le 0.305$ in T2 or $-0.324\le c_{H_{3}VV}\le 0.331$ in TY. Meanwhile, the effective couplings to fermions can be significantly larger than in the SM; for instance, in T2 and TY, the bottom-quark coupling $c_{H_{3}bb}$ can be as low as about $-13.5$ and as high as about $10.1$, respectively. Regarding the pseudoscalar A, the allowed mass range in T1 spans from about 98 GeV to about 1460 GeV; in T2, it starts around 446 GeV; in TX, around 350 GeV; and in TY, around 470 GeV. The corresponding $tan\beta$ values can be as small as around 0.8 or as large as 12, depending on the type.

Finally, the overall constraints from the Higgs data show the vital role played by searches for additional Higgs bosons. It is important to keep up-to-date with results from recent and future LHC runs. These are expected to be systematically included in updates by the $\mathtt{HT}$ group, and our interface/analysis code can be readily utilized to inspect any further effects on the four types of N2HDMs.