A Golden Probe of Nonlinear Higgs Dynamics

The most salient generic feature of a composite Higgs boson resides in the nonlinearity of its dynamics, which arises from degenerate vacua associated with the pseudo-Nambu-Goldstone (PNGB) nature of the Higgs boson. It has been shown that the nonlinear Higgs dynamics is universal in the IR and controlled only by a single parameter $f$, the decay constant of the PNGB Higgs. In this work we perform a fit, for the first time, to Wilson coefficients of ${\cal O}(p^4)$ operators in the nonlinear Lagrangian using the golden H $\to$ 4L decay channel. By utilizing both the"rate"information in the signal strength and the"shape"information in the fully differential spectra, we provide limits on the Goldstone decay constant $f$, as well as ${\cal O}(p^4)$ Wilson coefficients, using Run 2 data at the LHC. In rate measurements alone, the golden channel prefers a negative $\xi =v^2/f^2$ corresponding to a non-compact coset structure. Including the shape information, we identify regions of parameter space where current LHC constraint on $f$ is still weak, allowing for $\xi \lesssim 0.5$ or $\xi \gtrsim -0.5$. We also comment on future sensitivity at the high-luminosity upgrade of the LHC which could allow for simultaneous fits to multiple Wilson coefficients.

The Higgs boson plays a central role in our understanding of physics at the electroweak scale, the energy scale being probed by the Large Hadron Collider (LHC). Precise measurement of its property is among the top priorities of experimental programs at the LHC as well as any possible future electron-positron and hadron colliders. A major goal of these efforts is to understand the microscopic nature of the 125 GeV Higgs boson: • Is the Higgs boson a fundamental particle like the electron or a composite particle such as the pion?
Such a question could carry far reaching implications in our understanding of the Universe at the most fundamental scale, especially given the fact that no other fundamental scalar particles have been observed in nature. While current data at the LHC is consistent with the expectation of a Standard Model (SM) Higgs boson, the experimental uncertainty is still sizeable, on the order of 10 -20% or more [1][2][3]. Such a large uncertainty can hardly be characterized as "precise" and the question of whether the 125 GeV Higgs boson is indeed the SM Higgs remains open.
There is a long history in the idea of a composite Higgs boson, dating back to the classic papers [4,5] several decades ago, where the Higgs boson arises as a pseudo-Nambu-Goldstone boson (PNGB) like the pion in low-energy QCD. Models with a PNGB Higgs were revived and refined much later, via the little Higgs theories [6][7][8] and the holographic Higgs models [9,10]. By now they are collectively referred to as the composite Higgs models.
There are numerous composite Higgs models [11] which differ in several aspects. For example, the choice of symmetry breaking pattern G/H where G is the broken group in the UV and H is the unbroken group containing the electroweak SU (2) L × U (1) Y gauge groups.
The Higgs arises as a PNGB when G is spontaneously broken to H at a scale Λ ∼ 10 TeV. One other aspect that varies greatly is the implementation of additional fermions associated with the third generation quarks in the SM, which are introduced to reduce the UV sensitivity in the Higgs mass originating from the SM top quark. In some cases the new fermions do not even have to carry SM color charge [12] or electroweak quantum numbers [13,14].
In spite of all these variations in model-building, there is one salient prediction that is generic to the entire class of composite Higgs models: • Nonlinear dynamics in Higgs interactions with the electroweak gauge boson, as well as self-interactions carrying derivatives.
The origin of nonlinear Higgs interactions goes to the essence of a composite Higgs: the PNGB nature of the Higgs boson. It is the same nonlinear dynamics appearing in interactions of pions in chiral symmetry breaking, or any other Nambu-Goldstone bosons observed in nature. The theoretical tool employed to construct effective Lagrangians of Nambu-Goldstone bosons was developed half a century ago, in the seminal papers by Callan, Coleman, Wess and Zumino (CCWZ) [15,16]. In the CCWZ approach each symmetry breaking pattern G/H results in a seemingly different effective Lagrangian, each with its own set of nonlinear interactions and experimental predictions.
Only in recent years was it realized that the nonlinear interaction of a PNGB has very little to do with the details of the UV group G that is being spontaneously broken. This can be seen either by imposing shift symmetries on the Nambu-Goldstone bosons in the IR [17,18], or by soft bootstrapping tree-level amplitudes of nonlinear sigma model (nlσm) [19,20]. More specifically, the nonlinear interaction of Nambu-Goldstone bosons owes its presence to the existence of degenerate vacua in the deep IR and knows very little about the details of the broken Group G [21]. The only information on G resides in the overall normalization of the Goldstone decay constant f , which can be taken as an input parameter in the low-energy.
Following this progress, the complete list of modifications to the couplings of one Higgs boson to two electroweak gauge bosons (HVV), two Higgs bosons with two electroweak gauge bosons (HHVV), one Higgs coupled to three electroweak gauge bosons (HVVV), as well as triple gauge boson couplings (TGC), were studied in Refs. [22,23] up to four-derivative order and all orders in 1/f . These corrections are the universal predictions of a composite Higgs boson. In particular, a set of "universal relations" among the couplings were proposed which depend on only one input parameter f . Experimental confirmation of these universal relations would be a striking signal of the PNGB nature of the 125 GeV Higgs.
In this work we continue with the exploration of universal Higgs nonlinearity and its implications in Higgs coupling measurements. In particular, we will study the possibility of measuring and constraining Wilson coefficients of four-derivative operators in the effective Lagrangian of a composite Higgs boson using the H → 4L decay channel, the so-called "Golden channel" because it is the clearest and cleanest among all Higgs decay channels.
The small background contamination combined with the availability of full kinematic distributions of decay products offers not only a powerful probe of the spin and CP property of the Higgs boson [24,25], but also a unique opportunity to employ advanced multivariate techniques to enhance the experimental sensitivity [26][27][28][29][30][31][32][33][34]. We will use the publicly available LHC analyses in the 4L channel to probe and constrain, for the first time, Higgs nonlinear dynamics and the associated Wilson coefficients. Furthermore, we demonstrate that the conventional approach of bounding the decay constant f using the signal strength (total rate) in HVV measurements is incomplete, especially when effects of O(p 4 ) operators are included.
This work is organized as the follows. In Section II we briefly review the universal Higgs nonlinearity, listing all operators modifying the HVV couplings up to O(p 4 ). We also map out the correspondence to the tensor structure basis used in the LHC analyses. Then in Section the H→ 4L channel, as well as in H→ Zγ two body decays. Future sensitivity projections at the high-luminosity (HL) LHC are also provided in this section as well as a brief study of precision electroweak constraints. Finally we conclude in Section IV.

II. UNIVERSAL HIGGS NONLINEARITY
The effective Lagrangian of a PNGB Higgs boson contains two expansion parameters, • Nonlinear expansion characterized by π/f , where π denotes a generic Nambu-Goldstone field like a composite Higgs boson. The expansion in π/f is highly nonlinear for a PNGB Higgs boson.
• Derivative expansion chararcterized by ∂ µ /Λ. If some of the unbroken symmetry is gauged, one replaces the ordinary derivative by the gauge covariant derivative, In this power counting the gauge field A µ then counts as one derivative.
This is similar to the chiral Lagrangian in low-energy QCD, which describes interactions of pions as PNGB's arising from spontaneously broken SU (2) L × SU (2) R chiral symmetry [35].
In the chiral Lagrangian f = f π ∼ 93 MeV is the pion decay constant and Λ = Λ χSB ∼ 1 GeV is the scale where QCD becomes strongly coupled and chiral symmetry is spontaneously broken. In naive dimensional analysis (NDA) [36], requiring quantum corrections from higher loops to be comparable in size to the derivative expansion in ∂/Λ leads to the relation which is saturated in a strongly interacting theories with no relevant small adjustable parameters. Then the effective Lagrangian of Nambu-Goldstone bosons based on a nonlinear sigma model (nlσm) is organized as follows where L (n) represent operators containing n derivatives. In this work we focus on n ≤ 4. While L (n) incorporates the derivative expansion at a well-defined order, the nonlinear expansion is resummed to all orders in π/f in L (n) . Conventional wisdom from the CCWZ approach [15,16] has it that the nonlinear expansion realizes the spontaneously broken symmetries in the UV. The summation of operators to all orders in π/f and only at a certain order in ∂/Λ stems from requiring the resummed operators to have well-defined transformation properties under the broken symmetry. The nonlinear Higgs dynamics then follows. However, the CCWZ perspective obscures the universality in the Nambu-Goldstone interactions. It turns out that these nonlinear interactions arise entirely from the presence of degenerate vacua in the IR and are insensitive to the coset structure G/H in the UV [17,18,21,22]. The only parameter dependent on the coset structure is the normalization of f . In the end, by probing and measuring the Higgs nonlinear dynamics, one could potentially gain insight on a broad range of composite Higgs models, regardless of the G/H coset.
Invoking the custodial invariance and focusing on the scenario where the 125 GeV Higgs transforms as a fundamental representation of an unbroken SO(4) group, the leading twoderivative Lagrangian is simple, upon gauging an SU (2) L ×U (1) Y subgroup of SO(4) [22,23], where h is the 125 GeV Higgs boson. In the above sin θ ≡ v/f , where v = 246 GeV, is the vacuum misalignment angle and W ± and Z are the electroweak massive gauge bosons, whose masses can be read off from Eq. (3): where ξ ≡ v 2 /f 2 . Given that f is an input parameter whose normalization is UV dependent, Eq. (5) is the basis for extracting f using the signal strength in HVV coupling measurements. above the scale M ρ [37]. This effective theory is often referred to as the SILH Lagrangian, which has the same nonlinear structure as in Eq. (2), but with Λ → M ρ . In particular, at the two-derivative level, Eq. (3) remains unchanged even after integrating out M ρ , and corrections to HVV and HHVV couplings are also unchanged. There are effects that could potentially spoil the nonlinear structure of the effective Lagrangian. These are related to "explicit" symmetry breaking effects and examples include the Higgs potential and the Higgs coupling to fermions. However, as argued in Ref. [22], they would modify the nonlinearity HVV and HHVV couplings only at the loop-level. The typical spectrum in a composite Higgs model is displayed in Fig. 1.

A. O(p 4 ) Operators
The two-derivative nonlinear Lagrangian in Eq. (3) is valid when where E is the typical energy scale probed by the experiment. In particular, resumming contributions to all orders in v/f allows for an f not too far above the weak scale set by v ≈ 246 GeV. We can further extend the range of validity of the effective action by including higher order terms in the derivative expansion, such as O(p 4 ) operators. Chiral Lagrangian operators that are O(p 4 ) and to all orders in π/f have been enumerated in Ref. [38]. In the context of composite Higgs bosons the relevant four-derivative operators for Higgs couplings were given in Refs. [22,23] in the unitary gauge.
More specifically, the SILH Lagrangian at O(p 4 ), including the full nonlinearity structure, can be written as where c i are expected to be order unity constants parameterizing the incalculable UV physics at the scale Λ ∼ 4πf . In some cases operators contributing to couplings of neutral particles and an on-shell photon are further suppressed by additional loop factors. In total there are 7 , each with the corresponding unknown Wilson coefficients c i . These 7 Wilson coefficients contribute to a dozen different observables that can be measured in HVV and HHVV couplings [22,23]. Focusing on those relevant for HVV couplings we have, in the unitary gauge, where D = ∂ µ ∂ ν − η µν ∂ 2 . Although there are six coefficients C h i , i = 1, · · · , 6 in the unitary gauge, they are secretly related by only five Wilson coefficients residing in Eq. (8), as well as the input parameter defined by sin θ = √ ξ = v/f , in a composite Higgs model: In the above c w , c 2w , t w denote cos θ W , cos 2θ W , tan θ W respectively, where θ W is the weak mixing angle.

B. H→ 4L Tensor Structure
The operators listed in Eq. (9) enter into H→ 4L decays and manifest themselves through kinematic distributions of the decay product. For our analysis we utilize the "Golden Channel" analysis framework developed in Refs. [26, 27, 29-31, 33, 34] which parametrizes the effective Higgs boson couplings to pairs of neutral vector boson pairs in terms of the Lorentz tensor structures, where V = (ZZ, Zγ, γZ, γγ). We assume massless leptons, using the notation and conventions defined in [31], but have also included the A V 4 tensor structures which were not included previously. Note that electromagnetic gauge invariance requires A γγ In general the A V i can be momentum dependent form factors which are functions of Lorentz invariant products of the external momenta. For our purpose it is sufficient to take them to be real and constant coefficients in which case we have A 4 =Ā 4 . Furthermore, in the For simplicity we have kept only CP-even terms and neglected effective couplings to pairs of photons in Eq. (16), which do not appear in universal Higgs nonlinearity, though it would be straightforward to include them in the fit. Note that as ξ → 0 we recover the SM value for the tree level HZZ coupling, A ZZ 1 = 2. The important observation is: is directly related to ξ; the other tensor structures are dependent on both At the leading order in derivative expansion, the signal strength in HVV coupling measurements is attributed entirely to A ZZ 1 , which is then used to constrain ξ [39]. We see this is not the case anymore when O(p 4 ) effects are included. It turns out that all couplings in Eq. (17) can be extracted from the fully differential spectra of H→ 4L decays. In what follows we use CMS H→ 4L data to constrain the A V i coefficients, which then translate into limits on ξ and the Wilson coefficients C h i .

III. EXPERIMENTAL CONSTRAINTS
A. The Golden Channel: H → 4L In this subsection we will use both the "rate information," which pertains to the signal strength measured in H→ 4L channel, and the "shape information," as contained in the differential spectra of final state leptons, to place constraints on the Wilson coefficients in the nonlinear Higgs Lagrangian in Eq. (9). In particular, we utilize the fully differential decay width analytically computed in [27,29] for H → 2e2µ, 4e, and 4µ assuming on-shell decay of the Higgs boson. Interference effects between the different tensor structures in Eq. (16), as well as among identical final states in the case of 4e and 4µ, have been fully accounted for. Before examining current CMS constraints, we briefly review the H→ 4L partial decay width and fully differential spectra.
For our purpose it is convenient to single out the A V i dependence in the differential decay width, which can be written schematically as where O represents all observables available in the H→ 4L decay channel [27,29] and i, j sum over all of the possible tensor structures in Eq. (16). The advantage of doing so is then the remaining quantities dΓ ij /dO can be calculated and integrated over a particular phase space. We can then define the 'sub-widths' for each combination of which is just a numerical constant once a selection cut over the phase space is chosen. It is worth emphasizing that the sub-widths could be negative for certain combinations of tensor structures when they interfere destructively. However, the partial width written as the sum must be positive and is now a function of the effective couplings and the phase space cuts.
Eq. (18) and Eq. (20) allow for a full reconstruction of the differential decay spectra and the partial width of H→ 4L, respectively.
It will be convenient to normalize Γ H→4L to the (tree level) SM expectation, which corresponds to A ZZ 1 = 2 and all other couplings set to zero, We show in Fig. 2 all possible combinations ofΓ ij /Γ SM H→4L in the 2e2µ (left panel) and 4e/4µ (right panel) final state with cuts and reconstruction corresponding to 'CMS-like' phase space selections [40][41][42]. Note that in this work we have included the A ZZ/Zγ 4 couplings, which were not included in previous studies [31,34]. From Fig. 2 the normalized partial    [27,29] for 2e2µ (left) and 4e/4µ (right). We assume an on-shell Higgs boson at m h = 125 GeV and a CMS like phase space [40][41][42] for the 4L final state.
width in Eq. (21) now becomes, for the 2e2µ channel, where for completeness we have included all of the operators in Eq. (16), though below we will focus on those relevant for the nonlinear Higgs dynamics we are interested in. Note terms linear in the CP-odd couplings (A V V 3 ) do not appear because they integrate to (nearly) zero when choosing CMS-like selection cuts, which reflects the fact that "rate" measurements are not sensitive to CP violation. In this regard, one could employ the shape information in the differential spectra [31,33,34] or the construction of forward-backward asymmetries [43] to probe CP violation, but we do not explore this possibility here.
In the CMS analyses in Refs. [40,41], the signal strength in H→ 4L is measured in two categories, depending on whether the production channels involve Higgs couplings to fermions (ggH and ttH channels) or Higgs couplings to electroweak bosons (VBF and VH channels), where i = F, V represents the ggH+ttH channels and VBF+VH channels, respectively. In the following we will focus on the fermionic production channels, ggH+ttH, for two reasons: 1) ggH is by far the dominant production channel of the 125 GeV Higgs at the LHC and the uncertainty is smaller and 2) any modification in HVV couplings will enter into both the production and decay amplitudes in the VBF+VH channels. For simplicity we assume the production cross-sections in ggH+ttH channels are equal to their SM values as well as the total Higgs width. Under these assumptions the deviation in "rate measurements" in µ 4L F arises entirely from the decay amplitudes: For completeness we briefly summarize the procedures taken by the CMS collaboration in Refs. [40,41], which we refer the reader to for details. For CP-even couplings, which we focus on in this work, CMS built three multi-variate likelihood functions, each optimized for a particular anomalous Higgs coupling in Eq. (16): A ZZ 2 , A ZZ 4 and A Zγ 4 . In each likelihood function all anomalous couplings, other than the one the likelihood function is specifically optimized for, are set to zero [44]. In other words, in these analyses the anomalous HVV couplings are turned on only one at a time. (See Refs. [29,30] for a framework that allows for simultaneous measurements of several anomalous Higgs couplings at the same time.) The likelihood function allows one to constrain and fit: 1) the rate (signal strength) in H→ 4L decays and 2) the "fraction" of the observed 4L events originated from the anomalous HVV coupling. It turns out that the best fit value for the three CP-even anomalous HVV couplings analyzed in Refs. [40,41] all have central values extremely close to zero, albeit with varying degrees of uncertainties. As a result, CMS provided four different fits to the signal strength in 4L channel, under the assumption of vanishing anomalous HVV couplings.
The relevant CMS results are summarized in Table I, where we have computed the 95% C.L.
from the 67% C.L. by assuming a Gaussian distribution. As for the shape measurements, CMS provided limits on the ratio 1 1 Note we use a slightly different normalization than in the CMS analysis [40,41].  [40,41]. Each measurement is optimized with respect to a particular anomalous HVV coupling as indicated. However, the fit is performed assuming A V i = 0.
from the differential spectra in the 4L decays. This is possible because different tensor structures in Eq. (16) result in different shapes in the differential distributions. We quote the results in Table II and give a visual representation of the bounds in Fig. 3.
In their Run 2 analyses CMS did not provide a fit to A γγ/Zγ 2 [40,41], because these are better constrained from direct H→ γγ and H→Zγ two body decays. It turns out that only A Zγ 2 is universal in nonlinear Higgs dynamics, as A γγ 2 requires shift-symmetry breaking effects. In Section III B we perform a separate fit to A Zγ 2 using the direct Zγ channel. TABLE II. Summary of CMS shape measurements on ratios of anomalous HVV couplings [40,41].
In the middle column we indicate the assumptions on the Higgs total width.
It is worth emphasizing that, in order to properly constrain nonlinear Higgs dynamics, one would need to adapt the present CMS analyses. In particular, we would like to 1. fit the magnitude of the A V i coefficient in Eq. (16), instead of just the fraction.
2. allow all anomalous HVV couplings to be present at the same time.
It would be of interest to study the impact on the signal strength (µ F ) measurement by relaxing the assumption of vanishing anomalous HVV couplings. In what follows we will  Table II. perform a limited analysis within the framework of the CMS analyses, and await a more comprehensive experimental study in the future.

Constraints From Rate Measurements
Using the CMS bounds given in Table I and Table II which allows us to perform two kinds of fits to ξ: 1. In the right-panel of Fig. 4, we present bounds on ξ by setting all R i = 0 and applying the measured µ 4L F in Table I, where the constraints are obtained assuming that the anomalous HVV coupling vanishes. The values of ξ for the different analysis vary because different categorization and observables are utilized.
2. In the left-panel of Fig. 4, we use the central value of µ 4L F from Table I in the left-hand side of Eq. (26) and plug in the experimental limits on the corresponding R i from Table II in the right-hand side. We then invert Eq. (26) to derive a limit on ξ. The moral of these two fits on ξ is very different, and neither is perfect.
In the first fit, we assume all the observed 4L events came from the g µν tensor structure in Eq. (16), which is similar in spirit to the conventional approach of using the signal strength to constrain ξ without including O(p 4 ) effects. We have seen in Section II that this holds only at leading order in the derivative expansion. One could justify this approach somewhat by pointing out that the central values of the anomalous HVV couplings from the shape information are all extremely close to zero [41]. However the uncertainties remain significant, as can be seen in Table II, and it is not clear how to interpret µ 4L F when the anomalous couplings are turned on.
In the second fit, we attempted to subtract out the 4L events originating from the anomalous HVV coupling, by using the experimental limit on R V i in Table II. The remaining events then can be interpreted as arising entirely from the g µν tensor structure, whose coefficient is given by 2 √ 1 − ξ in Eq. (17). However, it is not clear what "signal strength" one should use to subtract out the anomalous 4L events, as µ 4L F in Table I is extracted assuming the anomalous coupling vanishes.
Having made these qualifications on our fitting procedures, there are interesting features in Fig. 4 which are likely to survive even after a more rigorous fitting is adopted: • The most prominent feature in the right panel of Fig. 4 is that ξ is preferred to be negative, due to the fact that µ 4L F in Table I is larger than 1 in all three different likelihood fits.
• In the left panel, the uncertainty in ξ is still rather large, allowing for ξ 0.5 when R Zγ 4 is turned on or ξ −0.5 when R ZZ 4 is allowed. It turns out that the reason behind the rather loose limits on ξ is different between these two scenarios.
In the case of R Zγ 4 , the loose constraint is due to a combination to two factors: the rather large experimental uncertainty in Table II and a large numerical coefficient, from phase space integration, in front of the |R Zγ 4 | 2 term in Eq. (26). The large experimental uncertainty implies the differential spectra from A Zγ 4 are quite similar to those from the leading order A ZZ 1 . As a consequence, the likelihood fit is unable to separate the two tensor structures.
On the other hand, while the limit on R ZZ 4 is the strongest in Table II,  It is possible to further obtain constraints on the Wilson coefficients defined in Eq. (9), by using the mapping in Eq. (17). Again we can perform two different fits using the rate and the shape information, respectively.
For rate measurements, we re-write Eq. (26) in terms of ξ, g ρ and C h i , where we have plugged in M ρ = g ρ f . With this we can use the rate measurements in 4L channel to examine bounds in the (C i , ξ) and (g ρ , ξ) two dimensional planes, by fixing the third variable. In the top row of Fig. 5, we fix g ρ = 1.5 and examine the (C i , ξ) plane.
(Recall that in QCD α s (m b ) ≈ 0.22, which corresponds to g s ≈ 1.7.) We see that turning on C h i could have a non-negligible impact on the extraction of ξ. However, the impact gets diminished as g ρ becomes larger and larger because of the 1/g 2 ρ and 1/g 4 ρ dependence in  [41]. Bottom: Same as top, but in the (g ρ , ξ) plane.
Eq. (27). This feature is demonstrated explicitly in the second row of Fig. 5 where we fix C h i = 1 and plot the constraints in the (g ρ , ξ) plane.
As for the shape measurements, since the CMS constraints are presented in terms of the percentage of 4L events originating from the anomalous HVV couplings [40,41], we find it convenient to obtain bounds on the ratios C h i /A ZZ 1 , which are independent of the 4L signal  Table II. strength. The outcome is presented in Fig. 6. We see that current bounds are still week and, depending on which coupling is varied in the fit, still allow for large positive or negative values of the C h i .

Projections with Multi-dimensional Parameter Likelihood Fits
Currently at the LHC only one parameter fits are performed in the 4L channel [40,41].
In the future this channel offers a golden opportunity to conduct multi-parameter fits, as demonstrated in the matrix element method (MEM) framework developed in Refs. [27,[29][30][31]34]. In this framework all decay observables are utilized and a combined likelihood for the 2e2µ, 4e, 4µ final states is constructed from the normalized fully differential decay width. The dominant qq → 4L background, computed analytically in Refs. [26,27,29], is also included in the likelihood. This likelihood function is a function of all anomalous HVV couplings in Eq. (16), and allowing all them to vary simultaneously we can obtain projection curves for future sensitivity to the Wilson coefficients in the nonlinear Higgs Lagrangian by the mapping in Eq. (17). Details of the statistical analysis and likelihood maximization procedure can be found in [27,[29][30][31]34].
In the left panel of Fig. 7 we show projections of the 95% C.L. contours for C h i /A ZZ simultaneously while the dashed lines are obtained allowing only a single coupling to vary, as done in current CMS analyses. In the projection we have included C h 4 , although it is not currently included in the CMS shape measurements in the 4L channel using Run 2 data.
On the right-panel we present projections for M ρ from the shape measurements using We can see that, in both plots, at large statistics single parameter and multi-parameter fits converge to similar values. However, at low statistics we see that single parameter fits could lead to an overly optimistic sensitivity to the Wilson coefficients and in particular for C h 4 . Eventually at the HL LHC [45], the sensitivity to each coupling could improve by around an order of magnitude from current limits. As can also be seen, the strongest sensitivity corresponds to C h 4 where values |C h 4 /A ZZ 1 | ∼ 4 can eventually be probed or correspondingly, scales M ρ ∼ 750 GeV. This stronger sensitivity is due to the fact the differential spectra corresponding to the A Zγ 2 can be efficiently distinguished from the dominant tree level A ZZ 1 operator [31]. This is primarily due to the photon propagator which introduces a pole at low di-lepton invariant mass. However, in the case of the A Zγ 4 (corresponding to C h 3 ), the k 2 1,2 dependence in the tensor structure (see Eq. (16)) cancels this pole thus eliminating the distinguishing feature in the differential spectra.  [46] and from the prospective HL-LHC [48].

B. Direct H → Zγ Decay
The Wilson coefficient C h 4 contributes to on-shell H → Zγ decays, which can be parametrized as [37]: where the O(1) coefficient in front of C h 4 is due to the fact that the SM H → Zγ process arises at the one-loop level. Although present measurements from the LHC on the signal strength µ Zγ only gives 6.6 at 95% CL upper limit [46,47], it already puts strong constraint with C h 4 > 0(C h 4 < 0) respectively. The HL-LHC will be able to measure this channel with 20% uncertainty [48], which will in turn give the constraint on C h 4 : In Fig. 8, we show the 95% C.L. constraint on C h 4 with M ρ = 1 TeV from the present measurement with integrated luminosity L = 36.1 fb −1 [46] and from the HL-LHC prospective [48].
While we have demonstrated that precision Higgs measurements, in particular the golden H → 4L channel, can directly constrain nonlinear Higgs dynamics, it is well-known that there are other low-energy, indirect constraints on anomalous HVV couplings from precision electroweak test (EWPT), which we consider in this section. The S, T parameters are related to ξ as follows [49]: where Λ is the UV cutoff that regularizes the logarithmic divergence. In the following, we will choose Λ = 4πf . The fact that the IR contributions to the S, T parameters have different sign put a strong constraint on the value of ξ, as the EWPT gives strong positive correlation 92% among the S and T parameter [50]. The present fit to the EWPT gives [50]: which will be used in the following analysis.
In addition to the IR contribution, there are potential UV contributions arising from the presence of O(p 4 ) operators in the effective Lagrangian. In particular, O + 5 will contribute to the S parameter at the tree-level.
There are also contributions from the one-loop threshold corrections involving the vector or fermionic resonances, which are model-dependent [51][52][53]. Instead of going into the detail of a particular UV construction, we simply choose some benchmark values for the these contributions. Note that the one loop effects can give negative and positive contributions to the S and T parameters, respectively, which are in the opposite direction from the IR contribution in Eq. (31). This observation can potentially relax the bound on ξ from EWPT.
In addition, since O + 5 contribution to the S parameter arises at the tree-level, as can been seen from the large coefficient in Eq. (33), only a small negative c + 5 is needed to relax the bound. They can be achieved by noticing that there are potential cancellations between different contributions to c + 5 in the Lagrangians of the composite spin-1 resonances [54]: where α 2 is the coefficient of the operator Q 2 defined in Ref. [54]. There are also potential positive contribution to the T parameter coming from the loop of the fermionic resonances, for example, the electroweak singlet top partner [55,56]. We are not going to discuss in detail this possibility, but only take a benchmark value of T f = 0.08 as an illustration. In the left panel of Fig. 9, we show the ∆χ 2 as function of ξ under different three different assumptions: We can see that the constraint on ξ is very strong in the absence of other contributions: ξ ∈ [−0.11, 0.002] at 95% C.L., which again prefers a negative value of ξ. A small negative value of c + 5 relax the bound on the negative ξ to −0.19 and additional positive contribution from T f can relax the bound on positive value of ξ ∈ [−0.015, 0.16]. In the right panel we present the 95% C.L. allowed region in the ξ − c + 5 plane with g ρ = 1.5. We can see that the bound on ξ is very strong with only IR contribution and can be significantly relaxed in the presence of small negative value of c + 5 and positive contribution T f , which is consistent with the results discussed above.
We emphasize that the precision electroweak constraints are orthogonal to direct measurements of Higgs couplings, as they involve different sets of assumptions. Therefore it is important to pursue both of them independently.

IV. CONCLUSIONS
Nonlinear dynamics of the 125 GeV Higgs boson is the most salient feature of a composite Higgs boson. The nonlinear interaction realize the PNGB nature of the Higgs, in the same way pions in low-energy QCD manifest their PNGB nature through the nonlinear interactions. In fact, the nonlinear dynamics can serve as the defining property of a composite Higgs boson, independent of how the fermionic sector is implemented, even when the so-called top partner is neutral under all SM charges. In particular, recent theoretical advances showed the nonlinear interaction is an IR property of the composite Higgs boson that is insensitive to the symmetry-breaking pattern invoked in the UV. At the two-derivative level, the nonlinear Lagrangian is determined by a single parameter f , the Goldstone decay constant.
In the low-energy regime, f needs to be taken as an input from the experimental data. In order to disentangle effects of O(p 2 ) operators from those of O(p 4 ) operators, it is crucial to include shape information, as different tensor structures lead to different shapes in the fully differential spectra. In this regard, the H→ 4L "golden channel" decay is the ideal probe of these nonlinear Higgs dynamics. In this work we have performed, for the first time, experimental fits to the Wilson coefficients of O(p 4 ) operators. In addition, we demonstrated the limitation of using only the signal strength to constrain the nonlinear parameter ξ at leading order in derivative expansions. We showed that it is important to include both the rate information in the signal strength measurements and the shape information in the fully differential spectra, in order to constrain the nonlinear parameter ξ and the Wilson coefficients C h i . The fitting procedures we adopted are less than ideal and limited by the methods employed in current LHC analyses.
We found that in rate measurements ξ is preferred to be negative, pointing to a non-compact coset structure in the UV [17,18]. Moreover, using the shape measurements we identified scenarios where ξ could be as large as +0.5 or -0.5, corresponding to turning on A Zγ 4 and A ZZ 4 , respectively. Such a large ξ could be further constrained in a specific UV model, either by precision electroweak measurements or direct searches in the fermionic sector. It would be an interesting model-building challenge to devise a concrete UV model realizing such a large ξ while avoiding other experimental constraints.
Our work points to a new experimental frontier in using precision Higgs measurements to probe nonlinear dynamics of a composite Higgs boson, beyond the conventional signal strength measurements. Given that the Higgs boson remains a top priority in current and future experimental programs in high-energy colliders, it is desirable to introduce new tools and techniques to further the experimental analysis in measurements of Higgs properties.
We leave this direction for future works.