Constraining new physics from Higgs measurements with Lilith: update to LHC Run 2 results

Lilith is a public Python library for constraining new physics from Higgs signal strength measurements. We here present version 2.0 of Lilith together with an updated XML database which includes the current ATLAS and CMS Run 2 Higgs results for 36/fb. Both the code and the database were extended from the ordinary Gaussian approximation employed in Lilith-1.1 to using variable Gaussian and Poisson likelihoods. Moreover, Lilith can now make use of correlation matrices of arbitrary dimension. We provide detailed validations of the implemented experimental results as well as a status of global fits for reduced Higgs couplings, Two-Higgs-doublet models of Type I and Type II, and invisible Higgs decays. Lilith-2.0 is available on GitHub and ready to be used to constrain a wide class of new physics scenarios.


Introduction
The LHC runs in 2010-2012 and 2015-2018 have led to a wealth of experimental results on the 125 GeV Higgs boson. From this emerges an increasingly precise picture of the various Higgs production and decay processes, and consequently of the Higgs couplings to the other particles of the Standard Model (SM), notably gauge bosons and third generation fermions. With all measurements so far agreeing with SM predictions, this poses severe constraints on scenarios of new physics, in which the properties of the observed Higgs boson could be affected in a variety of ways.
Assessing the compatibility of a non-SM Higgs sector with the ATLAS and CMS results requires to construct a likelihood, which is a non-trivial task. While this is best done by the experimental collaborations themselves, having at least an approximate global likelihood is very useful, as it allows theorists to pursue in-depth studies of the implications for their models. For this reason, the public code Lilith [1] (see also [2]) was created, making use of the Higgs signal strength measurements published by ATLAS and CMS, and the Tevatron experiments. 1 In this paper, we present version 2.0 of Lilith together with an updated database which includes the current published ATLAS and CMS Run 2 Higgs results for 36 fb −1 (status early June 2019).
Compared to HiggsSignals [4], which is written in Fortran90 and uses the signal strengths for individual experimental categories with their associated efficiencies as well as Simplified Template Cross Section (STXS) [5,6] results, Lilith is a light-weight Python library that uses as a primary input signal strength results in which the fundamental production and decay modes are unfolded from experimental categories. Here, the main production mechanisms X are: gluon fusion (ggH), vector-boson fusion (VBF), associated production with an electroweak gauge boson (WH and ZH, collectively denoted as VH) and associated production with top quarks, mainly ttH but also tH. The main decay modes Y accessible at the LHC are H → γγ, H → ZZ * → 4 , H → W W * → 2 2ν, H → bb and H → τ τ (with ≡ e, µ). 2 The signal strength framework is based on the narrow-width approximation and on the assumption that new physics results only in the scaling of SM Higgs processes. This makes it possible to combine the information from various measurements and assess the compatibility of given scalings of SM production and/or decay processes from a global fit to the Higgs data. This framework is very powerful as it can be used to constrain a wide variety of new physics models, see for example [7] and references therein.
For a proper inclusion of the recent Run 2 results from ATLAS and CMS, several improvements were necessary in Lilith. Concretely, in order to treat asymmetric uncertainties in a better way, we have extended the parametrisation of the likelihood to Gaussian functions of variable width ("variable Gaussian") as well as generalised Poisson functions. Moreover, Lilith can now make use of correlation matrices of arbitrary dimension. We have also added the tH and the gluon-initiated ZH production modes, and corrected some minor bugs in the code.
Results given in terms of signal strengths can be matched to new physics scenarios with the introduction of factors C X and C Y that scale the amplitudes for the production and decay of the SM Higgs boson, respectively, as for the different production modes X ∈ (ggH, VBF, WH, ZH, ttH, . . .) and decay modes Y ∈ (γγ, ZZ * , W W * , bb, τ τ , . . .), where the sum runs over all decays that exist for the SM Higgs boson. The factors C X and C Y can be identified to (or derived from) reduced couplings appearing in an effective Lagrangian. Following [8] and subsequent publications, we employ the notation 1 For a discussion of the use and usability of signal strength results, and recommendations on their presentation, see [3]. 2 When µ(X, Y ) is not directly available, µ = effX,Y × µ(X, Y ) is used, introducing appropriate efficiency factors effX,Y . For inclusive combinations of production channels for the same Y , the efficiencies become effX = σ SM (X)/ σ SM (X).
where C W,Z and C f (f = t, b, c, τ, µ) are bosonic and fermionic reduced couplings, respectively. In the limit where all reduced couplings go to 1, the SM case is recovered. In addition to these tree-level couplings, we define the loop-induced couplings C g and C γ of the Higgs to gluons and photons, respectively. If no new particles appear in the loops, C g and C γ are computed from the couplings in Eq. (3) following the procedure established in [9]. Alternatively, C g and C γ can be taken as free parameters. Apart from the different notation, this is equivalent to the so-called "kappa framework" of [9]. Finally note that often a subset of the C's in Eq. (3) is taken as universal, in particular C V ≡ C W = C Z (custodial symmetry), C U ≡ C t = C c and C D ≡ C b = C τ = C µ like in the Two-Higgs-doublet model (2HDM) of Type II, or C F ≡ C U = C D as in the 2HDM of Type I.
Last but not least, while the signal strength framework in principle requires that the Higgs signal be a sum of processes that exist for the SM Higgs boson, decays into invisible or undetected new particles, affecting only the Higgs total width, can be accounted for through without spoiling the approximation. The rest of the paper is organised as follows. We begin the discussion of novelties in Lilith-2.0 by presenting the extended XML format for experimental input in Section 2. This is followed by details on the calculation of the likelihood in Section 3. The inclusion of new production channels is described in Section 4. The Run 2 results included in the new database and their validation are presented in Section 5. In Section 6 we then give an overview of the current status of Higgs coupling fits with Lilith-2.0. We conclude in Section 7. Appendix A contains additional material illustrating the importance of various improvements discussed throughout the paper.
It is important to note that this paper is not a standalone documentation of Lilith-2.0. Instead, we present only what is new with respect to Lilith-1.1. For everything else, including instructions how to use the code, we refer the reader to the original manual [1].

Extended XML format for experimental input
In the Lilith database, every single experimental result is stored in a separate XML file. This allows to easily select the results to use in a fit, and it also makes maintaining and updating the database rather easy.
So far, the likelihood information could be specified in one or two dimensions in the form of [1]: 1D intervals given as best fit with 1σ error; 2D likelihood contours described as best fit point and parameters a, b, c which parametrise the inverse of the covariance matrix; or full likelihood information as 1D or 2D grids of −2 log L. The first two options, 1D intervals and 2D likelihood contours, declared as type="n" in the <expmu> tag, employ an ordinary Gaussian approximation; in the 1D case, asymmetric errors are accounted for by putting together two one-sided Gaussians with the same mean but different variances, while the 2D case assumes symmetric errors. This does does not always allow to describe the experimental data (i.e. the true likelihood) very well.
In order to treat asymmetric uncertainties in a better way, we have extended the XML format and likelihood calculation in Lilith to Gaussian functions of variable width ("variable Gaussian") as well as generalized Poisson functions [11]. The declaration is type="vn" for the variable Gaussian or type="p" for the Poisson form in the <expmu> tag. Both work for 1D and 2D data with the same syntax. Moreover, in order to make use of the multi-dimensional correlation matrices which both ATLAS and CMS have started to provide, we have added a new XML format for correlated signal strengths in more than two dimensions. This can be used with the ordinary or variable Gaussian approximations for the likelihood. In the following we give explicit examples for the different possibilities.
To ensure backwards compatibility, type="n" however still requires the tags <a>, <b>, <c> to give the inverse of the covariance matrix instead of <uncertainty> and <correlation>, see [1].
We are aware that having different formats for 2 and more than 2 dimensions is not necessary in principle. Nonetheless we chose to treat the 2D case separately (with axis tags "x" and "y" instead of "d1" and "d2") to stay as close as possible to what was done in Lilith-1.1. This may change in future versions.

Likelihood calculation
The statistical procedure used in Lilith was described in detail in [1]. The main quantity given as an output is the log-likelihood, −2 log L, evaluated in computelikelihood.py from the information given in the XML database. In this section, we explain how −2 log L is computed for type="vn" (variable Gaussian) and type="p" (Poisson) introduced in the previous section. For the old implementation of the ordinary Gaussian (type="n"), we refer the reader to [1].

Variable Gaussian
As shown in [11], a Gaussian function of variable width can be a good choice to deal with asymmetric uncertainties. We use the version linear in the variance, described as "Variable Gaussian (2)" in Section 3.6 of [11]. In the 1D case, the likelihood is then given by whereμ denotes the best-fit signal strength, and σ − and σ + are absolute values of the left and right uncertainties at 68.3% CL, respectively. If not stated otherwise, these notations are used for the entire section. For σ + = σ − , the symmetric case is obtained. The variable Gaussian form however has a singularity at µ =μ − (σ + σ − )/(σ + − σ − ), which can lead to numerical issues, although in practice this usually happens for µ's outside the range of interest (or reduced couplings outside their physically meaningful range). The numerical stability is however problematic when σ − → 0; in this case it can be better to use the ordinary 2-sided Gaussian implementation.
In the case of n-dimensional data (n > 1), we use the correlations given by the experimental collaboration, if available, together with the best fit points and the left and right uncertainties at 68.3% CL. When results are given in terms of 2D contour plots, we can also use the variable Gaussian form to numerically determine the best-fit point, uncertainties and their correlation, if not given explicitly by the experimental collaboration. For the n dimensional signal strength vector µ = (µ 1 , . . . , µ n ), the likelihood reads where the best fit pointμ = (μ 1 , . . . ,μ n ) and the covariance matrix is constructed from the correlation matrix ρ as with Here the σ − i and σ + i are the left and right uncertainties at 68.3% CL of the ith combination of production and/or decay channels, respectively. For multi-dimensional data in the ordinary Gaussian approximation, the relation between covariance matrix and the correlation matrix becomes where σ + = diag(σ + 1 , . . . , σ + n ) and σ − = diag(σ − 1 , . . . , σ − n ) .

Generalised Poisson
As an alternative to the variable Gaussian, a generalised Poisson form can be used for 1D and 2D results. For the 1D case, the likelihood is implemented according to Section 3.4, "Generalised Poisson", of [11] as where γ and ν are determined numerically from the equations .
More concretely, γ is determined from the expression on the left by bifurcation between γ = 0 and γ = 1/σ − and then inserted in the expression on the right to compute ν.
For the 2D case, we use the conditioning bivariate Poisson distribution described in [13], that has no restriction on the sign and magnitude of the correlation ρ. Here the joint distribution is a product of a marginal and a conditional distribution. The decision of which channel belongs to the marginal or the conditional distribution is based on the validation plots. To illustrate the method, we assume that the data of the channel X follows the marginal distribution, while data of the channel Y belongs to the conditional distribution. The joint log-likelihood is then and The function f reads where α is solved numerically from the correlation expression and the γ X , ν X and γ Y , ν Y are solutions of Eq. (11) for the X and Y channels, respectively.

New production channels
As mentioned in the introduction, we have also included a few new production channels. These are ZH production via gluon-gluon fusion (ggZH), Higgs production in association with a single top quark (tH), and Higgs production in association with a pair of bottom quarks (bbH). For the ZH production mode, the original implementation in Lilith included only the qq → ZH channel (qqZH). However, the loop-induced gluon-gluon fusion is not so small, about 14% of the total pp → ZH cross section at √ s = 13 TeV, and should hence be taken into account. Indeed, both ATLAS and CMS have been always including the ggZH contribution in their fits. From version 2 onwards, the ZH signal in Lilith is also the combination of the qq and gg initiated processes, with relative weights σ SM X /( X σ SM X ). In terms of the scale factors of Eq. (2), this gives For the SM cross sections the values given in [6,14] are used. The definitions of VH (ZH and WH) and VVH (ZH, WH and VBF) follow straightforwardly. The tH production also includes two contributions: t-channel tHq production and tHW production. The s-channel tHq cross section is much smaller and hence not included. Interference effects between these channels are also neglected. At √ s = 13 TeV, tHq is dominant, contributing about 80% of the tH cross section. As above, tHq and tHW are combined to tH production with efficiencies calculated from the SM cross sections. Moreover, following the usage in the ATLAS analysis [15], a combination of tH and ttH named 'top' production is defined in an analogous way. 5 Thus The values for σ SM ttH and σ SM tHq are again taken from [6,14], while σ SM tHW has been calculated at leading order using MadGraph [16] with factorization and renormalization scales µ F = µ R = (m t + m W + m H )/2 and the NNPDF30 lo as 0130 PDF set [17]. Other input parameters are the same as in [6] section I.6 for tH production. It is noted that the definition of a next-to-leading order cross section for the tHW channel is not straightforward because of the interferences with the ttH channel, see [18] for a discussion on this issue.
For the sake of completeness and to prepare for future data, the bbH production mode has also been added. The implementation is straightforward, anlogous to the ttH case.
Finally, since Lilith accepts reduced couplings as user input, the relations between the scaling factors and the reduced couplings are needed. These relations read, following the notation of [1], where the coefficients a i , e i , and f i are provided in Tables  5 Accordingly, prod="tHq", "tHW", "tH" and "top" can be used in the XML files in the database.  The ATLAS Run 2 results included in this release are summarised in Table 3   H → γγ (HIGG-2016-21): The ATLAS analysis [15] provides in Fig. 12 signal strengths for H → γγ separated into ggH, VBF, VH and "top" (ttH+tH) production modes. No correlations are given for the signal strengths, but we can use instead the correlations for the stage-0 STXS provided in Fig. 40a of the ATLAS paper, which should be a close enough match. It turns out, however, that the µ values rounded to one decimal do not allow to reproduce very well the ATLAS coupling fits for (C V , C F ) or (C γ , C g ). We have therefore extracted the best-fit points and uncertainties from the 1D profile likelihoods, which are provided as Auxiliary Figures 23a-d Fig. 12 of [15], but using more digits improves the coupling fits as shown in Fig. 1.
H → ZZ * → 4l (HIGG-2016-22): A similar issue as discussed for H → γγ above arises for H → ZZ * . In order to reasonably reproduce the C F vs. C V fit of ATLAS (Fig. 8b of [19]), we fit the 1D profile likelihoods for µ(ggH, ZZ * ) and µ(VBF, ZZ * ) shown in Aux. Figs. 7a and 7b of [19] as Poisson distributions. This gives µ(ggH, ZZ * ) 1.12 +0.25 −0.22 and µ(VBF, ZZ * ) 3.88 +1.75 −1.46 , which we implement as a bivariate Poisson distribution with correlation ρ = −0.41 (from Aux. Fig. 4c of [19]). For the VH and ttH production modes, lacking more information, we convert the given 95% CL limits into µ(VH, ZZ * ) = 0 ± 1.89 and µ(ttH, ZZ * ) = 0 ± 3.83 using a 2-sided Gaussian (assuming 1-sided limits gives a less good agreement with the ATLAS C F vs. C V fit). The validation is shown in Fig. 2 (see also Fig. 14 in Appendix A). This is a case where the variable Gaussian approximation performs less well than the Poisson likelihood.  Table 9 Lilith-2.0, DB 19.06 ATLAS official Figure 2: Fit of C F vs. C V for data from the ATLAS H → ZZ * analysis, using µ(ggH, ZZ * ) and µ(VBF, ZZ * ) as fitted from Aux. Figs. 7a and 7b of [19]; the ggH vs. VBF likelihood is then approximated as a bivariate Poissonian with correlation −0.41 (see text for more details). The 68%, 95% and 99.7% CL regions obtained with Lilith are shown as red, orange and yellow areas, and compared to the 68% and 95% CL contours from ATLAS (in blue). The best-fit point from Lilith is marked as a white star and the SM as a +.  [20], on the right for the H → τ τ channel from [21]. The 68% and 95% CL regions obtained with Lilith are shown in dark and light gray, respectively, and compared to the 68% and 95% CL contours from ATLAS (in blue). The best-fit points from Lilith and ATLAS are marked as white stars and blue dots, respectively.
H → W W * → 2l2ν (HIGG-2016-07): Ref. [20] focusses on the measurement of the inclusive ggH and VBF Higgs production cross sections in the H → W W * → eνµν channel. The paper quotes on page 13 signal strengths of µ(ggH, W W ) = 1.10 +0. 21 −0.20 and µ(VBF, W W ) = 0.62 +0. 36 −0.35 . We implemented these as a 2D result with a correlation of ρ = −0.08 using the variable Gaussian approximation; the correlation was fitted from the σ × BR plot, Fig. 9, of [20]. As no other validation material is available, we show in Fig. 3 (left) our reconstruction of the experimental likelihood in the µ(ggH, W W ) vs. µ(VBF, W W ) plane, comparing to the rescaled contours of Fig. 9 of the ATLAS paper.
H → τ τ (HIGG-2017-07): This ATLAS cross section measurement in the H → τ τ channel [21] provides as Aux. Fig. 5 the 68% and 95% CL contours in the µ(ggH, τ τ ) vs. µ(VBF, τ τ ) plane. A fit of a bivariate variable Gaussian to the 95% CL contour in this plot gives µ(ggH, τ τ ) 1.0 +0.72 −0.59 and µ(VBF, W W ) 1.20 +0.62 −0.56 with ρ −0.45, which are the values implemented in the database. As for H → W W above, no coupling fits are available which could be used for validation. We therefore show in Fig. 3 (right) our reconstruction of the experimental likelihood in the µ(ggH, τ τ ) vs. µ(VBF, τ τ ) plane. Note that a fit to the 68% CL contour of ATLAS gives a less good result.  [23] and µ(VBF, bb) = 3.0 +1.7 −1.6 [22]. No correlation data is available, so we implemented each of these as a 1D result; a Poisson likelihood is assumed per default but can easily be changed to a variable Gaussian if the user wishes to do so.
ttH production (HIGG-2017-02): The ATLAS paper [25], reporting evidence for ttH production, provides in Fig. 16 Fig. 4 the C F vs. C V fit from the implementation in Lilith to the official one from [25].
A few comments are in order here. First, the measurement of µ(ttH, γγ) given in [25] actually comes from [15] (HIGG-2016-21, see above) and is also included in the HIGG-2016-21 XML file; to avoid overlap when using both the HIGG-2016-21 and HIGG-2017-02 datasets, we provide a 3D XML file for the latter which includes only the V V , τ τ and bb, but not the γγ decay mode. The important point however is that the value given by ATLAS is not µ(ttH, γγ) but µ(ttH + tH, γγ). 7 This makes a big difference in the validation plot. Second, the individual measurement [26] gives µ(ttH, bb) to two decimals (0.84 +0.64 −0.61 ) instead just one (0.8 ± 0.6) in [25]. Again this makes a visible difference in Fig. 4, improving the quality of the fit, so we use the more precise numbers from [26]. The relevance of these points, and of the fitted correlations, is illustrated in Fig. 15 in Appendix A. Third, for µ(ttH, V V ) the contribution from H → W W * should dominate, but the concrete weights of the ZZ * and W W * decay modes are not given in [25]. (We use 95% for W W * and 5% for ZZ * as a rough estimate.) This is not a problem as long as C Z = C W ≡ C V , but one should not use the HIGG-2017-02 XML file for any other case.
H → invisible (HIGG-2016-28): Results from the search for invisibly decaying Higgs bosons produced in association with a Z boson are presented in [24]. A 95% CL upper limit of BR(H → inv.) < 0.67 is set for m H = 125 GeV assuming the SM ZH production cross section. In the Lilith database, we use a likelihood grid as function BR(H → inv.) extracted from Aux. Fig. 1c on the analysis' webpage.

CMS Run 2 results for 36 fb −1
The CMS Run 2 results included in this release are summarised in Table 4 Table 4: Overview of CMS Run 2 results included in this release. Note that we use the full 24 × 24 correlation matrix for the signal strengths for each production and decay mode combination provided in [12].
Combined measurements (HIG-17-031): CMS presented in [12] a combination of the individual measurements for the H → γγ [29], ZZ [30], W W [31], τ τ [32], bb [33,34] and µµ [35] decay modes as well as the ttH analyses [36][37][38]. We use the best fit values and uncertainties for the signal strengths for each production and decay mode combination presented in Table 3 of [12] together with the 24 × 24 correlation matrix provided as "Additional Figure 1" on the analysis webpage. Implemented as a variable Gaussian likelihood, this allows to reproduce well the coupling fits of the CMS paper as shown in Figs. 5 and 6.
V H, H → τ τ (HIG-18-007): The above data from [12] is supplemented by the results for the τ τ decay mode from the W H and ZH targeted analysis [28]. These are implemented in the form of 1D intervals for µ(ZH, τ τ ) and µ(W H, τ τ ) taken from Fig. 6 of [28].
H → invisible (HIG-17-023): In [27], CMS performed a search for invisible decays of a Higgs boson produced through vector boson fusion. We use the profile likelihood ratios for the qqH-, Z(ll)H-, V(qq')H-and ggH-tag categories extracted from their Fig. 8b together with the relative contributions from the different Higgs production mechanisms given in Table 6 of that paper. This assumes that the relative signal contributions stay roughly the same as for SM production cross sections. For validation, we reproduce in Fig. 7 the C g vs. C γ fit of [12], where the branching ratios of invisible and undetected decays are treated as free parameters.  Figure 5: Fit of C F vs. C V using best-fit values and uncertainties for the signal strengths for each production (ggH, VBF, WH, ZH, ttH) and decay (γγ, ZZ, W W , τ τ , bb, µµ) mode combination together with the 24 × 24 correlation matrix from the CMS combination paper [12]. The 1σ, 2σ and 3σ regions obtained with Lilith are shown as red, orange and yellow areas, and compared to the 1σ and 2σ contours from CMS (blue dots). The best-fit point from Lilith (CMS) is marked as a white star (blue dot), and the SM as a +.  Figure 6: Fit of tan β vs. cos(β − α) for the 2HDMs of Type I (left) and Type II (right) using the data from the combined CMS measurement [12]. The beige, orange and red filled areas show the 68%, 95% and 99.7% CL regions obtained with Lilith, while the blue dots mark the 95% CL contours from CMS.  Figure 7: Fit of C g vs. C γ using the data from the combined CMS measurement [12] and the search for invisible decays of a Higgs boson [27]. The branching ratios of invisible and undetected decays are treated as free parameters in the fit. The 1σ, 2σ and 3σ regions obtained with Lilith are shown as red, orange and yellow areas, and compared to the 1σ and 2σ contours from CMS (in blue). The best-fit point from Lilith (CMS) is marked as a white star (blue dot), and the SM as a +.

Status of Higgs coupling fits
In this section we give a brief overview of the current status of Higgs coupling fits. In all that follows, we use m H = 125.09 GeV. We begin by showing in Fig. 8 fits of C F vs. C V (left panel) and C g vs. C γ (right panel) using either the ATLAS (in blue) or the CMS (in green) Run 2 results in the current Lilith database, DB 19.06. As can be seen, the two experiments agree at the level of about 1σ, the ATLAS results being slightly closer to the SM (marked as a black + in all plots).
The situation when combining the results from both experiments is shown in Fig. 9. Using the Run 2 (Run 2 + Run 1) results of DB 19.06, 9 we find with the help of Minuit with a correlation of 0.31 (0.31). This assumes that contributions from new particles to the loop-induced couplings to gluons and photons as well as invisible or undetected decays are absent. Comparing the SM to the (C F , C V ) best fit gives −2 log(L SM /L C F ,C V max ) = 4.25 (5.51), corresponding to a p-value of 0.24 (0.14) based on Run 2 (Run 2 + Run 1) results.
Taking instead C g and C γ as free parameters with C F = C V = 1 (still assuming that invisible or undetected decays are absent), gives  It is also interesting to consider the couplings to up-type and down-type fermions as independent parameters. In this case, we find where we fitted separately for the two possible solutions of same-sign or opposite-sign C D with respect to C U , C V > 0. With −2 log L C D >0 max = 38.57 (43.40) compared to −2 log L C D <0 max = 39.07 (44.03), neither solution is clearly preferred by the data. Contours of constant CL in the C D vs. C U plane with C V profiled over can be seen in Fig. 10.
The (C F , C V ) and (C U , C D , C V ) fits above have their correspondence in the 2HDM of Type I and Type II, albeit C V is restricted to C V ≤ 1 in 2HDMs (and generally in models with only Higgs doublets and singlets). The couplings of the lighter scalar h are C F = cos α/ sin β in Type I, and C U = cos α/ sin β and C D = − sin α/ cos β in Type II; C V = sin(β − α) in both models. The fit results in the tan β vs. cos(β − α) plane are shown in Fig. 11. Note that for Type II the banana-shaped second branch corresponds to the "opposite-sign" solution for the bottom Yukawa coupling [39].
Before concluding, let us turn to invisible Higgs decays. Figure 12 (left) shows the 1D profile likelihood of BR(H → inv.) for two cases, SM couplings (in red) and C F and C V as free parameters (in blue). We see that the Run 2 results in DB 19.06 constrain BR(H → inv.) 5% at 95% CL for the SM-like case, and to BR(H → inv.) 17% when C F and C V are treated as free parameters; the case of free C U , C D , C V is not shown but gives almost the same result: BR(H → inv.) 18% at 95% CL. For Run 2+Run 1 results, these values tighten to BR(H → inv.) 4%, 15%, 16% for SM couplings, free C F , C V , and free C U , C D , C V , respectively. For completeness, Fig. 12 (right) shows the 68%, 95% and 99.7% CL regions from a 2D fit of BR(H → inv.) vs. C V with C F profiled over.

Conclusions
We presented Lilith-2.0, a light and easy-to-use Python tool for constraining new physics from signal strength measurements of the 125 GeV Higgs boson. The main novelties include • a better treatment of asymmetric uncertainties through the use of variable Gaussian or Poisson likelihoods where appropriate; • the use of multi-dimensional correlations whenever available; • a new database (DB 19.06) including the published ATLAS and CMS Run 2 Higgs results for 36 fb −1 .
We provided detailed validations of the results included in DB 19.06 and discussed the consequences of the available Run 2 results for fits of reduced Higgs couplings, 2HDMs of Type I and Type II, and invisible Higgs decays. Our analysis shows that the ATLAS and CMS results well agree with each other. The data is perfectly compatible with the SM, putting very tight constraints on any deviations. Indeed, our combination of the ATLAS and CMS results in global fits of (C F , C V ), (C g , C γ ) or (C U , C D , C V ) leads to a determination of these couplings to better than 10%. In particular, the uncertainty on C V shrinks to about 3-4%, and we observe a slight preference for C V > 1. In the context of 2HDMs, where C V ≤ 1, this forces one even deeper into the alignment limit [40,41]. Finally, the global fit also tightly constrains invisible Higgs decays -for SM-like couplings, to BR(H → inv.) 5% at 95% CL. Lilith-2.0 with its latest database is publicly available at http://lpsc.in2p3.fr/projects-th/lilith/ or directly on GitHub at https://github.com/sabinekraml/Lilith-2 and ready to be used to constrain a wide class of new physics scenarios. Lilith is also interfaced from micrOMEGAs [42] (v4.3 or higher). Readers who already have Lilith-1.1 in their micrOMEGAs installation can simply replace it with the new version 2.0. Given the high interest and ease of use of the signal strength framework, we kindly ask the experimental collaborations to continue to provide detailed signal strength results for the pure Higgs production×decay modes, including their correlations. The CMS combination paper [12] is an example of good practice in this respect. We want to stress, however, that when results are given for a combination of production and/or decay modes, it is important that the relative contributions be given and all assumptions be clearly spelled out. This is currently not the case in all publications. Moreover, as we have shown, it is crucial for a good usage of the experimental results, that the numbers quoted in tables and plots be precise enough. Such issues made the validation of some of the results in DB 19.06 seriously difficult, and we dearly hope that this will improve in the future.
Last but not least, we want to note that extracting results by digitizing curves from a plot, or typing the numbers for a large correlation matrix which is only available as an image, is painful, prone to errors, and should not be necessary in the modern information age. We therefore implore that all results be made available in digitized form, be it on HEPData or on the collaboration twiki page, by the time an analysis is published.

A Comparison with alternative implementations of the experimental results
To illustrate the importance of various improvements discussed throughout the paper, we here show versions of validation plots for alternative implementations of the experimental results. These can be compared with the corresponding validation plots for the final DB 19.06 release in Section 5.
To start with, we show in Fig. 13 fits of C F vs. C V , where we used the ordinary Gaussian approximation (i.e. with fixed width) instead of the variable Gaussian. The plot on the left is for the ATLAS H → γγ [15] data, which has a 4 × 4 correlation matrix for the ggH, VBF, VH and ttH+tH production modes. The plot on the right is for the CMS combination [12], which has a 24 × 24 correlation matrix. It is obvious that the fixed-width Gaussian does not correctly approximate the true likelihood. The variable Gaussian implementation, on the other hand, gives a satisfying result.
In some cases with large asymmetries, a Poisson distribution is more appropriate than a variable Gaussian. This is the case for the ATLAS H → ZZ * → 4l analysis [19] as shown in Fig. 14. Note, however, that although the Poisson form allows to better reproduce the 68% CL contour of ATLAS than the variable Gaussian, the 95% CL contour is still quite off. This is much improved by using the 1D profile likelihoods for µ(ggH, ZZ * ) and µ(VBF, ZZ * ) from Aux. Figs. 7a,b of [19] (also parametrised as Poisson likelihoods) instead of the best-fit and uncertainty values from Table 9 of [19], see Fig. 2.
Finally, Fig. 15 details the steps we made to achieve a good implementation and validation of the ATLAS ttH combination. The labels "best fits & uncertainties as given in HIGG-2017-02 (v1)" and "... (v2)" refer to identifying the measurement in the γγ final state as µ(ttH, γγ) or µ(ttH + tH, γγ). See also the discussion related to Fig. 4 Table 9 of [19] with correlation ρ = −0.41 from Aux. Fig. 4c on the analysis twiki page. On the left, the likelihood of µ(ggH, ZZ * ) vs. µ(VBF, ZZ * ) is parametrised as a 2D variable Gaussian, on the right as a 2D Poisson distribution. To be compared with Fig. 2 [19], are used to construct a 2D Poisson likelihood. VH and ttH production are treated as in Fig. 2. for Y = τ τ , γγ, bb, V V as quoted in Fig. 16 of ATLAS-HIGG-2017-02 [25]. Top right: same as on the left but including tH production for the γγ channel according to the second bullet point on p. 37 of [25]. Bottom left: adding the correlations fitted from Figs. 17a,b of [25]. Bottom right: using µ(ttH, bb) from [26] instead of the value from [25]. The bottom right plot is already very close to the final version shown in Fig. 4. The only difference is that for the latter µ(ttH + tH, γγ) fitted from the 1D profile likelihood in [15] was used.