Extracting Effective Higgs Couplings in the Golden Channel

Kinematic distributions in Higgs decays to four charged leptons, the so called `golden channel, are a powerful probe of the tensor structure of its couplings to neutral electroweak gauge bosons. In this study we construct the first part of a comprehensive analysis framework designed to maximize the information contained in this channel in order to perform direct extraction of the various possible Higgs couplings. To that end we first complete an earlier analytic calculation of the leading order fully differential cross sections for the golden channel signal and background to include the 4e and $4\mu$ final states with interference between identical final states. We also examine the relative fractions of the different possible combinations of scalar-tensor couplings by integrating the fully differential cross section over all kinematic variables as well as show various doubly differential spectra for both the signal and background. From these analytic expressions we then construct a `generator level' analysis framework based on the maximum likelihood method. We demonstrate the ability of our framework to perform multi-parameter extractions of all the possible effective couplings of a spin-0 scalar to pairs of neutral electroweak gauge bosons including any correlations. This framework provides a powerful method for study of these couplings and can be readily adapted to include the relevant detector and systematic effects which we demonstrate in an accompanying study to follow.


I. INTRODUCTION
With the recent discovery of the Higgs boson at the LHC [1,2] the focus now shifts to the determination of its detailed properties and in particular whether or not it possesses any anomalous couplings not predicted by the Standard Model. Ideally, the constraining or measuring of these couplings should be done through direct parameter extraction with minimal theoretical assumptions. The vast literature  on Higgs decays to four charged leptons (electrons and muons) through neutral electroweak gauge bosons, the so called 'golden channel', suggests that it can be a powerful channel in accomplishing this goal.
In addition, the high precision with which this channel is measured allows for one of the best opportunities to use analytic methods to analyze data. As has already been suggested for the golden channel [13,17,27] and to be further emphasized here, analytic methods are optimal for performing direct multi-parameter extraction within a minimal amount of computing time. Furthermore, as we show in an accompanying study [30,31], within an analytic framework one can also include the relevant detector effects in order to obtain a 'detector level' likelihood in terms of the full set of observables available in the four lepton final state. Of course other frameworks have also been recently constructed to study the golden channel * yi.chen@cern.ch † roberto.vega@th.u-psud.fr robertovegamorales2010@u.northwestern.edu (see for example recent Madgraph [21,32] or JHU generator [27] based implementations which also include the possibility to study other Higgs decay and production channels). In this study we construct the first part of a comprehensive analysis framework, based on a largely analytic implementation, designed to maximize the information contained in the golden channel in order to perform direct extraction of the various effective Higgs couplings.
We begin by extending our previous leading order analytic calculations [33], for both the signal and background in the 2e2µ final state, to now also include the 4e final state. We include the interference between identical final states as well as interference between all intermediate states. Explicitly we calculate for the signal process ϕ → ZZ+Zγ+γγ → 4e/4µ where ϕ is a spin-0 scalar and we have allowed for all possible tensor structures. This covers all possible couplings of a spin-0 scalar to ZZ, Zγ, or γγ pairs. For the dominant irreducible background we compute qq → 4e/4µ including both the t and s-channel process mediated by Z and γ vector bosons. All vector bosons are allowed to be on or off-shell and we do not distinguish between them in what follows.
After presenting the calculation of the analytic fully differential cross sections, we then examine various aspects of the golden channel in more detail. First, we isolate the individual contributions to the golden channel signal by obtaining the 'partial fractions' for each possible combination of tensor structures which can contribution to the ϕ → ZZ + Zγ + γγ → 4 (where 4 = 2e2µ, 4e, 4µ) process. This is done by integrating the differential cross section over the set of kinematic variables for a given arXiv:1310.2893v2 [hep-ph] 26 Mar 2014 phase space. These partial fractions give an indication of the relative contributions of each component to the golden channel and a rough picture of the potential sensitivity to the various tensor structures. As part of this integration we also show a number of doubly differential spectra for signal and background in the Appendix.
We then construct a maximum likelihood analysis using the analytic expressions of the fully differential cross sections to build the probability density functions (pdfs). This framework builds upon and extends recent studies which first introduce using analytic expressions to perform parameter extraction in the golden channel [13,17,27]. Using these analytic pdfs, we study the ability of the golden channel to directly extract the couplings between a spin-0 scalar and ZZ, Zγ, and γγ pairs. We validate our analysis framework by performing a number of simplified 'generator level' studies. To do this we choose an example parameter point in which all possible operators are simultaneously 'turned on' in order to demonstrate the validity of our maximization procedure as well as our ability to simultaneously extract the various couplings as well as their correlations.
Of course a proper treatment of the golden channel requires careful study of detector resolution and acceptance effects. This also includes an adequate treatment of the production variables for both signal and background as well as taking into account higher order contributions. We leave these issues to an accompanying paper [30] where we construct a 'detector level' analysis which includes a treatment of all these issues as well as systematic uncertainties while retaining the flexibility and speed in parameter extraction which we present at 'generator level' in this study.
The organization of this paper is as follows: in Sec. II we briefly review the kinematics of the four lepton final state. In Sec. III we describe the calculation of the signal fully differential cross section while in Sec. IV we describe the calculation of the background fully differential cross section. In Sec. V we examine the relative fractions of all the possible operators which might contribute to ϕ → ZZ + Zγ + γγ → 4 . We then present our analysis framework and perform an example parameter extraction to motivate the possibility of extracting the various couplings directly. We also comment on ongoing and future studies before concluding in Sec. VI. In the Appendix in Sec. VII we also show various 2D projections for both the signal and background in the 4e channel as well as the relative fractions for a second set of phase space cuts.

II. KINEMATIC VARIABLES
In this section we briefly discuss the set of observables used to parameterize the ϕ → ZZ +Zγ +γγ → 4 (where 4 = 2e2µ, 4e, 4µ) and qq → 4 fully differential cross sections. The kinematics of four lepton events are described in detail in [13] and are illustrated in Fig. 1. The invariant masses in the system are defined as, Here Z1 and Z2 can be either Z or γ.
• √ s ≡ m ϕ -The invariant mass of the four lepton system or equivalently the Higgs mass for the signal case.
• M 1 -The invariant mass of the lepton pair system which reconstructs closest to the Z mass.
• M 2 -The invariant mass of the other lepton pair system.
These variables are all independent subject to the constraint (M 1 + M 2 ) ≤ √ s. Note also that the 4e/4µ final state can be reconstructed in two different ways due to the identical final state interference. This is a quantum mechanical effect that occurs at the amplitude level and thus both reconstructions are valid. The definitions M 1 and M 2 remained unchanged however. The angular variables are defined as, • Θ -The 'production angle' between the momentum vectors of the lepton pair which reconstructs to M 1 and the total 4 system momentum.
• θ 1,2 -Polar angle of the momentum vectors of e − , µ − in the lepton pair rest frame.
• Φ 1 -The angle between the plane formed by the M 1 lepton pair and the 'production plane' formed out of the momenta of the incoming partons and the momenta of the two lepton pair systems.
• Φ -The angle between the decay planes of the final state lepton pairs in the rest frame of the 4 system.
We have ignored the off-set angle φ, defining a global rotation of the event which is 'flat' and thus not shown. We can group the angular variables as follows Ω = (Θ, cos θ 1 , cos θ 2 , Φ 1 , Φ).
There are also in principal the 'production' variables associated with the initial partonic state four momentum. This four momentum defines the invariant mass of the CM frame ( √ s), as well as the rapidity (Y ) defined as the motion along the longitudinal direction, and the momentum in the transverse direction ( p T ). In principal the inclusion of Y and p T as observables would increase the discriminating power of the golden channel, but as we are interested primarily in parameter extraction and these variables introduce additional systematic uncertainties we do not consider them in our set of observables. When including detector effects, however, these production variables must be properly accounted for as we will do in [30].

III. SIGNAL
In this section we present the calculation of the signal fully differential cross section for the process ϕ → ZZ + Zγ + γγ → 4e/4µ. We take ϕ to be a general spin-0 scalar and consider all possible couplings to any combination of Z and γ pairs allowing for mixtures of both CP even and odd interactions. We follow closely, with a slight variation in strategy and notation, the method used in [33] for the calculation of the 2e2µ final state and refer the reader there for many of the details. Here the only additional calculation needed is that for the identical final state interference in the 4e/4µ channels. Various validations of the calculation can be found in the Appendix as well as [33], [30], and [31].

A. Parametrization of Scalar-Tensor Couplings
The general couplings of a scalar ϕ to ZZ, Zγ or γγ pairs can be parametrized by the following vertex 1 , where ij = ZZ, Zγ, or γγ and k and k represent the four momentum of the intermediate vector bosons with v the Higgs vacuum expectation value (vev) which we have chosen as our overall normalization. The A 1ij,2ij,3ij are dimensionless arbitrary complex (momentum dependent) form factors. For the purposes of this study however, we will approximate the couplings as constant as is done in other similar analysis [13,17,21,27,32] though our framework can easily be made to include the full momentum dependence of the form factors. For the case of a scalar coupling to Zγ or γγ, electromagnetic gauge invariance requires A 1 = 0, while for ZZ it can be generated at tree level as in the SM or by higher dimensional operators. We can also write Eq.(1) as, where the coefficients A nij and Lorentz structure V µν n are those found in Eq. (1). Although it is more general, the parametrization in Eq.(1) can for example be mapped onto the Lagrangian 2 given by, where we have allowed only up to dimension five operators and Z µ is the Z field while V µν = ∂ µ V ν − ∂ ν V µ is the usual bosonic field strengths. The dual field strengths are defined as V µν = 1 2 µνρσ V ρσ . Thus for this Lagrangian we would have, A 1ZZ ≡ g h , A 2ZZ ≡ g Z , A 3ZZ ≡g Z and similarly for Zγ and γγ. This makes Eq.(1) a convenient parametrization for fitting to Lagrangian parameters that might be generated in various models at dimension five or less. If ϕ is purely the Standard Model Higgs, then A 1ZZ = 2, while all other coefficients are taken as approximately zero 3 . Note also that in this parameterization we have not made any theoretical assumptions about the nature of ϕ such as imposing that the couplings are related by SU (2) L ⊗U (1) Y gauge invariance for example.
We note that it is important to include all possible Higgs couplings including the Zγ and γγ contributions in the signal differential cross section. This is because since the Higgs appears to be mostly 'standard model like' [36] this means we are primarily searching for small anomalous deviations from the standard model prediction. Thus when attempting to extract specific couplings we must be sure that one small effect is not being mistaken for another (this is also why it is important to include the interference effects between the identical final state leptons). This is especially true because as we will see many of the couplings are correlated. Including all possible couplings and doing a simultaneous fit ensures we minimize the possibility of introducing a bias when attempting to extract these couplings.
We allow for all vertex structures in Eq. (1) to contribute simultaneously including all possible interference effects. Of course Eq.(1) can be mapped onto Lagrangians with dimension greater than five with appropriate translation of the parameters, but we work explicitly with the vertex in Eq.(1) and Eq.(2) when calculating the fully differential cross section for ϕ → ZZ + Zγ + γγ → 4e/4µ. Below we summarize the details of the calculation which is performed using the Tracer [37] package in Mathematica [38] to perform the necessary Dirac algebra.

B. Calculation
To compute the process ϕ → ZZ + Zγ + γγ → 4e/4µ we include the diagrams shown in Fig. 2 where i, j = Z, γ and parameterize the scalar coupling to gauge bosons as in Eq. (2). For any ij intermediate state, the amplitude M OF ij exists for both the opposite flavor (OF) 2e2µ final state as well as for the same flavor (SF) 4e/4µ final state. The amplitude M SF ij which is obtained by exchanging the four momentum of the particles (or anti-particles), is only present for the 4e/4µ final state. The total amplitude for any particular intermediate state is the sum of the two diagrams and can be written as, Assuming the final state leptons to be massless, we can write the OF amplitude as, where i, j label Z or γ while 1 and 2 label the final state leptons and can in principal be e or µ. The vector boson four momenta are given by k xy = (p x + p y ) where p x are the four momentum of the final state leptons. Note that we have also set k = k 11 and k = k 22 in the vertex function Γ µν ij . The SF amplitude can be obtained from the OF amplitude by swapping u 1 ↔ u 2 as well as p 1 ↔ p 2 and can be written as, where note an overall minus sign is included to account for the swapping of identical fermions and now k = k 12 and k = k 21 in the vertex Γ µν ij . Upon squaring Eq.(4) this gives for the amplitude squared, The M OF ij M * OFīj term is equivalent to the 2e2µ matrix element squared which was calculated in [33]. We repeat this part of the calculation here for clarity and consistency of notation. After summing over final state lepton polarizations we can obtain a general amplitude squared which encompasses any combination of interme-diate states and is given by, where Γ µν ijn are given in Eq.(1) and Eq.(2) and we have defined the objects, for the Dirac strings while for the propagators we have, OF SF Feynman diagrams [39] contributing to ϕ → ViVj → 4 where 4 = 2e, 4e/4µ and Vi,j = Z, γ. The arrows are to indicate the direction of momentum flow and 1, 2 label the lepton momenta. On the left we have the opposite (OF) flavor diagram present in both the 2e2µ and 4e/4µ channels. On the right we have the same flavor (SF) flavor diagram present only in the 4e/4µ channel. Note also that the diagram on the right hand side implicitly comes with an overall minus sign to account for the switching of identical fermions (1 ↔ 2).
The g i R,L represent the lepton couplings to Z and γ, but are in fact at this point general left and right handed couplings of a 'Z-like' spin-1 vector boson to a pair of fermions. The bars on Lorentz, i, j, and n indices are to indicate that the corresponding index belongs to the conjugated amplitude and are distinct indices from the un-bared ones. We treat all couplings at every vertex encountered when tracing over the Dirac strings as distinct as well as all Breit-Wigner factors so for any amplitude squared term there can in principal be four different vector bosons as intermediate states. In the case of the photon we have of course g γ R = g γ L = −e em and m γ = Γ γ = 0.
After expanding Eq.(8) we obtain, where a, b = (±, ±) with a and b corresponding to the fermion pairs labeled 1 and 2 respectively in the OF diagram of Fig. 2 and have defined, The T σσ 1± are the Dirac traces found in Eq.(9) and ± indicates whether the trace ends with a γ 5 (−) or not (+). From the objects in Eq. (12) we can go on to obtain the full amplitude squared for the 2e2µ channel as done in [33].
For the 4e/4µ final state we also have the second squared term M SF ij M * SFīj , but this is obtained easily from M OF ij M * OFīj by swapping p 1 ↔ p 2 . Thus the only new term left to calculate in the 4e/4µ case is the interference term M OF ij M * SFīj . Note also that the amplitudes in the 4e/4µ case come with a symmetry factor of 1/2 for the identical final states, which we explicitly add at a later step. After squaring the amplitude we find for the interference term, T ijīj (p 2 , γ, p2,γ, p 1 , σ, p1,σ) where Γ µν ijn are given in Eq.(1) and Eq.(2) and we have defined, Expanding out the terms in Eq.(13) we can write the interference term as, where the coefficients and Lorentz structure are now, The T γγσσ

12±
are the Dirac traces found in Eq. (14) and again ± indicates whether the trace ends with a γ 5 (−) or not (+). Note that again the vector boson momentum in V µν n of Eq.(16) is given by k = k 11 and k = k 22 , but now in V * μν n we have k = k 12 and k = k 21 . We can now take advantage of the fact that L ±± nn and L ± nn are independent of the intermediate state vector bosons to perform the sum over i, j = Z, γ and obtain general coefficients for the Lorentz structure which include all contributions from Z and γ gauge bosons, The full amplitude squared for ϕ → ZZ + Zγ + γγ → 4e/4µ can then be built out of the objects in Eqs. (12), (16), and (17) as follows 4 , where we have included the 1/4 symmetry factor for the identical final state fermions. One can also easily obtain the amplitude squared for any combination of vertex structures in Eq.(1) by not taking the sum over n andn and choosing the desired n,n combination. We will take advantage of this property when performing integration and when we examine the interference effects between different operators below. The final fully differential cross section (which is treated at fixed √ s) can then be obtained via, where d Ω = dc Θ dc θ1 dc θ2 dΦdΦ 1 (c θ = cos θ) and Π 4 is the final state massless lepton four body phase space derived following [41] and given by, 4 Analytic expressions may be obtained by emailing the authors or at a website which is currently under construction [40].
Unlike the 2e2µ final state, the coefficients C a nn in the interference term of Eq.(18) depend on the polar angles cos θ 1,2 and in particular through the denominators of the vector boson propagators (see Eq. (6)). This makes analytic integration difficult. Thus analytic expressions for the doubly differential mass spectra are not obtained in the 4e/4µ channel as they were for 2e2µ [33]. In Fig. 16 of the Appendix we show plots for the differential mass spectra after performing the angular integration numerically as well as various other doubly differential distributions for the SM signal hypotheses. Again, details of the validation procedure can be found in [30,31,33].

IV. BACKGROUND
The dominant irreducible background to the golden channel comes from qq annihilation into gauge bosons. At energies ∼ 125 GeV the dominant contribution comes from t-channel Zγ production [33]. However, as was seen in [33] contributions from s-channel process diagrams can effect the angular distributions, such as the distribution of the azimuthal angle between the lepton decay planes Φ defined in Sec.II. Furthermore, we include the ZZ and γγ contributions since in principal these are always present and may have observable interference effects due to the fact that they add at the amplitude level when decaying to charged leptons and can mimic some of the effects of the signal tensor structures. Of course higher order effects, including the gg initiated process [42][43][44] will contribute as well, but these are expected to be subdominant and mainly only effect the 'input' invariant mass (and overall normalization) for the fully differential cross sections. Since we are not including production variables in our set of observables and are not concerned with the overall normalization, neglecting these contributions has a minimal effect on our analysis framework, but as mentioned previously should properly be taken into account when including detector effects.
In this section we extend a previous calculation of the 2e2µ channel to include the 4e/4µ final state. The calculation follows in the same manner as for 2e2µ (with some slight changes in notation) except that now one must include the contribution from interference between the final state identical particles, which in some kinematic regimes can have non-negligible effects [21,42]. In this section we describe the calculation of this interference, while the parts of the calculation which are identical to the 2e2µ case can be found in [33]. Again we use the Tracer [37] package in Mathematica [38] to perform the necessary algebra.

A. Calculation
The background calculation is significantly more involved than the signal calculation due to a much larger number of Feynman diagrams (48 in total as opposed to 8 for signal) in addition to a more complicated Lorentz structure. As in the signal case the amplitude can be written as sum of opposite flavor (OF) amplitude and a same flavor (SF) amplitude. Thus the amplitude squared can again be written as, The first term M OF ij M * OFīj is exactly equivalent to the 2e2µ amplitude squared calculated in [33] to which we refer the reader for details. The second term M SF ij M * SFīj can be easily computed from the first by the simple exchange p 1 ↔ p 2 as was done in the signal case. Thus the only new term left to calculate is the identical final state interference term M OF ij M * SFīj . Following the strategy in [33] we organize the diagrams after 'twisting' them into the form shown in Fig. 3 where we allow the intermediate vector bosons to take on any combination of Z or γ, but once chosen are treated as fixed. We use the conventions indicated in [33] and in particular refer to all of the diagrams in Fig 3 as 'tchannel' type diagrams while the 'u-channel' diagrams are obtained by switching the vertex at which the vector bosons are attached. This is not to be confused with the typical vocabulary for this process which refers to diagrams (a) and (b) as t-channel and diagrams (c) − (f ) as s-channel. We find this re-naming convenient for organizing and reducing the many terms which need to be computed for the differential cross section (see [33] for a more detailed explanation). The diagrams on the left are labeled by OF, while those on the right are labeled by SF. Note also that the diagrams on the right hand side implicitly come with an overall minus sign to account for the switching of identical fermions (p 1 ↔ p 2 ).
The Lorentz structure for all of these amplitudes is clearly the same. One needs only to keep proper track of how the various momentum are routed through each diagram. We can see this by considering the amplitude explicitly. Using the massless initial quark and final state lepton approximation we can write any of the OF amplitudes on the left hand side in Fig. 3 as, where we label the amplitude by the 'long' Dirac string, in this case X. The labels X/Y /Z = 1, 2, q where 1, 2 are for final state lepton pairs while q is for the initial state quarks. The i, j = Z, γ label the vector bosons and n = t, u labels the t and u-channel diagrams in our new vocabulary. The internal vector boson momenta are again defined as k xy = (p x + p y ), while the internal fermion momentum are given by, For any of the SF amplitudes a similar formula as in Eq.(22) applies except we take p 1 ↔ p 2 and multiply by an overall minus sign in the corresponding OF amplitude with the quark string in the same position (this simply corresponds to diagrams in the same row of Fig. 3). Thus we have for the SF amplitude, while the internal fermion momentum are now given by, . To obtain any of the physical amplitudes one simply assigns the appropriate labels to Eq. (22) or Eq.(24) as well as the appropriate momenta. Thus for example for diagram (c) we have X → 1, Y → q, Z → 2, and n → t. To switch from t-channel type to u-channel diagrams one simply takes t → u and γ σ ↔ γ γ while to obtain the corresponding SF diagram simply take OF → SF and 1 ↔ 2 and multiply by an overall sign. Note that for the Z propagators we drop the momentum dependent terms since they do not contribute in the massless lepton approximation.
As in the case of the signal, the next step is to find a generalized amplitude squared for any two of the six diagrams. Since we are only concerned with obtaining the interference term M OF ij M * SFīj we need only consider the terms coming from multiplying the amplitudes on the left hand side (OF) with those on the right hand side (SF) [39] contributing to qq → ViVj → 4e/4µ and qq → Vi,j → 4e/4µ where Vi,j = Z, γ and 1, 2 label the lepton momenta. Note that although we define all diagrams as 't-channel' type, diagrams (c) − (f ) are in fact s-channel type in the usual convention so the fermions labeled by 1 and 2 are not to be confused as being in the initial state. This is taken into account in how the various momenta are assigned as indicated by the arrows. The diagrams on the left hand side are labeled by OF, while those on the right are labeled by SF. Note also that the diagrams on the right hand side implicitly come with an overall minus sign to account for the switching of identical fermions (1 ↔ 2).
in Fig. 3. These organize themselves into three distinct types of Lorentz structure. The first type is found when multiplying the two diagrams in the top row of Fig. 3 (corresponding to t-channel di-boson production in the conventional language). These give, where the D xyi are defined in Eq.(10) and the T X ijīj are defined in Eq. (14). Again the bars on Lorentz, i, j, and n indices are to indicate that the corresponding index belongs to the conjugated amplitude and are distinct indices from the un-bared ones. Expanding out the terms in Eq.(26) we can organize in a manner similar to Eqs. (11) and (15) writing the amplitude squared as, where again a, b = (±, ±) with a and b corresponding to the quark and lepton strings and we have defined the Lorentz structure coefficients, and Lorentz structure, where the T objects are the traces found in Eq. (26). The next type of Lorentz structure is found for any OF/SF pair of diagrams in (c)−(f ) (interference between s-channel diagrams in the usual language). For those in the same row we can write, where here Y /Z = 1, 2 while the T X jj are defined in Eq. (9) and we have also defined, Expanding out Eq.(30) we can write the amplitude squared as, where again a, b = (±, ±) and we have defined the Lorentz structure coefficients, (33) and Lorentz structure, where the T objects are the traces found in Eq. (30). For products of diagrams in different rows in (c) − (f ) we obtain the following, × T ijiījī (p Y , ν, P OF Y n , µ, pȲ ,γ, p Z , σ, pZ,μ, P SF Zn ,ν). Again expanding out Eq.(35) we can write the amplitude squared as, where again a, b = (±, ±) and we have defined the Lorentz structure coefficients, (37) and Lorentz structure, where the T objects are the traces found in Eq. (35). The final type of Lorentz structure occurs when a diagram from the first row (t-channel quark exchange diagram) interferes with one of the diagrams in (c) − (f ) (s-channel process in the usual language). For these we can write, where Y = 1, 2 and we have defined, as well as, After expanding out Eq.(39) we can write the amplitude squared as, where again a, b = (±, ±) and we have defined the Lorentz structure coefficients, and Lorentz structure, where the T objects are the traces found in Eq. (39). As in the signal case we take advantage of the fact that the Lorentz structures in Eqs. (29), (34), (38), and (44) are independent of the intermediate vector bosons to perform the sum over i, j in the Lorentz coefficients defined in Eqs. (28), (33), (37), and (43) to obtain, In this way we easily take into account all possible combinations of intermediate vector bosons.
We now have all of the pieces 5 necessary to build the total interference term M OF ij M * SFīj in Eq. (21) including all contributions from the intermediate vector 5 Expressions for the various coefficients and Lorentz structure can be obtained by emailing the authors or at [40].
bosons. Explicitly we have, where the sum over intermediate vector bosons has already been implicitly performed in Eq.(45) while the sum over n,n which includes the t and u channel contributions is shown explicitly (note that this also factors from the vector boson sum). The C ab XY coefficients are in general complex due to the factor of i multiplying the decay width in the massive vector boson propagators. The Lorentz structure is either purely real or purely imaginary depending on whether the term contains an even or odd number of traces ending in γ 5 . These traces give an overall factor of i (and an epsilon tensor). Thus if L ab XY nn contains an even number of these traces, then it is purely real and if it contains an odd number it is purely imaginary. Organizing in this manner allows for easier integration when obtaining the various projections (as well as when performing convolution to include detector effects [30,31]).
Plugging Eq.(46) into Eq.(21) and using the results from [33] as well as the fact that, we can obtain the complete amplitude squared for the qq → 4e/4µ background process, where we have included a symmetry factor of 1/4 and implicitly included a color factor of 1/3 as well as a 1/4 for averaging over initial state quark spins. Again the fully differential cross section is found by combining with the lepton four body phase space in Eq. (20) to give, This expression can now be combined with the result for the signal differential cross section to perform detailed analysis of the golden channel. As in the case for signal, one also finds in the interference terms a dependence on cos θ1,2 in the propagator denominators, making it difficult to perform analytic integration over the angular variables to obtain the doubly differential mass spectrum as was done in the 2e2µ case [33]. We thus perform this integration numerically and show in Fig. 17 of the Appendix the doubly differential mass spectra as well as various other doubly differential distributions. Again details of the validation procedure can be found in [30,31,33].

V. SCRUTINIZING THE GOLDEN CHANNEL
In this section we explore the potential of the golden channel to elucidate the nature of the couplings of a spin-0 scalar to neutral electroweak gauge bosons. We begin by examining the relative contributions of all the possible combinations of tensor structures in Eq.(2) to the total ϕ → 4 decay width. We then perform a 'toy' generator level analysis to demonstrate our parameter extraction procedure via maximization of the likelihood. We present various parameter fits to show the flexibility of our framework and its ability to extract the effective couplings including their correlations. We only focus on 'toy' parameter extractions in this study, since a proper study of the Higgs couplings requires careful inclusion of the relevant detector effects as well as an adequate treatment of production variables. We leave a more detailed investigation of the Higgs couplings in the golden channel including detector effects to an accompanying study [30].

A. Relative 'Partial Fractions'
The total decay width for ϕ → ZZ + Zγ + γγ → 4 can be decomposed into the various 'partial widths' formed out of pairs of tensor structures in Eq.(2) (or operators if interpreted in terms of Eq. (3)). Since each term will be quadratic in the couplings, we can label each partial width by the appropriate combination of couplings AnijA * nīj . They are obtained by integrating the fully differential decay width in Eq.(19) over the kinematic variables defined in Sec. II. We then normalize these partial widths to the standard model value to form the various 'partial fractions'.
We show in Fig. 4 a table of these partial fractions for every possible combination of AnijA * nīj which can contribute to the 2e2µ decay width. For these tables we take as our phase space 4 GeV < M1,2 and √ s = 125 GeV as well as p T > 2 GeV and |η | < 2.4 for the transverse momentum and rapidity respectively of the final state leptons. The couplings Anij have been separated into their real and imaginary components as Anij = AnijR + iAnijI and we have set all AnijR,I = 1. All of the |AnijR,I | 2 terms sit along the diagonal with the various interference terms making up the off-diagonal terms. Note that many of the interference terms are negative indicating destructive interference between the corresponding tensor structures (or operators).
In Fig. 5 we show the same plot for the 4e final state. One can see the change in the partial fractions and in particular the ZZ/γγ interference terms are significantly larger than in the 2e2µ channel. The blank entries indicate terms which are identically zero after integration. We can see that these entries are those for which CP violation in the form of interference between A1,2 and A3 tensor structures would occur. This is indicative of the fact that after one integrates over the kinematic all information on CP violation is lost. Of course for the fully differential decay width many of these terms are non-zero in principal allowing for sensitivity to CP violation in the golden channel. To get a rough idea of the size of these CP violating terms, in Figs. 14 and 15 in the Appendix we show the integral of the absolute value of the differential decay width.
Since all couplings are set to one, these tables essentially show how much each combination of tensor structures contributes to the ϕ → 4 phase space relative to the standard model contribution for which we have set A1ZZ = 2 and all other couplings to zero. From these values of the relative partial widths, one can gain some insight into which combination of operators the golden channel might be most sensitive to. Furthermore, for a specific model one can take the prediction for the values of the various couplings and simply multiply by the numbers given in Fig. 4-5 to get a feel for whether those couplings might be probed in the golden channel. For most realistic models, all couplings apart from A1ZZ are generated by higher dimensional operators and are expected to be small. Of course, these rates do not contain information about the shapes in the various distributions so in principal the sensitivity is greater than might be inferred from these values. In Sec. VII A of the Appendix we also show the same partial fractions for a 'CMS-like' phase space as well as show the same tables for the standard model prediction. Of course for a scalar resonance with a mass much larger or smaller than 125 GeV these numbers can change significantly.

B. Simplified Analysis
In order to demonstrate the flexibility and potential of our framework, we perform a simplified generator level analysis neglecting any detector effects and at a fixed center of mass energy of √ s = mϕ = 125 GeV . To do this we construct a maximum likelihood analysis using the fully differential cross sections in Eqs. (19) and (49) to build the signal plus background pdf from which the total likelihood will be constructed. Thus we have, where O = (s, M1, M2, Ω) is our final set of observables and f is the background fraction, which we must also extract. The signal and background pdfs are given by, where they have been normalized over O (at fixed √ s). With the pdfs in hand we can now write the likelihood of obtaining a particular data set containing N events as, After constructing L(f, λ) we then maximize with respect to f and λ to extract the values which maximize the likelihood λ andf for a given data set. To asses the error we then repeat this for a large number of pseudo experiments to obtain distributions forλ andf with a corresponding spread. Below we show the results for an example parameter point. More details on this procedure can be found in [30] and [31].

C. Fit Definition
To examine the Higgs couplings to neutral gauge bosons, we take as our hypothesis the vertex in Eq.(1). We can use an overall phase rotation to make one of the parameters real. Furthermore, we can avoid the need for the absolute normalization if we instead fit to ratios of couplings. Which parameter to make real and which ratios to construct explicitly is a matter of choice the most convenient of which depends on the fit being performed. Thus, in terms of the vertex as defined in Eqs.(2), we are explicitly fitting to, where R ij n are complex ratios defined as R ij n = Anij/|A| where |A| is some normalization to be chosen for each fit. Since one of the R ij n can always be made real there are in principal twelve undetermined parameters to fit for when neglecting the overall normalization (note R Zγ 1 = R γγ 1 = 0). Fitting to ratios also makes any dependence on the production variables, pT and Y minimal since they mainly only affect selection efficiencies when detector effects are eventually included [30].

D. Example Parameter Extraction
As a demonstration of our ability to perform parameter extraction, we analyze the following example parameter point: Note that even though A2ZZ is zero we still fit for it and therefore it is floated when performing the maximization. Thus we allow for all operators in Eq.(3) to be 'turned on' simultaneously, but we assume all coefficients to be real. Our framework can easily also allow for non-zero phases, but we do not consider them here for simplicity. The pseudo-data set to which we fit is obtained by generating large samples from the analytic expressions using a simply constructed event generator 6 . We generate both signal and background events at fixed energy √ s = 125 GeV and M1,2 > 4 GeV . Since we seek only to demonstrate the validity of our parameter extraction framework, we focus on the 2e2µ final state for simplicity. It would be interesting, however, to perform a dedicated study and examine how the sensitivity of the 2e2µ final state compares to the 4e/4µ final state for different choices of phase space, but we leave this for future work. The parameter extraction is performed by maximizing the likelihood function as described above.
We first perform a simultaneous extraction of all parameters including the correlations assuming a pure signal sample. We show in Fig. 6 one dimensional results for a large set of pseudo experiments containing 1000 events each. We have explicitly fit to the ratios of couplings R ij n = Anij/|A| where here we take |A| = |A1ZZ | (thus fixing R ZZ 1 = 1). The distribution for the extracted parameters obtained for the set of pseudo experiments is shown in blue with the true value indicated by the red vertical line. One can see that the true value sits near the center of the distribution, an indication that the maximization procedure is working properly and that the global maximum of the likelihood function is in fact being obtained in each pseudo experiment. The efficiency of convergence in our maximization is 99% and takes on the order of a few minutes to complete [30,31].
Of course there are also correlations between the parameters. To see this we can examine the different parameters  in pairs as shown in Fig. 7 again for 1000 events for each pseudo experiment and assuming a pure signal sample. The true value is indicated by the intersection of the two solid black lines which as can be seen falls near the center of the distribution. The colors indicate the density of pseudo experiments returning a particular value for the extracted parameters as indicated on the x and y axis and we have fit to R ij n = Anij/|A| where |A| = |A1ZZ |. We can see that in the majority of pseudo experiments the fit returns values close to the true ones. One can also see in these plots some of the potential correlations between the various parameters though of course the full set of correlations between the six parameters which are contained in the fit can not be displayed easily. This also demonstrates the importance of including all possible couplings in the differential cross section.
We can also examine how the results change when the qq → 4 background is included. As discussed above, when including both signal and background we must also extract the background fraction f . In Fig. 8 we show our results including background in the likelihood. In the top left plot we show the distribution of the extracted background fraction for 1000 signal plus 250 background events for a large set of pseudo experiments. In the additional plots we compare the results assuming a pure signal sample shown in black to those which include both signal and background which are shown in red. We fit to R ij n = Anij/|A|, but now take the overall normalization to be |A| = n,ij |Anij| 2 . We can see that the couplings which are affected the most by the inclusion of background are the Zγ couplings. This can be understood by the fact that near the signal region of 125 GeV , the background is primarily composed of the Zγ intermediate state [33]. In general, however, one can see that the effect of including background is small, an indication that there is strong discrimination between signal and background as implied by the differences in the various doubly differential spectra shown in [33] and in Figs. 16-17. It is also important to study how the spread of the distribution changes as a function of the numbers of events. In Fig. 9 we have plotted the results for various number of signal events ranging from 30 to 3000 per pseudo experiment with a 20% background fraction. The color indicates the density of pseudo experiments which return a value of the parameter as indicated on the y-axis. Here we again fit to R ij n = Anij/|A|, and take the overall normalization to be |A| = n,ij We can easily see by the color that the spread roughly decreases with the expected 1/ √ N scaling as the number of events is increased. The true value is now indicated by the solid black line, which as can be seen sits within the red region indicating that in most of the pseudo experiments the fit procedure returns a value of the parameter close to the true value. Quantifying more precisely how the spread, or more accurately the error, changes as a function of number of events requires the inclusion of detector effects and is beyond the scope of this study, but a more detailed analysis (using CMS criteria) is left to ongoing/future work.
One of the interesting questions to ask, is whether the golden channel is sensitive to the Zγ and γγ couplings of ϕ assuming it is the recently discovered resonance at ∼ 125 GeV. Since it has been firmly established that this resonance couples to ZZ through the ZµZ µ operator with a strength consistent with the SM prediction [36] it may per-haps be difficult to extract the Zγ and γγ couplings since they only occur through higher dimensional operators and will have couplings ∼ O(10 −2 − 10 −3 ), thus suppressing the partial widths corresponding to those operators in Figs. 4-5. Determining whether this is in fact impossible requires a detailed analysis including detector effects which is beyond the scope of this paper and we leave it to a future study.

VI. CONCLUSIONS AND OUTLOOK
In this study we have completed the first stage in the construction of a comprehensive analysis framework which builds upon earlier studies [13,17,27] and is aimed at extracting as much information as possible from the Higgs golden channel. First we extended previous analytic calculations for both signal and background in the 2e2µ Higgs 'golden channel' to include the 4e/4µ final states with the interference between identical final states. We have presented an overview of the calculations of the expressions as well as showing various doubly differential projections and relative 'partial fractions' for every combination of tensor structures.
We have also shown the potential of using these analytic expressions to perform parameter extraction of the various couplings of a spin-0 scalar to neutral electroweak gauge bosons including any correlations between parameters by implementing them into a maximum likelihood analysis. In order to show the validity of our maximization procedure we have focused on a simplified generator level analysis which includes both signal and background at fixed √ s. As our example parameter point, we have performed a simultaneous extraction of all parameters assuming real couplings (and overall normalization) of our scalar to ZZ, Zγ, and γγ pairs as well as the background fraction. We have shown that our maximum likelihood analysis gives accurate extraction of the parameters as well as the background fraction.
A more accurate analysis of course requires the inclusion of detector and systematic effects. We have not addressed these issues here and instead have left them for a series of accompanying studies of the golden channel [30,31] where we also demonstrate the advantage of analytic expressions when including detector effects. We have also neglected the use of 'production variables' into our set of observables since this requires careful treatment of the production mechanism which is beyond the scope of this study and furthermore would introduce additional systematic uncertainties. Since we fit to ratios of couplings and do not attempt to extract the overall normalization however, our results and analysis procedure are not overly sensitive to the production mechanism. We hope to include a detailed description of the production mechanism in future studies. In addition, we hope to conduct a detailed comparison between the sensitivity of the 2e2µ and 4e/4µ final states for different choices of the phase space cuts in order to determine the optimal phase space for extracting particular couplings to neutral electroweak gauge bosons.
In summary, we have demonstrated the potential of using analytic expressions in the golden channel to extract the couplings of a spin-0 scalar to neutral electroweak gauge bosons and have completed the first stage in the construction of a comprehensive analysis framework aimed at maximizing the power of this channel. This framework can now readily be adapted to include the relevant detector effects as well as any systematic uncertainties.  Here we examine the correlations between pairs of parameters. We conduct a large set of pseudo experiments with 1000 events for each and assuming a pure signal sample. The true value is indicated by the intersection of the two solid black lines and again we have fit to the ratios R ij n = Anij/|A| and take the normalization to be |A| = |A1ZZ |. In the additional plots we compare the results of the parameter extraction assuming a pure signal sample (black) to those which include both signal and background (red). We fit to the ratios R ij n = Anij/|A| and take the normalization to be |A| = n,ij |Anij| 2 .

VII. APPENDIX
In this Appendix we examine the 'partial fractions' of the various pairs of tensor structures which are found in Eq. (1). We also display a number of doubly differential spectra for a standard model signal as well as the qq → 4 background. Finally, we also show our validation of the signal and background calculations for the matrix element squared.

A. Relative 'Partial Fractions' for CMS Cuts
The total decay width for ϕ → ZZ + Zγ + γγ → 4 can be decomposed into the various 'partial widths' formed out of pairs of tensor structures in Eq.(2) (or operators if interpreted in terms of Eq. (3)). Since each term will be quadratic in the couplings, we can label each partial width by the appropriate combination of couplings AnijA * nīj . They are obtained by integrating the fully differential decay width in Eq. (19) over the kinematic variables defined in Sec. II. We then normalize these partial widths to the standard model value to form the various 'partial fractions'.
We show in Fig. 10 a table of these partial fractions for every possible combination of AnijA * nīj which can contribute to the 2e2µ decay width. For these partial fractions we now take a 'CMS-like' phase space of 40 GeV M1, 12 GeV M2 and √ s = 125 GeV as well as p T > 20, 10, 7, 7 GeV for the ordering of final state lepton pT and |η | < 2.4 for their rapidity. The couplings Anij have been separated into their real and imaginary components as Anij = AnijR + iAnijI and we have set all AnijR,I = 1. All of the |AnijR,I | 2 terms sit along the diagonal with the various interference terms making up the off-diagonal terms. Note that many of the interference terms are negative indicating destructive interference between the corresponding tensor structures (or operators).
In Fig. 11 we show the same plot for the 4e final state. One can see the change in the partial fractions and in particular the ZZ/γγ interference terms are significantly larger than in the 2e2µ channel. The blank entries indicate terms which are identically zero after integration. We can see that these entries are those for which CP violation in the form of interference between A1,2 and A3 tensor structures would occur. This is indicative of the fact that after one integrates over the kinematic all information on CP violation is lost. Of course for the fully differential decay width many of these terms are non-zero in principal allowing for sensitivity to CP violation in the golden channel.
Since all couplings are set to one, these tables essentially show how much each combination of tensor structures contributes to the ϕ → 4 phase space relative to the contribution from the partial width for which we have set A1ZZ = 2 and all other couplings to zero. From these values of the relative partial fractions, one can gain some insight into which combination of operators the golden channel might be most sensitive to. Furthermore, for a specific model one can take the prediction for the values of the various couplings and simply multiply by the numbers given in Fig. 10-11 to get a feel for whether those couplings might be probed in the golden channel.
For most realistic models, all couplings apart from A1ZZ are generated by higher dimensional operators and are expected to be small. In Figs. 12-13 we also show the same tables for the standard model prediction including the Zγ and γγ couplings for which we have A1ZZ = 2, A2Zγ 0.007, A2Zγ −0.008 [45] 7 while all other couplings zero. These values are normalized the same as in Figs. 10-11. Of course, these rates do not contain information about the shapes in the various distributions so in principal the sensitivity is greater than might be inferred from these values. Whether or not the golden channel has sensitivity to these couplings in the standard model requires careful study, which we leave for ongoing work.
In Figs. 14 and 15 we show the integral of the absolute value of the differential decay width. This gives a better indication of the shape differences in the different combinations of operators since some of them can integrate to zero when the absolute value is not taken. Furthermore we can see in this table some of the potential sensitivity in the golden channel to CP violation. Note that there are two sources of CP violation which occur. One is due to the interference between the A1,2 and A3 tensor structures, while the other occurs in the interference between the real and imaginary components of the couplings from different tensor structures.
One could also imagine attempting to find different sets of cuts in order to maximize the contribution of a particular combination of operators. In addition, the sensitivity between the 2e2µ and 4e/4µ final states may differ depending on the phase space that is chosen. We leave a detailed investigation of this issue to future work. These tables, however, obviously only give a partial picture of the sensitivity to the different operator combinations and are meant to be used only as a guide. Of course when performing parameter extraction the full kinematic information of the differential decay width is used. Obviously, for a scalar resonance with a mass much larger or smaller than 125 GeV these numbers can change significantly.

B. Doubly Differential spectra
In Fig. 16-17 we show various combinations of the doubly differential spectra for both the signal and background in     [45] in the 4e/4µ final state. For these partial fractions we take a 'CMS-like' phase space of 40 GeV M1, 12 GeV M2 and √ s = 125 GeV as well as p T > 20, 10, 7, 7 GeV for the ordering of final state lepton pT and |η | < 2.4 for their rapidity. They have been normalized to the partial width where we take A1ZZ = 2 and all other couplings zero. the 4e/4µ final state. These are primarily for illustration purposes, but from these one can get an idea of the correlations between the different kinematic variables. One can also see from these spectra the strong discriminating power between the signal and background in the golden channel. For the signal plots in Fig. 16 we only show the standard model result for which only A1ZZ is non-zero 8 . The background spectra are shown in Fig. 17. For all distributions the phase space is defined as 4 GeV < M1 < 120 GeV and 4 GeV < M2 < 120 GeV with √ s = 125 GeV for signal and background. We also take |η | < 2.4 and p T > 2 GeV for the lepton rapidity and transverse momentum. For these distributions we show the (M1, M2), (M1, Φ1), (M1, Φ), (M2, Φ1), (M2, Φ), (Φ1, Φ) doubly differential spectra.

C. Validation of Calculations
In this section we show a validation of the analytic calculations for the golden channel signal and background. Both the signal and background are validated against the Madgraph result for the leading order matrix element squared for a large number of random phase space points. For these comparisons we have generated 100k random phase space points in the range 5 GeV ≤ √ s ≤ 1000 GeV so these expressions are valid for essentially any scalar mass and energy range. We show the validation for the 4e/4µ final state, but as discussed above this is also validates the 2e2µ final state (though it was also explicitly validated in [33]) which is simply one term in the 4e/4µ matrix element squared.
We first show in Fig. 18 the validation for the ϕ → 4 calculation of the matrix element squared obtain in Eq. (18). We show the validation in two ways. In the top plot we show the Log(|M | 2 ) for a large number of random phase space points and plot the two results on top of one another. The Madgraph result is shown in red while the analytic result is shown in yellow. The two results are indistinguishable from one another and thus the two distributions sit on top of each other leading to the orange color seen.
In the bottom plot we show the fractional difference in their matrix elements squared for the same set of phase space points. The agreement is perfect up to very tiny differences due to numerical precision when evaluating the matrix elements squared for specific phase space points. To obtain the matrix element squared from Madgraph we have implemented the Lagrangian in Eq.(3) (or equivalently vertex in Eq.(1)) into the FeynRules/Madgraph [34,35] framework. We have chosen all of the Anij couplings to have random non-zero values for both their real and imaginary parts. Thus the complete expression including all tensor structures in Eq.(1) and their interference has been validated. One can now easily obtain any expression which includes a subset of the possible tensor structures by simply setting the unwanted Anij to zero.
In Fig. 19 we show the same validations for the leading order qq → 4 background again validated against the Madgraph result. Again we see essentially perfect agreement. See also [33] for how the different components of the qq → 4 depend on √ s. We also provide there an analytic expression for In the bottom plot we show the fractional difference in their matrix elements squared for the same set of phase space points. The tiny differences seen are due to numerical precision when evaluating the matrix elements squared for specific phase space points. In the bottom plot we show the fractional difference in their matrix elements squared for the same set of phase space points. The tiny differences seen are due to numerical precision when evaluating the matrix elements squared for specific phase space points.