Measurements of $WH$ and $ZH$ production in the $H \rightarrow b\bar{b}$ decay channel in $pp$ collisions at 13 TeV with the ATLAS detector

Measurements of the Standard Model Higgs boson decaying into a $b\bar{b}$ pair and produced in association with a $W$ or $Z$ boson decaying into leptons, using proton-proton collision data collected between 2015 and 2018 by the ATLAS detector, are presented. The measurements use collisions produced by the Large Hadron Collider at a centre-of-mass energy of $\sqrt{s} = $13 TeV, corresponding to an integrated luminosity of 139 fb$^{-1}$. The production of a Higgs boson in association with a $W$ or $Z$ boson is established with observed (expected) significances of 4.0 (4.1) and 5.3 (5.1) standard deviations, respectively. Cross-sections of associated production of a Higgs boson decaying into bottom quark pairs with an electroweak gauge boson, $W$ or $Z$, decaying into leptons are measured as a function of the gauge boson transverse momentum in kinematic fiducial volumes. The cross-section measurements are all consistent with the Standard Model expectations, and the total uncertainties vary from 30% in the high gauge boson transverse momentum regions to 85% in the low regions. Limits are subsequently set on the parameters of an effective Lagrangian sensitive to modifications of the $WH$ and $ZH$ processes as well as the Higgs boson decay into $b\bar{b}$.


Introduction
signal and control regions, a significant increase in the effective number of simulated events and re-derived background modelling uncertainties, including using a multivariate approach to estimate the modelling uncertainty in the dominant backgrounds. A complementary analysis using the same final states, but focussing on regions of higher Higgs boson transverse momentum not accessible using the techniques outlined in this paper, has also been undertaken [36]. The same dataset was used, resulting in some overlap in the events analysed.

Data and simulated event samples
The data used in this analysis were collected using unprescaled single-lepton or missing transverse momentum triggers at a centre-of-mass energy of 13 TeV during the 2015-2018 running periods. Events are selected for analysis only if they are of good quality and if all the relevant detector components are known to have been in good operating condition, which corresponds to a total integrated luminosity of 139.0 ± 2.4 fb −1 [41,42]. The recorded events contain an average of 34 inelastic pp collisions per bunch-crossing.
Monte Carlo (MC) simulated events are used to model most of the backgrounds from SM processes and the V H, H → bb signal processes. A summary of all the generators used for the simulation of the signal and background processes is shown in Table 1. Samples produced with alternative generators are used to estimate systematic uncertainties in the event modelling, as described in Section 7. The same event generators as in Ref. [33] are used; however, the number of simulated events in all samples has been increased by at least the factor by which the integrated luminosity grew compared to the previous publication (∼ 1.75). In addition, processes which significantly contributed to the statistical uncertainty of the background in the previous publication benefited from a further factor of two increase in the number of simulated events produced.
All simulated processes are normalised using the most accurate theoretical cross-section predictions currently available and were generated to next-to-leading-order (NLO) accuracy at least, except for the gg → Z H and gg → VV processes, which were generated at LO. All samples of simulated events were passed through the ATLAS detector simulation [43] based on G [44]. The effects of multiple interactions in the same and nearby bunch crossings (pile-up) were modelled by overlaying minimum-bias events, simulated using the soft QCD processes of P 8.186 [45] with the A3 [46] set of tuned parameters (tune) and NNPDF2.3LO [47] parton distribution functions (PDF). For all samples of simulated events, except for those generated using S [48], the E G 1.6.0 program [49] was used to describe the decays of bottom and charm hadrons. Table 1: The generators used for the simulation of the signal and background processes. Samples are generated considering decays into all three lepton ( ) flavours. If not specified, the order of the cross-section calculation refers to the expansion in the strong coupling constant (α S ). The acronyms ME, PS and UE stand for matrix element, parton shower and underlying event, respectively. ( ) The events were generated using the first PDF in the NNPDF3.0NLO set and subsequently reweighted to the PDF4LHC15NLO set [50] using the internal algorithm in P -B 2. ( †) The NNLO(QCD)+NLO(EW) cross-section calculation for the pp → Z H process already includes the gg → Z H contribution. The qq → Z H process is normalised using the cross-section for the pp → Z H process, after subtracting the gg → Z H contribution. An additional scale factor is applied to the qq → V H processes as a function of the transverse momentum of the vector boson, to account for electroweak (EW) corrections at NLO. This makes use of the V H differential cross-section computed with H [51,52]. Contributions from photon-induced processes are also included for pp → W H [53]. ( ‡) For the diboson samples the cross-sections are calculated by the Monte Carlo generator at NLO accuracy in QCD. to have p T > 7 GeV and |η| < 2.47, to have small impact parameters,3 to fulfil a loose track isolation requirement, and to meet a 'LooseLH' quality criterion computed from shower shape, track quality and track-cluster matching variables. In the 1-lepton channel, tight electrons are selected using a 'TightLH' likelihood requirement and a calorimeter-based isolation in addition the the track-based isolation.
Muons are required to be within the acceptance of the muon spectrometer |η| < 2.7, to have p T > 7 GeV, and to have small impact parameters. Loose muons are selected using a 'loose' quality criterion [88] and a loose track isolation requirement. In the 1-lepton channel, tight muons fulfil the 'medium' quality criterion and a stricter track isolation requirement.
Hadronically decaying τ-leptons [89,90] are required to have p T > 20 GeV and |η| < 2.5, to be outside the transition region between the barrel and endcap electromagnetic calorimeters 1.37 < |η| < 1.52, and to meet a 'medium' quality criterion [90]. Reconstructed hadronic τ-leptons are not directly used in the event selection, but are utilised in the missing transverse momentum calculation and are also used to avoid double-counting hadronic τ-leptons as other objects.
Jets are reconstructed from the energy in topological clusters of calorimeter cells [91] using the anti-k t algorithm [92] with radius parameter R = 0.4. Jet cleaning criteria are used to identify jets arising from non-collision backgrounds or noise in the calorimeters [93], and events containing such jets are removed.
Jets are required to have p T > 20 GeV in the central region (|η| < 2.5), and p T > 30 GeV outside the tracker acceptance (2.5 < |η| < 4.5). A jet vertex tagger [94] is used to remove jets with p T < 120 GeV and |η| < 2.5 that are identified as not being associated with the primary vertex of the hard interaction. Simulated jets are labelled as b-, cor light-flavour jets according to which hadrons with p T > 5 GeV are found within a cone of size ∆R = 0.3 around their axis [95]. In the central region, jets are identified as b-jets (b-tagged) using a multivariate discriminant [95] (MV2), with the selection tuned to produce an average efficiency of 70% for b-jets in simulated tt events, which corresponds to light-flavour (u-, d-, s-quark and gluon) jet and c-jet misidentification efficiencies of 0.3% and 12.5% respectively.
Simulated V+jets events are categorised according to the two b-tagged jets that are required in the event: V + ll when they are both light-flavour jets, V + cl when there is one c-jet and one light-flavour jet, and V + HF (heavy flavour) in all other cases (which after the b-tagging selection mainly consist of events with two b-jets).
In practice, b-tagging is not applied directly to simulated events containing light-flavour jets or c-jets, because the substantial MV2 rejection results in a significant statistical uncertainty for these background processes. Instead, all events with c-jets or light-flavour jets are weighted by the probability that these jets pass the b-tagging requirement [87]. This is an expansion of the weighting technique compared to the previous analysis, where only jets in the V + ll, V + cl and WW processes were treated in this manner. Applying the same treatment to all light-flavour jets and c-jets significantly increases the number of simulated events present after the full event selection, reducing the statistical uncertainty of the V + HF (tt) background by ∼ 65%-75% (∼ 25%). When comparing the direct application of the b-tagging to the weighting technique, differences were observed in a particular subset of events with a small angular separation between the jets, but it was verified that this has a negligible impact on the result.
In addition to the standard jet energy scale calibration [96], b-tagged jets receive additional flavour-specific corrections to improve their energy measurement (scale and resolution): if any muons are found within a p T -dependent cone around the jet axis, the four-momentum of the closest muon is added to that of the jet. In addition, a residual correction is applied to equalise the response to jets with leptonic or hadronic decays of heavy-flavour hadrons and to correct for resolution effects. This improves the resolution of the dijet mass by up to ∼ 20% [87]. Alternatively, in the 2-lepton channel for events with two or three jets, a per-event kinematic likelihood uses the complete reconstruction of all final-state objects to improve the estimate of the energy of the b-jets. This improves the resolution of the dijet mass by up to ∼ 40%.
The missing transverse momentum, E miss T , is reconstructed as the negative vector sum of the transverse momenta of leptons, photons, hadronically decaying τ-leptons and jets, and a 'soft-term', p miss,st T . The soft-term is calculated as the vectorial sum of the p T of tracks matched to the primary vertex but not associated with a reconstructed lepton or jet [97]. The magnitude of E miss T is referred to as E miss T . The track-based missing transverse momentum, p miss T , is calculated using only tracks reconstructed in the inner tracking detector and matched to the primary vertex.
An overlap removal procedure is applied to avoid any double-counting between leptons, including hadronically decaying τ-leptons, and jets.

Event selection and categorisation
Events are categorised into 0-, 1-and 2-lepton channels (referred to as the n-lepton channels) depending on the number of selected electrons and muons, to target the Z H → vvbb, W H → νbb and Z H → bb signatures, respectively. In all channels, events are required to have exactly two b-tagged jets, which form the Higgs boson candidate. At least one b-tagged jet is required to have p T greater than 45 GeV. Events are further split into 2-jet or 3-jet categories depending on whether any additional, untagged jets are present. In the 0-and 1-lepton channels, only one such jet is allowed, as the tt background is much larger in events with four jets or more. In the 2-lepton channel any number of jets is accepted in the 3-jet category, which increases the signal acceptance in this category by 100%.
The reconstructed transverse momentum of the vector boson, p V T , corresponds to E miss T in the 0-lepton channel, the vectorial sum of E miss T and the charged-lepton transverse momentum in the 1-lepton channel, and the transverse momentum of the 2-lepton system in the 2-lepton channel. Since the signal-to-background ratio increases for large p V T values [98,99], the analysis focuses on two high-p V T regions defined as 150 GeV < p V T < 250 GeV and p V T > 250 GeV. In the 2-lepton channel, an additional fiducial measurement region is studied via the inclusion of a medium-p V T region with 75 GeV < p V T < 150 GeV. The event selection for the three lepton channels is outlined in Table 2 with details provided below.

0-lepton channel
The online selection uses E miss T triggers with thresholds that varied from 70 GeV to 110 GeV between the 2015 and 2018 data-taking periods. Their efficiency is measured in W+jets, Z+jets and tt events using single-muon triggered data, which effectively selects events with large trigger-level E miss T values as muons are not included in the trigger E miss T calculation. The resulting trigger correction factors that are applied to the simulated events range from 0.95 at the offline E miss T threshold of 150 GeV to a negligible deviation from unity at E miss T values above 200 GeV. A requirement on the scalar sum of the transverse momenta of the jets, H T , removes a small part of the phase space (less than 1%) where the trigger efficiency depends mildly on the number of jets in the event. Events with any loose lepton are rejected. High E miss T in multi-jet events typically arises from mismeasured jets in the calorimeters. Such events are efficiently removed by requirements on the angular separation of the E miss T , jets, and p miss T .

1-lepton channel
In the electron sub-channel, events are required to satisfy a logical OR of single-electron triggers with p T thresholds that started at 24 GeV in 2015 and increased to 26 GeV in 2016-2018. 4 The muon sub-channel uses the same E miss T triggers and correction factors as the 0-lepton channel. As these triggers effectively select on p V T , given that muons are not included in the trigger E miss T calculation, they perform more efficiently than the single-muon triggers in the analysis phase space, which have a lower efficiency due to the more limited coverage of the muon trigger system in the central region. Events are required to have exactly one tight muon with p T > 25 GeV or one tight electron with p T > 27 GeV and no additional loose leptons. In the electron sub-channel an additional selection of E miss T > 30 GeV is applied to reduce the background from multi-jet production.
Control regions High and low ∆R(b 1 , b 2 ) side-bands

2-lepton channel
The trigger selection in the electron sub-channel is the same as in the 1-lepton channel.
In the muon sub-channel, an OR of single-muon triggers is used, with lowest p T thresholds increasing from 2016-2018 and ranging from 20 GeV to 26 GeV. Events must have exactly two same-flavour loose leptons, one of which must have p T > 27 GeV, and the invariant mass of the lepton pair must be close to the Z boson mass.

Signal and control regions
The three n-lepton channels, two jet categories and two (0-lepton, 1-lepton) or three (2-lepton) p V T regions result in a total of 14 analysis regions. Each analysis region is further split into a signal region (SR) and two control regions (CRs), resulting in a total of 42 regions. The CRs are enriched in either V + HF or tt events and defined using a continuous selection on the ∆R between the two b-tagged jets, ∆R(b 1 , b 2 ), as a function of p V T , with the b-tagged jets labelled in decreasing p T as b 1 and b 2 . A lower and upper requirement on ∆R(b 1 , b 2 ) is applied, creating two CRs, referred to as the low and high ∆R CRs, shown in Figure 1. In the 1-lepton channel, the high ∆R selection was tuned such that the SR and low ∆R CR contain 95% (85%) of the signal in the 2-jet (3-jet) categories, whilst the low ∆R selection was tuned such that the SR contains 90% of the diboson yield, to ensure that a sufficient number of these events remain when conducting the diboson validation analysis. The same ∆R selection is applied in all three n-lepton channels and keeps over 93% of the signal in the 2-jet categories and over 81% (68%) of the signal in the 3-jet (≥ 3-jet) categories.5 The acceptances in the three n-lepton channels after the event selection, as well as the predicted crosssections times branching fractions for (W/Z)H with W → ν, Z → , Z → νν, and H → bb are given in Table 3. The non-negligible acceptance for the qq → W H process in the 0-lepton channel is mostly due to events with a hadronically decaying τ-lepton produced in the W decay, which are not explicitly vetoed and which could also be misidentified as a jet or subsequently decay to a low-p T electron or muon that fails to satisfy the selection criteria. The larger acceptance for the gg → Z H process compared with qq → Z H is due to the harder p V T spectrum of the gluon-induced process.

Simplified template cross-section categories
Cross-section measurements are conducted in the reduced V H, V → leptons stage-1.2 STXS region scheme [100, 101] described in Ref. [35] and summarised in Table 4. In this scheme, qq → Z H and gg → Z H are treated as a single Z H process, since there is currently not enough sensitivity to distinguish between them. The expected signal distributions and acceptance times efficiencies for each STXS region are estimated from the simulated signal samples by selecting events using information from the generator's 'truth' record, in particular the truth p V T , denoted by p V , t T . The signal yield in each reconstructed-event category for each STXS region is shown in Figure 2(a), with the corresponding fraction of signal events shown in Figure 2(b). The key improvement compared to the previous publication is the addition of a reconstructed-event category with p V T > 250 GeV. This region is more aligned with the STXS regions and significantly reduces the correlation between the STXS measurements in the two highest p V , t T bins. The acceptance times efficiency for W H events with p W , t T < 150 GeV or Z H events with p Z, t T < 75 GeV is at the level of 0.1% or smaller. Given the lack of sensitivity to these regions, the signal cross-section in these regions is constrained to the SM prediction, within the theoretical uncertainties. These regions contribute only marginally to the selected event sample and the impact on the final results is negligible.  aims to separate the V Z, Z → bb diboson process from the V H signal and other background processes, is used to validate the V H analysis. In each set, BDTs are trained in eight regions, obtained by merging some of the 14 analysis regions. In particular, the 150 GeV < p V T < 250 GeV and p V T > 250 GeV analysis regions in each lepton channel and jet category are merged for the training, as no increase in sensitivity was found when undertaking separate trainings in the two regions. The outputs of the BDTs, evaluated in each signal region, are used as final discriminating variables.
The BDT input variables used in the three lepton channels are detailed in Table 5. The separation of two b-tagged jets in pseudorapidity is denoted by |∆η(b 1 , b 2 )|. In 3-jet events, the third jet is labelled as jet 3 and the mass of the 3-jet system is denoted m bb j . The azimuthal angle between the vector boson and the system of the Higgs boson candidate formed from the two b-tagged jets is denoted ∆φ(V, bb), and their pseudorapidity separation is denoted ∆η(V, bb). In the 0-lepton channel, m eff is defined as the scalar sum of the transverse momenta of all jets and the E miss T (m eff = H T + E miss T ). In the 1-lepton channel, the angle between the lepton and the closest b-tagged jet in the transverse plane is denoted min(∆φ( , b)) and two variables are used to improve the rejection of the tt background: the rapidity difference between the W and Higgs boson candidates, |∆y(V, bb)| and, assuming that the event is tt, the reconstructed top quark mass, m top . The latter is calculated as the invariant mass of the lepton, the reconstructed neutrino and the b-tagged jet that yields the lower mass value. For both variables, the transverse component of the neutrino momentum is identified with E miss T , and the longitudinal component is obtained by applying a W-mass constraint to the lepton-neutrino system. The variable E miss T / √ S T , where S T is the scalar sum of transverse momenta of the charged leptons and jets in the event, is defined for use in the 2-lepton channel.
In addition to the above, which were all used in the previous iteration of the analysis [33], the following variables are also input to the BDTs: • Binned MV2 b-tagging discriminant: The MV2 discriminant for the two b-tagged jets is input to the BDT. The MV2 discriminant is grouped into two bins corresponding to efficiencies of 0-60% and 60%-70%, which are calibrated to data [95, 102,103]. This variable provides additional rejection against backgrounds where a c-jet or light-flavour jet has been misidentified as a b-jet, especially W → cq in the tt and Wt backgrounds. This improves the sensitivity in the 1-lepton (0-lepton) channel by ∼ 10% (∼ 7%). The binned MV2 discriminant does not provide any additional sensitivity in the 2-lepton channel, where the backgrounds are dominated by processes containing two b-jets.
• Magnitude of the track-based E miss T soft-term, p miss,st T : In the 0-lepton channel this provides additional rejection against the tt background, which may contain unreconstructed objects, such as leptons or b-jets, due to kinematic and detector acceptance. The presence of such objects in an event will result in a larger p miss,st T for tt events than for signal events. This improves the sensitivity in the 0-lepton channel by ∼ 2%-3%.
• Z boson polarisation, cos θ( − , Z): The cos θ( − , Z) is calculated as the cosine of the polar angle between the lepton ( − ) direction in the Z rest frame and the flight direction of the Z boson in the laboratory frame. The Z bosons from the Z H signal process are expected to have a different polarisation compared to those from the dominant Z+jets background [104], which provides additional background rejection in the 2-lepton channel. This improves the sensitivity in the 2-lepton channel by ∼ 7%.
The distributions of all input variables of the BDTs are compared between data and simulation, and good agreement is found within the uncertainties. The same training procedures and BDT output binning transformation as those detailed in Ref. [33] are used, with the exception that the training algorithm was updated to use gradient boosting in the TMVA [105] framework. Variable 0-lepton 1-lepton 2-lepton

Background modelling
The simulated event samples summarised in Section 3 are used to model all background processes, except for the tt background in the 2-lepton channel6 and the multi-jet background in the 1-lepton channel, which are both estimated using data-driven techniques, as discussed below.

Data-driven tt background estimation
In the 2-lepton channel a high-purity control region, over 99% pure in tt and single-top-quark Wt events (jointly referred to as the top background), is defined using the nominal event selection, but replacing the same-flavour lepton selection with a requirement of exactly one electron and one muon. This region is referred to as the eµ-control region, eµ-CR. As these top background events typically contain two W bosons which decay into leptons, they are symmetric in lepton flavour. The events in the eµ-CR are directly used to model the shape and normalisation of the same-flavour lepton top background in the nominal selection. Any bias caused from the lepton trigger, reconstruction, identification or acceptance, is determined by comparing the yield of simulated top background events in the nominal selection with that in the eµ control region. No significant bias in the shape or normalisation is observed for any of the important kinematic variables, including the BDT discriminant. A ratio of the top yield in the analysis region to that in the eµ-CR of 1.00 ± 0.01 (1.01 ± 0.01) is determined using simulation, for the 2-jet (≥ 3-jet) region, where the uncertainty in the ratio is the statistical uncertainty resulting from the simulated samples. As no evaluated theoretical or experimental uncertainties create any bias beyond the statistical uncertainty of the ratio, the latter is assigned as an extrapolation uncertainty. This method has the advantage that all the experimental and theoretical uncertainties are eliminated, resulting in the data statistics in the eµ-CR becoming the dominant uncertainty source for the data-driven top background estimate.

Multi-jet background estimation
Multi-jet (MJ) event production has a large cross-section and thus, despite not being a source of genuine missing transverse momentum or prompt leptons, has the potential to contribute a non-negligible amount of background. Using the same techniques detailed in Ref. [33], the MJ background was demonstrated to be negligible in both the 0-and 2-lepton channels.
In the 1-lepton channel, the MJ background is reduced to the percent level and is predicted using the same method as described in Ref. [33] with minor changes to account for the use of the MV2(b j ) variables in the BDT. The MJ background is modelled from data in an MJ-enriched control region (MJ-CR), from which all simulated backgrounds are subtracted. The MJ-CR is defined by applying the nominal event selection, except for the stricter lepton isolation requirement, which is inverted. The requirement on the number of b-tagged jets is relaxed from two (2-b-tag MJ-CR) to one (1-b-tag MJ-CR) to increase the statistical precision. To correctly estimate the 2-b-tag MJ BDT shape, the values of both the MV2(b 1 ) and MV2(b 2 ) BDT input variables in the 1-b-tag events, are replaced with values emulated from a joint MV2(b 1 ) and MV2(b 2 ) probability distribution derived from the 2-b-tag MJ-CR. The normalisation of the MJ background is then determined from a template fit to the m W T distribution after applying the nominal selection with a 2-b-tag requirement, using the MJ shape predicted from the 1-b-tag MJ-CR and the shapes of the other backgrounds from simulation.

Systematic uncertainties
The sources of systematic uncertainty can be broadly divided into three groups: those of an experimental nature, those related to the modelling of the backgrounds and those associated with the Higgs boson signal simulation. The estimation of the uncertainties closely follows the methodology outlined in Refs. [35,87] and is briefly summarised below.

Experimental uncertainties
The dominant experimental uncertainties originate from the b-tagging correction factors, jet energy scale calibration and the modelling of the jet energy resolution. The b-tagging correction factors, determined from the difference between the efficiencies measured in data and simulation, are evaluated in five MV2 discriminant bins and are derived separately for b-jets, c-jets and light-flavour jets [95, 102,103]. All of the correction factors for the three jet flavours have uncertainties estimated from multiple measurements, which are decomposed into uncorrelated components that are then treated independently. After removal of components that have a negligible impact, the total numbers of independent systematic components are 29 uncertainties for b-jets, 18 for c-jets and 10 for light-flavour jets. The uncertainties in the jet energy scale and resolution are based on their respective measurements [96,106]. The various sources of uncertainty in the calibration of the jet energy scale are characterised by 30 independent components and the jet energy resolution is characterised by 8 independent components. An additional specific uncertainty in the energy calibration of band c-jets is also considered.
Uncertainties in the reconstruction, identification, isolation and trigger efficiencies of muons [88] and electrons [107] are considered, along with the uncertainty in their energy scale and resolution. These are found to have only a small impact on the result. The uncertainties in the energy scale and resolution of the jets and leptons are propagated to the calculation of E miss T , which also has additional uncertainties from the modelling of the underlying event and momentum scale, momentum resolution and reconstruction efficiency of the tracks used to compute the soft-term [97,108]. An uncertainty is assigned to the E miss T trigger correction factors, determined from the ratio of the trigger efficiency in data and simulation, to account for the statistical uncertainty in the measured correction factors and for differences between the correction factors determined from W + jets, Z + jets and tt events. The uncertainty in the combined 2015-2018 integrated luminosity is 1.7%. It is derived following a methodology similar to that detailed in Ref.
[41], and using the LUCID-2 detector for the baseline luminosity measurements [42]. The average number of interactions per bunch crossing in the simulation is rescaled by 1.03 to improve agreement between simulation and data, based on the measurement of the visible cross-section in minimum-bias events [109], and an uncertainty, as large as the correction, is included.

Background uncertainties
Modelling uncertainties are derived for the simulated samples and broadly cover three areas: normalisations (referred to as normalisation uncertainties), acceptance differences that affect the relative normalisations between regions with a common underlying normalisation (referred to as relative acceptance uncertainties), and the shapes of the differential distributions of the kinematic variables (referred to as shape uncertainties).
The overall cross-sections and associated normalisation uncertainties for the background processes are taken from the currently most accurate calculations as detailed in Table 1, apart from the main backgrounds (Z + HF, W + HF, tt) whose normalisations are left unconstrained (floated) in the global likelihood fit.
The relative acceptance and shape uncertainties are derived from either particle-level or reconstruction-level comparisons between nominal and alternative simulated samples, or from comparisons with data in control regions. The alternative samples are produced either by different generators or by altering the nominal generator's parameter values. When relative acceptance uncertainties are estimated, the nominal and alternative samples are normalised using the same production cross-section. Shape uncertainties are estimated within a signal region, an analysis region or a set of analysis regions, depending on the distribution being varied, with the nominal and alternative samples scaled to have the same normalisation in the studied area. When the shape uncertainty acts over an analysis region or a set of analysis regions, it also controls the migration between the component regions via the relative acceptance changes from the shape variation (referred to as a shape plus migration uncertainty), as opposed to a shape uncertainty which only alters the shape within a single SR (referred to as just a shape uncertainty). Unless stated otherwise, the uncertainty is taken from the alternative sample that differs most in shape from the nominal sample.
Shape uncertainties for Z + HF, single-top and diboson backgrounds are derived for the m bb and p V T variables, as it was found sufficient to consider the changes induced in these variables to cover the overall shape variation of the BDT discriminant. For W + HF and tt backgrounds, a more sophisticated multidimensional parameterisation method is introduced to estimate the shape uncertainties of the final discriminant [110]. In this method, a BDT (referred to as BDT S ) is trained to discriminate the nominal sample from an alternative sample, using the kinematic variables from the BDT V H (Table 5) as input variables, except for the p V T . Before training, the p V T distribution of the nominal sample is reweighted to match that of the alternative sample. The p V T difference is considered as a separate, uncorrelated uncertainty, in a manner similar to that for the other backgrounds. The ratio of the BDT S distributions evaluated for the alternative and nominal samples provide a reweighting function (referred to as R BDT ), which can be used to correct the nominal sample to match the alternative sample. This method simultaneously maps the n-dimensional space formed by the kinematic variables of the two generators onto each other. It is verified that, after being reweighted by R BDT , the input variable distributions for the nominal sample are in good agreement with those of the alternative sample.
The systematic uncertainties affecting the modelling of the background samples are summarised in Tables 6  and 7, and key details of the treatment of the backgrounds are reported below.
V + jets production The V + jets backgrounds are subdivided into three different components based upon the jet flavour labels of the two b-tagged jets in the event. The main background contributions (V + bb, V + bc, V + bl and V + cc) are jointly considered as the V + HF background. Their overall normalisations are free to float in the global likelihood fit, separately in the 2-and 3-jet categories. For the Z + HF background, the normalisations are also floated separately in the 75 GeV < p V T < 150 GeV and p V T > 150 GeV regions. The remaining flavour components, V + cl and V + ll, constitute less than ∼ 1% of the background in each analysis region and only normalisation uncertainties are included.
Uncertainties are estimated for the relative normalisation of the four heavy-flavour components that constitute the V + HF background. These are taken as uncertainties in the bc, cc and bl yields compared with the dominant bb yield and are estimated separately in each lepton channel in a manner similar to the acceptance systematic uncertainties. Relative acceptance uncertainties for the W + HF background are estimated for the ratio of the event yield in the 0-lepton channel to that in the 1-lepton channel. For the Z + HF background, there is a relative acceptance uncertainty in the ratio of the event yield in the 0-lepton channel to that in the 2-lepton channel in the p V T > 150 GeV region. For both W + HF and Z + HF, relative acceptance uncertainties are estimated for the ratio of the event yield in the SR to that in the CRs.
For Z + HF, shape uncertainties are derived for m bb and p V T , which are evaluated from comparisons with data in the m bb side-bands (m bb < 80 GeV or m bb > 140 GeV), after subtracting backgrounds other than Z + jets. For W + HF, uncertainties are derived for p V T and the R BDT method from comparisons of the nominal sample (S ) with an alternative sample (M G 5_aMC@NLO+P 8 [111,112]).
t t production In the 0-and 1-lepton channels (jointly referred to as 0+1-lepton channel) separate floating normalisations are used for the 2-jet region and 3-jet region. Uncertainties are derived from comparisons between the nominal sample (P +P 8) and alternative samples corresponding to matrixelement (M G 5_aMC@NLO+P 8) and parton-shower (P +H 7 [113]) generator variations. Table 6: Summary of the systematic uncertainties in the background modelling for Z + jets, W + jets, tt, single-topquark and multi-jet production. 'ME' indicates a matrix element generator variation and 'PS' indicates a parton shower generator variation. An 'M+S' symbol is used when a shape uncertainty includes a migration effect that allows relative acceptance changes between regions, whilst 'S' indicates that the uncertainty only acts upon the shape in the signal region. Instances where an uncertainty is considered independently in different regions are detailed in parentheses. Where the size of an acceptance systematic uncertainty varies between regions, a range is displayed.
The dominant flavour component of the two b-tagged jets in tt is bb. However, there is a sizeable bc component which has a more signal-like topology. Uncertainties in the relative composition of three components, bb, bc, and any other flavour configuration (referred to as 'other') are estimated from the difference in the ratio of the bc or other components to the bb yield between the nominal sample and the alternative matrix element and parton shower generator samples. Shape uncertainties are derived for p V T and using the R BDT method in the 0+1-lepton channels from comparisons with the alternative parton shower and matrix element generator samples.
In the 2-lepton channel the tt background is estimated by a data-driven method as discussed in Section 6.1. The uncertainty in this background is dominated by the statistical uncertainty of the eµ control region data events.

Single-top-quark production
In the Wtand t-channels, uncertainties are derived for the normalisation, relative acceptance and shapes of the m bb and p V T distributions. For the Wt-channel, the estimated modelling uncertainties are applied independently according to the flavour of the two b-tagged jets, due to the different regions of phase space being probed when there are two b-jets (bb) present compared with events where there are fewer b-jets present (referred to as 'other'). Those uncertainties are evaluated from comparisons between the nominal sample (P +P 8 using the diagram removal scheme [114]) and alternative samples with parton-shower variations (P +H ++) and a different scheme to account for the interference between Wt and tt production (P +P 8 using the diagram subtraction scheme) [115]. Only a normalisation uncertainty is derived for the s-channel, since its contribution is at a very low level.

Diboson production
The diboson backgrounds are composed of three distinct processes: W Z, WW and Z Z production. Given the small contribution from WW production (< 0.1% of the total background) only a normalisation uncertainty is assigned. For the more important contributions from the W Z and Z Z backgrounds, uncertainties are considered in the overall normalisation, the relative acceptance between regions and the m bb and p V T shapes. These are derived following the procedure described in Ref. [87] and are outlined in Table 7, which includes comparisons of the nominal sample (S ) with alternative samples (P +P 8 and P +H ++).

Multi-jet background uncertainties
The systematic uncertainties in the multi-jet background estimate in the 1-lepton channel are derived by following the procedure outlined in Ref. [33]. Two different uncertainty components are considered, those which alter the normalisation and those which alter the multi-jet BDT template shape.

Signal uncertainties
The systematic uncertainties that affect the modelling of the signal are summarised in Table 8 and are estimated with procedures that closely follow those outlined in Ref. [27,35,116,117]. The systematic uncertainties in the calculations of the V H production cross-sections and the H → bb branching fraction7 are assigned following the recommendations of the LHC Higgs Cross Section Working Group [31,70,71,118,119].
Uncertainties in the m bb and p V T signal shape are estimated, as described in Ref. [33], from scale variations, PDF and α S (PDF+α S ) uncertainties, from varying the parton shower and underlying event (PS/UE) models using AZNLO tuning variations and from comparisons with alternative parton-shower generator samples (P +H 7). In addition, a systematic uncertainty from higher-order EW corrections effects is taken into account as a variation in the shape of the p V T distributions for qq → V H production. Acceptance uncertainties, evaluated according to STXS regions, correctly accounting for the migration and correlations between regions, are evaluated for the scale variations, PS/UE models and PDF+α S .
For the STXS measurement, the signal uncertainties are separated into two groups, uncertainties in the acceptance and shape of kinematic distributions which alter the signal modelling (theoretical modelling uncertainties) and the uncertainties in the prediction of the production cross-section for each of these regions (theoretical cross-section uncertainties). Whilst theoretical modelling uncertainties enter the STXS measurements, theoretical cross-section uncertainties only affect the predictions with which they are compared, and are therefore not included in the likelihood function.

Statistical analysis
The statistical procedure is based on a likelihood function L(µ, θ), constructed as the product of Poisson probability terms over the bins of the input distributions, with parameters of interest (POI) extracted by maximising the likelihood. The effects of systematic uncertainties enter the likelihood as nuisance parameters (NP), θ. Most of the uncertainties discussed in Section 7 are constrained with Gaussian or log-normal probability density functions. The normalisations of the largest backgrounds, tt, W + HF and Z + HF, can be reliably determined by the fit, so they are left unconstrained in the likelihood. The uncertainties due to the limited number of events in the simulated samples used for the background predictions are included using the Beeston-Barlow technique [120]. As detailed in Ref. [121], systematic variations that are subject to large statistical fluctuations are smoothed, and systematic uncertainties that have a negligible impact on the final results are pruned away region-by-region (treating signal and control regions separately).
The global likelihood fit comprises 14 signal regions, defined as the 2-and 3-jet categories in the two high-p V T (150 < p V T < 250 GeV and p V T > 250 GeV) regions for the three channels, and in the medium-p V T region (75 < p V T < 150 GeV) for the 2-lepton channel. The 28 control regions are also input as event yields in all fit configurations.
Three different versions of the analysis are studied, which differ in the distributions input to the fit.
• The nominal analysis, referred to as the multivariate analysis, uses the BDT V H multivariate discriminant output distributions as the inputs to the fit. Three different POI configurations are studied. Firstly, a single-POI fit measures µ bb V H , the signal strength that multiplies the SM Higgs boson V H production cross-section times the branching fraction into bb. Secondly, a two-POI fit is undertaken, which jointly measures the signal strengths of the W H and Z H components. Finally, a five-POI fit version measures the signal cross-section multiplied by the H → bb and V → leptons branching fractions in the five STXS regions (see Table 4).
• The dijet-mass cross-check analysis uses the m bb distributions, instead of the BDT V H distributions, as inputs to a single-POI fit to measure µ bb V H .
• The diboson validation analysis, a measurement of the signal strength of the W Z and Z Z processes, uses the BDT V Z output distributions. The SM Higgs boson is included as a background process normalised to the predicted SM cross-section with an uncertainty of 50%, which conservatively encompasses the previous measurement and uncertainty [33]. Two POI configurations are evaluated, firstly a single-POI fit to measure µ bb V Z , the signal strength of the combined W Z and Z Z diboson processes, and secondly a two-POI fit to simultaneously measure the W Z and Z Z signal strengths.
The background predictions in all post-fit distributions and tables are obtained by normalising the backgrounds and setting the nuisance parameters according to the values determined by the fit used to measure µ bb V H .

Signal strength measurements
The post-fit normalisation factors of the unconstrained backgrounds in the global likelihood fit are shown for the single-POI multivariate analysis in Table 9, the post-fit signal and background yields are shown in Tables 10 and 11, and Figure 3 shows the BDT V H output distributions in the high-p V T 2-jet SRs, which are most sensitive to the signal. For the V H production mode the background-only hypothesis is rejected with a significance of 6.7 standard deviations, to be compared with an expectation of 6.7 standard deviations [122].
The results of the combined fit when measuring signal strengths separately for the W H and Z H production processes are shown in Figure 4. The W H and Z H production modes reject the background-only hypothesis with observed (expected) significances of 4.0 (4.1) and 5.3 (5.1) standard deviations, respectively. The fitted values of the two signal strengths are:      Table 12. The impact of a set of systematic uncertainties is defined as the difference in quadrature between the uncertainty in µ computed when all NPs are fitted and that when the NPs in the set are fixed to their best-fit values. The total statistical uncertainty is defined as the uncertainty in µ when all the NPs are fixed to their best-fit values. The total systematic uncertainty is then defined as the difference in quadrature between the total uncertainty in µ and the total statistical uncertainty. For the W H and Z H signal strength measurements the total statistical and systematic uncertainties are similar in size, with the b-tagging, jet, E miss T , background modelling and signal systematic uncertainties all making important contributions to the total systematic uncertainty. The impact of the statistical uncertainty from the simulated event samples has been significantly reduced compared to the previous result [35], due to the measures taken to considerably enhance the number of simulated events.  Using the 'bootstrap' method [121], the dijet-mass and nominal multivariate analysis results are found to be statistically compatible at the level of 1.1 standard deviations. The observed excess rejects the background-only hypothesis with a significance of 5.5 standard deviations, compared to an expectation of 4.9 standard deviations. Good agreement is also found when comparing the values of signal strengths in the individual channels from the dijet-mass analysis with those from the multivariate analysis.
The m bb distribution is shown in Figure 5 summed over all channels and regions, weighted by their respective values of the ratio of fitted Higgs boson signal to background yields and after subtraction of all backgrounds except for the W Z and Z Z diboson processes.

Diboson validation
The measurement of V Z production using a multivariate approach, as a validation of the Higgs boson analysis, returns a signal strength of −0.14 = 0.93 +0.07 −0.06 (stat.) +0.14 −0.12 (syst.), in good agreement with the Standard Model prediction. Analogously to the nominal analysis, fits are also performed with separate signal strengths for the W Z and Z Z production modes, and the results are shown in Figure 6.

Cross-section measurements
The measured V H cross-sections times the H → bb and V → leptons branching fractions, σ × B, together with the SM predictions in the reduced STXS regions, are summarised in Table 13 and Figure 7. The cross-sections are all consistent with the Standard Model expectations and are measured with relative uncertainties varying from 30% in the highest p V T region to 85% in the lowest p V T region. The data statistical uncertainty is the largest single uncertainty in all regions, although in the lower p V T regions systematic uncertainties make a sizeable contribution to the total uncertainty. In all regions there are large contributions from the background modelling, b-tagging and jet systematic uncertainties. In the lowest p V T region in both the W H and Z H measurements, the E miss T uncertainty is one of the largest uncertainties. For the Z H measurements, the signal uncertainties also make a sizeable contribution due to the limited precision of the theoretical calculations of the gg → Z H process.

Constraints on effective interactions
The strength and tensor structure of the process V H, H → bb are investigated using an effective Lagrangian approach. Extra terms are added to the SM Lagrangian (L SM ) to obtain an effective Lagrangian (L SMEFT ) following the approach in Refs. [124,125]:  The STXS measurements are used to constrain the coefficients of the operators in the 'Warsaw' formulation [126], which provides a complete set of independent operators when considering those allowed by the SM gauge symmetries. Thirteen operators directly affect the V H cross-section [127]. This analysis has significant sensitivity to the six operators detailed in Table 14, in addition to the operator which directly affects the H → bb decay width.
Following methodologies similar to those outlined in Refs. [127,128], a parameterisation of the STXS production cross-section and Higgs boson decay rates in terms of the SMEFT parameters is derived, in this case based upon leading-order predictions made using the SMEFTsim package [125]. The interference terms between the SM and BSM amplitudes are linear in the coefficients and of order 1/Λ 2 , while BSM contributions are quadratic in the coefficients and of order 1/Λ 4 . Linear terms from D = 8 operators are suppressed by the same 1/Λ 4 factor as the quadratic D = 6 terms. However, it is currently not possible to include such terms, so results for both the linear and linear plus quadratic D = 6 terms are studied to provide some indication of the effect D = 8 linear terms could have on the result. Modifications of the gg → Z H production cross-section are only introduced by either higher-dimension (D ≥ 8) operators or corrections that are formally at NNLO in QCD, and are not included in this study. The expected gg → Z H contribution is fixed to the SM prediction within uncertainties. The dependence of the experimental acceptance in each analysis region on the Wilson coefficients is not accounted for in this study, although it was verified that the impact on the acceptance from the EFT operators was at most 10%.
Maximum-likelihood fits across the STXS regions are performed to determine the Wilson coefficients. All coefficients but one are assumed to vanish, and one-dimensional confidence level (CL) intervals are inferred for the coefficient under study both with and without the quadratic terms. An example negative-log-likelihood one-dimensional projection is shown in Figure 8 for c H q3 , and the 68% and 95% CL intervals are summarised in Figure 9 for the four coefficients to which the analysis has greatest sensitivity, in addition to the c dH coefficient which directly affects the H → bb decay width. As detailed in Table 14, the O Hu , O H d and O H q1 operators have a similar impact and as such are found to be highly degenerate, so only a representative result for O Hu is shown. The coefficient c H q3 is constrained at 68% CL to be no more than a few percent, whilst the constraints on the other three coefficients range from 10%-30% to order unity and c dH has much weaker constraints. In most cases the observed constraints are found to significantly depend on the presence of the quadratic terms, indicating that D = 8 linear terms could also have a non-negligible effect.
These limits were also produced using the full likelihood and using only the STXS measurement central values and covariance matrix. It was found that the two methods produced results that are consistent with each other within ∼ 10%-20% for the majority of operators and to within ∼ 30% for the two operators with the weakest constraints, O dH and O HW B . As there are only five STXS regions, attempting to simultaneously extract constraints on multiple coefficients, some of which have similar effects, leads to unmanageable correlations. An alternative approach is to fit an orthogonal set of linear combinations of the Wilson coefficients of the Warsaw-basis operators. This removes the assumption, inherent in the one-dimensional limits, that only one operator acts at a time. Based upon the procedure outlined in Ref. [127], eigenvectors are determined from the Hessian matrix of the STXS likelihood fit to data, after it has been re-expressed in terms of the Wilson coefficients. This approach only considers the linear terms and the H → bb partial width, with a dedicated independent parameter added to account for the modifications to the total width. The resulting five eigenvectors are shown in Table 15. They are labelled as E0-E4 and ordered in terms of experimental sensitivity, with E0 having the greatest and E4 the least. The eigenvectors contain information about the sensitivity of the analysis to degenerate deformations of the SM. The leading eigenvector, E0, consists almost exclusively of c H q3 , which is also the coefficient most constrained in the one-dimensional limits, with similar limits obtained in both cases. The second eigenvector, E1, is dominated by c Hu , but has sizeable contributions from c H d and c H q1 , suggesting only a linear combination of these coefficients can be constrained given the degeneracy between them. The eigenvector E2 demonstrates sensitivity to a combination of the branching ratio and c HW , whilst E3 has limited sensitivity to a combination of c HW B and c H q1 . The analysis has negligible sensitivity to the fifth eigenvector. Figure 10 shows the impact on the STXS cross-section measurements when varying the coefficients for the four leading eigenvectors within their 1σ bounds. The analysis has greatest sensitivity to coefficients which predominantly increase the cross-section in the higher p V T STXS regions (E0 and E1), with lower sensitivity to those which predominantly impact the lower p V T STXS regions (E2 and E3).  [126]. All modifications that alter the branching ratio are absorbed into an additional independent term (∆BR/BR SM ), which linearly alters the branching ratio and all contributions with a coefficient below 0.2 are omitted. The full composition of the eigenvectors is available in the HEPData repository [123].  Figure 10: The impact of the leading four eigenvectors on the STXS cross-section measurements. The change to the cross-section is indicated at the +1σ (solid) and −1σ (dashed) limits of the corresponding Wilson coefficients, extracted from a simultaneous fit to data of all five eigenvectors.

Conclusion
Measurements are presented of the Standard Model Higgs boson decaying into a bb pair and produced in association with a W or Z boson, using data collected by the ATLAS experiment in proton-proton collisions from Run 2 of the LHC. The data correspond to an integrated luminosity of 139 fb −1 collected at a centre-of-mass energy of Cross-sections of associated production of a Higgs boson decaying into bottom quark pairs and an electroweak gauge boson, W or Z, decaying into leptons are measured as a function of the gauge boson transverse momentum in kinematic fiducial volumes in the simplified template cross-section framework. The uncertainties in the measurements vary from 30% in the highest p V T regions to 85% in the lowest, and are in agreement with the Standard Model predictions.
Limits are also set on the coefficients of effective Lagrangian operators which affect the V H production and H → bb decay. Limits are studied for both the variation of a single coefficient and also the simultaneous variation of a set of linear combinations of coefficients. The allowed range of the individual or linear combinations of the coefficients, to which the analysis has the greatest sensitivity, is limited to a few percent.                      The ATLAS Collaboration