Multi-Higgs Production and Unitarity in Vector-Boson Fusion at Future Hadron Colliders

We study multi-Higgs final states in vector boson fusion (VBF) processes at the LHC and at future proton-proton colliders, focusing on the prospects for measurements at 27 TeV and at 100 TeV. We use an effective Lagrangian which includes higher-dimensional operators in the mass eigenstates which are relevant to VBF processes and relate this to specific parameterizations and models for new physics in the Higgs sector. We derive theoretical constraints on the parameter space from the unitarity of 2 → n scattering amplitudes and apply the results to V V → hh and hhh processes. We numerically compute cross sections with appropriate VBF cuts for 14 TeV, 27 TeV, and 100 TeV respectively. The results describe the potential for constraining new physics in the Higgs sector at the LHC and at future hadron colliders. ar X iv :1 80 8. 05 53 4v 2 [ he pph ] 2 9 A ug 2 01 8


Introduction
After the discovery of the Higgs boson with mass m h ≈ 125 GeV at the LHC [1,2], the measurements of all of its properties have become essential for the search for new physics beyond the Standard Model (SM). In the SM, the Higgs boson has three types of interaction at tree level: (1) the Yukawa interactions with fermions; (2) the interaction with electroweak gauge bosons (W ± and Z ); (3) the triple and quartic Higgs self-interactions. Establishing the last type of interaction is a crucial test of our current understanding of electroweak symmetry breaking (EWSB). However, any direct measurement of the triple and quartic Higgs self-coupling involves producing two or more Higgs bosons in a single elementary process. For all possible multi-Higgs production processes, the rates are very small.
In analogy to Higgs pair-production, the triple-Higgs final state and thus the quartic Higgs self-coupling can also be studied in the VBF production channel. The VBF topology, which implies forward jets with suppressed QCD activity in the central region, improves the signalto-background ratio considerably. This process is also sensitive to a hhhV V interaction which does not exist in the SM but should be expected for a strongly interacting Higgs sector [63]. Furthermore, an anomalous hV V coupling would have a strong impact on triple Higgs production in VBF [64].
In this paper, we study multi-Higgs production in VBF processes in an effective-field theory (EFT) approach. In this framework, anomalous effects are parameterized by the coefficients of higher-dimensional operators. By investigating the consequences of S-matrix unitarity for the amplitudes V V → hh and V V → hhh, we constrain the energy-dependent parameter region where the EFT yields a physically meaningful parameterization. We use the packages WHIZARD [65] and Madgraph [66] to compute the tree-level cross sections with these anomalous couplings taken into account. Evaluating numerical results at 14 TeV,27 TeV and at 100 TeV, we obtain projections for bounds on the parameter space at future hadron colliders and compare them to the potential of the LHC. This paper is organized as follows. In Sec. 2, we establish the framework and introduce the effective Lagrangian. In Sec. 3, we consider the unitarity constraints that arise in this class of processes. Sec. 4 contains the calculation of multi-Higgs production in the context of the EFT, and we present numerical results. We conclude this paper with a discussion of our findings in Sec. 5.

Effective Lagrangian for VBF Higgs production
As a theoretical framework for studying multiple Higgs production in VBF, we consider a generic set of anomalous Higgs interactions, in analogy with our earlier paper where we investigated multiple Higgs production in gluon-gluon fusion, Ref. [61]. We introduce a parameterization in terms of a local Lagrangian which incorporates all fields that may contribute, to leading order in a derivative expansion. This effective Lagrangian respects QED and QCD gauge invariance and conserves the discrete symmetries of the SM. For the VBF class of processes, we do not need to consider anomalous contact interactions to light quarks or gluons, and we can ignore couplings to the third generation of fermions.
In the SM at the tree level, we have the relations λ 3 = λ 4 = g W,a1 = g W,a2 = g Z ,a1 = g Z ,a2 = 1 and κ 5 = κ 6 = g V,b1 = g V,b2 = g V,b3 = g W,a3 = g Z ,a3 = 0, where the subscript V denotes W, A, X , Z . It is understood that the corresponding terms have been removed from L SM , such that they are not double-counted.
The higher-order operators in the kinetic-energy term (proportional to ξ 5 and ξ 6 ) are redundant and can be eliminated by applying the equation of motion of the Higgs field or by a non-linear transformation [67]. To eliminate ξ 5 , we may replace h → h + a 2v h 2 and get parameter shifts such as λ 3 → λ 3 + a λ 4 → λ 4 + a 2 + 6aλ 3 (2.5)  4 . To facilitate the comparison between parameterizations, we choose to retain these parameters in Table 2 below (cf. also Ref. [61]), but do not include ξ 5,6 in the numerical calculations.
symmetry violation in present electroweak precision data. The assumption of custodial symmetry should be independently confirmed or discarded by analyzing high-precision electroweak and Higgs data that will be collected at the HL-LHC and at later experiments. The phenomenological Lagrangian (2.1) provides a robust parameterization of new physics in the Higgs-electroweak sector, under the condition that no new on-shell states appear in the kinematically accessible range. In Table 1 we summarize the dependency of various Higgs production processes on the coefficients in the effective Lagrangian, as can be read off from the contributing Feynman diagrams.
The Lagrangian (2.1) does not manifestly exhibit electroweak gauge invariance. We may express all anomalous Higgs interactions in terms of a SU (2) L × U (1) Y invariant Lagrangian. The power series in physical fields and derivatives becomes an expansion in terms of canonical dimension for gauge-invariant operators. We expect a gauge-invariant Lagrangian to result from integrating out manifestly gauge-invariant new physics at higher scales, for an example, cf. Ref. [68]. If we truncate the gauge-invariant expansion at some fixed order, we obtain relations among the operator coefficients in (2.1). In the following subsection, we present the relations for the SILH parameterization [67] of the gauge-invariant EFT truncated at dimension six, cf. also Refs. [69][70][71][72]. Allowing for further higher-dimensional operators (dimension eight) in the EFT would relax those relations, so that the full phenomenological Lagrangian (2.1) is recovered with a different power-counting for the coefficients.

Relation to the SILH effective Lagrangian
There are various versions of the SILH effective-Lagrangian [67] parameterization. We refer to the following definition: Alternative parameterizations (operator bases) which are equivalent to the above Lagrangian at this order of the expansion, are obtained by applying linear or non-linear field redefinitions or adding total derivatives, and truncating the result again at dimension six. For instance, we may trade covariant derivatives of boson fields for flavor-universal contact terms if we apply the classical equations of motions, equivalent to specific non-linear field redefinitions: Boson-fermion contact terms yield subleading effects in dedicated VBF data analyses.
To reduce background, we apply cuts that enhance the quasi on-shell contribution for intermediate vector bosons, and we optimize the analysis for resonant final-state vector bosons. In this context, it is preferable to remove extra derivatives by the above relations. We also account for electroweak symmetry breaking and select the unitarity gauge. Finally, we can omit any contact terms and arrive at the phenomenological Higgs Lagrangian (2.1-2.4). The derivation, including conventions not listed here, is detailed in Appendix A.
We conclude this part with a remark on oblique corrections. In Ref. [73], the parameterŜ is defined byŜ where c W B /v 2 g g is the coefficient of the operator H † σ i HW i µν B µν /g g for non-canonical gauge fields, and c H is the coefficient of the operator |H † D µ H | 2 . If we translate the basis of Ref. [73] to our version of the SILH effective Lagrangian, we havê We recall thatŜ is constrained by data at the 10 −3 level; the precise value depends on the variation of the other electroweak parameterT . For our purposes, we have set c T to zero.

Relation to models of Higgs-inflation
Higgs-inflation models [74][75][76] provide an interesting example of a scenario where new physics is associated with the Higgs sector, with little impact on other SM particles. Such models are notoriously difficult to identify, and any possible probe of Higgs interactions should be investigated. We briefly review the derivation of the phenomenological Higgs Lagrangian for this model.

Higgs Inflation
c HW Consider the Higgs field coupled to gravity in a non-minimal way. The model is originally formulated as a Lagrangian in the Jordan frame, (2.20) The value of ξ can vary between 1 ξ 10 17 , corresponding to M M P .
For investigating the phenomenology, we apply the conformal transformation from the Jor-dan frame to the Einstein frameĝ This transformation leads to a non-minimal kinetic term for the Higgs field. In the unitary gauge H = 1 2 (0, h) T , we may introduce a scalar field χ as a transformed Higgs field, The action in the Einstein frame is whereR is calculated using the metricĝ µν . We neglect any renormalization-group running effect. The effective Higgs potential is In the context of collider physics, we are looking at small field values h χ and Ω 2 1, so the potential for the field χ is close to that of the initial Higgs field. Inflation physics is described by the large-field behavior of the Higgs field, the Higgs thus acting as an inflation, where h M P / ξ (or χ 6M P ). In this range, we can approximate The potential is exponentially flat at large h, as appropriate for a model of inflation. We are interested in a collider study and thus assume small h field values, so we replace χ by h again. We plug Eq. (2.22) into Eq. (2.23) and omit higher-order terms. Also re-instating the Higgs doublet notation H , we arrive at 1 Details regarding the gauge interaction can be found in Ref. [77]. We obtain corrections to the coefficients of the following operators: In Table. 2 we list the coefficient expressions for the Higgs inflation model and relate them to the SILH operator basis and to the Higgs Lagrangian that we use for our study. It is evident that the SILH operator basis, which is appropriate for a generic strongly interacting model, incorporates directions in parameter space which are absent in the more specific model of Higgs inflation. Obviously, specific measurements of Higgs self-interactions become essential if such a class of model is realized.

Constraints on parameters from the unitarity of S matrix
The importance of unitarity in quantum mechanics and quantum field theory has been emphasized since the very first papers on the subject. In the context of relativistic quantum field theory, unitarity is a constraint on the S operator as the most generic dynamic observable. While the basic idea, its formulation as the optical theorem, and its most elementary application to spinless 2 → 2 scattering are textbook knowledge, the concrete formulation for high-energy 2 → n scattering in the context of the SM and its extensions, is not as familiar. In this section, we provide an explicit introduction without the restriction to elastic scattering, and proceed to apply the generic formalism to the processes that we are interested in. Further details of the derivation are given in Appendix B. The goal is to determine energy-dependent constraints on the free parameters of the EFT.
An EFT, i.e., higher-dimensional operators in a local Lagrangian, is a method to generate and parameterize a low-energy expansion of the scattering amplitudes for unknown highenergy dynamics. A valid quantum field theory underlying any physics beyond the SM will result in a unitary scattering matrix, but this need not hold for the approximation generated by the corresponding EFT. In fact, a term of dimension d in the Lagrangian generically generates uncancelled factors that rise proportional to E d −4 , where E is the overall energy scale in a scattering process [78].
The unitarity requirement for the complete S matrix, evaluated for the set of m → n scattering amplitudes, relates the forward n → n elastic scattering amplitude to the interference of all m → n scattering amplitudes, integrated over phase space. In a simplified treatment where we neglect all n, m > 2, we can use angular-momentum conservation to simplify the scattering matrix to a finite N × N matrix which may be diagonalized in terms of partial-wave amplitudes, cf., e.g., Ref. [79]. This relation implies a strict upper bound on each partial wave, which can be exploited to derive energy-dependent constraints on the parameters of an EFT. For parameter values which violate those constraints, the EFT is invalid as a useful approximation, independent of the true underlying theory.
If we include 2 → n scattering with n > 2, we may apply similar methods and obtain comparatively simple constraints if we introduce extra assumptions [80,81] or neglect spins [82]. Extending these ideas, below we consider the generic formalism in its consequences for V V → hh, V V → hhh, and express the results in the form of inequalities for the EFT parameters.

General unitarity constraints
Unitarity is the conservation of probability in a quantum theory, applied to the S operator that encodes the scattering of observable particles: S † S = 1. Its nontrivial part T , defined by S = 1 + i T , thus satisfies the universal relation We are interested in unitarity conditions for matrix elements between asymptotic states which consist of a finite number n a of particles with well-defined masses. We denote multi-particle states collectively by |α, Φ a 〉, where Φ a is a shorthand for the kinematical configuration of n a on-shell four-momenta (the phase-space point), and α denotes the set of discrete quantum numbers such as helicity and color. Furthermore, we fix the total momentum of a multiparticle state a to p a . With this constraint, the manifold of configurations (α, Φ a ) becomes a compact manifold for each fixed n a . Momentum conservation allows us to introduce the matrix elements of the scattering amplitude operator M between the initial state |α, Φ a 〉 and the final state |β, Φ b 〉, We take matrix elements on both sides of Eq. (3.1) and insert a complete set of multi-particle states |γ, Φ c 〉, where dΦ c denotes the canonical Lorentz-invariant measure on the phase-space {Φ c } constrained by p c = p a = p b . For convenience, we may introduce a bijective mapping between the unit hypercube in d a = 3n a − 4 dimensions, {x a ∈ R d a ; 0 < (x a ) i < 1} and the manifold {Φ a }, for each fixed n a . For instance, we may factorize phase space as a tree consisting of n a − 1 momentum splittings of type 1 → 2, with p a at the root. There are 2(n a − 1) angular variables and n a − 2 invariant-mass variables. This mapping introduces a Jacobian J a (x a ) = dΦ a /dx a , which should incorporate symmetry factors where appropriate. This construction corresponds to a common method of evaluating phase-space integrals. If we introduce amplitude functions which include the Jacobian factors as follows, If massless particles are involved, the sum over intermediate states is infinite, and the matrix elements contain non-integrable infrared, collinear, and Coulomb singularities, so the integrals do not converge. To remedy this, we may introduce some version of phase-space slicing and sum over nearly degenerate states, which introduces indefinite particle numbers [83,84]. However, in the present context where we are studying the production of neutral massive bosons, with photons, quarks and gluons acting as spectators which are treated in standard QED and QCD perturbation theory, we may ignore this complication and assume all relevant states to be massive. The sum over intermediate states then is a finite sum, the matrix elements and the Jacobians are finite, and the integration manifold (the union of the unit hypercubes for all (n a , α) is compact.
In such a situation, it is possible to introduce a scalar product of square-integrable functions on the integration manifold and to find a complete basis of functions which are mutually orthonormal with respect to this scalar product. For instance, choosing the canonical scalar product, we could take a straightforward Fourier expansion. A more physical choice could involve spherical harmonics for the normalized angular variables and an arbitrary basis for the invariant-mass variables. In the two-particle case where there are no free invariant masses, this becomes the standard partial-wave expansion. We note that for each particle combination a, we may choose a different kind of expansion for the corresponding phase space Φ a (x a ).
We adopt, for simplicity, the canonical scalar product and a corresponding orthonormal where A is an appropriate (multi-)index which labels the basis functions. We expand the amplitudes as and thus express the scattering in terms of an actual matrix with elements a αβ AB . (The factor 2 has been inserted for consistency with the standard two-particle partial-wave expansion.) Explicitly, the coefficients are (3.8) They take complex values and depend only on the total momentum, a αβ AB = a αβ AB (p a ), where p a = p b . If we choose a phase-space parameterisation which preserves Lorentz invariance, the coefficients depend only on s = p 2 a . We obtain a discrete version of Eq. (3.5) [85,86], where all coefficients are finite and the sums are convergent if the simplifications regarding massless states are applied, as described above.
Eq. (3.9) encodes all unitarity relations of the scattering matrix in question. To derive constraints on individual amplitudes, we need a positivity condition. We may diagonalize the scattering matrix and obtain exact relations for superpositions of states. Alternatively, we may derive less comprehensive but phenomenologically more useful relations by focusing on diagonal matrix elements, i.e., α = β and A = B , To cast this in the intuitive geometry of the Argand circle, we may express the diagonal amplitude in terms of its real and imaginary parts and write That is, each complex-valued elastic amplitude a αα A A (s) must lie on a circle with radius r around i /2, where the elastic radius r = 1/2 is reduced by the total contribution of all inelastic channels.
The exact relation (3.12) yields strict upper bounds for the elastic amplitude as well as for the total inelastic contribution, which trivially translates into a bound for each individual final state in this representation. We read off Examples for the application of these bounds, referring also to the treatments in Refs. [79][80][81][82], can be found in Appendix B. The last inequality in Eq. (3.13) gives the unitarity constraints on inelastic scattering, and we note that it is independent of the basis H γ . To see this, we define the coefficients b αγ A as follows: which is clearly independent of H γ . Using the expansion in Eq. (3.7), we find The unitarity constraint for inelastic scattering can be written as

Unitarity Constraints from V V → hh
We now apply the generic formalism to the scattering process W + W − → hh. We assume that the on-shell approximation is justified for the purpose of deriving unitarity bounds, i.e., we can safely treat the incoming vector bosons as on-shell with a pair invariant mass M (W W ) = M (hh) =ŝ. (In the actual process, the incoming propagators are space-like with a virtuality of O(m W ).) We are thus looking at an inelastic channel. If we expand in a discrete basis as described in the preceding subsection, we deduce the bounds for the Higgs-pair production amplitudes where A and C are (multi-)indices for the initial-state and final-state basis, respectively. We note that the initial-state particles carry spin as well as momentum, while the final-state phase space manifold is trivially given by the unit sphere, for fixed energy ŝ. Figure 1. Four types of Feynman diagrams are shown which contribute to the processes W + W − → hh.
As shown in Fig. 1, there are four Feynman diagrams which contribute to the W + W − → hh process in the SM, and this breakdown remains valid in the EFT, (3.18) We refer to these as the s-channel, t -channel, u-channel, and contact-interaction amplitudes, respectively.
In the high energy limit s m 2 W , m 2 H , the leading contribution in the EFT is proportional to s. We thus write a series expansion as follows, in terms of the rescaled energy s/v 2 as a dimensionless expansion parameter, where m i are the coefficients in the expansion. In Table 3 we list the prefactors of the leading contribution for each amplitude and each one of the four independent helicity modes. The amplitudes of all other helicity modes are related to the four modes that we include in the table. We find that with the exception of one term (t /u-channel +−), all leading contributions are independent of the kinematics, so the table entries translate directly into bounds for amplitude coefficients once a suitable basis has been chosen. We also observe that only the ++, +−, and 00 modes lead to amplitudes rising proportional to s, so we may focus on those when considering unitarity bounds. Actually, the leading contribution to the +0 helicity amplitude is proportional to g W,a1 g W,b1 s/v 2 . If g W,a1 , which is constrained via the 00 mode, does not deviate grossly from its SM value, the +0 mode leads to a bound on g W,b1 which has the same s dependence but is weaker than the constraint that we get from the +− amplitude.
Angular-momentum conservation directs the choice of phase-space basis. The final state is described by a straightforward partial-wave expansion. For the initial state, we should couple helicity configuration helicity with orbital angular momentum to total angular momentum j , i.e., adopt the Wigner Dmatrix formalism (cf. Appendix B). We thus derive individual bounds for amplitude coefficients The strongest bounds on the EFT coefficients that we obtain for this process, are the following ones: In particular, the +− mode contributes a bound on g W,b1 , i.e., the hW + T W − T interaction, which is independent of the other EFT parameters. The helicity amplitudes of W + W − → hhh processes correspond to seven types of Feynman diagrams. Similar to W + W − → hh, in the high energy limit, the amplitude can be expanded as a series in powers of s/v 2 , We list the leading term m 0 in Table 4, for each helicity combination. Where the coefficient is phase-space dependent, we denote it as O (C ), where C is a combination of coupling constants. We find that the +0 helicity mode does not contribute to m 0 , and that the unitarity bounds resulting from the m 1 terms are weaker than the remaining ones, as long as g W,a1 , g W,a2 , κ 5 are not far from their respective SM values.

Unitarity Constraints from
Since this is also an inelastic channel, we obtain unitarity bounds on the b-coefficients defined in Eq. (3.14), Note that the b-coefficients are independent of the phase-space parameterisation and of the basis functions for the triple-Higgs system; only the phase-space parameterisation and the basis functions for the W -boson pair do matter. As discussed in Appendix B, after we choose the Wigner D-matrix as our basis for the W + W − state, the b-coefficients are diagonal and can be denoted as b j (h 1 h 2 ), where j represents the total angular momentum, and h i = + − 0 are the helicities of the two W bosons. We calculate the (reduced) b-coefficients directly according to Eq. (3.14). Although the result is independent of the phase-space parameterisation for the triple-Higgs system, an explicit expression is required for phase-space integration; we adopt the form given in App. B.4. We give the results for the three helicity modes below: 1. For the 00 helicity mode, the amplitude is independent of phase-space parameters. The optimal choice is given in Appendix B, and since the initial state is a two-particles state, the Wigner D-matrix as a basis yields this optimal bound, 2. For the ++ helicity mode, the type-f contribution is phase-space dependent, and it yields a non-zero b j for j > 0. However, we checked that the dependence is of minor importance, and the bounds from b j ≤ 1 4 with j > 0 turn out to be much weaker than the bounds (3. 27) with f 1 = 7.49994 ± 0.00005 and f 2 = 0.0658 ± 0.0006 computed by numerical integration. The negligible f 2 reflects the fact that the dependence of g 3 W,b1 on phase-space is small.
3. For the +− helicity mode, only b j with j = 2, 4, . . . are non-zero, and among them the largest one is b 2 , which is given by

Multi-Higgs production via VBF processes with dimension-6 operators
To evaluate the sensitivity of future colliders to new effects in multi-Higgs production, we compute the cross sections for the processes pp → hh j j and pp → hhh j j including the full dependence on the higher-dimensional operator coefficients, represented by the anomalous couplings in the phenomenological Lagrangian (2.1). To enhance the contribution of the VBF subprocess, we apply standard VBF cuts, as listed in Table ( Choosing M = 0 and g Shhh = −1 and restricting the calculation to the triple-Higgs production process, g w,a3 and g w,b3 become equivalent to the parameters in our convention, and the resulting amplitude expression is identical to the one that follows from using (2.1) directly, cf. Fig. 3. We have cross-checked numerical results from both implementations against each other. As another cross-check, we have validated selected results against the package VBFNLO [87,88], with good agreement. For the pure SM, we obtain the cross sections after VBF cuts as listed in Table 6. All numerical results are computed at leading order in the strong and electroweak perturbative expansions.
The VBF cuts in Table 5 force the remnant jets to a back-to-back configuration, with high energy and momentum, as it is expected from q → W q splitting in the VBF signal region. Due to the finite mass of the vector bosons, the VBF contribution is maximised for transversal momenta of the order of the W mass. We require p T ( j ) > 20 GeV for 14 and 27 TeV, and 30 GeV for 100 TeV, respectively. Regarding the transition from LHC kinematics to a 100 TeV collider, Figure 3. Triple-Higgs production diagram with a five-point vertex W W hhh effectively generated by an auxiliary field S.
Fraction of Events Fraction of Events > 500 GeV > 500 GeV > 800 GeV Table 5. Acceptance cuts used for the calculation of VBF Higgs production in pp collision (VBF cuts), for three different collider energies.
our numerical results demonstrate that the forward jets can acquire significantly larger rapidity than at lower energy (Fig. 4). Likewise, the produced Higgs bosons are distributed over a broader η range (Fig. 5). Therefore, we assume a better rapidity coverage for the detector at 100 TeV and have adapted our cuts in Table 5 accordingly.

Higgs pair production
For the on-shell process W + W − → hh, the effective Lagrangian (2.1) results in an expression of the form

4)
M u = U 1 g 2 W,a1 +U 2 g W,a1 g W,b1 +U 3 g 2 W,b1 , (4.5) where we have used the breakdown in terms of basic Feynman graphs as given in Fig. 1. The coefficients S i , T i , U i , and X i describe the expansion in terms of anomalous couplings. They include the SM parts and account for the polarisation wave functions of vector bosons, external momentum, propagators of internal particles, and the SM interactions. The full process pp → hh j j admits a similar expansion with new coefficients which incorporate the integration over the PDF of incoming partons and the complete phase space for off-shell intermediate vector bosons: Here, A i ,hh denotes second-order polynomials of the parameters g V,a1 , g V,b1 , λ 3 , κ 5 , g V,a2 and g V,b2 . The prefactors F hh i denote the integrated form factors which we compute numerically. We may simplify the general expression by making use of phenomenological information from expected precision data. For instance, the W W h vertex should be severely constrained by data from the Higgs decay to W W as well as VBF single-Higgs production. Higgs factories such as the CEPC, the ILC, or the CLIC collider allow for an absolute model-independent measurement of the W W h interaction. We account for this expectation and fix g W,a1 and g W,b1 , the Higgs couplings to longitudinal and W s, to their SM values, respectively. Furthermore, we assume the custodial-symmetry relations g W,a 2 = g Z ,a 2 and g W,b 2 = g Z ,b 2 whenever contributions of Z boson are considered. The resulting expression takes the form +K 6 κ 5 + K 7 g W,a2 κ 5 + K 8 g W,b2 κ 5 + K 9 κ 2 5 + K 10 λ 3 +K 11 g W,a2 λ 3 + K 12 g W,b2 λ 3 + K 13 κ 5 λ 3 + K 14 λ 2 3 (4.8) We note that the parameters λ 3 and κ 5 are accessible in the gluon-fusion process g g → hh, which may be measured with greater precision than Higgs pairs in VBF. In effect, VBF data will primarily constrain g W,a2 and g W,b2 , i.e., the anomalous contact interactions of a longitudinal or transversal W pair with a Higgs pair, respectively. The numerical calculation has been performed with Madgraph 5. We tabulate the results in terms of the coefficients K 0 −K 14 in Table 7. By using the cross section expression (4.8), we can derive bounds on the parameters which can be achieved at the LHC and at a future hadron collider. In Fig. 6, we display the dependence of the cross section on the two parameters g W,a2 and g W,b2 , respectively, for 14 TeV, 27 TeV, and for 100 TeV. One parameter is varied at a time. We mark the SM value, which in all cases is close to the minimum of the cross section. In the LHC (HE-LHC) plots, the horizontal lines indicate cross section values 10 fb and 30 fb (5/15 fb) which may be expected as realistic bounds if no signal is observed, respectively. In the 100 TeV plots, we show the bounds that result from constraining the cross section to ±10 %. (This is a conservative estimate; the study in Ref. [37] argues that a precision of 1% could be achieved.) We verify the expectation that the sensitivity to new physics that results from putting upper bounds on the cross section, increases considerably between the LHC and the HE-LHC. At 100 TeV, we expect an actual measurement of this process within the SM region.
We also show unitarity bounds for the parameter on the horizontal axis, for each plot. We have chosen the values 3 TeV, 5 TeV, and 14.5 TeV for the effective energy value that enters the inequalities (3. 26-3.27), to apply to a collider energy of 14 TeV, 27 TeV, and 100 TeV, respectively. We expect that for any realistic model, the cross section curves outside those bounds will flatten out in order to remain consistent with the optical theorem, such that constraining the cross section is not meaningful unless the experimental resolution reaches a certain threshold. Clearly, inserting a fixed cutoff value Λ when dealing with unitarity constraints for energy distributions can only yield a crude estimate of the EFT limitations. An appropriate framework for dealing with this problem has been described in Refs. [89][90][91].
In Fig. 7, we show projections for the correlated inclusion plots in terms of g W,a 2 -g W,b 2 , g W,a 2 -λ 3 , and g W,a 2 -κ 5 , respectively. (We recall that where numerical results are at the margin  of saturating unitarity limits, as discussed above, they should be interpreted with some care.) We observe that from the measurement of Higgs pairs in VBF alone, at the 14 TeV LHC we can determine g W,a 2 and g W,b 2 with similar precision, while the parameter λ 3 is only weakly constrained. The constraint on λ 3 is comparable with a recent parton level analysis [38], which shows that the HL-LHC can exclude the regions λ 3 > 3 and λ 3 < 0 by the VBF mode. Increasing the energy to 27 TeV and to 100 TeV, the precision improves as expected. The gain in sensitivity is particularly striking for g W,b 2 , the coupling to of a Higgs pair to transversal W gauge bosons. This property is evident from the huge value of K 5 in Table 7. Physically, the emission of transversally polarized W s from jets, which couple directly to the Higgs pair via the anomalous coupling g W,b 2 , is unsuppressed at asymptotic energy while the coupling to longitudinally polarized W s involves a W mass-mixing interaction and therefore becomes subleading.
The ultimate precision on the parameters clearly depends on the resolution power of the experimental analysis, which we do not consider in detail in this work. For a ±10 % measurement, we can read off the sensitivity on g W,a2 and g W,b2 from Fig. 7. If a boosted-Higgs analysis can constrain the cross section to ±3 fb and thus g W,a2 to 1 % precision [37], we analogously obtain a constraint on g W,b2 of the order ±0.008. This is within the range of the loop-induced interaction strength for this 'vertex.   . Projections for correlated bounds in the planes g a 2 -g b 2 , g a 2 -λ 3 , and g a 2 -κ 5 , from the process pp → hh j j at three different collider energies.
The numerical results clearly reflect the strong gauge cancellation which occurs between individual terms of (4.8) in the SM limit. Some of coefficients K i , such as K 1 , K 2 , K 3 , are one order of magnitude larger than the cross section in the SM. Furthermore, in our conventions, all linear terms, i.e. terms K 1 g W,a 2 , K 3 g W,b 3 , K 6 κ 5 , and K 10 λ 3 , are negative. The sign of these interference terms is retained when the collision energy increases from 14 TeV to 100 TeV.

Triple Higgs production
For the hhh j j final state, the breakdown of the amplitude in terms of Feynman graphs, Fig. 2, reads as follows: (4.10) M 3 = C 1 g W,a2 λ 3 +C 2 g W,a2 κ 5 +C 3 g W,b2 λ 3 +C 4 g W,b2 κ 5 (4.12) + E 3 g W,a2 g W,b1 + E 4 g W,b2 g W,b1 (4.14) For the full off-shell process, which involves the integration over the PDF of incoming partons, we can write Taking into account the results of the previous section, we may simplify the analysis. We can assume that g W,a1 , g W,a2 , λ 3 and g W,b1 , g W b 2 , κ 5 are known from the measurements of single-Higgs and double Higgs production to sufficient precision. We set them to their SM values 1 and 0, respectively, within this section. Keeping the remaining parameters, the cross section of the full process can be parameterized as σ(pp → V V j j → hhh j j ) = C 0 +C 1 g W,a3 +C 2 g 2 W,a3 +C 3 g W,b3 +C 4 g W,a3 g W,b3 +C 5 g 2 W,b3 + C 6 κ 6 +C 7 g W,a3 κ 6 +C 8 g W,b3 κ 6 +C 9 κ 2 6 +C 10 λ 4 + C 11 g W,a3 λ 4 +C 12 g W,b3 λ 4 +C 13 κ 6 λ 4 +C 14 λ 2 4 (4.18) We have computed the C i coefficients by numerical integration using WHIZARD. The results are listed in Table 8.
The numerical results enable us to derive projected bounds on the parameters. In Fig. 8, we display the dependence of the cross section on g W,a 3 and g W,b 3 , the direct coupling of longigudinal and transversal W bosons to three Higgs bosons, respectively. We only show the result for the 100 TeV collider. For reference, we include the unitarity bounds derived from Eqs. (3.26-3.27), evaluated for an effective energy of s = Λ 2 = (9 TeV) 2 . Since the cross section of the triple-Higgs process pp → hhh j j in the SM is very small (4.50 × 10 −2 fb), we do not expect data to provide a measurement, but rather an upper bound on the cross section, and thus exclusion limits for the anomalous couplings. If we assume an experimental sensitivity to a cross section of 10 fb, it turns out that in the 00 mode, the result already exceeds the unitarity bound for the chosen  effective energy. The perspective for a meaningful measurement of g W,a 3 alone in this channel becomes questionable. By contrast, for the ++ mode which probes the transversal coupling g W,b 3 , the measurement does provide a constraint. As discussed above, a more quantitative statement near the margin of unitarity saturation would require leaving the straightforward EFT approximation [89][90][91]. This is beyond the scope of the present paper.
In Fig. 9 we show EFT projections for correlated bounds in terms of g W,a 3 -g W,b 3 and g W,a 3λ 4 , respectively. We observe again that the coupling to transversal W bosons g W,b 3 , can be constrained more strongly than the other couplings. This is due to the fact that the value of C 5 in Table (8) is around 20 times larger than C 2 , and C 2 is around 4 times larger. As we have noted above for the double-Higgs case, the emission of transversal gauge bosons from the incoming partons at high energy is unsuppressed. In fact, comparing the 14 TeV LHC to a 100 TeV collider, the enhancement factors of the coefficients of C 5 and C 9 are 2×10 4 and 3×10 3 , respectively. The leading on-shell amplitudes grow proportional to s 3 , eqs. (3. 26-3.27). The terms that depend on λ 4 are subleading.
For a sensitivity limit of 10 fb at 100 TeV, we read off from Fig. 9c/d that the parameter g W,b 3 can be bounded to around 2%, and g W,a 3 to around 20%. This has to be understood with the caveat of unitarity constraints, as discussed above. The parameter λ 4 is also constrained, Fig. 9d, but this constraint is not expected to be improve on the measurement of g g → hhh. For 27 TeV, the expected bounds are accordingly weaker.

Multi-Higgs boson production with a strongly-interacting Higgs sector
In the preceding sections, we based the analysis on the hypothesis that all couplings that can be precisely determined, assume values close to their SM prediction. We could thus derive bounds  on further anomalous effects iteratively, starting from single-Higgs couplings and repeating the analysis up to triple-Higgs couplings.
One of the characteristic features of composite-Higgs models is the property that the Higgsvector boson couplings g W,a1 and g W,a2 as well as the Higgs-potential parameters λ 3 , λ 4 can deviate from the values of the SM by a sizable amount. In this section, we therefore consider a model where these couplings may take arbitrary values, but instead we keep the higher-order coefficients at their respective SM values, g a 3 = g b 1 = g b 2 = g b 3 = ξ 5 = ξ 6 = 0, to simplify the interpretation. (The parameters ξ 5 and ξ 6 can be eliminated in any case, cf. Eq. (2.5-2.10).) Measurements of multi-Higgs production are thus interpreted as supplementary data that improve on the Higgs-parameter determination from lower-order processes. This model corresponds to the truncated EFT in the linear representation (2.11), taken at face value without restricting the parameters to small values. In Sec. 2.2, we have introduced a Higgs-inflation model as a physics scenario which leads to an effective model of this class.
The parameterization of the tree-level cross section for pp → hh j j takes the form where D i can be computed numerically, and we have made use of the structure of the tree-level amplitude to limit the powers of the couplings that can appear. . Projections for correlated bounds in the planes g a 3 -g b 3 , g a 3 -λ 4 , from the process pp → hhh j j at two different collider energies.
Analogously, the cross section for pp → hhh j j is parameterized as below As before, the coefficients D i and T i are computed numerically using WHIZARD. The values of D i are presented in Table 9. The numerical values are large and rise rapidly with energy. However, when computing the complete amplitude, we observe a cancellation among D 0 , D 1 , and D 2 , independent of the value of λ 3 . Another cancellation occurs between D 3 and D 4 . The numerical values of T i are displayed in Table 10. Again, there are cancellations between terms.
In the SM limit, all terms proportional to positive powers of s have to vanish. To eliminate such partial cancellations from the beginning, we re-organize the calculation and redefine the couplings as δg 2 = g 2 g 2 1 − 1. We parameterize the cross section in the following form We do not assume the δg values to be small. The cancellations that we absorb by redefining the coefficients are numerically relevant: The primed coefficients are given in Table ( 11). Cancellations occur between (T 0 ,T 1 ,T 2 ) which yields T 0 , and between (T 1 ,T 2 ) which yields T 1 . At the energy of the LHC, the original coefficient T 2 is 8×10 2 larger than T 0 . At 100 TeV, the cancellation removes a factor of 4.5×10 4 , due to the s 3 enhancement of T 2 . Similar cancellations occur for (T 3 , T 4 , T 5 ), (T 6 , T 7 , T 8 ), (T 9 , T 10 ), (T 12 , T 13 ) and (T 14 , T 15 ). Furthermore, T 2 is 25 times larger than T 1 at 14 TeV, and 600 times larger than T 1 at 100 TeV.
Due to the s 3 dependence of the matrix element of the process pp → hhh j j given in Eq. (3.26) and Eq. (3.27), the enhancement of the ratio of leading coefficients T 0 − T 2 at a 100 TeV and the LHC 14TeV is around 3000, which is one order magnitude larger than that of D 0 − D 2 in the process pp → hh j j , as demonstrated in Tables 9 and 10, respectively. The coefficients T 3 − T 5 represent the next-to-leading contribution, which is one order of magnitude smaller than T 0 − T 2 .
This exercise demonstrates the exceptional behavior of the SM as a gauge theory. All positive powers of s in the amplitude cancel and disappear in the SM limit. The lesson is that for phenomenological estimates at ultra-high energies, it is important to split off the SM part in such amplitudes and cross sections, even if the deviation from the SM is not small.
As a final result for this model, in Fig. 10, we show the contours of constant cross section for pp → hh j j at the LHC and a 100 TeV collider in the g a 1 − g a 2 plane, respectively. Similarly, in Fig. 11 we show the exclusion regions in the g a 1 −g a 2 and g a 1 −λ 3 planes that can be derived from a cross-section measurement for pp → hhh j j . We note that there are two points in parameter space where the cross section is as small as in the SM: the first point corresponds to g a 1 ∼ 0;  Table 9. Coefficients D 0 − D 5 (in fb) in the expression (4.19) for the process pp → hh j j at three different collider energies.    Table 11. Redefined coefficients T 0 − T 17 (in fb) in the expression (4.21) for the process pp → hhh j j at two different collider energies.
the second point corresponds to the case g 2 a 1 − g a 2 ∼ 0. The first point is inconsistent with the observation of single-Higgs boson production in VBF, so only the second point is allowed.
For the specific case of the Higgs-inflation model of Sec. 2.2, we compare the ratio of cross sections of pp → hh j j and pp → hhh j j to their SM prediction at the LHC and a 100 TeV collider. The result may be expressed in terms ofx = 6ξ 2 v 2 /M 2 p , it is shown in Fig. 12.   Figure 10. Projections for correlated bounds in the planes g a 1 -g a 2 and g a 1 -λ 3 , from the process pp → hh j j with a strongly interacting Higgs sector at two different collider energies.

Obtaining and relating bounds on a generic Higgs sector
To put our results on multi-Higgs boson in their proper context, in this section we first briefly review the phenomenology of the dominant final states that need to be analyzed in an experiment. We also discuss the overall prospects for constraining the generic effective Lagrangian, comparing potential results from different channels at various future colliders.
The decay channels of a triple-Higgs state have been considered in Ref. [58]. In particular, the decay channels hhh → bbbbγγ and hhh → bbW W * W W * have been studied for the gluongluon fusion process [57,58,61]. For concreteness, we focus on the signal and background of the 6b, 4b2τ, and 4b2W final states, referring to SM values of branching ratios and cross sections.
The 6b final state, hhh → bbbbbb, has a branching ratio 20.30%. One of the main background for this final state is pp → htt , which can decay to 4b+jets. After applying the VBF cuts, we find that the σ × B R of this background is 3.48 × 10 4 ab. It is challenging to further reduce  Figure 11. Projections for correlated bounds in the planes g a 1 -g a 2 and g a 1 -λ 3 , from the process pp → hhh j j at a 100 TeV collider. Similarly, for the 4b2τ final state, the process pp → htt is also the main background, whose branching ratio is 7.16%. This can be suppressed by using 3b-tagging and τ-tagging techniques.
The channel hhh → bbbbW W * has the largest branching ratio [58]. The preferred decay channels of the W W * system are the semi-leptonic and fully-leptonic decays, whose branching ratios are 3.21% and 1.01%, respectively. Assuming that at a 100 TeV collider, an integrated luminosity of 30 ab −1 can be collected, we obtain 44 and 14 signal events, respectively. The dominant backgrounds of the hhh → bbbbW W * channel include pp → htt +jets, pp → Z tt +jets, pp → bbtt +jets, etc. VBF cuts reduce these backgrounds by two orders of magnitude, but they are still six orders of magnitude larger than the signal.
We conclude that in the SM, an unambiguous discovery of a triple-Higgs signal in the VBF mode remains a real challenge. In the presence of new physics, particularly if the Higgs sector is strongly interacting, cross sections can be enhanced by two orders of magnitude, so in that case the situation is more promising.
Considering the complete effective Lagrangian (2.1), we may summarize our results and results known from the literature as follows: • Constraints on L h The parameters of this part can only be constrained by measurements of multi-Higgs processes. Current LHC data can only give a weak bound λ 3 ∈ [−8.82, 15.04] [3]. Regarding the prospects at the HL-LHC, we refer to Ref. [24] which quotes a precision of 40% for this coupling; further studies are currently under way.
At the ILC, λ 3 can be measured via double Higgs-strahlung, e + e − → hhZ . Assuming a high-luminosity running scenario, this coupling can be constrained within accuracy 27% [92]. For the CLIC collider, Ref. [53] reports a precision of 19 % with a polarized electron beam for s = 1.4 TeV and 3 TeV combined, where the VBF production mechanism of Higgs pairs dominates.
At a lower energy e + e − collider, such as the CEPC, constraints on λ 3 can be derived indirectly via the loop corrections to the process e + e − → hZ [93,94]. Analogously, an indirect determination of λ 4 was considered in Ref. [54,95]. Direct constraints on this parameter can be derived from high-energy data on triple Higgs production at CLIC [54], where the enhancement of the cross section in the presence of anomalous couplings follows a similar pattern as for the pp processes that we consider in this work.
The current bounds on the remaining parameters are rather weak. At a 100 TeV machine, we expect that λ 3 and ξ 5 can be determined with a precision of about 10%. The parameters λ 4 and ξ 6 contribute only to triple-Higgs production, at leading order. In our previous study of triple-Higgs production in ggF, we obtained λ 4 ∈ [−13, 20] and ξ 6 ∈ [−2.3, 1.5] [61].

• Constraints on L V V h
Current LHC data strongly constrain the h → γγ decay channel. The measured signal strength is consistent with the SM prediction [96].
The couplings of a single Higgs to vector bosons can be measured in the decay process h → V V * . Recently, ATLAS presented results for h → Z Z * from an analysis of 13 TeV data [32]. Translating from the convention of Ref. [32], we have g Z ,b1 ∈ [0. 8, 4.7].
Bounds on g W,b1 as well first limits on the double-Higgs interactions g W,a2 ,g W,b2 , and g Z ,b2 can be inferred from current and upcoming LHC data. In the absence of a dedicated experimental analysis, we may use the projections for the bounds as given in Refs. [37,97]. Note that Ref. [97] used the SILH EFT basis where at the dimension-six level, we have g W,b1 = g W,b2 , g Z ,b1 = g Z ,b2 .
Lepton-collider data for the Higgs-production processes e + e − → hZ or e + e − → hνν give direct access to these parameters. For the coupling of the Higgs to transversal W bosons, the expectation is g W,b1 ∈ [−0.30, 0.30] [98]. The coupling to transversal Z bosons is naturally measured with much greater precision. For instance, the 240 GeV CEPC can achieve g Z ,b1 ∈ [−4.2, 4.2] × 10 −4 [99].
The parameter g X ,b1 is involved in the decay process h → Z γ. So far, LHC data do not show evidence for this process [100]; the upper limit on the cross section is 6.6 times the SM prediction. At the CEPC collider, the coupling hZ γ can be contrained by measuring e + e − → hγ [101] or measuring the angular observables of e + e − → hZ [99]. The projected bounds for CEPC are g X ,b1 ∈ [−8.4, 8.4] × 10 −4 .

• Constraints on L V h
Recent cross-section measurement of h → W W * and h → Z Z * can be found in [32,33]. The observed signal strengths are 1.05 ∼ 1.29 times the SM prediction, which gives a strong bound on the parameter g W,a1 .
At a lepton collider, the hV V couplings can be determined with high precision and independent of any model assumptions. Studies for CEPC and ILC imply that at both colliders, g W,a1 can be constrained down to 1% precision [92,98]. The projection for a global modelindependent of CLIC data yields similar results for achievable accuracy on the hV V couplings [53].
We summarize these values in Table 12. The new results from this work are included as projected bounds for the parameters g W,b 2 , g W,a 3 , and g W,b 3 , cf. Figs. 6 and 8.

Summary and Discussion
We have studied multi-Higgs production in VBF at the LHC and at future hadron colliders. While our emphasis lies on the rare triple-Higgs production mode, we have performed our calculation for a generic Higgs-sector effective Lagrangian and treated double-and triple-Higgs production processes in a common framework.
It is not surprising that even at a 100 TeV proton-proton collider, observing the multi-Higgs final state is very difficult and may not be possible at all in VBF if the SM is valid. Nevertheless, our study shows how such a measurement can nevertheless improve our overall knowledge about the interactions of the Higgs field. If all anomalous couplings are allowed to vary freely, constraining the direct couplings of three Higgs fields to themselves and to vector bosons is the logical next step after the analogous study of Higgs pairs. Anomalous effect in such couplings can spoil the delicate gauge cancellations of the SM, and therefore can lead to an increase in rates of several orders of magnitude. Turning this around, putting bounds on a rare process does yield meaningful constraints on parameter space, as shown explicitly in this paper.
Such large enhancement effects inevitably raise the issue of unitarity violation within the assumed model, which cannot occur in reality and therefore would limit the actual sensitivity to parameters. Applying the generic formalism of unitarity to multi-Higgs production processes, we have computed the energy-dependent region of parameter space where the model can still be considered as valid. We found that while couplings of multiple Higgs fields to transversal  , g W,b2 and g Z ,b2 , we take the projected bounds of Ref. [37,97]. ILC projections are taken from Ref. [92], which includes initial phase with 500 fb −1 at 250 GeV, 200 fb −1 at 350 GeV, and 500 fb −1 at 250 GeV, and a luminosity-upgraded phase with 3500 fb −1 at 500 GeV and 1500 fb −1 at 250 GeV. The CEPC projections are taken from Ref. [93,99], which assumes s = 240 GeV and 5 ab −1 integrated luminosity. The bounds for a 100 TeV hadron collider are obtained by assuming 30 ab −1 integrated luminosity, from Ref. [61] and from this work. vector bosons can be constrained within that region, constraints on couplings to longitudinally polarized vector fields become only marginally useful. In that case, for a detailed and quantitative assessment of the experimental resolution power, the straightforward effective-theory approach should be replaced by parameterized unitary models as a physically more reliable source of simulated data.
For an alternative interpretation of the same measurements, we have investigated two more specific models: the SILH effective Lagrangian imposes gauge-invariance constraints on the Higgs-electroweak interactions. We truncate the operator series expansion at dimension six but do not restrict the study to small values of the higher-dimensional coefficients. This Lagrangian does not constitute a consistent expansion if applied to the high-energy processes of this paper over the full kinematical range. Nevertheless, it provides a simple model framework for a coherent analysis of Higgs physics at colliders. An even more restricted model is obtained as the low-energy limit of a Higgs-inflation scenario, which as such recently has drawn much attention as a possible connection between Higgs physics and cosmology. For both models, we re-interpret our results on multi-Higgs production, and show appropriate projections for the achievable bounds in parameter space at future proton-proton colliders, namely at the proposed 27 TeV HE-LHC and at a future 100 TeV collider.

A Relating the SILH parameterization to the Higgs EFT Lagrangian
In our notation, the field strength tensors of the U (1) and SU (2) gauge groups are defined as respectively. The mass eigenstates of the gauge bosons are In unitary gauge, the Higgs doublet is given by By using the equation of motion of W µν and B µν (cf. [72]), we obtain the following expressions for the operators with coefficients c W and c B : To obtain expressions for the operators with coefficients c HW and c H B , we write the relations Here we used that . Also note that {D µ , D ν } vanishes when it is being contracted with an anti-symmetric tensor W µν , B µν . Eq. A.12 and Eq. A.13 correspond to the analogous relations in Ref. [71]. Expanding in unitary gauge, we obtain We arrive at the following kinetic part of the Lagrangian which includes the anomalous contributions, The fields may be rescaled by to rewrite the Lagrangian in terms of normalized fields as To eliminate the Z A mixing term, we introduce a linear shift as follows, This leads to In the final result, all electroweak gauge bosons are canonically normalized, and we may omit the primes from the redefined fields. The factors ζ h , ζ W , ζ Z ,ζ Z ζ A and ζ AZ are introduced for convenience.
From Eq. A.10 and Eq. A.11, we also have the mass terms There are shifts in the W mass and Z mass given by After rescaling the fields, we read off the parameter relations that are listed in Table. 2

B.1 2 → 2 scattering
The application of unitarity conditions to elastic 2 → 2 scattering is well known. In this subsection, for completeness, we provide the explict derivation and its connection to the generic formulas in Sec. 3. The derivation is not restricted to elastic scattering; it applies to any combination of two-particle initial and final states a and b, respectively. For a two-particle state vector |α, Φ a 〉, working in the center of mass frame, it is convenient to choose the polar angle θ a and azimuthal angle φ a as phase-space parameters, or correspondingly, the normalized kinematics variables are x a = ( 1 2 (cos θ a +1), φ a 2π ). The Jacobian determinant is given by where m a1 , m a2 are the masses of particles in a. S α is the symmetry factor that accounts for identical particles in a with quantum-number combination α: if the two particles are identical then S α = 2, otherwise S α = 1. Following Ref. [102], in the center of mass frame, the scattering matrix from the two-particle state |α, Φ a 〉 to |β, Φ b 〉 can be expressed as follows: 2 where θ a (θ b ) are the polar angles and φ a (φ b ) are the azimuthal angles for the states |α, Φ a 〉 (|β, Φ b 〉), respectively. ζ 1 , ζ 2 , ζ 3 denote corresponding Euler angles which represent the rotation from direction (θ a , φ a ) to direction (θ b , φ b ). We use the standard convention for parameterizing four-momenta in terms of polar and azimuthal angles, p µ = (E , | p| sin θ cos φ, | p| sin θ sin φ, | p| cos θ) (B.3) If the particle is a massive vector boson, we define the polarization states as follows: |p, +〉 = where m = E 2 − | p| 2 . This expansion suggests that we choose the Wigner D-matrix as an orthonormal basis for the 2-particle phase spaces, As a result, in the scattering amplitude between two-particle states the corresponding amplitude a becomes diagonal and depends only on one index: which by construction are independent of the phase-space parameterization pertaining to Φ c . We keep the dependency on the discrete quantum numbers γ of the intermediate state c. In analogy to the 2 → 2 case above, we may use any orthonormal basis for the initial twoparticle state a. Choosing the same Wigner D-matrix expansion is most convenient, however, since due to angular-momentum conservation the b coefficients only depend on one index, At this point, we may discuss the connection to previous literature on the subject [80][81][82].
• In Refs. [80,81], unitarity constraints are formulated for the total cross section of 2 → n scattering under the assumption that the j = 0 partial wave (s-wave) dominates. This assumption applies to some subset of the states that we consider here, but clearly is not justified for the generic case of polarized vector-boson scattering.
In fact, with our notation, the cross section for a → c with discrete quantum numbers α, γ is given by: where b j are the reduced b-coefficients after choosing the Wigner D-matrix as basis.
Assuming that the j = 0 partial wave dominates in the high-energy limit, we obtain which is equivalent to the result of Ref. [80,81]. This inequality applies to any polarized cross section and could provide a stronger bound than its equivalent for an unpolarized cross section.
• Ref. [82] considers the more generic situation of 2 → n scattering without s-wave dominance, but restricts the derivation to spin-less particles. In that case, the Wigner D-matrix formalism collapses to the familiar formalism of Legendre polynomials and spherical harmonics. By the general relation With these relations, it is easy to verify that our formulas reduce to the ones of Ref. [82] in the spin-less case.

B.3 Generalized s-wave
For some helicity combinations, the unitarity condition for 2 → n scattering becomes independent of phase-space parameters in the high-energy limit. This is the situation which was considered in Refs. [80,81]. In this subsection, we work out the details for our application. In the high-energy limit, effectively we can treat all external particles as massless, p 2 i = 0. The generalized s-wave condition for scattering a → c takes the form 〈γ, Φ c |M |α, Φ a 〉 ≈ C , (B.20) where C is a constant with respect to the kinematical parameters, for fixed total four-momentum. In fact, in the EFT approximation, this situation occurs naturally for some of the terms since the leading contributions become polynomials of the Lorentz invariants.
(a) For the case of inelastic scattering α = γ, the b-coefficients with (multi-)index A take the form where the total phase-space volume ∆ γ is given by [105,106] ∆ γ ≡   To realize the optimal bound within a given phase-space parameterization, the following condition should be satisfied: