Search for the decay of the Higgs boson to a 𝒁 boson and a light pseudoscalar particle decaying to two photons

A search for the decay of the Higgs boson to a 𝑍 boson and a light, pseudoscalar particle, 𝑎 , decaying respectively to two leptons and to two photons is reported. The search uses the full LHC Run 2 proton–proton collision data at √ 𝑠 = 13 TeV, corresponding to 139 fb − 1 collected by the ATLAS detector. This is one of the first searches for this specific decay mode of the Higgs boson, and it probes unexplored parameter space in models with axion-like particles (ALPs) and extended scalar sectors. The mass of the 𝑎 particle is assumed to be in the range 0.1–33 GeV. The data are analysed in two categories: a merged category where the photons from the 𝑎 decay are reconstructed in the ATLAS calorimeter as a single cluster, and a resolved category in which two separate photons are detected. The main background processes are from Standard Model 𝑍 boson production in association with photons or jets. The data are in agreement with the background predictions, and upper limits on the branching ratio of the Higgs boson decay to 𝑍𝑎 times the branching ratio 𝑎 → 𝛾𝛾 are derived at the 95% confidence level and they range from 0.08% to 2% depending on the mass of the 𝑎 particle. The results are also interpreted in the context of ALP models.


Introduction
Light pseudoscalars that couple to Higgs bosons appear in many well-motivated new physics scenarios.
A few examples of such scenarios are the next-to-minimal supersymmetric standard model [1], axion models [2,3], and models with dark or extended Higgs sectors [4][5][6].In addition, they may be related to a possible muon  − 2 deviation with respect to the Standard Model (SM) predicted value, e.g. as in Ref. [7].The search for elusive particles in the decays of the 125 GeV Higgs boson [8, 9] is a highly promising avenue for their discovery.Due to the narrow decay width of the Higgs boson [10], even minor contributions from physics beyond the SM can result in final states with significant branching fractions.The most recent combination of the ATLAS measurements of the Higgs boson production rates does not rule out decays to yet unobserved states with a branching ratio of up to 12% [11] at 95% confidence level (CL).The corresponding result from the CMS experiment is 16% [12].These results indicate that there may be sizeable exotic decays of the Higgs bosons that have not been observed yet.The LHC experiments have looked for such particles in Higgs boson decays, usually assuming that the Higgs boson, , decays to a pair of narrow-width, pseudoscalar particles  ( → ), which subsequently decay to fermions or photons.Such searches have been performed in the  [13, 14], 4 [15,16], 4 [17,18],  or 4 [19][20][21][22][23][24],  [25], 4 [26-28] and 2 + 2 jets [29] final states.In addition to the  →  channel, the decay mode of  →  is well-motivated, but remains less studied.Such decays are particularly motivated by axion models [2].This decay has been searched for by ATLAS in the case where the  boson decays to leptons and  decays to jets [30], but this is one of the first searches for the  →  with  → .In this letter, results from the search for the final state in which the  boson decays to leptons and the  particle decays to a photon pair are presented using the full LHC Run 2 proton-proton ( ) collision data collected by the ATLAS detector during 2015-2018 and corresponding to an integrated luminosity of 139 fb −1 .The search explores  particle masses,   , between 0.1 GeV and 33 GeV.When the invariant mass of the two photons from the  decay is below 2 GeV and the transverse momentum of the diphoton system is high, their electromagnetic showers typically overlap in the detector and a single photon is reconstructed.The search presented here considers both the merged and the resolved  →  decays.This letter is structured as follows.A description of the ATLAS detector is given in Section 2, followed by a presentation of the data and simulated samples in Section 3. The physics object reconstruction is discussed in Section 4, followed by an explanation of the analysis in Section 5, the presentation of the search results in Section 6 and the conclusion in Section 7.

ATLAS detector
The ATLAS detector [31] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer.The inner tracking detector covers the pseudorapidity range || < 2.5.It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. 1ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe.The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards.Cylindrical coordinates (, ) are used in the transverse plane,  being the azimuthal angle around the -axis.

Data and simulated event samples
The search for  →  is performed using the full LHC Run 2 (2015-2018) √  = 13 TeV   collision dataset recorded by the ATLAS detector and corresponding to an integrated luminosity of 139.0 fb −1 [34], which includes only data-taking periods where all relevant detector subsystems were operational [35].The data events are collected using a set of triggers that require the presence of one or two electrons [36] or one or two muons [37] with thresholds in the transverse energy for electrons and the transverse momentum for muons that are in the range 12-26 GeV, depending on the lepton flavour, the number of leptons and the data-taking period [38].Data-driven corrections are applied to the event-level trigger efficiencies.
The Higgs bosons that are generated undergo subsequent decay into a  boson and a pseudoscalar particle, , through Pythia 8.244.3 [49].The decay of the  boson is isotropic in its decay frame.The simulated samples are generated with the  particle mass in the range 0.1-33 GeV.
The main backgrounds for this search come from  boson production in association with jets (+jets) or a photon ( + ).The production of +jets is simulated with the Sherpa 2.2.1 [65] generator using NLO matrix elements for up to two partons, and leading-order (LO) matrix elements for up to four partons calculated with the Comix [66] and OpenLoops [67][68][69] libraries.They are matched with the Sherpa parton shower [70] using the MEPS@NLO prescription [71][72][73][74] using the set of tuned parameters developed by the Sherpa authors.The NNPDF3.0nnlo set of PDFs [75] is used and the samples are normalised to a NNLO prediction [76].The production of ( + ) final states is simulated with the Sherpa 2.2.8 generator.Matrix elements at NLO QCD accuracy for up to one additional parton and LO accuracy for up to three additional parton emissions are matched and merged with the Sherpa parton shower based on Catani-Seymour dipole factorisation [66,70] using the MEPS@NLO prescription.The virtual QCD corrections for matrix elements at NLO accuracy are provided by the OpenLoops 2 library [67][68][69]77].The generation uses the NNPDF3.0nnloPDF set [75], along with the dedicated set of tuned parton-shower parameters developed by the Sherpa authors.Finally, the Higgs boson decay  →  [78] is found to be negligible, and thus not included.
The effects of multiple   interactions in the same bunch crossing as the hard scatter, and in neighbouring bunch crossings, (pile-up) are included using simulated events generated with Pythia 8.186 [79] with the A3 minimum-bias tune.Simulated events are weighted to reproduce the distribution of the average number of interactions per bunch crossing observed in data.
All the generated events are passed through the simulation of the ATLAS detector [80] based on Geant4 [81] and the same algorithms as for data.

Object reconstruction
Electron candidates are created from an energy deposit in the EM calorimeter matched with a track in the inner detector.They are required to be within the fiducial region of the detector, || < 2.47, excluding the transition region between the barrel and endcap calorimeters (1.37 < || < 1.52).The associated track must have | 0 |/  0 < 5 and | 0 | sin  < 0.5 mm, where  0 ( 0 ) is the transverse (longitudinal) impact parameter relative to the primary vertex, and   0 is the uncertainty in  0 .The reconstructed electrons are required to have  T > 20 GeV and their quality is ensured by using a likelihood method with the requirement to satisfy "medium" identification criteria [82].The energy of the electron is calibrated through a series of successive steps, utilising a combination of simulation-based and data-driven correction factors [82].
Muon candidates are reconstructed in the range || < 2.5 by combining tracks in the inner detector with tracks in the muon spectrometer.Muon candidates must have | 0 |/  0 < 3 and | 0 | sin  < 0.5 mm.In addition they must have  T > 20 GeV and satisfy "medium" quality requirements [83].The efficiency of muon reconstruction and identification, along with the momentum calibration and associated systematic uncertainties, have been estimated using the methods described in Refs [83,84].
Both muon and electron candidates are required to satisfy isolation requirements based on tracking and calorimeter measurements in order to suppress hadronic and non-prompt lepton backgrounds [82,84].In addition, electrons reconstructed within Δ = 0.2 from a muon are not considered in the analysis.
Data-driven corrections are applied at the event level to ensure accurate reconstruction, identification and isolation efficiencies of muons and electrons.
Photon candidates are reconstructed from topological clusters of energy deposits in the EM calorimeter and calibrated as described in Ref. [82].Photon candidates are required to have  T > 10 GeV and || < 2.37, excluding the transition region between the barrel and endcap calorimeters.Backgrounds from jets are reduced by requiring them to fulfil "loose" identification criteria based on shower shapes in the ATLAS EM calorimeter [82].Finally, photon candidates are not considered if they are reconstructed within Δ < 0.3 from an electron or muon.

Analysis
The signature that is searched for consists of two high- T leptons from the  boson decay and a photon pair from the  particle decay.The boost of the photon pair determines whether it will be reconstructed as a single photon or as two resolved photons.This motivates a categorisation of events depending on whether the  pair is reconstructed as one ("merged" category) or two ("resolved" category) photons.The main sources of background for both categories is +jets and  +  production.
The events used in this search are required to have exactly two electrons or two muons (hereafter referred to as the two leptons), with one of them satisfying  T > 27 GeV and the other  T > 20 GeV.Subsequently, the two leptons are required to have opposite electric charge and their invariant mass must be in the range 81-101 GeV to be compatible with the mass of the  boson.Finally, the transverse momentum of the two-lepton system is required to be greater than 10 GeV, a criterion that is optimised to improve both background rejection and to reduce uncertainties in the background modelling.
The events are further required to contain at least one photon.If the event contains two or more photons, for all possible pairs for which the angular distance between the two photons satisfies Δ  < 1.5, the quantity  = Δ    T /(2  ) is calculated.In this variable,   T and   denote the transverse momentum and the invariant mass of the photon pair, respectively. 2If there is a photon pair with 0.96 <  < 1.2 then the event is included in the "resolved" category.In case of multiple photon pairs satisfying this requirement, the pair with the  value closest to unity is considered as the  →  pair candidate.If the event fails the resolved category requirement, then it becomes part of the "merged" category if it contains at least one photon with  T > 20 GeV.In this case, the highest- T photon is selected as the  →  decay candidate.The resolved category is sensitive to relatively high  particle masses, whereas the merged category is designed for the low mass regime with typically   ≲ 2 GeV.
After event categorisation, the  →  decay candidate photons are required to pass isolation criteria [82].For the merged category, only a track based isolation criterion is used by requiring the sum of the  T of the reconstructed tracks within a cone of Δ = 0.2 around the photon to be less than 5% of the photon  T .For the resolved category photons both track and calorimetric isolation are used.For the track isolation, if Δ(, ) > 0.22 the same track-based photon isolation is used as for the merged category; otherwise no track isolation is applied.For the calorimeter isolation, the sum of the topological clusters  T in a cone size Δ = 0.2 around the photon is required to be less than 6.5% of the photon  T .The isolation cuts efficiency exceeds 94% and 51% for isolated photons that pass the identification criteria discussed in Section 4 and the merged and resolved selection respectively.The total event selection acceptance times efficiency for signal events in the full mass range in the merged and resolved categories is in the range 4-11% and 0.4-2%, respectively for ( →  → ℓℓ) generated events.Overlayed on (c) for comparison is a signal distribution for   = 9 GeV, and the signal plus background fit (solid blue) with its background component (dashed blue).For (b) and (c) the branching ratio of the Higgs boson decay to  times the branching ratio  →  is assumed to be 50%.

Resolved category
The events in the resolved category are required to have the invariant mass of the  system,    , in the range 110-140 GeV, so that they are compatible with the 125 GeV Higgs boson decay.The    distributions for data and simulation are shown in Fig. 1(a).Here simulation is used to provide insight into the background composition but is not employed for the final background estimation.The main source of background is +jets production, amounting to about 90% of all background events.Photons in this background mostly arise from  0 and other light meson decays.The remaining 10% of the background events are due to photons from  +  production.
The final result of the search in the resolved category is extracted by fitting the invariant mass of the  pair,   .For the signal, the   distributions are modelled with a parametric function of which the parameters are interpolated to include signal masses for which simulation is not available.A double-sided crystal ball function [85] is found to be an excellent choice for the   distribution shape.Its parameters are interpolated to cover any  particle mass using linear functions.Examples of signal   distributions in the mass range used for this category are shown in Fig. 1(b), assuming a branching ratio of the Higgs boson decay to  times the branching ratio  →  of 50%.
The background distribution corresponds to a smooth and falling spectrum for   ≥ 2.0 GeV.This allows the background estimation to be performed with a data-driven method in which continuum shape is parametrised by an analytic function, which is derived in a control region obtained by inverting the    requirement.This method defines the region of interest for the resolved category to start from   = 2.0 GeV.
The parametric form of the background is:

𝑑
, with four parameters , ,  and .It has two components: a Fermi-Dirac probability function to describe the bulk of the distribution, and an exponential function to capture the tail end of the distribution.This functional form is chosen based on its ability to adequately model the background shape in several control regions across the full invariant mass range of the search.
The bias arising from the choice of the background model is evaluated by fitting the signal and background contributions to data that are known to have no signal [86].In this way, the biases inherent to this method, which is known as "spurious signal", are estimated.Templates are created using data or simulation in the control region, and are used as probability distribution functions to generate distributions with the expected number of events in the signal region.These distributions are fitted with a signal and background model and the "spurious signal" is extracted.The median values of the extracted signal for each   assumption is used across all templates to define an envelope that approximates the bias induced by the method.This bias is included as a systematic uncertainty in this category.Its value depends on the mass range, and the impact is limited as the size is well below the statistical uncertainties in the signal extraction fit for the entire mass range.
Post-fit   distributions with a binning of 200 MeV in the signal region and including all the systematic uncertainties on the background, as described earlier, and for signal (see Section 5.3) are shown in Fig. 1(c).
In this figure, a representative axion signal with a mass of   = 9 GeV and a branching ratio of the Higgs boson decay to  times the branching ratio  →  of 50% is shown.

Merged category
The events in the merged category are required to have the mass of the  system,    in the range 110-130 GeV.In addition, the photon must satisfy  ratio > 0.8.The quantity  ratio is the ratio of the energy difference between the maximum energy deposit and the energy deposit in the secondary maximum in the photon cluster to the sum of these energies [87].This variable is very efficient in discriminating photons against jets mis-identified as photons.The  ratio requirement reduces the background from +jets by 60%, leaving  +  as the dominant background, corresponding to about 75% of the expected background events in the signal region.The    pre-fit distribution after the  ratio requirement is shown in Fig. 2(a).
The background estimation in the merged category uses simulation with shape corrections derived from a control region that is composed of events that fail the    requirement.This control region is referred to as the "sideband" region.The shape corrections are applied to  +  and +jets events to correct for  differences between simulation and data for a number of chosen variables.The sideband region primarily consists of an enriched population of  +  events.Conversely, a sample enriched in +jets events is obtained by reversing the  ratio selection within the sideband.The distribution of  ratio in this control region before the fit is shown in Fig. 2(c).
The +jets and  +  modelling is tested and corrected using data from the sideband control region.In particular, corrections on the shape of the two-lepton system  T ,   T , the two-lepton-plus-photon system  T ,    T , the angular distance between the photon and the closest lepton, Δ ℓ  , and the photon  T are individually derived by comparing the shapes of these distributions in the sideband region in data and in simulation.The most relevant correction is for the shape of   T and is shown pre-fit after the  ratio requirement in Fig. 2(b).
After applying all the corrections, the signal is extracted through a simultaneous fit of two key observables: the angular distance between the two-lepton system and the  →  candidate photon, Δ   , in the signal region; and the  ratio distribution in the sideband region.This approach effectively constrains the normalisation of the +jets and  +  backgrounds using the available data in the sideband region.To ensure enhanced purity of the regions and minimise correlation between normalisation factors, the  ratio distribution is divided into two bins.The first bin represents  ratio < 0.8, while the second bin corresponds to  ratio > 0.8.In the signal region, as shown in Fig. 2(d) post-fit, the background exhibits a pronounced peak for large Δ   values, while the signal shape demonstrates a long tail for lower Δ   values.

Systematic uncertainties
Systematic effects from both experimental and theoretical sources affect the signal and background estimation used in deriving the results of the analysis.These may change the number of expected events and the shape of the fitted distributions.The following uncertainties are included in the signal distributions of both the merged and resolved categories and in the background estimation of the merged category.Dedicated systematic uncertainties coming from the background estimation in the resolved category and the +jets or  +  modelling in the merged category are described in Sections 5.1 and 5.2, respectively.
Background shape uncertainties of the MC distributions in the merged category are evaluated from the shape corrections applied to the samples.The full size of the shape correction is taken as an estimate of the associated systematic uncertainty.The impact of these uncertainties on the final result applied to the  +  process in the Δ   distribution, is found to be less than 0.5%.The most significant impact is from the shape of   T in the +jets process.However the impact on the final result is small as the +jets is a sub-dominant background.
Normalisation uncertainties related to modelling corrections of electrons and muons in the simulation with respect to data as well as the trigger efficiencies are considered.The largest impacts are from the electron identification and the muon reconstruction, which are found, however, to be less than 1.5% and 0.5% respectively across all simulated samples.Uncertainties on the data to simulation correction factors for photon identification and isolation efficiencies are evaluated as the difference between the correction factors applied or not, and they are found to be at maximum 0.3% and 6.5% respectively for the signal, and 0.3% and 2.0% for the background processes.The photon energy scale and resolution uncertainties are found to be negligible and thus not included in the analysis.
A pile-up modelling uncertainty is assigned to account for the difference between the predicted and measured inelastic   cross-sections [88].This uncertainty is evaluated to have a less than 1.5% effect on the event yield across all signal samples and less than 15% for the background processes.
Modelling uncertainty arises from various sources, including limited MC sample statistics, uncertainty in the renormalisation and factorisation scales, the choice of MC generator for both the signal and backgrounds, and the uncertainty in the parton distribution functions.These are evaluated to be less than 10% and 2% for the gluon-gluon fusion and the VBF processes, respectively.
The uncertainty in the integrated luminosity of the combined dataset from 2015 to 2018 is 1.7%, using the same methodology as in Ref. [34], and obtained using the LUCID-2 detector [89] for the primary luminosity measurements, complemented by measurements using the inner detector and calorimeters.
The background parametrisation in the resolved category is data-driven and extracted from the fit in the sideband region, thus only systematic uncertainties related to this fit and the signal modeling are considered.The systematic variation of the fit parameters is found to have a negligible impact on the fit result.The bias from the choice of the background model is evaluated with the "spurious signal" method, as discussed in Section 5.1.
All systematic uncertainties are implemented through nuisance parameters in the likelihood function used to fit the data and are profiled in the final fit, as described in more detail in Section 6.A summary of the impact on the limit for the merged category of each systematic uncertainty is reported in Table 1 for a signal-plus-background fit assuming a signal model with an axion mass of 1 GeV.These are computed by repeating the fit procedure after fixing the nuisance parameter under evaluation to its best-fit value, with all other parameters floating.The impact is then obtained by comparing the quadratic difference between the regular fit and the alternative fit with the parameter fixed.The contributions to the 68% confidence interval of the fitted signal strength from different sources of systematic uncertainties, for a signal-plus-background fit for the merged category assuming a signal model with an axion mass of 1 GeV.The evaluation is performed by fixing the parameter related to the tested source of systematic uncertainty to the best-fit value and redoing the fit with all other parameters floating.The impact is then obtained by comparing the quadratic difference between the default fit and the new model with the parameter fixed.The individual uncertainties can be correlated, so do not necessarily add up quadratically to the total uncertainty.The percentages show the size of the uncertainty relative to the parameter of interest.

Results
The statistical analysis of the data uses binned maximum-likelihood fits of the   and Δ   distributions in the resolved and merged categories, respectively.The search in the observed   distribution is performed in the 2 to 33 GeV mass range with a 200 MeV step, smaller than the  particle invariant mass resolution.The systematic uncertainties discussed in Section 5.3 are accounted for in the fit by means of nuisance parameters constrained by Gaussian penalty terms in the likelihood function.In the absense of signal, expected and observed 95% CL exclusion limits on the branching ratio of the Higgs boson decay to  times the branching ratio  →  are computed as a function of   using the CL  [90] modified frequentist approach.The exclusion limits are computed using the asymptotic approximation [91].
The 95% CL upper limits for both categories are shown in Fig. 3(a).The validity of the resolved category is limited to masses larger than 2 GeV and is complementary to the merged category.The search excludes a range of branching ratios of the branching ratio of the Higgs boson decay to  times the branching ratio  →  from 0.08% to 2%, depending on the assumed   .No significant excesses are observed.
The upper limits are also interpreted in the context of ALP models [3], in terms of the ALP mass and its effective coupling to photons, |  |/Λ.These are presented in Fig. 3

Conclusion
Data recorded by the ATLAS experiment at the LHC corresponding to an integrated luminosity of 139 fb −1 from proton-proton collisions at a centre-of-mass energy of 13 TeV, are used to search for a rare decay of the Higgs boson to a  boson and an axion particle , with a mass between 0.1 GeV and 33 GeV.The  boson is reconstructed using an electron or muon pair, while the particle  candidate is reconstructed from a pair of photons.The analysis accounts for the cases in which both photons are close enough to be reconstructed as a single photon in the detector, or topologies where they are reconstructed as two separate photons.No significant deviations with respect to the SM predictions are observed and upper limits are set on the branching ratio of the Higgs boson decay to  times the branching ratio  → , ranging from from 0.08% to 2% depending on the mass of the  particle.This is one of the first direct limits on the Higgs boson branching fraction into a  boson and an ALP decaying to a pair of photons, reaching sensitivity to parameter space that was not accessible by previous searches.

Figure 1 :
Figure 1: Distributions related to the resolved category: (a)    distribution along with the range, indicated by vertical dashed lines, in which the resolved category events are found; (b) signal   distributions for a selection of   values; (c)   post-fit distribution.Overlayed on (c) for comparison is a signal distribution for   = 9 GeV, and the signal plus background fit (solid blue) with its background component (dashed blue).For (b) and (c) the branching ratio of the Higgs boson decay to  times the branching ratio  →  is assumed to be 50%.

Figure 2 :
Figure 2: Distributions related to the merged category: (a)    distribution along with the range, indicated by vertical dashed lines, in which the merged category events after the  ratio cut are found; (b)   T distribution after categorisation, the  ratio cut and the corrections derived from the sideband control region; (c) pre-fit distribution of the  ratio variable in the sideband control region; (d) post-fit final discriminating variable Δ   in the signal region.Signal distributions for   values used in this category are overlayed for comparison, assuming a branching ratio of the Higgs boson decay to  times the branching ratio  →  of 100%.All the systematic uncertainties discussed in Section 5.3 are included in these plots.
(b) for different values of the Higgs coupling to , |   |/Λ), ranging from 0.05 TeV −1 to 0.4 TeV −1 .Depending on the |  |/Λ coupling values the ALP may have a sizeable lifetime resulting in displaced decays.The results of the analysis are fully accurate only for prompt axions decays, and only models with ALPs decaying up to the first ATLAS pixel layer (   = 33 mm) are considered.In the same figure, 90% CL exclusion regions on the ALP parameter space from direct experimental searches are also shown for comparison, with these contours being adapted from Refs.[2, 92].

Figure 3 :
Figure 3: (a) Expected and observed 95% CL upper limits on the branching ratio of the Higgs boson decay to  times the branching ratio  →  as a function of the  particle mass in the merged (  ≤ 2 GeV) and the resolved (  > 2 GeV) categories.(b) ATLAS observed 95% CL exclusion contours limits in terms of the ALP mass and its effective coupling to photons, |  |/Λ, for different values of the Higgs coupling to , |   |/Λ.The overlaid contour limits from other direct experimental searches are shown as well.The collider bounds (LHC, LEP, CDF) are displayed at 95% CL, while the remaining bounds (SN1987a, Cosmology and Beam Dump) are presented at 90% CL.The red band shows the preferred parameter space where the ( − 2) anomaly can be explained at 95% CL.These contours are adapted from Refs.[2, 92].