Leptoquark-assisted Singlet-mediated Di-Higgs Production at the LHC

At the LHC, the gluon-initiated processes are considered to be the primary source of di-Higgs production. However, in the presence of a new resonance, the light-quark initiated processes can also contribute significantly. In this paper, we look at the di-Higgs production mediated by a new singlet scalar. The singlet is produced in both quark-antiquark and gluon fusion processes through loops involving a scalar leptoquark and right-handed neutrinos. With benchmark parameters inspired from the recent resonant di-Higgs searches by the ATLAS collaboration, we examine the prospects of such a resonance in the TeV-range at the High-Luminosity LHC (HL-LHC) in the $b\bar{b} \tau^{+}\tau^{-}$ mode with a multivariate analysis. We obtain the $5\sigma$ and $2\sigma$ contours and find that a significant part of the parameter space is within the reach of the HL-LHC.


I. INTRODUCTION
The current measurements of the Higgs boson cross section and its couplings agree with the Standard Model (SM) predictions. Still, new physics (NP) could be hiding in the Higgs sector if it primarily affects the Higgs couplings with the third generation quarks and vector bosons up to a level of 10%-20%. The prospects are much more open if the effects are exclusive to the first-and second-generation quark couplings, as the bounds on these are much more relaxed. NP could also show up in the Higgs self interaction. At the LHC, a direct probe of the Higgs self interaction is the di-Higgs production through gluon fusion, gg → hh. At the leading order (LO), the self-coupling appears in a virtual-h-mediated process where the virtual Higgs is produced through quark triangles. There are also box diagrams, independent of the Higgs self coupling. In the SM, the di-Higgs cross section is almost three orders of magnitude smaller than the single Higgs production. This is because the box diagrams interfere destructively with the triangle diagrams [1][2][3][4][5], and producing two Higgses requires more energy than producing one (phase-space suppression). Ignoring the finite topquark-mass effects, one gets σ (gg → hh) ≈ 33 fb with √ s = 13 TeV for m h = 125 GeV [6][7][8][9][10][11][12] at the next-to-next-to-leading order (NNLO) in α s . This accidental suppression of the di-Higgs production cross section makes the process a good candidate to probe for NP in the Higgs sector.
In this letter, we extend this idea and consider the pp → hh process mediated by an s-channel SM-Singlet scalar resonance. Except for the Higgs, the singlet scalar couples to the SM fields only through its interaction with a scalar LQ. There are mainly three motivations for considering this setup. First, such a boson is well-motivated in many NP scenarios. Second, since, unlike the Higgs, a singlet scalar cannot couple with the left-handed fermions, producing them at the colliders is not straightforward. Hence, studying the LHC phenomenology and prospects of a coloured-scalar-assisted production of a singlet scalar that dominantly decays to a Higgs pair is interesting. Third, the recent ATLAS searches for the resonant di-Higgs production (pp → X → hh) show a small excess of about 3.2σ local significance around M X ≈ 1 − 1.1 TeV [99,100]. The analyses were done at √ s = 13 TeV with 139 fb −1 integrated luminosity for different final states: bbbb, bbτ + τ − , bbγγ, and a combined one. Of these, the data in the bbγγ mode go only up to a TeV. The enhancement is prominently visible in the bbτ + τ − and combined modes. The data put an upper bound on resonant hh production cross section. The CMS collaboration has also put similar upper limits on the hh cross section from various searches with 35.9 fb −1 and 138 fb −1 of data [101,102]. 1 Since the excess is not statistically significant, it could be just a statistical fluctuation. However, we take it as a motivation for choosing our model parameters to study the LHC phenomenology of the singlet scalar.
In principle, LQs can be either scalars or vectors in local quantum field theories. Here we consider a weak-singlet charge −1/3 scalar LQ (commonly known as S 1 ) for our purpose. The basic framework is the same as that in Ref. [81] where we considered a simple extension of the SM augmented with an S 1 LQ and three generations of RH neutrino ν R . There we showed that the heavy neutrinos can induce a boost to the down-typequark Yukawa interactions through a diagonal coupling with the quarks and S 1 . The enhanced Yukawa couplings can also be realised through a radiative generation of the dimension-6 operator of the form f d (H † H/Λ 2 ) (q L Hd R ) + H.c. (with Λ ∼ TeV) which can enhance the single Higgs production rate through the qq → h processes. In this framework, the coefficient f d is restricted by LHC limits on the masses and the couplings of new fields (e.g. LQs, ν R ) to the SM ones. However, the limits become weaker, if the SM Higgs boson is replaced by a new scalar φ . Hence, a large production rate can be obtained for the qq → φ processes. Through the same process, a resonant φ can generate the hh final state at the LHC, which is further enhanced at the resonance. We will see if any significant excess in the 2h cross-section can be observed if the gluon fusion process is supplemented by the quark fusion. We will also investigate the prospects of this process at the HL-LHC.
The plan of the letter is as follows. In the next section, we review the model with the S 1 LQ; in Sec. III, we discuss the di-Higgs production in the SM and via the new scalar resonance; in Sec. IV, we discuss the prospects of the channel at the HL-LHC and conclude in Sec V.

II. A RECAP OF THE LQ MODEL
The model in Ref. [81] is an extension of the SM with chiral neutrinos and an S 1 -type scalar LQ. The LQ transforms under the SM gauge group as 3 , 1, 1/3 with Q EM = T 3 +Y . For our purpose, we introduce a real SM singlet scalar φ . The general fermionic interaction Lagrangian for S 1 in addition to the kinetic term can be written in the notation of Ref. [82] as, where D µ = ∂ µ + ig S T a G a µ + iY B µ is the covariant derivative; T a and G a µ are the colour generators and the gluon fields, respectively; B µ is the gauge mediator of hypercharge Y and g S is the strong coupling. We have suppressed the colour indices. The superscript C denotes charge conjugation; {i, j} and {a, b} are flavour and SU(2) indices, respectively. The SM quark and lepton doublets are denoted by Q L and L L , respectively. We add the scalar interaction terms to the Lagrangian in Eq. (1) as, Here, H denotes the SM Higgs doublet, andM S 1 and M φ define the bare mass parameters for S 1 and φ respectively, and λ and µ are dimension-one couplings proportional to some NP scale. For simplicity, we shall assume µ → 0. We denote the physical Higgs field after the electroweak symmetry breaking as h ≡ h 125 . The physical masses can be obtained via where v 246 GeV is the vacuum expectation value (VEV) of the SM Higgs. The physical mass of S 1 is given as In the subsequent discussion, we set the quark and neutrino mixing matrices to identity, assume the LQ couplings to be flavour diagonal, and put y RR 1 = 0 to avoid flavour bounds [103]. We write y LL 1 = g L and y RR 1 = g R for simplicity. For a closer look, we can expand Eq. (2) for the first generation as where we have simplified y X 1 ii as y X i . Since the flavour of the neutrino is irrelevant for the LHC, we simply write ν to denote a neutrino. The LQ-gluon couplings are generated from the covariant derivative term: where, G µ = T a G a µ . The second term of Eq. (6) produces the gS 1 S 1 vertex [vertex factor ∼ −g S (p + p ) µ ] while the last term is responsible for ggS 1 S 1 four point interaction [vertex factor ∼ g 2 S ]. Here p and p are the incoming and outgoing momenta of the of S 1 .

A. Mixing between φ and the SM Higgs
After the electroweak symmetry breaking, h and φ mix through the term, µ H † H φ . In the (φ h) T basis, the mass matrix can be written as, We can rotate (φ h) T to a physical basis (H 1 H 2 ) T as: In the new basis, the mass matrix is diagonal and it is defined as, where U is an orthogonal rotation matrix and The mixing angle is given as, For a small mixing angle, we can assume H 1 to be mostly singlet-like while H 2 will show a strong doublet-like nature. Using Eq. (8), we can recast the important interaction terms as, where s θ = sin θ and c θ = cos θ .
In the subsequent analysis, we assume c θ → 1 and s θ → 0. In this limit, we assign m h SM M H 2 = m h and M φ = M H 1 in the following sections. The above assignment is valid only if µ → 0 or M φ m h > µ . In our case, the second condition is applicable as we are primarily interested in a new scalar resonance at ∼ 1.1 TeV in the di-Higgs production rate.

III. DI-HIGGS PRODUCTION THROUGH A SCALAR RESONANCE IN THE LQ MODEL
We first look at the leading-order (LO) SM contributions to the di-Higgs production at the LHC.
A. Di-Higgs production in the SM Here we start with reviewing the dominant production processes for the gg → hh [1][2][3][4][5], as shown in Fig. 1. The amplitude for where a, b are the colour indices,ŝ,t,û are the parton-level Mandelstam variables. The independent tensor projections P µν and Q µν are defined as, where p T = (ût − m 4 h )/ŝ is the transverse momentum of an outgoing Higgs. The triangle diagram contributes only to F (ŝ,t,û, m 2 t ), while box diagrams contribute to both F (ŝ,t,û, m 2 t ) and G (ŝ,t,û, m 2 t ). The triangle diagram has no angular momentum dependence, thus, it is an s-wave contribution. The form factors F box and G box can be attributed to the spin-0 and spin-2 parts of the box amplitude, respectively. The angular dependencies of the form factors by considering the partial wave decomposition of the scattering amplitude are discussed in Refs. [9,104]. We can split the SM contribution to Here we have used the following short hands for three-point and four-point functions, We have ignored all other quark loops except the t-quark loop since it gives the dominant contribution. Similarly, G SM (ŝ,t,û, m 2 t ) can be expressed as, As shown, the leading order (LO) amplitude includes the topquark-mass-dependent terms. The contribution at the next-toleading order (NLO) in the strong coupling constant both in the infinite top-mass limit and with the full top-mass dependence are considered in Refs. [8,9,11,12,62,105]. For higher-order results in different top mass limits, see Refs. [6,[106][107][108][109].
B. The gg → hh process in presence of φ The diagrams for the gg → hh through the exchanges of a scalar s i where s i = φ , h are shown in Fig. 2. The complete invariant amplitude for the gg → hh processes can be obtained by considering the contributions from both the triangle [ Fig. 2(a)] and bubble [ Fig. 2(b)] diagrams: The loop function in the gg s i effective vertex is given by, For the SM Higgs, λ i = λ v and µ i = 3m 2 h /v, while for the singlet scalar φ , λ i = λ and µ i = µ /2. The total gg → hh differential cross section including all the resonance channel diagrams (SM+NP) and the box diagrams (SM only) can be calculated as, C. The qq → hh process Fig. 3(a) shows a generic diagram for the qq → hh processes, where s i = {φ , h} and q is a down-type quark. The blob represents the s i qq effective vertex. The effective couplings are available in Ref. [81] which we repeat below for completeness.
• For s i = φ : For the singlet scalar, there is only one possible vertex, as shown in Fig. 3(b). The corresponding effective coupling is given by, where we use a compact notation, C q = g L g R y ν .  (31)]. Hence, we compare it with the upper limit on σ (X → hh) (dashed line) as observed by the ATLAS collaboration in the combined channel [100]. The black diamond signifies the parameter choice for which σ (pp → φ ) ≈ 14.5 fb at M φ = 1.1 TeV.
• For s i = h: For the SM Higgs, there are three possible one loop diagrams along with the tree-level contribution [81]. The total effective coupling for the qqh vertex is given by, where D 0 , C 0 and B 0 are the Passarino-Veltman fourpoint, triangle, and bubble integrals, respectively. The y SM q term is the SM tree-level contribution. The differential cross section for qq → hh at the leading order is given by, (28)

D. The resonant φ contribution and choice of parameters
By replacing the s i propagator in Eqs. (24) and (28) by the Breit-Wigner shape, we obtain the gluon and quark contributions to the di-Higgs process including the SM-NP interference. However, since we focus on a TeV-range φ , the φ resonance appears far from the Higgs contribution, making the interference negligible. Hence, in this case, it is possible to look at the resonant φ contribution separately. This is also justified by the fact that in the parameter range of our interest, the decay width of φ is quite narrow, Γ φ /M φ < 1% (for example, for M S 1 = 1 TeV, M ν R = 0.75 TeV, C q = 3.5, λ = 3.5 TeV and µ = 200 GeV, we get Γ φ 1 GeV). As a result, we could use the narrow-width approximation to estimate the resonant φ contribution to the di-Higgs production.
In other words, where BR stands for branching ratio. For µ 100 GeV, the tree-level decay of φ dominates, making BR(φ → hh) ≈ 1. We take a conservative µ = 200 GeV to keep θ small. Hence, we can further simplify the above expression as We show σ (pp → φ ) at the √ s = 13 TeV LHC in Fig. 4. The combination C q = g L g R y ν parametrises the quark-initiated contribution. The quark contribution switches off for C q = 0. Since σ qq is proportional to C 2 q , σ (pp → φ ) increases rapidly with increasing C q . For this plot, we use a benchmark set of parameters, M ν R = 0.75 TeV, M S 1 = 1 TeV and λ = 3.5 TeV. Our choice of the right-handed neutrino mass is motivated by Ref. [110,111]. The current direct search limits on S 1 lie between 1.7 − 1.8 TeV [112]. These limits are obtained assuming that the S 1 decays exclusively to one or two final states (like, e.g., e j or e j and ν e j). However, in our case, the S 1 can decay to ue, cµ,tτ, dν, sν, bν, dν R , sν R , sν R final states, making the BR in each light mode (i.e., those without ν R -there is no direct limit available from these exotic decay modes) to be less than 1/6 [81]. As a result, a TeV S 1 is allowed. Moreover, the combination C q = g L g R y ν can be order one even when g L is less than one. Thus, it is possible to take M S 1 = 1 TeV and still avoid the constraints from low-energy observables and the singleproduction searches. Finally, the dimension-one coupling λ comes from a NP scale which we assume to be heavier than the entire particle spectrum of our model. The particular choice is motivated by the fact that for λ ∼ 3.5 TeV, σ (pp → φ ) at 13 TeV fits the small excess seen by ATLAS [100] at M φ = 1.1 TeV with perturbative couplings.

IV. LHC PHENOMENOLOGY
To study the prospects of our model at the LHC, we consider the bbτ + τ − channel where one of the Higgs boson decays to a bb pair and the other to a pair of τ's. We use MADGRAPH5 [120] to generate the signal and background events at the tree level. We implement the following effective Lagrangian in FEYN- where, µ hφ and Λ g are the dimension-one couplings of the singlet φ to the SM Higgs and gluons, respectively. To estimate the signal strength for our LHC analysis, we use the Kfactors, K qq (gg) = 1.3 (1.72) from Ref. [62]. We incorporate the higher order corrections to the SM background processes as well. We multiply the LO cross sections with the appropriate QCD K-factors known in the literature (see Table I). We use the NNPDF2.3LO [123] parton distribution functions with the default renormalisation and factorisation scales. For showering and hadronisation, these parton level events are passed through PYTHIA6 [124]. The detector simulations are performed in DELPHES 3.5.0 [125] using the ATLAS card, incorporating the latest b quark and τ lepton identification efficiencies 3 [126,127]. We use the FastJet package [128] for jet clustering with the anti-k T algorithm [129] with the radius parameter R = 0.4. 2 Since we consider the φ to be about as heavy as the LQ running in the loops, one must be careful in using the effective Lagrangian for computing cross sections. Therefore, we use it only to simulate the kinematic distributions of the Higgs pair, not the signal cross sections. 3 We take the b-quark tagging efficiency to be a constant beyond the p T regime mentioned in Ref. [126].

FIG. 5. Pearson's linear correlation coefficients in the signal for the inputs chosen for BDT analysis.
Feature 0443   TABLE III. Summary of method-unspecific and method-specific ranking for the input features at the benchmark mass point.
A simple kinematic-cuts-based analysis is not very promising in this case, as the SM backgrounds are quite significant. Hence, we employ a multivariate analysis that uses Boosted Decision Trees (BDTs) to estimate the signal significance at the HL-LHC.
These cuts lead to a significant drop in the background events and a moderate drop in the signal events. After these cuts, the efficiencies of quark mediated and gluon mediated signal events are similar (≈ 4.5%). Demanding two b quarks essentially eliminates the W backgrounds. We find that only the tt backgrounds remain significant after these cuts. Using the event-level information and the kinematic features of identified objects, we define the following basic set of quantities: where m T (b 1 , b 2 ) is the transverse mass of the bb-system (both transverse and invariant mass of the bb-system have similar distributions and are highly correlated, we pick m T since it has slightly better separation visibly), ∆R(x, y) is the separation between x and y in the η − φ plane,M φ is the estimator of the invariant mass of the singlet φ constructed from the tagged τ hleptons and b-quarks and H T is the sum of the p T of all jets. We picked these features from a larger set of kinematic variables created from the identified objects after dropping the highly correlated ones (linear correlation > 70% in the signal)-we pick the one with the greatest separation out of a highly co-related pair. The Pearson correlation coefficients between the final set of features are shown in Fig. 5.
We use these features as inputs for the BDT analysis with the adaptive boosting (AdaBoost) optimisation. After obtaining statistically independent event samples for all the signals and background processes considered, we take 30% of the dataset for testing and the rest for training the algorithm. Since there are multiple processes in each of the signal and background classes, we weigh each by the ratio of the expected number of events from the process at the HL-LHC to the number of events fed into the BDT algorithm. We choose the benchmark mass, M φ = 1 TeV to optimise the BDT algorithm, which we implement using the TMVA [130] package in ROOT. The Kolmogorov-Smirnov test ensures that the signal and background distributions are well separated-see Fig. 6 for the normalised BDT response and classifier efficiency. The set of optimised BDT hyperparameters are given in Table II. To show the discriminating power of this set of features in the context of BDTs, we list two measures for each feature: method-unspecific separation and method-specific ranking. The method-unspecific separation which is defined as, whereŷ S andŷ B are the probability density functions of the signal and background respectively for a particular feature y.
For signal and background distributions with no overlap, it is unity-features with higher method-unspecific separation values are better at resolving the signal events from the background. Method-specific ranking is derived by counting how often the variables are used to split decision tree nodes, and by weighting each split occurrence by the separation gain-squared it has achieved and by the number of events in the node. We show both these quantities in decreasing order of discriminating power of the features in Table III. The significance is computed using the formula, is the number of signal (background) events surviving at the HL-LHC after the BDT cut. We use the optimised analysis setup to study the collider reach for different masses of φ in the range 0.8−1.6 TeV, for the benchmark λ and C q . Based on our parametrisation, the number of signal events scales with λ and C q as, whereN gg(qq) s is the number of gg(qq)-initiated events surviving the BDT cut at 3 ab −1 when λ = 1 TeV (and C q = 1). We present the 5σ and 2σ limits in terms of λ and C q in Fig. 7. One can easily translate these limits in terms of the signal cross section, correcting for the efficiency of each process. We see that for C q = 3.5, λ = 3.5 TeV, a 1.1 TeV φ (i.e., the parameters to accommodate the current ATLAS excess) can be discovered at the HL-LHC-the signal significance for this point is slightly above 5σ .

V. SUMMARY AND CONCLUSIONS
In this letter, we have studied the resonant di-Higgs production mediated by a SM-singlet scalar, φ at the LHC. In our model, the only SM particle the singlet directly couples to is the Higgs. It couples to S 1 , the charge −1/3 weak-singlet scalar Leptoquark, through which it gains an effective coupling with gluon pairs. It also couples to down-type quark-antiquark pairs through loops involving S 1 and ν R . We showed in an earlier paper [81] that the quark fusion contribution to σ (pp → φ ) could be comparable, or even more significant, than the gluon fusion channel for a suitable choice of parameters. Moreover, since the singlet φ has tree-level couplings to the Higgs, it dominantly decays to a pair of Higgses. As a result, we have studied an interesting setup where the φ can be substantially produced at the LHC through both quark and gluon fusion channels and decays to the di-Higgs final state. The parameter ranges for our study were motivated from the (3.2σ local significance) excess in the di-Higgs spectrum around M hh ≈ 1.1 TeV observed by the ATLAS collaboration [99,100] at the 13 TeV LHC.
We have investigated the prospects of resonant di-Higgs production in our model at the HL-LHC. For this purpose, we have chosen the bbτ + τ − final state (in which the ATLAS excess was observed [99]). To overcome the huge SM background in this channel, we performed a multivariate analysis using Boosted Decision Trees to estimate the signal significance at the HL-LHC. Our study has shown that if we choose the model parameters to account for the excess at the 13 TeV LHC, the HL-LHC can discover such a resonance. However, there are other variants of multivariate analysis (e.g., BDTs with different boosting techniques, deep-neural-network-based analyses, etc., for e.g. see [131,132]) that could lead to a better signal yield at the HL-LHC. Moreover, with better b-and τ-tagging efficiencies in the future, one could improve this further. A similar study can also be performed in the bbbb mode, which we postpone for a future study. We have presented the results of our study for a wide range of model parameters that can be easily translated to cross sections and, hence, useful for other similar models.

VI. ACKNOWLEDGMENTS
Our computations were supported in part by SAMKHYA: the High Performance Computing Facility provided by the Institute of Physics (IoP), Bhubaneswar, India. C. N. is supported by the DST-Inspire Fellowship.