The measurement of the Higgs self-coupling at the LHC: theoretical status

Now that the Higgs boson has been observed by the ATLAS and CMS experiments at the LHC, the next important step would be to measure accurately its properties to establish the details of the electroweak symmetry breaking mechanism. Among the measurements which need to be performed, the determination of the Higgs self-coupling in processes where the Higgs boson is produced in pairs is of utmost importance. In this paper, we discuss the various processes which allow for the measurement of the trilinear Higgs coupling: double Higgs production in gluon fusion, vector boson fusion, double Higgs-strahlung and associated production with a top quark pair. We first evaluate the production cross sections for these processes at the LHC with center-of-mass energies ranging from the present \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \sqrt{s}=8 $\end{document} TeV to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \sqrt{s}=100 $\end{document} TeV, and discuss their sensitivity to the trilinear Higgs coupling. We include the various higher order QCD radiative corrections, at next-to-leading order for gluon and vector boson fusion and at next-to-next-to-leading order for associated double Higgs production with a gauge boson. The theoretical uncertainties on these cross sections are estimated. Finally, we discuss the various channels which could allow for the detection of the double Higgs production signal at the LHC and estimate their potential to probe the trilinear Higgs coupling.


Introduction
A bosonic particle with a mass of about 125 GeV has been observed by the ATLAS and CMS Collaborations at the LHC [1][2][3][4] and it has, grosso modo, the properties of the long sought Higgs particle predicted in the Standard Model (SM) [5][6][7][8][9]. This closes the first chapter of the probing of the mechanism that triggers the breaking of the electroweak symmetry and generates the fundamental particle masses. Another, equally important chapter is now opening: the precise determination of the properties of the produced particle. This is of extreme importance in order to establish that this particle is indeed the relic of the mechanism responsible for the electroweak symmetry breaking and, eventually, to pin down effects of new physics if additional ingredients beyond those of the SM are involved in the symmetry breaking mechanism. To do so, besides measuring the mass, the total decay JHEP04(2013)151 width and the spin-parity quantum numbers of the particle, a precise determination of its couplings to fermions and gauge bosons is needed in order to verify the fundamental prediction that they are indeed proportional to the particle masses. Furthermore, it is necessary to measure the Higgs self-interactions. This is the only way to reconstruct the scalar potential of the Higgs doublet field Φ, that is responsible for spontaneous electroweak symmetry breaking, with v = 246 GeV. Rewriting the Higgs potential in terms of a physical Higgs boson leads to the trilinear Higgs self-coupling λ HHH , which in the SM is uniquely related to the mass of the Higgs boson, This coupling is only accessible in double Higgs production [10][11][12][13][14][15][16][17][18]. One thus needs to consider the usual channels in which the Higgs boson is produced singly [19][20][21][22], but allows for the state to be off mass-shell and to split up into two real Higgs bosons. At hadron colliders, four main classes of processes have been advocated for Higgs pair production: a) the gluon fusion mechanism, gg → HH, which is mediated by loops of heavy quarks (mainly top quarks) that couple strongly to the Higgs boson [24][25][26][27]; b) the W W/ZZ fusion processes (VBF), qq ′ → V * V * qq ′ → HHqq ′ (V = W, Z), which lead to two Higgs particles and two jets in the final state [24,[28][29][30][31][32]; c) the double Higgs-strahlung process, qq ′ → V * → V HH (V = W, Z), in which the Higgs bosons are radiated from either a W or a Z boson [33]; d) associated production of two Higgs bosons with a top quark pair, pp → ttHH [34].
As they are of higher order in the electroweak coupling and the phase space is small due to the production of two heavy particles in the final state, these processes have much lower production cross sections, at least two orders of magnitude smaller, compared to the single Higgs production case. In addition, besides the diagrams with H * → HH splitting, there are other topologies which do not involve the trilinear Higgs coupling, e.g. with both Higgs bosons radiated from the gauge boson or fermion lines, and which lead to the same final state. These topologies will thus dilute the dependence of the production cross sections for double Higgs production on the λ HHH coupling. The measurement of the trilinear Higgs coupling is therefore an extremely challenging task and very high collider luminosities as well as high energies are required. We should note that to probe the quadrilinear Higgs coupling, λ HHHH = 3M 2 H /v 2 , which is further suppressed by a power of v compared to the triple Higgs coupling, one needs to consider triple Higgs production processes [10][11][12][35][36][37]. As their cross sections are too small to be measurable, these processes are not viable in a foreseen future so that the determination of this last coupling seems hopeless.

JHEP04(2013)151
In ref. [15][16][17], the cross sections for the double Higgs production processes and the prospects of extracting the Higgs self-coupling have been discussed for the LHC with a 14 TeV center-of-mass (c.m.) energy in both the SM and its minimal supersymmetric extension (MSSM) where additional channels occur in the various processes.
In the present paper, we update the previous analysis. In a sense, the task is made easier now that the Higgs boson mass is known and can be fixed to M H ≈ 125 GeV. However, lower c.m. energies have to be considered such as the current one, √ s = 8 TeV.
In addition, there are plans to upgrade the LHC which could allow to reach c.m. energies of about 30 TeV [38] and even up to 100 TeV. These very high energies will be of crucial help to probe these processes.
Another major update compared to ref. [15][16][17] is that we will consider all main processes beyond leading order (LO) in perturbation theory, i.e. we will implement the important higher order QCD corrections. In the case of the gluon fusion mechanism, gg → HH, the QCD corrections have been calculated at next-to-leading-order (NLO) in the low energy limit in ref. [39]. They turn out to be quite large, almost doubling the production cross section at √ s = 14 TeV, in much the same manner as for single Higgs production [40][41][42][43][44][45][46].
In fact, the QCD corrections for single and double Higgs productions are intimately related and one should expect, as in the case of gg → H, a further increase of the total cross section by ≈ 30% once the next-to-next-to-leading (NNLO) corrections are also included [47][48][49].
It is well known that for single Higgs production in the vector boson fusion process qq ′ → Hqq ′ there is no gluon exchange between the two incoming/outgoing quarks as the initial and final quarks are in color singlet states at LO. Then the NLO QCD corrections consist simply of the known corrections to the structure functions [50][51][52]. The same can be said in the case of double Higgs production qq ′ → HHqq ′ , and in this paper we will implement the NLO QCD corrections to this process in the structure function approach. The NNLO corrections in this approach turn out to be negligibly small for single Higgs production [53,54] and we will thus ignore them for double Higgs production.
In the single Higgs-strahlung process, qq ′ → V * → V H, the NLO QCD corrections can be inferred from those of the Drell-Yan process qq ′ → V * [55][56][57]. This can be extended to NNLO [47,58,59] but, in the case of ZH production, one needs to include the gg initiated contribution, gg → ZH [59][60][61] as well as some additional subleading corrections [62]. The same is also true for double Higgs-strahlung and we will include in this paper the Drell-Yan part of the corrections up to NNLO. In the case of ZHH final states, we will determine the additional contribution of the pentagon diagram gg → ZHH which turns out to be quite substantial, increasing the total cross section by up to 30% at √ s = 14 TeV.
In the case of the pp → ttHH process, the determination of the cross section at LO is already rather complicated. We will not consider any correction beyond this order (that, in any case, has not been calculated) and just display the total cross section without further analysis. We simply note that the QCD corrections in the single Higgs case, pp → ttH, turn out to be quite modest. At NLO, they are small at √ s = 8 TeV and increase the cross section by less than ≈ 20% at √ s = 14 TeV [63][64][65][66]. We also note that this channel is plagued by huge QCD backgrounds.

JHEP04(2013)151
In addition, the electroweak radiative corrections to these double Higgs production processes have not been calculated yet. Nevertheless, one expects that they are similar in size to those affecting the single Higgs production case, which are at the few percent level at the presently planned LHC c.m. energies [67][68][69][70][71][72][73][74][75][76] (see also ref. [77] for a review). They should thus not affect the cross sections in a significant way and we will ignore this issue in our analysis.
After determining the K-factors, i.e. the ratios of the higher order to the lowest order cross sections consistently evaluated with the value of the strong coupling α s and the parton distribution functions taken at the considered perturbative order, a next step will be to estimate the theoretical uncertainties on the production cross sections in the various processes. These stem from the variation of the renormalization and factorization scales that enter the processes (and which gives a rough measure of the missing higher order contributions), the uncertainties in the parton distribution functions (PDF) and the associated one on the strong coupling constant α s and, in the case of the gg → HH process, the uncertainty from the use of an effective approach with an infinitely heavy virtual top quark, to derive the NLO corrections (see also ref. [78,79]). This will be done in much the same way as for the more widely studied single Higgs production case [77,80].
Finally, we perform a preliminary analysis of the various channels in which the Higgs pair can be observed at the LHC with a c.m. energy of √ s = 14 TeV assuming up to 3000 fb −1 collected data, and explore their potential to probe the λ HHH coupling. Restricting ourselves to the dominant gg → HH mechanism in a parton level approach, 1 we first discuss the kinematics of the process, in particular the transverse momentum distribution of the Higgs bosons and their rapidity distribution at leading order. We then evaluate the possible cross sections for both the signal and the major backgrounds. As the Higgs boson of a mass around 125 GeV dominantly decays into b-quark pairs with a branching ratio of ≈ 60% and other decay modes such as H → γγ and H → W W * → 2ℓ2ν are rare [89,90], and as the production cross sections are already low, we will focus on the three possibly promising detection channels gg → HH → bbγγ, bbττ and bbW + W − . Very high luminosities, O(ab −1 ) would be required to have some sensitivity on the λ HHH coupling. The rest of the paper is organized as follows. In the next section, we discuss the QCD radiative corrections to double Higgs production in the gluon fusion, vector boson fusion and Higgs-strahlung processes (the ttHH process will be only considered at tree-level) and how they are implemented in the programs HPAIR [91], VBFNLO [92] and a code developed by us to evaluate the inclusive cross sections in Higgs-strahlung processes. In section 3, we evaluate the various theoretical uncertainties affecting these cross sections and collect at M H = 125 GeV the double Higgs production cross sections at the various LHC energies. We also study the sensitivity in the different processes to the trilinear Higgs self-coupling.
1 Early and more recent parton level analyses of various detection channels have been performed in refs. [81][82][83][84][85][86] with the recent ones heavily relying on jet-substructure techniques [87]. However, a full and realistic assessment of the LHC to probe the trilinear coupling would require the knowledge of the exact experimental conditions with very high luminosities and a full simulation of the detectors which is beyond the scope of this paper. Only the ATLAS and CMS Collaborations are in a position to perform such detailed investigations and preliminary studies have already started [88].
Associated production with top-quarks: qq/gg → ttHH Section 4 will be devoted to a general discussion of the channels that could allow for the detection of the two Higgs boson signal at a high-luminosity 14 TeV LHC, concentrating on the dominant gg → HH process, together with an analysis of the major backgrounds. A short conclusion is given in the last section.

Higgs pairs at higher orders in QCD
Generic Feynman diagrams for the four main classes of processes leading to double Higgs production at hadron colliders, gluon fusion, W W/ZZ fusion, double Higgs-strahlung and associated production with a top quark pair, are shown in figure 1. As can be seen in each process, one of the Feynman diagrams involves the trilinear Higgs boson coupling, λ HHH = 3M 2 H /v, which can thus be probed in principle. The other diagrams involve the couplings of the Higgs boson to fermions and gauge bosons and are probed in the processes where the Higgs particle is produced singly.

JHEP04(2013)151
In this section we will discuss the production cross sections for the first three classes of processes, including the higher order QCD corrections. We will first review the gluon channel and then we will move on to the higher-order corrections in the weak boson fusion and Higgs-strahlung channels.

The gluon fusion process
The gluon fusion process is -in analogy to single Higgs production -the dominant Higgs pair production process. The cross section is about one order of magnitude larger than the second largest process which is vector boson fusion. As can be inferred from figure 1a) it is mediated by loops of heavy quarks which in the SM are mainly top quarks. Bottom quark loops contribute to the total cross section with less than 1% at LO.
The process is known at NLO QCD in an effective field theory (EFT) approximation by applying the low energy theorem (LET) [40][41][42][43][44][45][46][93][94][95] which means that effective couplings of the gluons to the Higgs bosons are obtained by using the infinite quark mass approximation. The hadronic cross section at LO is given by with s being the hadronic c.m. energy, τ 0 = 4M 2 H /s, and f g the gluon distribution function taken at a typical scale µ F specified below. The partonic cross section at LO,σ LO , can be cast into the form withŝ andt denoting the partonic Mandelstam variables. The triangular and box form factors F △ , F and G approach constant values in the infinite top quark mass limit, The expressions with the complete mass dependence are rather lengthy and can be found in ref. [27] as well as the NLO QCD corrections in the LET approximation in ref. [39]. The full LO expressions for F △ , F and G are used wherever they appear in the NLO corrections in order to improve the perturbative results, similar to what has been done in the single Higgs production case where using the exact LO expression reduces the disagreement between the full NLO result and the LET result [19][20][21][22][40][41][42][43][44][45][46].
For the numerical evaluation we have used the publicly available code HPAIR [91] in which the known NLO corrections are implemented. As a central scale for this process we choose

JHEP04(2013)151
where M HH denotes the invariant mass of the Higgs boson pair. This is motivated by the fact that the natural scale choice in the process gg → H is µ 0 = M H . Extending this to Higgs pair production naturally leads to the scale choice of eq. (2.5). The motivation to switch to µ 0 = 1/2 M H in single Higgs production comes from the fact that it is a way to acccount for the ∼ +10% next-to-next-to-leading logarithmic (NNLL) corrections [77,[96][97][98] in a fixed order NNLO calculation. It also improves the perturbative convergence from NLO to NNLO [99]. Still NNLO and NNLL calculations for gg → HH process are not available at the moment, not to mention an exact NLO calculation that would be the starting point of further improvements. It then means that there is no way to check wether these nice features appearing in single Higgs production when using µ 0 = 1/2 M H would still hold in the case of Higgs pair production when using µ 0 = 1/2 M HH . We then stick to the scale choice of eq. (2.5). The K-factor, describing the ratio of the cross section at NLO using NLO PDFs and NLO α s to the leading order cross section consistently evaluated with LO PDFs and LO α s , for this process is

The vector boson fusion process
The structure of the Higgs pair production process through vector boson fusion [28][29][30][31] is very similar to the single Higgs production case. The vector boson fusion process can be viewed as the double elastic scattering of two (anti)quarks with two Higgs bosons radiated off the weak bosons that fuse. In particular this means that the interference with the double Higgs-strahlung process qq ′ → V * HH → qq ′ HH is negligible and this latter process is treated separately. This is justified by the kinematics of the process with two widely separated quark jets of high invariant mass and by the color flow of the process. This leads to the structure function approach that has been applied with success to calculate the QCD corrections in the vector boson fusion production of a single Higgs boson [50][51][52][53][54]. Generic diagrams contributing at NLO QCD order are shown in figure 2. For simplicity only the diagrams with the QCD corrections to the upper quark line are shown. The calculation involving the second quark line is identical. The blob of the vertex V V HH is a shortcut for the diagrams depicted in figure 3, which include charged currents (CC) with W ± bosons and neutral currents (NC) with a Z boson exchange. As can be seen only one of the three diagrams involves the trilinear Higgs coupling. The other diagrams act as irreducible background and lower the sensitivity of the production process to the Higgs self-coupling. We have calculated the NLO QCD corrections in complete analogy to the single Higgs VBF process [51]. The real emission contributions are given by a gluon attached to the quark lines either in the initial or the final state and from the gluon-quark initial state. As we are working in the structure function approach, the corrections of the upper and lower quark lines do not interfere and are simply added incoherently. The amplitudes have the following structure, Figure 2. Generic diagrams contributing to the NLO corrections to qq ′ → HHqq ′ . Shown are the LO diagram (upper left) and the NLO corrections for the upper quark line. The blob of the V V HH vertex is a shortcut for the three diagrams shown in figures 1b) and 3. where T µν V * V * stands for the tensor structure of the diagrams depicted in figure 3 and J q,q ′ µ are the quark currents of the upper and lower lines, respectively, with four-momenta q, q ′ . The calculation is done numerically using the Catani-Seymour dipole subtraction method [100] to regularize the infrared divergencies. The formulae for the subtraction terms as well as the finite corrections are identical to the ones for single Higgs VBF production as only the quark currents are involved. They can be found in ref. [51].
We have implemented this calculation in the VBFNLO code [92] in which we have provided the tensor structure depicted in figure 3 which has been calculated with Mad-Graph [101]. Up to now the VBFNLO implementation only involves on-shell Higgs pairs. We have found an increase of ∼ +7% of the total cross section compared to the LO result when using the central scale with Q V * being the momentum of the exchanged weak bosons (V * = W * , Z * ). 2 This result is in agreement with a previous calculation done in the context of the two Higgs doublet model [102].

The Higgs-strahlung process
The production of a Higgs pair in association with a vector boson has been calculated for the first time quite a while ago [33] and shares common aspects with the single Higgsstrahlung process. The NLO corrections can be implemented in complete analogy to single Higgs-strahlung [55][56][57]. We will update in this paper the former results and present the NNLO corrections to the W HH and ZHH inclusive production cross sections. These calculations have been implemented in a code which shall become publicly available. At LO the process pp → V HH (V = W, Z) is given by quark-antiquark annihilations in s-channel mediated processes involving three Feynman diagrams, see figure 1c). As can be seen only one of the three diagrams involves the trilinear Higgs coupling. The sensitivity to this coupling is then diluted by the remaining diagrams. After integrating over the azimuthal angle we are left with the following partonic differential cross section withŝ being the partonic c.m. energy (see also ref. [15][16][17]), where we use of the following notation, and the reduced couplings of the quarks to the vector bosons, a q = v q = √ 2 for V = W and any quark q, a u = 1 and v u = 1 − 8/3 sin 2 θ W for q = u, s and V = Z, a d = −1 and v d = −1 + 4/3 sin 2 θ W for q = d, c, b and V = Z. The coefficients f i as well as C HHH are The coefficient C HHH includes the trilinear Higgs coupling λ HHH . Figure 4. Feynman diagrams contributing to the NLO QCD corrections for Drell-Yan production.

JHEP04(2013)151
In order to obtain the full hadronic section, the differential partonic cross section of eq. (2.9) is convoluted with the quark parton distribution functions, f q , f q ′ taken at a typical scale µ F specified below: where s stands for the hadronic c.m. energy and τ 0 = (2M H + M V ) 2 /s. The total partonic cross sectionσ V HH has been obtained after the integration of eq. (2.9) over x 1 , x 2 . The calculation of the NLO QCD corrections is similar to the single Higgs-strahlung case. In fact, this process can be viewed as the Drell-Yan production pp → V * followed by the splitting process V * → V HH. The off-shell vector boson can have any momentum k 2 with (2M H + M V ) 2 ≤ k 2 ≤ŝ. This factorization is in principle valid at all orders for the Drell-Yan like contributions and leads, after folding with the PDF, to (2.14) In the expressions above ij stands for any initial partonic subprocess, ∆ ij is the Drell-Yan correction, z = k 2 /ŝ andσ LO stands for the LO partonic cross section of the process figure 4 the generic diagrams which contribute at NLO to the Drell-Yan process qq ′ → V * are depicted. The NLO QCD corrections increase the total cross section by ∼ +17% at 14 TeV for M H = 125 GeV.
We have calculated the NNLO corrections, which have not been available so far, in the same way except for the process involving a Z boson. In fact there are additional contributions that are specific to the case of a Z boson, involving an effective Zgg vertex. Similar to what is stated in ref. [59] for the single Higgs production case, only the specific gluon fusion initiated process will be of non-negligible contribution and will be described below. Let us start with the NNLO QCD Drell-Yan contribution. Some generic diagrams contributing to the NNLO corrections to qq ′ → V * are shown in figure 5. We apply the procedure as described by eq. (2.13) and the expression is then given by Figure 5. Some Feynman diagrams contributing at NNLO QCD to Drell-Yan production.
The expressions for the coefficients ∆ (i=1,2) (z) refer to the NLO and NNLO corrections, respectively. As they are too lengthy to be reproduced here, we refer the reader to the appendix B of ref. [58] and to ref. [47]. The expressions given there have to be rescaled by a factor of (π/α s ) i , and M ≡ µ F , R ≡ µ R . In our calculation we have included the full CKM matrix elements in the quark luminosity as well as the initial bottom quark contribution. Figure 6. Some generic diagrams contributing to gg → ZHH. For the triangle+box topologies, only those involving the trilinear Higgs couplings are depicted.

JHEP04(2013)151
We use the central scale where M V HH denotes the invariant mass of the V HH system. The Drell-Yan NNLO QCD corrections eq. (2.16) turn out to be very small. They typically increase the cross section by a few percent at 14 TeV.
The last contribution ∆σ gg→ZHH , see diagrams in figure 6, is only present in the case of Higgs pair production in association with a Z boson. It stems from the process gg → ZHH which is loop-mediated already at LO. Being of order α 2 s it contributes to the total cross section pp → ZHH at NNLO QCD. The process is mediated by quark loops in triangle, box and pentagon topologies. In the latter two topologies, only top and bottom quarks contribute as the Yukawa couplings to light quarks can be neglected. At the LHC the contribution of the gluon fusion channel is substantial in contrast to the single Higgs production case. Indeed, while in the latter the contribution is of order ∼ +5% compared to the NNLO QCD Drell-Yan contribution, in the case of Higgs pair production it contributes with ∼ +20 · · · + 30% depending on the c.m. energy. This enhancement can be explained by the additional pentagon topology which a) involves two top Yukawa couplings and b) softens the destructive interference between the triangle and box diagrams that is present in the single Higgs production case. Furthermore, the suppression by a power (α s /π) 2 is partly compensated by the increased gluon luminosity at high energies. This explains why this channel, which has been calculated using FeynArts/FormCalc [103][104][105][106], should be taken into account. It also implies that the scale variation in pp → ZHH will be worse than in pp → W HH because of the O(α 2 s ) gluon fusion mechanism appearing at NNLO.

Cross sections and sensitivity at the LHC
In this section we will present the results for the calculation of the total cross sections including the higher-order corrections discussed in the previous section as well as the various related theoretical uncertainties. We will use the MSTW2008 PDF set [107] as our reference set. We choose the following values for the W , Z and top quark masses and for the strong coupling constant at LO, NLO and NNLO, The electromagnetic constant α is calculated in the G µ scheme from the values of M W and M Z given above. For the estimate of the residual theoretical uncertainties in the various Higgs pair production processes we considered the following uncertainties: 1. the scale uncertainty, stemming from the missing higher order contributions and estimated by varying the renormalization scale µ R and the factorization scale µ F in the interval 1 2 µ 0 ≤ µ R , µ F ≤ 2µ 0 with some restrictions on the ratio µ R /µ F depending on the process; 2. the PDF and related α s errors. The PDFs are non-perturbative quantities fitted from the data and not calculated from QCD first principles. It is then compulsory to estimate the impact of the uncertainties on this fit and on the value of the strong coupling constant α s (M 2 Z ) which is also fitted together with the PDFs; 3. in the case of the gluon fusion process there is a third source of uncertainties which comes from the use of the effective field theory approximation to calculate the NLO QCD corrections, where top loops are taken into account in the infinite top mass approximation and bottom loops are neglected.
In the following we will present results for M H = 125 GeV. Note that the results for the total cross sections and uncertainties are nearly the same for M H = 126 GeV. The total cross sections at the LHC for the four classes of Higgs pair production processes are shown in figure 7 as a function of the c.m. energy. For all processes the numerical uncertainties are below the permille level and have been ignored. The central scales which have been As can be inferred from the figure and also seen in table 1 the largest cross section is given by the gluon fusion channel which is one order of magnitude larger than the vector boson fusion cross section. All processes are ∼ 1000 times smaller than the corresponding single Higgs production channels, implying that high luminosities are required to probe the Higgs pair production channels at the LHC.  3.1 Theoretical uncertainties in the gluon channel 3.1.1 Theoretical uncertainty due to missing higher order corrections

JHEP04(2013)151
The large K-factor for this process of about 1.5 − 2 depending on the c.m. energy shows that the inclusion of higher order corrections is essential. An estimate on the size of the uncertainties due to the missing higher order corrections can be obtained by a variation of the factorization and renormalization scales of this process. In analogy to single Higgs production studies [77,80] we have estimated the error due to missing higher order corrections by varying µ R , µ F in the interval As can be seen in figure 8 we find sizeable scale uncertainties ∆ µ of order ∼ +20%/−17% at 8 TeV down to +12%/−10% at 100 TeV. Compared to the single Higgs production case the scale uncertainty is twice as large [77,80]. However, this should not be a surprise as there are NNLO QCD corrections available for the top loop (in a heavy top mass expansion) in the process gg → H while they are unknown for the process gg → HH. be to compare different parameter sets, such as MSTW [107], CT10 [108], ABM11 [109], GJR08 [110], HERA 1.5 [111] and NNPDF 2.3 [112]. This is exemplified in figure 9 where the predictions using the six previous PDF sets are displayed. As can be seen there are large discrepancies over the whole considered c.m. energy range. At low energies the smallest prediction comes from ABM11 which is ∼ 22% smaller than the prediction made with the MSTW set while at high energies ABM11 and MSTW lead to similar results whereas the result obtained with the GJR08 set deviates by ∼ −15%. The CT10 predictions show about −5% difference over the whole energy range with respect to the cross section obtained with MSTW while the HERA prediction starts from lower values and eventually reaches the CT10 result. Finally the cross sections calculated with the NNPDF set decrease over the energy range, starting from being very similar to the MSTW result to reach at √ s = 100 TeV the one calculated with CT10.
Another source of uncertainty due to the PDF sets comes from the experimental uncertainties on the fitted data. The so-called Hessian method, used by the MSTW collaboration, provides additional PDF sets next to the best-fit PDF. Additional 2N P DF sets reflect the ±1σ variation around the minimal χ 2 of all N P DF parameters that enter the fit. Using the 90% CL error PDF sets provided by the MSTW collaboration a PDF error of about 6% is obtained for √ s = 8 TeV. The uncertainty shrinks to ∼ 2% for √ s = 100 TeV.
In addition to the PDF uncertainties, there is also an uncertainty due to the errors on the value of the strong coupling constant α s . The MSTW collaboration provides additional PDF sets such that the combined PDF+α s uncertainties can be evaluated [113]. At NLO the MSTW PDF set uses α s (M Z ) = 0.12018 +0.00122 −0.00151 (at 68% CL) or +0.00317 −0.00386 (at 90% CL) . The combined PDF and α s error is much larger than the pure PDF error. At 8 TeV the PDF error of +5.8%/−6.0% rises to a combined error of +8.5%/−8.3%. At 33 TeV the rise is even larger -from the pure PDF error of +2.5%/−2.7% to the combined PDF+α s error of +6.2%/−5.4%.
There is also a theoretical uncertainty on α s stemming from scale variation or ambiguities in the heavy flavour scheme definition. The MSTW collaboration estimates this uncertainty for α s at NLO to ∆α s (M Z ) = ±0.003 [113]. However this uncertainty is already included in the scale uncertainty on the input data sets used to fit the PDF and has been taken into account in the definition of the PDF+α s uncertainty. Thus, it does not have to be taken into account separately and the combined PDF+α s error calculated with the MSTW 2008 PDF set will be our default PDF+α s uncertainty.
However, even if all these uncertainties for the MSTW PDF set are taken into account, the different PDF set predictions do not agree. There might be agreement if also uncertainties of the other PDF sets are taken into account, as done in ref. [80] for the case of single Higgs production. This means that the PDF uncertainty might be underestimated, but this issue is still an open debate (see for example ref. [114] for a new discussion about theoretical issues in the determination of PDFs) and improvements may come with the help of new LHC data taken into account in the fits of the various PDF collaborations.

3.1.3
The effective theory approach The last source of theoretical errors that we consider is the use of the LET for the calculation of the NLO corrections. At LO it was found in ref. [ Such an expansion works very well for single Higgs production, since √ŝ = M H (at LO) for the production of an on-shell Higgs boson whereas in Higgs pair production we have This means that for Higgs pair production m top ≫ √ŝ is never fulfilled for M H = 125 GeV so that the LET approximation is not valid at LO [78,79].
At NLO, however, the LET approximation works much better in case the LO cross section includes the full mass dependence. The reason is that the NLO corrections are dominated by soft and collinear gluon effects. They factorize in the Born term and in the NLO correction contributions, meaning that the K-factor is not strongly affected from any finite mass effects. Based on the results for the single Higgs case [40][41][42][43][44][45][46] where the deviation between the exact and asymptotic NLO results amounts to less than 7% for M H < 700 GeV, we estimate the error from applying an effective field theory approach for the calculation of the NLO corrections to 10%.

Total uncertainty
In order to obtain the total uncertainty we follow the procedure advocated in ref. [115]. Since quadratic addition is too optimistic (as stated by the LHC Higgs Cross section Working Group, see ref. [77]), and as the linear uncertainty might be too conservative, the procedure adopted is a compromise between these two ways of combining the individual theoretical uncertainties. We first calculate the scale uncertainty and then add on top of that the PDF+α s uncertainty calculated for the minimal and maximal cross sections with respect to the scale variation. The LET error is eventually added linearly. This is shown in figure 10 where we display the total cross section including the combined theoretical uncertainty. It is found to be sizeable, ranging from ∼ +42%/ − 33% at 8 TeV down to ∼ +30%/−25% at 100 TeV. The numbers can be found in table 2.

VBF and Higgs-strahlung processes
We will not repeat the detailed description of the previous uncertainties in this subsection and only summarize how they affect the VBF and Higgs-strahlung inclusive cross JHEP04(2013)151 sections. In both channels, only the scale uncertainties and the PDF+α s errors are taken into account, the calculation being exact at a given order.

The VBF channel
As stated in section 2 we use the central scale µ 0 = µ R = µ F = Q V * , that is the momentum transfer of the exchanged weak boson. Note that a cut of Q V * ≥ 2 GeV has to be applied as stated in the previous section. The scale uncertainty is calculated in exactly the same way as for the gluon fusion mechanism, exploring the range µ 0 /2 ≤ µ R , µ F ≤ 2µ 0 . We have checked that imposing the restriction 1/2 ≤ µ R /µ F ≤ 2 does not modify the final result. We obtain very small scale uncertainties ranging from ∼ ±2% at 8 TeV down to ∼ ±1% and even lower at 33 TeV as can be seen in figure 11 (left). This is in sharp contrast with the ±8% uncertainty obtained at LO at 8 TeV for example, which illustrates the high level of precision already obtained with NLO QCD corrections.
The 90% CL PDF+α s uncertainties are calculated following the recipe presented in the gluon fusion subsection. The PDF+α s uncertainty dominates the total error, ranging from +7%/−4% at 8 TeV down to +5%/−3% at 100 TeV.
The total error has been obtained by adding linearly the scale and PDF+α s uncertainties, given the small variation of the cross section with respect to the choice of the scale. This process has a total theoretical uncertainty which is always below 10%, from +9%/−6% at 8 TeV to +6%/−4% both at 33 and 100 TeV as can be read off table 3. The total uncertainty is displayed in figure 11 (right) as a function of the c.m. energy. The QCD corrections drastically reduce the residual theoretical uncertainty.   Table 3. The total Higgs pair production cross section at NLO in the vector boson fusion process at the LHC (in fb) for given c.m. energies (in TeV) at the central scale µ F = µ R = Q V * for M H = 125 GeV. The corresponding shifts due to the theoretical uncertainties from the various sources discussed are also shown as well as the total uncertainty when all errors are added linearly.

The associated Higgs pair production with a vector boson
The cross section is calculated with the central scale µ 0 = µ R = µ F = M V HH which is the invariant mass of the W/Z + Higgs pair system. The scales are varied in the interval µ 0 /2 ≤ µ R = µ F ≤ 2µ 0 . The factorization and renormalization scales can be chosen to be the same as the impact of taking them differently is totally negligible, given the fact that the scale µ R only appears from NLO on and that we have a NNLO calculation which then reduces any non-negligible contribution arising from the difference between renormalization and factorization scales. As noticed in section 2, the scale uncertainty is expected to be worse in the ZHH channel because of the significant impact of the gluon fusion contribution. This is indeed the case as we obtain a scale uncertainty below ±0.5% in the W HH channel whereas the uncertainty in the ZHH channel is ∆ µ ∼ ±3% at 8 TeV and slightly more at higher energies to reach ∼ ±5% at 33 TeV, as can be seen in figure 12 (left). The total PDF+α s error that we obtain is very similar for the two channels W HH and ZHH. It varies from ∼ ±4% at 8 TeV down to ∼ ±3% at 33 TeV, with a slightly higher uncertainty at 100 TeV.   The total error has been obtained exactly as in the VBF case, given the very small variation of the cross section with respect to the scale choice. It is dominated by the PDF+α s uncertainty and amounts to +5%/ − 4% (+4%/ − 4%) at 8 (100) TeV in the W HH channel and +7%/−5% (+8%/−8%) at 8 (100) TeV in the ZHH channel. The total theoretical uncertainty in the Higgs-strahlung channels is less than 10%. The numbers are given in tables 4 and 5. The total uncertainty bands for the W HH and ZHH channels are displayed in figure 12 (right).

Sensitivity to the trilinear Higgs coupling in the main channels
We end this section by a brief study of the sensitivity of the three main channels to the trilinear Higgs coupling that we want to probe. Indeed, as can be seen in figure 1, all processes do not only involve a diagram with the trilinear Higgs couplings but also additional contributions which then dilute the sensitivity. In order to study the sensitivity within the SM, we rescale the coupling λ HHH in terms of the SM trilinear Higgs  as λ HHH = κλ SM HHH . This is in the same spirit as the study done in ref. [15][16][17] and its goal is to give a way to estimate the precision one could expect in the extraction of the SM trilinear Higgs coupling from HH measurements at the LHC. In particular the variation of the trilinear Higgs coupling should not be viewed as coming from some beyond the SM physics model and it should be noted that quite arbitrary deviations of the trilinear Higgs couplings emerge from non-vanishing higher-dimension operators starting with dimension 6.
In figure 13 the three main Higgs pair production cross sections are displayed as a function of κ for three different c.m. energies √ s = 8, 14 and 33 TeV. The left panels show the total cross sections while the right panels show the ratio between the cross sections at a given λ HHH = κλ SM HHH and the SM cross section with κ = 1. As can be seen, the most sensitive channel is by far the VBF production mode, in particular for low and high values of κ. The shapes of the cross sections with respect to a variation of κ are the same in all channels and at all energies with a minimum reached at κ ∼ −1, 2 and 3 for Higgs-strahlung, VBF and gluon fusion, respectively. The right panels of figure 13 also show that the sensitivity decreases when √ s increases. This is to be expected as the diagrams involving the trilinear Higgs self-coupling are mediated by s-channel propagators which get suppressed with increasing energy, so that the relative importance of these diagrams with respect to the remaining ones is suppressed. Again it is important to note that the sensitivity tested here does not give information on the sensitivity to Higgs self-couplings in models beyond the SM. It only tests within the SM how accurately the respective Higgs pair production process has to be measured in order to extract the SM trilinear Higgs self-coupling with a certain accuracy. For example the gluon fusion cross section has to be measured with an accuracy of ∼ 50% at √ s = 8 TeV in order to be able to extract the trilinear Higgs self-coupling with an accuracy of 50%, as can be inferred from figure 13 upper left. Similar discussions can be found in refs. [15][16][17]39].

Prospects at the LHC
As shown in the previous section, inclusive Higgs boson pair production is dominated by gluon fusion at LHC energies. Other processes, such as weak boson fusion, qq ′ → qq ′ HH, associated production with heavy gauge bosons, qq ′ → V HH (V = W, Z), or associated production with top quark pairs, gg/qq → ttHH, yield cross sections which are factors of 10 -30 smaller than that for gg → HH.  Figure 13. The sensitivity of the various Higgs pair production processes to the trilinear SM Higgs self-coupling at different c.m. energies. The left panels display the total cross sections, the right panels display the ratio between the cross sections at a given κ = λ HHH /λ SM HHH and the cross sections at κ = 1.
following, we examine channels where one Higgs boson decays into a b quark pair and the other Higgs boson decays into either a photon pair, gg → HH → bbγγ, into a τ pair, gg → HH → bbττ , or into an off-shell W boson pair, gg → HH → bbW * W * . Following the Higgs Cross section Working Group recommendations [77], we assume a branching ratio JHEP04(2013)151 of 57.7% for a 125 GeV Higgs boson decaying into b quarks, 0.228% for the Higgs boson decaying into a photon pair, 6.12% for the Higgs boson decaying into a τ pair and 21.50% for the Higgs boson decaying into off-shell W * bosons.
At the time of the analysis, no generator existed for the signal process, but the matrix element for Higgs pair production in the gluon fusion channel was available in the Fortran code HPAIR [39,91]. In parallel to the approach used by the program described in [116][117][118], the HPAIR matrix element was added to Pythia 6 [119]. It has been checked that the cross sections produced by HPAIR and by the Pythia 6 implementation are consistent.
All tree-level background processes are calculated using Madgraph 5 [120]. Signal and background cross sections are evaluated using the MSTW2008 parton distribution functions [107] with the corresponding value of α s at the investigated order in perturbative QCD. The effects of QCD corrections are included in our calculation via multiplicative factors which are summarized in the following subsections and have been introduced in section 3 for the signal cross sections.

Kinematical distributions of gg → HH
In this subsection the characteristic distributions of the gluon fusion process gg → HH are studied for several observables. In figure 14, we show for each of the two final state Higgs bosons the normalized distributions of the transverse momentum P T,H and the pseudorapidity η H , as well as the invariant mass M HH , the helicity angle θ ⋆ HH which is the angle between the off-shell Higgs boson, boosted back into the Higgs boson pair rest frame, and the Higgs boson pair direction, and the rapidity y HH of the Higgs boson pair. The distributions of each observable are shown for different values of the trilinear Higgs coupling λ/λ SM = 0 (green curve), 1 (red curve) and 2 (blue curve). We also include in the plots the typical background qq → ZH coming from the Higgs boson itself (black curve), the Z boson faking a Higgs boson.
The Higgs bosons from inclusive Higgs pair production are typically boosted, as we can see in the upper left plot of figure 14 where the P T,H distributions reach their maximum for P T,H ∼ 150 GeV. For λ/λ SM = 2, the triangle diagram interferes destructively with the box diagram, which explains the dip in the P T,H distribution. This high transverse momentum spectrum can also be interpreted in terms of the low pseudorapidity of the two Higgs bosons which have a typical symmetric distribution with the maximum around zero, see figure 14 upper right. The qq → ZH background has a completely different topology with less boosted Higgs and Z bosons, P T,H/Z ∼ 50 GeV, with pseudorapidity of order |η H/Z | ∼ 2 as can be seen in the upper left and right plots of figure 14. The middle left plot of figure 14 displays the distributions of the invariant mass of the Higgs boson pair. For the signal the typical value is M HH 400 GeV to be compared to a lower value of M ZH 250 GeV for the ZH background. Also note that an important depletion appears in the signal for λ/λ SM = 2 caused by the destructive interference between the triangle and box contribution. Due to the Higgs boson scalar nature a known discriminant observable is the angle θ ⋆ HH [121]. The middle right plot in figure 14 shows that signal and ZH background have similar distributions thus making this observable less discriminant than others described before but still efficient for some specific backgrounds, as we will see   in the following. From the bottom plot of figure 14, it can be inferred that the rapidity distribution of the Higgs pair is narrower for the signal than for the ZH background.
In the following the different decay channels HH → bbγγ, HH → bbττ and HH → bbW + W − will be investigated in more detail.

The bbγγ decay channel
In this subsection, the bbγγ final state for the production of two Higgs bosons with a mass of 125 GeV at √ s = 14 TeV is investigated. Earlier studies can be found in ref. [81][82][83][84].
The calculation of the signal, pp → HH → bbγγ, is performed as described above by incor-  Table 6. K-factors for gg → HH, bbγγ, ttH and ZH production at √ s = 14 TeV [77]. The Higgs boson mass is assumed to be M H = 125 GeV.
porating the matrix element extracted from the program HPAIR into Pythia 6. We include the effects of NLO QCD corrections on the signal by a multiplicative factor, K N LO = 1.88, corresponding to a 125 GeV Higgs boson and a c.m. energy of 14 TeV. Here we set the factorization and renormalization scales equal to M HH . The generated background processes are the QCD process bbγγ and the associated production of a Higgs boson with a tt pair, ttH, with the Higgs boson subsequently decaying into a photon pair and the top quarks decaying into a W boson and a b quark, as well as single Higgs-strahlung ZH with the Higgs boson decaying into γγ and the Z boson decaying into bb. The QCD corrections have been included by a multiplicative K-factor applied to the respective LO cross section. All K-factors are taken at NLO except for the single Higgs-strahlung process which is taken at NNLO QCD and the case of bbγγ in which no higher order corrections are taken into account. The various K-factors are given in table 6 and taken from ref. [77]. The factorization and renormalization scales have been set to M HH for the signal and to specific values for each process for the backgrounds.
In this analysis, the signal and background processes are generated with exclusive cuts. The basic acceptance cuts are motivated by the fact that the b quark pair and the photon pair reconstruct the Higgs mass according to the resolutions expected for ATLAS and CMS. Note that starting from this section all the plots include the decays and the acceptance cuts specific to each final state.
In detail, we veto events with leptons having soft transverse momentum p T,ℓ > 20 GeV and with a pseudorapidity |η ℓ | < 2.4 in order to reduce the ttH background. Furthermore we also veto events with QCD jets p T,jet > 20 GeV and with a pseudorapidity |η jet | < 2.4 to further reduce the ttH background. We require exactly one b quark pair and one photon pair. The b quark pair is restricted to have p T,b > 30 GeV, |η b | < 2.4 and ∆R(b, b) > 0.4, where ∆R(b, b) denotes the isolation of the two b quarks defined by the distance ∆R = (∆η) 2 + (∆φ) 2 in the pseudorapidity and azimuthal angle plane (η, φ). We consider the b-tagging efficiency to be 70%. The photon pair has to fulfill p T,γ > 30 GeV, |η γ | < 2.4 and ∆R(γ, γ) > 0.4. The two reconstructed Higgs bosons, from the b quark pair and from the photon pair, have to reproduce the Higgs boson mass within a window of 25 GeV, 112.5 GeV < M bb < 137.5 GeV, and a window of 10 GeV, 120 GeV < M γγ < 130 GeV, respectively. We require additional isolations between the b quarks and the photons being ∆R(γ, b) > 0.4.
Based on the distributions shown in figure 15, apart from the acceptance cuts we have applied more advanced cuts for this parton level analysis. We first require the reconstructed invariant mass of the Higgs pair to fulfill M HH > 350 GeV. Furthermore we remove events which do not satisfy P T,H > 100 GeV. We also constrain the pseudorapidity of the two reconstructed Higgs bosons, |η H | < 2, and the isolation between the two b jets to be ∆R(b, b) < 2.5.

JHEP04(2013)151
[GeV]  The results are collected in table 7. The local decrease of the sensitivity between the cut on M HH and the cut on P T,H is explained by the fact that we accept to have a reduced sensitivity locally during the chain of cuts in order to enhance the final significance. In the case described in this section a cut on P T,H alone reduces the sensitivity as does a cut on η H alone, but the first cut actually improves the discrimination between the signal and the background in the pseudorapidity distribution, hence allowing for a larger improvement when applying the η H cut just after the P T,H cut. Eventually all the cuts allow for an improvement of the significance by two orders of magnitude, that is the ratio of signal events S over the square root of background events B, S/ √ B. The final value for S/ √ B is 16.3 for an integrated luminosity of L = 3000 fb −1 , corresponding to 51 signal events. Therefore this channel seems promising.
A realistic assessment of the prospects for measuring the signal in the bbγγ final state depends mostly on a realistic simulation of the diphoton fake rate due to multijet production, which is the dominant background in such an analysis. Our first parton level study gives a rough idea of how promising the bbγγ channel is.
In the following we perform a full analysis including fragmentation and hadronization effects, initial and final state radiations by using Pythia 6.4 in order to assess more reliably whether the promising features of this channel survive in a real experimental environment. All the events are processed through Delphes [122], the tool which is used for the detector  simulation. For the jet reconstruction we use the anti-k T algorithm with a radius parameter of R = 0.5. We still consider a b-tagging efficiency of 70%. We keep the same acceptance cuts as before except for the transverse momentum of the reconstructed b jet and photon which we increase up to p T,b/γ > 50 GeV. We also enlarge the window for the reconstructed Higgs boson coming from the b quark pair, by requiring 100 GeV < M bb < 135 GeV. We select events with exactly two reconstructed b jets and two photons.
The final result is displayed in the last line of table 7. The final significance S/ √ B for this simulation has decreased to 6.5 for an integrated luminosity of 3000 fb −1 , corresponding to 47 events. Though low, the significance nevertheless is promising enough to trigger a real experimental analysis as can be performed only by the experimental collaborations and which is beyond the scope of this work.

The bbττ decay channel
The bbττ decay channel is promising for low mass Higgs boson pair production at the LHC and has been previously studied in ref. [81][82][83][84][85]. An important part of this analysis will depend on the ability to reconstruct the b quark pair and the τ pair. As the real experimental assessment of such a reconstruction is beyond the scope of our work we will perform in the following a parton level analysis, assuming a perfect τ reconstruction. 3 The analysis thus represents an optimistic estimate of what can be done at best to extract the trilinear Higgs self-coupling in the bbττ channel.
We consider the two QCD-QED continuum background final states bbττ and bbτ ν τ τν τ calculated at tree-level. The bbτ ν τ τν τ final state background dominantly stems from tt production with the subsequent top quark decays t → bW and W →τ ν τ . We also include backgrounds coming from single Higgs production in association with a Z boson and the subsequent decays H → ττ and Z → bb or H → bb and Z → ττ . in table 8. All K-factors are taken at NLO except for the single Higgs-strahlung process which is taken at NNLO QCD. The factorization and renormalization scales have been taken at M HH for the signal and at specific values for each background process.
Concerning the choice of our cuts, we demand exactly one b quark pair and one τ pair. The b quark pair is required to fulfill p T,b > 30 GeV and |η b | < 2.4. We assume the b-tagging efficiency to be 70% and the τ -tagging efficiency to be 50% and we neglect fake rates in both cases. The τ pair has to fulfill p T,τ > 30 GeV and |η τ | < 2.4. The reconstructed Higgs boson from the b quark pair is required to reproduce the Higgs mass within a window of 25 GeV, 112.5 GeV < M bb < 137.5 GeV. The reconstructed Higgs boson from the τ pair needs to reproduce the Higgs mass within a window of 50 GeV, 100 GeV < M ττ < 150 GeV or within a window of 25 GeV, 112.5 GeV < M ττ < 137.5 GeV, in a more optimistic scenario. In addition to these acceptance cuts we also add more advanced cuts for this parton level analysis, based on the distributions shown in figure 16 and in a similar way as what has been done in the previous bbγγ analysis. We first demand the invariant mass of the reconstructed Higgs pair to fulfill M HH > 350 GeV. In addition, we remove events which do not satisfy P T,H > 100 GeV. We do not use a cut on the pseudorapidity of the reconstructed Higgs bosons in this analysis as it would reduce the significance.
The different results of our parton level analysis are summarized in table 9. The cuts allow for an improvement of two orders of magnitude in the significance S/ √ B. The final S/ √ B is 6.71 for an integrated luminosity of L = 3000 fb −1 , corresponding to 330 signal events. We then conclude that this channel is promising. In the last line we reproduce our result for the optimistic requirement of 112.5 GeV < M τ τ < 137.5 GeV leading to the final significance S/ √ B = 9.36 for an integrated luminosity of 3000 fb −1 . Already for a planned mid-term integrated luminosity of 300 fb −1 at 14 TeV the expectations are promising with 33 signal events and a significance S/ √ B = 2.96 in the optimistic scenario.

The bbW + W − decay channel
The analysis in this channel is difficult as the leptonic W boson decays lead to missing energy in the final state. Consequently, one of the two Higgs bosons cannot a priori be reconstructed equally well as the other Higgs boson, thus reducing our capability to efficiently remove the background with the canonical acceptance cuts previously applied in the other decay channels. This channel with one lepton plus jets final state has been studied in [85,86]. We only consider here the decay W → ℓν ℓ (ℓ = e, µ) with a branching ratio of 10.8%. We take into account the continuum background which contains all processes with the bbℓν ℓ ℓν ℓ final states at tree-level, for example qq/gg → b * b * → γbZb → bbℓℓZ with    Table 9. Cross section values of the of HH signal and the various backgrounds expected at the LHC at √ s = 14 TeV, the signal to background ratio S/B and the significance S/ √ B for L = 3000 fb −1 in the bbττ channel after applying the cuts discussed in the text. the subsequent splitting Z → ν ℓνℓ . We proceed in a similar manner as in the previous analyses. We generate the signal and the backgrounds with the following parton-level cuts. We require that the b quarks fulfill p T,b > 30 GeV and |η b | < 2.4. We consider the btagging efficiency to be 70%. The leptons have to fulfill p T,ℓ > 20 GeV and |η ℓ | < 2.4. The reconstructed Higgs boson from the b quark pair has to reproduce the Higgs boson mass within a window of 25 GeV, 112.5 GeV < M bb < 137.5 GeV. We also require that the missing transverse energy respects E miss We first require the transverse mass of the lepton pair to be M T < 125 GeV. We then remove events which do not satisfy ∆φ ℓ 1 ℓ 2 < 1.2 and we also add a constraint on the angle between the two leptons, ∆θ ℓ 1 ℓ 2 < 1.0. We demand the missing transverse energy to fulfill E miss T > 120 GeV and the projected energy to satisfyẼ miss T < 80 GeV. Note that thẽ E miss T distribution displayed in figure 17 is obtained after the acceptance cuts having been  Table 10. Cross section values of the HH signal and the considered background expected at the LHC at √ s = 14 TeV, the signal to background ratio S/B and the significance S/ √ B for L=3000 fb −1 in the bbW + W − channel after applying the cuts discussed in the text. applied but before the advanced cuts. The cuts on M T , ∆φ ℓ 1 ℓ 2 , ∆θ ℓ 1 ℓ 2 andẼ miss T modify this distribution and explain why theẼ miss T cut, which would seem not to be efficient, actually improves the significance.
The results for the LHC at √ s = 14 TeV are summarized in table 10. While the cuts allow for an improvement of the significance S/ √ B by about one order of magnitude, we are still left with a very small signal to background ratio. Thus, this channel using the final states considered here is not very promising.

Conclusions
In this paper we have discussed in detail the main Higgs pair production processes at the LHC, gluon fusion, vector boson fusion, double Higgs-strahlung and associated production with a top quark pair. They allow for the determination of the trilinear Higgs self-coupling λ HHH , which represents a first important step towards the reconstruction of the Higgs potential and thus the final verification of the Higgs mechanism as the origin of electroweak symmetry breaking. We have included the important QCD corrections at NLO to gluon fusion and vector boson fusion and calculated for the first time the NNLO corrections to double Higgs-strahlung. It turns out that the gluon initiated process to ZHH production which contributes at NNLO is sizeable in contrast to the single Higgs-strahlung case. We have discussed in detail the various uncertainties of the different processes and provided numbers for the cross sections and the total uncertainties at four c.m. energies, i.e. 8, 14, 33 and 100 TeV. It turns out that they are of the order of 40% in the gluon fusion channel while they are much more limited in the vector bosons fusion and double Higgs-strahlung processes, i.e. below 10%. Within the SM we also studied the sensitivities of the double Higgs production processes to the trilinear Higgs self-coupling in order to get an estimate of how accurately the cross sections have to be measured in order to extract the Higgs self-interaction with sufficient accuracy.
In the second part of our work we have performed a parton level analysis for the dominant Higgs pair production process through gluon fusion in different final states which are JHEP04(2013)151 bbγγ, bbττ and bbW + W − with the W bosons decaying leptonically. Due to the smallness of the signal and the large QCD backgrounds the analysis is challenging. The bbW + W − final state leads to a very small signal to background ratio after applying acceptance and selection cuts so that it is not promising. On the other hand, the significances obtained in the bbγγ and bbττ final states after cuts are ∼ 16 and ∼ 9, respectively, with not too small event numbers. They are thus promising enough to start a real experimental analysis taking into account detector and hadronization effects, which is beyond the scope of our work. Performing a first simulation on the detector level for the bbγγ state shows, however, that the prospects are good in case of high luminosities. Taking into account theoretical and statistical uncertainties and using the sensitivity plot, figure 13, the trilinear Higgs self-coupling λ HHH can be expected to be measured within a factor of two.