Unitarization effects in EFT predictions of WZ scattering at the LHC

Effective field theories are an incredibly powerful tool in order to study and understand the true nature of the symmetry breaking sector dynamics of the Standard Model. However, they can suffer from some theoretical problems such as that of unitarity violation. Nevertheless, in order to interpret experimental data correctly a fully unitary prescription is needed. To this purpose, unitarization methods are addressed, but each of them leads to a different (unitary) prediction. Because of this, there is an inherent theoretical uncertainty in the determination of the effective field theory parameters due to the choice of one unitarization scheme. In this work, we quantify this uncertainty assuming a strongly interacting electroweak symmetry breaking sector, described by the effective electroweak chiral Lagrangian. We focus on the bosonic part of this effective Lagrangian and choose in particular the WZ scattering as our main VBS channel to study the sensitivity to new physics at the LHC. We study the different predictions of various well known unitarization methods, considering the full coupled system of helicity amplitudes, and construct the 95\% confidence level exclusion regions for the most relevant electroweak chiral Lagrangian parameters, given by the two anomalous quartic gauge couplings $a_4$ and $a_5$. This provides a consistent analysis of the different constraints on EChL parameters that can be achieved by using different unitarization methods in a combined way.


Introduction
Although the discovery of the Higgs boson by the ATLAS and CMS experiments supposed a great success of the Standard Model (SM), it also posed a lot of new questions about the symmetry breaking sector of the Electroweak (EW) Theory. Questions such as why is the Higgs boson so light, being its mass so similar to that of the EW gauge bosons, is the Higgs boson an elementary or a composite particle, what mechanism generates its potential, and others. All this can be summarized in the fact that the dynamical generation of electroweak symmetry breaking (EWSB) is still a mystery to be solved.
A very efficient way to try to understand the true nature of the EWSB sector of the SM is to use effective field theories (EFTs). Effective theories allow us to describe in a model independent way the relevant beyond the SM (BSM) physics that might be responsible for the dynamical generation of EWSB. In these theories the ultraviolet (UV) dynamics are, in principle, unknown, but their effects at low energies remain present encoded in a finite set of low energy parameters. If these low energy parameters were measured, we would have a hint towards the UV completion that might describe best the true dynamics of the EWSB sector.
With the aim of obtaining such a measurement, many experimental searches at the LHC are devoted to look for signals predicted by these effective theories. The most characteristic of these signals are those coming from vector boson scattering (VBS) processes (for recent reviews on VBS physics see, for instance, [1][2][3][4] and references therein) since it is there where the interactions among the longitudinal EW gauge bosons will appear dominantly, due to their relation to the scalar Goldstone boson modes which are associated to the spontaneous symmetry breaking taking place in the EW Theory.
Nevertheless, predictions of observable rates computed with effective theories can carry some theoretical problems. It is typical from such theories, due to the energy structure of the operators involved, to suffer from unitarity violation problems at high energies, such as the ones being now probed by the LHC. However, predictions that are to be tested at colliders must be fully unitary to be consistent with the underlying quantum field theory. Therefore, a prescription is needed to translate these non-unitary predictions into reliable, unitary ones with which interpret the experimental data. These prescriptions are called unitarization methods or procedures, that drive unitary the non-unitary EFT predictions (for some illustrative reviews on different unitarization methods in the context of VBS see, for instance, [1,5,6]). The problem with these methods is that the various manners of unitarizing the computation of an observable lead to different final results [1,[5][6][7][8][9][10][11][12][13][14][15]. Thus, a theoretical uncertainty arises when computing unitarized EFT predictions due to the fact that there is a variety of ways of achieving such a unitary outcome.
Current constraints imposed on some of the mentioned low-energy parameters by LHC experiments do not take this theoretical uncertainty into account. They are all basically based on searches for anomalous quartic gauge couplings that are then interpreted using the theoretical EFT predictions in different ways, i.e., using different unitarization methods or no unitarization method at all. For instance, the most recent constraints given in [16][17][18][19] provide a model independent experimental analysis, do not rely upon any particular method or employ one unitarization method only.
In this work, we quantify the uncertainty due to the choice of unitarization scheme present in the determination of some of the most relevant low-energy constants for VBS processes. To this aim, we assume a strongly interacting EWSB sector, properly described by the effective electroweak chiral Lagrangian (EChL) (nowdays called also Higgs effective field theory (HEFT)) (see, for instance, [20][21][22][23][24][25][26][27][28][29][30][31][32][33], before the Higgs boson discovery, and [8][9][10][11][12][13][14][15][34][35][36][37][38][39][40][41][42] after this discovery) and we focus on the bosonic sector of this Lagrangian, choosing to study the particular VBS process given by the WZ channel, as an example which is also interesting from the experimental detection perspective. Within this framework, we characterize the unitarity violation that arises in the predictions of the WZ → WZ cross sections, and we analyze the impact that a variety of well stablished unitarization methods have on them. We pay special attention to the fact that all helicity states of the incoming and outgoing gauge bosons might play a relevant role in the unitarization process, and consider them all at once as a coupled system. Then, we move on to the LHC scenario. We use the Effective W Approximation (EWA) [43,44] to give predictions of pp →WZ+X events at the LHC for different unitarization schemes. In order to check that the EWA works for our purpose here, we compare the EWA predictions for the cases of the SM and the EChL with the corresponding full results from the Monte Carlo MadGraph v5 (MG5) [45,46], and we find very good agreement in both cases. Finally, in order to provide a quantitative analysis of the implications of our study on the LHC searches, we choose in particular to compare our results with those in [16]. Concretely, we translate the ATLAS constraints from [16] to construct the 95% exclusion regions in some of the EChL parameter space for each of the considered unitarization methods, giving, at the same time, the total theoretical uncertainty driven by the variety of these methods, which represents our main conclussions in this work.
The paper is organized as follows. In section 2 we summarize the main features of the EChL and the relevant operators for describing deviations with respect to the SM in VBS processes. We also introduce the issue of unitarity violation in effective field theories. In section 3 we present the different methods we will use to deal with this problem and the corresponding predictions for the WZ → WZ scattering process within the EChL framework. We will also comment on the importance of taking into account the whole coupled system of helicity amplitudes, contrary to what is usually done in the literature. Section 4 is devoted to the presentation of our main results, in which we show the different predictions of the various unitarization methods in VBS processes at the LHC. In this section we display the impact of using different unitarization methods on the constraints that can be imposed on the relevant EChL parameters. The final section summarizes our main conclusions.

The Electroweak Chiral Lagrangian and the violation of unitarity
As already introduced in the previous section, we will work under the assumption of a strongly interacting EWSB sector. Even though the description of physics beyond the SM is unavoidably model dependent, we will employ a technology that is as model independent as possible. This technology is based on the effective theory parametrization of the possible BSM interactions, and the appropriate approach to a strongly interacting EWSB sector in this framework is to use the effective EChL, also known as HEFT in the literature. In this context, the information of the short-range theory is encoded in a certain number of coefficients of local operators, often called low-energy parameters.
The EChL is a gauged non-linear effective field theory based on the SU (2) L × SU (2) R chiral symmetry of the EWSB sector that is spontaneously broken down to the subgroup SU (2) L+R , usually called the custodial symmetry group or EW isospin group. It contains the EW gauge bosons, W ± , Z and γ, their corresponding would-be Goldstone bosons, w ± , z, and the Higgs scalar boson, H as dynamical fields, being the latter a singlet under the EW chiral symmetry. For the sake of simplicity, and since their contribution to VBS processes in the strongly interacting EWSB case should be negligible, we will not discuss the fermion sector in this article.
The w ± , z are introduced in a matrix field U (w ± , z) = 1 + iw a τ a /v + O(w 2 ) that takes values in the SU (2) L × SU (2) R /SU (2) L+R coset, and transforms as U → g L U g † R under the SU (2) L × SU (2) R global group. The EW gauge bosons are introduced through the following covariant derivative and field strength tensors: Here we will asume that custodial symmetry is preserved in the EWSB sector, except for the explicit breaking due to the gauging of the U (1) Y symmetry. This assumption is based on experimental measurements of the ρ parameter and of the effective couplings between the Higgs and the EW gauge bosons, that disfavor custodial breaking other than that induced from g = 0.
According to the usual counting rules within the chiral Lagrangian approach 1 , the different EChL operators are organized by means of their 'chiral dimension'. This chiral dimension can be found by following the scaling with the external momentum, p, of the various contributing basic functions, since, after all, the EChL structure is based on a momentum expansion. Derivatives and masses are considered as soft scales of the EFT and of the same order in the chiral counting, i.e., of O(p). Furthermore, in order to have a power counting consistent with the loop expansion, one needs ∂ µ , (gv) and (g v) ∼ O(p) or, equivalently, ∂ µ , m W , m Z ∼ O(p). The typical energy scale that controls the size of the various contributing terms in this chiral expansion is provided by 4πv, where v = 246 GeV is the vacuum expectation value of the Higgs field, in close analogy to the typical scale of 4πf π , with f π = 94 MeV, for the case of the chiral Lagrangian in QCD. In the scenarios where there are resonances that emerge typically from the assumed strongly interacting underlying UV theory, then there are additional mass scales given by the masses of the resonances to account for in the EChL. However, in this work we will assume that there are not emergent resonances below roughly 4πv ∼ 3 TeV and, therefore, this will be our unique energy scale parameter in the EW chiral expansion. All the other masses involved, m H , m W and m Z , are soft masses, as we have already said.
With all these considerations in mind, one then constructs the EChL up to a given order in the chiral expansion. This Lagrangian must be CP , Lorentz and SU (2) L × U (1) Y gauge invariant. Furthermore, with our simplifying assumption, it should be as well custodial symmetry preserving. For the present work, we include the terms with chiral dimension up to O(p 4 ), thus, the EChL can be generically written as: where L 2 refers to the terms with chiral dimension O(p 2 ), L 4 refers to the terms with chiral dimension O(p 4 ), and L GF and L FP are the gauge-fixing (GF) and the corresponding Fadeev-Popov (FP) terms. Nevertheless, not all the operators that can be included a priori in L 2 and L 4 have the same relevance in VBS processes. We conclude that for our purpose of describing the most relevant deviations from the SM in VBS it will be sufficient to work with just a subset of EChL operators, i.e., the most relevant ones. These operators are the following 2 : In this framework, as we have said, the chiral parameters a and b intervene in the interactions between two EW gauge bosons and one or two Higgs bosons, respectively. The other parameters, the ones controlling the strength of the O(p 4 ) operators, appear in the self interactions of the EW gauge bosons. The SM prediction is recovered for a = b = 1 and a i = 0.
The fact that in the context of this strongly interacting dynamics, operators, and thus, interactions among gauge bosons, scale directly with the external momentum, leads to a scenario in which predictions of observables can behave pathologically with energy from a certain energy scale upwards. This pathology translates into a violation of unitarity of the S matrix, which basically implies an unphysical leak in the interaction probability among EW gauge bosons. The energy at which this violation of unitarity occurs can be easily computed by studying the unitarization condition, implemented at the level of the partial waves, defined in this work as: where J is the total angular momentum of the system, λ = λ 1 − λ 2 , λ = λ 3 − λ 4 , being λ i the helicity states of the external gauge bosons, and where d J λ,λ (cos θ) are the Wigner functions. In this framework, the violation of unitarity occurs when the expression: is not fulfilled. Therefore, the violation of unitarity can take place for any value of the angular momentum J, and for every helicity channel, in principle. Furthermore, this expression implies that the unitarity condition of a particular helicity amplitude might depend on the amplitudes corresponding to other helicity channels too. This is an important statement and, therefore, when studying the possible unitarity violation of the various channels involved in WZ scattering, all helicity states should be considered consistently as a coupled system. Eq. (2.7) can be rewritten in a more friendly way in the following manner: This way, the value of the energy at which the J th partial wave crosses the unitarity limit, i.e., |a J (s)| = 1, defines the unitarity violation scale. With these considerations in mind, we now want to understand the relevance of each EChL parameter in the violation of unitarity. Due to the fact that every of these low energy parameters has a different role in the scattering of EW gauge bosons, each of them will have a different impact on this issue. In order to characterize the violation of unitarity induced by each of the EChL parameters presented in Eqs. (2.4) and (2.5), we compute the total cross section of WZ → WZ scattering in the EChL at the tree level for different representative values of one parameter at a time, setting the rest of them to their SM value. In Fig. 1 we show these cross sections as a function of the center of mass energy of the process and we mark the unitarity-violating predictions with dashed lines. The value of the energy at which each cross section overcomes the unitarity limit is chosen as the lowest one at which any of the corresponding J and/or helicity partial wave crosses the unitarity bound defined in Eq. (2.8). In these plots it can be clearly seen that in this scattering process the parameters a, a 1 and a 2 (upper left, upper right and middle left panels respectively) do not play a relevant role in the violation of unitarity, since in the whole energy range that has been studied in this work there is no unitarity violation driven from these coefficients. Notice that the b parameter, which controls the interaction between two EW gauge bosons and two Higgs bosons, does not appear in this scattering at tree level. When the parameter a 3 is considered (middle right panel), however, cross sections show a unitarity violating behavior in this same energy range. This happens only for large values of a 3 , of the order of 0.1, for which unitarity is violated at around 2 TeV. However, this size of 0.1 is already at the border of being in conflict with the EW precision data and, therefore, a realistic choice for this a 3 parameter should assume a smaller value than this, leading in consequence to non-violation of unitarity in the energy range relevant for VBS at the LHC. Overall, it is clear that a 4 and a 5 are the most relevant parameters regarding the issue of the violation of unitarity in this channel. If one takes a look  at the two lower panels of Fig. 1 it is manifest that for values of these two parameters between 0.1 and 10 −3 , the violation of unitarity occurs well inside the energy range considered in this work. Actually, for values of the order of 0.1, the crossing of the unitarity limit takes place at really low values of the energy √ s ∼ 800 GeV. Another interesting feature to pay attention to is that, in this particular channel, a 4 has a bigger impact in the cross section values than a 5 .
Based on these results, from now on we will consider only a 4 and a 5 , since their contribution to the violation of unitarity is, by far, the most relevant one in this context. This is indeed intuitive, as the effective operators related to these two parameters are the only ones that remain present in L 4 if the EW gauge interactions were switched off, i.e., in the limit g, g → 0 that corresponds to just keeping the self-interactions among the scalar modes. This is similar to the situation in chiral perturbation theory (ChPT) of low energy QCD if the electromagnetic gauge interactions were switched off by taking the limit e → 0 in the chiral Lagrangian which leads to just keeping the self-interactions among pions that provide the dominant contributions in pion-pion scattering. In summary, this means that within the EChL, if the underlying UV theory is strongly interacting, the dominant contributions from BSM physics to the scattering of longitudinal EW bosons are expected to be provided mainly by a 4 and a 5 . Thus, and although we of course work in the framework in which the electroweak interactions are still present (i.e. we consider g and g non-vanishing), relying on the assumption that the EW chiral coefficients a 4 and a 5 parameterize most relevantly the deviations from BSM physics in VBS, and specially in terms of the violation of unitary, is well justified.
At this point, it is important to make a comment on the experimental constraints imposed on these parameters. Nowadays, several experimental studies have been devoted to set bounds on the values of these and other EFT coefficients. Here, we will discuss those concerning a 4 and a 5 . For a summary on the bounds imposed upon other EChL parameters see, for instance, [42]. Regarding the a 4 and a 5 constraints, there is not a unique value for them given by the LHC experimental collaborations. Although all focus their searches on VBS observables, and on the search for anomalous quartic EW gauge boson coupling signals, they differ in the interpretation of the data to extract the bounds on the EFT parameters.
The most recent experimental searches at the LHC with √ s = 13 TeV aimed to constraint the parameter space of effective field theories for EWSB are explained in [17,18]. In [18] a maximum total cross section of various VBS processes, and, therefore, a model independent experimental study is reported, whereas in [17] direct bounds on the linear counterparts of some EChL parameters are provided. Concerning the results in [17], the translation to the a 4 and a 5 95% C.L. constraints corresponds to: (2.9) These are obtained without unitarizing the EChL (or, in those references the linear EFT 3 ) predictions at all, through a combined study of different VBS channels and analyzing the effect of each parameter at a time. One should keep in mind that these values for the a 4 and a 5 bounds might be overestimated, since the issue of the violation of unitarity has been neglected in the corresponding study. The same happens in [19], where only the WZ channel is considered in order to provide experimental bounds on a 4 and a 5 . We shall not directly use these bounds as reference, since they do not take into account the issue of the violation of unitarity, but we are planning to do it in a future work. Another interesting bound on a 4 and a 5 is the one provided in [16] for the LHC run at √ s = 8 TeV. There, a K-matrix unitarization analysis, following the procedure proposed in [7], is performed, and the EChL [a 4 , a 5 ] parameter space is constrained, as it is shown in Fig. 2, borrowed from [16]. We will rely mainly upon this experimental search of [16] as a first example in order to concrete quantitatively our final conclusions. Besides, as the overall constraints imposed in the EChL parameters in this study are of the order of a 4 ∼ a 5 ∼ 0.01, we will use these values as reference to illustrate different VBS features without loss of generality.
In this section we have defined and studied the violation of unitarity in WZ → WZ scattering in the context of the EChL. This violation of unitarity takes place at values of the energy of O(TeV) that are accesible at the LHC for different values of the most relevant parameters, a 4 and a 5 , which are the ones we will base our study on. The fact that unitarity is not fulfilled at those energies is not compatible with the correct interpretation of experimental data in order to test the EFT, since we need unitary predictions to be consistent with the underlying quantum field theory. To obtain these unitary predictions something must be done in the proper way. This is precisely to what the next section is devoted.

Restoring Unitarity in WZ → WZ scattering
In the previous section we have stated that the EChL, and specially the operators governed by a 4 and a 5 , lead to unitarity violation predictions for WZ → WZ scattering cross sections in the energy range accesible by the LHC. However, in order to make the EFT testable at colliders, we need to solve this problem and obtain fully unitary results for the relevant observables. To this aim, unitarization methods are addressed: prescriptions to construct unitary scattering amplitudes from the raw, non-unitary, EFT predictions. This is what we will do in this section, but, before entering in the specific details of these unitarization methods, some general considerations have to be commented.
First of all, it is important to have in mind that relying on a particular unitarization method for the EFT implies to make some assumptions about the UV complete theory. There is therefore a trade between obtaining unitary predictions for observables and losing some of the model independence inherent to EFTs. Nevertheless, there is a caveat in this statement.
When the EFT includes by construction the presence of resonant heavy states in the spectrum the various unitarization methods for VBS usualy provide comparable results, since the main features of the resonances (mass and width) are present in all cases. However, when the resonances are instead generated dynamically by the unitarization method itself (as it is the case of the Inverse Amplitude Method) this is not the case anymore and the results may vary substancially from one method to another one. Nevertheless, it is important to notice that if the unitarization methods provide amplitudes with the proper analytical structure, they can all accommodate dynamically generated resonances whose mass and width are predicted to be more or less the same independently of the employed method. Furthermore, in the different case of non-resonant scenarios, i.e., when there are not clear emergent peaks in VBS and one searches for smooth deviations from the SM continuum, different unitarization methods can lead to outstandingly different predictions for diverse observables. This suggests that, in order not to lose the appealing model independence of EFTs in the non-resonant case, the predictions given from the different unitarization methods available have to be contrasted, and a quantitative estimate of their differences should be provided. This inevitably introduces a theoretical uncertainty in the unitarized EFT predictions, which is precisely the one we want to quantify in this work. Therefore, we will focus in the case in which new resonant states do not manifest in the energies we are going to explore at the LHC via VBS. Besides, if present, they would also suppose a completely different experimental setup and search strategy. For a recent study of these emergent resonances at the LHC via WZ scattering within the EFT approach see, for instance, [42].
Second of all, if we recall the unitarity condition given in Eq. (2.7) that all unitarized amplitudes must fulfil, we see once again that the unitarity of a particular helicity channel does not depend just on itself but in other helicity amplitudes as well. This implies that considering only the most pathological of these amplitudes in terms of the violation of unitarity, could mean that we are neglecting important effects given the fact that the helicity system is coupled. In general, the most worrying helicity channel regarding the violation of unitarity will be the purely longitudinal one W L Z L → W L Z L . This is easily understood since the longitudinal modes of the EW gauge bosons are directly related to the strongly interacting Goldstone bosons. Thus, the bigger the number of longitudinally polarized gauge bosons involved in the scattering, the lower the energy at which unitarity will be violated. In any case, the fact that the purely longitudinal helicity channel dominates at high energies and dominates the violation of unitarity depends on the particular setup that one is considering.
By studying the partial waves with the lowest values of angular momenta, J = 0, 1, 2, for the 81 helicity channels independently and for different values of a 4 and a 5 and at different center of mass energies, one can disentangle the relevance of the purely longitudinal case with respect to the other helicity channels. These three lowest-order partial waves are the ones that should contain all the unitarity violating effects. This fact can be understood through the Equivalence Theorem (ET) that relates de EW gauge boson scattering amplitudes with the scalar scattering amplitudes at energies well above the EW scale. In the scattering involving just scalars in the external legs, since we have a polynomial expansion in energy up to order s 2 once we compute with the EChL at order O(p 4 ), all partial waves with J > 2 project to 0. Thus, the unitarity violation arising from the strongly interacting character of the interactions among scalars must be encoded in just the three mentioned partial waves even if we consider full gauge bosons in the external legs of our computations.
With this in mind, we have calculated the absolute value of the three lowest-order partial waves for all the helicity channels at a certain center of mass energy and for a particular value   Predictions are shown for a fixed center of mass energy of √ s = 1 TeV and for a 4 = a 5 = 0.01 (with the other parameters set to their SM value) as reference. Incoming and outgoing states can be interpreted indistinctly since the results are presented in a symmetric way due to time-reversal invariance. The included labels of these 9 incoming WZ and 9 outgoing WZ states with two polarized gauge bosons, longitudinal (L) and/or transverse (T +,− ), are ordered and denoted here correspondingly by: LL, of a 4 and a 5 , in order to understand the implication of the different helicity amplitudes in the total cross section and in the unitarization process. In Fig. 3 we present an example of this for the reference values of a 4 = a 5 = 0.01 and for a representative center of mass energy of 1 TeV. Looking at this figure, one can observe various interesting features. The first one is that, in general terms, the J = 0 partial wave (left panel) is around one order of magnitude bigger than the other two, J = 1, 2 (middle and right panels, respectively), as it is already well known in the literature. The second one is that only for that same value of the angular momentum, J = 0, the purely longitudinal scattering (displayed in the (1,1) entry of these "matrices", where incoming and outgoing states can be interpreted indistinctly since the results are presented in a symmetric way due to time-reversal invariance) dominates, being it a factor 5 larger than the next contributing helicity channel and thus becoming practically the only relevant amplitude to take into account. In the other two cases, J = 1, 2, the LL → LL case is no longer dominating the picture and other helicity channels become important. In particular we see in this figure that T + T + → T + T + and T − T − → T − T − play a relevant role in J = 1 and T + T − → T − T + and T − T + → T + T − do it in J = 2. This points towards the fact that, in some setups and for determined values of the relevant chiral parameters, neglecting the unitarity-violating effects of channels other than the purely longitudinal one could lead to incomplete predictions. This is the reason why we will consider the whole coupled system of the 81 helicity amplitudes when applying the mentioned unitarization methods.
The third comment one has to make regarding unitarization methods is somehow obvious, but important: all unitarization schemes have to provide similar predictions in the low energy region, i.e., above but not far from the WZ threshold production. This is a well known feature in the context of ChPT where the scattering amplitudes from the chiral Lagrangian, unitarized with the various methods, do recover the ChPT prediction at low energies, in agreement with the well known low-energy theorems.
Having stated all these considerations, we proceed to briefly explain the unitarization methods that we are going to consider in the present work. We have selected them based on the fact that they are the most used ones nowadays in the literature. Since these methods are the ones that are currently being used to interpret the experimental data in order to obtain information about the EFT, we find pertinent to contrast their predictions. They can be classified in two categories: 1) the ones that directly suppress by hand the pathological energy behavior of the amplitudes with energy (that we call here, as it is usual in the literature, Cut off, Form Factor and Kink), and 2) the ones that unitarize the first three partial waves from which then the total unitary amplitude is reconstructed (K-matrix and Inverse Amplitude Method (IAM)). Furthermore, they differ in their physical implications and motivation, and in their analytical properties, that we will discuss in the next paragraphs. Despite these differences, and the fact that some of them could be more physically justified than others, there is in principle no prior to choose a particular method with respect to the others.
We now list the five unitarization prescriptions considered in these work with a brief explanation of each of them (for some illustrative reviews on different unitarization methods in the context of VBS see, for instance, [1,5,6]) .
• Cut off: The Cut off is not a unitarization method per se but a way to obtain unitary amplitudes by just discarding those predictions given for energy values above the unitarity violation scale Λ, defined in the previous section as the lowest value of √ s at which any partial wave crosses the unitarity bound stated in Eq. (2.8). This would mean to reject the predictions of the cross sections marked with dashed lines in Fig. 1, sticking only to those that respect the unitarity condition (i.e., solid lines in these figures).
• Form Factor (FF): In this case, instead of obviating part of the results computed from the raw EFT, what is done is to suppress the pathological behavior of the amplitudes with energy above the scale at which each of them violate unitarity. To that purpose, a smooth, continuous function of the form: is employed. Here s is the center of mass energy squared, Λ i is the specific value of √ s at which the helicity channel i violates unitarity according to Eq. (2.8) and ξ i is the minimum exponent that is sufficient to fix the pathological behavior of the corresponding i th helicity amplitude with energy. Thus, every non-unitary helicity amplitude will be unitarized in the following manner: withÂ being the unitary amplitude and A the non-unitary EFT prediction. With all these unitarized amplitudes, then, one would be able to recover a unitary unpolarized, total cross section. In the present case, and for the values of the chiral parameters that are going to be probed in this work, the scales at which unitarity is violated for all helicity channels are above the maximum center of mass energy considered, except in the purely longitudinal case. We have checked that including the Form Factor suppression given in Eq. (3.2) for all helicity channels (notice that not only the scale is different in each channel, but also the exponent since they depend differently with energy) is equivalent to do it just in the LL → LL one for the energies and parameters we are considering, so, for simplicity, from now one we will apply Eq. (3.2) to the scattering of longitudinally polarized gauge bosons leaving the rest unchanged. In this way, our prescription to apply the Form Factor unitarization method can be summarized as: recalling that any other helicity amplitude is left unaffected. The exponent has been set to ξ LLLL = 2 since it is the minimum value necessary to repair the anomalous growth with energy of the LL → LL amplitude. The scale Λ LLLL has been computed with the VBFNLO utility to calculate Form Factors [47][48][49].
• Kink: The so called Kink unitarization method is very similar to the Form Factor. Conceptually, it is the same, and the only difference present between both prescriptions is that the suppression in the Kink method is not performed smoothly, but with a step function: Except for this fact, the rest of the discussion regarding the Form Factor is equally valid for the Kink, so, in this case, we will also apply the method only to the LL → LL amplitude with an exponent of ξ LLLL = 2.
• K-matrix: The K-matrix unitarization method has been extensively studied and implemented in the context of ChPT in QCD. This method is a prescription applied to the partial wave amplitudes and basically projects the non-unitary ones into the Argand circle through a stereographic projection. This means that it takes a real, non unitary partial wave amplitude to which an imaginary part is added ad hoc such that the unitarity limit is saturated. For each helicity partial wave amplitude, this is achieved by using the following simple formula: However, as we have already commented throughout the text, the unitarity condition implies that the whole coupled system of helicities has to be taken into account in our unitarization procedures. Thus, we solve this coupled system in terms of matrices, for which we construct a 9×9 matrix, whose entries correspond to the 81 possible helicity amplitudes of the elastic WZ scattering we are studying, and we unitarize it using the K-matrix method. This way we have: being α the 9×9 matrix containing the whole system of helicity partial wave amplitudes. Now, what we need is to reconstruct, from these unitary partial waves, the complete scattering amplitude. To this aim, we substitute from the initial, non-unitary amplitude, the unitarity violating partial waves by their unitarized versions. As we have already explained in the text, these partial waves are those that correspond to J = 0, 1, 2, so, what we do is to subtract these three partial waves from the total amplitude to then add the same partial waves after the K-matrix unitarization has been performed: Here we denote asα J;K−matrix [λ 1 λ 2 λ 3 λ 4 ] (s) (in the rest of the formulas it is implicit that all partial waves depend solely on s) the element of the 9×9 matrix that corresponds to the λ 1 λ 2 λ 3 λ 4 polarization state. In this way, we obtain a unitary amplitude in which we maintain all the fundamental properties introduced by all the partial wave amplitudes, including those with higher J > 2 that, since are not involved in the violation of unitarity, remain unaffected. The numerical computations in this K-matrix case and the next one, IAM, have been performed with a private mathematica code developed by us.
• Inverse Amplitude Method (IAM): The Inverse Amplitude Method is, probably, the most profoundly studied unitarization prescription considered in this work. It is very well known in the context of ChPT for pion-pion scattering, and its accuracy has been proved in various scenarios, like, for instance, in the prediction of the ρ meson as an emergent resonance in these scattering processes. It is based on the application of dispersion relations (bidirectional mathematical prescriptions allowing to relate the real and imaginary parts of complex functions) to the inverse of the partial wave amplitudes computed in the EFT framework. This unitarization procedure can be actually understood as the result of the first Padé approximant derived from the chiral expansion series provided by ChPT. In practice, this method implements an approximate re-summation of loops with bubbles in the s-channel of the given scattering process. Therefore in the present context of the EChL it accounts for re-scattering effects in the scattering of the two EW gauge bosons, i.e., WZ in our chosen example, which are not taken into account with the other unitarization methods. Notice that this makes sense in the context of a strongly interacting theory since these re-scattering contributions are not suppressed as in weakly interacting systems.
In summary, if one starts with the typical result for a given partial wave amplitude from the chiral Lagrangian, given by the sum of the two contributions in the chiral expansion, one of order O(p 2 ) and the other one of order O(p 4 ) , the corresponding prediction of the IAM leads to the following unitarized helicity partial wave amplitudes: where a (2) is the contribution to the partial wave amplitude computed with the operators from the L 2 Lagrangian (Eq. (2.4)) at the tree level, which is of order O(p 2 ) and a (4) is the contribution to the partial wave amplitude computed with the operators from the L 4 (Eq. (2.5)) at the tree level, plus the contribution computed with the operators from the L 2 Lagrangian at one loop level, which are both of O(p 4 ). In the present work, since the computation of the complete one loop level amplitudes that enter in a (4) has not been performed yet due to the difficulty of the task, we will evaluate this here in an approximate way. Following the usual features in ChPT, we take the imaginary part of this contribution to be a (2) 2 so that the unitarity condition is fulfilled perturbatively, and neglect the real contribution of the loops which are expected to provide a very small contribution, not being relevant for the present computation.
Once again we encounter ourselves in the scenario in which we have a prescription to unitarize each helicity amplitude independently. However, we want to take the whole coupled system of helicities in full generality, as explained above. We construct once more the 9×9 matrix α, this time splitting it into its O(p 2 ) and O(p 4 ) contributions, that contains the 81 helicity amplitudes, and we unitarize it using the IAM in the following matricial manner: At this point, to obtain a fully unitary amplitude, we use the same trick as in the K-matrix case, i.e., we replace the unitarity violating partial waves of the total amplitude by their IAM-unitarized version, following Eq. (3.7) with the only change of K-matrix→IAM. It is pertinent to make now some comments regarding important differences between the IAM unitarization method and the rest we have considered. The IAM does not provide just unitary predictions, but also succeeds to get partial wave amplitudes with the appropriate analytical structure (for more details on this, see, for instance, [6]). This implies that it is the only method, amongst the ones studied in this work, that can accommodate dynamically generated resonances, since these appear as complex poles in the second Riemann sheet of the partial wave with the corresponding J quantum number. This is in contrast to the unitarized partial waves with the K-matrix method that do not have such poles. These resonances are characteristic of strongly interacting theories, and appear naturally at high energies, such as in the case of low-energy QCD. Furthermore, it is worth commenting that, according to [11], similar results as those obtained with the IAM regarding the appearance of dynamical resonances are also provided by other alternative unitarization methods that lead to the proper analytical structure. Example of such methods are the N/D or the improved K-matrix, which for shortness, we have decided not to include here. Nevertheless, for the forthcoming study at the LHC, as we have already said, we are interested in studying the non-resonant case of the unitarized theory, so the differences among the various unitarization methods will come in terms of smooth deviations from the SM continuum via WZ scattering rather than from the appearance of peaks due to the emergence of resonances. It is important, though, to keep in mind that the IAM has some peculiarities regarding its structure and physical motivation, that differentiates it from the others.
We have already discussed briefly each of the unitarization procedures that we consider in this work and the specific way in which we implement them. Now, what we need is to study the different predictions they provide in regard of VBS observables. To that purpose, we apply each of them, in the above explained manner, to the WZ→WZ scattering amplitudes for different values of a 4 and a 5 . The final results of this computation can be seen in Fig. 4, where we present the total, unpolarized cross section of the WZ→WZ scattering as a function of the center of mass energy for the different unitarization methods used. We also display the SM prediction and the non-unitarized EChL prediction for comparison. We consider,  A great number of interesting facts can be extracted from these plots. Firstly, and most clearly, one can see that different unitarization methods lead, indeed, to very different predictions for this observable. These predictions differ, in general, from those of the raw EChL and the SM, as well. Therefore, one can expect that these differences might be very well seen experimentally. Secondly, another interesting issue arises from the comparison of both panels in the figure. It appears plainly that the value of the EChL prediction and of the K-matrix prediction are, in general, smaller in the case in which both parameters, a 4 and a 5 , have opposite sign. On the contrary, the Kink and the Form Factor provide larger results in this same case. This means that for the non-unitarized prediction and for the K-matrix, the regions of the parameter space in which a 4 and a 5 have the same sign will be more constrained, whereas for the Kink and the Form Factor the opposite-sided regions will be the most constrained ones. We have checked that the predictions for the scenario in which both parameters are negative gives the same results as the one in which they are both positive, and that in the case in which they have opposite sign the same result is obtained when either of the parameters is positive/negative. Thirdly, a comment has to be made regarding the Cut off procedure. The unitarity violation scale is not explicitly shown in these plots, but it can be inferred from the position of the "knee" in the Kink prediction. As it is clear, discarding the values of the cross section above this scale will imply to lose a lot of sensitivity, and will of course correspond to a very different prediction with respect to the other studied cases.
Regarding the IAM, we can clearly see that for the particular choice of parameters in the left panel of Fig. 4 its prediction lies very close to the SM one. In this case the IAM does not provide an emergent resonance in WZ scattering, since for these particular values of the EChL parameters there are not poles in the reconstructed total amplitude from the IAM partial waves and, as a consequence, the outcome provided by the IAM when applied to the LHC context will not show any departure from the SM continuum. In contrast, for the particular choice of parameters in the right panel, there is indeed an emergent resonance below 1 TeV, which we have decided not to include in this plot in the right panel since it is most probably already excluded by the present searches at the LHC.
When other particular values of the EChL parameters are chosen, different patterns in the predictions of the VBS cross sections from the various unitarization methods can appear. In general, the choice of smaller values of |a 4 | and |a 5 | than those in Fig. 4 typically lead, in the non-resonant case, to closer predictions for the various unitarization methods in the studied energy range, and also closer to the SM prediction. This can be clearly seen in the upper left pannel in Fig. 5 where the parameters have been set to a 4 = a 5 = 0.0001 and a = 0.9 (or equivalently, ∆a = a − 1 = −0.1). For this particular choice, a scalar resonance emerges close to 3 TeV in the IAM unitarized predictions, which does not manifest in the channel of our interest here W Z → W Z but in the W W → ZZ channel. This can be seen clearly in the plot of the upper right panel in Fig. 5, which we have included for comparison. In this case, studying this alternative VBS channel W W → ZZ at the LHC seems more appropriate in order to analyze the distortions with respect to the SM predictions due to BSM physics represented by this particular choice of parameters.
The other example included in Fig. 5, where the parameters are set to a 4 = 0.0004, a 5 = −0.0001 and again a = 0.9, displays the emergence of a vector resonance in the IAM prediction for W Z → W Z (lower left panel) close to 2500 GeV, and a scalar resonance close to 2800 GeV in the IAM prediction for W W → ZZ (lower right panel). This resonant behavior is only found in the predictions with the IAM but not in the predictions with the other unitarization methods.
In summary, regarding the IAM, the appearance of dynamically generated resonances in the energy range of a few TeV occurs indeed for a continuum set of a 4 and a 5 values of the order of O(10 −3 −10 −4 ) and its properties, mass and width, also depend on the other relevant parameters, particularly on a. These features of the IAM have been studied extensively in the literature and are not the main focus of the present paper which, as we have said, is mainly devoted to the non-resonant case. Thus, for the rest of this work we will focus on the other unitarization methods which will produce instead smooth distortions from the SM continuum.
Finally, a significant point has to be made concerning the K-matrix cross sections, since they are the ones we will use in the next section as a link to the experimental results. We have compared our estimates, obtained with the K-matrix procedure explained in the pages above, with the ones provided by the Wizard group [5,7]. In the given references, the authors construct unitarized four point functions that can be introduced in a Monte Carlo event generator. Their prescription is based in the T-matrix unitarization method, that they implement in a similar way than us: replacing the unitarity violating partial wave amplitudes of the total amplitude by their T-matrix unitarized version 4 . This prescription is used, actually, by the ATLAS collaboration in order to constrain the EChL parameter space [16]. Nevertheless, their work is based on the ET, and they unitarize all the helicity amplitudes using the ET calculation, valid only to describe the longitudinally polarized gauge bosons at high energies. Thus, given this difference between their method and ours, we consider pertinent to make some comments about the discrepancies we have found.
Our predictions match those of the Wizard group for all the LL → LL amplitudes we have considered, i.e., for all the studied energies and values of the chiral parameters. helicity channels differ. In the case in which the purely longitudinal scattering dominates at high energies, both procedures give rise to the same values for the cross sections. If other helicity channels have important contributions to the total cross section, we obtain different predictions. This can be the case if the values of a 4 and a 5 are very small, of the order of, for instance, 10 −4 . The authors in [5,7] themselves comment on the limitations of their approach in this regime, so we are proposing here a way to avoid these limitations. We have seen that different unitarization methods lead to very different predictions for the values of the cross section of the elastic WZ scattering. For this reason one can expect that the translation of these results to the LHC scenario would also show the different behaviors present at the subprocess level. Precisely because of this, the experimental measurements and constraints interpreted using one method or another will be different, and this difference can be understood as a theoretical uncertainty which is precisely the one that we want to quantify in this work. Thus, in the next section, we will present our results for the LHC, and we will give an estimate of this uncertainty in the experimental determination of a 4 and a 5 due to the unitarization scheme choice.

Parameter determination uncertainties at the LHC due to unitarization scheme choice
In the previous section we learnt that the predictions for WZ scattering observables computed in the EChL framework can be very different depending on the unitarization method we apply to them. This was manifest at the subprocess level, but now we want to study and quantify these deviations as they would be seen at the LHC. In order to compute the total cross section at the LHC we have first used the simple tool provided by the Effective W Approximation [43,44] and then we have compared this approximate result with the full result from MG5 [45,46]. The EWA is the translation to the massive EW gauge bosons case of the familiar Weiszäcker-Williams, or EPA, approximation for photons [50,51]. This framework has two important advantages: the first one is that has the intuitive physical interpretation of the distribution functions of the W and the Z as the PDFs in the parton model, and the second one is that it is computationally simple and, as we will see next, leads to very good results within the context we are working on here. The EWA provides probability functions, f W,Z (x), for the W and the Z that describe the probability of the EW gauge boson to be radiated collinearly from a fermion carrying a fraction x of its total momentum. In order to get the total cross section at the LHC for the full process that starts with protons, these functions, taking quarks as the mentioned fermions, are then convoluted with the PDFs of the quarks 5 and with the corresponding subprocess cross section for the scattering of on-shell EW gauge bosons, σ(W Z → W Z) in our present case. Furthermore, the computation with the EWA requires to separate the different polarizacions for the EW gauge bosons, and to use accordingly the corresponding probability function for the polarized W or Z. For the numerical computation of the cross section at the LHC with the EWA in this work we have developed our own private PYTHON code.
There are several studies in the literature that use the EWA to obtain reliable estimates. However, not all of them employ the same probability functions. For this work, we have considered and compared four of these implementations of the EWA. These four implementations are: 1) the original EWA functions given in [43], including first the Leading Log Approximation (LLA) ones (eqs. 2.19 and 2.29 in [43]); and second 2) the improved ones which go beyond the LLA by keeping O(M 2 V /E 2 ) corrections, with M V the EW gauge boson mass and E the energy of the initial quark (eqs. 2.18 and 2.28 in [43]); 3) the EWA functions derived from [44]; and 4) the simplified functions of the beyond LLA given in [7]. In principle, all should lead to similar results for the pp →WZ+X process, and they do at high invariant masses of the final diboson system. Nevertheless, they differ quite a lot at lower energies. It is worth mentioning that to compute the pp →WZ+X rates with the EWA, one has to consider the contributions from two different subprocesses: the intermediate state with a W and a Z radiated from the initial protons that then scatter, and, in addition, the case in which a W and a photon are radiated and then scatter. The latter is of great importance in the low energy region where it dominates indeed over the other one. For the photon case we have used the well established probability function of the Weiszäcker-Williams approximation [50,51]. In order to select the most accurate probability function for the EW gauge boson case among the ones available in the literature, we have compared the results of the above four mentioned approaches to the full results for the complete process pp →WZ+X obtained using MG5. Notice that for this comparison we have generated MG5 events of the exclusive process pp →WZjj, which automatically contain all the topologies, i.e., the VBS topologies and all the others contributing to the same order in perturbation theory. Besides, in order to compare properly both results, the MG5 one and the EWA one, one has to set particular kinematical cuts on the final state particles. In particular, as it is well known, in order to regularize the Coulomb singularity produced by the diagrams with a photon interchanged in a t-channel, some minimal cuts have to be imposed on the final particles. Concretely, for this quantitative comparison of the total cross sections we give the following cuts on the transverse momentum and pseudorapidity of the final gauge bosons V and jets j, and the angular separation among the jets: |p T j | > 5 GeV; |η j | < 10 ; ∆R jj > 0.1 , (4.1) both in the EWA and MG5 for the cuts concerning the gauge bosons, and in MG5 events only for the ones concerning the extra jets. With these considerations in mind, we have compared quantitatively the predictions of the pp →WZ+X processes in the SM and in the EChL for a 4 = a 5 = 0.01 within the EWA for the four probability functions considered against the MG5 computation of the pp →WZjj events. From our numerical comparison (not included here for shortness) we have reached the conclusion that the original, improved probability functions given in [43] are the ones that better match the MG5 prediction. The others overestimate the probability of radiating a EW gauge boson at low fractions of momentum of the initial quarks, thus missing the correct prediction of the cross section at low energies where most of events lie. In Fig. 6, we display the results of the differential cross section distribution with respect to the invariant mass of the final gauge bosons, computed in the SM (green) and in the EChL (gray) for a 4 = a 5 = 0.01 using the EWA and employing these improved probability functions. We also show the MG5 prediction for these same distributions as solid, darker lines of each corresponding color, as well as the total cross sections obtained with both procedures. Cuts in Eq. (4.1) have been required, if applicable, and center of mass energy has been set to √ s = 14 TeV, as it will be considered for the rest of the work. Regarding the comparison shown in this figure, it is manifest that the EWA works remarkably well, specially at high invariant masses. Not only the total MG5 cross section is recovered within a factor 1.5 at the worst in the SM case and 1.15 in the EChL case, but also the invariant mass distributions match considerably well. Now that we have checked that our computations obtained with the EWA employing the improved probability functions provide reliable predictions of pp →WZ+X observables, we move on to characterize the behavior of the different unitarization methods at the LHC. To that purpose, we have convoluted the subprocess cross sections of each of the studied unitarization methods, corresponding to the different curves in Fig. 4, with the EW gauge bosons probability functions and with the CT10 set of PDFs [52], evaluated at Q 2 = M 2 W . The results are displayed in Fig. 7, where we present the invariant mass distributions of the differential cross section of the process pp →WZ+X computed with the EChL for a 4 = a 5 = 0.01 (left) and a 4 = −a 5 = 0.01 (right) and unitarized with the diverse procedures we have described in the previous section. The non-unitarized EChL and the SM predictions are also shown, for comparison. The unitarity violation scale is marked with a dashed line in each case. The final gauge bosons are required to have |η V | < 2 and |p T V | > 20 GeV and the evaluation is performed at √ s = 14 TeV. From these curves we can see that the translation of the subprocess results to the LHC is direct, and the conclusions regarding the results are very similar. The different predictions among the various unitarization methods are still manifest, which clearly indicates that the experimental constraints imposed on the EChL parameters will strongly depend on the unitarization method used to analyze the data. Besides, the same pattern of the predictions concerning the relative sign of the chiral parameters is encountered: in the EChL and the K-matrix case, same sign a 4 and a 5 lead to larger predictions than in the opposite sign case. For the Form Factor and the Kink, the reverse setup is recovered. This still points towards the fact that same sign values of a 4 and a 5 will be more constrained in the EChL and the K-matrix case, opposite to the Form Factor and the Kink case. The IAM is not shown in these plots since, as we mentioned, it is more suitable for the resonant case. Besides, as we have seen before, in the present non-resonant case, for the chosen particular channel W Z → W Z, and with the simplified setup of just two non-vanishing chiral coefficients, a 4 and a 5 , the IAM predictions are very close to the SM ones. Notice that it will not be the case if other channels were considered (for instance, we have checked this explicitly for W W → ZZ) and other chiral coefficients (in particular, we have checked this for a = 0.9, and |a 4 |, |a 5 | ∼ O(10 −3 − 10 −4 )) were also non-vanishing. Regarding the Cut off, it is clear that integrating only up to the unitarity violation scale to obtain the total cross section will lead to much smaller predictions than in the rest of the cases. Finally, it is worth commenting that, as it should be, again all predictions match the EChL one at low invariant masses.
We have now characterized the different predictions of the studied unitarization methods at the LHC. The next step should be to translate these predictions into uncertainties in the extracted constraints on the parameter space of the EChL. In order to do that, we will base our approach upon the ATLAS results for √ s = 8 TeV given in ref. [16]. In the mentioned reference, a very sophisticated experimental analysis is performed, specially regarding triggers, background estimations, and event selection. Then, using the K-matrix (or T-matrix) unitarization prescription proposed in [5,7] , the 95% C.L. exclusion regions in the [a 4 , a 5 ] (sometimes called [α 4 , α 5 ] in the literature) parameter space are obtained. It is beyond the scope of this work to reproduce accurately the experimental analysis of the ATLAS searches. However, there is a consistent way in which we can use their results to obtain the experimental constraints corresponding to other unitarization methods apart from the K-matrix one. Our approach is the following: first, we take the a 4 and a 5 values lying on the contour of the WZ observed "elipse" provided by the ATLAS study. With those values, we evaluate the total cross section following our K-matrix unitarization procedure for the LHC case, that is, indeed, constant over the mentioned values. This should be equivalent to what ATLAS has performed, since we have checked that for these values of the parameters our prescription matches the one given by the Wizard group. The cross section that we obtain represents the equivalent cross section in our framework to the one that ATLAS has measured experimentally. It is, so to say, a translation between the experimental results and our naive results. Now, what we do is to find the values of a 4 and a 5 that lead to the same cross section for the other unitarization methods considered in the present work. In this way, we construct the 95% exclusion regions in the [a 4 , a 5 ] plane for the various unitarization schemes presented in the previous section, to see how they differ in magnitude and shape. By applying this procedure, we are assuming that the selection cuts required to be fulfilled by the ATLAS search affect all our predictions equivalently. This could not be the case, but we expect the differences to be small, so our prescription should be a good first approximation to the issue. Furthermore, it is worth commenting that, regarding the backgrounds, since they are the same to all of our signals, it is well justified to proceed in this way.
The final results of the present work, i.e., the 95% C.L. exclusion regions in the [a 4 , a 5 ] plane for different unitarization scheme choices, are presented in Fig. 8. There we show the corresponding limits for the case in which no unitarization is performed at all (EChL, in gray), for the K-matrix unitarization, matching, of course, the ATLAS results (purple), for the Form Factor prediction (blue) and for the Kink (yellow). We also show the total exclusion region, obtained by the overlap of the former ones.
Many interesting features can be extracted from this figure. First of all, and most importantly, it is indeed very clear that using one unitarization method at a time to interpret experimental data does not consider the full EFT picture. Since there are many unitarization prescriptions that lead to very different constraints, one should take them all into account in order to provide a reliable bound on the EFT parameters. These different constraints can vary even in an order of magnitude, as it is manifest in Fig. 8. For instance, the Form Factor prescription leads to bounds on [a 4 , a 5 ] of the order of [0.8,0.4], roughly speaking, whereas the case in which there is no unitarization performed leads to constraints of the order of [0.04,0.08]. Notice that these latter bounds, i.e., those obtained from the raw, non-unitarized EChL, are not directly comparable to those given in [19], since our results correspond to √ s = 8 TeV and the ones reported in the mentioned reference to √ s = 13 TeV. We leave the precise computation of the 13 TeV results for a future work. It is also obvious from this figure that the Kink method leads to more stringent constraints than the FF method (the corresponding pseudo-ellipse is smaller and oriented similarly to the FF one). Also the Kmatrix method leads to more stringent constraints than the FF one and, in this case, with a different orientation of the pseudo-ellipse (indeed, similar to the EChL one). Interestingly, for the present studied case of non-resonant pp → WZ +X events, there is not just a difference in the magnitudes of the bounds of the EChL parameters, but also in the role of each of them, a 4 and a 5 . This feature was already stated before, since in Fig. 7 we had already seen that points lying in the region of the plane in which a 4 and a 5 have the same sign should be more constrained in the EChL and the K-matrix case, just in the opposite direction to the Form Factor and the Kink case. At this point, two further comments have to be made. The first concerns the IAM, whose prediction is not present in this figure. It is due to the same argument we have been commenting throughout the text, that can be summarized in the fact that, for our particular setup (non-resonant case with deviation with respect to the SM coming only from the two considered O(p 4 ) operators, i.e., for a = b = 1 and just a 4 and a 5 non-vanishing) this method is not suitable to impose reliable constraints on the EChL parameters. Nevertheless, the IAM can be extremely useful when looking for new physics signals at the LHC in the resonant case, as it has been studied in [42]. The second concerns the Cut off, also not present in the figure. Since this procedure implies to sum events only up a to a determined invariant mass of the diboson system to obtain the total cross section, a problem arises concerning the backgrounds. In our approach, we are always integrating over the whole studied energy region, for all the unitarization method predictions. This means that the background is considered to be the same for all of our signals and we can use the translation from the ATLAS results safely. However, if we now change the picture and integrate over a smaller invariant mass region, such as in the Cut off procedure, we should take into account this same integration over the background, and the pure translation form the ATLAS results fails, since we don't know the background scaling with energy. For this reason, we have not included the Cut off prediction in our final results, but we really do believe that it should be also considered in proper experimental searches.

Conclusions
It is undoubtable that effective field theories constitute a remarkable, model independent tool to help us understand the true nature of the electroweak symmetry breaking sector. They typically suffer, however, from unitarity violation problems, due to the energy structure of the operators they contain. If the predictions of observables in such theories violate unitarity from some energy scale upwards, they are, in principle, not compatible with the underlying quantum field theory. Therefore, reliable, unitary predictions are needed to interpret experimental data in order to obtain information about the effective theory and thus about the dynamics it describes. With the aim of obtaining these predictions, unitarization methods are addressed. There are, nevertheless, many available options to drive non-unitary observables computed with the raw effective theory into unitary ones. This ambiguity supposes then a theoretical uncertainty that has to be taken into account when constraining the parameter space of such theories, that, up to now, has been made using one of these prescriptions at a time or no prescription at all. In this work, we provide a first approximation to quantifying this uncertainty by studying the non-resonant case of the elastic WZ scattering at the LHC.
In order to do so, we use the electroweak chiral Lagrangian or Higgs effective field theory, that is the most approriate EFT if one assumes a strongly interacting EWSB system. With this EChL, we study the violation of unitarity in the WZ→WZ scattering and select the most relevant operators concerning it. These correspond to the ones controlled by the chiral coefficients a 4 and a 5 , which parametrize in this context the anomalous interactions among four massive electroweak gauge bosons. Furthermore, at this point, we also analyze the relevance in the unitarization procedure of each of the helicity channels participating in the scattering. Although nowadays most unitarization studies assume that the purely longitudinal scattering is sufficient to understand the unitarity violation processes, we consider of much importance to take into account the whole coupled system of helicity states, since the unitarity conditions relates them. Therefore, we implement this coupled system analysis in our final results.
Once we have characterized our framework, we study the predictions of different unitarization methods at the subprocess level. We choose five of these methods: Cut off, Form Factor, Kink, K-matrix and Inverse Amplitude Method, which are the most used ones in the literature, to illustrate how different their predictions can be. Their concrete implementation as well as a brief explanation of each of them can be found in Section 3. When analyzing each of these methods' predictions at the subprocess level, compared to those of the raw effective theory and of the SM, one is convinced that they lead to very different results, and that this fact should be taken into account at the time to impose constraints on the parameter space of the effective theory.
Moving on to the LHC case, we use the Effective W Approximation to give estimates of the predictions of the various unitarization methods considered for the pp → WZ+X process. In order to be sure that the EWA works properly for our purpose here, we have first compared its predictions at the LHC of the total cross sections and differential cross sections with the invariant mass M W Z , both in the SM and in the EChL case, with the corresponding full predictions provided by MadGraph. In this comparison, we study various probability functions available in the literature for the massive electroweak gauge bosons and select the ones that better reproduce the MadGraph simulation of the total pp → WZ+X process. Concretely, we find that the improved EWA functions in [43] are the most accurate providing predictions which are in very good agreement with the MadGraph result. Afterwards, we employ these most accurate EWA functions to obtain the predictions of the invariant mass M ZW distributions of the differential cross section of the pp → WZ+X events for the different unitarization methods discussed in this work. We conclude again that the various unitarization methods provide very different predictions not only for the subprocess but also for the total process at the LHC.
Finally, we construct, based on the ATLAS results for √ s = 8 TeV given in [16], the 95% exclusion regions in the [a 4 , a 5 ] plane for the various unitarization schemes. The main results of the work are contained in Fig. 8, from which very interesting features can be extracted. The most important of them is that it is indeed very clear that using one unitarization method at a time to interpret experimental data does not consider the full effective theory picture. Since there are many unitarization prescriptions that lead to very different constraints, one should take them all into account in order to provide a reliable bound on the EFT parameters. These different constraints can vary even in an order of magnitude. As an example, the Form Factor method leads to bounds on [a 4 , a 5 ] of the order of ∼ [0.8, 0.4] whereas the pure EChL prediction, without unitarization, leads to constraints of the order of ∼ [0.04, 0.08]. Furthermore, the differences do not lie just in the magnitudes of the bounds, but in the role of a 4 and a 5 , what can be seen in the shapes of the different exclusion regions.
The main conclusion of this work, is, therefore, that there is a theoretical uncertainty present in the experimental determination of effective theory parameters due to the unitarization scheme choice. A first approximation to this uncertainty has been quantified in the present work analyzing the predictions of pp → WZ+X events at the LHC from the EChL in terms of a 4 and a 5 and with different unitarization methods. We believe that it is important to take these uncertainties into account when relying upon experimental values of the constraints of effective theory parameters, in order to consider the full effective theory properties correctly.