{'fraction_non_alphanumeric': 0.032412110019128026, 'fraction_numerical': 0.009735505573510982, 'mean_word_length': 4.564151497357604, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decoupling Theoretical Uncertainties from Measurements of the Higgs Boson

We develop a technique to present Higgs coupling measurements, which decouple the poorly defined theoretical uncertainties associated to inclusive and exclusive cross section predictions. The technique simplifies the combination of multiple measurements and can be used in a more general setting. We illustrate the approach with toy LHC Higgs coupling measurements and a collection of new physics models.


Introduction
The discovery of the Higgs boson [1][2][3][4] has initiated a vigorous program of precision Higgs measurements at the LHC. The primary focus is on the couplings of this new state in terms of an effective Lagrangian based on the Standard Model gauge structure. The key question is if this Lagrangian is fully renormalizable, as predicted in the Standard Model. Deviations in the properties of the newly found particle from the Standard Model predictions would indicate the first hint for physics beyond the Standard Model.
These precision studies provide one of the key motivation for new accelerators [5,6]. Complicating the interpretation of both the current results and future prospects of the LHC is the presence of uncertainties on the theoretical link between Lagrangian parameters and LHC observables. In this paper we outline a strategy to decouple the theoretical uncertainties from the experimental results while retaining the ability to incorporate those uncertainties in a subsequent stage we refer to as recoupling.
In reporting results, experimentalists strive to find a representation of the results that is as free of theoretical assumptions as possible while still being convenient for addressing specific theories of interest. For example, the measurement of a cross-section in a well defined fiducial region allows experiments to report results that are not tied to a specific theory at the cost of requiring the reader to estimate a model-dependent acceptance. This approach has advantages, but is difficult to generalize. For example, extracting the Higgs couplings requires simultaneous inference of multiple production and decay modes from the combination of multiple searches with correlated experimental systematic uncertainties [7][8][9][10]. A convenient approach to combining multiple measurements is the Best Linear Unbiased Estimator (BLUE) technique [11]; however, that approach is limited to measurements that are well within the Gaussian regime. The methodology presented here addresses both of these challenges.
The first wave of Higgs coupling measurements [2,13,14] involves signal strength modifiers µ pd that scale the total rate of events for a given combination of production (p) and decay (d) modes, The ATLAS experiment has made such likelihood functions available [15]. If we focus on measuring the Higgs couplings in the Standard Model Lagrangian we can link the event rates to the set of shifted Higgs couplings g x /g SM x to any Standard Model particle x. Kinematic distributions (and thus selection efficiencies and detector acceptances) will not change as long as the Higgs couplings are roughly in their Standard Model ranges. While σ p and Γ d can usually be directly linked to a specific coupling g x , the appearance of the total width forces us to make non-trivial assumptions and induces strong correlations. The LHC Cross Section Working Group has defined a number of benchmark scenarios that specify how the production cross section and branching ratios are modified with respect to the standard model values [16], and the ATLAS and CMS collaborations are reporting their results in terms of these benchmark scenarios.
Thus far, results for signal strengths and Higgs couplings include the theoretical uncertainty as part of the Standard Model prediction. Large uncertainties on the production appear because of unknown higher orders in the perturbative QCD expansion. The problem with these theoretical uncertainties is two-fold: first, the size of the associated uncertainty in not known. Second, the uncertainties are not statistical in nature -there is no random variable associated to missing higher orders. While there has been effort to assign a Bayesian degree of belief to these uncertainties [17] and an effort to complete the perturbative series [18], we lack an objective probabilistic interpretation of these uncertainties. While the field may be able to settle on a consensus for both the size and shape of these uncertainties, it is an area of debate and will certainly change with time as theoretical progress is made. Therefore, our view is that we should develop a technique in which the theoretical uncertainties are factorized from the experimental result.
In Section 2 we describe the approach allowing us to decouple and recouple the theoretical uncertainties from the experimental measurement. This approach is independent of our specific application to the Higgs coupling measurements at the LHC. In Section 3 we briefly discuss the issues with theoretical uncertainties at hadron colliders and how they affect the Higgs couplings measurement. Section 4 demonstrates the procedure with a toy example based on the ATLAS results presented in Ref. [14]. Next, we consider some specific new physics scenarios that can be tested through Higgs coupling measurements in Section 5. A simple example is worked in detail in Appendix A and more details on the new physics models in Appendix B.

The Statistical Model
In this Section we outline briefly the statistical modeling approach used by the LHC experiments following Ref. [19]. Once the statistical model has been constructed, the LHC experiments employ the profile likelihood ratio to define confidence intervals [20] on the parameters of interest.
The coupling measurements require defining several disjoint categories of events, indexed by c, which satisfy specific selection criteria designed, in part, to be particularly sensitive to a particular production or decay mode. Each category has associated to it an expected (observed) number of events ν c (n c ). Each category may also have some discriminating variable(s) x, such as an invariant mass, and a corresponding probability density function f (x). The data associated to the c th category is denoted D c = {x 1 , . . . , x nc }. In general the expected number of events and their distribution will depend on both the signal strength parameters µ and nuisance parameters α. The nuisance parameters α parametrize both theoretical and experimental uncertainties. Typically, the expected number of events is written ν c (µ, α) = p,d µ pd s cpd (α) + b c (α) , (2.1) where b c (α) is the background in this category, and s cpd (α) = L(α) σ SM p (α) BR SM d (α) cpd (α) (2.2) is the expected Standard Model signal for production mode p and decay mode d predicted by the product of the integrated luminosity, cross section, branching ratio, and selection efficiency, each of which may depend on theoretical uncertainties parametrized by α. Constraint terms associated with systematic uncertainties are described as f i (a i |α i ), where α i are nuisance parameters and a i are auxiliary or control measurements designed to estimate those nuisance parameters. In the case of experimental uncertainties, there are often real auxiliary measurements that are summarized by f i (a i |α i ). However, in the case of most theoretical uncertainties, the auxiliary measurement does not truly exist and an ad hoc f i (a i |α i ) is introduced for convenience. The full likelihood function used by the experiments [21] can be written as a product of the main experimental measurement and the constraint terms L full (µ, α) = c∈category Pois(n c |ν c (µ, α)) nc e=1 f c (x e |µ, α) ≡L main (µ,α) i∈syst f i (a i |α i ) ≡Lconstr (α) . (2.3) Typically, confidence intervals are then defined by contours of the profile likelihood ratio whereα(µ) is the conditional maximum likelihood estimate andμ,α are the unconditional maximum likelihood estimates [20].

The Effective Signal Strength
We are interested in inferring the values of the signal strength parameters µ, which scale the signal expectation s cpd ; however, the presence of experimental and theoretical uncertainties mean that the signal and background expectations are functions of the nuisance parameters as in Eq.(2.1). Alternatively, we can introduce effective scale factor with respect to some fixed reference scenario α 0 , so the expected number of events can be re-written The key conceptual jump is to realize that we can think of µ eff cpd not as a function, but as a well-defined parameter free of theoretical uncertainty that we can infer directly.
While the signal strength parameters µ pd we ultimately want to infer are independent of the details the individual analysis categories, the effective signal strength µ eff cpd is specific to the c th category due to the selection efficiency (and, more generally, the distributions f c (x|µ, α)). In particular, selection requirements that leverage exclusive or differential properties of a specific production and decay will introduce category-specific theoretical uncertainties.
State-of-the art Higgs property measurements can include hundreds of categories of events. The full likelihood L full defined in Eq.(2.3) encodes a detailed description of the correlated effect of common experimental systematic uncertainties. In practice, we want to find some coarse graining of the many categories into a few groups so that we can ignore the category index within each group. By taking a common category-weighted signal strength factor we can suppress the category index and make the replacement µ eff cpd → µ eff pd . By doing so, we can define The goal is to show that by providing L eff (µ eff ), the reparametrization µ eff (µ, α), and the constraint terms f i (a i |α i ) that we can recouple these ingredients and approximate the full likelihood In the case of an inclusive cross section uncertainty µ eff cpd is the same for all c -in which case we say the effect of that uncertainty is category-universal. If all uncertainties are category-universal, then it is possible for this approach to be exact. More generally the grouping of categories will lead to µ eff (µ, α) encoding some weighted effect from the individual categories. We will discuss some examples in the next section.

Reparametrization Templates
The art of this approach lies in choosing a template for the reparametrization in which the coefficients of the template can be effectively deduced from the likelihood. We treat the likelihood as a "black box" since the diversity and complexity of statistical models created by experimentalists and encompassed by Eq.(2.3) is so diverse. In the case that the reparametrization µ eff (µ, α) is category-universal, which can be trivially achieved if the ingredients are explicitly provided for each category, this reformulation of the likelihood can be exact.
For example, a natural way to parametrize the dependence of the expected signal due to uncertainties that modify inclusive production cross sections is In this situation, µ eff pd (µ, α) is not linear in the full set of parameters, but is bi-linear in (µ pd , α). This µ scaling is important for capturing the behavior of the likelihood away from the maximum likelihood estimate and distinguishes this approach from techniques such as BLUE [11].
In the case of uncertainties that only effect the background through the equivalent form of the effective signal strength is 11) which is linear in µ pd and α. Because b c (α 0 )/s cpd (α 0 ) is a constant, this pre-factor can be absorbed into φ ci and the category-weighted effect would simply be written φ i .
Another example is motivated by the large uncertainty associated to gluon-fusion (ggF) Higgs production with two additional jets, which populates the categories meant to isolate weak boson fusion (VBF). We would expect the uncertainty to modify the µ eff c,p=VBF,d signal strength for weak boson fusion, but be proportional to µ ggF . Thus we should anticipate templates of the form Combining these three situations, a fairly general template would be where we can identify η pi = η p pi from Eq.(2.8) and b c (α 0 )φ ci /s cpd (α 0 ) = φ i . This general template involves n α (n 2 p + 1) coefficients (per grouping of categories).
The LHC experiments must cope with additional complications. First is the fact that uncertainties on α i can be large, and linear extrapolation of the effect of this uncertainty on a signal or background rate can lead to unphysical expectations like s cpd < 0 or b c < 0. To cope with this, the LHC experiments typically implement log-normal priors (constraint terms), which are implemented via reparametrization so that f i (a i |α i ) is a Gaussian distribution and s cpd (α i ) is an exponential response function. This approach has two advantages: it ensures s cpd (α i ) > 0 and it allows for multiple signal expectations with different sensitivities to a common source of uncertainty to have a similar log-normal behavior parametrized by α i .
Second, experiments often have a few numbers with which to parametrize the signal (and background) expectations. Typically, this is based on a nominal s cpd (α i = 0) and "±1σ" variations on s cpd (α i = ±1) (using a conventional scaling of α i ). Often the variation of s cpd (α i = ±1) with respect to s cpd (α i = 0) is asymmetric, which requires some assumptions about the intermediate behavior and the use of ad hoc interpolation algorithms, such as second degree polynomial or higher degree polynomials that match the exponential extrapolation up to the second derivative [22,23].
Finally, to ensure the positivity of s cpd (α) under the joint effect of several sources of uncertainty, it is common that the interpolation/extrapolation template is multiplicative over the nuisance parameters. In generic terms, this often leads to signal expectations parametrized as where I(α i ) is some interpolation/extrapolation function based on s cpd evaluated at several points in α i . These ad hoc choices influence the resulting inference and further strengthen the motivation to decouple theoretical uncertainties from the presentation of experimental results.

Determining the Coefficients of the Reparametrization Template
Following Eq.(2.7), our goal is to show that by providing L eff (µ eff ), the reparametrization µ eff (µ, α), and the constraint terms f i (a i , α i ) we can "recouple" the ingredients to approximate the original likelihood L full (µ, α). Below we develop two approaches to determine the coefficients of these reparametrization templates.
Via the local covariance matrix A direct path towards this goal is to use a reparametrization template like the ones introduced in Sec. 2.3 and determine the coefficients that reproduce the local covariance structure around the maximum likelihood estimate (μ,α). The maximum likelihood estimate is defined by −∇ ln L = 0, while the Hessian of − ln L captures the local covariance structure and is referred to as the observed Fisher information matrix or precision matrix. In particular, the observed Fisher information is an estimate of the inverse of the covariance matrix It is helpful to factorize the full likelihood in Eq.(2.3) as a product of the main measurement and the constraint terms for the nuisance parameters α. This is equivalent to decomposing the full information matrix into two parts By fixing α = α 0 the experiments can provide the (profile) likelihood for the effective signal strength L eff (µ eff ) via Eq. (2.6). We denote the Fisher information for this effective signal strength as V −1 eff . Reparametrizing L eff (µ eff ) via µ eff → µ eff (µ, α) implies the following transformation to the Fisher information matrix 17) where the Jacobian of the transformation is defined by Even though V −1 eff may be a constant matrix, the Jacobian in Eq.(2.17) depends on (µ, α) so that the resulting information matrix need not be constant. This is an important point as the analysis below is only sensitive to the linear approximation of µ eff (µ, α) at (μ,α) although the template may encode important non-linear behavior. Combining Eq. (2.16) and Eq.(2.17) gives us The final stage of this procedure is either to check if Eq.(2.19) holds for a particular reparametrization µ eff → µ eff (µ, α) or to determine the coefficients for a template by imposing this equality.
Perhaps the most intuitive way to think of the effect of a given uncertainty is the shift it induces in the best fit value of the signal strengths µ. This way of thinking requires keeping the α fixed and considering the likelihood is only a function of µ. For example, Fig. 1(b) shows the shift in the L(µ) contour due to a shift in the gluon fusion inclusive cross section. We denote the best fit µ with fixed α aŝ µ fix (α) ≡ argmax µ L full (µ, α) (2.20) withμ fix (α) =μ. Because of the definition L full (µ, α 0 ) ≡ L eff (µ) in Eq. (2.6), this implies thatμ fix (α 0 ) =μ eff . Finally, if the main measurement is not able to measure the nuisance parameters, i.e. there is flat direction in the likelihood or a degeneracy between µ and α, then α 0 =α. By using a specific template -in this case Eq.(2.9) -we can equate µ eff p (μ fix p , α i ) =μ fix p (1 + η ip α i ) =μ p and then explicitly evaluate the partial derivative that quantifies the shift to the best fit value of µ These partial derivatives are visualized as vectors in the signal strength plane in Fig. 2.
Another way to arrive at Eq. (2.21) is to approximate the likelihood in the neighborhood of the maximum likelihood estimate as a multivariate Gaussian G(µ, α|μ,α, Σ). The conditional distribution of µ given α is also a multivariate Gaussian with mean given bŷ 22) where Σ c is the upper-right sub-block of the full covariance matrix Σ (2. 23) In situations that the main measurement does not constrain or pull on the nuisance parameters, Σ α is just the covariance matrix associated to the constraint term defined in Eq. (2.3). In general Σ c will depend both on the constraint terms and the main measurement. Clearly, the conditional likelihood with α fixed as in Eq. (2.20) is independent of the constraint term, thus the product (Σ c Σ −1 α ) can only depend on the details of the main measurement. Through the techniques developed for regression in the general linear model [24] one can show that ∂μ fix where the right-most equivalence is specific to the template of Eq.(2.9) and explicit computation requires an assumption about the constraint terms even though the result is independent of those assumptions.This is shown explicitly for a simple example in Appendix A. Note, one can eliminate the constraint terms and only consider the main measurement, in which case Σ in Eq.(2.24) is a singular matrix and one must use the pseudo-inverse. While this approach is mathematically cumbersome, it is equivalent to Eq.(2.21) and Eq. (2.25).
A third approach is to impose Eq.(2.17) with V −1 constr subtracted from both sides (2.25) This allows us to proceed without assumptions on the constraint terms, which we are trying to decouple from the procedure. The upper-right sub-bloc of this matrix leads to a system of equations that can be used to determine the coefficients of the template. Again in the case of template of Eq.(2.9), these linear equations provide the same solutions for η ip as Eq. (2.21) and Eq. (2.24). These equations are solved explicitly for a simple example in Appendix A.
It is worth noting that in the case of an arbitrary template µ eff (µ, α; η), these equations fix the linear behavior at α 0 (i.e. ∂µ eff /∂α i | α=α 0 ), which may be non-trivially related to the coefficients η. For example, this is the case for the most common parametrizations encompassed by Eq. (2.14) and used by the LHC experiments and to implement asymmetric uncertainties [22,23].
The information matrix from the full likelihood function has (n 2 + n)/2 independent components, where n = n p + n α is the sum of the number of parameters of interest and nuisance parameters. An n p × n p sub-block describes Σ µ or V −1 eff and a n α × n α sub-block describes Σ α or V −1 constr. . Thus, there are a remaining n p ×n α numbers in the local covariance encoded in Σ c or V −1 main that can be used to determine the coefficients of the template. This counting matches precisely for the templates of Eq.(2.8) and Eq.(2.3), but leaves ambiguity for the more general template of Eq. (2.12).
The art of choosing the template As mentioned above, the local information provided by the observed information matrix evaluated at the maximum likelihood estimate is only sensitive to the linear behavior of µ eff (µ, α) at (μ,α). For example Eq.(2.8) and Eq. (2.3) can have the same linear behavior at (μ,α) though the effect of an uncertainty α i scales with µ in the former and not in the latter. This difference is illustrated in Fig. 3. Thus, while both templates will be equivalent locally, if one wants to capture the higher-order corrections encoded in the different templates some additional information will be required. This information can either be injected by hand or by using information about the likelihood away from (μ,α), as described below. For instance, one may know that a particular nuisance parameter α i is associated to the uncertainty in the inclusive cross section for the p th production mode, thus the coefficients of the general template can be restricted by hand to η p pi = δ pp η ip and φ ip = 0. In another extreme case, the uncertainty on gluon fusion production with two jets primarily affects the inference of µ eff p=VBF signal strength for vector boson fusion, but the size of the shift is proportional to µ ggF . Here the general template would be constrained so that η p pi = 0 unless p = ggF and p = VBF or vice versa. We explore a simple model for this situation in Scenario C of Appendix A.
Via a global learning approach Ideally we would have a formalism that would work with the black box likelihood L full (µ, α) and a general template like Eq.(2.12) without having to introduce by hand restrictions on that template as described above. In order to do that we must introduce information about the likelihood away from (μ,α). A flexible approach to that problem is based on the ideas of machine learning and function approximation in which one aims to minimize a loss function with respect to some model parameters (in this case the template coefficients). The loss function needs to be a scalar evaluated over the (µ, α) space, which in frequentist terms has no measure. However, from the point of view of decision theory, one can introduce some weighting over the parameter space (without regarding it as a Bayesian prior) and evaluate The choice of the weighting function π(µ, α) is arbitrary, but a reasonable choice is the Bayesian posterior with respect to some baseline constraint terms interpreted as a prior on α, which leads to This approach will put the highest weight for L recouple to approximate L full near the best fit point (μ,α) and lesser weight as one moves away from it. Importantly, this minimum loss solution can be found numerically and is well defined even when the number of parameters in the template is larger than n p × n α . Furthermore, this approach may be more robust in the case of very complicated likelihood functions where the numerical accuracy of the covariance matrix, information matrix, and partial derivatives needed in Eqs. (2.21), (2.24), and (2.25) may be poor. In situations where additional experimental uncertainties α exp have been profiled in providing L eff (µ eff ) -a situation discussed in more detail in Sec. 2.5 -one must take care that the loss function makes the comparison for equivalent values of the profiled nuisance parameters. For example, when creating L eff (µ eff ), one can keep track of the profiled valueŝ α exp (µ eff ) and then in Eq.(2.26) make the replacement L full (µ, α) → L full (µ, α,α exp (µ eff (µ, α)) ). (2.28) We demonstrate the effectiveness of this approach in Sec. 4 and Scenarios B and C of Appendix A.
Software The software implementation of the reparametrization templates described in Sec. 2.3 as well as the three strategies for determining the coefficients of those templates from the local covariance matrix and the learning approach described in described in Sec. 2.4 is available at Ref. [25]. Experiments can use this software on their full RooFit/RooStats [26] models and obtain the effective likelihood L eff (µ eff ) as well as the reparametrization µ eff (µ, α) for publication. These ingredients can be supplied in a technology independent format enabling others to perform the recoupling stage, modify constraint terms associated to theoretical uncertainties, combine multiple results, and create likelihood scans in benchmark models.

Grouping of Categories and Combinations
In Section 2.2 we discuss the coarse graining of the many categories into a few groups. By taking a common category-weighted signal strength factor we are able to suppress the category index and make the replacement µ eff cpd → µ eff pd for each group. In practice, this grouping is often based on the Higgs decay mode indexed by d. For example, ATLAS provided the profile likelihood λ(µ ggF , µ VBF ) for the three decay modes d = γγ, W W, ZZ [15]. These profile likelihoods are colloquially referred to as "likelihoods"; however, they have eliminated several nuisance parameters via profiling as defined in Eq. (2.4). In addition to the common theoretical uncertainties, these three likelihoods share common experimental systematic uncertainties. Naïvely combining these likelihoods will thus double count the common constraint terms f i (a i |α i ) and lead to an artificial reduction in the uncertainty. This effect can be seen by comparing a naïve combination of the three individual profile likelihoods provided by ATLAS with the official ATLAS combined result (which avoids the erroneous double counting). The same effect is demonstrated with a toy example in Fig. 4.
The decouple/recouple procedure described here has focused on theoretical uncertainties for reasons that will be elaborated in Sec. 3; however, the same technique can be used for experimental uncertainties. In particular, sources of systematic uncertainties that are anticipated to be common to other analysis should not be profiled or it will not be possible to avoid double counting of constraint terms. The luminosity uncertainty is an example of an experimental uncertainty that should not be profiled, but instead the effect of the uncertainty should be parametrized in the template µ eff (µ, α). As a rule of thumb, one can safely profile uncertainties that are specific to the categories in a given analysis, but should include common (correlated) sources of uncertainty in the reparametrization.

Anticipating Higher Dimension Operators and Other Effects
At this point, the Higgs coupling analysis is primarily focused on small deviations in the coefficients of the SM operators. This justifies the scaling of the SM Higgs expectations as µ pd s cpd , since the efficiency and acceptance of the modified signal in the c th category is not effected. The presence of non-SM operators will generically change the pattern of the signal across the categories in a way specific to the operators under consideration [27]. In particular, operators that modify kinematic distributions will affect the cut efficiencies and acceptances cpd (α) of Eq.(2.2) One approach to anticipate these studies is to characterize the effect of a deviation in the number of events in each individual category. In particular, by considering a perturbation of the form where α c is some addition or reduction to the number of events in that specific category. This would add one nuisance parameter for each category, each of which would have the effect of shiftingμ eff . These α c can be seen as a basis for possible new physics effects and can be incorporated into µ eff (µ, α) using the same formalism. To use this information in the context of a specific new physics model, one would need to estimate the perturbation to the number of events in each category (α c ) through knowledge of the details of the selection for that category. If that is possible, then one can parametrize each α c in terms of the coefficients of these new operators. This may involve parametrizing several α c in terms of just a few operator coefficients.

Theoretical Uncertainties
In this section we focus on the dominant theoretical uncertainties plaguing Higgs coupling measurements at the LHC. There are several sources of uncertainties, which can be associated with the theoretical description of LHC processes [4]: 1. the inclusive rate for a given production and decay, for example the total cross section σ tot (gg → H + X).
2. efficiency and acceptance of kinematic cuts for the signal process, for example cuts on the transverse momentum of the Higgs boson or a central jet veto.
3. the parton densities describing the initial state, which have both experimental and theoretical uncertainties.
4. the parton shower, hadronization, underlying event, and other key parts of the Monte Carlo description of signal and background events.
Linking event rates in a given kinematic region to perturbative theoretical predictions might be considered the critical problem for the LHC's precision Higgs program. According to the above list it can be separated into two parts: first, the total Higgs production cross section parametrically depends on the Higgs coupling in a given renormalization scheme. This relation is based, for example, on collinear factorization, includes a small number of scales, and can be computed in perturbative QCD. However, this connection cannot be replaced or checked experimentally, which means that Higgs coupling measurements will eventually be limited by theory predictions [5,6]. Kinematic cuts linked to Higgs-plus-jet production [28], tagging jets in weak boson fusion [29,30], or advanced analysis methods lead to additional complications. They do not automatically obey collinear factorization and induce a large number of energy scales and possibly large logarithms [31][32][33][34][35]. Already for the relatively safe inclusive Higgs production cross section, the next-next-to-leading order corrections more than double the predicted number of produced Higgs bosons [36][37][38][39]. The source of these problems is the choice of protons as colliding particles, combined with the large value of the strong coupling constant, which leads to a poor convergence of the perturbative QCD description.
Aside from the poor convergence in perturbative QCD the corresponding theoretical uncertainties face a more fundamental problem: there does not exist any well-defined estimate for an uncertainty on a production cross section computed in perturbative QCD. Traditionally, we derive a range of allowed cross section values using a variation of the factorization and renormalization scales. This is based on the fact that these scales are artifacts of the perturbative expansion, so the scale dependence has to vanish once we include all orders in perturbation theory. This recipe captures some of the effects of the theoretical uncertainties, but we know from the Drell-Yan process or Higgs production in gluon fusion that it does not give a conservative estimate. For Higgs production at the LHC we know that equal variation of the factorization and renormalization scales leads to a cancellation of the two scale dependencies as one possible reason for the underestimate of the theoretical uncertainties [37]. On the other hand, a separate variation of the factorization and renormalization scales ruins the theoretical description of the scale-dependent parton densities in terms of resummed large logarithms of transverse jet momenta [4]. Independently varying the two scales makes sense in a phenomenologically, but not in terms of a formal QCD description. A recent, promising approach might be to combine the scale variation with an extrapolation of the perturbative series [18]. However, no matter what recipe we choose, it is clear that the size of any theoretical uncertainty is poorly defined.
In addition to the overall size of the uncertainty, most statistical analysis techniques require the uncertainty to be quantified in the likelihood, i.e. the constraint terms of Eq. (2.3). Frequentist interpretation of the likelihood requires one to be able to identify some random observable who's probability is defined with a corresponding ensemble. No traditional ensemble defined by repeated observation exists for theoretical uncertainties. Missing higher order terms in a perturbative prediction are not random in nature, there is no ensemble, thus L constr (α) is ill-defined. Attempts to provide a meaningful degree of belief for the perturbative uncertainties [17] requires a Bayesian interpretation and attributes a fundamental meaning to the perturbative QCD expansion.
In a practical sense we must make some choice for the likelihood L constr even if it fundamentally ill-defined. Unfortunately, the various choices have their own unique pathologies. The often-used Gaussian and log-normal distributions inject information by ascribing a preferred (peak) value and allowing for long tails into what are untenable values for the theoretical prediction. Alternatively, we can try to avoid introducing a peak by adopting a constant likelihood function, as proposed in the RFit scheme [10,40]. The choice of distributions might not appear numerically relevant, but in combination with a profile likelihood, it can lead to significant differences. For instance, the combined effect of uncertainties described in the RFit scheme add linearly rather than in quadrature [4].
Given that the size as well as the shape of the theoretical uncertainties entering the Higgs couplings measurements are subject to variations in time, geography, and personal taste, the best solution is to decouple them from the experimental result as proposed in Section 2. This allows for efficient tests of different assumptions on the theoretical uncertainties as well as an efficient implementation of a perceived or actual improvement.

A Toy Example
We now consider a toy example that is representative of the current ATLAS results for H → γγ, W W, ZZ [14]. The statistical model here is based purely on number of events in various categories without including discriminating variable distributions for m γγ , m 4 , or m T (the terms f c (x|µ, α) in Eq. (2.3)). Each decay mode groups together several categories of events that together provide sensitivity to the underlying production modes. We model the uncertainty on the signal expectation from luminosity, parton distribution functions, the inclusive gluon fusion cross section, and the uncertainty on the cross section for gluon fusion in association with two or more jets.
The H → γγ likelihood includes a simplified version of the 14 categories considered by ATLAS including the low-and high-p Tt categories, the low-and high-mass 2-jet categories, the high-E / T significance category, and the lepton-tagged category. The H → ZZ → 4 likelihood includes three ggF-like categories (for 4µ, 2e2µ, and 4e) as well as a VBF-like and a VH-like category. The H → W W → ν ν likelihood includes 0-, 1-, and 2-jet categories. The HistFactory script and RooFit/RooStats workspace for this toy model can be found at Ref. [41]. Standard Model 68% CL full model 68% CL w/o theory uncert. Figure 1(a) shows the likelihood contours for the three different decays with and without theoretical uncertainties, which are modeled using Gaussian constraint terms and a linear response as in Eq. (2.8). Fig. 1(b) shows the shift to the contours without theory uncertainty due to fixing the inclusive gluon fusion cross section to its "+1σ" value estimated from QCD scale variation. The larger gluon fusion cross section leads to a smaller inferred value for µ ggF . This can be repeated for each of the nuisance parameters α i as in Eq. (2.21). The corresponding partial derivatives are visualized in Fig. 2  Next we demonstrate the recoupling stage based on the local covariance structure at (μ,α). If we choose the reparametrization template in Eq.(2.9) we can solve for the η pi coefficients that reproduce the local covariance structure using either Eq.(2.21), Eq.(2.24), or Eq. (2.25). Having determined the η pi coefficients, the reparametrization µ eff (µ, α) is specified and we can create the decoupled likelihood via Eq.(2.7). Fig. 3(a) compares the full likelihood to the recoupled likelihood using this template. The result is not bad, the effect of the theoretical uncertainties have been recovered; however, there is a significant discrepancy for H → ZZ → 4 at the top of the contour.
The source of this discrepancy between L full and L recouple contours is the non-category universal nature of the theoretical uncertainty associated to gluon-fusion+≥2 jets, which primarily affects the VBF-like categories. Correspondingly, the effect of this uncertainty is a shift in the µ VBF direction, as illustraited in Fig. 2. The reparametrization template in Eq.(2.9) implies that this uncertainty would scale with µ VBF , though in the full model this source of uncertainty scales with gluon fusion. At the top of the contour, µ VBF ≈ 5, while µ ggF ≈ 1, which leads to a large inflation on the effect of this uncertainty.
As discussed in Sec. 2.4 the local covariance structure of L recouple is only sensitive to the linear behavior of the template, so this physical insight must either be put in "by hand" or one should instead use more global information about the likelihood as in the learning approach. To illustrate the "by hand" approach, we can consider the more general template in Eq. (2.13). The general template can be constrained so that η p pi = 0 unless p = ggF and p = VBF (or vice versa), where i is the index for the nuisance parameter associated to the gluon fusion + ≥ 2 jet cross section. In addition, we set φ i = 0 for all i. The same procedures are followed to determine the coefficients that reproduce the local covariance structure for this template. The resulting contour is shown in Fig. 3(b), where the agreement is improved, particularly near the top of the H → ZZ → 4 contour. The alternative learning approach works with the general template in Eq.(2.13) without restricting the terms by hand. The optimized values of the coefficients leads to even better agreement of the contours compared to those shown in Fig. 3(b).   Figure 3. Comparison of effective model using different templates with same local covariance Next, we combine these three likelihoods. As discussed in Sec. 2.5, the experiments typically present results grouped by decay mode. In this example, the three likelihoods share common sources of systematic uncertainty. Naïvely combining these likelihoods double counts the common constraint terms f i (a i |α i ) and leads to an artificial reduction in the uncertainty. This effect can be seen in the top plots of Fig. 4, which compare a naïve combination to the the full combined result that avoids the erroneous double counting. There is no meaning to a combined contour in the signal strength plane µ ggF,d − µ VBF,d due to the decay index d. Thus, two simple 2-parameter benchmark models are used to present the result. Fig. 4(a) shows a benchmark model that dictates all µ p,d based on the scaling of the fermionic couplings g f = κ f g SM f and weak boson couplings g V = κ V g SM V . Similarly, Fig. 4(b) considered a simple 2-parameter benchmark in which scale the effective Hgg by a factor κ g and the effective Hγγ coupling by a factor κ γ [16]. The contours from the naive combination are considerably smaller in the κ V and κ g directions, leading to poor agreement with the full combination (based on all 22 categories of events without double counting constraint terms).
In contrast, the recoupling approach allows one to avoid double counting constraint terms and for a consistent profiling over the common sources of uncertainty (both theoretical and experimental), which leads to an improved agreement. In addition, the learning approach of Eq.(2.26) has been used for L recouple in Fig. 4. Note, in this example the sources of experimental uncertainty that are unique to one of the decays were profiled in L eff (µ eff ), as would be done by the experiments. Thus, there is a small effect that is neglected due to the fact that the profiling of these analysis-specific uncertainties is slightly affected by identifying the common α in the combination.
Perhaps the greatest benefit of this approach is that the theoretical uncertainties have been decoupled from the experimental result encoded in L eff (µ eff ). Moreover, the recoupling procedure in Eq.(2.7) allows one to replace the constraint terms in L constr (α) with some  other constraint terms L constr (α). The bottom plots of Fig. 4 demonstrate the change by using the same L eff and reparametrization template, but with uncertainties inflated by 30%. This recouped approach is compared to the full model where the same modification is made to the constraint term. In addition, an example of replacing Gaussian constraint terms with the Rfit scheme is given in Appendix A.

New Physics Effects vs Theoretical Uncertainty
Any coupling measurement based on a Lagrangian description of underlying field theory assumes a set of operators, specifically those describing the Higgs interactions. In the Standard Model this includes the renormalizable Higgs Lagrangian with the tree-level interactions to massive gauge bosons and fermions. The effective Higgs couplings to photons and gluons are loop-induced, but they avoid the decoupling theorem and are only suppressed by 1/v. This means that at LHC energies they have to be included in the description of the Higgs signatures. Additional higher-dimensional operators can be included in the analysis once the corresponding measurements separate them from deviations in the renormalizable Higgs operators [27]. The prefactors of each operator, combined with a renormalization scheme, constitute the measurable couplings. When extracting the couplings of a light Higgs boson with mass 126 GeV we follow an effective theory approach. This allows us to only consider modifications of the light Higgs boson while keeping the ultraviolet properties of a functioning, partly decoupled Higgs sector [42]. We analyze extensions of the SM Higgs sector based on Higgs coupling modifications, written as [10,[42][43][44] For loop-induced Higgs couplings it is important that we first define the deviation due to a shift in the SM-like loops and then separately treat additional particles in the loop. The Higgs signal strength is given in Eq.(1.1). The three ratios of the production rate, the decay width, and the total width depend on one or more Higgs couplings. Because the total Higgs width cannot be measured at the LHC we assume that it is given by the sum of all observable partial widths, where we align the relevant second generation modes with the corresponding third generation decay measurements.
In different new physics models we describe the leading deviations from the Standard Model in terms of one free parameter ξ. For more details and a complete set of references we refer to Appendix B and to Ref. [42]. Because new physics effects which violate custodial symmetry are unlikely to be discovered in the Higgs sector, we simplify our example analysis by assuming ∆ W = ∆ Z ≡ ∆ V . In all models considered below we also find ∆ V < 0. For example two-Higgs-doublet models can violate both of these features at the loop level, so eventually we should release this simplification. We choose ξ such that in the custodial limit the (effective) Higgs coupling to gauge bosons gets modified as

2)
Dark singlet Dark singlet models include an additional scalar S which couples to the Higgs through dimension-4 portal interactions [45]. The extra scalar does not form a VEV, and its decays are precluded by a discrete Z 2 symmetry. Because it couples to the Higgs it leads to an invisible decay width and hence impacts the Higgs measurements through an invisible width, In this case the scaling pattern of Eq.(5.2) does not hold for the actual coupling g V , but for its apparent value from the rate measurement.
Additional singlet In the presence of an additional SU (2) L singlet with non-zero VEV, its mixing with the Higgs boson is described by an angle sin θ ≡ ξ [46]. All Higgs couplings to fermions and gauge bosons are rescaled by the common factor The phenomenological equivalence of the dark singlet, the singlet mixing, and the simplest strongly interacting single form factor models can only be broken by an observation of invisible Higgs decays.

Composite Higgs Minimal composite Higgs models (MCHM) describe the Higgs boson
as a pseudo-Nambu-Goldstone boson in a new strongly interacting sector with a spontaneously broken global symmetry [47,48]. In the Randall-Sundrum picture this global symmetry has to include the local gauge groups of the Standard Model and can in addition include a global custodial symmetry of the Standard Model. This way one of the states in the Higgs-sector will be light, while the others reside at f m H . The light Higgs couplings are shifted proportional to the ratio * ξ = v f .

(5.5)
Depending on the symmetry structure this Goldstone-protected strongly interacting Higgs sector predicts different coupling patterns for fermions and gauge bosons. In the MCHM5 setup the ratio of production rates scales like Additional doublet In the most general setup with an additional Higgs doublet [49], the Yukawa-aligned 2HDM, the different Higgs couplings to the heavy fermions vary independently. However, the light Higgs couplings to the massive gauge bosons are universally modified by the mixing angle cos(β − α) ≡ ξ, The decoupling parameter ξ parameterizes the distance from the Standard Model limit. In the fermion sector four setups accommodate the flavor symmetry of the Standard Model [50] -type-I, where all fermions couple to just one Higgs doublet Φ 2 * Note that for consistency reasons this definition of ξ differs from the original ξ = (v/f ) 2 .
-type-II, where up-type (down-type) fermions couple exclusively to Φ 2 (Φ 1 ) -lepton-specific, with a type-I quark sector and a type-II lepton sector -flipped, with a type-II quark sector and a type-I lepton sector These natural flavor conserving models correspond to particular cases of the aligned 2HDM with specific alignment angles. The modification factors ∆ t,b,τ depend on the mixing angles α and β. We vary sin(β − α) while fixing tan β = 1.5. We also fix the charged Higgs contribution to g γ via m H ± = 500 GeV andλ = 1 for the relevant self coupling defined in Ref. [42]. This choice leads to a small charged Higgs contribution with the same sign as the top loop. It sharpens the destructive interference with the W loop and yields µ p,γγ < 1 in the limit ξ → 0.
MSSM The Higgs sector of the MSSM is a particular case of a type-II 2HDM with a supersymmetric Higgs potential [51]. This identifies the Higgs self couplings with gauge couplings, so neither the additional heavy Higgs masses nor the mixing angle α are free quantities anymore. The decoupling limit can be described in terms of the general 2HDM, as defined in Eq.(5.7). The only difference is that for the MSSM the gauge couplings are related to the heavy Higgs masses as entering the gauge boson couplings just as in Eq.(5.7). These tree-level Higgs mass and coupling patterns may be strongly modified by the inclusion of quantum effects [52]. In the computation of the decoupling parameter ξ we therefore use α eff , defined as the mixing angle of the two scalar Higgs states with the loop-corrected mass matrix. The ratio of the two vacuum expectation values is fixed to tan β = 10 . The resulting size of ξ is limited by the range m A 0 = 200 − 1000 GeV, which will mean that we never reach ξ = 0.2. We use the MSSM Higgs boson cross sections and branching ratios given by Feyn- Higgs [53]. The Higgs signal strengths are defined after identifying the (lightest) Higgs masses in both models m SM H = m MSSM h 0 = 126 GeV. We show a benchmark scenario with a maximum light Higgs mass generated through large stop mixing (m max h ). Compared to the general 2HDM, the possible departures from the linear correlations are milder. This can be understood from the more constrained 2HDM potential in the MSSM case, which implies a fast transition to the decoupling regime [42]. The largest deviations arise in the low-m A 0 regime or for light stop and stau masses.
Signal strengths For the above described modifications to a Standard-Model-like Higgs sector we find simple patterns in the two-dimensional plane of coupling strengths. In Figure 5 we compare the correlated modifications as functions of the decoupling parameter ξ. In general, the behavior in the H → V V and H → γγ decay planes should be similar, as long as the loop-induced Higgs-photon coupling is dominated by the W loop. We focus on ξ < 0.4, corresponding to a modification of g V by 8%, which the LHC is sensitive to with sufficient luminosity. In addition, we mark a deviation by ξ = 0.2, equivalent to a 2% coupling deviation. The latter could be considered the target of a linear collider analysis.
For the simplest models, the dark singlet and the singlet mixing, all correlations follow a straight diagonal line towards reduced coupling strengths µ p,d < 1. This is due to the simple mixing pattern and the net scaling of the LHC event rates as σ×BR ∝ g 2 . The same pattern appears for the simplest strongly interacting models with a single Higgs form factor. For the more complex strongly interacting model MCHM5, we find both correlated (µ GF,τ τ − µ VBF,τ τ ) and anti-correlated (µ GF,V V −µ VBF,V V ) patterns. It predicts an increased number of weak boson fusion events whenever the couplings ∆ V in the production process decrease more slowly than the coupling ∆ f in the total width, as shown in Eq. (5.6).
Larger departures from the Standard Model are possible in the 2HDM, thanks to its more flexible coupling structure. In the type-I setup there is essentially no mechanism to increase the number of weak boson fusion events as compared to the Standard Model, because of the increase in the Higgs width combined with the reduced gauge boson coupling. This is described in more detail in Appendix B. For gluon fusion production combined with a fermionic decay the suppression by the total width can be compensated by the production and decay couplings. In the type-II setup both signal strengths can, unexpectedly, be enhanced for bosonic Higgs decays. The reason is a strongly decreased partial Higgs width to bottoms and taus which cannot be generated in the type-I model. For fermionic decays the direct link between the bottom and tau down-type Yukawas leads to a systematically decreased event rate. The fact that for ξ → 0 the 2HDM rates do not match the Standard Model is linked to the finite contribution of the charged Higgs to the effective photon-Higgs coupling.
Finally, the MSSM as a constrained type-II 2HDM shows limited signal strength variations because of the supersymmetric constraints. Unlike the general 2HDM even in the type-II setup the MSSM does not allow for a free variation of the two parameters α and β, which affect the Yukawa couplings in a complicated manner. Departures from the decou-pling limit in this case lie below ξ 0.2 for the considered parameter space configurations. In the MSSM deviations from the Standard Model in the limit ξ → 0 arise through contributions of the sfermions and the charged Higgs to both the effective gluon-Higgs and photon-Higgs couplings. In that sense the parameter ξ does not fully track down the decoupling limit for the loop-induced Higgs couplings, similarly to the 2HDM case. The effect of a shifted bottom Yukawa is not sufficient to overcome the reduction in g V , which means that unlike for the 2HDM both signal strength deviations for the V V and γγ decays are negative, leading to an (almost linear) correlated suppression. Moreover, we see that the typical deviations in the signal strengths can be achieved for small values ξ < 0.2 in the MSSM, because quantum effects dominate over the mere tree-level rescaling ∆ V .
A Heuristic for Robustness to Theory Uncertainty One of the great challenges of the Higgs coupling program is to understand what type of deviation from the Standard Model prediction would be compelling enough to make a claim for new physics. It is clear that the pattern of deviations in the various production and decay modes carries much more information than considering them individually. Furthermore, the inability to measure the total width of the Higgs necessitates either assumptions on the total width or consideration of various ratios in which the total width cancels. Both approaches lead to strong correlations in the inferred couplings. If the theoretical uncertainties were well defined and statistical in nature, the significance of any given deviation could be readily assessed by standard statistical methods. However, the ill-defined nature of theoretical uncertainties is beyond the scope of rigorous statistical procedures.
For example, if one were to see a 4σ deviation from the Standard Model that could be reduced to a 2σ deviation by inflating a theoretical uncertainty by a factor of two or by changing from a Gaussian constraint to an Rfit constraint, then the deviation would most likely be met with healthy skepticism by the community. However, a deviation that is orthogonal to the effect of a theoretical uncertainty is much more robust.
This motivates a heuristic to evaluate the robustness of an observed deviationμ to a theoretical uncertainty parametrized by α i , which we denote R i (µ). We want a larger deviation from the Standard Model to be reflected as a larger value for this robustness heuristic, so we begin with the length of the vector |µ − 1|. We also want the robustness to be larger as the deviation from the Standard Model (µ − 1) becomes more orthogonal to the shift in the signal strength induced from varying α i as defined in Eq.(2.21) and denoted as ∂ α i µ fix . Since the magnitude of the theoretical uncertainty is poorly defined, it is natural that we only consider the angle between (µ − 1) and ∂ α i µ fix . This leads to the robustness heuristic This same heuristic can be used to assess an expected departure from the Standard Model based on a new physics model parametrized by ξ via composition R i (µ(ξ)). Figure 6 shows the heuristic for the new physics effects in Fig. 5 combined with the theoretical uncertainties in the toy model illustrated in Fig. 2. As expected from the discussion of the different models, the different two-Higgs-doublet models can be distinguished from the  Figure 6. The sensitivity heuristic R i (ξ) evaluated for various new physics models and the theoretical uncertainties i associated to the gluon fusion cross section for ≥ 0-jets and ≥ 2-jets.
Standard Model even for ξ ∼ 0. While for the type-I model the VBF topology is more robust with respect to potential QCD effects, deviations due to a type-II model are more robust to the inclusive Higgs production rate. The most reliable signature for the strongly interacting MCHM5 model would be observed in weak-boson-fusion Higgs production with a decay H → W W . For the MSSM the ξ axis is rescaled, because deviations ξ 0.1 are hardly generated in our scan over MSSM spectra. The robustness of a supersymmetric Higgs sector is roughly equal for the inclusive and VBF Higgs topologies.

Conclusion
Motivated by the fact that ill-defined theoretical uncertainties will eventually be the limiting factor in Higgs coupling measurements at the LHC, we have developed a technique to decouple the theoretical uncertainties from the experimental results while retaining the ability to incorporate those uncertainties in a subsequent stage we refer to as recoupling.
This approach is amenable to simultaneously measuring multiple quantities, such as a vector of signal strength parameters µ for different Higgs production and decay signatures. Moreover, the technique lends itself well to combinations with several common sources of uncertainty that induce correlations among the contributing measurements. In that respect it is similar to the Best Linear Unbiased Estimator (BLUE) [11] technique, but not restricted to Gaussian measurements or linear response to the source of uncertainty. We considered a toy example modeled after the current ATLAS Higgs coupling measurements where the measurements are not in the Gaussian regime and these non-linear effects are important for approximating the full likelihood function.
One of the most powerful features of this approach is that it allows one to change the assumptions on both the magnitude and the shape of the uncertainty in the recoupling stage, which may occur long after the experimental groups have released their results. This includes the ability to introduce a priori correlations in the source of the systematics, which might have been neglected originally. These same capabilities would be possible if the experiments published the full statistical model using RooFit/RooStats [26] as suggested in Refs. [54]; however, the approach outlined here is less technology dependent. The conceptual picture of these uncertainties leading to shifts in the inferred values of the parameters is intuitive and provides convenient visualizations as in Fig. 2.
In order for the experiments to present their results in this way, they would need to • publish the effective likelihood L eff (µ eff ) profiling only parameters that are not anticipated to be common to other measurements, • publish the reparametrization template µ eff (µ, α), and • document the conventions that establish meaning to the nuisance parameters α.
In return, our approach allows for a flexible treatment of systematic uncertainties and removes the burden of choosing a description of theoretical uncertainties from the experimental groups. It also allows for future improvements in the theoretical description of Higgs processes at the LHC to be easily incorporated. We discussed multiple strategies to determine the reparametrization template in Sec. 2.
In the absence of a reliable measure of theoretical uncertainties the key question becomes how easily an apparent deviation of experimental measurements from the Standard Model description can be explained by a change in the assumed theoretical uncertainties. If we describe new physics effects in the Higgs sector as one-parameter deviations from the Standard Model decoupling limit in the signal strength planes, the directions of these deviations can be compared to the effects of a change in the theoretical uncertainties in the same planes. The direct comparison of these two possible explanations for an experimental observation leads us to the robustness heuristic presented in Eq.(5.9).
While this paper focuses on the application of our approach to Higgs coupling measurements and the theoretical uncertainties associated with QCD, the technique is quite general and may find broad applications.

A A Worked Example
In this Appendix we work through a simple example that illustrates explicitly the procedure described in Sec. 2. A Mathematica notebook carrying out these calculations can be found in Ref. [41]. We consider two Gaussian measurements (indexed by c = 1, 2), which can be thought of as approximating the Poisson distribution for number counting analyses or approximating the maximum likelihood estimator of some more complicated analysis. We consider two signal processes (indexed by p = 1, 2) with nominal (e.g. Standard Model) expectations s cp and signal strength modifiers µ p . Finally, we consider one systematic effect parametrized by α that shifts the signal expectation from the nominal α 0 = 0 so that the expectation is Including a Gaussian constraint term for the parameter α leads to the full likelihood The maximum likelihood estimators arê The full Fisher information matrix is straightforward to calculate, but cumbersome to write explicitly. Results below only require the off-diagonal block elements In order to decouple the uncertainty, the experiment would provide the likelihood for an effective signal strength with respect to the nominal prediction α 0 = 0. The effective likelihood would be based on which has the same maximum likelihood estimates above and the following information matrix The effective likelihood has the form: Now we must choose a reparametrization template to be used to recouple the the uncertainty due to α. First let us choose the template of Eq.(2.9), so that µ eff p = µ p (1+η p α). Note, that the original model in Eq.(A.2) had four η cp while the template only has two η p . If the effect of α is category universal so that η c=1,p = η c=2,p , then recoupling based on the this template will reproduce the full model exactly; however, in the more general situation η c=1,p = η c=2,p it will not. Now we solve for the coefficients of the template based on the local covariance matrix as described in Sec. 2.4. We outlined three equivalent approaches based on Eqs. (2.21), (2.24), and (2.25). Let us demonstrate the last of these three approaches. Based on the template of Eq.(2.9), the reparametrization will lead to the Jacobian As in Eq.(2.25), we use this Jacobian to relate the information matrix of the effective likelihood and the main measurement The (i = µ p , j = α) sub-block of this matrix leads to the following system of linear equationŝ which can easily be inverted to provide solutions for η p . In practice one would work with numerical representations of the maximum likelihood estimators and Fisher information matrices, but here we present the result symbolically. Note, in the category-universal situation, the solution simplifies to η p = η 1p = η 2p as expected. In the non-category-universal situation it is possible for η p to lie outside of the range [η 1p , η 2p ].
Three example scenarios To make these examples more explicit, Table 1 specifies three scenarios for the coefficients in Eq.(A.1). Scenarios A and B are meant to be representative of LHC Higgs measurements in which the first category is gluon-fusion-like with significant background and signal dominated by the first production mode (p = 1 = ggF), while the second category is VBF-like with negligible background and signal dominated by the second production mode (p = 2 = VBF). The only difference between Scenarios A and B are the η cp that quantify the response to the source of uncertainty. Since Scenario A is category-universal we achieve an exact reproduction of the full likelihood, while for Scenario B we don't expect exact results. Scenario C is meant to probe the extreme case in which the uncertainty for the first production mode only affects the expectation in the second category (which is representative of the theory uncertainty of gluon-fusion + ≥ 2jets.  Figure 7 shows a comparison of the full and recoupled likelihood using the 'aligned' template of Eq.(2.9) and the solutions to the coefficients in Eq.(A.11). Scenario A is reproduced exactly, the non-category-universal property of Scenario B leads to a slight discrepancy, and the extreme non-category-universal property of Scenario C leads to a substantial discrepancy. In particular, while the effect of the uncertainty in Scenario C is to change the inferred value of µ VBF , the size of this effect should scale with µ ggF . The agreement can be improved dramatically by moving to the more general template of Eq.(2.13) that includes η p pi . One approach is to fix 'by hand' the coefficients of the template η ggF ggF = 0, η VBF VBF = 0 and φ = 0 and determine the two remaining coefficients using the local covariance matrix (which proceeds as above with a different Jacobian transformation).
An alternate approach is to use the unrestricted template of Eq.(2.13) and utilizing the 'learning' approach of Eq. (2.26).
Finally, the second row of Fig. 7 shows a comparison of the full and decoupled likelihood with a modified constraint term. In particular, the Gaussian constraint is replaced with an Rfit constraint term: Gaus(0|α, 1) → Uniform(−1, 1). The coefficients for the templates are the same as for both rows of Fig. 7

B New Physics Models
In this appendix we will give a more detailed picture of the new physics models and their features briefly discussed in Section 5. In particular, we will motivate and discuss the description of new physics effects by a single parameter ξ, defined as the modification of the Higgs couplings to massive gauge bosons, i.e. ∆ V −ξ 2 /2. Note that this unified definition of ξ differs from Ref. [42] for some of the new physics models. All signal strength deviations we compute by rescaling the SM production cross section, branching ratio and total width [12], while for the MSSM case we use FeynHiggs [53].
Dark singlet A dark singlet is defined as a model with an additional scalar particle S which does not have a vacuum expectation value and hence cannot mix with the Higgs boson. In addition, we assume that its only interaction with the Standard Model will be the dimension-4 portal interactions in the combined scalar potential [45], This interaction with strength λ 3 can, if kinematically allowed, lead to an invisible Higgs decay width [55] directly linked to a possible dark matter agent [56] Such an invisible Higgs width contributes to the total Higgs width and hence to the number of predicted Higgs events at the LHC. This leads to an apparent reduction of all couplings, including the Higgs coupling to massive gauge bosons, shown in Eq. (5.3). Such a universal modification of all predicted event numbers corresponds to a diagonal pattern, for example in the two-dimensional µ VBF,d vs µ GF,d plane. In the upper panels of Figure 8 we show these identical diagonal correlations for a set of decay-diagonal channels (left panels) and non-diagonal channels (right panels). The correlations are all identical. In the lower panels we show the expected signal strengths relative to the Standard Model as a function of the invisible width.
Additional singlet If the additional SU (2) L singlet acquires a finite VEV v S the combined Higgs potential [46,57] V with the portal interaction λ 3 leads to singlet-doublet mixing. The rotation to mass eigenstates h and H defines the angle Dark singlet µ GF,γγµ GF,ττ µ VBF,γγµ VBF,ττ µ HV,bbµ VBF,ττ µ p,d Figure 8. Dark singlet: correlated Higgs signal strengths for the decay-diagonal channels (left) and non-diagonal channels (right). The first of the two signal strengths in the notation µ p1,d1 − µ p2,d2 is shown on the vertical axis, the second on the horizontal axis. In the lower panels we give the different signal strengths as a function of Γ inv /Γ SM = ξ 2 .
All Higgs couplings to fermions and gauge bosons are rescaled by cos θ < 1. In addition, the mostly Higgs state can decay into two lighter, mostly singlet states [58] Γ The signature for such a decay depends on the lifetime and the decay channels of the lighter state h. If sizeable, this additional decay channel contributes to the total Higgs width, entering the predicted number of LHC events as a second universal modification. This universal coupling modification again predicts diagonal lines for all µ p 1 ,d 1 − µ p 2 ,d 2 correlations, as shown in Figure 9. If we do not observe the new Higgs decay modes for each of the additional singlet scenarios, the degeneracy in the signal strength deviations makes it impossible to distinguish a dark singlet, an additional singlet, and the simplest strongly interacting form factor models (as discussed in the next paragraph) [59].
Composite Higgs Depending on the symmetry structure this strongly interacting but Goldstone-protected Higgs sector predicts different coupling patterns for fermions and gauge bosons [48]. The MCHM4 setup with ∆ V = ∆ f = 1 − ξ 2 − 1 is experimentally equivalent to a mixing term in a singlet extension. The reduced couplings of the lightest Higgs state reflect the fact that at the energy scale f there exist many strongly interacting Higgs fields which share the unitarization of the usual 2 → 2 scattering processes. The MCHM5 setup is more interesting. While ∆ V is identical to the MCHM4 case the fermions follow a different pattern, 1 + ∆ f = (1 − 2ξ 2 )/ 1 − ξ 2 . If we ignore the heavy states' contributions to the the effective Higgs-photon and Higgs-gluon couplings we find that the Singlet µ GF,γγ -µ GF,ττ µ VBF,γγ -µ VBF,ττ µ HV,bb -µ VBF,ττ µ p,d Figure 9. Singlet mixing: correlated Higgs signal strengths for the decay-diagonal channels (left) and non-diagonal channels (right). The first of the two signal strengths in the notation µ p1,d1 −µ p2,d2 is shown on the vertical axis, the second on the horizontal axis. In the lower panels we give the different signal strengths as a function of cos θ = 1 − ξ 2 . ratio of production rates scales like Unlike all previous models the MCHM5 setup accommodates µ p,d > 1, for example The sharper suppression in the fermion couplings still leads to depleted gluon fusion channels. In Figure 10 we Additional doublet An even richer structure in particular in separating the gauge and Yukawa couplings of the light Higgs appears when we add a second Higgs doublet [49,50]. Up to terms of O(ξ 3 ) the modification to the LHC signal strengths (1 + ∆ p ) 2 (1 + ∆ d ) 2 − 1 reads  Figure 10. Composite Higgs: correlated Higgs signal strengths for the decay-diagonal channels (left) and non-diagonal channels (right). The first of the two signal strengths in the notation µ p1,d1 − µ p2,d2 is shown on the vertical axis, the second on the horizontal axis. In the lower panels we give the different signal strengths as a function of These expressions hold for the general Yukawa-aligned model [42,60]. For example a type-I model corresponds to γ b = γ τ = π/2, while a type-II model appears if we set γ b = γ τ = 0. The main feature is that the leading signal strength deviations involving fermions arise at order ξ, while gauge couplings only vary with ξ 2 [61]. The detailed 2HDM signal strength patterns are to a large extend model-dependent, as they are tied to the specific Yukawa structures. In the type-I setup all Yukawas couplings and hence the effective Higgs-gluon coupling are shifted by a common factor cos α/ sin β. This way they are not suppressed for sin(β − α) 1. The Higgs couplings to gauge bosons scale like sin(β − α). As a result, all gluon fusion channels show an increased signal strength as compared to their weak boson fusion counterparts. Type-II models link up-like and down-like fermions to different Higgs doublets, implying separate Yukawa modifications. The bottom and tau final states are suppressed for a wide range of sin(β −α) -eventually leading to the unphysically large signal strength deviations visible in the plots. For an even larger deviation from sin(β − α) 1 models of type-II, lepton-specific and flipped feature sign-inverted Yukawas. These sign ambiguities are visible in the ellipsoidal correlated variations in the signal strength plane. The additional charged Higgs-mediated contribution ∆ γ is responsible for the mild offset µ p,γγ = 1 in the decoupling limit.
MSSM A supersymmetric Lagrangian requires two Higgs doublets to give mass to uptype and down-type fermions. At tree level the MSSM Higgs sector is a type-II 2HDM, but with the different quartic Higgs couplings fixed to gauge couplings. If we assume the lighter of the scalar Higgs bosons to lie at 126 GeV the heavy Higgs masses are almost 2HDM -Type II µ HV,bb -µ VBF,ττ µ HV,bb µ GF,γγ -µ GF,ττ µ VBF,γγ -µ VBF,ττ µ VBF,ττ µ GF,γγ µ GF,ττ µ VBF,γγ Figure 11. 2HDM: correlated Higgs signal strengths for decay-diagonal channels (left two panels) and non-diagonal channels (right two panels), each for type-I and type-II setups. The first of the two signal strengths in the notation µ p1,d1 − µ p2,d2 is shown on the vertical axis, the second on the horizontal axis. In the lower panels we give the different signal strengths as a function of cos(β − α) = ξ.
degenerate at m H 0 m A 0 m H ± up to O(m 2 Z /m 2 A 0 ) corrections. Just like for the general 2HDM the decoupling is described by ξ = cos(β − α).
At tree level the entire MSSM Higgs sector is fully described by m A 0 and tan β. Loop corrections in the top-stop sector lead to significant corrections and yield an additional parameter dependence for example on the stop trilinear coupling A t . The leading quantum corrections to the Higgs couplings can be accounted for through an effective mixing angle α → α eff which we use in the definition of ξ. Additional contributions may further modify the effective Higgs coupling pattern. One-loop triangle corrections to the Yukawa couplings are governed by gluino-sbottom loops [62] and shift the bottom quark Yukawas by ∆ b ∼ O(m 2 Z tan 2 β/m 2 A 0 ). They give rise to a delayed decoupling [51], meaning that the decoupling limit now demands m 2 Z tan 2 β m 2 A 0 for large tan β. These effects depend on the detailed supersymmetric mass spectrum and have to be included in a full MSSM parameter study [63]. In this discussion of the light Higgs couplings we omit any constraint from the new particles outside the Higgs sector. This also means that we ignore the possibility of light neutralinos 2m χ 0 1 < m h 0 contributing to the invisible Higgs width. In specific supersymmetric models the coupling of the lightest neutralino to the light Higgs boson can be important, leading to an accelerated dark matter annihilation in the early universe. Nevertheless, even the couplings required for this light Higgs funnel are unlikely to give a measurable invisible Higgs decay rate at the LHC.
Technically, we use the MSSM Higgs cross sections and branching ratios given by FeynHiggs 2.9.5 [53]. The corresponding Higgs signal strengths are obtained by normalizing the individual production and decay rates to their SM counterparts [12], after identifying the (lightest) Higgs mass in both models m SM H = m MSSM h 0 . We consider two benchmarks relevant to LHC searches [64], all compatible with the observed ∼ 126 GeV resonance: the m max h case with maximum stop mixing, heavy squarks, and tan β = 10 and the light stop case with maximum stop mixing, but generally light squarks, and tan β = 35. Both of them are numerically similar. The difference is that only the light stop scenario  Figure 12. MSSM: correlated Higgs signal strengths for decay-diagonal channels (left two panels) and non-diagonal channels (right two panels), each for two MSSM benchmark spectra. The first of the two signal strengths in the notation µ p1,d1 − µ p2,d2 is shown on the vertical axis, the second on the horizontal axis. In the lower panels we give the different signal strengths as a function of m A , linked to ξ via Eq.(5.8).
generates sizeable O(10%) corrections to ∆ g and at the same time gives relatively large negative ∆ b corrections to the bottom Yukawa. In Fig. 12 we see that compared to the general 2HDM, the possible departures from the linear correlations are milder. The two benchmarks predict very similar signal strengths. All this can be understood from the more constrained Higgs potential and the moderate negative corrections to the total width. The largest deviations obviously arise in the low-m A 0 regime.