Search for Higgs Boson Decays into a Z Boson and a Light Hadronically Decaying Resonance Using 13 TeV pp Collision Data from the ATLAS Detector

A search for Higgs boson decays into a Z boson and a light resonance in two-lepton plus jet events is performed, using a pp collision dataset with an integrated luminosity of 139 fb − 1 collected at ﬃﬃﬃ s p ¼ 13 TeV by the ATLAS experiment at the CERN LHC. The resonance considered is a light boson with a mass below 4 GeV from a possible extended scalar sector or a charmonium state. Multivariate discriminants are used for the event selection and for evaluating the mass of the light resonance. No excess of events above the expected background is found. Observed (expected) 95% confidence-level upper limits are set on the Higgs boson production cross section times branching fraction to a Z boson and the signal resonance, with values in the range 17 – 340 pb ( 16 þ 6 − 5 – 320 þ 130 − 90 pb) for the different light spin-0 boson mass and branching fraction hypotheses, and with values of 110 and 100 pb ( 100 þ 40 − 30 and 100 þ 40 − 30 pb) for the η c and J= ψ hypotheses, respectively.

The structure of the standard model (SM) scalar sector is the subject of intense scrutiny by the ATLAS [1] and CMS [2] Collaborations at the CERN Large Hadron Collider (LHC) [3]. At the current level of precision, all of the measured properties of the Higgs boson (H) [4,5] are found to be consistent with their SM predictions [6][7][8][9][10], and no additional Higgs boson has been observed to date. However, given the small natural decay width of the Higgs boson, even small additional contributions from physics beyond the SM can lead to final states with substantial, and thus possibly detectable, branching fractions (B) [11]. This Letter presents a search for Higgs boson decays into a Z boson and a hadronically decaying light resonance in events with a same-flavor lepton pair (electrons or muons) and a jet in the ATLAS detector. Hadronic decays of an η c or of a J=ψ charmonium resonance (Q), or of a light spin-0 boson from an extended Higgs sector with a mass up to 4 GeV, are considered and are reconstructed as a single jet.
The Yukawa sector of the SM [12] does not provide an explanation for the observed fermion mass hierarchy. As a result, a wide range of new physics scenarios have been proposed, including the Froggatt-Nielsen mechanism [13] and the ; for a recent overview, see Ref. [15]. The couplings of the Higgs boson to the third-generation fermions [16][17][18][19][20][21] have been observed, and a program to probe its couplings to the first-and second-generation charged leptons has been established [22][23][24][25]. For its couplings to first-and secondgeneration quarks, several approaches are being explored. Focusing on the Higgs boson's coupling to the charm quark, direct searches have been performed for Higgs boson decays into charm quarks [26,27] and for exclusive decays into a J=ψ and a photon [28,29], with no excess observed. Constraints from differential cross section measurements of Higgs boson production versus transverse momentum (p T ) have also been derived [30,31]. Higgs boson decays into a gauge boson and a charmonium state, including an η c or a J=ψ, have been proposed as another way to access the coupling of the Higgs boson to the charm quark [32][33][34] and to probe the nature of the Higgs boson [35]. This search follows the last approach and maximizes the signal acceptance by focusing on inclusive hadronic final states of the mesons in H → Zη c and H → ZJ=ψ decays, which have SM branching fractions of 1.4 × 10 −5 and 2.2 × 10 −6 [35], respectively.
While the SM posits a single complex Higgs doublet field [36,37], extended Higgs sectors are motivated [38] and provide a rich phenomenology of additional scalars. Two such models discussed here are the two-Higgs-doublet model (2HDM) [11,39] and the 2HDM with an additional scalar singlet (2HDM þ S) [11,40]. These represent two of the simplest extensions of the scalar sector, and with their type-II fermion couplings they are necessary to generate the masses in the minimal supersymmetric SM and the next-tominimal supersymmetric SM, respectively [41]. Both of these models can include additional light pseudoscalars (a) with significant BðH → ZaÞ or BðH → aaÞ [11]. In the 2HDMðþSÞ, these two B values can be adjusted independently, therefore searches for H → aa do not constrain BðH → ZaÞ, so that searches for the latter decay are required [11,34]. Despite the Yukawa nature of the a to fermion couplings, there are large regions of parameter space depending on the mass of a and the ratio of the vacuum expectation values of the two Higgs-doublet fields (tan β) [11], where these pseudoscalars decay mainly to gluons and light up-type quarks, as the decays into downtype fermions are suppressed. These experimental signatures are also relevant in axion models [42][43][44], models of electroweak baryogenesis [45], neutrino mass models [46], dark-matter models [46,47], and models of grand unification [48]. Previous searches for Higgs boson decays into light scalars have been performed at the Tevatron [49] and the LHC [50][51][52][53][54][55][56][57][58][59]. However, these were mostly focused on searches for H → aa, in final states including leptons, photons, or bottom quarks. By targeting the H → Za, a → hadrons decay channel, this search accesses new, previously unexplored regions of the parameter space.
Searches for hadronic decays of light resonances are challenging at the LHC due to the large multijet background. However, substantial progress has been made in the use of jet substructure techniques in boosted final states [60], typically in searches or measurements involving heavy resonances [61,62]. In this Letter, jet substructure variables enable the reconstruction of a light, boosted, hadronic final state. Information from the individual substructure variables is combined using machine learning techniques. Specifically, for event selection, a multilayer perceptron (MLP) [63] classifier is employed. Given the range of masses considered, the classifier is provided with resonance-mass-related information from a separate MLP-based mass estimator, which results in improved classification performance over the full mass range.
This search is performed using the complete run 2 pp collision dataset, produced between 2015 and 2018 at a center-of-mass energy ffiffi ffi s p ¼ 13 TeV by the LHC. The data were collected by the ATLAS detector [1] and correspond to an integrated luminosity of 139 fb −1 .
Monte Carlo (MC) samples of simulated events are used to model the signal selection efficiency. The signal samples were generated via the gluon-gluon fusion process using POWHEG-BOX v2 [64][65][66], with the CT10 next-to-leading order (NLO) parton distribution function (PDF) set [67]. Particle decays, hadronization, parton showers, and the underlying event were modeled using PYTHIA v8.212 [68] and Evt Gen v1.6.0 [69], interfaced to the AZNLO [70] set of tuned parameters and the CTEQ6 L1 PDF set [71]. Next-tonext-to-leading order (NNLO) corrections are applied to the p T distribution of the Higgs boson. The a branching fractions were determined using PYTHIA 8 [68] with a 2HDM tan β value of 1, which predicts a → gg to be the dominant decay mode until a → cc becomes kinematically accessible. The signal MC samples used in this analysis have a masses of 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, and 4 GeV. The Z boson is required to decay into pairs of electrons, muons, or τ leptons.
The background is dominated by Z þ jets events, modeled using SHERPA 2.2.1 [72] interfaced to the NNPDF 3.0 (NNLO) PDF set [73]. The inclusive production cross sections are known to NNLO in QCD [74]. The ZZ, ZW, and tt processes contribute < 1% of the total background in this search. The diboson backgrounds were modeled using SHERPA 2.2.1 interfaced to the NNPDF 3.0 (NNLO) PDF set, except for gluon-induced ZZ production, which was modeled using SHERPA 2.2.2 [72]. All of the SHERPA samples used a set of tuned parameters developed by the SHERPA authors. The tt process was modeled using POWHEG-BOX v2, while the subsequent decay, hadronization, parton shower, and underlying event were modeled using PYTHIA v8.230 and EvtGen v1.6.0. The  Events are selected by a combination of single electron or muon triggers for each data-taking period [82][83][84][85], and the online lepton reconstructed by the trigger is required to be within ΔR ¼ 0.1 [86] of an off-line reconstructed lepton. Events are required to have at least one reconstructed primary interaction vertex [87]. Electron candidates are reconstructed by matching tracks in the inner detector to topological energy clusters in the electromagnetic calorimeter [80] and must pass a likelihoodbased selection, which requires the shower profile to be compatible with that of an electromagnetic shower. Muons are reconstructed using tracks in the muon spectrometer, matched to tracks in the inner detector where available [88]. Electrons and muons are each required to have p T > 18 GeV, and at least one must have p T > 27 GeV. Electrons (muons) are required to be reconstructed within jηj < 2.47 (jηj < 2.7), but electrons within 1.37 < jηj < 1.52 are excluded. The transverse energy sum in a cone of size ΔR ¼ 0.2 around the electron [muon] in the calorimeter must be less than 20% (30%) of the lepton's p T , and the summed p T of tracks within a cone of variable size ΔR ¼ minð0.2; 10 GeV=p T Þ [ΔR ¼ minð0.15; 10 GeV=p T Þ] around the electron [muon] must be less than 15% of its p T . Contributions from nearby electrons and muons are removed from these cones. If an inner detector track is present, muons must also have a longitudinal impact parameter jz 0 sin θj < 0.5 mm and a transverse impact parameter jd 0 j < 1 mm relative to the primary interaction vertex. At least two same-flavor opposite-sign electrons or muons are required to pass this selection and have an invariant mass compatible with the mass of the Z boson: 81 < m ll < 101 GeV. If multiple same-flavor opposite-sign lepton pairs fulfill this requirement, the pairing with an invariant mass closest to that of the Z boson is chosen. Z → ττ decays are reconstructed through the leptonic decays of the τ leptons.
The hadronically decaying resonance is reconstructed as a single jet using the anti-k t jet algorithm [89,90] with a radius parameter of 0.4, formed from topological calorimeter energy clusters [91,92] and calibrated to the electromagnetic energy scale. Jet energies are corrected for contributions from simultaneous inelastic pp interactions (pileup) using a jet-area-based technique [93,94] and calibrated [95,96] using p T -and η-dependent correction factors determined from simulation, with residual corrections from in situ measurements applied to data and internal jet properties. Jets are required to have p T > 20 GeV and jηj < 2.5 and satisfy a jet cleaning requirement [97]. To reject jets from pileup interactions, jets with p T < 60 GeV and jηj < 2.4 are required to pass a "jet vertex tagger" [79] requirement. An overlap removal procedure resolves cases in which multiple electrons, muons, or jets are reconstructed from the same detector signature. Higgs boson candidates are reconstructed from the lepton pair and jet system, which is required to have an invariant mass passing a loose preselection requirement: m llj < 250 GeV. If multiple jets satisfy these requirements, the jet with the highest p T is selected. The acceptance for this preselection, evaluated using generator-level MC samples, varies between 28% and 29% for the different Q=a signal hypotheses.
MLPs [63] are used to select signal events passing this preselection. The MLP input variables are built using tracks matched to the calorimeter jet by ghost association [93], in which the tracks are included in the jet clustering process as with negligible energy and their angles from the jet axis. This allows the MLP to benefit from the high resolution of the tracking detector. These tracks must have p T > 500 MeV and jηj < 2.5 and pass loose quality and track-tovertex association requirements [98] to reject fake tracks from the reconstruction and tracks from pileup, respectively. Six dimensionless variables are constructed using these tracks: the ratio of the p T of the highest p T track to the p T of the ghost-associated track system; the angular separation ΔR between the highest-p T track and the calorimeter jet axis; NSubJettiness 2 [99], using exclusive-k t subjet axes with radius parameters of 0.2, and a jet axis radius parameter of 0.4; angularity (2) [100]; and U 1 ð0.7Þ and M 2 ð0.3Þ, which are modified energy correlation functions [101] designed for quark-gluon discrimination and to target two-pronged substructure, respectively. These variables primarily capitalize on the presence of a narrow resonance or two-pronged substructure in the track system. Initially, a regression MLP [63], using four hidden layers of 12 nodes, is trained using the above input variables and the a signal samples to estimate the mass of a, as shown in Fig. 1 FIG. 1. Output of (a) the regression and (b) the classification MLPs, for data, background, and three signal hypotheses. Events are required to pass the complete event selection, including the 120 < m llj < 135 GeV requirement, but not the requirement on the classification MLP output variable. The background normalization is set equal to that of the data, and the signal normalizations assume the SM Higgs boson inclusive production cross section and BðH → ZaÞ ¼ 100%, and in (a) the signal normalization is scaled up by a factor of 100. The error bars (hatched regions) represent the data (MC) sample statistical uncertainty, in both the histograms and the ratio plots. In (b) the region to the right of the dashed line is the signal region.

(a). This estimated mass is then
provided alongside the six input variables to a classification MLP [63], to inform the classifier about the part of the hadronic resonance mass spectrum where the specific event lies. This classification MLP has two hidden layers of six and five nodes and is trained using the a signal samples and the background samples. The 0.75 GeV a signal sample is excluded from the training of the classification MLP to ensure an even spacing between the a mass hypotheses, so the training is not biased toward lower masses. Both MLPs use sigmoidal response functions with summed inputs and are trained using backpropagation with a mean-square estimator [63], as these resulted in optimal discrimination without overtraining. The addition of the regression MLP was found to result in about a 13% improvement in the S= ffiffiffi ffi B p of the classification MLP, where S and B are the expected numbers of signal and background events passing the MLP requirement, respectively. The classification MLP output variable (M) is shown in Fig. 1(b).
The signal region (SR) for this search is defined by the requirements 120 < m llj < 135 GeV and M > 0.0524, chosen to maximize the expected S= ffiffiffi ffi B p , averaged over the various a mass hypotheses. The efficiency of this MLP requirement for events passing the preselection is ð0.761 AE 0.020Þ% for the background, ð5.89 AE 0.24Þ% and ð6.66 AE 0.26Þ% for H → Zη c and H → ZJ=ψ, respectively, and between ð1.88 AE 0.15Þ% and ð45.9 AE 0.8Þ% for H → Za. The efficiencies for the complete selection are estimated using MC samples and are ð0.545 AE 0.022Þ% and ð0.560 AE 0.022Þ% for H → Zη c and H → ZJ=ψ, respectively, and range between ð0.140 AE 0.011Þ% and ð3.27 AE 0.06Þ% for H → Za. The efficiencies are highest for the lowest a mass hypotheses, due to higher probabilities to pass the MLP requirement. The efficiency for H → Zη c events to pass the MLP requirement is lower than that of H → ZJ=ψ events, as J=ψ decays tend to have a lower charged hadron multiplicity. Using the predicted cross section for inclusive SM Higgs boson production of 55.7 þ3.0 −3.9 pb [102], and B½H → ZðQ=aÞ ¼ 100%, gives expected signal yields of 4260 and 4370 for H → Zη c and H → ZJ=ψ, respectively, and between 1090 and 25600 for H → Za.
A "modified ABCD estimate" of the total background in the SR is derived using four regions: A, defined by 0.0341 < M < 0.0524, expected to contain about 10% of the total background, and 155 < m llj < 175 GeV; B, defined by the m llj requirement of the SR and the M requirement of region A; C, defined by the M requirement of the SR and the m llj requirement of A; and D, which is the SR. An initial data-driven background estimate in the SR is calculated as D ¼ BC=A, then MC samples, reweighted to match data, are used to correct this estimate for the 13% correlation between the m llj and M variables. This reweighting is performed in the p T of the calorimeter jet, the number of ghost-associated tracks and U 1 ð0.7Þ. This background estimate is 82400 AE 2900 events in the SR, where the uncertainty is due to the limited data and MC sample statistics. The background estimation method is found to be consistent with data within 1.7 times the total statistical and systematic uncertainty in 14 validation regions, defined in regions of the m llj and M variables.
A measure of σðpp → HÞB½H → ZðQ=aÞ is extracted for a given signal hypothesis using a maximum-likelihood fit [103] to the number of events observed in the SR. The systematic uncertainties are included in the likelihood fit as nuisance parameters, which modify the signal efficiencies or the simulation-based correction used to calculate the expected background yield. These systematic uncertainties include uncertainties in the signal and background modeling and experimental uncertainties. The sources of modeling uncertainty include the limited MC sample statistics, renormalization scale and choice of MC generator for the signal and background, and a signal uncertainty to account for the extrapolation from gluon-gluon fusion signal samples to the inclusive Higgs boson production cross section. The effects of factorization scale and PDF uncertainties are found to be negligible. reconstruction. The total uncertainty on the extracted signal yield is dominated by the background modeling uncertainties, the largest being due to limited MC sample statistics. The total uncertainty on the background in the SR is 3700 events, where the uncertainty due to the limited data and MC sample statistics is 2900 and the modeling uncertainty is 2300. The data statistical uncertainty corresponds to approximately 8% of the total uncertainty on the extracted signal yield.
The SR contains 82 908 data events. This result is compatible with the SM background-only expectation, and the three-body mass distribution is shown in Fig. 2 Table I. In the absence of systematic uncertainties, these limits would range between 1.9 and 55 pb for the different signal hypotheses. To simplify the interpretation, the upper limits are quoted for Bða → ggÞ ¼ 100% and Bða → ssÞ ¼ 100%. Because of the Yukawa ordering of the decays of Higgs bosons, only decays into gluon and strange quark pairs are considered. The tighter limits for the a → ss decays are due to a higher MLP selection efficiency. The systematic uncertainties for a → gg and a → ss decay hypotheses are estimated using the inclusive decays as modeled in PYTHIA 8, which is a good approximation due to the dominance of the background modeling uncertainties. This assumption allows a limit to be set on the decay of a into gg or ss final states, in any ratio, by using a weighted sum of the two limits. Higgs boson decays to a Z boson and a quarkonium state other than the considered signal process are not included in the statistical interpretation.
In conclusion, a search has been performed for Higgs boson decays into a Z boson and either a η c or J=ψ charmonium state, or a light spin-0 boson. No excess is found, and 95% CL upper limits are set on σðpp → HÞB½H → ZðQ=aÞ, with values of 110 and 100 pb for the H → Zη c and H → ZJ=ψ hypotheses, respectively, and with values in the range 17-340 pb for the H → Za signal hypotheses. Assuming the SM prediction for inclusive Higgs boson production, the limits on charmonium decay modes correspond to branching fraction limits in excess of 100%. This is the first direct limit on decays of the observed Higgs boson to light scalars, decaying to light quarks or gluons. Because of the large value of Bða → hadronsÞ over the entire 2HDMðþSÞ parameter space, these limits represent tight, direct constraints for low (high) tan β in the type-II and type-III (