Quantum interference among heavy NMSSM Higgs bosons

In the Next-to-Minimal Supersymmetric Standard Model (NMSSM), it is possible to have strong mass degeneracies between the new singlet-like scalar and the heavy doublet-like scalar, as well as between the singlet-like and doublet-like pseudoscalar Higgs states. When the difference in the masses of such states is comparable with the sum of their widths, the quantum mechanical interference between their propagators can become significant. We study these effects by taking into account the full Higgs boson propagator matrix in the calculation of the production process of $\tau^+\tau^-$ pairs in gluon fusion at the Large Hadron Collider (LHC). We find that, while these interference effects are sizeable, they are not resolvable in terms of the distributions of differential cross sections, owing to the poor detector resolution of the $\tau^+\tau^-$ invariant mass. They are, however, identifiable via the inclusive cross sections, which are subject to significant variations with respect to the standard approaches, wherein the propagating Higgs bosons are treated independently from one another. We quantify these effects for several representative benchmark points, extracted from a large set of points, obtained by numerical scanning of the NMSSM parameter space, that satisfy the most important experimental constraints currently available.


I. INTRODUCTION
Supersymmetry (SUSY) predicts the existence of at least two Higgs doublets of opposite hypercharge, which are necessary to cancel chiral anomalies, unlike the standard model (SM), where only one Higgs doublet is sufficient. In the minimal supersymmetric standard model (MSSM), after the breaking of the electroweak (EW) symmetry, these two doublets result in a total of three physical neutral Higgs states, two scalars (h and H, with m h < m H ) and a pseudoscalar (A), as well as a charged pair (H AE ). In principle, either of the h and H could be the Higgs boson, h obs , observed at the LHC [1,2]. However, the model's parameter space region where H has a mass near 125 GeV and almost SM-like couplings to EW gauge bosons is very tightly constrained by experimental data [3]. Alternatively, the condition on h to be a candidate for h obs , while a phenomenologically much more favored scenario, pushes m H and m A upwards into the so-called decoupling regime (see, e.g., [4]), where they are nearly identical.
Thus, if a bump appears in the LHC data near a certain (large) invariant mass of a fermion-antifermion or photon pair above the SM expectation, it could very well be due to these two states combined, unless their mass difference is large enough to enable the experiment to resolve them individually. A formidable difference in the statistical significances of the peaks observed in the fermionic channels versus the W þ W − and ZZ channels could also be interpreted as a hint of the fH; Ag mass-degeneracy. However, if the Higgs sector of the MSSM is CP-violating, the two heavy mass-degenerate states could mutually interfere quantum mechanically. This interference could significantly alter the cross section expected for a given SM final state that these Higgs bosons decay into [5][6][7].
In the NMSSM [8], which originally aimed at addressing the μ-problem of the MSSM [9,10] by invoking an additional SUð2Þ L singlet Higgs field, there exists an extra scalar, h s , and an extra pseudoscalar, a s , in the Higgs sector. The masses of these two new Higgs bosons are essentially free parameters of the model and either of them could very well lie near 125 GeV, along with the MSSM-like h. The scenario with m h s near m h is in fact well motivated by naturalness considerations, since the doublet-singlet mixing can enhance the tree-level mass of h appreciably [11][12][13]. The potentially strong effects of interference in the scenario where these Higgs bosons are too close in mass for the peaks to be separately identifiable, given the current diphoton mass resolution at the LHC, have therefore been recently analysed in [14]. In this work, we consider the alternative scenarios wherein one (or possibly even both) of h s and a s could be mass-degenerate with the heavier H and A instead of h. Our main objective here is to investigate the phenomenological implications of such a scenario, which is a viable one in general extended Higgs sectors, without delving into its theoretical motivations in the specific model considered.
For our analysis, we first perform numerical scans of the NMSSM parameter space to identify regions where m h s ≈ m H ð≈m A Þ or m a s ≈ m A ð≈m H Þ. We then quantify the impact of quantum interference on the process where nearly mass-degenerate heavy Higgs bosons are produced in gluon fusion at the 14 TeV LHC and decay into τ þ τ − pairs. This is done by comparing the cross section obtained by including the full Higgs propagator matrix in the expression for amplitude and the one obtained by assuming two (or more) individual Breit-Wigner (BW) Higgs propagators with nearly identical poles. We restrict ourselves to a CP-conserving NMSSM Higgs sector here, so that the scalar and pseudoscalar interaction eigenstates do not mix. Hence, the off-diagonal elements in the Higgs propagator matrix contribute to the interference affects only when the two mass-degenerate states have the same CP-identities. Furthermore, employing the diphoton final state, as in our previous study pertaining to 125 GeV Higgs bosons in a similar context [14], is not feasible here. This is because the cross section, which is already phase-space suppressed, gets further depleted considerably due to the extremely small partial decay widths-and consequently branching ratios (BRs)-of the heavy Higgs states into photon pairs.
During the numerical scanning, each randomly generated parameter space point is tested against the most crucial and recent experimental constraints, including those from the LHC, from B-physics measurements and from dark matter (DM) searches. Among the ones passing all the constraints, we then identify a few benchmark points (BPs) and carry out the cross section comparison noted above for them, using a Monte Carlo (MC) integration code developed in house. We find that, as concluded by the previous studies for the MSSM [7,15], the interference is (almost) always destructive and can result in a sizeable reduction of the cross section for the considered process. However, the heavy Higgs boson masses for the successful points collected in our scans are always quite large (>800 GeV). As a result, the effects of interference on the shape of the distribution of the differential cross section may not be detectable even with an integrated luminosity of 6000 fb −1 at the LHC, owing to the low τ þ τ − invariant mass resolution.
The article is organised as follows. In the next section we will briefly revisit the Higgs sector of the NMSSM and the analytical expression for the cross section that includes the full Higgs propagator matrix. In Sec. III we will discuss in some detail our methodology for the parameter space scanning. We will present the results of our comparative numerical analysis of event rates in Sec. IV and our conclusions in Sec. V.

II. THE τ + τ − SIGNAL FROM HEAVY NMSSM HIGGS BOSONS AT THE LHC
A. The NMSSM Higgs sector The Higgs potential of the Z 3 -symmetric NMSSM is written in terms of the two SUð2Þ L doublets H u and H d , with Y ¼ AE1, 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 H u and m S are the soft Higgs masses, and g 1 and g 2 are the Uð1Þ Y and SUð2Þ L gauge coupling constants, respectively, with g 2 ¼ g 2 1 þg 2 2 2 . By taking the second derivative of V 0 after developing the fields H d , H u and S around their respective vacuum expectation values (VEVs), v d , v u and v s , as one obtains the tree-level 5 × 5 neutral Higgs mass-squared matrix in the H T ¼ ðH dR ; H uR ; S R ; H I ; S I Þ basis, from which the Goldstone state, G, has been rotated away. In the CP-conserving limit, where all the Higgs sector coupling parameters are real, after the inclusion of the higher order corrections, this mass matrix can be expressed in the block-diagonal form where the 3 × 3 submatrix M 2 S corresponds to the CP-even states and the 2 × 2 matrix M 2 P to the CP-odd states, while all the CP-mixing terms vanish. The CP-even and CP-odd Higgs mass eigenstates can be obtained from the interaction states through the rotations ðh s ; h; HÞ T ¼ R H ðH dR ; H uR ; S R Þ T and respectively, where R H and R A are orthogonal matrices. Note that we have chosen the above notation for the physical Higgs states in order to distinguish between the MSSM-like and the NMSSM-specific ones. Throughout our analysis below, we identify the lighter MSSM-like (or doublet-like) scalar h with the h obs , so that the H is always heavier than 125 GeV, and we disregard the alternative possibility of H playing the role of the h obs . The mass of the singlet-like h s (a s ) can be smaller or larger than m h and/or m H (m A ).

B. Gluon fusion production of Higgs bosons
The differential cross section for the process pp → H → τ þ τ − (with H collectively denoting the five neutral Higgs bosons of the model, and assuming vanishing contribution from production modes other than gluon fusion) can be written as where gðx 1 Þ and gðx 2 Þ are the parton distribution functions (PDFs) of the two incoming gluons having squared centerof-mass (CM) energyŝ, and τ ¼ŝ=s, with s being the total CM energy of the pp system. In the limit of negligible nonfactorizable loop contributions, the amplitude-squared in the above equation can be cast into the most general form where λ ¼ AE1 and σ ¼ AE1 are the helicities of the incoming gluons and outgoing τ AE , respectively. The matrix element for the production part, when the Higgs sector is CP-conserving, is given as where the scalar and pseudoscalar form factors, S g i and P g i , respectively, can be found in, e.g., [16,17]. The amplitude for the decay part is similarly given as . In principle, for calculating the differential cross section given in Eq. (5), one should include in the expression for A 2 gg→τ þ τ − the full NMSSM Higgs boson propagator matrix, written as where ImΠ ij ðŝÞ are the absorptive parts of the Higgs selfenergies, for i, j ¼ 1-5, and M ii ≡ŝ − m 2 H i þ iImΠ ii ðŝÞ, with m H i being the renormalised mass of the ith of the five neutral Higgs bosons. The explicit expressions for ImΠ ij ðŝÞ can be found in the Appendix of [14]. The total cross section can then be obtained by integrating Eq. (5) over ffiffi f s p . In general, however, the off-diagonal absorptive terms in the propagator are assumed to be negligible compared to all m H i . In that case, the amplitudes due to all H i can be summed incoherently as In fact, if a given Higgs boson's width, Γ H i , defined through iImΠ ii ≡ im H i Γ i from the corresponding diagonal term M ii in Eq. (9), is additionally much smaller than m H i , one can further apply the narrow width approximation (NWA), to the amplitude above. This approximation implies onshell production of each H i , so that one can simply factorise its contribution to the total cross section as In the following, we define the total cross section for τ þ τ − production calculated using the above NWA as for all mass-degenerate H n . This is the most commonly adopted approach, and we aim to examine to what extent this cross section might differ from the one obtained by the incoherent summing of amplitudes with individual BW propagators, as in Eq. (10), which we refer to as σ BW . We then assess the further impact of the interference effects, by calculating the cross section, σ Int , calculated using the complete amplitude expression given in Eq. (6).

III. PARAMETER SPACE: SCANS, CONSTRAINTS, AND FAVORED REGIONS
In order to find solutions with nearly mass-degenerate fh s ; Hg and/or fa s ; Ag pairs that also satisfy experimental constraints from a variety of sources, we first performed numerical scanning of the NMSSM parameter space. As the model contains a large number of free parameters at the EW scale, we fixed the soft masses of the sfermions as M Q 1;2;3 ¼ M U 1;2;3 ¼ M D 1;2;3 ¼ 3 TeV and M L 1;2;3 ¼ M E 1;2;3 ¼ 2 TeV, and of the gauginos as 2M 1 ¼ M 2 ¼ 1 3 M 3 ¼ 1 TeV, since variations in these parameters are not expected to have a significant influence on the overall findings of our particular case study. Note also that the parameters A λ and A κ appearing in the potential in Eq. (1) can be traded for the pseudoscalar masses m P ð∼m a s Þ and m A as inputs, plus we assume universal trilinear couplings of the charged sfermions: A 0 ≡ Aũ ;c;t ¼ Ad ;s;b ¼ Aẽ ;μ;τ . The Higgs mass spectra and BRs for each randomly generated set of input parameters, over the wide initial scanned ranges of which are given in the second column of Table I, were calculated using the public code NMSSMTools-v5.3.0 [18,19].
As for the experimental constraints, a scanned input point was rejected if it did not predict m h lying within the 123-127 GeV bracket (thus allowing for up to AE2 GeV uncertainty in the theoretical prediction of its mass, given the experimental measurement of 125.09 AE 0.32 GeV [20]). In addition, we required the neutralino DM relic abundance, calculated by NMSSMTools via an interface to MicrOmegas [21,22], to satisfy Ω˜χ0 1 h 2 ≤ 0.131. This difference of the upper limit enforced on Ω˜χ0 1 h 2 from the actual PLANCK measurement of Ω DM ¼ 0.119 [23] is to accommodate up to a þ10% possible error in its theoretical evaluation. A point was also discarded during the scan if the spin-independentχ 0 1 -proton scattering cross section, σ SI p , did not satisfy the 95% confidence level (CL) limits from the XENON1T direct detection experiment [24]. The theoretical estimate of this cross section is also written out in the NMSSMTools output, as are those of the B-physics observables. The points collected in the scan were further filtered by the requirement on the most constraining of these observables to lie within 2σ of their latest measurements, which read (i) BRðB → X s γÞ × 10 4 ¼ 3.32 AE 0.15 [25], [26]. The successful points were then run through HiggsBounds-v4.3.1 [27] to test each of the additional NMSSM Higgs bosons against the exclusion bounds from TABLE I. Wide and narrowed-down (for a given scenario) scanned ranges of the seven NMSSM parameters considered free in this study.

Parameter
Initial wide scanned range Narrow range for scenario 1 Narrow range for scenario 2 with LEP, TeVatron and LHC as an added precaution, since this is done by NMSSMTools itself also. These points were further subjected to the 95% CL exclusion limits from the combined analysis ofχ 0 1χ 0 1 → bb annihilation in dwarf spheroidal galaxies performed by the Fermi-LAT and MAGIC collaborations [28], which are currently the strongest of the DM indirect detection bounds. Next, we calculated the theoretical predictions of the signal strengths, μ X , of h in the SM decay channels, X ¼ γγ, ZZ Ã , W þ W −Ã , τ þ τ − , bb, for each point using the public program HiggsSignals-v1.4.0 [29]. Since the discovery of h obs , the CMS and ATLAS collaborations have frequently updated the measurements of μ X independently from each other, and have also released their combined results based on the ffiffi ffi s p ¼ 7 and 8 TeV data for each channel in [30]. However, no such combined analysis for the ffiffi ffi s p ¼ 13 TeV data has been published so far and, while the measurements by the two groups have generally been in agreement with each other, there are also non-negligible differences in at least one of the channels (see, e.g., [31,32]). Given that these results have increasingly favored the SM predictions, instead of choosing one of the two results over the other or performing fits to both, which would be beyond the scope of this study, we simply enforced μ X ¼ 1 AE 0.34 for each X on all points, i.e., the theoretical signal rates were required to lie within 1σ of the SM value.
Finally, it has previously been established in the literature [15,33] that the interference effects between two Higgs states grow as the mass-splitting between them drops off compared to the sum of their widths. We therefore defined where i ¼ 1 implies the lighter and i ¼ 2 the heavier of the two nearly mass-degenerate scalars (X ¼ H) or pseudoscalars (X ¼ A), and retained only the points for which Λ X 1 > 1 or Λ X 2 > 1. These points were then split into two categories, (i) scenario 1: m h s ≈ m H , (ii) scenario 2: m a s ≈ m A . Note, however, that for a vast majority of the points left after applying the above filters, the doubletlike H and A were found to be highly mass-degenerate and thus lying in the decoupling regime of the MSSM. This is understandable, in particular in scenario 1, since the requirement for h s to have a mass close to that of H forces h to be almost entirely doubletlike, with maximal tree-level mass and SM-like couplings. Furthermore, for all the points the lightest neutralino turned out to be a pure higgsino, with the difference in its mass from that of the next-to-lightest neutralino always lying within 10-13 GeV. A recent ATLAS search [34] has ruled out such a DM for a mass up to ∼145 GeV (which effectively translates into μ eff ≥ 145 GeV). We therefore removed all the points with m˜χ0 1 < 145 GeV from the scanned set. We point out here that the limits from a CMS search [35] which could also be of relevance here have already been incorporated in NMSSMTools-v5.3.0, so that each scanned point was intrinsically tested against them.
The points in scenario 2 can be further divided into two distinct sets, one with h s lighter than h and the other with h s heavier than h. Columns 3-5 of Table I show the ranges of the input parameters covered by the points corresponding to the two scenarios. In order to find more solutions with a possibly enhanced mass-degeneracy in a given scenario, we performed secondary scans of each of these narroweddown parameters ranges. All the constraints and conditions noted above were reapplied to the points collected in these scans and the successful ones were combined with the corresponding set from the initial scan.
In the top left panel of Fig. 1, we show Λ H i for the final set of points belonging to scenario 1 (where H 1 , and likewise H 2 , can be either one of h s or H). We see that Γ H 1 can be up to 8 times larger than Δm H , as illustrated by the color map. According to the top right panel, the maximum width of H 1 for these points is about 4 GeV (while that of H 2 , not shown here, is much smaller, as can be deduced from the overall smaller Λ H 2 in the left panel). We point out again that Λ H 2 > 1 whenever Λ H 1 < 1, and vice versa. In the bottom left panel of Fig. 1, which similarly shows Λ A i for scenario 2, one notices Γ A 1 being as much as 50 times larger than Δm A . In this figure (and in the subsequent figures for this scenario) circles correspond to the subset of points for which m h s > m h and triangles to those with m h s < m h . Hence, while in the former case Λ A 2 never exceeds Λ A 1 , there are a few of the latter points for which both Λ A 1 , Λ A 1 ∼ 40. The bottom right panel illustrates that in the m h s > m h case large values of Λ A 1 are a consequence of large, ∼10 GeV, values of the A 1 width. In the m h s < m h case, in contrast, Γ A 1 stays low generally (between 2-4 GeV), and large Λ A 1 results mainly from the fact that Δm A can reach values lower than in the m h s > m h case, according to the colour map in the bottom left panel. The points encircled in green in Fig. 1 are the BPs we identified for our cross section analysis, to be explained in the next section.
We now briefly discuss the parameter combinations that lead to strong mass-degeneracies between Higgs bosons in the two scenarios. Figure 2 corresponds to scenario 1, and shows that a smaller Δm H favours lower values of both λ and κ but larger values of tan β (within the considered range). The reason is that the mass-degeneracy condition in this scenario effectively implies the decoupling of h s from h too. The smaller values of λ essential for that in turn put a virtual cut-off on μ eff ≡ λv s of about 250 GeV, below which it can vary relatively freely. Scenario 2 shows a much stronger dependence on the correlations between the input parameters, as seen in Fig. 3.
In particular, the parameter space regions with smaller Δm A are rather distinct for the m h s < m h and m h s > m h cases. Recall that the m h s < m h case forces the model into the "natural" limit, so that both h and h s are states with a large singlet-doublet mixing, made possible by larger λ and κ, and smaller tan β, as noted for the triangles in the various panels of the figure. Such values of these parameters, with some fine tuning, allow μ eff to reach up to ∼400 GeV. Note, however, that in the m h s > m h case also it is possible for h s to be a highly mixed state lying very close in mass to the h, and for such points the parameter values can be similar to the ones in the m h s < m h case. Generally though, the parameter combinations in the m h s > m h case follow a somewhat similar trend to that seen for scenario 1, with relatively small (large) λ and κ (tan β), and μ eff restricted to below 250 GeV.

IV. CROSS SECTION ANALYSIS
As noted earlier, we selected some BPs from the two scenarios to study the impact of quantum interference between the Higgs bosons on the cross section for the pp → τ þ τ − process at the LHC with ffiffi ffi s p ¼ 14 TeV. The calculation of the cross sections σ BW and σ Int defined in Sec. II B was carried out using a FORTRAN program developed in-house, which performs Higgs propagator matrix inversion and numerical integration through interfaces with the LAPACK package [36] and a locally modified version of the VEGAS routine [37], respectively. The Higgs boson masses and couplings for a given NMSSM parameter space point, which the program takes as inputs for the calculation of the propagator matrix as well as the form factors for Higgs boson production via gluon fusion, were written out by a suitably modified version of NMSSMTools. For comparison, we also calculated σ H 1 …H n by multiplying the σðgg → H i Þ obtained from the public code SusHi v1.6.0 [38] with BRðH i → τ þ τ − Þ given by NMSSMTools. We note here that we required SusHi to evaluate these σðgg → H i Þ at the next-to-next-toleading order (NNLO) in QCD, while our program calculates the production process only at the leading order (LO). 1 For a more accurate comparison, we therefore estimated the k NNLO factor by obtaining also the LO σðgg → H i Þ from SusHi for each H i in a given BP and multiplied our LO amplitude-squared with it before summing over all H i 's to obtain the total σ BW . In the case of σ Int though, since the contribution due to every H i is not computed separately, we simply multiply the total amplitude-squared with an average of the k NNLO factors for the (three) Higgs bosons near a given τ þ τ − invariant mass (or, equivalently, ffiffi f s p ). For all the cross sections, the CT10 [39] PDF sets were used for gluons and the default values of the SM input parameters inside NMSSMTools were retained, except for the t-quark mass, which was fixed to 172.5 GeV.
The values of the input parameters, the masses and widths of the Higgs bosons, as well as the total cross sections calculated using each of the approaches explained above are given in Table II for all the selected BPs. BP1-BP4 correspond to scenario 1, and the first one of these has been selected such that m h s and m H lie relatively close to their lowest observed values, in order to minimise phase-space suppression of the cross section, while also requiring Λ H 1;2 > 1. However, these masses are still only slightly smaller than 900 GeV, and thus, while the σ H 1 …H n for this BP is the highest among all, it nevertheless stays below 1 fb. Also noted down in the table is the quantity Δσ BW ≡ σ BW −σ H…H n σ H…H n × 100 and its value for BP1 indicates that the total cross section assuming BW propagators is enhanced by 16% over the one obtained using the NWA. Δσ Int ≡ σ Int −σ BW σ BW × 100 in the last row of the table similarly quantifies the impact of including the full propagator matrix in the amplitude. We note that σ Int is 11% smaller than σ BW , implying that the interference between the h s and H propagators is destructive, reducing the total cross FIG. 2. Distributions of the input parameters, showing the correlations among them that lead to a strong mass-degeneracy between h s and H, for the points corresponding to scenario 1 obtained from the numerical scans. 1 Including higher order (HO) corrections for the production process is a highly involved task, which would be beyond the scope of this work. Their inclusion is, however, tantamount to rescaling the cross section for a given H i by a "k factor," which is defined as k HO ≡ σ HO =σ LO , with HO implying the perturbative order at which the cross section is to be evaluated. section somewhat, although it is still larger than σ h s HA (i.e., σ H…H n , defined in Eq. (12), calculated assuming the NWA for H 1 ; …; H n ¼ h s ; H; A).
For all the points with m h s and m H similar to those for BP1, Λ H 1;2 stays very close to 1. Larger values of these ratios are typically obtained for higher masses, and BP2 was selected such that it has the former slightly increased, while the latter are still below 900 GeV. Even this slight increase in mass reduces each cross section by a few percent compared to the BP1 case, while enhancing the interference effects only marginally. The interference for this BP, and in fact also for the other two BPs for scenario 1, is again destructive. The selection criteria for BP3 were maximal Λ H 1 only, while for BP4 maximal Λ H 1;2 , with their respective m h s and m H being very similar. This effectively implies that for BP3 one has Γ H 2 ð¼ Γ h s Þ considerably smaller than Γ H 1 ð¼ Γ H Þ, while for BP4 both these widths are much closer in size. This difference leads to a clear rise in the magnitude of interference for BP4, which leads to a Δσ Int of 19%, compared to the BP3, for which the increase in the cross section is 11%, identically to BPs 1 and 2. Δσ BW for BPs 3 and 4 are also considerably smaller than for BP1, despite a large hike in m h s and m H , as is the respective absolute value of each type of the cross section itself. Note that μ eff for both these high-mass BPs lies just above the effective exclusion limit from ATLAS discussed earlier, and could thus be within the reach of subsequent analyses at the LHC in the near future. BRðH i → τ þ τ − Þ, with H i ¼ h s , H, A, for all the above BPs is ∼10%.
BP5 belongs to the m h s > m h case and BP6 to the m h s < m h case of scenario 2. As with BP3 and BP4, for BP5 also one finds that Γ a s is somewhat smaller than Γ A , while for BP6 the two widths are much closer in magnitude. Furthermore, σ BW shows an enhancement of about 22% over σ Ha s A for BP5, but σ Int is exactly the same as σ BW , implying that interference is nonexistent between the a s and A propagators. This is also the case for BP6, despite Λ a s and Λ A both being very large due to a very strong massdegeneracy between a s and A. σ h s HA for BP6 is sizeably smaller than for BP5, which is a consequence of the relatively small BRðA → τ þ τ − Þ, even though the respective pseudoscalar masses are only about 4 GeV larger. Notice that for both of these BPs, even though Γ a s are of order 1 GeV, the τ þ τ − partial decay widths, and hence the BRða s → τ þ τ − Þ, are extremely small, resulting in the vanishing of the interference effects.
Next, we also calculated the differential σ BW and σ Int with respect to the ffiffi f s p , in order to examine if the interference effects can lead to visible differences between the shapes of their distributions and may therefore be possible to probe at the LHC. For these we employ a binning template closely replicating the one used by the ATLAS collaboration in the searches for heavy resonances in the τ þ τ − decay channel [40], based on an expected detector mass resolution of ∼15% − 20% of M τ þ τ − . Our template thus assumes bins of width 50 GeV, 100 GeV and 150 GeV for ffiffi f s p ¼ 0-500 GeV, 500-800 GeV and 800-1400 GeV, respectively. The differential distributions are shown in Fig. 4 for all the six BPs, with the red lines corresponding to σ BW and the blue ones to σ Int . Since the mass-degenerate Higgs bosons in both the scenarios investigated are always heavier than 800 GeV and the remaining two Higgs bosons typically much lighter, these distributions have a lower cutoff at ffiffi f s p ¼ 500 GeV. In fact, for these distributions we included only H i with i ¼ 3-5 in the calculation of the two differential cross sections.
Importantly though, the two differential distributions do not reveal much beyond what can be inferred from the total cross sections, owing to the poor M τ þ τ − resolution at the LHC. We noticed earlier that Δm H in scenario 1 never exceeds 4 GeV, and even this maximum value is about 1=38th of the adopted widths of the bins around ffiffi f s p ¼ 1 TeV, despite ∼0.15 × M τ þ τ − itself being on the conservative side of the expected detector resolution. Thus, the red boxes in each of the first four panels of Fig. 4 rise higher than the blue ones (recall that the interference is always destructive) only in the bins within which the three mass-degenerate Higgs bosons lie. In the remaining bins, the boxes in the two histograms almost entirely overlap. Evidently, for BP5 and BP6 from scenario 2, a complete TABLE II. Masses and decay widths of the Higgs bosons, and the cross sections obtained using the three approaches discussed in the text, for the six selected benchmark points. BPs 1-4 correspond to scenario 1 and BPs 5 and 6 to scenario 2. Blank space in front of a quantity implies that it is not relevant for the given BP. overlapping of the two lines occurs even in the bin containing the three heavy Higgs bosons (so that the red lines, plotted first, are entirely invisible behind the blue ones), since no interference effects appear there. Finally, we took the distributions for the differential σ BW and σ Int for the four BPs of scenario 1 and convolved them with a Gaussian of 150 GeV width, using the ListConvolve function [41] in Mathematica. The purpose of this exercise is to emulate the detector effect of smearing the cross sections over a few adjacent bins, to see whether the resulting distributions can be visually distinguishable from each other, despite the poor mass resolution. Figure 5 shows these convolved distributions. We note that the convolution spreads out the differences in the heights of the red and blue boxes to the bins around the ones containing the resonant Higgs masses. However, the shapes of these convolved distributions for σ Int do not show any peculiar features. They could very well be the distributions of the differential σ BW for slightly different parameter space points, with the peaks arising from a single heavy Higgs resonance or even from nearly mass-degenerate H and A (with the underlying model being simply the MSSM). Therefore, even with an integrated luminosity as large as 6000 fb −1 (as assumed in this figure, in order to reduce the sizes of the error bars, which account for the statistical error only), the LHC will not be able to exploit the interference effects in order to identify two Higgs resonances with highly identical masses.

V. CONCLUSIONS
The commonly adopted approach of calculating the cross section for a given 2 → 2 process by factorising it into the production and decay parts, assuming a narrow width of the mediator, by construction cannot account for the possible quantum interference among the propagators of several mass-degenerate states. In this study, we have considered the specific example of the NMSSM, wherein nearly identical-mass pairs of CP-even or CP-odd Higgs bosons are viable over substantial regions of the parameter space. These regions were found by numerical scanning of broad ranges of the model parameters, while imposing the most important experimental constraints, including those from the LHC pertaining to the Higgs, exotic and flavour sectors, as well as those from the DM searches.
By analysing six illustrative benchmark points from the scanned set, we have highlighted the importance of taking the interference effects into account. This was done by including the full propagator matrix in the calculation of the cross section for the process of production of τ þ τ − in gluon fusion at the 14 TeV LHC. We have shown that this cross section can deviate considerably from the one obtained by employing the NWA, and even from the one obtained assuming BW propagators, an approach most experimental searches are based on. This deviation, in fact, implies a reduction in the cross section in the case of two massdegenerate CP-even Higgs bosons, as the interference is always destructive. In the case of CP-odd states, on the other hand, no interference effects appear. We have also reasserted the fact that the smaller the mass-splitting between two nearby Higgs bosons compared to the sum of their widths, the larger the interference effects.
The reason for considering the τ þ τ − decay channel of the Higgs bosons is that the γγ channel, while cleaner, has a prohibitively small decay rate for Higgs bosons with masses close to 1 TeV, The BR for the τ þ τ − is significantly larger, but the price to pay is a poor experimental resolution and mass reconstruction. As a result of which, we have concluded that the LHC will be unable to disentangle the two resonances, even if they have a mass splitting of a few GeV and the integrated luminosity is as large as 6000 fb −1 . Thus, even though the gg → τ þ τ − channel is an effective search instrument for heavy Higgs bosons, in the particular scenario when the latter are highly degenerate in mass, the former is unable be resolve them. Alternative Higgs boson production and decay modes that might be more suitable for this purpose are therefore the subject of a future study.   5. Distributions of the differential cross sections for the four selected BPs of scenario 1, after convolution with Gaussians of width 150 GeV. The color convention for the lines is the same as in Fig. 4, and the error bars on them correspond to an assumed integrated luminosity of 6000 fb −1 .