Direct constraint on the Higgs–charm coupling from a search for Higgs boson decays into charm quarks with the ATLAS detector

A search for the Higgs boson decaying into a pair of charm quarks is presented. The analysis uses proton–proton collisions to target the production of a Higgs boson in association with a leptonically decaying W or Z boson. The dataset delivered by the LHC at a centre-of-mass energy of and recorded by the ATLAS detector corresponds to an integrated luminosity of 139 fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{ fb}^{-1}$$\end{document}. Flavour-tagging algorithms are used to identify jets originating from the hadronisation of charm quarks. The analysis method is validated with the simultaneous measurement of WW, WZ and ZZ production, with observed (expected) significances of 2.6 (2.2) standard deviations above the background-only prediction for the (W/Z)Z(→cc¯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(W/Z)Z(\rightarrow c{\bar{c}})$$\end{document} process and 3.8 (4.6) standard deviations for the (W/Z)W(→cq)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(W/Z)W(\rightarrow cq)$$\end{document} process. The (W/Z)H(→cc¯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(W/Z)H(\rightarrow c {\bar{c}})$$\end{document} search yields an observed (expected) upper limit of 26 (31) times the predicted Standard Model cross-section times branching fraction for a Higgs boson with a mass of , corresponding to an observed (expected) constraint on the charm Yukawa coupling modifier |κc|<8.5(12.4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\kappa _c| < 8.5~(12.4)$$\end{document}, at the 95% confidence level. A combination with the ATLAS (W/Z)H,H→bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(W/Z)H, H\rightarrow b{\bar{b}}$$\end{document} analysis is performed, allowing the ratio κc/κb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa _c / \kappa _b$$\end{document} to be constrained to less than 4.5 at the 95% confidence level, smaller than the ratio of the b- and c-quark masses, and therefore determines the Higgs-charm coupling to be weaker than the Higgs-bottom coupling at the 95% confidence level.


Introduction
Since the discovery of a new particle, H , with a mass of approximately 125 GeV by the ATLAS [1] and CMS [2] collaborations at the LHC [3], studies of its properties have indicated that it is consistent with the Standard Model (SM) Higgs boson [4][5][6][7]. The interactions between the Higgs boson and the charged fermions of the third generation have been observed by both the ATLAS [8][9][10] and CMS [11][12][13] collaborations, and CMS has reported evidence for the decay e-mail: atlas.publications@cern.ch of the Higgs boson into a pair of muons [14], while ATLAS reported a 2σ excess over the background-only prediction [15]. Direct searches by the ATLAS and CMS collaborations for Higgs boson decays into a charm quark-antiquark pair, H → cc [16,17], decays into an electron-positron pair [18,19], exclusive decays into mesons [20][21][22][23][24][25], and reinterpretations of the Higgs p T spectrum [26,27], have not yet found experimental evidence for Higgs boson couplings to the first-generation fermions or second-generation quarks. Taken together, the results of these measurements and searches are consistent with the prediction that the coupling strength of the Higgs boson to each fermion scales proportionally to the fermion's mass.
In the SM the branching fraction of H → cc is 2.89% [28], approximately 20 times smaller than the branching fraction of the Higgs boson to a bottom quark-antiquark pair, H → bb. Physics beyond the Standard Model can significantly enhance or reduce the coupling of the Higgs boson to the charm quark, and in turn the H → cc branching fraction [29][30][31][32][33][34][35]. Direct searches for H → cc in protonproton ( pp) collisions have set upper limits on the crosssection times branching fraction for the production of a W or Z boson in association with a Higgs boson decaying into cc. The ATLAS Collaboration has performed a search in the Z (→ )H (→ cc) channel, where = e, μ, using 36.1 fb −1 of pp collision data recorded at √ s = 13 TeV [16], setting an observed (expected) upper limit at 110 (150) times the SM prediction, at 95% confidence level (CL). The CMS Collaboration has also performed a search using 35.9 fb −1 of pp collision data recorded at 13 TeV [17]; the search was conducted in three channels based on the number of charged leptons, namely the 0-, 1-and 2-lepton channels, targeting the Z H → ννcc, W H → νcc and Z H → cc signatures, respectively. These were combined to set an observed (expected) upper limit of 70 (37) times the SM prediction, at 95% CL. This paper presents a new search for V H(→ cc), where V = W or Z , using 139 fb −1 of pp collision data collected at a centre-of-mass energy of 13 TeV by the ATLAS detec-tor from 2015 to 2018. Events are selected in the 0-, 1-and 2-lepton channels and are categorised according to the transverse momentum, p T , of the vector boson and the number of jets.
Higgs boson candidates are constructed from the two jets with the highest p T . One of the main challenges in searching for H → cc is to recognise jets originating from the hadronisation of charm quarks. To identify these jets, a multivariate charm-jet tagging algorithm is used. Additionally, a bottom-jet identification algorithm is used to veto bottom jets, ensuring this analysis is orthogonal to the recent ATLAS V H(→ bb) measurement [36]. Events are selected if one or both of the two highestp T jets are c-tagged.
In order to search for the H → cc signal the distributions of the dijet invariant mass, m cc , in all event categories are used simultaneously in a binned maximum-likelihood fit, which allows the signal yield and the main background normalisations to be extracted. The analysis strategy is validated by the simultaneous measurement of the diboson processes in which one of the bosons decays to at least one charm quark, V W (→ cq) and V Z(→ cc), where q is a down-type quark. The result is interpreted in the kappa framework [37,38] in terms of κ c , the modifier of the coupling between the Higgs boson and the charm quark. The analysis is combined with the ATLAS V H, H → bb measurement [36] and the results are interpreted in the kappa framework in terms of both κ b and κ c , and in terms of the ratio κ c /κ b .

ATLAS detector
The ATLAS experiment [39] 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 (ID) surrounded by a thin superconducting solenoid providing a 2T axial magnetic field, electromagnetic and hadron 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. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|η| < 1.7). The endcap and forward regions are instrumented with LAr calorimeters for both 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the zaxis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of R ≡ ( η) 2 + ( φ) 2 . the EM and hadronic energy measurements up to |η| = 4.9.
The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer includes a system of precision chambers for tracking and fast detectors for triggering. A two-level trigger system is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz [40]. This is followed by a software-based trigger that reduces the accepted event rate to 1kHz on average depending on the data-taking conditions. An extensive software suite [41] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Dataset and simulated event samples
This analysis uses data recorded by the ATLAS detector during Run 2 of the LHC, which took place from 2015 to 2018 at a centre-of-mass energy of 13 TeV. Data were collected using a combination of missing transverse momentum triggers [42], in the 0-and 1-lepton channels, and single-lepton triggers [43,44], in the 1-and 2-lepton channels. Events are required to be of good quality and recorded while all relevant detector components were in operation [45]. The dataset corresponds to an integrated luminosity of 139.0 ± 2.4 fb −1 [46]. The Monte Carlo (MC) simulation samples used in this analysis are largely the same as those used in the ATLAS V H(→ bb) analysis [36], and are summarised in Table 1. Samples of simulated events were generated for V H production with a Higgs boson mass, m H , of 125 GeV, for both H → cc and H → bb decays, with branching fractions of 2.89% and 58.2%, respectively, and for the main background processes (tt, single-top, V + jets and diboson). The samples are used to optimise the analysis and perform the statistical analysis of the data.
The background from multi-jet events is negligible in the 0-and 2-lepton channels after applying the selection criteria detailed in Sect. 4. In the 1-lepton channel, it is estimated using a data-driven method. All samples of simulated events are initially normalised to the most accurate theoretical cross-section predictions currently available. Samples produced using alternative event generators are used to assess systematic uncertainties in the modelling of the signal and background processes, and are discussed in Sect. 5.
All samples of MC events were passed through the ATLAS detector response simulation [47] based on Geant4 [48] and were reconstructed with the same algorithms as used for data. The effect of multiple interactions in the same and neighbour- Table 1 Signal and background processes and their corresponding MC generators used in this analysis.The acronyms ME and PS stand for matrix element and parton shower, respectively. The crosssection order refers to the order of the cross-section calculation used for process normalisation in QCD, unless otherwise stated, with ((N)N)LO and ((N)N)LL standing for ((next-to-)next-to-)leading order and ((next-to-)next-to-)leading log, respectively Process ing bunch crossings (pile-up) was modelled by overlaying the simulated hard-scattering event with inelastic pp events generated with PYTHIA8.186 [49] using the NNPDF2.3lo set of parton distribution functions (PDF) [50] and the A3 set of tuned parameters (A3 tune) [51]. The simulated events were weighted to reproduce the distribution of the average number of interactions per bunch crossing seen in data. The EvtGen1.6.0 program [52]  . Electrons are required to have p T > 7 GeV and |η| < 2.47. They must satisfy the loose identification criterion, based on a likelihood discriminant combining observables related to the shower shape in the calorimeter and to the track matching the energy cluster, and are required to be isolated in both the ID and calorimeter using p T -dependent criteria. In the 1-lepton channel, more stringent requirements are placed on the identification and isolation of electrons. These electrons, called tight electrons, are required to satisfy the tight likelihood criterion and a stricter calorimeter-based isolation.
Muons are reconstructed within the acceptance of the muon spectrometer, |η| < 2.7 [89]. They are required to have p T > 7 GeV, to satisfy the loose identification criteria and to be isolated in the ID using p T -dependent criteria. As with electrons, in the 1-lepton channel more stringent requirements are placed on the identification and isolation of muons. These muons, called tight muons, must satisfy the medium identification criteria and a stricter track-based isolation, and have |η| < 2.5.
Hadronically decaying τ -leptons [90,91], identified with a medium quality criterion [91], are required to have p T > 20 GeV and |η| < 2.5, excluding the transition region of 1.37 < |η| < 1.52 between the barrel and endcap sections of the electromagnetic calorimeter. Reconstructed τ -leptons are not directly used in the event selection, but are used in the calculation of missing transverse momentum and to avoid double-counting reconstructed τ -leptons as other objects.
Jets are reconstructed from topological clusters of energy deposits in the calorimeter [92-94] by using the anti-k t algorithm [95,96] with a radius parameter of R = 0.4. The jets are classified as central or forward jets depending on their pseudorapidity. Jets are classified as forward jets if they have 2.5 < |η| < 4.5 and p T > 30 GeV. Jets are classified as central jets if they have p T > 20 GeV and |η| < 2.5. Additionally, central jets with p T < 120 GeV are required to be identified as originating from the primary vertex using a jetvertex tagging algorithm [97]. To improve the measurement of each jet's energy and direction, and consequently the measurement of m cc , if any muons are found within a cone of jetp T -dependent size around the jet axis, the four-momentum of the muon closest to the jet is added to the jet four-momentum, following the procedure described in Ref. [36]. An overlap removal procedure is applied to avoid double-counting between electrons, muons, hadronically decaying τ -leptons and jets.
Central jets are tagged as containing either b-or c-hadrons using two discriminants resulting from multivariate tagging algorithms, MV2 and DL1 [98]. Jets are b-tagged using the MV2 discriminant, configured to select b-jets with 70% efficiency in simulated tt events. A c-tagging configuration of the DL1 discriminant, DL1 c , was optimised for this analysis, and includes a veto on jets b-tagged by the MV2 algorithm. This configuration gives an average efficiency of 27% to tag c-jets in simulated tt events, and b-and light-jet misidentification rates of 8% and 1.6%, respectively. The efficiencies in simulation are calibrated to match those in data using control samples of tt and Z+jets events with a precision of 5%-10%, using methods identical to those applied to b-tagging algorithms [98-100]. Jets in simulated events are labelled using information from the MC generator's event 'truth' record, exclusively as b-, c-, or τ -jets, in this order, according to whether they contain a b-hadron, c-hadron, or τ -lepton with p T > 5 GeV within a cone of size R = 0.3 around their axis. Jets not labelled as b-, c-or τ -jets are labelled as light jets. Diboson, V + jets and top-quark backgrounds are classified according to the flavour labels of the jets that form the Higgs boson candidate in those selected events.
To maximise the statistical power of the available MC samples, the c-tagging requirement is not applied to the diboson, V + jets or top-quark samples. Instead, events are weighted by the probabilities for each jet to be c-tagged, based on its flavour label and as a function of the jet p T and |η|, to obtain predictions for events with either one or two ctagged jets. This is referred to as truth-flavour tagging since it uses information from the MC generator's event 'truth' record. In the V + jets samples, differences of up to 20% are observed between the two methods and are corrected for by weights assigned to each jet, dependent on the R to the closest other jet and on the flavour labels of the jet and the closest other jet. Finally, to correct for any residual non-closure in the truth-flavour tagging procedure, small normalisation corrections are applied to the diboson, V + jets and top-quark predictions such that the number of events for each process in each analysis region (defined in Sect. 4.2) matches that obtained when directly applying c-tagging. These normalisation corrections vary between 0.9 and 1.05.
The missing transverse momentum, E miss T is reconstructed as the negative of the vector sum of the transverse momenta of electrons, muons, hadronically decaying τ -leptons, jets, and a 'soft term' which is constructed from tracks associated with the primary vertex but not with any reconstructed object [101]. The magnitude of the E miss T is referred to as E miss T . The track-based missing transverse momentum, p miss T , is constructed using all ID tracks associated with the primary vertex and satisfying the quality criteria detailed in Ref. [102], with its magnitude denoted by p miss T .

Event selection and categorisation
Events are categorised into 0-, 1-and 2-lepton channels based on the number of loose electrons and muons they contain. Events with at least two central jets are selected, and they are further categorised as 2-or 3-jet events according to the total number of jets. Events with more than three jets are rejected in the 0-and 1-lepton channels to reduce the tt background.
In the 2-lepton channel, events with more than three jets are included in the 3-jet category.
Since the signal-to-background ratio increases for large transverse momentum of the vector boson, p V T , events with reconstructed p V T > 75 GeV are selected [103]. Two p V T regions are used: 75 GeV < p V T < 150 GeV (only in the 2-lepton channel) and p V T > 150 GeV. The p V T corresponds to E miss T in the 0-lepton channel, the magnitude of the vector sum of the E miss T and the lepton p T in the 1-lepton channel, and the magnitude of the vector sum of the two lepton transverse momenta in the 2-lepton channel.
The main discriminating variable in this analysis is the invariant mass, m cc , of the two central jets with the highest p T , hereafter referred to as signal jets. At least one signal jet must have p T > 45 GeV. Signal regions are composed of events in which one or both of these jets are c-tagged, with the two cases defining separate categories, referred to as 1-c-tag and 2-c-tag, respectively. Furthermore, any additional nonsignal jet must not be b-tagged. This requirement means that events in the signal regions can contain at most one b-tagged jet. Combined with an identical jet selection, this ensures that selected events are orthogonal to those selected in the ATLAS V H(→ bb) analysis [36]. Events selected in the control regions are not completely orthogonal with those selected in the V H(→ bb) analysis, and the impact of this is discussed in Sect. 7. In total, 16 signal regions are defined, arising from the combination of three lepton channels, two p V T categories (in the 2-lepton channel), two number-of-jets categories and two number-of-c-tagged-jets categories.
To reduce the background contamination in all channels, the R between the two signal jets is required to be R < 2.3 in events with 75 < p V T < 150 GeV, R < 1.6 in events with 150 < p V T < 250 GeV and R < 1.2 in events with p V T > 250 GeV. The R selection criteria, optimised for each p V T range, are approximately 80% efficient for signal events. For each signal region, a corresponding control region is defined as containing events failing the R selection, up to a maximum R of 2.5. These control regions, hereafter referred to as the high-R control regions, are designed to constrain the V + jets background normalisations and shapes.
In addition to the high-R control regions, further control regions are defined to constrain the modelling of the most important other background processes. Top control regions, enriched in tt and single-top events, are defined in all lepton channels. In the 0-and 1-lepton channels, events with three jets are selected, and in these events exactly one of the signal jets is c-tagged and the non-signal jet is b-tagged, resulting in one control region for each of these lepton channels. In the 2-lepton channel a pure sample of tt events is selected by requiring the two leptons to have different flavours (eμ) and opposite electric charges. One control region is defined for each number-of-jets category and p V T region combination, resulting in a total of four 2-lepton top control regions, each containing exactly one c-tagged jet.
Finally, in the 1-and 2-lepton channels, events in which neither of the two signal jets are c-tagged and no non-signal jets are b-tagged are used to constrain the normalisation of the V + light-jets backgrounds. One control region is defined for each number-of-jets category and p V T region combination, for a total of six 0-c-tag control regions.
A total of 44 analysis regions are defined: 16 signal regions, 16 high-R control regions, 6 top control regions and 6 0-c-tag control regions. The signal-region selection criteria specific to each channel are described below, and summarised in Table 2. 0-lepton channel Data were collected using E miss T triggers with thresholds ranging from 70 GeV in 2015 to 110 GeV in 2018 [42]. Events must not contain any loose electrons or muons, and are required to have E miss T > 150 GeV. At 150 GeV the E miss T triggers are approximately 75-90% efficient, depending on the year, reaching a full efficiency plateau at about 200 GeV. A requirement on the scalar sum of jet transverse momenta, H T , of 120 (150) GeV in 2-jet (3-jet) events is imposed to remove a small region of phase space where the trigger efficiency depends on the number of jets. To remove non-collision backgrounds, p miss T is required to exceed 30 GeV. Background multi-jet events with high E miss T typically arise from mismeasured jet energies in the calorimeter and can be rejected using angular separation requirements (detailed in Table 2) between the jets, E miss T and p miss T . 1-lepton channel Events must contain exactly one loose lepton, that is then required to also be tight. If the lepton is an electron (muon), it must have p T > 27 (25) GeV  In the muon sub-channel, data were collected with the same E miss T triggers as in the 0lepton channel. The online E miss T calculation does not include muons, so these triggers effectively select on p V T and perform more efficiently than single-muon triggers in the analysis phase space. In the electron subchannel, single-electron triggers were used to collect data with thresholds starting at 24 GeV for data collected in 2015 and 26 GeV for data collected between 2016 and 2018 [43]. In the electron sub-channel, there is a significant background from events with jets misidentified as electrons at low E miss T . These events are rejected by requiring E miss T > 30 GeV. The transverse mass of the reconstructed W boson, 2 m W T , is required to be less than 120 GeV. The number of multi-jet background events that survive the 1-lepton channel event selection is estimated in each analysis region by performing a template fit to the m W T distribution, which offers good discrimination between multi-jet and simulated backgrounds. The multi-jet m W T templates are extracted from control regions defined by inverting the tight isolation requirements on the leptons, after subtraction of the simulated backgrounds. The shapes of the multi-jet m cc distributions are obtained using the same procedure. More details of the template-fit method can be found in Ref. [9]. 2-lepton channel The 2-lepton channel must contain exactly two same-flavour loose leptons. At least one of the leptons must have p T > 27 GeV to be consistent with the online single-lepton trigger selection. In the dimuon case, they must have opposite electric charges. This requirement is not imposed in the dielectron subchannel due to a higher probability of charge misidentification. The invariant mass of the two leptons, m , is required to be consistent with the mass of the Z boson, 81 < m < 101 GeV.
Following the event selection, the Z + jets, W + jets and tt processes constitute the main backgrounds in the 0-lepton channel. In the 1-lepton channel, the main backgrounds arise from the W + jets and tt processes. For both the 0-and 1-lepton channels, the relative background composition depends substantially on the analysis region. In the 2-lepton channel the main background is Z + jets in all regions. The efficiency to select the V H(→ cc) signal, in which the V decays to leptons, is ≈ 1-2%, and the expected signal-to-background ratio in the mass range 100 GeV < m cc < 150 GeV is (1-7) × 10 −4 in the 1-ctag signal regions and (0.6-8) × 10 −3 in the 2-c-tag signal regions.

Systematic uncertainties
The sources of systematic uncertainties affecting this analysis can be broadly divided into two groups: those related to experimental effects and those due to the theoretical modelling of signal and background processes. The estimation of these uncertainties closely follows the procedures outlined in Ref. [36].
, where E miss T is used as an approximation for p ν T .

Experimental uncertainties
The leading experimental uncertainties in this analysis are due to imperfect calibration of the c-tagging efficiency, jet energy scales and jet energy resolutions. Correction factors for c-and b-tagging are determined from the difference between tagging efficiencies in data and simulation and are derived separately for c-jets, b-jets and light-flavour jets [98-100]. The uncertainties in the correction factors originate from multiple sources and are decomposed into independent components that are correlated between different analysis regions. Two additional uncertainties are included in MC samples where truth-flavour tagging, described in Sect. 4, is used. An uncertainty is included in the V + jets samples, equal to the R correction that is applied to improve agreement between the truth-tagged and direct-tagged simulation samples, and for each MC prediction an uncertainty is included in the overall normalisation correction between direct and truth-flavour tagging. The uncertainties in the calibration of jet energy scales and resolutions are estimated from multiple measurements [92]. Uncertainties in the jet energy scale and resolution are combined into independent components that are correlated between different analysis regions. An additional uncertainty in the calibration of b-and c-jets is also included.
Uncertainties in the reconstruction, identification, isolation and trigger efficiencies of electrons and muons, and the uncertainties in their energy scale and resolution, have been measured in data and found to be negligible compared to other uncertainties [88,89]. These uncertainties, along with the jet energy scale and resolution uncertainties, are propagated to the calculation of E miss T following the method described in Ref. [36]. The E miss T calculation has additional uncertainties associated with the p T scale, p T resolution and reconstruction efficiency of the tracks used to build the E miss T soft term, and with the modelling of the underlying event [101]. An uncertainty in the E miss T trigger efficiency is also included.
The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [46], obtained using the LUCID-2 detector [104] for the primary luminosity measurements. The average number of interactions per bunch-crossing in simulation is scaled by 1.03 to improve agreement with data, with an uncertainty corresponding to the full size of the correction (±3%).

Signal and background modelling uncertainties
Modelling uncertainties are evaluated using samples of simulated events. For each process, four categories of uncertainty are considered: cross-section and acceptance uncertainties, which account for the overall normalisation of backgrounds that are not allowed to float freely in the global like- Table 3 Summary of the background modelling systematic uncertainties considered

Top(b) normalisation Floating
Top(other) normalisation Floating

Multi-jet (1-lepton)
Normalisation 20%-100% The values given refer to the size of the uncertainty affecting the yield of each background. Where the size of an acceptance systematic uncertainty varies between analysis regions, a range is displayed. Uncertainties in the shapes of the m cc distributions are not shown below, but are taken into account for all backgrounds. CR and SR stand for control region and signal region lihood fit; flavour-fraction uncertainties, which account for uncertainties in the make-up of subcomponents of each background; relative acceptance uncertainties, which account for the relative normalisations of backgrounds in cases where the overall normalisation is considered correlated across two or more analysis regions; and shape uncertainties, accounting for uncertainties in the shape of the m cc distribution in each region. Unless stated otherwise, normalisation and shape effects of each systematic uncertainty are treated as uncorrelated with one another to ensure that modelling uncertainties are not artificially constrained. The decision to correlate or decorrelate normalisation and shape effects is made for each modelling uncertainty considered and is based on studies of different correlation models, and in cases where the impact of the choice is non-negligible, the model giving the most conservative uncertainty is used. Background modelling uncertainties considered in this analysis are summarised in Table 3.
VH Uncertainties in the V H(→ cc) signal are evaluated following the recommendations of the LHC Higgs Working Group [105,106], and include uncertainties in the cross-section of V H production and in the H → cc branching fraction. In addition, acceptance and shape uncertainties are evaluated by comparing the nominal V H(→ cc) samples with alternatives generated using Powheg+Herwig7 [107], and by independently varying the renormalisation (μ r ) and factorisation (μ f ) scales by factors of one-half and two in the nominal generator.
Comparisons are made separately for the three production processes; qq → Z H, qq → W H, and gg → Z H. Uncertainties in the total acceptance are found to be similar between lepton channels, 4-6% in the quark-initiated processes and 31-35% for gg → Z H. Similarly, uncertainties in the ratio of 3-jet to 2-jet events are found to be similar between channels for the quark-initiated processes, 6-12%, while uncertainties in the gg → Z H processes are larger, 19-56%. In the 2-lepton channel an uncertainty is included to account for differences in acceptance between the two p V T regions, and is 2% for qq → Z H and 5% for gg → Z H. Uncertainties in the shapes of the m cc distributions for each of the three production processes are evaluated in a similar way, and good agreement is found between lepton channels, allowing the use of one shape uncertainty per production process. Uncertainties in the normalisation of the W H(→ bb) and Z H(→ bb) background are taken from the recent ATLAS measurements [36]. Uncertainties in the number of jets and p V T acceptance ratios are set to be the same as those derived for the H → cc signal. Shape uncertainties in the H → bb background are evaluated in an equivalent way to the H → cc signal shape uncertainties.
Diboson Uncertainties in the diboson prediction are evaluated by comparing the nominal diboson MC samples, generated using Sherpa, with alternative samples generated using Powheg+Pythia8 and by independently varying μ r and μ f in the nominal samples by factors of one-half and two. Inclusive acceptance uncertainties are assigned to the W W , W Z and Z Z processes, and uncertainties in the ratios of 3-jet to 2-jet events and the ratios of events in different p V T regions are also included. Uncertainties in the m cc shape of the diboson signal, V W (→ cq) and V Z(→ cc), and background components are evaluated separately by comparing MC events from the two generator set-ups in each lepton channel. These comparisons are made inclusively in the number of jets and, in the 2-lepton channel, split into p V T regions. The shape uncertainties are found to be consistent between channels.
V+jets production The largest backgrounds in this analysis originate from V + jets, with W + jets and Z + jets being the largest background in the 1-lepton and 2-lepton channels, respectively, and a combination of the two making up the majority of the background in the 0lepton channel. In all channels, the V + jets background is divided into three subsamples based on the flavours of the two signal jets: V + hf, V + mf and V + lf, where hf, mf and lf stand for heavy-flavour, mixed-flavour and light-flavour, respectively. The V + hf background consists of the V + bb and V + cc contributions, the V + mf background consists of the V + bl, V + cl and V + bc contributions, where l refers to a light jet that doesn't contain a heavy-flavour hadron or τ -lepton, and the V + lf background consists of the remaining contribution where neither of the signal jets contain a heavy-flavour hadron. In the 0-lepton channel a non-negligible component of events contains a W → τ ν decay in which the τ -lepton decays hadronically and is selected as one of the two signal jets. These jets are considered as light jets for the purpose of assigning these events to the W + mf and W + lf background components. Uncertainties in the V + jets backgrounds due to normalisation, acceptance, flavour-fractions and shapes are evaluated using alternative Monte Carlo generator setups. These include taking the same Sherpa set-up used to generate the nominal V + jets samples but varying μ r and μ f by factors of one-half and two independently, and separate samples generated using Mad-Graph5_aMC@NLO [108] at leading order in QCD with up to four additional partons in the matrix element calculation, interfaced to PYTHIA8. The normalisations of all the W + jets and Z + jets components are free to float in the likelihood fit to data. The W + lf and Z + lf components float independently in the 2-jet and 3-jet signal regions. The Z + hf, Z + mf and Z + lf components float independently in the two p V T regions. The normalisations of the V + lf backgrounds are constrained by the 0-c-tag control regions, while the V + hf and V + mf components are constrained by the signal regions and high-R control regions. Uncertainties are included to account for acceptance effects between number-of-jet categories and lepton channels. Uncertainties in the relative contributions of the components of the V + hf, V + mf and V + lf backgrounds are found to be consistent between lepton channels, so one uncertainty per component is used, taken from the lepton channel which offers the most precise estimate. Shape uncertainties are derived in each analysis region and channel, separately for W + jets and Z + jets, with an uncertainty being included for each of the three comparisons performed, namely the comparisons between the nominal Sherpa sample and the μ r and μ f variations, and the comparison with MadGraph5_aMC@NLO+ Pythia8. The comparison between Sherpa and MadGraph5_aMC@NLO+Pythia8 is performed in a two-step process. First, a comparison is made in the high-R control region, with differences in both the shape and normalisation between the two models propagated to the signal regions in a correlated way. Second, the Sherpa model is weighted such that the two models agree in the high-R control region and any residual difference between the two models in the signal regions is included as a shape-only uncertainty.

Top-quark background
The top-quark background comprises tt and single-top-quark events. In the 2-lepton channel, the normalisation of this background in each signal region is determined using the corresponding top control region described in Sect. 4.2. Due to the small size of this background, no shape uncertainties are considered. In the 0-and 1-lepton channels, where this background is larger, tt and single-top W t-channel events are divided into two components based on the flavours of the two signal jets: top(b) events, in which at least one of the signal jets originates from a b-quark; and top(other) events. For the latter component, the two signal jets are mostly produced in the decay of a W boson and therefore their invariant mass peaks at the W -boson mass, while for the former component it does not. The normalisations of these two components are free to float separately in the global likelihood fit, with information from the top control regions contributing significantly. The background from t-and schannel single-top-quark production is small and is considered separately, with uncertainties in the t-channel and s-channel production cross-sections included.
Acceptance and shape uncertainties are derived separately for each component by comparing the nominal Powheg+Pythia8 tt and single-top MC samples with alternative samples, generated using Powheg+Herwig7 and MadGraph5_aMC@NLO+Pythia8 [108]. In addition, the impact of additional radiation is assessed using Powheg+Pythia8 samples with modified parameter values. Uncertainties are included to cover differences in the normalisation between lepton channels, between number-of-jets categories, and between the signal regions and top control regions. The dominant single-top contribution comes from W t production. The estimated uncertainty in the relative contributions of tt and W t events to the total top-quark background is included, as is an uncertainty due to the interference between tt and W t events, evaluated using an alternative W t MC sample in which the tt/W t interference is dealt with using the diagram subtraction scheme instead of the diagram removal scheme used in the nominal W t sample [109,110]. Shape differences in each of the various MC sample comparisons are considered as separate shape uncertainties. Shape uncertainties are derived for each component of the total top-quark background, top(b) and top(other), separately for tt and W t events.
Multi-jet Uncertainties in the multi-jet background are evaluated separately in the electron and muon subchannels. Normalisation and shape uncertainties are derived by changing the definition of the multi-jet control region and by modifying the normalisation of the tt and W+jets backgrounds in this control region by up to 25%. Additionally, the impact of using an alternative variable instead of m W T , namely the azimuthal angle between the charged lepton and E miss T , φ( , E miss T ), in the template fit is considered as an uncertainty.

Statistical analysis and results
A binned maximum-likelihood fit to the m cc distribution is performed across the 44 analysis regions to extract three parameters of interest (POI), μ. The parameters of interest, μ V H(cc) , μ V W (cq) and μ V Z(cc) , correspond to signal strengths that multiply the SM cross-sections times branching fractions of the V H(→ cc), V W (→ cq) and V Z(→ cc) processes, and are extracted by maximising the likelihood function with respect to both μ and nuisance parameters, which account for the systematic uncertainties discussed in Sect. 5. Uncertainties are constrained with Gaussian or log-normal distributions in the likelihood function with the exception of the normalisations of the V + jets and top-quark backgrounds, which are allowed to float freely in the fit. The uncertainties due to the limited number of events in the simu-lated samples used in the fit are included using the Beeston-Barlow technique [111] for the total MC prediction, excluding the V H(cc) signal. Systematic uncertainties that exhibit large fluctuations are smoothed and uncertainties with negligible impact on the final results are 'pruned' following the procedures outlined in Ref. [112].
The V H(→ bb) background, expected to be about eight (two) times larger than the SM H → cc signal in the 1c-tag (2-c-tag) signal regions, is included in the likelihood function with uncertainties; however, at the present level of signal sensitivity it does not significantly impact the search for V H(→ cc).
The m cc resolution is studied using simulation in the 2lepton channel and its value is 10-20 GeV depending on the signal region. The resolution is better in the 2-jet signal regions than those with more than two jets, and better in the 2-c-tag signal regions than the 1-c-tag signal regions. In the 2-lepton channel, the resolution in the p V T > 150 GeV signal regions is better than in the 75 GeV < p V T < 150 GeV signal regions. The following m cc ranges and uniform binnings are used in the various signal and control regions: Selected post-fit m cc distributions, where all normalisations and nuisance parameters are adjusted by the likelihood fit, are shown in Fig. 1 for the 0-, 1-, and 2-lepton channels. Post-fit distributions for the remaining analysis regions can be found in the Appendix. Table 4 shows the values of the free-floating background normalisation parameters obtained from the likelihood fit to data. The fitted signal strengths are: The correlation between μ V H(cc) and μ V W (cq) (μ V Z(cc) ) is 17% (16%), while μ V W (cq) and μ V Z(cc) are 17% anticorrelated. The probability of compatibility with the SM, defined as all three POIs being equal to unity, is 84%. The observed (expected) significances of the V W (→ cq) and V Z(→ cc) signals are 3.8 (4.6) and 2.6 (2.2) standard deviations, respectively. For the μ V H(cc) signal strength,   Events / 10 GeV   is observed (expected) at 95% CL using a modified frequentist CL s method [113], with the profile-likelihood ratio as the test statistic [114], using the RooFit/RooStats framework [115,116]. The limits for the three lepton channels are summarised in Fig. 2.

ATLAS
The effects of systematic uncertainties are summarised in Table 5. For each POI, the statistical uncertainty is obtained from a fit in which all nuisance parameters are fixed to their post-fit values. The total systematic uncertainty is found by subtracting the squared statistical uncertainty from the squared total uncertainty, i.e. σ syst. = (σ 2 tot. − σ 2 stat. ) 1/2 . Similarly, the impact of a subset of the systematic uncertainties is assessed by performing the fit with only their nuisance parameters fixed to their post-fit values. For each POI, the impact is then computed as the square root of the decrease in the squared uncertainty of that POI between the nominal fit and the fit with the nuisance parameters fixed. Despite the additional uncertainties it introduces, the use of truthflavour tagging improves the expected limit on μ V H(cc) by about 10% due to the improved statistical precision in the MC predictions.
The improvements in this analysis relative to the previous ATLAS search for Z H, H → cc [16] are quantified by performing a fit in the 2-lepton channel to the 2015-2016 data, corresponding to 36 fb −1 . Using the same signal regions as the previous analysis a 36% improvement in the expected limit is found, with most of the improvement due to better flavour-tagging performance. After also including the new 2lepton signal and control regions introduced in this analysis, a 43% improvement in the expected limit is found. Adding the full Run-2 dataset, along with the 0-and 1-lepton channels, the expected limit is improved by a factor of five in this analysis, relative to the previous ATLAS search.
The m cc distributions for events with either one or two c-tagged jets, summed over all channels and regions, after background subtraction, and using the fitted signal strengths, are shown in Fig. 3.
The best-fit value of the V H(→ cc) signal strength is interpreted within the kappa framework [37,38], by reparam- where B SM H →cc is the H → cc branching fraction predicted in the SM.
Constraints on κ c are set using the profile-likelihood ratio test statistic and are shown at 95% CL for each of the three channels and for the combined likelihood fit in Fig. 4. The combination allows an observed (expected) constraint of |κ c | < 8.5 (12.4) to be set at the 95% CL.

Combination with V H, H → bb
A combination of the analysis presented in this paper with the ATLAS V H, H → bb measurement [36] is performed by creating a likelihood function that is the product of the individual likelihood functions of the two analyses.     Fig. 6a and b, respectively. The likelihood function is symmetric relative to the sign of κ c but not to the sign of κ b due to the inclusion of κ b in the parameterisation of σ (gg → Z H) [117], resulting in two minima in the expected likelihood scan. For most values of κ b , a value of κ c is allowed at 95% CL that compensates for the effect of κ b via the width of the Higgs boson and vice versa. The observed best-fit value is (κ b , κ c ) = (−1.02, 0). The differ-ence in the value of the log-likelhood function between the best-fit value and (κ b , κ c ) = (+1.02, 0) is 0.02. These constraints complement those from measurements of the Higgs boson p T spectrum [27,118]. The ratio |κ c /κ b | is constrained to be less than 4.5 at 95% CL (5.1 expected). The observed value is smaller than the ratio of the b-and c-quark masses, m b /m c = 4.578 ± 0.008 [119], and therefore constrains the coupling of the Higgs boson to the charm quark to be weaker than the coupling of the Higgs boson to the bottom quark at 95% CL. The profile likelihood scan, parameterised in terms of κ c /κ b , with κ b as a free parameter, is shown in Fig. 7.

Conclusion
A direct search for the decay of a Higgs boson to a charm quark-antiquark pair has been performed using 139 fb −1 of pp collision data recorded at √ s = 13 TeV by the ATLAS detector at the LHC. The search uses three channels, Z H → ννcc, W H → νcc and Z H → cc. Signal events are identified using a multivariate charm tagging algorithm.
To enhance the signal sensitivity, events are categorised according to the p T of the reconstructed vector boson, the number of jets and the number of c-tagged jets. The m cc observable is used as the main discriminant in the likelihood fit to extract the signal. The analysis strategy is validated with the study of diboson production, which is found to be consistent with the SM prediction, with observed (expected) significances of 2.6 (2.2) standard deviations for the V Z(→ cc) process and 3.8 (4.6) standard deviations for the V W (→ cq) process.
The analysis yields an observed (expected) limit of 26 (31 +12 −8 ) times the predicted SM cross-section times branching fraction for a Higgs boson, with a mass of 125 GeV, decaying into a charm quark-antiquark pair, at the 95% confidence level, the most stringent limit to date. The expected limit is a factor of five more stringent than in the previous ATLAS search for Z H, H → cc, due to the larger dataset, improved c-tagging, and inclusion of the 0-and 1lepton channels and additional signal and control regions. The result is interpreted in the kappa framework, considering effects on the Higgs boson width and setting all other couplings to their SM values, which results in an observed (expected) constraint on the charm Yukawa coupling modifier strength |κ c | < 8.5 (12.4), at the 95% confidence level.
A combination with the ATLAS H → bb measurement is performed. Interpreted in the kappa framework the combination constrains the observed ratio |κ c /κ b | to be < 4.5 at the 95% confidence level. This is less than the ratio of the b-and c-quark masses, m b /m c , and thus constrains the coupling of the Higgs boson to the charm quark to be weaker than the coupling of the Higgs boson to the bottom quark at the 95% confidence level. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes. All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN. Data associated with journal publications are also made available: tables and data from plots (e.g. cross section values, likelihood profiles, selection efficiencies, cross section limits, ...) are stored in appropriate repositories such as HEPDATA (http:// hepdata.cedar.ac.uk/). ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models. For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.org/). This information is taken from the ATLAS Data Access Policy, which is a public document that can be downloaded from http:  Fig. 14. Figures 15, 16  and 17 show the post-fit background composition, including in control regions, for the 0-, 1-and 2-lepton channels, respectively.    Events / 10 GeV Events / 10 GeV   Events / 10 GeV Events / 10 GeV       Events / 20 GeV Events / 20 GeV