Constraining new physics from Higgs measurements with
[1mm] Lilith: update to LHC Run 2 results
Sabine Kraml^{1*}, Tran Quang Loc^{2,3}, Dao Thi Nhung^{2}, Le Duc Ninh^{2}
1 Laboratoire de Physique Subatomique et de Cosmologie, Université GrenobleAlpes,
CNRS/IN2P3, 53 Avenue des Martyrs, F38026 Grenoble, France
2 Institute For Interdisciplinary Research in Science and Education, ICISE,
590000 Quy Nhon, Vietnam
3 VNUHCMUniversity of Science, 227 Nguyen Van Cu, Dist. 5,
Ho Chi Minh City, Vietnam
[2mm] *
August 7, 2021
Abstract
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 Lilith1.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, TwoHiggsdoublet models of Type I and Type II, and invisible Higgs decays. Lilith2.0 is available on GitHub and ready to be used to constrain a wide class of new physics scenarios.
1 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 nonSM Higgs sector with the ATLAS and CMS results requires to construct a likelihood, which is a nontrivial 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 indepth 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}^{1}1For a discussion of the use and usability of signal strength results, and recommendations on their presentation, see [3]. 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.
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 lightweight Python library that uses as a primary input signal strength results
(1) 
in which the fundamental production and decay modes are unfolded from experimental categories. Here, the main production mechanisms are: gluon fusion (), vectorboson fusion (), associated production with an electroweak gauge boson ( and , collectively denoted as ) and associated production with top quarks, mainly but also . The main decay modes accessible at the LHC are , , , and (with ).^{2}^{2}2When is not directly available, is used, introducing appropriate efficiency factors . For inclusive combinations of production channels for the same , the efficiencies become .
The signal strength framework is based on the narrowwidth approximation and on the assumption that new physics results only in the scaling of SM Higgs processes.^{3}^{3}3In other words, the Lagrangian has the same tensor structure as the SM, see e.g. the discussion in [3]. 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 and the gluoninitiated 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 and that scale the amplitudes for the production and decay of the SM Higgs boson, respectively, as
(2) 
for the different production modes
(3) 
where and () 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 treelevel couplings, we define the loopinduced couplings and of the Higgs to gluons and photons, respectively. If no new particles appear in the loops, and are computed from the couplings in Eq. (3) following the procedure established in [9]. Alternatively, and can be taken as free parameters. Apart from the different notation, this is equivalent to the socalled “kappa framework” of [9]. Finally note that often a subset of the ’s in Eq. (3) is taken as universal, in particular (custodial symmetry), and like in the TwoHiggsdoublet model (2HDM) of Type II, or 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
(4) 
without spoiling the approximation.
The rest of the paper is organised as follows. We begin the discussion of novelties in Lilith2.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 Lilith2.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 Lilith2.0.
Instead, we present only what is new with respect to Lilith1.1.
For everything else, including instructions how to use the code, we refer the reader to the
original manual [1].
2 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.
The root tag of each XML file is <expmu>, which has two mandatory attributes, dim and type to specify the type of signal strength result. Production and decay modes are specified via prod and decay attributes either directly in the <expmu> tag or in the efficiency <eff> tags. The latter option allows for the inclusion of different production and decay modes in one XML file. Additional (optional) information can be provided in <experiment>, <source>, <sqrts>, <CL> and <mass> tags. Taking the result from the combined ATLAS and CMS Run 1 analysis [10] as a concrete example, the structure is
<expmu decay="gammagamma" dim="2" type="n"> <experiment>ATLASCMS</experiment> <source type="publication">CMSHIG15002; ATLASHIGG201507</source> <sqrts>7+8</sqrts> <mass>125.09</mass> <CL>68%</CL> <eff axis="x" prod="ggH">1.</eff> <eff axis="y" prod="VVH">1.</eff> <! (...) > </expmu>
where <! (...) > is a placeholder for the actual likelihood information. For a detailed description, we refer to the original Lilith manual [1]. In the following, we assume that the reader is familiar with the basic syntax.
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 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 . 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 onesided 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 multidimensional 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.
1D likelihood parametrisation
For 1D data, the format remains the same as in [1]. For example, a signal strength is implemented as
<eff prod="ZH">1.</eff> <bestfit>1.12</bestfit> <param> <uncertainty side="left">0.45</uncertainty> <uncertainty side="right">0.50</uncertainty> </param>
The <bestfit> tag contains the bestfit value, while the <uncertainty> tag contains the left (negative) and right (positive) errors.^{4}^{4}4The values in the <uncertainty> tag can be given with or without a sign. The choice of likelihood function is done by setting type="n" for an ordinary, 2sided Gaussian (as in Lilith1.1); type="vn" for a variable Gaussian; or type="p" for a Poisson form in the <expmu> tag.
2D likelihood parametrisation
For type="vn" and type="p", signal strengths in 2D with a correlation are now described in an analogous way as 1D data. For example, and with a correlation of can be implemented as
<expmu decay="WW" dim="2" type="vn"> <eff axis="x" prod="ggH">1.0</eff> <eff axis="y" prod="VBF">1.0</eff> <bestfit> <x>1.10</x> <y>0.62</y> </bestfit> <param> <uncertainty axis="x" side="left">0.20</uncertainty> <uncertainty axis="x" side="right">+0.21</uncertainty> <uncertainty axis="y" side="left">0.35</uncertainty> <uncertainty axis="y" side="right">+0.36</uncertainty> <correlation>0.08</correlation> </param> </expmu>
Here, the <eff> tag is used to declare the x and y axes, specified by their production and/or decay channels together with the corresponding efficiencies. The <bestfit> tag specifies the location of the bestfit point in the (x,y) plane. The <uncertainty> tags contain the left (negative) and right (positive) errors for the x and y axes, and finally the <correlation> tag specifies the correlation between x and y. The choice of likelihood function is again done by setting type="vn" or type="p" in the <expmu> tag.
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].
Multidimensional data
For correlated signal strengths in more than 2 dimensions, a new format is introduced. We here illustrate it by means of the CMS result [12], which has signal strengths for 24 production and decay mode combinations plus a correlation matrix. First, we set dim="24" and label the various signal strengths as axes d1, d2, d3, … d24:^{5}^{5}5The <experiment>, <source>, <sqrts>, etc. tags are omitted for brevity.
<expmu dim="24" type="vn"> <eff axis="d1" prod="ggH" decay="gammagamma">1.0</eff> <eff axis="d2" prod="ggH" decay="ZZ">1.0</eff> <eff axis="d3" prod="ggH" decay="WW">1.0</eff> ... <eff axis="d24" prod="ttH" decay="tautau">1.0</eff>
The bestfit values for each axis are specified as
<bestfit> <d1>1.16</d1> <d2>1.22</d2> <d3>1.35</d3> ... <d24>0.23</d24> </bestfit>
The <param> tag then contains the uncertainties and correlations in the form
<param> <uncertainty axis="d1" side="left">0.18</uncertainty> <uncertainty axis="d1" side="right">+0.21</uncertainty> <uncertainty axis="d2" side="left">0.21</uncertainty> <uncertainty axis="d2" side="right">+0.23</uncertainty> ... <uncertainty axis="d24" side="left">0.88</uncertainty> <uncertainty axis="d24" side="right">+1.03</uncertainty> <correlation entry="d1d2">0.12</correlation> <correlation entry="d1d3">0.16</correlation> <correlation entry="d1d4">0.08</correlation> ... <correlation entry="d23d24">0</correlation> </param> </expmu>
This will also work for type="n", see Eq. (9) in the next section.
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 Lilith1.1. This may change in future versions.
3 Likelihood calculation
The statistical procedure used in Lilith was described in detail in [1]. The main quantity given as an output is the loglikelihood, , evaluated in computelikelihood.py from the information given in the XML database. In this section, we explain how 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].
3.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
(5) 
where denotes the bestfit signal strength, and and are absolute values of the left and right uncertainties at % 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 can also be problematic when , in which case it may be better to use the ordinary 2sided Gaussian implementation, in particular for 1D data; this has to be decided casebycase upon validation. In any case, Lilith issues an error message whenever a numerical problem occurs.
In the case of dimensional data (), we use the correlations given by the experimental collaboration, if available, together with the best fit points and the left and right uncertainties at % CL. When results are given in terms of 2D contour plots, we can also use the variable Gaussian form to numerically determine the bestfit point, uncertainties and their correlation, if not given explicitly by the experimental collaboration. For the dimensional signal strength vector , the likelihood reads
(6) 
where the best fit point and the covariance matrix is constructed from the correlation matrix as
(7) 
with
(8) 
Here the and are the left and right uncertainties at % CL of the th combination of production and/or decay channels, respectively. For multidimensional data in the ordinary Gaussian approximation, the relation between covariance matrix and the correlation matrix becomes
(9) 
where and .
3.2 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
(10) 
where and are determined numerically from the equations
(11) 
More concretely, is determined from the expression on the left by bifurcation between and 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 follows the marginal distribution, while data of the channel belongs to the conditional distribution. The joint loglikelihood is then
(12) 
with
(13) 
and
(14) 
The function reads
(15) 
where is solved numerically from the correlation expression
(16) 
and the and are solutions of Eq. (11) for the and channels, respectively.
4 New production channels
As mentioned in the introduction, we have also included a few new production channels. These are ZH production via gluongluon 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 channel (qqZH). However, the loopinduced gluongluon fusion is not so small, about of the total cross section at 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 and initiated processes, with relative weights . In terms of the scale factors of Eq. (2), this gives
(17) 
For the SM cross sections, the values given in [6, 14, 15] are used, where the qqZH cross section is calculated at the NexttoNexttoLeading Order (NNLO) QCD + NexttoLeading Order (NLO) electroweak level, and the ggZH cross section is calculated at the NLO QCD level. The definitions of VH (ZH and WH) and VVH (ZH, WH and VBF) follow straightforwardly. The WH cross section is calculated at NNLO QCD and NLO electroweak accuracies, while the VBF cross section is calculated at approximate NNLO QCD and NLO electroweak accuracies. Further details are provided in [6, 14, 15].
The production also includes two contributions: channel production and production. The channel cross section is much smaller and hence not included. Interference effects between these channels are also neglected. At TeV, is dominant, contributing about of the cross section. As above, tHq and tHW are combined to production with efficiencies calculated from the SM cross sections. Moreover, following the usage in the ATLAS analysis [16], a combination of tH and ttH named ‘top’ production is defined in an analogous way.^{6}^{6}6Accordingly, prod="tHq", "tHW", "tH" and "top" can be used in the XML files in the database. Thus
(18)  
(19) 
The value for is calculated at NLO QCD and NLO electroweak accuracies, while is calculated at NLO QCD level as given in [6, 14, 15]. has been calculated at leading order using MadGraph [17] with factorization and renormalization scales and the NNPDF30_lo_as_0130 PDF set [18]. Other input parameters are the same as in [6] section I.6 for production. It is noted that the definition of an NLO cross section for the channel is not straightforward because of the interferences with the ttH channel, see [19] 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, analogous 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],
(20)  
(21)  
(22)  
(23) 
where the coefficients , , and , being the SM predictions, are provided in Tables 1 and 2. Note that for the time being fixed values for GeV are used. Moreover, for the case TeV, the values at TeV are used; the differences are negligible in the current approximations.
[TeV]  

[TeV]  

Finally, we note that the accuracies of the remaining SM predictions, including ggH cross sections and the Higgs branching fractions, are unchanged compared to the previous version of Lilith; in particular (N)NLO QCD corrections are included as explained in [1].
5 ATLAS and CMS results included in the database update
In this section, we discuss the ATLAS and CMS Run 2 results included in the Lilith database. Most of this is based on the database release DB 19.06 (June 2019). Three ATLAS results (HIGG201714, HIGG201610 and HIGG201854) were added to this during the peerreview process, so that the current latest version is DB 19.09 (September 2019). There is no other difference between DB 19.06 and DB 19.09. The validation plots in this section carry their original DB version number, that is 19.06 or 19.09, while the fit results in Section 6 are all for the complete DB 19.09.
5.1 ATLAS Run 2 results for 36 fb
The ATLAS Run 2 results included in this release are summarised in Table 3 and explained in more detail below.
mode  inv.  

ggH  [16]  [21]  [22]  [23]  –  [24]  – 
VBF  [16]  [21]  [22]  [23]  [25]  [24]  [26] 
WH  [16]  [21]  [27]  –  [28]  [24]  – 
ZH  [27]  –  [28]  [29]  
ttH  [16]  [21, 30]  [30]  [30]  [30, 31]  –  – 
(HIGG201621):
The ATLAS analysis [16] provides in Fig. 12 signal strengths for 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 stage0 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 or .
We have therefore extracted the bestfit points and uncertainties
from the 1D profile likelihoods, which are provided as Auxiliary Figures 23ad on the analysis webpage,
as^{7}^{7}7In the XML file, we use the exact numbers from the fit to the 1D profile likelihoods.
,
,
and
(using a Poisson likelihood).
These numbers are consistent with the rounded values in Fig. 12 of [16], but using more digits
improves the coupling fits as shown in Fig. 1.
(HIGG201622): A similar issue as discussed for above arises
for . In order to reasonably reproduce the vs. fit of ATLAS (Fig. 8b of [21]),
we fit the 1D profile likelihoods for and shown in Aux. Figs. 7a and 7b
of [21] as Poisson distributions. This gives and
, which we implement as a bivariate Poisson distribution with
correlation (from Aux. Fig. 4c of [21]).
For the VH and ttH production modes, lacking more information, we convert the given 95% CL limits into
and using a 2sided Gaussian
(assuming 1sided limits gives a less good agreement with the ATLAS vs. 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.
(HIGG201607 and HIGG201714): Ref. [22] focusses on the measurement of the
inclusive ggH and VBF Higgs production cross sections in the channel. The paper quotes
on page 13 signal strengths of and .
We implemented these as a 2D result with a correlation of using the variable Gaussian approximation;
the correlation was fitted from the plot, Fig. 9, of [22].
In addition, Ref. [27] presents the measurement of the () channel for Higgs boson production in association with a vector boson. Using again the variable Gaussian approximation, we extracted
and
with correlation from Fig. 8 of that paper.
As no other validation material is available,
we show in Fig. 3 (top left and top right plots) our reconstruction of the experimental likelihoods in the
vs. and vs. planes,
comparing respectively to the rescaled contours of Fig. 9 of [22] and Fig. 8 of [27].^{8}^{8}8Here and in the following, whenever the experimental paper gives only measured cross sections but no signal strengths, we use the recommended SM cross section values from the LHC Higgs Cross Section Working Group [15] for appropriate rescaling.
(HIGG201707): This ATLAS cross section measurement in the channel [23]
provides as Aux. Fig. 5 the 68% and 95% CL contours in the vs. plane.
A fit of a bivariate variable Gaussian to the 95% CL contour in this plot
gives and with
, which are the values implemented in the database.
As for above, no coupling fits are available which could be used for validation.
We therefore show in Fig. 3 (bottom plot) our reconstruction of the experimental likelihood in the
vs. plane.
Note that a fit to the 68% CL contour of ATLAS gives a less good result.
(HIGG201610): In Ref. [24], ATLAS reports a measured overall signal strength
of , from which a 95% CL limit of is computed using the prescription.
For consistency with the 95% CL limit, that is to avoid being overconstraining,
we implement this as in the Lilith database.
The relative contributions of ggH (90%), VBF (7%) and VH (3%) production are estimated from Aux. Table 3
on the analysis’ webpage, using the sum of all eight orthogonal categories.
(HIGG201629 and HIGG201630): For the decay mode, ATLAS gives
, [28]
and [25].
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 (HIGG201702): The ATLAS paper [30], reporting evidence for production, provides in Fig. 16 the signal strength results broken down into , , and decay modes from a combined analysis of all searches. Correlations are not given explicitly but can be estimated from Figs. 17a and 17b in [30] as for the correlation between the and decay modes and for that between the and decay modes. For validation, we compare in Fig. 4 the vs. fit from the implementation in Lilith to the official one from [30].
A few comments are in order here. First, the measurement of given in [30]
actually comes from [16] (HIGG201621, see above) and is also included in the HIGG201621 XML file;
to avoid overlap when using both the HIGG201621 and HIGG201702 datasets, we provide a 3D XML file for the latter
which includes only the , and , but not the decay mode.
The important point however is that the value given by ATLAS is not but .^{9}^{9}9Lacking more precise information, ttH and tH production are combined according to Eq. (19).
This makes a big difference in the validation plot.
Second, the individual measurement [31] gives to two decimals
() instead just one () in [30]. Again this makes a visible difference
in Fig. 4, improving the quality of the fit, so we use the more precise numbers from [31].
The relevance of these points, and of the fitted correlations, is illustrated in Fig. 15
in Appendix A.
Third, for the contribution from should dominate, but the concrete weights of the
and decay modes are not given in [30]. (We use 95% for and 5% for as a rough estimate.)
This is not a problem as long as , but one should
not use the HIGG201702 XML file for any other case.
invisible (HIGG201628 and HIGG201854): Results from the search for invisibly decaying Higgs bosons produced in association with a boson are presented in [29]. A 95% CL upper limit of BR is set for GeV assuming the SM production cross section. In the Lilith database, we use a likelihood grid as function BR extracted from Aux. Fig. 1c on the analysis’ webpage. The combination of Run 2 searches for invisible Higgs boson decays in [26] tightens BR at 95% CL; here, we use the likelihood grid for VBF production of invisibly decaying Higgs bosons extracted from Fig. 1a.
5.2 CMS Run 2 results for 36 fb
The CMS Run 2 results included in this release are summarised in Table 4 and explained in more detail below.
mode  inv.  

ggH  [12]  [12]  [12]  [12]  [12]  [12]  [32] 
VBF  [12]  [12]  [12]  [12]  –  [12]  [32] 
WH  [12]  [12]  [12]  [33]  [12]  –  [32] 
ZH  [12]  [12]  [12]  [33]  [12]  –  [32] 
ttH  [12]  [12]  [12]  [12]  [12]  –  – 
Combined measurements (HIG17031):
CMS presented in [12] a combination of the individual measurements for the
[34], [35], [36],
[37], [38, 39] and [40]
decay modes as well as the analyses [41, 42, 43].
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 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.
, (HIG18007): The above data from [12] is supplemented by the results
for the decay mode from the and targeted analysis [33]. These are implemented in the
form of 1D intervals for and taken from Fig. 6 of [33].
invisible (HIG17023):
In [32], CMS performed a search for invisible decays of a Higgs boson produced through vector boson fusion,
setting a limit of BR at 95% CL.
We use the profile likelihood ratios for the qqH, Z(ll)H, V(qq’)H and ggHtag 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 vs. fit of [12], where the branching ratios of invisible and undetected decays
are treated as free parameters.^{10}^{10}10The profiling in Fig. 7 was done with Minuit.
Since Minuit does not allow conditional limits, in this case ,
we demanded that both BR and BR be less than 50%.
6 Status of Higgs coupling fits
In this section we give a brief overview of the current status of Higgs coupling fits. We remind the reader that these fits rely on the specific assumptions mentioned in the Introduction and do not represent general, modelindependent statements on “the Higgs couplings”. Most importantly, new physics effects are assumed to follow the same Lorentz structure as the SM Higgs couplings, such that new physics contributions to the Higgs production and decay processes can be parameterised in terms of reduced coupling factors and some coefficients depending only on SM predictions, as seen e.g. from Eqs. (20)–(23) in Section 4. Details about the accuracy of the SM predictions are also provided in Section 4. In all that follows, we use GeV. For other global fits to Run 2 Higgs results see, e.g., [44, 45, 46, 47].
We begin by showing in Fig. 8 fits of vs. (left panel) and vs. (right panel) using either the ATLAS (in blue) or the CMS (in green) Run 2 results in the current Lilith database, DB 19.09. As can be seen, the two experiments agree at the level of about , 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.09,^{11}^{11}11For Run 1, we use the results from the official ATLAS+CMS combination [10] available in DB 19.09, plus the individual ATLAS and CMS results. We also note here that for the SM, we get using the Run 2 results in DB 19.09 (53 measurements) and from the combination of Run 2 and Run 1 results (66 measurements). we find with the help of Minuit
(24) 
with a correlation of (). This assumes that contributions from new particles to the loopinduced couplings to gluons and photons as well as invisible or undetected decays are absent. Comparing the SM to the best fit gives (), corresponding to a value of () based on Run 2 (Run 2 + Run 1) results.
Taking instead and as free parameters with (still assuming that invisible or undetected decays are absent), gives
(25) 
with correlation () from Run 2 (combining Run 2 and Run 1) results; here the SM point has a value of ().
It is also interesting to consider the couplings to uptype and downtype fermions as independent parameters. In this case, we find

(26) 
where we fitted separately for the two possible solutions of samesign or oppositesign with respect to . With () compared to (), neither solution is clearly preferred by the data. Contours of constant CL in the vs. plane with profiled over can be seen in Fig. 10.
The () and () fits above have their correspondence in the 2HDM of Type I and Type II, albeit is restricted to in 2HDMs (and generally in models with only Higgs doublets and singlets). The couplings of the lighter scalar are in Type I, and and in Type II; in both models. The fit results in the vs. plane are shown in Fig. 11. Note that for Type II the bananashaped second branch corresponds to the “oppositesign” solution for the bottom Yukawa coupling [48].
Before concluding, let us turn to invisible Higgs decays. Figure 12 (left) shows the 1D profile likelihood of for two cases, SM couplings (in red) and and as free parameters (in blue). We find that the Run 2 results in DB 19.09 constrain at 95% CL for the SMlike case,^{12}^{12}12This strong bound is in fact primarily driven by the signal strength measurements in the SM final states, as invisible decays reduce the branching ratios to (visible) SM particles [49]. and to when and are treated as free parameters; the case of free , , is not shown but gives the same result. For Run 2+Run 1 results, these values tighten to , 15%, 15% for SM couplings, free , , and free , , , respectively. For completeness, Fig. 12 (right) shows the 68%, 95% and 99.7% CL regions from a 2D fit of BR() vs. with profiled over.
7 Conclusions
We presented Lilith2.0, a light and easytouse 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 multidimensional correlations whenever available;

a new database (DB 19.09) including the published ATLAS and CMS Run 2 Higgs results for 36 fb.
We provided detailed validations of the results included in DB 19.09 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 (, ), (, ) or (, , ) leads to a determination of these couplings to better than 10%. In particular, the uncertainty on shrinks to about 3–4%, and we observe a slight preference for . In the context of 2HDMs, where , this forces one even deeper into the alignment limit [50, 51]. Finally, the global fit also tightly constrains invisible Higgs decays — for SMlike couplings, to at 95% CL.
Lilith2.0 with its latest database is publicly available at
or directly on GitHub at https://github.com/sabinekraml/Lilith2 and ready to be used to constrain a wide class of new physics scenarios. Lilith is also interfaced from micrOMEGAs [52] (v4.3 or higher). Readers who already have Lilith1.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 productiondecay 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.09 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.
Acknowledgements
S.K. thanks R. Schöfbeck, W. Waltenberger and N. Wardle for helpful discussions. This work was supported in part by the IN2P3 theory project “LHCitools: methods and tools for the interpretation of the LHC Run 2 results for new physics”. D.T.N and L.D.N. are funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.012017.78. D.T.N. thanks the LPSC Grenoble for hospitality and financial support for a research visit within the LHCitools project. T.Q.L. acknowledges the hospitality and financial support of the ICISE Quy Nhon.
Appendix 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 official release in Section 5.
To start with, we show in Fig. 13 fits of vs. , 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 [16] data, which has a 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 correlation matrix. It is obvious that the fixedwidth 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 analysis [21] 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 and from Aux. Figs. 7a,b of [21] (also parametrised as Poisson likelihoods) instead of the bestfit and uncertainty values from Table 9 of [21], 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 HIGG201702 (v1)” and “… (v2)” refer to identifying the measurement in the final state as or . See also the discussion related to Fig. 4 in Section 5.1.
References
 [1] J. Bernon and B. Dumont, Lilith: a tool for constraining new physics from Higgs measurements, Eur. Phys. J. C75(9), 440 (2015), doi:10.1140/epjc/s1005201536459, arXiv:1502.04138.
 [2] J. Bernon, B. Dumont and S. Kraml, Status of Higgs couplings after run 1 of the LHC, Phys. Rev. D90, 071301 (2014), doi:10.1103/PhysRevD.90.071301, arXiv:1409.1588.
 [3] F. Boudjema et al., On the presentation of the LHC Higgs Results, In Workshop on Likelihoods for the LHC Searches Geneva, Switzerland, January 2123, 2013 (2013), arXiv:1307.5865.
 [4] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, : Confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC, Eur. Phys. J. C74(2), 2711 (2014), doi:10.1140/epjc/s1005201327114, https://higgsbounds.hepforge.org/, arXiv:1305.1933.
 [5] J. R. Andersen et al., Les Houches 2015: Physics at TeV Colliders Standard Model Working Group Report, In 9th Les Houches Workshop on Physics at TeV Colliders (PhysTeV 2015) Les Houches, France, June 119, 2015 (2016), arXiv:1605.04692.
 [6] D. de Florian et al., Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector (2016), doi:10.23731/CYRM2017002, arXiv:1610.07922.
 [7] G. Belanger, B. Dumont, U. Ellwanger, J. F. Gunion and S. Kraml, Global fit to Higgs signal strengths and couplings and implications for extended Higgs sectors, Phys. Rev. D88, 075008 (2013), doi:10.1103/PhysRevD.88.075008, arXiv:1306.2941.
 [8] G. Belanger, B. Dumont, U. Ellwanger, J. F. Gunion and S. Kraml, Higgs Couplings at the End of 2012, JHEP 02, 053 (2013), doi:10.1007/JHEP02(2013)053, arXiv:1212.5244.
 [9] J. R. Andersen et al., Handbook of LHC Higgs Cross Sections: 3. Higgs Properties (2013), doi:10.5170/CERN2013004, arXiv:1307.1347.
 [10] G. Aad et al., Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at and 8 TeV, JHEP 08, 045 (2016), doi:10.1007/JHEP08(2016)045, arXiv:1606.02266.
 [11] R. Barlow, Asymmetric statistical errors, In Statistical Problems in Particle Physics, Astrophysics and Cosmology (PHYSTAT 05): Proceedings, Oxford, UK, September 1215, 2005, pp. 56–59 (2004), arXiv:physics/0406120.
 [12] A. M. Sirunyan et al., Combined measurements of Higgs boson couplings in proton–proton collisions at TeV, Eur. Phys. J. C79(5), 421 (2019), doi:10.1140/epjc/s100520196909y, https://cmsresults.web.cern.ch/cmsresults/publicresults/publications/HIG17031/, arXiv:1809.10733.
 [13] P. Berkhout and E. Plug, A bivariate Poisson count data model using conditional probabilities, In Statistica Neerlandica, pp. 349–364 (2004).
 [14] https://twiki.cern.ch/twiki/bin/view/LHCPhysics/CERNYellowReportPageAt8TeV.
 [15] https://twiki.cern.ch/twiki/bin/view/LHCPhysics/CERNYellowReportPageAt13TeV.
 [16] M. Aaboud et al., Measurements of Higgs boson properties in the diphoton decay channel with 36 fb of collision data at TeV with the ATLAS detector, Phys. Rev. D98, 052005 (2018), doi:10.1103/PhysRevD.98.052005, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201621/, arXiv:1802.04146.
 [17] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, The automated computation of treelevel and nexttoleading order differential cross sections, and their matching to parton shower simulations, JHEP 07, 079 (2014), doi:10.1007/JHEP07(2014)079, arXiv:1405.0301.
 [18] R. D. Ball et al., Parton distributions for the LHC Run II, JHEP 04, 040 (2015), doi:10.1007/JHEP04(2015)040, arXiv:1410.8849.
 [19] F. Demartin, B. Maier, F. Maltoni, K. Mawatari and M. Zaro, tWH associated production at the LHC, Eur. Phys. J. C77(1), 34 (2017), doi:10.1140/epjc/s1005201746017, arXiv:1607.05862.
 [20] https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHXSWG2KAPPA.
 [21] M. Aaboud et al., Measurement of the Higgs boson coupling properties in the decay channel at = 13 TeV with the ATLAS detector, JHEP 03, 095 (2018), doi:10.1007/JHEP03(2018)095, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201622/, arXiv:1712.02304.
 [22] M. Aaboud et al., Measurements of gluongluon fusion and vectorboson fusion Higgs boson production crosssections in the decay channel in collisions at TeV with the ATLAS detector, Phys. Lett. B789, 508 (2019), doi:10.1016/j.physletb.2018.11.064, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201607, arXiv:1808.09054.
 [23] M. Aaboud et al., Crosssection measurements of the Higgs boson decaying into a pair of leptons in protonproton collisions at TeV with the ATLAS detector, Phys. Rev. D99, 072001 (2019), doi:10.1103/PhysRevD.99.072001, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201707, arXiv:1811.08856.
 [24] M. Aaboud et al., Search for the dimuon decay of the Higgs boson in collisions at TeV with the ATLAS detector, Phys. Rev. Lett. 119(5), 051802 (2017), doi:10.1103/PhysRevLett.119.051802, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201610/, arXiv:1705.04582.
 [25] M. Aaboud et al., Search for Higgs bosons produced via vectorboson fusion and decaying into bottom quark pairs in collisions with the ATLAS detector, Phys. Rev. D98(5), 052003 (2018), doi:10.1103/PhysRevD.98.052003, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201630, arXiv:1807.08639.
 [26] M. Aaboud et al., Combination of searches for invisible Higgs boson decays with the ATLAS experiment, Phys. Rev. Lett. 122(23), 231801 (2019), doi:10.1103/PhysRevLett.122.231801, arXiv:1904.05105.
 [27] G. Aad et al., Measurement of the production cross section for a Higgs boson in association with a vector boson in the channel in collisions at TeV with the ATLAS detector (2019), https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201714/, arXiv:1903.10052.
 [28] M. Aaboud et al., Evidence for the decay with the ATLAS detector, JHEP 12, 024 (2017), doi:10.1007/JHEP12(2017)024, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201629, arXiv:1708.03299.
 [29] M. Aaboud et al., Search for an invisibly decaying Higgs boson or dark matter candidates produced in association with a boson in collisions at TeV with the ATLAS detector, Phys. Lett. B776, 318 (2018), doi:10.1016/j.physletb.2017.11.049, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201628/, arXiv:1708.09624.
 [30] M. Aaboud et al., Evidence for the associated production of the Higgs boson and a top quark pair with the ATLAS detector, Phys. Rev. D97(7), 072003 (2018), doi:10.1103/PhysRevD.97.072003, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201702, arXiv:1712.08891.
 [31] M. Aaboud et al., Search for the standard model Higgs boson produced in association with top quarks and decaying into a pair in collisions at TeV with the ATLAS detector, Phys. Rev. D97(7), 072016 (2018), doi:10.1103/PhysRevD.97.072016, https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HIGG201703, arXiv:1712.08895.
 [32] A. M. Sirunyan et al., Search for invisible decays of a Higgs boson produced through vector boson fusion in protonproton collisions at TeV (2018), doi:10.1016/j.physletb.2019.04.025, https://cmsresults.web.cern.ch/cmsresults/publicresults/publications/HIG17023/, arXiv:1809.05937.
 [33] A. M. Sirunyan et al., Search for the associated production of the Higgs boson and a vector boson in protonproton collisions at TeV via Higgs boson decays to leptons, JHEP 06, 093 (2019), doi:10.1007/JHEP06(2019)093, https://cmsresults.web.cern.ch/cmsresults/publicresults/publications/HIG18007/, arXiv:1809.03590.
 [34] A. M. Sirunyan et al., Measurements of Higgs boson properties in the diphoton decay channel in protonproton collisions at 13 TeV, JHEP 11, 185 (2018), doi:10.1007/JHEP11(2018)185, arXiv:1804.02716.
 [35] A. M. Sirunyan et al., Measurements of properties of the Higgs boson decaying into the fourlepton final state in pp collisions at TeV, JHEP 11, 047 (2017), doi:10.1007/JHEP11(2017)047, arXiv:1706.09936.
 [36] A. M. Sirunyan et al., Measurements of properties of the Higgs boson decaying to a W boson pair in pp collisions at TeV, Phys. Lett. B791, 96 (2019), doi:10.1016/j.physletb.2018.12.073, arXiv:1806.05246.
 [37] A. M. Sirunyan et al., Observation of the Higgs boson decay to a pair of leptons with the CMS detector, Phys. Lett. B779, 283 (2018), doi:10.1016/j.physletb.2018.02.004, arXiv:1708.00373.
 [38] A. M. Sirunyan et al., Evidence for the Higgs boson decay to a bottom quark–antiquark pair, Phys. Lett. B780, 501 (2018), doi:10.1016/j.physletb.2018.02.050, arXiv:1709.07497.
 [39] A. M. Sirunyan et al., Inclusive search for a highly boosted Higgs boson decaying to a bottom quarkantiquark pair, Phys. Rev. Lett. 120(7), 071802 (2018), doi:10.1103/PhysRevLett.120.071802, arXiv:1709.05543.
 [40] A. M. Sirunyan et al., Search for the Higgs boson decaying to two muons in protonproton collisions at TeV, Phys. Rev. Lett. 122(2), 021801 (2019), doi:10.1103/PhysRevLett.122.021801, arXiv:1807.06325.
 [41] A. M. Sirunyan et al., Evidence for associated production of a Higgs boson with a top quark pair in final states with electrons, muons, and hadronically decaying leptons at 13 TeV, JHEP 08, 066 (2018), doi:10.1007/JHEP08(2018)066, arXiv:1803.05485.
 [42] A. M. Sirunyan et al., Search for production in the decay channel with leptonic decays in protonproton collisions at TeV, JHEP 03, 026 (2019), doi:10.1007/JHEP03(2019)026, arXiv:1804.03682.
 [43] A. M. Sirunyan et al., Search for H production in the alljet final state in protonproton collisions at TeV, JHEP 06, 101 (2018), doi:10.1007/JHEP06(2018)101, arXiv:1803.06986.
 [44] V. Sanz and J. Setford, Composite Higgs Models after Run 2, Adv. High Energy Phys. 2018, 7168480 (2018), doi:10.1155/2018/7168480, arXiv:1703.10190.
 [45] D. Chowdhury and O. Eberhardt, Update of Global TwoHiggsDoublet Model Fits, JHEP 05, 161 (2018), doi:10.1007/JHEP05(2018)161, arXiv:1711.02095.
 [46] J. Ellis, C. W. Murphy, V. Sanz and T. You, Updated Global SMEFT Fit to Higgs, Diboson and Electroweak Data, JHEP 06, 146 (2018), doi:10.1007/JHEP06(2018)146, arXiv:1803.03252.
 [47] A. Biekötter, T. Corbett and T. Plehn, The GaugeHiggs Legacy of the LHC Run II, SciPost Phys. 6, 064 (2019), doi:10.21468/SciPostPhys.6.6.064, arXiv:1812.07587.
 [48] P. M. Ferreira, J. F. Gunion, H. E. Haber and R. Santos, Probing wrongsign Yukawa couplings at the LHC and a future linear collider, Phys. Rev. D89(11), 115003 (2014), doi:10.1103/PhysRevD.89.115003, arXiv:1403.4736.
 [49] G. Belanger, B. Dumont, U. Ellwanger, J. F. Gunion and S. Kraml, Status of invisible Higgs decays, Phys. Lett. B723, 340 (2013), doi:10.1016/j.physletb.2013.05.024, arXiv:1302.5694.
 [50] J. Bernon, J. F. Gunion, H. E. Haber, Y. Jiang and S. Kraml, Scrutinizing the alignment limit in twoHiggsdoublet models: m=125 GeV, Phys. Rev. D92(7), 075004 (2015), doi:10.1103/PhysRevD.92.075004, arXiv:1507.00933.
 [51] J. Bernon, J. F. Gunion, H. E. Haber, Y. Jiang and S. Kraml, Scrutinizing the alignment limit in twoHiggsdoublet models. II. m=125 GeV, Phys. Rev. D93(3), 035027 (2016), doi:10.1103/PhysRevD.93.035027, arXiv:1511.03682.
 [52] D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. Da Silva, S. Kraml, U. Laa and A. Pukhov, Collider limits on new physics within micrOMEGAs_4.3, Comput. Phys. Commun. 222, 327 (2018), doi:10.1016/j.cpc.2017.08.028, arXiv:1606.03834.