Four of a kind: exploring quadruple-Higgs boson production at the HL-LHC

More than a decade after its discovery, the Higgs boson still has much to tell us. Many of its properties have been measured with increasing precision, but we know far less about how Higgs bosons interact with one another. These self-interactions offer a way to probe the shape of the Higgs potential, which governs electroweak symmetry breaking. Measuring them could reveal whether the Higgs sector behaves as the Standard Model predicts or contains new particles and interactions.

The Higgs potential may also hold clues about the early Universe. Changes to its shape could have affected the electroweak phase transition, when particles acquired mass through their interaction with the Higgs field. If this transition was strongly first order, it could have helped create the imbalance between matter and antimatter that we observe today.

At the LHC, one way to explore the Higgs potential is to produce several Higgs bosons in a single collision. Two interactions are especially important. The trilinear coupling involves three Higgs bosons, while the quartic coupling involves four. Their strengths relative to the Standard Model predictions are written as κ3\kappa_3 and κ4\kappa_4, with both equal to one in the Standard Model.

Higgs-boson pair production (HH) provides the most direct route to κ3\kappa_3. The process has not yet been observed, but the ATLAS and CMS combination of Run 2 searches [1] sets a 95% confidence-level upper limit of 2.5 times the predicted Standard Model rate. Assuming that the other relevant couplings have their Standard Model values, the result constrains κ3\kappa_3 to between −0.71 and 6.1.

To probe κ4\kappa_4 directly, we need to produce three Higgs bosons (HHH). This process is extremely rare. At 13 TeV, its predicted Standard Model cross section is only about 0.08 fb, roughly 400 times smaller than that of Higgs-boson pair production. Its rate also depends on both κ3\kappa_3 and κ4\kappa_4, making their effects difficult to separate. Despite these challenges, triple-Higgs boson production has now entered the experimental programme.

In 2024, ATLAS presented the first LHC search for triple-Higgs boson production [2], using 126 inverse fb of Run 2 data at 13 TeV. The search focused on HHH → 6b, in which all three Higgs bosons decay into bottom-quark pairs. In earlier theoretical studies, we proposed exploring this signature in Standard Model triple-Higgs production at a future 100 TeV proton collider [3] and in resonantly enhanced triple-Higgs production at the LHC [4].

The large branching fraction helps. A 125 GeV Higgs boson decays into bb̄ about 58% of the time, so roughly one in five HHH events produces six bottom quarks. Even so, only around two Standard Model events of this kind would have been produced in the full dataset before acceptance and selection. Finding them among the enormous QCD multijet background is therefore a formidable challenge.

ATLAS searched for the signal in events with at least six b-tagged jets, using events with exactly four or five b-tagged jets to estimate the background. The jets were paired to reconstruct three Higgs-boson candidates, while neural networks helped distinguish potential signal events from the background.

No significant excess was found. ATLAS set a 95% confidence-level upper limit of 59 fb on the production rate, corresponding to 760 times the Standard Model prediction. Varying one self-coupling at a time gave:

11<κ3<17(κ4=1),230<κ4<240(κ3=1).-11 < \kappa_3 < 17 \qquad (\kappa_4 = 1), \qquad -230 < \kappa_4 < 240 \qquad (\kappa_3 = 1).

These were the first direct experimental constraints on the quartic Higgs self-coupling.

CMS has since reported a more sensitive search [5], using 138 inverse fb of Run 2 data. At the heart of the analysis is SPANet, a neural network designed to assign the jets to the most likely Higgs-boson candidates while separating signal from background. The search considered both resolved events, in which the bottom quarks from Higgs-boson decays produce separate jets, and events in which the decay products of a boosted Higgs boson are reconstructed as a single large-radius jet. The resolved events provided most of the sensitivity.

CMS also found no significant excess. Its observed upper limit of 44 fb, or 588 times the Standard Model prediction, led to the 95% confidence-level intervals

7.4<κ3<12.4(κ4=1),177<κ4<185(κ3=1).-7.4 < \kappa_3 < 12.4 \qquad (\kappa_4 = 1), \qquad -177 < \kappa_4 < 185 \qquad (\kappa_3 = 1).

For the first time, a direct non-resonant HHHHHH search excluded part of the Higgs self-coupling parameter space allowed by perturbative unitarity bounds from Higgs-boson scattering.

ATLAS also explored a different route: ppXSHHHHpp\rightarrow X\rightarrow SH\rightarrow HHH. Here, two new scalar particles decay in succession, producing three Higgs bosons through a resonant cascade. Such processes occur, for example, in the Two Real Singlet Model, which extends the Standard Model with two additional real scalar fields that are singlets under its gauge symmetries. These processes can greatly increase the rate. The search tested masses up to about 1.5 TeV for X and 1 TeV for S, but found no significant excess.

Together, these searches constrain scenarios in which triple-Higgs boson production is strongly enhanced. The bounds on the self-couplings still allow large departures from the Standard Model and depend on assumptions about other interactions, including a Standard Model-like top-quark Yukawa coupling. Could an even higher Higgs-boson multiplicity provide additional information?

Adding a fourth Higgs boson comes at a steep price in the Standard Model. For gluon-fusion production at 14 TeV, our study [6] gives the leading-order prediction:


σ(ggHHHH)1.24×104fb\sigma(gg \to HHHH) \simeq 1.24 \times 10^{-4}\,\mathrm{fb}

about 350 times smaller than the corresponding leading-order triple-Higgs boson production rate. Even with the projected integrated luminosity of 3000 inverse fb at the High-Luminosity LHC (HL-LHC), fewer than half an event would be produced before the Higgs bosons decay.

We focused on the case in which all four Higgs bosons decay into bottom-quark pairs, producing a final state with eight bottom quarks. Assuming Standard Model decays, this happens in approximately 11% of HHHHHHHH events and is the largest individual decay channel once the subsequent decays of W and Z bosons are specified. Yet even this favourable branching fraction leaves only about 0.04 Standard Model events before detector acceptance and event selection.

Finding the signal would also pose a reconstruction challenge. With eight resolved jets, there are already 105 different ways to arrange them into four pairs, even before mass or kinematic information is considered. The task is to identify which pairing, if any, is compatible with four Higgs-boson decays.

These numbers put an observation of Standard Model four-Higgs boson production beyond the reach of the HL-LHC. That does not make the process irrelevant. Its dependence on κ3\kappa_3 and κ4\kappa_4 differs from that of HHHH and HHHHHH production, and interference among the contributing diagrams can make the rate less suppressed in some regions of parameter space. New scalar resonances could increase it much more dramatically. We examined both possibilities.

The ATLAS six-b signal region in HHHHHH searches asks for at least six b-tagged jets. Those two words leave room for extra company: an HHHH8b HHHH \rightarrow 8b event can enter the same region when six, seven, or all eight of its bottom-quark jets are reconstructed and tagged.

How often could this happen? Comparing the two all-bottom decay channels, we find that

σ(ggHHHH8b)σ(ggHHH6b)\frac{\sigma(gg \to HHHH \to 8b)}{\sigma(gg \to HHH \to 6b)}

lies between 1% and 10% across substantial regions of the (κ3\kappa_3, κ4\kappa_4) plane away from the Standard Model, as shown in Figure 1. In these regions, adding a fourth Higgs boson does not reduce the rate as dramatically as the Standard Model numbers might suggest.

Figure 1: Comparing four-Higgs-boson and three-Higgs-boson events. The contours show the ratio of production cross sections, including decays to eight and six bottom quarks, respectively, at 14 TeV. The axes use c3=κ31c_3 = \kappa_3 – 1 and d4=κ41d_4 = \kappa_4 – 1, placing the Standard Model at (0, 0).

The additional bottom-quark jets can also help distinguish the signal from the large QCD background, although mistagged jets and other background processes must still be considered. With a loose selection inspired by the ATLAS analysis, we find that HHHH events could contribute between 10% and 100% of the selected signal rate from HHHHHH and HHH+bbHHH+b\bar{b} in some regions with strongly modified self-couplings.

This means that existing triple-Higgs boson searches may already contain useful information about four-Higgs boson production. Including the HHHHHHHH contribution could provide initial constraints without requiring a dedicated event category. It could also be important when translating the results of these searches into bounds on the Higgs self-couplings.

An eight-bottom final state does not necessarily come from four Higgs bosons. Two or three Higgs bosons produced with additional bottom quarks can leave the same signature. We therefore combined HHHHHHHH, HHH+bbHHH + b\bar{b} , and HH+4bHH+4b in a single signal model. Our analysis reconstructs four Higgs-boson candidates from separate jets and uses a machine-learning classifier to distinguish signal from the QCD multijet background.

At the Standard Model point, only about 0.4 signal events remain after selection, compared with approximately 22 background events. Around 98% of the signal comes from HH+4bHH+4b. The genuine HHHHHHHH contribution is tiny, but can grow substantially when the self-couplings move away from their Standard Model values.

Allowing both couplings to vary, the projections of our simultaneous 95% confidence region reaches

10.3<c3<11.2,222<d4<234,-10.3 < c_3 < 11.2, \qquad -222 < d_4 < 234,

along the two axes, where c3=κ31c_3 = \kappa_3 – 1 and d4=κ41d_4 = \kappa_4 – 1. The result is shown in Figure 2 and includes statistical uncertainties only.

Different combinations of the two self-couplings can produce similar event yields. Another final state may respond differently and help separate these possibilities. The shape of the allowed region can therefore be as informative as its overall size.

Figure 2: Projected constraints on the Higgs self-couplings at the HL-LHC. The red curve shows the simultaneous 95% confidence contour from our eight-b signal model, while the blue curve shows the ATLAS HHH6bHHH \rightarrow 6b projection [7]. Neither includes systematic uncertainties. The pale-red band shows the effect of varying our background estimate by a factor of four; it is not a systematic-uncertainty band. The black dashed curve marks the perturbative unitarity bound from HHHH scattering, and the star marks the Standard Model point.

Our projected reach for the trilinear coupling is comparable to that of the ATLAS analysis, while ATLAS obtains a narrower range for the quartic coupling. The different contour shapes suggest that the eight-b-jet analysis could provide complementary information in a combined determination of the two self-couplings.

These results give a first look at what an eight-b-jet search could achieve. A more realistic assessment will require detailed detector simulation and a careful study of systematic uncertainties in both the signal and background.

Heavy scalar particles offer another route to four Higgs bosons. Their direct or cascade decays could increase the rate by many orders of magnitude. In some models, events with several Higgs bosons could even provide the clearest sign of these new particles.

A recent study [8] shows how this can happen near the alignment limit. In this limit, the observed Higgs boson retains its Standard Model-like behaviour even though additional scalars are present. Small mixing can suppress the new particles’ decays into fermions and gauge bosons, while other interactions or decay chains can still produce several Higgs bosons. Suitable, although somewhat fine-tuned, choices of couplings can also suppress decays into Higgs-boson pairs.

We considered two examples, shown in Figure 3. A new scalar S can decay directly into four Standard Model-like Higgs bosons, ppSHHHHpp \to S \to HHHH or a heavier scalar h3h_3 can decay into two lighter scalars h2h_2, each producing a Higgs-boson pair, pph3h2h2(HH)(HH)pp \to h_3 \to h_2 h_2 \to (HH)(HH).

The cascade provides an extra clue when reconstructing the event. Each pair of Higgs bosons comes from the same lighter scalar and can reveal its mass. Combining all four Higgs bosons can then reveal the heavier scalar. In the direct decay, only the common parent of all four can be reconstructed.

Figure 3: Two possible routes to four Higgs bosons. A new scalar can decay directly into four Higgs bosons (left), or a heavier scalar can decay through two intermediate scalars, each producing a Higgs-boson pair (right). The labels h and h1 denote the observed Higgs boson in the two scenarios.

The masses of the new particles determine how the events appear in the detector. For lighter resonances, the bottom quarks are more likely to form separate, resolved jets. Heavier resonances can produce boosted Higgs bosons whose decay products merge into a single jet. We combined resolved and boosted reconstruction methods to cover both possibilities.

For the direct decay, our projected 95% confidence-level upper limits range from about 1.3 to 11 fb over the masses studied. At 1.5 TeV, for example, we obtain:

σ(ppS)×BR(SHHHH)<2.5fb\sigma(pp \to S)\times BR(S \to HHHH) < 2.5\,\mathrm{fb}

For the cascade, at Mh2=625GeVM_{h_2}=625\,\mathrm{GeV} and Mh3=1.5TeVM_{h_3}=1.5\,\mathrm{TeV}, the corresponding result is

σ(pph3)BR(h3h2h2)[BR(h2HH)]2<2.6fb.\sigma(pp\to h_3)\, \mathrm{BR}(h_3\to h_2h_2)\, \left[\mathrm{BR}(h_2\to HH)\right]^2 <2.6\,\mathrm{fb}.

Here, BR denotes a branching fraction, and both results refer to four Higgs bosons before they decay. These values provide useful benchmarks; testing a specific model also requires its predicted rates and branching fractions, together with an assessment of the relevant uncertainties.

Observing Standard Model quadruple-Higgs boson production is out of reach at the HL-LHC. Four-Higgs-boson events can nevertheless matter when modified self-couplings enhance their rate: some may enter inclusive triple-Higgs boson signal regions and affect how those searches are interpreted. This gives a practical reason to include HHHHHHHH contributions when testing large departures from the Standard Model.

An eight-b-jet search offers another way to explore these departures. Its signal can contain HH+4bHH+4b, HHH+bbHHH + b\bar{b}, and HHHH, so the three contributions need to be considered together. The different sensitivity contours of the eight- and six-b-jet searches suggest that the two channels could provide complementary information, although a combined analysis would be needed to establish the benefit.

The prospects become especially exciting if new scalar particles exist. Their direct or cascade decays could dramatically enhance triple- and quadruple-Higgs boson production, particularly when more conventional decay channels are suppressed. Multi-Higgs boson final states may then provide the first sign of a richer scalar sector. Our projected rate limits provide benchmarks for testing such scenarios, although establishing their experimental reach will require more detailed detector simulations and a fuller treatment of uncertainties.

These challenging but promising signatures are compelling targets not only for the LHC and its high-luminosity phase, but also for future higher-energy colliders, where larger production rates could bring the detailed exploration of higher Higgs-boson multiplicities within reach. Four Higgs bosons are unlikely to be produced often—but when they are, they may have an unusually interesting story to tell.

[1] ATLAS and CMS Collaborations, Combination of ATLAS and CMS searches for Higgs boson pair production at √s = 13 TeV (2026). arXiv:2602.23991.

[2] ATLAS Collaboration, Search for triple Higgs boson production in the 6b final state using pp collisions at √s = 13 TeV with the ATLAS detector, Phys. Rev. D 111, 032006 (2025). arXiv:2411.02040.

[3] A. Papaefstathiou, G. Tetlalmatzi-Xolocotzi and M. Zaro, Triple Higgs boson production to six b-jets at a 100 TeV proton collider, Eur. Phys. J. C 79, 947 (2019). arXiv:1909.09166.

[4] A. Papaefstathiou, T. Robens and G. Tetlalmatzi-Xolocotzi, Triple Higgs Boson Production at the Large Hadron Collider with Two Real Singlet Scalars, JHEP 05, 193 (2021). arXiv:2101.00037.

[5] CMS Collaboration, Search for nonresonant triple Higgs boson production in the final state with six bottom quarks in proton-proton collisions at √s = 13 TeV (2026). arXiv:2607.05145.

[6] A. Papaefstathiou and G. Tetlalmatzi-Xolocotzi, Quadruple-Higgs boson production at the high-luminosity LHC (2026). arXiv:2608.16998.

[7] ATLAS Collaboration, HL-LHC prospects for the measurement of triple-Higgs production in the 6b final state at the ATLAS experiment (2025). ATL-PHYS-PUB-2025-003.

[8] S. Roy and C. E. M. Wagner, Alignment and Enhanced Multi-Higgs Production (2026). arXiv:2605.29345.