Search for exotic Higgs boson decays H $\to$ $\mathcal{A}\mathcal{A}$ $\to$ 4$\gamma$ with events containing two merged diphotons in proton-proton collisions at $\sqrt{s}$ = 13 TeV

We present the first direct search for exotic Higgs boson decays H $\to$ $\mathcal{A}\mathcal{A}$, $\mathcal{A}$ $\to$ $\gamma\gamma$ in events with two photonlike objects. The hypothetical particle $\mathcal{A}$ is a low-mass spin-0 particle decaying promptly to a merged diphoton reconstructed as a single photonlike object. We analyze the data collected by the CMS experiment at $\sqrt{s}$ = 13 TeV corresponding to an integrated luminosity of 136 fb$^{-1}$. No excess above the estimated background is found. We set upper limits on the branching fraction $\mathcal{B}$(H $\to$ $\mathcal{A}\mathcal{A}$ $\to$ 4$\gamma$) of (0.9-3.3) $\times$ 10$^{-3}$ at 95% confidence level for masses of $\mathcal{A}$ in the range 0.1-1.2 GeV.


1
Despite current constraints on the Higgs boson coupling to standard model (SM) particles from experiments at the CERN LHC [1,2], the Higgs sector remains an important area for physics searches beyond the SM (BSM).This is possible because the Higgs sector can potentially access BSM states that do not directly participate in the SM gauge interactions, called SM neutral.The simplest extensions of the Higgs sector include exotic decays of the form H → AA, where H is the 125 GeV boson with a non-SM decay, and A is a hypothetical new spin-0 particle decaying to a pair of SM particles [3].Such decays are found in BSM models containing an additional SMneutral singlet [3,4] and are of interest in searches for axionlike particle (ALP) production [5][6][7], including analyses performed at CMS [8][9][10][11][12].In ALP models, A is identified as a CP-odd spin-0 particle where CP is the charge-conjugation and parity operator.The experimental search presented in this Letter is insensitive to the CP quantum numbers of A, since its polarization is not measured.
For particle A masses m A ≲ 1 GeV, decays to many of the heavier fermions become inaccessible and the diphoton decay mode A → γγ becomes particularly relevant [3,13].Generic A → γγ searches are prevalent in astrophysics, cosmology [14][15][16], and particle collider experiments [17,18] because of the potential impact of A on stellar formation and the evolution of the early Universe.The particle A is also a potential dark matter candidate [19][20][21][22].
A measurement of the diphoton invariant mass spectrum at the LHC in this A mass regime is challenging if A is produced from the decay H → AA.In the A mass range of m A < 0.4 GeV, A will have a Lorentz boost γ L ≳ 150, where γ L = E A /m A and E A is the A energy.In the CMS experiment, this boost corresponds to a distance between the photons from the A → γγ decay at the electromagnetic calorimeter (ECAL) that is equal to or less than the Molière radius (r M ) of the ECAL material.The two photons will be predominantly reconstructed as one photonlike object (labeled Γ) by the standard CMS photon reconstruction algorithm, as described in Ref. [23].If A decays promptly, the H → AA → 4γ decay will lead to a measured two-photon invariant mass m ΓΓ peak degenerate with that of the SM H → γγ decay.The H → AA → 4γ signal will then be hidden by reconstructed SM H → γγ events [1,24].With A decays to more-massive particles being inaccessible, A can also be long-lived [3].
To separate a possible signal from the SM background, we present the first search that directly measures the invariant mass spectrum of merged diphoton candidates Γ reconstructed in events resembling an SM H → γγ final state.We assume that each A in the H → AA decay has the same mass and decays promptly to photons, though we provide estimated results where these assumptions are relaxed.The direct search is made possible by a novel particle reconstruction technique that is able to measure the invariant mass, m Γ , of merged diphoton candidates, something not possible with standard CMS reconstruction software.The technique utilizes deep learning algorithms trained directly on ECAL energy deposits to estimate m A in a so-called end-to-end m Γ reconstruction [23].It is applied for the first time at CMS to probe masses from 0.1 GeV (corresponding to γ L ≈ 600 with the diphotons collimated enough to enter the same ECAL crystal), to 1.2 GeV (corresponding to γ L ≈ 50 and diphotons separated by around 3r M at the ECAL).Tabulated results for this analysis can be found in [25].Signal events with larger m A have at least one of the A → γγ decays reconstructed as two distinct photons and are not studied.We note that the technique requires access to the full CMS event content which is not readily accessible because of storage constraints.No special event triggers are otherwise required.
The A → γγ decay is primarily detected in the ECAL, which contains a barrel section covering the pseudorapidity range |η| < 1.48 and two endcap sections on either end extending the range to |η| < 3. The crystal r M is comparable to the barrel crystal width of 2.2 cm, resulting in more than 90% of the energy of photons converting only in the ECAL barrel to be laterally contained within a 3×3 crystal matrix [31].
A sample of events containing two reconstructed photonlike objects (ΓΓ) is selected for analysis.We assume that each reconstructed photonlike object Γ corresponds to a single merged diphoton candidate.Events selected for study are required to pass trigger, photon reconstruction, and kinematic selection requirements similar to those used in the CMS SM H → γγ analysis [24].Key differences in the photon reconstruction criteria are described in the following paragraph.The SM background processes contributing to the selected sample are composed of nonresonant quantum chromodynamic (QCD) dijet and γ + jet production, prompt-diphoton production, and H → γγ decays.
Since the m A range we study spans a wide range of A boosts, the opening angles of the A → γγ decays vary as well [23].We thus apply a looser requirement on the output of the multivariate photon tagger [24], which then accepts more photons than for the CMS H → γγ search.This increases the contribution from QCD background processes containing hadronic jets, which may be reconstructed as photonlike objects if they radiate an energetic photon or contain one or more neutral-meson decays (π 0 /η → γγ).To compensate, we set a more restrictive requirement on the transverse momenta of charged particles in a cone around the selected photon, I ch [24].Since neutral mesons are mainly produced inside hadronic jets, they are more likely to be accompanied by charged particle tracks.Thresholds for the above requirements are chosen to maximize the significance of a possible H → AA → 4γ signal over the m A range of interest.After all selection criteria are applied, signal events are selected with an estimated efficiency of 8%-24%, decreasing with m A , while background events are reduced by >99%, as determined from simulation.
To simplify the application of the new m Γ reconstruction technique, only events with both ΓΓ reconstructed in the barrel section of the ECAL are analyzed.These account for about two thirds of the total expected signal yield.Events with more than two selected photons passing the selection criteria are not used to maintain orthogonality with a complementary CMS search at higher m A where the A → γγ photons are fully resolved [41].
To discriminate H → AA → 4γ events from background, the end-to-end m Γ reconstruction technique is used to measure the two-dimensional (2D) invariant mass spectrum, m Γ 1 versus m Γ 2 (labeled 2D-m Γ ), where Γ 1(2) is the higher-(lower-) energy reconstructed photon.Each Γ is assumed to correspond to a single A → γγ leg of the presumed H → AA → 4γ decay.We construct signal (S) 2D-m Γ templates for each signal mass hypothesis m A = 0.1-1.2GeV, in 0.1 GeV steps, using simulation.This defines our signal model at each m A hypothesis.The relative m Γ resolution varies from approximately 100% to 20% for m A = 0.1-1.2GeV [23].Background (B) 2D-m Γ templates are constructed for the SM processes contributing to our selected sample, derived from data and simulation.The sum of these background templates defines our background model, representing the SM-only hypothesis.A scan over different signal m A hypotheses is performed.At each point, maximum likelihood estimation (MLE) [42] is used to extract the best fit between the observed 2D-m Γ spectrum and the signal-plus-background model, µS + B, for some signal strength parameter µ.The detection of a potential H → AA → 4γ contribution at a given m A is expressed in terms of the significance of the extracted best fit signal strength.
Simulated H → AA → 4γ events are generated with MADGRAPH MC@NLO [43] at leading order using the SM + dark vector + dark Higgs model [3,4] and a Higgs boson mass of 125 GeV.A rounded mass value of 125 GeV is chosen as this search is insensitive to the exact value of the Higgs boson mass.The Higgs boson is produced with up to one associated jet, and the A has a negligible lifetime.The events are fully simulated in the CMS detector and account for effects from the underlying event and additional pp collisions in the same or nearby bunch crossings.Signal samples are generated separately for each m A hypothesis.
The simulated H → AA → 4γ events are used to construct the signal 2D-m Γ templates.Theoretical estimates of the H → AA production cross section vary depending on the model assumptions.Our final results, however, are expressed in terms of the signal branching fraction B(H → AA → 4γ), which is independent of the choice of cross section used to normalize the signal templates.
Background 2D-m Γ templates are derived to account for the H → γγ contribution and all other nonresonant processes.The SM H → γγ 2D-m Γ template is found from simulation.The H → γγ events are generated with MADGRAPH MC@NLO [43] at next-to-leading order, with up to two extra jets in the matrix element calculation using FXFX merging [44].The generated events are fully simulated in the CMS detector in the same way that the signal events described earlier are simulated.The normalization of the template is determined by fixing the template yield to the integrated luminosity of the data sample times the SM theoretical prediction for the total inclusive H production cross section [45] times the H → γγ branching fraction and efficiency to pass selection criteria.
An overall background 2D-m Γ template for all other nonresonant processes is estimated from data.We divide the selected events in data into regions based on m ΓΓ and 2D-m Γ .Since the two A particles come from the decay of a Higgs boson with a mass of 125 GeV, we define an m ΓΓ signal region (SR) around the Higgs boson mass (labeled m H -SR) using 110 < m ΓΓ < 140 GeV and sideband (SB) regions below (labeled m H -SB low ) and above (labeled m H -SB high ) the Higgs boson mass, by requiring 100 < m ΓΓ < 110 GeV and 140 < m ΓΓ < 180 GeV, respectively.Since we assume that the two A particles from the H → AA decay have equal mass, we also divide the 2D-m Γ distribution into an SR along the diagonal (labeled m A -SR), by requiring 3 GeV.The final SR in which the MLE test is performed corresponds to the intersection of the above two SRs, m A -SR ∩ m H -SR. A binned likelihood fit of µS + B versus the observed 2D-m Γ distribution with bin sizes of 0.05 × 0.05 GeV 2 is used.The boundaries between the SR and SB regions are tuned using the simulation to maximize the signal enrichment in the SR, while maintaining a sufficient number of events in the SB regions to give a good background estimation.Figure 1 (center) shows the observed 2D-m Γ distribution of the selected events in data from m H -SR and the respective 1D projections of m Γ 1 (left) and m Γ 2 (right) in the m A -SR. Illustrative contours in 2D-m Γ of the signal template (center) and its 1D projections (left and right) are also provided.
An estimate of the nonresonant background shape of the 2D-m Γ template in the m H -SR is obtained from data SBs by taking the event-weighted average of the shapes from the two m H -SB  The statistical uncertainties in the former are negligible and the total uncertainties in the latter are barely visible as green bands.The spectra of simulated H → AA → 4γ events for m A = 0.1 (purple dashed curve), 0.4 (gray dotted curve), and 1.0 GeV (orange dash-dotted curve) are also provided.They are each normalized to the value of B(H → AA → 4γ) that is expected to be excluded by the background model (described under the CL s criterion in our results) times 10 3 .The black points in the lower panels of the left and right plots give the ratios of the data to the predicted background distributions.The vertical bars represent the statistical uncertainties in the former, and the green bands represent the total uncertainty in the latter.
regions (together labeled m H -SB). To determine the normalization of the template, we assume that the ratio of the number of events along the m A -SR diagonal, from that observed in the m H -SRs to that estimated from the m H -SB, N(m A -SR ∩ m H -SR)/N(m A -SR ∩ m H -SB), is the same as for the m A -SB off diagonal, N(m A -SB ∩ m H -SR)/N(m A -SB ∩ m H -SB).This is justified since the ratio of the number of events in the m H -SR to the m H -SB as a function of 2D-m Γ is estimated to be constant after the reweighting procedure described below.The normalization for the template shape is then determined by solving for the expression for N(m A -SR ∩ m H -SR).
For hadronic jets passing the event selection, there is an energy dependence, and thus an m ΓΓ dependence, that violates the assumption described in the previous paragraph.With increasing jet energy, more energy becomes available to produce additional hadrons in the jet, which increases the effective mass of the jet and thus its reconstructed m Γ .To correct for this effect, prior to deriving the 2D-m Γ shape, events from the two m H -SB regions are first reweighted so that their transverse momentum distributions match that in the target m H -SR.
The normalized H → γγ and nonresonant background templates are then added together.Their combined yield is renormalized so that the ratio of the predicted to the observed yields in the m A -SB region remains unity.The H → γγ template accounts for about 0.4% of the total background template yield.To account for residual differences in the 2D-m Γ shape between the background estimate and the observed distribution in the m A -SB region, the total background template is multiplied by a 2D polynomial function pol(m Γ 1 , m Γ 2 ) over the full 2D-m Γ range.A linear polynomial is used with parameters chosen to maximize the likelihood with respect to the observed 2D-m Γ shape in the m A -SB region.No further improvement in the goodness of fit is found with higher polynomial orders.The resulting corrected 2D-m Γ background template defines the background model used in the MLE fit.The left and right plots in Fig. 1 show that the SM background in the observed 2D-m Γ distribution is dominated by single photonlike objects, which exhibit a smoothly falling m Γ spectrum, rather than by neutral-meson decays, which would be manifested as peaks [23].Neutral-meson decays from γ + jet production are more likely to be reconstructed in the lower-energy Γ 2 .The background estimation is validated using an orthogonal data sample obtained by inverting the I ch requirement on Γ 2 , to ensure negligible signal contamination.No significant bias is observed when the signal extraction procedure is performed on this sample.
Uncertainties in the predicted signal and background templates are treated as nuisance parameters in the MLE procedure used to determine the best fit signal strength μ.The dominant uncertainties impacting the extracted μ are those from statistical uncertainties in the background template's shape.Their largest impact on the relative uncertainty in μ varies between 15% and 20%, depending on m A .For m A = 0.1 GeV, where the m Γ resolution is poorest and the background contribution is largest, systematic uncertainties affecting the background template normalization are also important.These include systematic uncertainties associated with the best fit parameters of pol(m Γ 1 , m Γ 2 ) and the relative contribution of m H -SB low versus m H -SB high events in the nonresonant background 2D-m Γ template.These systematic uncertainties impact the relative uncertainty in μ by about 25% (10%) for m A = 0.1 (1.2) GeV.
The systematic uncertainty in the signal strength from using the 2D-m Γ template determined from simulation is estimated from a sample of electrons in events with Z → e + e − decays, selected from both data and simulation with the "tag-and-probe" technique [46].Electrons are preferred over decays of neutral mesons in jets because of the complicating effects of accompanying hadrons in the same jet [23].Specifically, uncertainties associated with a relative m Γ scale shift and an increase in the smearing of the mass peak are estimated and found to have a negligible impact on μ.
The best fit background estimate determined from the MLE procedure is shown by the blue solid curves in Fig. 1 (left and right), together with its associated best fit total statistical plus systematic (stat + syst) uncertainties.We find no statistically significant excess in the data over the SM background predictions for m A masses in the range 0.1-1.2GeV.
The CL s criterion [47,48] is used to interpret this result in terms of excluded B(H → AA → 4γ) values.The observed upper limit on B(H → AA → 4γ) at 95% confidence level (CL) as a function of m A in the range 0.1-1.2GeV is shown in Fig. 2, and varies between (0.9-3.3)×10 −3 for m A values 0.1-1.2GeV.The expected 95% CL limits and their associated 68 and 95% confidence intervals (CIs) are determined by simulating SM background-only pseudo-experiments.
The LHC measurements of B(H → γγ) [1,2] give an effective upper bound on a possible measurement of B(H → AA → 4γ) because of the degeneracy between the final states.The constraint from the CMS measurement [1] is shown in Fig. 2. It is relevant for values of m A ≈ 0.1 GeV where the A → γγ decay resembles a single photon and increases at larger m A .Our observed upper limits thus set the best constraints for this decay mode in the m A range that we study.
We estimate the upper limits for long-lived A decays by comparing the signal yield in the m A -SR ∩ m H -SR for different simulated A decay lengths compared with that for prompt decays.For m A = 0.1 (0.4) GeV, the 95% CL upper limit on B(H → AA → 4γ) is 1.6 (0.9) times the prompt-decay upper limit for an A decay length of cτ 0 = 1 mm, and 30 (3) times larger for cτ 0 = 10 mm, with a linear interpolation between the two limits in both cases.Better upper limits result when the increased merging from long-lived A decays improves photon reconstruction more than it degrades m Γ resolution.The prompt-decay upper limits are also relevant for models with dissimilar A masses, H → A 1 A 2 , with m A 1 ̸ = m A 2 , for mass differences less than the m A -SR window, |m A 1 − m A 2 | ≲ 0.3 GeV.For larger mass differences, the signal mass peak would fall outside of the m A -SR and be absorbed into the m A -SB, making a measurement impossible.
In summary, the results of a search for the exotic Higgs boson decay H → AA → 4γ for a Higgs boson mass of 125 GeV are presented.Events reconstructed with two photonlike objects are used, where each photonlike object is assumed to be a merged A → γγ candidate.A method is developed to measure the invariant mass of merged diphoton candidates to discriminate a potential signal from the standard model background.This is the first search of its kind at CMS and made possible by a novel end-to-end reconstruction technique of merged diphotons.No excess of events above the estimated background is found.An upper limit on the branching fraction B(H → AA → 4γ) of (0.9-3.3)×10 −3 is set at 95% C.L. for masses of A in the range 0.1-1.2GeV, assuming prompt A decays.These are the current best constraints on H → AA → 4γ in this m A range.
Tabulated results for this analysis are provided in the HEPData record [25].

Figure 1 :
Figure 1: Mass distributions from selected events in data.Center: the 2D-m Γ distribution for data events in the m H -SR. The red dashed lines indicate the m A -SR boundaries.The contours of simulated H → AA → 4γ events for m A = 0.4 GeV are plotted for 75% (solid contour) and 50% (dotted contour) of the distribution maximum.The corresponding m Γ 1 (left) and m Γ 2 (right) projections for the overlap of the m H -SR and m A -SR are also shown.The data distributions (black points) are plotted against the total predicted background distributions (blue curves).The statistical uncertainties in the former are negligible and the total uncertainties in the latter are barely visible as green bands.The spectra of simulated H → AA → 4γ events for m A = 0.1 (purple dashed curve), 0.4 (gray dotted curve), and 1.0 GeV (orange dash-dotted curve) are also provided.They are each normalized to the value of B(H → AA → 4γ) that is expected to be excluded by the background model (described under the CL s criterion in our results) times 10 3 .The black points in the lower panels of the left and right plots give the ratios of the data to the predicted background distributions.The vertical bars represent the statistical uncertainties in the former, and the green bands represent the total uncertainty in the latter.

Figure 2 :
Figure2: Observed (black solid curve with points) and median expected (blue dashed curve) 95% CL upper limit on B(H → AA → 4γ) as a function of m A for prompt A decays.The 68% (green band) and 95% (yellow band) CIs are plotted around the expected limit.The 95% CL upper limit from the CMS measurement[1] of B(H → γγ) is also shown (red band, where the width represents the uncertainty in the measurement).