1 Introduction

Since the discovery of a new particle, H, with a mass of approximately 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 of the Higgs boson into a pair of muons [14], while ATLAS reported a \(2\sigma \) 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\rightarrow c{\bar{c}}\) [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_{\text {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\rightarrow c{\bar{c}}\) is 2.89% [28], approximately 20 times smaller than the branching fraction of the Higgs boson to a bottom quark–antiquark pair, \(H \rightarrow b{\bar{b}}\). 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 \rightarrow c{\bar{c}}\) branching fraction [29,30,31,32,33,34,35]. Direct searches for \(H \rightarrow c{\bar{c}}\) in proton–proton (pp) collisions have set upper limits on the cross-section times branching fraction for the production of a W or Z boson in association with a Higgs boson decaying into \(c{\bar{c}}\). The ATLAS Collaboration has performed a search in the \(Z(\rightarrow \ell \ell )H (\rightarrow c{\bar{c}})\) channel, where \(\ell =e,\mu \), using \(36.1~\text{ fb}^{-1}\) of pp collision data recorded at [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~\text{ fb}^{-1}\) of pp collision data recorded at [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 \(ZH\rightarrow \nu \nu c{\bar{c}}\), \(WH\rightarrow \ell \nu c{\bar{c}}\) and \(ZH \rightarrow \ell \ell c{\bar{c}}\) 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 \(VH (\rightarrow c{\bar{c}})\), where \(V=W\) or Z, using \(139~\text{ fb}^{-1}\) of pp collision data collected at a centre-of-mass energy of by the ATLAS detector from 2015 to 2018. Events are selected in the 0-, 1- and 2-lepton channels and are categorised according to the transverse momentum, \(p_{\text {T}}\), of the vector boson and the number of jets.

Higgs boson candidates are constructed from the two jets with the highest \(p_{\text {T}}\). One of the main challenges in searching for \(H \rightarrow c{\bar{c}}\) 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 \(VH(\rightarrow b{\bar{b}})\) measurement [36]. Events are selected if one or both of the two highest-\(p_{\text {T}}\) jets are c-tagged.

In order to search for the \(H \rightarrow c{\bar{c}}\) 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, \(VW(\rightarrow cq)\) and \(VZ(\rightarrow c{\bar{c}})\), where q is a down-type quark. The result is interpreted in the kappa framework [37, 38] in terms of \(\kappa _c\), the modifier of the coupling between the Higgs boson and the charm quark. The analysis is combined with the ATLAS \(VH, H \rightarrow b{\bar{b}}\) measurement [36] and the results are interpreted in the kappa framework in terms of both \(\kappa _b\) and \(\kappa _c\), and in terms of the ratio \(\kappa _c / \kappa _b\).

2 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\pi \) coverage in solid angle.Footnote 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 \(|\eta | < 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 (\(|\eta | < 1.7\)). The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to \(|\eta | = 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.

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 cross-section 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

3 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\pm 2.4~\text{ fb}^{-1}\) [46].

The Monte Carlo (MC) simulation samples used in this analysis are largely the same as those used in the ATLAS \(VH(\rightarrow b{\bar{b}})\) analysis [36], and are summarised in Table 1. Samples of simulated events were generated for VH production with a Higgs boson mass, \(m_H\), of , for both \(H \rightarrow c{\bar{c}}\) and \(H \rightarrow b{\bar{b}}\) decays, with branching fractions of 2.89% and 58.2%, respectively, and for the main background processes (\(t{\bar{t}} \), 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 neighbouring 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] was used to describe the decays of bottom and charm hadrons in all samples, except those generated with Sherpa [53,54,55].

4 Object and event selection

4.1 Object selection

Interaction vertices are reconstructed from tracks in the ID. The vertex with the highest sum of squared transverse momenta of associated tracks is used as the primary vertex [87].

Electrons are reconstructed by matching ID tracks with energy clusters in the EM calorimeter [88]. Electrons are required to have and \(|\eta | < 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_{\text {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, \(|\eta | < 2.7\) [89]. They are required to have , to satisfy the loose identification criteria and to be isolated in the ID using \(p_{\text {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 \(|\eta | < 2.5\).

Hadronically decaying \(\tau \)-leptons [90, 91], identified with a medium quality criterion [91], are required to have and \(|\eta | < 2.5\), excluding the transition region of \(1.37< |\eta | < 1.52\) between the barrel and endcap sections of the electromagnetic calorimeter. Reconstructed \(\tau \)-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 \(\tau \)-leptons as other objects.

Jets are reconstructed from topological clusters of energy deposits in the calorimeter [92,93,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< |\eta | < 4.5\) and . Jets are classified as central jets if they have and \(|\eta | < 2.5\). Additionally, central jets with are required to be identified as originating from the primary vertex using a jet-vertex 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 jet-\(p_{\text {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 \(\tau \)-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 \(t\bar{t}\) events. A c-tagging configuration of the DL1 discriminant, \(\mathrm {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 \(t{\bar{t}}\) 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 \(t{\bar{t}} \) and \(\text {Z}\,+\,\text {jets}\) events with a precision of 5%–10%, using methods identical to those applied to b-tagging algorithms [98,99,100]. Jets in simulated events are labelled using information from the MC generator’s event ‘truth’ record, exclusively as b-, c-, or \(\tau \)-jets, in this order, according to whether they contain a b-hadron, c-hadron, or \(\tau \)-lepton with within a cone of size \(\Delta R = 0.3\) around their axis. Jets not labelled as b-, c- or \(\tau \)-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_{\text {T}}\) and \(|\eta |\), to obtain predictions for events with either one or two c-tagged 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 \(\Delta 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, \({\varvec{E_{\mathrm {T}}^{\mathrm {miss}}}}\) is reconstructed as the negative of the vector sum of the transverse momenta of electrons, muons, hadronically decaying \(\tau \)-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 \({\varvec{E_{\mathrm {T}}^{\mathrm {miss}}}}\) is referred to as \(E_{\mathrm {T}}^{\mathrm {miss}}\). The track-based missing transverse momentum, \({\varvec{p^{\mathrm {miss}}_{\mathrm {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^{\mathrm {miss}}_{\mathrm {T}}\).

4.2 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 \(t{\bar{t}}\) 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_{\text {T}} ^V\), events with reconstructed are selected [103]. Two \(p_{\text {T}} ^V\) regions are used: (only in the 2-lepton channel) and . The \(p_{\text {T}} ^V\) corresponds to \(E_{\text {T}}^{\text {miss}}\) in the 0-lepton channel, the magnitude of the vector sum of the \({\varvec{E^{\mathrm {miss}}_{\mathrm {T}}}}\) and the lepton \({\varvec{p_{\mathrm {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_{\text {T}}\), hereafter referred to as signal jets. At least one signal jet must have . 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 non-signal 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 \(VH(\rightarrow b{\bar{b}})\) analysis [36]. Events selected in the control regions are not completely orthogonal with those selected in the \(VH(\rightarrow b{\bar{b}})\) 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_{\text {T}} ^V\) 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 \(\Delta R\) between the two signal jets is required to be \(\Delta R < 2.3\) in events with , \(\Delta R < 1.6\) in events with and \(\Delta R < 1.2\) in events with . The \(\Delta R\) selection criteria, optimised for each \(p_{\text {T}} ^V\) range, are approximately \(80\%\) efficient for signal events. For each signal region, a corresponding control region is defined as containing events failing the \(\Delta R\) selection, up to a maximum \(\Delta R\) of 2.5. These control regions, hereafter referred to as the high-\(\Delta R\) control regions, are designed to constrain the \(V+\) jets background normalisations and shapes.

In addition to the high-\(\Delta R\) control regions, further control regions are defined to constrain the modelling of the most important other background processes. Top control regions, enriched in \(t{\bar{t}}\) 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 \(t{\bar{t}}\) events is selected by requiring the two leptons to have different flavours (\(e\mu \)) and opposite electric charges. One control region is defined for each number-of-jets category and \(p_{\text {T}} ^V\) 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_{\text {T}} ^V\) region combination, for a total of six 0-c-tag control regions.

Table 2 Summary of the signal region event selection in the 0-, 1- and 2-lepton channels. Jet1 and jet2 refer to the two signal jets and H refers to the jet1–jet2 system

A total of 44 analysis regions are defined: 16 signal regions, 16 high-\(\Delta 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_{\text {T}}^{\text {miss}}\) triggers with thresholds ranging from in 2015 to in 2018 [42]. Events must not contain any loose electrons or muons, and are required to have . At the \(E_{\text {T}}^{\text {miss}}\) triggers are approximately 75–90% efficient, depending on the year, reaching a full efficiency plateau at about . A requirement on the scalar sum of jet transverse momenta, \(H_{\mathrm {T}}\), of 120 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^{\mathrm {miss}}_{\mathrm {T}}\) is required to exceed . Background multi-jet events with high \(E_{\text {T}}^{\text {miss}}\) typically arise from mismeasured jet energies in the calorimeter and can be rejected using angular separation requirements (detailed in Table 2) between the jets, \({\varvec{E^{\mathrm {miss}}_{\mathrm {T}}}}\) and \({\varvec{p^{\mathrm {miss}}_{\mathrm {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 and . In the muon sub-channel, data were collected with the same \(E_{\text {T}}^{\text {miss}}\) triggers as in the 0-lepton channel. The online \(E_{\text {T}}^{\text {miss}}\) calculation does not include muons, so these triggers effectively select on and perform more efficiently than single-muon triggers in the analysis phase space. In the electron sub-channel, single-electron triggers were used to collect data with thresholds starting at for data collected in 2015 and 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_{\text {T}}^{\text {miss}}\). These events are rejected by requiring . The transverse mass of the reconstructed W boson,Footnote 2\(m_{\mathrm {T}}^W\), is required to be less than . 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_{\mathrm {T}}^W\) distribution, which offers good discrimination between multi-jet and simulated backgrounds. The multi-jet \(m_{\mathrm {T}}^W\) 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 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 sub-channel due to a higher probability of charge misidentification. The invariant mass of the two leptons, \(m_{\ell \ell }\), is required to be consistent with the mass of the Z boson, .

Following the event selection, the \(\text {Z}\,+\,\text {jets}\), \(\text {W}\,+\,\text {jets}\) and \(t{\bar{t}}\) processes constitute the main backgrounds in the 0-lepton channel. In the 1-lepton channel, the main backgrounds arise from the \(\text {W}\,+\,\text {jets}\) and \(t{\bar{t}}\) 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 \(\text {Z}\,+\,\text {jets}\) in all regions. The efficiency to select the \(VH(\rightarrow c{\bar{c}})\) signal, in which the V decays to leptons, is \(\approx 1\text{-- }2\)%, and the expected signal-to-background ratio in the mass range is \((1\text{-- }7) \times 10^{-4}\) in the 1-c-tag signal regions and \((0.6\text{-- }8) \times 10^{-3}\) in the 2-c-tag signal regions.

5 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].

5.1 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,99,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 \(\Delta 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_{\text {T}}^{\text {miss}}\) following the method described in Ref. [36]. The \(E_{\text {T}}^{\text {miss}}\) calculation has additional uncertainties associated with the \(p_{\text {T}}\) scale, \(p_{\text {T}}\) resolution and reconstruction efficiency of the tracks used to build the \(E_{\text {T}}^{\text {miss}}\) soft term, and with the modelling of the underlying event [101]. An uncertainty in the \(E_{\text {T}}^{\text {miss}}\) 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 (\(\pm 3\%\)).

Table 3 Summary of the background modelling systematic uncertainties considered

5.2 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 likelihood 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 \(VH(\rightarrow c{\bar{c}})\) signal are evaluated following the recommendations of the LHC Higgs Working Group [105, 106], and include uncertainties in the cross-section of VH production and in the \(H\rightarrow c{\bar{c}}\) branching fraction. In addition, acceptance and shape uncertainties are evaluated by comparing the nominal \(VH(\rightarrow c{\bar{c}})\) samples with alternatives generated using Powheg+Herwig7 [107], and by independently varying the renormalisation (\(\mu _\mathrm {r}\)) and factorisation (\(\mu _\mathrm {f}\)) scales by factors of one-half and two in the nominal generator. Comparisons are made separately for the three production processes; \(q{\bar{q}} \rightarrow ZH\), \(q{\bar{q}} \rightarrow WH\), and \(gg \rightarrow ZH\). 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 \rightarrow ZH\). 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 \rightarrow ZH\) processes are larger, 19–56%. In the 2-lepton channel an uncertainty is included to account for differences in acceptance between the two \(p_{\text {T}} ^V\) regions, and is 2% for \(q{\bar{q}} \rightarrow ZH\) and 5% for \(gg \rightarrow ZH\). 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 \(WH(\rightarrow b{\bar{b}})\) and \(ZH(\rightarrow b{\bar{b}})\) background are taken from the recent ATLAS measurements [36]. Uncertainties in the number of jets and \(p_{\text {T}} ^V\) acceptance ratios are set to be the same as those derived for the \(H\rightarrow c{\bar{c}}\) signal. Shape uncertainties in the \(H \rightarrow b{\bar{b}}\) background are evaluated in an equivalent way to the \(H\rightarrow ~c{\bar{c}}\) 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 \(\mu _\mathrm {r}\) and \(\mu _\mathrm {f}\) in the nominal samples by factors of one-half and two. Inclusive acceptance uncertainties are assigned to the WW, WZ and ZZ processes, and uncertainties in the ratios of 3-jet to 2-jet events and the ratios of events in different \(p_{\text {T}} ^V\) regions are also included. Uncertainties in the \(m_{cc}\) shape of the diboson signal, \(VW(\rightarrow cq)\) and \(VZ(\rightarrow c{\bar{c}})\), 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_{\text {T}} ^V\) 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 \(\text {W}\,+\,\text {jets}\) and \(\text {Z}\,+\,\text {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 0-lepton 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 \(\tau \)-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 \rightarrow \tau \nu \) decay in which the \(\tau \)-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 set-ups. These include taking the same Sherpa set-up used to generate the nominal \(V+\) jets samples but varying \(\mu _\mathrm {r}\) and \(\mu _\mathrm {f}\) by factors of one-half and two independently, and separate samples generated using MadGraph5_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 \(\text {W}\,+\,\text {jets}\) and \(\text {Z}\,+\,\text {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_{\text {T}} ^V\) 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-\(\Delta 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 \(\text {W}\,+\,\text {jets}\) and \(\text {Z}\,+\,\text {jets}\), with an uncertainty being included for each of the three comparisons performed, namely the comparisons between the nominal Sherpa sample and the \(\mu _\mathrm {r}\) and \(\mu _\mathrm {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-\(\Delta 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-\(\Delta 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 \(t{\bar{t}}\) 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, \(t{\bar{t}}\) and single-top Wt-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 s-channel 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 \(t{\bar{t}}\) 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 Wt production. The estimated uncertainty in the relative contributions of \(t{\bar{t}}\) and Wt events to the total top-quark background is included, as is an uncertainty due to the interference between \(t{\bar{t}}\) and Wt events, evaluated using an alternative Wt MC sample in which the \(t{\bar{t}}/Wt\) interference is dealt with using the diagram subtraction scheme instead of the diagram removal scheme used in the nominal Wt 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 \(t{\bar{t}}\) and Wt events.

  • Multi-jet Uncertainties in the multi-jet background are evaluated separately in the electron and muon sub-channels. Normalisation and shape uncertainties are derived by changing the definition of the multi-jet control region and by modifying the normalisation of the \(t{\bar{t}}\) and \(\text {W}\,+\,\text {jets}\) backgrounds in this control region by up to 25%. Additionally, the impact of using an alternative variable instead of \(m_{\mathrm {T}}^W\), namely the azimuthal angle between the charged lepton and \({\varvec{E_{\mathrm {T}}^{\mathrm {miss}}}}\), \(\Delta \phi (\ell , {\varvec{E_{\mathrm {T}}^{\mathrm {miss}}}})\), in the template fit is considered as an uncertainty.

Fig. 1
figure 1

Post-fit distributions for six selected signal regions out of 44 analysis regions, with two jets and for the 0-lepton (left), 1-lepton (centre), and 2-lepton (right) channels, with one c-tag (top) and two c-tags (bottom). The total signal-plus-background prediction is shown by the solid black line and includes the \(H\rightarrow c{\bar{c}}\) signal scaled to the best-fit value of \(\mu _{VH(c{\bar{c}})} = -9\). The \(H\rightarrow c{\bar{c}}\) signal is also shown as an unfilled histogram scaled to 300 times the SM prediction. The post-fit uncertainty is shown as the hatched background. The ratio of the data to the sum of the post-fit signal plus background is shown in the lower panel

6 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), \({\varvec{\mu }}\). The parameters of interest, \(\mu _{VH(c{\bar{c}})}\), \(\mu _{VW(cq)}\) and \(\mu _{VZ(c{\bar{c}})}\), correspond to signal strengths that multiply the SM cross-sections times branching fractions of the \(VH(\rightarrow c{\bar{c}})\), \(VW(\rightarrow cq)\) and \(VZ(\rightarrow c{\bar{c}})\) processes, and are extracted by maximising the likelihood function with respect to both \({\varvec{\mu }}\) 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 simulated samples used in the fit are included using the Beeston–Barlow technique [111] for the total MC prediction, excluding the \(VH(c{\bar{c}})\) 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 \(VH(\rightarrow b{\bar{b}})\) background, expected to be about eight (two) times larger than the SM \(H \rightarrow c{\bar{c}}\) signal in the 1-c-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 \(VH(\rightarrow c{\bar{c}})\).

The \(m_{cc}\) resolution is studied using simulation in the 2-lepton channel and its value is 10–20  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 signal regions is better than in the signal regions. The following \(m_{cc}\) ranges and uniform binnings are used in the various signal and control regions:

  • 16 bins from 50 to in the signal regions and 0- and 1-lepton top control regions, with the exception of the 2-lepton, 2-c-tag, signal regions.

  • 9 (10) bins from 50 to in the 2-lepton, 2-c-tag, , 2-jet (3-jet) signal region.

  • 13 bins from 80 to in each of the high-\(\Delta R\) control regions, with the exception of the 2-lepton, 2-c-tag, high-\(\Delta R\) regions, where 9 bins from 80 to are used.

  • A single bin from 50 to in each of the 0-c-tag control regions and 2-lepton top 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:

Table 4 Values of the free-floating background normalisation parameters obtained from the likelihood fit to data. The uncertainties represent the combined statistical and systematic uncertainties. Unless otherwise stated, normalisation parameters are correlated across all \(p_{\text {T}} ^V\) and number-of-jets analysis regions
$$\begin{aligned} \begin{array}{lllll} \mu _{VH(c{\bar{c}})}= &{} -9 &{} \pm 10 \mathrm{(stat.)} &{} \pm 12 \mathrm{(syst.)}\\ \mu _{VW(cq)}= &{} 0.83 &{}\pm 0.11 \mathrm{(stat.)}&{} \pm 0.21 \mathrm{(syst.)} \\ \mu _{VZ(c{\bar{c}})}= &{} 1.16 &{}\pm 0.32 \mathrm{(stat.)}&{} \pm 0.36 \mathrm{(syst.)}. \\ \end{array} \end{aligned}$$

The correlation between \(\mu _{VH(c{\bar{c}})}\) and \(\mu _{VW(cq)}~(\mu _{VZ(c{\bar{c}})})\) is 17% (16%), while \(\mu _{VW(cq)}\) and \(\mu _{VZ(c{\bar{c}})}\) are 17% anti-correlated. The probability of compatibility with the SM, defined as all three POIs being equal to unity, is 84%. The observed (expected) significances of the \(VW(\rightarrow cq)\) and \(VZ(\rightarrow c{\bar{c}})\) signals are 3.8 (4.6) and 2.6 (2.2) standard deviations, respectively. For the \(\mu _{VH(c{\bar{c}})}\) signal strength, an upper limit of \(26~(31^{+12}_{-8})\) is observed (expected) at 95% CL using a modified frequentist \(\mathrm {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.

Fig. 2
figure 2

The observed and expected 95% CL upper limits on the cross-section times branching fraction normalized to its SM prediction in each lepton channel and for the combined fit. The single-channel limits are obtained using a five-POI fit, in which each channel has a separate \(VH(\rightarrow c{\bar{c}})\) parameter of interest

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. \(\sigma _{\mathrm {syst.}} = (\sigma _{\mathrm {tot.}}^2 - \sigma _{\mathrm {stat.}}^2)^{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 truth-flavour tagging improves the expected limit on \(\mu _{VH(c{\bar{c}})}\) 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 \(ZH, H \rightarrow c{\bar{c}}\) [16] are quantified by performing a fit in the 2-lepton channel to the 2015–2016 data, corresponding to \(36~\text{ 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 2-lepton 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 \(VH(\rightarrow c{\bar{c}})\) signal strength is interpreted within the kappa framework [37, 38], by reparameterising \(\mu _{VH(c{\bar{c}})}\) in the likelihood function in terms of the Higgs–charm coupling modifier, \(\kappa _c\), assuming that the coupling modifier only affects the Higgs boson decays. Including effects in both the partial and full width, considering only SM decays and setting all other couplings to their SM predictions, \(\mu _{VH(c{\bar{c}})}\) is parameterised as a function of \(\kappa _c\)

$$\begin{aligned} \mu _{VH(c{\bar{c}})}(\kappa _c) = \frac{\kappa _c^2}{1 + B^{\mathrm {SM}}_{H \rightarrow c{\bar{c}}}(\kappa _c^2 - 1)}\,, \end{aligned}$$

where \(B^{\mathrm {SM}}_{H \rightarrow c{\bar{c}}}\) is the \(H \rightarrow c{\bar{c}}\) branching fraction predicted in the SM.

Table 5 Breakdown of contributions to the uncertainty in the fitted values of \(\mu _{VH(c{\bar{c}})}\), \(\mu _{VW(cq)}\) and \(\mu _{VZ(c{\bar{c}})}\). The sum in quadrature of uncertainties from different sources may differ from the total due to correlations. In cases where the upward and downward systematic variations have different values, the mean of the absolute values is shown
Fig. 3
figure 3

The post-fit \(m_{cc}\) distribution summed over all signal regions after subtracting backgrounds, leaving only the \(VH(\rightarrow c{\bar{c}})\), \(VW(\rightarrow cq)\) and \(VZ(\rightarrow c{\bar{c}})\) processes, for events with a one c-tag and b two c-tags. The red filled histogram corresponds to the \(VH, H\rightarrow c{\bar{c}}\) signal for the fitted value of \(\mu _{VH(c{\bar{c}})} = -9\), while the open red histogram corresponds to the signal expected at the 95% CL upper limit on \(\mu _{VH(c{\bar{c}})}~(\mu _{VH(c{\bar{c}})}=26)\). The hatched band shows the uncertainty of the fitted background

Constraints on \(\kappa _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 \(|\kappa _c| < 8.5~(12.4)\) to be set at the 95% CL.

7 Combination with \(VH, H\rightarrow b{\bar{b}}\)

A combination of the analysis presented in this paper with the ATLAS \(VH, H\rightarrow b{\bar{b}}\) measurement [36] is performed by creating a likelihood function that is the product of the individual likelihood functions of the two analyses. Two parameters of interest are used, \(\mu _{VH(c{\bar{c}})}\) and \(\mu _{VH(b{\bar{b}})}\) for the \(VH, H \rightarrow c{\bar{c}}\) and \(VH, H \rightarrow b{\bar{b}}\) signal strengths, respectively, and are included in both of the input likelihood functions. The importance of including both signal strengths in a combined likelihood was pointed out by Refs. [29, 30]. Experimental systematic uncertainties that are common to both analyses, detailed in Sect. 5.1, are considered correlated between the two analyses, with the exception of the flavour-tagging systematic uncertainties due to the different calibration procedures used for b- and c-tagging. Background normalisations and modelling uncertainties are considered uncorrelated between the two analyses.

$$\begin{aligned} \begin{array}{lll} \mu _{VH(c{\bar{c}})}= &{} -9 &{}\pm 10 \mathrm{(stat.)} \pm ~11 \mathrm{(syst.)} \\ \mu _{VH(b{\bar{b}})}= &{} 1.06 &{}\pm 0.12 \mathrm{(stat.)} ^{+0.15}_{-0.13} \mathrm{(syst.)} \\ \end{array} \end{aligned}$$

with a correlation coefficient of \(-12\%\). The fitted signal strengths are consistent with those found in the individual analyses. The expected and observed best-fit values and their 68% and 95% CL contours are shown in Fig. 5.

Although the signal regions of the two analyses are orthogonal due to the b-tagging veto used in the c-tagging definition, a small overlap of events occurs in the control regions used in the two analyses. To test the impact of this, events that appear in both analyses are removed from the \(VH, H\rightarrow c{\bar{c}}\) control regions. The results are unchanged. Treating the normalisations of the backgrounds as correlated between the two analyses is also tested and does not affect the expected sensitivity.

Fig. 4
figure 4

Expected and observed values of the negative profile log-likelihood ratio as a function of \(\kappa _c\). The single-channel likelihoods are obtained using a five-POI fit, in which each channel has a separate \(VH(\rightarrow c{\bar{c}})\) parameter of interest. The order of the lines in the plot matches the order of the lines in the legend

Fig. 5
figure 5

The observed and expected best fit values of \(\mu _{VH(c{\bar{c}})}\) and \(\mu _{VH(b{\bar{b}})}\) with their 68% and 95% CL contours

The best-fit values of \(\mu _{VH(b{\bar{b}})}\) and \(\mu _{VH(c{\bar{c}})}\) are interpreted in the kappa framework by parameterising the likelihood function in terms of both \(\kappa _b\) and \(\kappa _c\), while setting all other couplings to their SM predictions and considering only SM Higgs boson decays. Constraints on \(\kappa _b\) and \(\kappa _c\) are set using the profile-likelihood ratio test statistic. The expected and observed constraints are shown in Fig. 6a and b, respectively. The likelihood function is symmetric relative to the sign of \(\kappa _c\) but not to the sign of \(\kappa _b\) due to the inclusion of \(\kappa _b\) in the parameterisation of \(\sigma (gg \rightarrow ZH)\) [117], resulting in two minima in the expected likelihood scan. For most values of \(\kappa _b\), a value of \(\kappa _c\) is allowed at 95% CL that compensates for the effect of \(\kappa _b\) via the width of the Higgs boson and vice versa. The observed best-fit value is \((\kappa _b, \kappa _c) = (-1.02, 0)\). The difference in the value of the log-likelhood function between the best-fit value and \((\kappa _b, \kappa _c) = (+1.02, 0)\) is 0.02. These constraints complement those from measurements of the Higgs boson \(p_{\text {T}} \) spectrum [27, 118]. The ratio \(|\kappa _c/\kappa _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 \pm 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 \(\kappa _c/\kappa _b\), with \(\kappa _b\) as a free parameter, is shown in Fig. 7.

Fig. 6
figure 6

The a expected and b observed constraints on \(\kappa _c\) and \(\kappa _b\) at 68% and 95% confidence levels

8 Conclusion

A direct search for the decay of a Higgs boson to a charm quark–antiquark pair has been performed using \(139~\text{ fb}^{-1}\) of pp collision data recorded at by the ATLAS detector at the LHC. The search uses three channels, \(ZH \rightarrow \nu \nu c{\bar{c}}\), \(WH \rightarrow \ell \nu c{\bar{c}}\) and \(ZH \rightarrow \ell \ell c{\bar{c}}\). Signal events are identified using a multivariate charm tagging algorithm.

To enhance the signal sensitivity, events are categorised according to the \(p_{\text {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 \(VZ(\rightarrow c{\bar{c}})\) process and 3.8 (4.6) standard deviations for the \(VW(\rightarrow cq)\) process.

Fig. 7
figure 7

Expected and observed values of the combined \(VH, H\rightarrow c{\bar{c}}\) and \(VH, H\rightarrow b{\bar{b}}\) negative profile log-likelihood ratio as a function of \(\kappa _c/\kappa _b\), where \(\kappa _b\) is a free parameter. The vertical green lines correspond to the values of \(|\kappa _c / \kappa _b|\) for which the Higgs–charm and Higgs–bottom couplings are equal, where each coupling strength \(|\kappa _i y_i|\) is the product of the \(\kappa _i\) modifier and the Yukawa coupling, \(y_i\), for \(i = b,c\), and is equal to \(m_b / m_c = 4.578 \pm 0.008\) [119]

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 , 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 \(ZH, H\rightarrow c{\bar{c}}\), due to the larger dataset, improved c-tagging, and inclusion of the 0- and 1-lepton 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 \(|\kappa _c| < 8.5~(12.4)\), at the 95% confidence level.

A combination with the ATLAS \(H \rightarrow b{\bar{b}}\) measurement is performed. Interpreted in the kappa framework the combination constrains the observed ratio \(|\kappa _c / \kappa _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.