Quantum interference effects in Higgs boson pair-production beyond the Standard Model

New physics frameworks like the Next-to-Minimal Supersymmetric Standard Model and the Next-to-2-Higgs-doublet Model contain three neutral CP-even Higgs bosons. It is possible for the heavier two of these states to have masses identical to each other, which can result in a sizeable quantum interference between their propagators in processes they mediate. For both these models, we study the impact of such interference on the pair-production of the lightest of the three scalars, which we identify with the observed 125 GeV Higgs boson, in the gluon-fusion channel at the Large Hadron Collider (LHC). We find that the inclusion of these effects can substantially alter the cross section, compared to its value when they are ignored, for this process. Our results illustrate the importance of taking possible quantum interference effects into account not only when investigating the phenomenology of extended Higgs sectors at the future Run(s) of the LHC, but also when imposing its current exclusion bounds on the parameter spaces of these models.


Introduction
Pair-production of the Higgs boson state, h obs , discovered in 2012 [1,2] is a key process for measuring the Higgs self-coupling at the Run 3 of the Large Hadron Collider (LHC), as well as at its now approved high-luminosity upgrade (HL-LHC). This process represents a direct probe of the mechanism of electroweak symmetry breaking (EWSB), since the Higgs self-coupling enters the Higgs potential directly. Accessing it experimentally will thus be of extreme importance in order to understand whether mass generation in nature occurs within the Standard Model (SM) or in some scenario incorporating new physics.
In a beyond-the-SM (BSM) framework containing an extended Higgs sector, the phenomenology of the pair-production process of the h obs candidate, i.e., the Higgs boson with mass lying near 125 GeV, can deviate significantly from that in the SM due to two main reasons. First, the h obs self-coupling gets modified from its predicted value in the SM owing to the mixing between various interaction eigenstates. Secondly, the additional Higgs states also enter the resonant channel, so that the other Higgs trilinear couplings appearing in the Lagrangian of the model also come into play. While there exists plenty of literature on h obs pair-production in BSM scenarios at the various energy and luminosity stages of the LHC, most often this is limited to frameworks wherein only one CP-even companion to the SM-like Higgs state exists, like the 2-Higgs Doublet Model (2HDM) [3] or the Minimal Supersymmetric Standard Model (MSSM) [4,5].
In the MSSM, the requirement for the lighter of its two scalar states, H 1 , to be a good h obs candidate pushes the heavier scalar, H 2 , as well as the solo pseudoscalar, A, up into the so-called decoupling regime [6], where they have identical masses. If the NMSSM Higgs sector is CP-violating, all the interaction eigenstates can mix together to yield three CP-indefinite physical states, with the two nearly mass-degenerate heavy states now labelled H 2 and H 3 . When the mass-splitting between these two is comparable to the sum of their widths, a description of the intervening propagators which takes into account the imaginary parts of the one-loop self-energies, alongside the customary real parts, becomes necessary [7,8]. This is because the imaginary off-diagonal entries of the Higgs propagator matrix can induce quantum interference between these states, so that the one produced in, for example, gluon-fusion can potentially oscillates into the other one before decaying into a given SM final state. This can significantly alter not only the total production cross section but also the shape of the differential cross section distribution for that final state [9,10].
In the 2-Higgs-doublet model (2HDM), obtained by simply adding an additional Higgs doublet to the SM, which results again in two scalar and one pseudoscalar Higgs states, a mass-degeneracy between H 2 and A is not a precondition for the H 1 to have properties identical to the h obs . It is nonetheless a possibility not ruled out by any experimental results, and the aforementioned interference effects can become significant in this model also if it has a CP-violating Higgs sector with m H 2 ≈ m H 3 . In a BSM scenario containing three or more CP-even Higgs bosons, the quantum interference effects can appear in processes involving Higgs propagators without the need to invoke CP-violation. A minimal realisation of such a scenario would be the extension of the two models mentioned above by a singlet Higgs field, resulting in an extra scalar state in their Higgs sectors.
In the context of Supersymmetry, adding a complex singlet Higgs field to the MSSM can address some of its theoretical and phenomenological shortcomings, resulting in the so-called Next-to-MSSM (NMSSM) [11,12,13,14]. In this model, some particular configurations of the free parameters can yield a SM-like H 1 along with H 2 and H 3 that are nearly mass-degenerate. We have previously investigated the aforementioned interference effects in the NMSSM, in the scenario where H 1 and H 2 are mass-degenerate [15], as well as in the alternative scenario with m H 2 ≈ m H 3 [16]. The first study pertained to the production process for the γγ final state and the second to that of τ + τ − at the LHC. Both these studies found the results from the calculation embedding the full Higgs propagator matrix to be notably different from the ones using the standard approximation where only one term containing a Breit-Wigner (BW) propagator corresponding to each of the Higgs bosons appears in the amplitude expression. It was also shown in those papers that the expected mass resolutions of the respective final states at the LHC may, however, not allow it to disentangle the two Higgs states from each other.
In this article, we investigate the implications of the quantum interference on the gluon-initiated pair-production of the SM-like H 1 state of the NMSSM at the LHC with √ s = 14 TeV, and also of its non-Supersymmetric counterpart, the Next-to-2HDM (N2HDM). The latter model is obtained by introducing a real singlet Higgs field into the 2HDM, and while it is phenomenologically similar to the NMSSM, a crucial advantage the N2HDM has is that the physical Higgs boson masses can themselves be the input parameters. This grants us the freedom of setting the H 2 and H 3 masses exactly equal and assessing the impact of this maximal mass-degeneracy on the said process. This model additionally allows us to analyse how the various Higgs couplings govern the relative sizes of the interference effects, so that the general inferences can be extended to other multi-Higgs BSM scenarios.
The article is organised as follows. In the next section we briefly revisit the Higgs pair-production process at the LHC. In section 3 we discuss some details of the NMSSM and the N2HDM, as well as of our numerical computational tool. In section 4 we present our analysis and discuss its results. In section 5 we conclude our findings.
2 Higgs pair-production at the LHC The cross section for the (inclusive) process pp → H i H j , where i, j = 1, ..., N for a model with N CP-even Higgs bosons but without any additional particle content beyond the SM, can be written at the leading order (LO) as where g(x 1 ) and g(x 2 ) are the parton distribution functions (PDFs) of the two incoming gluons having squared centre-of-mass (CM) energyŝ = x 1 x 2 s, given in terms of the total CM energy, s, of the pp system, and by defining τ ≡ŝ s = x 1 x 2 . The amplitude-squared in Eq. (1) can be written, following the notation of Ref. [17], as where denotes the contribution from the Higgs-mediated triangle loop diagram, Fig. 1 (left), and 2 refers to the quark-box diagram, Fig. 1 (right).
The coefficient corresponding to the box contributions in Eq. (2) is written in terms of the Yukawa couplings as The form factor F 2 corresponds to the case when the gluons have a combined total spin of S z = 0 along the proton beam, while G 2 refers to the case with S z = 2. The full expressions for F 2 and G 2 within the SM can be found in the appendix of Ref. [17]. The Higgs-mediated triangle loop diagram contributes only to the S z = 0 case. The corresponding form factor, for state H l attached to the triangle, is written as H j H i g g Figure 1: The Feynman diagrams contributing to Higgs boson pair-production in a model with an extended Higgs sector, but no additional particle content besides the SM.
where the scalar and pseudoscalar components, S g l and P g l , respectively, can be found in, e.g., Refs. [18,19]. In case of a single Higgs boson, as in the SM, the triangle coefficient in Eq. (2) is given as where λ hhh is the Higgs trilinear self-coupling. In multi-Higgs models like the ones we intend to explore here, the above coefficient is generalised to Here, λ H i H j H k are the Higgs trilinear couplings and D kl (ŝ), with k, l = 1, ..., N , are the entries of the Higgs propagator matrix. This modified C l allows the possibility of interference between two different Higgs intermediate states, induced by higher order quantum effects, as illustrated by Fig. 1 (left).
The main focus of this study is to investigate the above-mentioned quantum effects in the specific scenario with N = 3, which permits the resonant pair-production of the lightest Higgs state via the two, mutually interfering, heavier states. In this scenario, the (symmetric) propagator matrix is written as where the ImΠ kl (ŝ)'s are the absorptive parts of the Higgs self-energies, and m H k is the renormalised mass of the k-th Higgs boson. The explicit expressions for ImΠ kl (ŝ) can be found in the Appendix of Ref. [15]. In general, however, the off-diagonal absorptive terms in the propagator are assumed to be negligible, in which case the D(ŝ) becomes a diagonal matrix and C can, to a good approximation, be reduced to a sum over three terms containing BW propagators corresponding to each H l , as in Eq. (5).

Models and computational tools
Two new physics models that are consistent with the N = 3 scenario are the NMSSM and N2HDM.
In both these models, we identified the lightest of the three scalars, H 1 , with h obs , and analysed the impact of the interference between the heavier H 2 and H 3 on its pair-production at the √ s = 14 TeV LHC. It has previously been emphasised in literature [8,15,20,21] that these effects are more pronounced for large (combined) total widths of the intermediate Higgs states compared to the mass splitting between these. Using this as a guideline, we numerically scanned the parameter spaces of the two models to find their potentially relevant configurations. Below we discuss some details of the two models as well as of our calculation of the H 1 H 1 production cross section.

NMSSM
As a follow-up of our previous analyses, the first model that we investigate is the Z 3 -symmetric NMSSM. The Higgs potential in this model is written in terms of the two SU (2) L doublets H u and H d , with Y = ±1, and the singlet S as Here λ and κ are dimensionless Higgs trilinear couplings and A λ and A κ are their respective soft SUSY-breaking counterparts, m H d , m Hu and m S are the soft Higgs masses, while g 1 and g 2 are the U (1) Y and SU (2) L gauge coupling constants, respectively. The neutral components of the fields H d , H u and S are developed around their respective Vacuum Expectation Values (VEVs) v d , v u and v S , when EW symmetry is broken, as By taking the second derivative of V NMSSM , one then obtains the tree-level 3 × 3 neutral CP-even Higgs mass-squared matrix, M 2 H , in the (H dR , H uR , S R ) T basis. The orthogonal matrix R rotates these interaction eigenstates into the physical states as The matrix M 2 H thus gets diagonalised as with Higgs boson masses in the ascending order, i.e., m H 1 < m H 2 < m H 3 . Our current analysis pertains to the 'phenomenological' NMSSM, wherein all the free parameters, including the above Higgs sector ones, are input at the EW scale. Since variations in non-Higgs sector parameters are expected to have little impact on our particular phenomenological scenario, we fixed the soft squark masses as M Q 1,2,3 = M U 1,2,3 = M D 1,2,3 = 3 TeV, the slepton masses as M L 1,2,3 = M E 1,2,3 = 2 TeV, the soft gaugino masses as 2M 1 = M 2 = 1 3 M 3 = 1 TeV and the trilinear couplings of the charged sfermions as A 0 ≡ Aũ ,c,t = Ad ,s,b = Aẽ ,μ,τ . This resulted in A 0 , tan β (≡ vu v d ), µ eff (≡ λv s ), λ, κ, m P and m A as the complete set of inputs. The final two parameters are the bare masses of the two pseudoscalars, which are a trade-off for A λ and A κ using the minimisation conditions of the Higgs potential.
For numerically generating the Higgs boson mass spectra and Branching Ratios (BRs) corresponding to each set of values of the 7 model input parameters, each randomly selected from a broadly-defined range, we used the public code NMSSMTools-v5.5.2 [22, 23,24]. Each parameter space point was required to satisfy all the theoretical and experimental constraints defined in  NMSSMTools, which include limits from the Higgs searches at the Large Electron-Positron (LEP) collider, the TeVatron and the LHC, from the direct and indirect searches for neutralino Dark Matter (DM) and estimates of its relic abundance and from B-physics measurements. Output points satisfying all these constraints were further run through HiggsBounds-v5.7.0 [25,26,27,28,29] to test the Higgs sector observables of the model against the latest exclusion bounds from the LHC that might not (yet) have been included in NMSSTools itself.
In Fig. 2 we show the {m H 2 , m H 3 } plane for the resulting successful points, where one notices that the limits from the direct searches at the LHC rule out a mass below ∼ 405 GeV for the (predominantly doublet-like) H 3 , over the entire parameter space explored. Moreover, our initial scan with wide input ranges of the parameters yielded only one point (out of nearly two thousand violet points in the figure) with ∆m H ≡ m H 3 − m H 2 less than 5 GeV, lying just above the LHC exclusion bound for m H 3 . In order to find solutions with larger H 2 -H 3 mass-degeneracy, we therefore performed another scan of the narrowed-down parameter space region around the said point. Indeed, several points with ∆m H < 1 GeV were obtained with this secondary scan, which are plotted in blue colour in the figure. The coordinates of the point with the smallest ∆m H are which result in the following Higgs mass spectrum: The total widths of the three scalars yielded by this parameter space point are Γ H 1 = 4.76 MeV, Γ H 2 = 535.4 MeV and Γ H 3 = 24.78 MeV.

N2HDM
Since the scans for the NMSSM did not generate any points with a H 3 lighter than ∼ 405 GeV, we extended our analysis to the N2HDM also. In this model the Higgs boson masses are input parameters, as opposed to the NMSSM, where they are derived in terms of input variables, which allows greater freedom in the selection of the other free parameters relevant to the process under investigation. The N2HDM is obtained by adding a real singlet scalar field, S, to the (CP-conserving) 2HDM, and its Higgs potential reads where H u and H d are doublet fields similar to the NMSSM ones. This potential has a generic form and observes two symmetries , which is softly broken by the term containing m 2 12 , and ii) a spontaneously broken The charge assignments of the fermions under the Z 2 symmetry define the four types of the underlying (N)2HDM. Our adopted notation for the doublet Higgs fields is intended to indicate the Type-II N2HDM specifically, wherein the fermions have Z 2 charges such that the doublet H u couples only to the up-type quarks and H d to the down-type quarks and charged leptons. Upon EWSB, the two doublet fields are expanded around their respective VEVs according to Eq. (9), while the real singlet is expanded in this model as S = v S + S R . After minimisation of the potential and rotation of the scalar mass matrix, as in the NMSSM, the masses of the three physical CPeven Higgs states are obtained, with m H 1 < m H 2 < m H 3 . Besides these, the Higgs sector of the model also contains a CP-odd Higgs boson A. The Type-II N2HDM is thus essentially the non-Supersymmetic counterpart of the NMSSM, with fewer symmetries impinging on the properties of the CP-even Higgs sector, which makes it a more suitable fit for our comparative investigation than the other N2HDM types. For details of the Higgs sector of the N2HDM, we refer the reader to Ref. [30].
There are twelve free parameters in the potential in Eq. (14): λ 1,··· , 7,S , m 2 Hu , m 2 H d , m 2 S , m 2 12 . Relations between these parameters and the VEVs, arising from the minimisation conditions of the Higgs potential, allow us to trade m 2 and v S . Moreover, the eight quartic couplings can be traded for the physical masses, m H 1,2,3 , m H ± , m A , and the three independent parameters of the mixing matrix R in Eq. (10). These parameters, taken to be R 11 , R 12 and R 23 , can then further be replaced by the top-Yukawa and gauge couplings of the H 1 , defined in units of the corresponding couplings of the Higgs boson in the SM as Thus, for the purpose of this study, the following independent real parameters representing the N2HDM were randomly scanned in the given ranges using the public tool ScannerS-2 [31,32]: where sign(R 13 ) takes into account the sign ambiguity in the neutral scalar mixing. While the above ranges were mostly guided by existing literature on the model (see, e.g., Refs. [30,33,34]), the ones of g 2 H 1 V V and g 2 H 1 tt were based roughly on the current 2σ error-bar on the measurements of the corresponding couplings for the h obs at the LHC [35]. Several scans were performed for this model, in all of which we fixed v = 246 GeV and m H 1 = 125 GeV. The values of m H 2,3 , in contrast, were set to certain different values of interest in different scans (as will be explained in the next section). The purpose of the numerical scanning was to find configurations of the parameters in Eq. (16) that satisfied theoretical conditions such as unitarity and vacuum stability, and were at the same time consistent with precision EW and B-physics measurements. In addition to these checks performed internally by ScannerS, testing of the Higgs sector observables against the exclusion bounds from direct collider searches was also performed for each scanned point, by interfacing it with N2HDECAY [36] and HiggsBounds. Finally ScannerS was also interfaced with the program HiggsSignals-2 [37], which performs a χ 2 -fit of the h obs properties for a given model point against the LHC measurements, and rules it out if ∆χ 2 = χ 2 N2HDM − χ 2 SM > 6.18 (assuming a 2σ Gaussian error on the best-fit value).

Cross section calculation
For the output points from the scans, we proceeded to calculate the inclusive pp → H 1 H 1 cross section, using a FORTRAN code prepared in-house. For evaluating σ LO given in Eq. (1), the expressions corresponding to the triangle and box form factors were formulated following the public code HPAIR-v2.00 [17,38,39] and includes only the SM and the MSSM. The numerical computation of the next-to-LO (NLO) corrections to σ LO , which can be expressed as [38] ∆σ = ∆σ virt + ∆σ gg + ∆σ gq + ∆σq q , were also imported from HPAIR, since they are generic to all models. Besides catering to models beyond the MSSM, another significant way that our cross section calculator differs from HPAIR, which evaluates individual BW propagators for each intermediate Higgs boson in the triangle diagram, is in the incorporation of the full propagator matrix of Eq. (7). This allows us to estimate the magnitude of the effects resulting from the off-diagonal terms in the matrix, by including or neglecting these during the cross section computation for a given point by our code.
Since the input parameters as well as the particle contents, and hence the Higgs self-energy contributions, of the NMSSM and N2HDM are mutually rather different, we prepared a separate code for each of these models. In order to check the accuracy and consistency of our base code, we compared the pp → H 1 H 1 cross sections calculated in the MSSM limit of the NMSSM for a few test points with the ones obtained from HPAIR. We found the two sets of results to be in very good (within 1%) agreement. Note here that the higher order QCD corrections for this process have now been evaluated up to the next-to-next-to-next-to-LO [40,41,42,43,44,45] in the SM. We presume that these can be extended straightforwardly to the multi-Higgs models discussed here, and their overall impact would amount to a simple rescaling of our NLO calculations.

Analysis results
To quantify the magnitude of the triangle-box interference arising in the S z = 0 channel and, additionally, the full-propagator effects within the triangle diagram, we calculated the integrated cross sections corresponding to the following cases for each successful point from the scans for the two models: a) without triangle-box interference, with diagonal-only propagator matrix, b) with triangle-box interference, with diagonal-only propagator matrix, c) with triangle-box interference, with full propagator matrix. Below, these three cross sections will be referred to as σ a , σ b , and σ c , respectively.

The NMSSM
The top half of Fig. 3 (left) shows the cross sections σ a (top half) and σ b (bottom half) with diagonal-only propagator matrix, as functions of the H 3 mass. One sees a large negative impact of the triangle-box interference, reducing the cross section uniformly by ∼ 35 fb for all the points. We note here that, in models with Supersymmetry, the box and triangle diagrams in principle include loops from squarks also. Here we take the view that the squarks are always too heavy to contribute significantly to either of these production processes (recall that we fixed the soft squark masses to 3 TeV in our parameter space scans, to prevent the physical sparticle masses from conflicting with the direct search results from the LHC), and thus retain only the quark loops. A detailed study of the impact of the inclusion of squarks in the MSSM and the NMSSM (without the Higgs propagator interference effects) can be found in Refs. [46,47,48]. The small blue and red islands near the lowest allowed m H 3 and with overall largest cross sections in the top and bottom halves, respectively, are the points with ∆m H < 5 GeV obtained from the secondary scan.
In the numerical calculation of the propagator matrix, in contrast, the (one-loop) Higgs selfenergies due to all the relevant NMSSM particles were included. The right panel of Fig. 3, however, shows negligible impact of introducing the full propagator matrix. This figure, restricted only to the points with ∆m H < 5 GeV, shows the ratio of the cross sections corresponding to cases b and c, R σ ≡ σ b /σ c , against the ratio of the sum of H 2 and H 3 widths, Γ H , and ∆m H . Note that Γ H ranges between 535 MeV and 565 MeV for all the points, implying that when Γ H 2 , depicted by the colour map in the figure, reaches its maximum value, Γ H 3 is at its minimum, and vice versa. The fact that the lowest ∆m H obtained is 0.8 GeV, according to Eq. (13), implies that Γ H /∆m H is always smaller than 1 and hence the above mentioned condition of larger Γ H than ∆m H for a sizeable enhancement in the propagator effects is never met. Still, one can notice a small gradual increase in R σ , meaning an increasing negative effect of the full propagator, as Γ H rises with respect to ∆m H . This effect is more pronounced for points with H 2 and H 3 widths closer to each other in magnitude, as illustrated by the violet/red points in the top left quadrant of the figure. A larger gap between these two widths, in contrast, generally tends to slightly increase σ c compared to σ b (the points in the bottom left quadrant).
The overall smallness of R σ in the NMSSM can be attributed partly to the large squark and slepton masses, so that their contribution to the Higgs self-energies is diminished, and partly to the specific Yukawa and gauge coupling combinations of the H 2 and H 3 in the narrow parameter space region yielding large mass degeneracy between these. Fig. 4 shows R σ as a function of these coupling combinations. One notices in these figures that the colour-mapped points, which correspond to the parameter space region with ∆m H < 5 GeV, mark the boundaries of the (black) points from the extended scan. Thus, the search results from the LHC, besides directly constraining the mass of the H 3 to lie above ∼ 405 GeV, also restrict its top-Yukawa coupling to fairly small values, with either sign. The condition of mass degeneracy with H 3 then also dictates the signs and sizes of the H 2 couplings.
According to the left panel of Fig. 4, while the H 2 and H 3 top-Yukawa couplings can take up three different sign combinations in general, for points with ∆m H < 5 GeV, the sign of g H 2 tt is always negative, while that of g H 3 tt can be both negative or positive. However, only positive g H 3 tt values appear for large negative values of g H 2 tt . As the magnitude of the latter drops, that of the former increases, with R σ also rising slowly, until both reach equal values (with opposite signs). At that point, the sign of g H 3 tt flips to negative, giving the largest R σ according to the colour map. A further increase in its magnitude, however, along with a decrease in the size of g H 2 tt , leads to a lowering of R σ again. In short, largest (allowed) values of one of the two top-Yukawa couplings, whether positive or negative, coupled with the smallest value of the other, results in σ c > σ b and, as the two tend towards each other, σ c starts to lower towards σ b and eventually below it.
The   Higgs boson is always opposite to that of its top-Yukawa coupling, so that g H 2 bb is positive only, conversely to g H 2 tt . Furthermore, R σ shows a similar trend with the variation in the sizes of g H 2 bb and g H 3 bb as with the top-Yukawa couplings -the largest (allowed) value of one bottom-Yukawa coupling paired with the smallest value of the other yields σ c > σ b , while σ c ≤ σ b results from their comparable magnitudes. The dependence of R σ on the relative signs and magnitudes of g H 2 V V and g H 3 V V follows the behaviour of the top-Yukawa couplings exactly, as seen in the right panel of the Fig. 4. Their allowed values are, however, much smaller than even those of the top-Yukawa couplings, pointing towards the decoupling regime of the (N)MSSM. As for the remaining couplings of the H 2 and H 3 , even when the corresponding (s)particles have sufficiently low masses, including A 1 as well as χ 0 1,2 and χ ± 1 (which are higgsino-like and thus have masses ∼ µ eff ∼ 150 GeV, see Eq. (13)), their influence on R σ is too small to merit a discussion here.

The Type-II N2HDM
As indicated earlier, the Higgs boson masses are input parameters in the N2HDM, which allows us to investigate H 2 and H 3 with exactly equal masses, that can also be much lower than those obtained in the NMSSM. For a direct comparison with the NMSSM though, in our first scan for this model we set m H 2 = m H 3 = 410 GeV, and the σ a and σ b for the 100 successful points thus obtained are shown in blue and red, respectively, in the left panel of Fig. 5 against the width of H 3 . Contrary to the NMSSM, triangle-box interference does not reduce the cross section uniformly for all the points. While for most of the points σ b is smaller by a few tens of fb than σ a , the former is larger than the latter by upto 10 fb for a few points. This is owing to the wider ranges of magnitudes as well as sign combinations for the H 2 and H 3 couplings being available in this model, as will be explained later. Notice also that Γ H 2 can reach a few GeVs and, in fact, Γ H 3 can simultaneously be quite large, as illustrated by the horizontal axis of the right panel of the figure.
The vertical axis of this panel shows the impact of including the full propagator matrix in the cross section calculation.
For a number of points from the first scan with m H 2 = m H 3 = 410 GeV, seen in red in Fig. 5  (right), σ b is a few times larger than σ c , but for some points it gets reduced by upto 35%, implying a net positive contribution from the off-diagonal terms in the matrix. For two of the red points, though, R σ exceeds 100, meaning a two orders of magnitude reduction in σ b . (We point out here that these two points were omitted from the left panel to keep the y-axis scale visually interesting, but the corresponding cross sections will be provided below.) To assess the effect of reducing the mass degeneracy, the points from scans with ∆m H = 5 GeV and ∆m H = 10 GeV, with m H 3 still fixed to 410 GeV, are also plotted in this figure in orange and green, respectively. Evidently, a larger ∆m H results in smaller fluctuations in σ b , as the R σ value lies very close to 1 for all the 100 green points. The violet points in the figure correspond to the scan with m H 2 = m H 3 = 300 GeV. While in general R σ can deviate substantially from 1 for these points also, its maximum value does not exceed 3. The main reason for this is that the widths of H 2 and H 3 are always lower than 1 GeV in this case, unlike the m H 2 = m H 3 = 410 GeV case, owing to the fact that their masses lie below the tt production threshold. Lowering m H 2 and m H 3 even further to under the H 1 H 1 threshold expectedly results in a vanishing impact of the off-diagonal propagator matrix terms, as demonstrated by the cyan points in the figure, which are all clustered together below Γ H 200 MeV. For a detailed investigation, we selected six Benchmark Points (BPs) from our main scan with m H 2 = m H 3 = 410 GeV. The input parameters, the widths and couplings of H 2 and H 3 as well as the four cross sections corresponding to these points are given in Table 1. σ 2×2 in the table implies the cross section obtained by setting m H 3 → ∞, with all the other input parameters fixed to their exact values for a given BP, and is quoted for reference. BP1 and BP2 are the two points with the highest R σ in Fig. 5, for BP3 and BP4 the R σ value lies very close to 1 while BP5 and BP6 are chosen from amongst the points for which σ c is slightly enhanced compared to σ b .
In the 2HDM, and the N2HDM by extension, of the type-II, the B-physics measurement strongly constrain m H ± [49,30] and therefore the latter lies above 600 GeV for all the successful points from the scans, while tan β is also pushed to smaller values, as can be noted in the table. m A is then also restricted to values close to m H ± by the EW precision constraints. One feature distinguishing the points with the largest R σ (BP1 and BP2) from the rest of the BPs are the larger m H ± and g H 1 tt values and relatively small tan β. Such parameter configurations result in specific combinations of the couplings of H 2 and H 3 for BP1 and BP2, which in turn lead to very high R σ for these. For these two points, g H 2 tt and g H 3 tt are both positive and large while g H 2 bb , g H 3 bb , g H 2 V V and g H 3 V V are all negative. In the case of BP3, g H 2 tt and g H 3 tt have signs opposite to each other while g H 2 bb and g H 3 bb are both negative. We note here that, again in contrast with the NMSSM, g H 2 bb is negative for all the 100 N2HDM points for the m H 2 = m H 3 = 410 GeV scenario, and also that, for a majority of these points, three out of the four Yukawa couplings had the same signs. BP4 and BP5 are very similar points in that the two top-Yukawa couplings have signs that are opposite not only to each other but also to the signs of the corresponding bottom-Yukawa couplings. BP6 is the only point of its kind found in the scan, with g H 2 tt , g H 3 tt and g H 2 bb all having negative signs, and it therefore uniquely exhibits a constructive triangle-box interference as well as constructive propagator interference, so that σ a < σ b < σ c .
In Fig. 6 we show the four cross sections as functions of the most important H 2 couplings in this context, for all the six BPs. Columns 1, 2, 3 and 4 of the figure correspond to σ 2×2 , σ b , σ c and R σ , respectively, and rows 1, 2, 3, 4 and 5 to g H 2 tt , g H 2 bb , g H 2 V V , g H 2 H 1 H 1 and g H 2 ,A,Z , respectively. The plotted ranges of these couplings are indicative of those observed across all the 100 points obtained for the m H 2 = m H 3 = 410 GeV scenario. Once again, for each BP in a given panel, all the remaining couplings are fixed to the values described in Table 1. The red lines in the figure correspond to BP1, green to BP2, olive to BP3, violet to BP4, blue to BP5 and orange to BP6. The point on a line in a given panel marks the actual value of the plotted coupling for that BP. The horizontal black lines in the panels in columns 1, 2 and 3, indicate the current experimental limit on h obs pair-production cross section [50], which we approximate to be 1 pb for m H 2 = m H 3 = 410 GeV considered here. Note also that, since only the product of the corresponding H 2 and H 3 couplings enters the iImΠ 23 (ŝ) element of the Higgs propagator matrix, the behaviour of σ c with varying H 3 couplings should by and large mimic that with varying H 2 couplings.
In the top row of the figure, one sees that the presence of an additional Higgs boson degenerate in mass with H 2 eliminates the peaks appearing at specific values of g H 2 tt in σ 2×2 , so that the variations in σ b in column two are relatively smooth. As expected, σ b shows a very similar behaviour for BP1 and BP2, reaching values even higher than the true ones for slightly different positive g H 2 tt (recall that g H 3 tt is also positive for these two points). Thus g H 2 tt 0.1 would be ruled out by the LHC h obs h obs -production limits. When g H 2 tt switches sign to negative, σ b drops to much smaller values. The introduction of the full propagator matrix then largely mitigates the very strong dependence of the cross section on positive g H 3 tt , as seen in the third column, bringing it down to values consistent with experimental bounds. And since σ c shows little variation with g H 2 tt , the shapes of the red and green lines in the top row of column 4 (and also in the rows beneath it) are very similar to those in column 2, with R σ reaching about 900 for the BP1.
Cross sections for BP3 and BP4, both of which have g H 3 tt with mutually opposite signs but very similar magnitudes, show similar trends to each other with the variations in g H 2 tt across the four top panels. For these two points, the peaks in σ 2×2 are the tallest, while σ b and even σ c violates the experimental bound for large negative g H 2 tt . BP4 and BP5, likewise mimic each other's behaviour for positive g H 2 tt , but since BP4 has a negative g H 3 tt larger in magnitude than that in BP5, its dependence on negative g H 2 tt is much more pronounced for both σ b and σ c . The second row of the figure shows negligible dependence of σ b on g H 2 bb for all the BPs expect 1 and 2 and, conversely, the least variation in σ c for these two BPs. This is due to the fact that for these points g H 2 bb and g H 3 bb both have negative signs, opposite to signs of the two top-Yukawa couplings which have a much more dominant effect.
The third row in the figure illustrates that the g H 2 V V coupling plays a role as crucial as the top-Yukawa couplings. Similarly to the NMSSM, these couplings originally have generally quite small magnitudes, as a consequence of the very SM-like properties of the H 1 . For this coupling, the two peaks appearing in σ 2×2 are replaced by a very tall narrow peak in σ b , especially for BP1 and BP2, close to g H 2 V V = 0. The introduction of the full propagator matrix smoothens the behaviour of the cross section with respect to this coupling, except for BP3, for which σ c hits a sharp peak when the negative g H 2 V V has the same magnitude as the positive g H 3 V V . The shapes of all the lines are therefore largely dictated by the interplay between the signs and sizes of the top-Yukawa and gauge couplings of H 2 and H 3 . In row 4 is depicted the dependence of the cross sections on g H 2 H 1 H 1 , which is the only coupling of significance other than the ones discussed above. Here, σ 2×2 shows a sharp dip at the zero of this coupling, since it also enters the H 2 → H 1 H 1 decay besides the self-energies. This sharp dip shifts away from zero for σ b , according to the relative sign of the diagonal H 3 contributions to the propagator. It returns to zero when the off-diagonal terms are also turned on. Around the minimum, σ c shows a fairly symmmetric behaviour in both signs of g H 2 H 1 H 1 , as do σ b and σ b . Unlike these two cross sections, however, σ c increases rather smoothly.
Finally, the bottom row of the figure shows a negligible dependence of each of the cross sections on the g H 2 AZ couplings, even though the AZ pair is the lightest of all the non-SM pairs contributing to the Higgs self-energies, besides the H + W − pair. Evidently then, the plots for any of the remaining couplings noted in Table 1 look exactly the same as the corresponding ones for g H 2 AZ in this row and can thus be safely ignored.   Table 1. The point on a line marks the actual value of the plotted coupling for the corresponding BP. See text for more details.

Conclusions
The commonly adopted approach of calculating the cross section for a given 2 → 2 process by factorising it into its production and decay parts cannot account for possible quantum interference amongst the propagators of several mass-degenerate states. This is because such an approach by construction assumes narrow widths for the resonant mediators. In some previous papers we explored such interference effects in the case of the gluon-fusion production of certain SM final states, via two highly mass-degenerate Higgs mediators. The mass-splitting between the two intermediate Higgs states being comparable to or smaller than the sum of their widths is a precondition for such effects to be sizeable, in both the integral and the differential cross section. The reason for their onset is that the imaginary off-diagonal elements of the Higgs propagator matrix become comparable to the imaginary parts of the diagonal elements, irrespectively of whether they are taken into account coherently or incoherently. These studies were performed within the illustrative theoretical framework of the NMSSM.
In this article, we have extended our investigation to the pair-production of the lightest of the three NMSSM neutral Higgs scalars, H 1 , at the 14 TeV LHC, taking into account the contributions of the triangle as well as the box diagram that this process proceeds through. We have investigated the impact of not only the interference between these two topologies, but also of the aforementioned propagator interference between H 2 and H 3 within the triangle topology on the cross section for H 1 H 1 production. Furthermore, since the lowest H 3 mass attainable in the NMSSM is 405 GeV, owing to the bounds from the LHC searches, we have also included the N2HDM in our analysis. In this framework, since the physical Higgs boson masses are input parameters, we could choose any desired (unique) value for m H 2 and m H 3 , which allowed us to study also the scenario where they contribute non-resonantly to the triangle topology for the studied process.
In the case of the NMSSM, we have found the effects of the inclusion of the full Higgs propagator matrix in the triangle topology to be similar in size to those established in our previous studies. However, given the various constraints imposed, since the minimal mass-splitting between H 2 and H 3 is obtained in a very narrow region of the parameter space where, however, the sum of their widths never exceeds it, the effect is largely subdued. In this region the box diagram and the triangle diagram with (an off-shell) H 1 in the propagator contribute much more dominantly to the H 1 pair-production process. The narrowness of this region also means that a nearly constant negative interference is always observed between the two topologies.
In the N2HDM, on the other hand, we have seen that the propagator interference effects can modify the cross section by more than two orders of magnitude. Of particular importance is the observation that these effects tend to 'regulate' the behavior of the total H 1 H 1 cross section, smoothing the peaks that appear in it for certain specific values of the H 2 and H 3 couplings, and generally bringing it down to values consistent with the current LHC limits. Moreover, herein the interference between the box and triangle topologies can be positive or negative, as a consequence of the relatively wider ranges of the magnitudes and sign combinations of the Higgs boson Yukawa couplings. Clearly, such a disparity between the results obtained for this model and those for the NMSSM is due to the fact that supersymmetry imposes strong limitations on the masses and couplings of the heavy Higgs states. In the N2HDM, these quantities are essentially free parameters. But even in this model, when m H 2 and m H 3 lie below the H 1 H 1 production threshold, the propagator interference effects tend to vanish. In this case the one-loop two-point functions corresponding to the off-diagonal elements in the Higgs propagator matrix are too small to be able to overcome the kinematic suppression.
Finally, we emphasise that we have reached the above conclusions on the basis of a detailed analysis at the level of the total cross section. As for their phenomenological relevance, we believe that the LHC may develop sensitivity to all such interesting dynamics already at its upcoming Run 3, at least in the N2HDM. Hence, for the purpose of aiding experimental efforts in establishing all the effects studied here, we have proposed some BPs, compliant with the latest theoretical and experimental constraints, that are amenable to dedicated probes by the ATLAS and CMS collaborations.