The correlation matrix of Higgs rates at the LHC

The imperfect knowledge of the Higgs boson decay rates and cross sections at the LHC constitutes a critical systematic uncertainty in the study of the Higgs boson properties. We show that the full covariance matrix between the Higgs rates can be determined from the most elementary sources of uncertainty by a direct application of probability theory. We evaluate the error magnitudes and full correlation matrix on the set of Higgs cross sections and branching ratios at $\sqrt{s}=7$, $8$, $13$ and $14$ TeV, which are provided in ancillary files. The impact of this correlation matrix on the global fits is illustrated with the latest $7$+$8$ TeV Higgs dataset.


Introduction
The measurement of the production and decay rates of the Higgs boson at the LHC [1,2] has provided a new set of precision tests for the Standard Model (SM) of particle physics, giving access to the couplings of the Higgs to the other particles of the SM. Such information on the Higgs couplings can in turn be recast into constraints on theories extending the SM. The SM Higgs couplings being all fixed (now that there is a measurement of the Higgs boson mass), any significant deviation would reveal the existence of a new physics lying beyond the Standard Model.
This characterisation of the Higgs couplings relies on the comparison between the observed event numbers and those predicted by the SM in a given detection channel. These theoretical predictions of the SM event rates are not infinitely precise. Rather, they are plagued by a number of uncertainties, including the Parton Distribution Functions (PDF) uncertainties and the error inherent to the perturbative computation of amplitudes for the partonic processes.
The uncertainties on the LHC Higgs production rates can reach about 10% in relative magnitude. This is enough to substantially influence the outcome of the Higgs fits, so that the effect of these uncertainties should be carefully taken into account. As a matter of fact, the uncertainties on the cross sections have been found to have a non-negligible impact already on the fits of the 7 + 8 TeV data (see e.g. Ref. [3]). Besides, the statistical uncertainties will decrease with more data becoming available, hence the uncertainties on Higgs rates predictions are expected to become more and more important in the upcoming Higgs analyses.
The uncertainties on the SM Higgs production cross sections and partial decay widths are frequently taken into account in the global fits available in the literature. However, the correlations among these uncertainties are poorly known, and are usually approximated or neglected. For example in the case of the 7+8 TeV CMS-ATLAS global fit of the Higgs couplings [3], the correlation matrix among cross sections is approximated with either 0%, 100% or −100% elements. A similar approximation has been adopted in Ref. [4]. Yet, these Higgs rates correlations can in principle have a substantial impact on the fit [4], so that a complete treatment of the correlations is desirable.
The present paper is dedicated to address this shortcoming of the Higgs analyses, by providing a complete set of uncertainties on the Higgs rates at the LHC, including the correlations among all the rates. Our study covers the uncertainties on inclusive cross sections and branching ratios. No LHC-related experimental systematic uncertainties (on efficiencies or luminosity for example) are included.
We first show in Section 2 that the complete correlation matrix of the Higgs rates can be directly obtained from a straightforward application of probability theory to the set of elementary sources of uncertainty. A necessary condition for this approach to work is that the relative magnitude of the elementary uncertainties be 100% in order to allow error propagation, which will turn out to be a very good approximation. The subsequent formalism for error combinations is described in Section 3. In Section 4, we then enumerate all the elementary sources of errors affecting the LHC Higgs rates predictions and estimate their magnitudes, including the error on the Higgs mass determination. The full correlation matrix and combined errors on the Higgs rates are given in Section 5. These results are independent of the distributions associated to the elementary uncertainties. Finally in Section 6, we derive the analytical marginal likelihood including all correlations, using the fact that a central limit theorem is at work in presence of the large number of independent Higgs uncertainties [5]. Fits of the full 7+8 Higgs dataset are also performed.

Covariance and correlation matrix of expected Higgs rates
The likelihood associated to the Higgs events has in general the form P ( Higgs data | θ, δ) ≡ L(θ, δ) . (2.1) Here θ = (θ 1 , θ 2 , . . .) denotes the set of parameters of interest, which in the context of Higgs fits typically parametrises the deviations of the observed Higgs rates from the Standard Model predictions. The nuisance parameters δ = (δ 1 , δ 2 , . . .) model any kind of systematic uncertainties. This paper is focused on the systematic uncertainties on the inclusive Higgs rates, i.e. the SM predictions for production cross sections and decay widths with no phase space cuts. 1 By definition, these uncertainties enter in the likelihood via the cross sections σ X and the decay widths Γ Y , so that the Higgs likelihood has in general the form Here and below, we do not display the experimental nuisance parameters, which are irrelevant for this work. Inferring information about the θ parameters necessitates to marginalise the unwanted parameters δ, either through an integration in case of a Bayesian treatment, or a maximisation (i.e. profiling) in case of a frequentist treatment, see Ref. [4] for a review. In either case, a "prior" distribution π(δ) is associated to the nuisance parameters. In practice, there can be a large number of nuisance parameters δ, so that marginalisation methods can be technically difficult to perform. However, in the situation of small relative magnitudes (see Eq.(3.1)), a simplification is available: it is in principle possible to combine all sources of uncertainty together before marginalising.
Such combinations of uncertainties in both Bayesian and frequentist frameworks have been discussed in Ref. [4], where it has been found that Bayesian combinations are both better defined and simpler than the frequentist ones. In the particular case of Gaussian priors, Bayesian and frequentist combinations are found to be equivalent. In the Bayesian case, the combination is defined by where the new dimensionful δσ X , δΓ Y nuisance parameters encapsulate the combination of all the elementary sources of uncertainty, and σ 0 X and Γ 0 Y are the nominal values when no uncertainty is taken into account. The practical interest of this combination is that the dimension of the integral on the right-handed side is lower than on the left-handed side, rendering the operation of marginalisation more economic to achieve.
The key point to correctly apply the method of preliminary combination lies in the determination of the combined priorπ. This might seem a very challenging task at first sight, as theπ prior arises from a complicated combination of all elementary priors, and in general does not factorise. However, the covariance matrix generated byπ, which is a central object for inference, is independent of the shape of the combined prior. Moreover it is approximately independent of the elementary priors shape if elementary uncertainties have small relative magnitude (see Sec. 3). 2 The covariance matrix arising fromπ can thus be determined without further reference to prior shapes. This is what we are going to compute throughout this paper. Moreover, in case of a somewhat large number of sources of uncertainty with similar magnitudes, the shape ofπ automatically tends to a multivariate Gaussian by virtue of the Lyapunov central-limit theorem (see Ref. [5]). In this caseπ is completely determined by the covariance matrix and the mean value -the latter will always be set to zero with no loss of generality in our conventions. Finally, the frequentist correlation matrix matches the Bayesian one in case of priors with Gaussian shapes.
The goal of this work is to provide an accurate and consistent determination of the covariance matrix generated byπ, and the distribution of uncertainties on expected Higgs rates. In our analysis, for the combined uncertainty on decay rates, we are going to work at the level of branching ratios B Y instead of partial widths. The error on partial widths are propagated into branching ratios following For further purpose, it is convenient to put the nuisance parameters under a standardised form, where the normalised nuisance parameters δ X, , and the numbers ∆ 2 X , ∆ 2 Y correspond to the relative variances for σ X and B Y . The interest of this notation is that the relative magnitudes of the uncertainties are explicitly factored out. In turn, the covariance matrix Cov[δ X , δ X ] is precisely the correlation matrix between the production modes. This applies similarly to the decay modes. 3 2 The covariance matrix of the combined uncertainties is the second central moment of theπ distribution. We recall that it is given by V αβ = dδZ δα δ βπ (δZ ) − dδZ δαπ(δZ ) dδ Z δ βπ (δ Z ). The combination of elementary uncertainties is prior-independent when they are linearly combined (see also [4]). 3 An important conceptual remark is that the linear parametrisation in Eq. (2.5) requires the domain Finally we join the production and decay labels by adopting a unified notation Z = (X, Y ), so that The complete correlation matrix of the expected Higgs rates generated by theπ(δ Z ) distribution is then given by 8) and the complete covariance matrix is given by We emphasise that the primary purpose of this paper is to provide reliable covariance and correlation matrices associated with the expected Higgs rates. The choice of the exact shape forπ is left to the user -although a Gaussian shape is certainly well-motivated [4,5].

Covariance matrix from combination of elementary uncertainties
Here we derive a general expression for the covariance matrix of quantities that are affected by arbitrary elementary uncertainties. This simple formalism will be applied to the Higgs rates in the next sections. The reader only interested in the final results for the Higgs can safely jump to Section 5.
By definition the covariance relates two quantities at once. It is thus sufficient to consider only the case of two quantities to obtain a fully general result. Let us consider two quantities A, B, subject to elementary sources of uncertainties respectively labelled by n = (1, . . . , p) and n = (1, . . . , p ). The nuisance parameters are denoted δ A n , δ B n , so that . At this point, we make the crucial assumption that the magnitude of all relative uncertainties (i.e. the ∆'s) is small enough so that a Taylor expansion of A and B can be performed and stays valid up to δ = O(1). The most general form for the uncertainties is then for δ to satisfy δ > −1/∆ to forbid negative event numbers. For positive quantities, an exponential parametrisation σ = σ0e ∆δ is in principle much more natural. However, for the purpose of determining the covariance matrix, using the linear parametrisation is enough.
For the purpose of evaluating the covariance of A and B, it is enough to stop the ∆ expansion at linear order. In the most general case, one should assume correlations among all the sources of errors. The correlation matrix within the two groups of nuisance parameters δ A n , δ B n and the one between the two groups are respectively given by In particular, if a given source of uncertainty s affects both observables A and B, one has ρ ss AB = ±1. The correlation matrices ρ The next step is the combination of these elementary uncertainties. The global uncertainties on A and B are parametrised as Their correlation coefficient is denoted 5) and the correlation coefficient ρ AB is given by Equations (3.5) and (3.6) fully specify the covariance matrix of A and B, and thus determine entirely the errors on A and B induced by the elementary sources of error. The covariance matrix of (A, B) is expressed as Of particular interest is the correlation matrix between the elementary uncertainties of the two observables, ρ nn AB . From Eq. (3.4), it is clear that the more coefficients ρ nn AB are close to one, the more A and B get correlated. Similarly, the more coefficients ρ nn AB are close to zero, the more A and B get uncorrelated. This shows explicitly the competition between the uncertainties that correlate A and B and the ones that decorrelate A and B, which results in general in a non-trivial global correlation between A and B. Another observation is that the sources of uncertainty that are considered as elementary tend typically to be independent of each other. 4 Thus one can expect ρ nm A = 0, ρ nm B = 0 and ρ nn AB = 0 for n = n . The case n = n corresponds instead to the same nuisance parameter entering in both A and B, in which case the correlation is total, ρ nn AB = ±1, as already mentioned. In practice, the uncertainties that are not provided at the level of hadronic cross sections need to be propagated using a Taylor expansion with respect to ∆ 1. This will be the case for most of the errors in this study. For a cross section σ(Q) depending on a quantity Q subject to an uncertainty described as Q(δ) = Q 0 × (1 + δ Q ∆ Q ), the error gets propagated through the cross section as where σ 0 = σ(Q 0 ) and the error on the cross section corresponds therefore to When necessary, the partial derivative in Eq. (3.9) is obtained by varying the cross section with respect to the quantity Q.

Inventory of the uncertainties on Higgs rates
At the LHC, the Higgs boson can be produced on-shell and its decay products can be detected. The process of inclusive Higgs production followed by its decay is parametrised as pp where the ellipses denote extra states produced in association with the Higgs. The SM Higgs production mechanisms accessible at the LHC are i) gluon-gluon fusion (ggH), ii) vector boson fusion (VBF), iii) associated production with an electroweak gauge boson V = W, Z (VH), iv) associated production with a top quark pair tt (ttH), and v) associated production with a bottom quark pair bb (bbH). Moreover, the cross sections at different energies should be considered separately -although a high correlation between the errors of a given cross section taken at different energies can be expected. The production modes X will be therefore taken in the following list, tainty as correctly as possible, as the more one delves into the origin of uncertainty, the more its description becomes a set of elementary sources unrelated to each other. The fact that most of elementary uncertainties are independent is not a crucial feature for this paper, which is focussed on the combined covariance matrix. In contrast, independence is very important to claim that the shape of the combined prior converges to a Gaussian (see [5]).
The main decay modes are The index of all Higgs rates Z = (X, Y ) = (X 7 TeV , X 8 TeV , X 13 TeV , X 14 TeV , Y ) takes therefore its values in a list of dimension 34. Following Section 2, the Higgs production cross sections are denoted as σ X and the partial decay width as Γ Y . The nuisance parameters are written δ n X for a given cross section σ X , and are respectively associated with relative magnitude ∆ n X . Therefore one has The combined uncertainty is defined as the correlation matrix is defined as and the complete covariance matrix as Cov The notations are completely identical for branching ratios. For the latter, Eq. (2.4) has to be used to translate partial width uncertainties into branching ratio uncertainties.
Having defined the formalism, we now turn to the inventory of all the elementary sources of uncertainty. Our aim is to specify the relative magnitudes of all the elementary systematic uncertainties ∆ Z , as well as their possible correlations, ρ nm Z , ρ n m Z , ρ nn ZZ , with Z = (X, Y ). As discussed in the previous section, only the simultaneous variations of the rates are needed for the purpose of evaluating the correlation matrix among production modes. The details of these uncertainties have been discussed in Ref. [4] and references therein, apart from the uncertainty from the Higgs mass which is, to the best of our knowledge, taken into account for the first time in this paper. We give here a brief summary of the Higgs rate uncertainties.
• The experimental uncertainty on the parton distribution functions are released by the PDF4LHC15 working group in Ref. [6]. We use the PDF4LHC15 100 error set, recommended for precision physics, which includes a set of 100 independent elementary nuisance parameters. These errors propagate into the cross sections with a relative magnitude ∆ PDF,i [7], and on the Higgs mass m h are given in Tab. 1. For the Higgs mass we combined in quadrature the statistical and systematic uncertainties reported in Ref. [8].
• The uncertainties from Effective Field Theory (EFT) approximations in the ggH matrix elements are evaluated from Ref. [9]. For the EFT approximation applied  Finally, we need to specify the correlations among this set of elementary systematic uncertainties. Most of them are independent. The only non-trivial case is the one of the µ R and µ F scale dependence. For both µ R and µ F , some degree of correlation is expected among the production modes. To make sure we cover all the possibilities, we will consider two extreme cases, bearing in mind that the truth is somewhere in between. We consider 5 i) A fully independent case where there is an independent µ X R and µ X F for each production mode, and any µ X R is independent of any µ X F so that = 0 for any (X, X ) . ii) A fully correlated case where a universal µ R and µ F is assumed for all production modes, and µ R and µ F are further assumed to be 100% correlated so that there is a single nuisance parameter µ R = µ F ≡ µ. The correlation matrices are therefore 6    partial derivatives (see Eq.(3.9)), one obtains the ∆'s and ρ's that enter in the combination formulas Eqs. (3.5), (3.6). This fully determines the covariance matrix of the Higgs rates. To compute the cross sections we use SusHi v1.5.0 [11] for ggH and bbH, modified versions of VV2H v1.10 [12,13] and V2HV v1.10 [13,14] interfaced with LHAPDF [15] for VBF and VH, and an NLO version of HQQ [13,16] for ttH, based on Ref. [17]. For the ZH mode we also included the gg → ZH contribution, computed with VH@NNLO v1.2.1 [18]. The decay widths have been computed using HDECAY v6.51 [19]. The nominal values for σ 0 X and Γ 0 Y are given in Tabs. 2 and 3, respectively. These generators have been chosen in order to provide reproducible results. All hypotheses and cuts used to generate these cross sections are in principle publicly available.
The complete 34 × 34 correlation matrix ρ ZZ and the vector of relative magnitudes ∆ Z = (∆ X 7 TeV , ∆ X 8 TeV , ∆ X 13 TeV , ∆ X 14 TeV , ∆ Y ), all in both scale correlation cases i) and ii), are included in ancillary files attached to the arXiv version of this manuscript. The relative magnitudes ∆ Z are shown in Tab. 4. The 7,8,13,14 TeV blocks of the correlation matrices are shown in the Appendix, in Tabs. 5 to 12. The correlations of cross sections between various energies, which are encoded in the non-diagonal blocks of the complete ρ ZZ matrix given in the ancillary files, can deviate by ∼ 10% with respect to the correlations of cross sections at same energy. 7

On the correlations among ggH and ttH production modes
Our result on the correlation between the ggH and ttH modes depends strongly on the assumptions for scale correlations. The correlations we obtain are either very small or 7 The choice of generators is expected to affect only marginaly the correlation matrices. We checked that results given by, for example, HIGLU v4.34 [20,21] and SusHi give very similar uncertainty on the ggF cross-section. The discrepancy between the two estimations of the error is of order 0.1%. Moreover, the central values from HIGLU and SusHi are well compatible within these uncertainties. These cross-sections are in good agreement with NNLO values of the LHCXSWG.  positive. In case i) we obtain while for case ii) we get ρ ggH 7TeV ttH 7TeV = 82% , ρ ggH 8TeV ttH 8TeV = 82% , ρ ggH 13TeV ttH 13TeV = 80% , ρ ggH 14TeV ttH 14TeV = 80% .

(5.2)
These numbers arise from a non-trivial competition among the various sources of uncertainties. Indeed, the PDF uncertainties and the α s error tend to induce a negative correlation, while the m h and m t errors induce a positive correlation. Without the contribution from scale errors, the total correlations would be O(−10%), which is in reasonable agreement with the result of Ref. [6] (c.f. Table 3) [see also Ref. [22] (c.f. 2)). 9 The ρ ggH ttH case is a good example of a non-trivial combination of uncertainties, that requires the formalism presented in Sec. 3 in order to be treated correctly. 8 The correlation of ∼ −60% between ggH and ttH reported in Ref. [10] has been revised. 9 In addition to these extreme cases, one can also consider an intermediate one, close to i), where the scales µ X R and µ X F are 100% correlated for a given cross-section X, but are independent between each different cross-sections. This corresponds to the configuration ρ µ R µ F

Application to 7+8 TeV Higgs data
In this section we compute and display the simplified marginal likelihood, following closely the derivation of Ref. [5]. 10 Here we perform a Bayesian marginalisation, i.e. an integration over the nuisance parameters δ Z , with a multivariate normal prior whose covariance matrix is the ρ ZZ matrix obtained through the previous sections. This marginal likelihood is analytical and involves explicitly ρ ZZ and the relative magnitudes ∆ Z . One introduces the observed and expected signal strengthsμ I and µ I (θ), the latter being a function of a set of parameters of interest θ. The absolute statistical error on the signal strength µ I is written ∆µ I . The experimental efficiencies for a cross section X in a detection channel I are denoted I X . 11 = 0 for any (X, X ). The ggF-ttH correlations are then These numbers turn out to be close to those of Eq. (5.1). 10 The derivation of the simplified marginal likelihood consists of steps of error propagation and combination, and uses the fact that a central limit theorem applies to the combined errors. 11 The signal strengths are defined asμI =NI /N SM I , µI (θ) = N BSM I (θ)/N SM I , whereNI is the observed event number in channel I, N SM is the expected event number in the Standard Model, and similarly N BSM (θ)I is the expected event number in a modelisation of physics beyond the SM with parameters θ. In this section, it is assumed selection efficiencies such that, BSM I X = SM I X ≡ I X , following Refs. [23,24]. This relation is a good approximation when using Higgs signal strengths [24], is exact in the absence of higher-order derivative operators, and allows us to conveniently estimate the N BSM (θ)I event numbers.

The marginal likelihood is given bȳ
and where δ ZZ is the Kronecker delta, and ∆ I Y is the branching ratio uncertainty of the final state Y selected in the channel I. The L stat (θ) term is the Poisson likelihood containing the statistical error only. The L sys (θ) term encodes the effects of the systematic uncertainties. In particular, it contains all the correlations among signal strengths induced by the systematic uncertainties.
This computation is in fact slightly simpler than the general case presented in Ref. [5]. Although the number of elementary uncertainties is large, they get combined in only a few nuisance parameters δ X,Y at the level of cross sections and branching ratios, which have been checked in Ref. [4] to approximately follow a multivariate normal distribution, as expected from the central limit theorem. Because the number of nuisance parameters is smaller than the number of available Higgs channels, it is convenient to marginalise over these nuisance parameters instead of propagating the errors up to the event numbers.
As a final application, we present a global fit of the full set of 7 + 8 TeV data, following closely Ref. [4]. For ATLAS data, the diphoton final state results are taken from Ref. [25], the ZZ channel is from Ref.
[26], the W W channel from Ref.
[27], the bb from Ref. [28] and the ττ from Ref. [29]. The combined channels are studied in Ref. [31]. For CMS data, the diphoton final state has been presented in Ref. [32], the ZZ channel measurements are provided in Ref. [33], the W W ones in Ref. [34], the bb in Ref. [35] and the ττ in Ref. [36] (see also the combined channel analysis [38]). As a simple model of potential new physics effects, the deviations of the expected signal strengths from the SM are parametrised as where y t,b,c,τ are the SM Yukawa coupling constants (in mass eigenbasis), the subscript L/R indicates the fermion chirality, v is the Higgs vacuum expectation value, g hW W = 2M 2 W /v and g hZZ = M 2 Z /v are the electroweak gauge boson couplings. The c V,f parameters are defined such that the limiting case c V,f → 1 corresponds to the SM.
The Bayesian regions at 68, 95, 99% credibility level are presented in the c V −c f plane in Fig. 1, assuming flat prior for c V,f . The shapes of the regions are rather similar between the two extreme cases for renormalisation/factorisation scale correlations. The main difference between the two extreme cases turns out to be the shift of the best-fit regions with respect to the fit without Higgs rate uncertainties (shown with dotted regions in Fig. 1). This discrepancy between the two extreme cases of scale correlations illustrates the importance of properly including the correlation matrices in the analyses of the Higgs rates. The impact of correlations can also be seen by comparing Fig. 1 with the best-fit regions of Fig. 6 from Ref. [4], where schematic correlations have been used.
As the LHC will accumulate more and more data on the Higgs processes, the statistical error bars will decrease, letting the uncertainties on the Higgs rates be among the dominant ones. Therefore, the correct treatment of these uncertainties is expected to become more and more crucial for the global Higgs fits.

Conclusion
We have evaluated the full covariance matrix of the expected Higgs production and decay rates, through a combination of their elementary sources of uncertainty. Our calculation follows from a direct application of probability theory, and is completely consistent with the framework of Bayesian statistics. The obtained covariance matrix is prior-independent, and is also consistent with a frequentist combination of Gaussian elementary uncertainties. A frequentist combination of non-Gaussian uncertainties is prior-dependent and beyond the scope of this work.
We provide the error magnitudes and correlation matrices on the set of Higgs cross sections and branching ratios at √ s = 7, 8, 13 and 14 TeV. Most of the elementary uncertainties are evaluated via error propagation in the SusHi, VV2H, V2HV, HQQ and HDECAY codes. This error propagation is equivalent to truncating the expansion in the relative magnitude of the errors (∆) at first order, which we systematically checked to be a good approximation. Our set of elementary uncertainties includes the error on the Higgs mass determination by CMS and ATLAS.
The error magnitudes and the correlation matrices are provided under various formats with the arXiv version of this paper, and can be readily used for further global fits of the Higgs data. The lack of knowledge about the correlations among the elementary uncertainties coming from the various renormalisation and factorisation scales has been taken into account by considering two extreme cases for these correlations. In these extreme cases, we find for example either a very small or a large positive correlation between the ggH and ttH rates.
Finally we present global fits of the latest 7 + 8 TeV dataset, showing the impact of the full Higgs rate covariance matrix on the Higgs coupling determination. In doing so we use a simplified likelihood framework [5], relying on the fact that the distribution of the combined nuisance parameters is approximately Gaussian.

A Blocks of the Higgs rate correlation matrix
Here we display the 7, 8, 13, and 14 TeV blocks of the correlation matrix ρ ZZ , assuming either scale correlations i) or ii), as defined in Eqs. (9) and (4.9). ggH VBF ZH WH ttH bbH W W ZZ γγ Zγ gg bb cc ss ττ µμ                 Table 12: Correlation matrix between the expected Higgs production and decay rates at √ s = 14 TeV, expressed in percent, assuming the scale correlations ii).