Tripartite entanglement and Bell non-locality in loop-induced Higgs boson decays

In this article, we study quantum entanglement properties of the three-body $H\to\gamma l\bar{l}$ decays (for $l=e,\mu,\tau$) within the context of the Standard Model augmented with CP-violating interactions in the lepton Yukawa sector. Our aim is to elucidate the distribution of entanglement between the final photon, lepton and antilepton across the phase-space. These rare Higgs boson decays occur at 1-loop level, presenting a unique opportunity to scrutinize quantum correlations of fundamental interactions in tripartite systems by computing concurrence measures and investigating Bell non-locality. Moreover, we explore post-decay and autodistillation phenomena. Multipartite entanglement measures have much richer structure than those in the bipartite case, thus deserve more attention in collider phenomenology. In this line, we analyze here novel observables for these three-body Higgs boson decays, which can be extended to other multiparticle systems within the high-energy regime. We found that entanglement manifests among final particles, occasionally achieving a maximally entangled state in specific kinematical configurations. Also, these decay channels are promising for Bell non-locality tests but CP-effects are suppressed by lepton masses in this kind of observables.


Introduction
Quantum entanglement stands as a pivotal resource for tasks which can not be performed via classical resources.Quantum Information theory develops the manipulation, control and distribution of the entanglement in a given system, with applications ranging from cryptography [1] to teleportation [2] and quantum computation [3].In particular, entanglement can be generated when two systems interact and elementary particle collisions, described by Quantum Field Theory (QFT), provide a natural framework for studying such properties of fundamental interactions.However, it receives a very recent attention in the high-energy physics (HEP) community, see for instance review [4] and references therein.Notably, the ATLAS Collaboration observed entanglement in t t production [5][6][7][8][9][10][11] with a significance exceeding 5σ [12], despite collider detectors were not initially designed for probing such properties.
The S-matrix formulation in QFT allows to compute decay and scattering amplitudes at a given order in perturbation theory.Radiative corrections, arising from closed loops in the propagation of virtual particles, represent genuine quantum effects and the inhered correlations between initial and final states warrant attention.
On the other hand, tripartite entanglement within the HEP context is in development: QED scattering processes with spectator particle were analyzed in [26][27][28], Ref. [29] studies heavy fermion decaying into three fermions via generic (pseudo)scalar, (pseudo)vector and (pseudo)tensor interactions, Ref. [30] explores three-flavor entanglement in neutrino oscillations, and Ref. [31] addresses entanglement among the two spins and total angular momentum in H → ZZ, W W .General properties of multipartite systems were presented in [32,33].Furthermore, the tripartite Higgs boson decays considered here can be related to bipartite ones, in order to investigate post-decay entanglement [34] and autodistillation phenomena [35,36].
Beyond the entanglement due to correlations among constituents of a system, is the concept of non-locality.The advantage of quantum mechanics over classical theories for certain information processing tasks lies precisely in the non-locality of quantum correlations.On the contrary, local realistic (LR) or hidden variable (LHV) theories are described by local objective properties that are independent of observation.LR assumption has experimental consequences providing constraints on the statistics of two or more physically separated systems through Bell inequalities [37].These inequalities can be violated just by predictions of quantum mechanics.The structure of non-local correlations is much richer (but also less understood) for multipartite systems than for bipartite ones.In particular, there exist different notions of non-locality by extending the bipartite definition, see for instance [38].Bell inequality tests to these three-body Higgs decays are also explored in this work.
Since the Higgs boson discovery in 2012, the determination of its properties is part of the major experimental program of ATLAS and CMS Collaborations.In the Standard Model (SM) context, the H → γl l decays were discussed and related to two-body H → γγ, γZ decays in [39][40][41][42].These decays were also examined within various beyond SM theories [43][44][45][46][47].There are proposals to study CP properties of the Higgs boson in these three-body decays via the forward-backward asymmetry [48][49][50] and other polarization-dependent observables were introduced in [51,52].This work analyzes novel observables for these three-body decays and extends the understanding of quantum interactions within such systems.Concerning the experimental searches, the 13 TeV data analyses for the decay of a Higgs boson in the γl l channel were performed by CMS [53] and ATLAS [54].
Despite the fact that we do not delve into experimental aspects, this work provides a new conceptual analysis of entanglement and non-locality properties of multipartite states produced in loop-induced Higgs decays.Predictions of related quantities from Monte-Carlo simulations in complete collider events lie beyond the scope of this work.Consequently, the aim of this paper is to identify kinematical regions of the phase-space where interesting quantum mechanical measurements might be performed.To our knowledge, this represents the first analysis of tripartite entanglement in a full 1-loop SM computation.It is worth noting that the application of this analysis at the detector level would require polarization measurements of the final high-energetic photons (see for instance a related discussion in [13]).Although this kind of measurement is not currently available in ATLAS and CMS detectors, in contrast to the case of massive gauge bosons, the LHCb Collaboration performed analysis for photon polarization in b-baryon decays [55].There are also proposals to study CP properties of the Higgs boson through the di-photon decay [56,57].
The paper is structured as follows: Section 2 provides a comprehensive overview of the 1-loop computation of the helicity amplitudes corresponding to H → γl l decays.In Section 3, we introduce the tripartite density matrix formalism along with the measures associated with entanglement and Bell non-locality.Analytical results can be found in this section.The Section 4 is devoted to a comparative analysis of entanglement properties of these three-body decays with the related two-body H → l l, γγ, γZ decays.The numerical results of this work are collected in Section 5 for each taulepton, muon and electron cases.We then summarize the main findings and outline future perspectives in Section 6. Appendices contain detailed information for the helicity amplitude computation and show the main numerical results using Dalitz plots representation.

Three-body Higgs boson decays
The present computation is performed in the context of the SM with additional CP-violating terms in the lepton Yukawa sector.Concretely, we consider a generic interaction between the Higgs boson to each lepton l = e, µ, τ as where Y l = m l /v is the Yukawa coupling, with m l as the lepton mass and the vev of the Higgs field v = 246 GeV.The magnitude of this interaction and the CP-phase are parametrized by κ l CP ∈ ℜ + and δ l CP ∈ [0, 2π].The SM is recovered for κ l CP = 1 and δ l CP = 0. We are interested in quantum entanglement properties of the spin degrees of freedom in the rare H → γl l decays, then we computed the corresponding helicity amplitudes.We denote the four-momenta of the photon, lepton and antilepton by k, p − and p + and their helicities, along the direction of motion, as s 1 , s 2 and s 3 .In addition, ε s1 = ε s1 (k) is the polarization vector of the photon and u s2 = u s2 (p − ), v s3 = v s3 (p + ) are the lepton and antilepton spinors (conventions collected in Appendix A).The generic diagrams corresponding to H → γl l are depicted in Fig. 1.Diagrams (a) and (b), in the first row of this figure, correspond to the photon emission process at O(ℏ 0 ), i.e. leading order (LO) in perturbation theory.This tree level contribution is suppressed by one power of the Yukawa coupling and the resulting helicity amplitude is where the global factor A Tree is em l /v (e as the electromagnetic coupling constant).On the other hand, diagrams corresponding to the electroweak 1-loop O(ℏ) contribution, i.e next-to-leading order (NLO), are schematically represented in the second row of Fig. 1.The gray blobs denote the one-particle-irreducible (1PI) Green functions renormalized in the on-shell scheme.Diagrams (c) and (d) correspond to photon emission, diagram (e) represents the 4-legs 1PI (usually called boxes), and (f) has H → γV * as two-body intermediate Higgs boson decay (with V = γ, Z).In this O(ℏ) computation, the usual linear covariant R ξ -gauge is implemented for the bosonic loops.For the fermionic loops, only the top-quark Yukawa coupling is considered and the rest are neglected.In particular, O(ℏ) lepton mass effects are relevant just close to the lepton-pair production threshold and we avoid this region considering dilepton invariant mass above 0.1m H [41,42,49]. Hence, there is not suppression with Y l nor CP-effects in this contribution (in contrast to the tree level).Also, we have vanishing renormalized 1PI corresponding to 2-legs Green functions Hγ , HZ , Hϕ Z and γZ , and 3-legs Green function Γ Hγϕ Z in this setup1 .Hence the resulting helicity amplitude is written as where P L,R = (1 ∓ γ 5 )/2 and the form factors a 1,2 and b 1,2 are functions of the momenta k, p − and p + .Notice that electroweak (EW) 1-loop corrections involve the Higgs coupling to bosons and to the top-quark, then yield a nonzero decay amplitudes even for vanishing lepton Yukawa.The full computation using the R ξ -gauge of all diagrams represented in Fig. 1 in the setup described previously, was performed in [49].In particular, the form factors were provided in the corresponding ancillary file, which is implemented for this work.
The diagrams (f) in the second row of Fig. 1 correspond to the two-body intermediate Higgs decay H → γV * → γl l.In that case, we have a 1 = a 2 and b 1 = b 2 since lepton momenta are combined in order to get the intermediate gauge boson momentum q = p − + p + .The results for the renormalized 3-legs Green function Γ HγV * in the R ξ -gauge can be found in [58][59][60].Notice that the predictions corresponding to H → γV * → γl l are very similar for the three flavors because lepton masses are only present in the spinors (leading to negligible effects far to the threshold), there is no tree level Hl l interaction, lepton Yukawas in the loops are neglected and we have universality in V l l interaction.Particular attention is devoted to the resonant production of the Z boson, for which a non-vanishing decay width Γ Z is considered by means of a Breit-Wigner distribution in the subprocesses type (f) [42,49].Other resonant productions decaying into a lepton-pair correspond to the quarkonium states J/ψ and Υ(nS), but they are also rejected by imposing a lower bound (0.1m H ) to the dilepton invariant mass [53,54].
Furthermore, it is interesting to consider the hybrid computation [48,50] in which the tree level and the two-body intermediate Higgs decay are combined, i.e. diagrams (a), (b) and (f) of Fig. 1.In summary, the three considered computations in this work are: (2.4) Helicity amplitudes in Eqs.(2.2)-(2.3)are Lorentz invariant and they can be written in terms of the Mandelstam variables s = (p − + p + ) 2 , t = (k + p − ) 2 and u = (k + p + ) 2 , which satisfy the relation s + t + u = m 2 H + 2m 2 l .However, a convenient choice of the reference frame simplifies the resulting expressions.For the present computation, we consider the rest frame of the lepton-pair where the z-axis is along the direction of the lepton, the y-axis is perpendicular to the decay plane and the photon momentum has positive x-component.Then the two relevant kinematical variables are the dilepton invariant mass (m l l = √ s) and the polar angle between photon and lepton (θ γl ).Details of the kinematics are gathered in Appendix A.
The differential decay width for H → γl l is computed as (2.5) The last sum involves the 8 helicity amplitudes (see Appendix B).The polar angle θ γl belongs to [0, π] and the Mandelstam variable s has relevant limits where a lower cut E γ cut to the photon energy is imposed in order to avoid infrared (IR) divergences [49].We apply E γ cut = 1 GeV, in the Higgs boson rest frame, for all numerical evaluations in this work and no additional cuts over the final particles are applied.In addition, we consider the following input SM parameters (2.7) The resulting differential decay width respect to the dilepton invariant mass for the three flavors are shown in Fig. 2. It is assumed that κ l CP = 1 and δ l CP = 0, i.e. the SM values.Solid lines correspond to the computations in Eq. (2.4): full (black), hybrid (orange) and two-body intermediate decay (green).Blue and yellow dashed lines account for the tree level in Eq. (2.2) and complete 1-loop (diagrams (c)-(f)) in Eq. ( 2.3), respectively.In addition to the relevant range of Eq. (2.6), the range 2m l ≤ m l l ≤ 0.1m H is included in shaded region for illustrative purposes (just EW 1-loop corrections are considered and we are excluding the contributions associated to J/ψ and Υ(nS) resonances).See also Table 1 in Appendix B for numerical estimations within this setup.
Taulepton case A general conclusion from Fig. 2 of the full computation is that the interference between tree level and 1-loop contributions is negligible, or in other words, either dominates tree level or dominates 1-loop.The only exception is the tau-lepton case (upper plot) near the Z-pole peak, where the boson resonance of the 1-loop contribution interferes constructively with the tree-level.In particular, the tree level contribution strongly dominates for energies above 30 GeV in this tau-lepton case.Regarding the muon case (left-lower plot), the 1-loop controls the behavior up to 100 GeV and the tree level does in the high dimuon invariant mass.The electron case (right-lower plot) is clearly dominated by the 1-loop contribution in the whole range.By construction in Eq. (2.4), the hybrid computation (orange solid lines) are very close to the blue dashed lines when the tree level dominates but they are close to the two-body intermediate (green solid) when the tree level is negligible.

3-qubit formalism
We are interested in quantum entanglement properties of the spin degrees of freedom corresponding to the final particles in the process H → γl l.Since photons have two transverse polarizations and leptons have two helicities along the momentum direction, they correspond to qubits in the Quantum Information language.This 3-qubit system is described by the pure state |ψ⟩, which is expanded using the helicity amplitude basis {+, −} ⊗ {+, −} ⊗ {+, −}.Therefore, the associated 8×8 density matrix ρ = |ψ⟩⟨ψ| for a 3-qubit system, in terms of the helicity amplitudes of the process, is where the first factor in the r.h.s is the total unpolarized square amplitude and determines the normalization Tr[ρ]=1.
An important question is to recognize and quantify the entanglement in a given quantum state.Concurrence is one of the well-defined quantitative measures for entanglement/separability criteria [61,62].However multipartite systems have richer structure than bipartite ones.In particular, just the so-called genuine separability is a direct generalization of the bipartite separability and there are many types of partial separability.Following [29,32] to this 3-qubits system, just one bipartition (among 1-2 or 2-1) is relevant since C i|jk = C jk|i .Then we only consider the one-to-other concurrences where ρ jk is the reduced density matrix of subsystem jk by tracing over particle i, i.e.
The relevance of this quantifier is that a state described by ρ is biseparable if and only if C jk|i = 0.
The one-to-other concurrences define the concurrence triangle and the corresponding area represents a measure of the genuine entanglement of this 3-qubit system (GTE) by computing where S = (C 23|1 + C 31|2 + C 12|3 )/2 is the semiperimeter of the concurrence triangle [63].Genuine multipartite entanglement arises in a quantum state when it cannot be expressed as a convex combination of biseparable states, or equivalently, the system is entangled respect to all bipartitions of the parties.
In addition, the entanglement between two individual particles is evaluated by the one-to-one concurrences C jk , which are obtained from the eigenvalues η jk 's (in decreasing order) of the matrix Next, we can analytically compute the previous quantifiers for the final state of the H → γl l decay.The 8 helicity amplitudes are written in terms of the generic 1-loop form factors and are collected in Eq. (B.1).Remember that we impose a lower limit in Eq. (2.6) to the dilepton invariant mass, then m l ≪ √ s regime is satisfied for electron and muon cases, whereas tiny tau-lepton mass effects are present near this lower limit.Hence, in the high dilepton invariant mass regime respect to m l , we can neglect lepton masses (except in the A Tree factor) and a compact expression combining Eqs.(2.2)-(2.3)can be written for the pure final state: where the normalization factor is In the vanishing lepton mass limit, the tree level contribution has an IR divergence when the photon is collinear to lepton or antilepton (θ γl = 0, π), and the upper limit in Eq. (2.6) must be also imposed.Keeping lepton mass, the collinear configuration is well-defined, as shown in Eq. (B.2).Notice that each helicity amplitude receive contributions either from tree level or 1-loop in the massless lepton regime.In particular, interference terms appear proportional to lepton masses and can be considered negligible.
Furthermore, the previous quantifiers in Eqs.(3.2)-(3.4)can be computed analytically in some particular cases.Considering just the tree level, i.e.Eq. (2.2), a direct computation yields to with c γl = cos(θ γl ), c CP = cos δ l CP and we keep m l terms.Notably, the tree level concurrences of the lepton and antilepton respect to the other particles attain the maximal theoretical value 1 independently on the kinematics.In that case, the F 3 measure of Eq. (3.3) is very similar to C l l|γ since they have values in the interval [0, 1].
On the other hand, the photon-to-dilepton concurrence dependence on the CP-phase δ l CP is suppressed by the lepton mass, and then negligible in the relevant energy range of Eq. (2.6).Also, the global factors A Tree and κ l CP are cancel out due to the normalization of the state ρ2 .In particular, C Tree l l|γ never vanishes but in the high energy regime it simplifies to which is very close to zero in the cut energy √ s cut and we almost have the biseparable state where the normalization factor is omitted.On the contrary, if the photon is collinear with lepton or antilepton, C Tree l l|γ reaches the maximal value 1.
Now considering just the 1-loop contribution and neglecting lepton mass terms, the one-to-other concurrences are We found that the 1-loop concurrences of the lepton and antilepton respect to the other particles never vanish and achieve the maximal value 1 if a 1 = a 2 and b 1 = b 2 (as in the two-body intermediate decay).In that case, the photon-to-dilepton concurrence vanishes when c γl = 0 and we have the biseparable state where the normalization factor is omitted.
In general, considering both tree level and 1-loop contributions, the conditions for genuine entanglement using the concurrence vector formalism to this 3-qubit system [32] are {q 0 , q 1 , q 2 } = {0, 0, 0}, where for which we neglected lepton mass terms and used Eq.(3.5).Of course, our previous findings in Eqs.(3.7)-(3.9)for either tree level or 1-loop contribution are recovered from these conditions.The q 0 = 0 and q 1 = 0 equations relate the form factors to each other, and q 2 = 0 establish a relation between them and the tree level factors.For each kind of computation (full or two-body intermediate decay), we have non-trivial dependence of the form factors with s and cos(θ γl ), and very particular kinematical configurations could correspond to biseparable states (the numerical analysis is developed in Section 5).
Concerning the one-to-one concurrences, they require a 4×4 diagonalization and analytical expressions result just by neglecting lepton mass terms.In that regime, both photon-to-lepton C γl and photon-to-antilepton C γ l exactly vanish for both tree level and 1-loop contributions (then lepton mass terms in the spinors will be relevant in the numerical analysis) and the lepton-to-antilepton C l l for the tree level is very compact (3.11) Observe that C Tree l l never vanishes and is very close to 1 in the cut energy √ s cut .The expression corresponding to the 1-loop contribution is not illuminating and it is omitted, however when the photon is collinear with the lepton or antilepton, it vanishes.The previous two classes of entanglement (one-to-other and one-to-one) are inequivalent.In fact, the i-to-jk entanglement limits the entanglements i-j and i-k by means of the Coffman-Kundu-Wootters (CKW) monogamy inequality [64] where the rhs corresponds to the three-tangle measure, which is the same for all permutations of subsystem indices.For the present rare decays, Eqs.(3.7)- (3.11) show that the three-tangle is essentially the same as the square of the C l l|γ concurrence.
Finally, regarding the Bell non-locality, local realism (LR) had experimental consequences that can be tested through the inequality of the expectation value of certain Bell operator B: where β LR is the locally realist bound and, if an experiment shows ⟨B⟩ > β LR , then there is a non-local communication between different particles of a composed system.For the 2-qubit case, the Clauser-Horne-Shimony-Holt (CHSH) operator [65] yields to β LR = 2 and it is the only relevant (optimal) Bell operator for this case.
In the multipartite case, non-locality displays a more complex structure, then much richer, and the characterization of multipartite non-local correlations results a challenging problem.In particular, there exist different notions of non-locality arising as extensions of the bipartite definition [38].In this work, we implement the 3-qubit Mermin operator [66] with âi = ⃗ a i • ⃗ σ, bi = ⃗ b i • ⃗ σ and ĉi = ⃗ c i • ⃗ σ are Hermitian operators acting on the Hilbert spaces of the photon, lepton and antilepton, respectively.The maximum value of this operator in a local realistic theory is ⟨M 3 ⟩ LR = 2, whereas the quantum maximum value is ⟨M 3 ⟩ QM = 4.There are other options for the Bell operator in the literature, as for example S 3 proposed by Svetlichny in the seminal work [67].Notice that Mermin inequality can be violated by biseparable states.This is in contrast to the S 3 case [38,68], which yields to a stronger inequality but it is not a necessary requirement for genuine tripartite non-locality.In addition, the ratio ⟨B⟩ QM /⟨B⟩ LR is greater in the Mermin case.We found that both operators essentially have the same behaviour for these three-body Higgs boson decays.We also verified that some final state configurations result in predictions compatible with LR using S 3 in Eq. (3.13) but violates this inequality using M 3 .Then we just present the numerical results of the Mermin operator in Section 5.
Observe that the expectation value of the Mermin operator in Eq. (3.14) involves a maximization over the six directions corresponding to ⃗ a 1,2 , ⃗ b 1,2 and ⃗ c 1,2 .As far as we know, there is no closed analytical way to determine this maximum for an arbitrary 3-qubit system, as the 2-qubit case using the Horodecki condition [69] for the CHSH operator.

Related bipartite Higgs boson decays
In this section, we compare the entanglement properties of our process of interest H → γl l respect to related bipartite Higgs boson decays, that were previously studied in the literature.On the one hand, the photon emission process, diagrams (a)-(b) in Fig. 1, can be connected with H → l l.On the other hand, diagrams type (f) are associated to the H → γV processes (V = γ, Z).We will relate the one-to-one and one-to-other concurrences of the previous section with those of the bipartite H → l l and H → γV decays, respectively.Some comments are in order: the dilepton decay channel is a well defined observable and H → γl l is part of its NLO corrections.However, H → γV * → γl l is unphysical (in particular, gauge-dependent) and it can be extracted as pseudo-observable by imposing kinematical cuts on the dilepton invariant mass [41,49,53,54].

H → l l and CP-effects
Entanglement properties of H → l l were analyzed for the tau-lepton case in [13,14].The CPviolating interaction, as in Eq. (2.1), was also treated in these references.The authors conclude that the concurrence of the ditau system is maximal regardless of the CP-phase and its determination was obtained by a direct fit of the entries in the correlation matrix of the 2-qubit final state.In addition, CP properties of the lepton Yukawa sector were measured for the H → τ τ decay by CMS [70] and ATLAS [71] and there are proposal to study them through the forward-backward asymmetry in H → γl l decays [48,50].Now we focus on the tree level contribution since the 1-loop does not depend on κ l CP nor δ l CP .The reduced density matrix ρ l l of the dilepton subsystem, after tracing over the photon helicity, can be decomposed as where the coefficients A and B are the spin polarizations and C is the spin correlation matrix of the resulting bipartite subsystem.For these rare decays, neglecting the lepton mass, the A and B coefficients vanish and the correlation matrix results The concurrence of this dilepton subsystem was presented in Eq. (3.11).The dependence on the CP-phase in the correlation matrix is missing in this entanglement measure, but sensitivity to δ l CP could be obtained by a direct fit of the entries in the correlation matrix, as in the case of H → τ τ .
We are analyzing novel observables for these tripartite Higgs boson decays.In the previous section, we found that CP-effects are suppressed by the lepton mass for the one-to-other concurrences in Eq. (3.7).In particular, for the relevant dilepton invariant mass range of Eq. (2.6), they are roughly independent of κ l CP and δ l CP .In other words, this kind of quantifiers are not sensitive to the new physics introduced by these parameters.Hence, we assume κ l CP = 1 and δ l CP = 0, i.e.SM computation, for the rest of this manuscript.

H → γZ and post-decay entanglement
Other related bipartite decays are H → γγ and H → γZ, associated to diagrams type (f) of Fig. 1.From the experimental point of view, the H → γγ was one of the golden decay channels for the Higgs boson discovery.The latest ATLAS and CMS measurements of the Higgs boson mass in this channel are presented in [72,73] and the corresponding cross-section can be found in [74,75].On the other hand, a recent combined analysis of ATLAS and CMS found evidence for H → γZ decay [76], which agrees with the SM theoretical expectation within 1.9 standard deviations.
Entanglement properties of the diphoton channel was previously studied in [13].This final state is maximally entangled, i.e. the concurrence attains the theoretical maximum 1.Furthermore, the corresponding CHSH operator saturates the Cirelson bound [77] and the Bell inequality is maximally violated.
On the other hand, requiring dilepton invariant mass close to the Z-pole peak, the pseudoobservable H → γZ * can be enhanced.This final state corresponds to a bipartite system composed by one qubit and one qutrit.This state was studied in [25] when it is coming from vector boson scattering and, as far as we know, this is the first time that entanglement properties of H → γZ are studied.The SM amplitude of this decay is expressed as where {k µ , s 1 } and {q ν , s1 } are the momentum and helicity of the photon and Z boson, ε's are their polarization vectors and the 3-legs form factor V HγZ accounts for the 1-loop contribution (that only depends on k and q momenta).The precise dependence of this form factor on the momenta is not relevant for the following discussion.In the {+, −} ⊗ {+, 0, −} basis, the resulting density matrix is 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/2 1/2 0 0 0 0 1/2 1/2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 which is independent of the form factor (since it is a global factor of the helicity amplitudes).The concurrence for a bipartite system with arbitrary dimension d 1 ⊗ d 2 was defined in [78] and the theoretical maximum is 2(d − 1)/d with d = min{d 1 , d 2 }.For this bipartite decay, the concurrence achieves the maximum 1 and the γZ state (coming from the Higgs boson) results maximally entangled.However, the 2⊗3 generalized CHSH operator [16,25,79] for this decay never exceeds 2.
The H → γZ and our process of interest H → γl l allow to test the post-decay entanglement [34] and autodistillation phenomena [35,36].Concretely, the concurrence C initial of an 'initial' bipartite state {a, ā} is computed.The particle ā decays along the process ā → b 1 b 2 and the entangled 'initial' state led to the entangled 'final' subsystems {a, b 1 } and {a, b 2 } with concurrences C final 1 and C final 2 , respectively.A priori, different amount of entanglement is expected and the autodistillation phenomena arises when final concurrences after decay are greater than the initial one.In our context, we consider the two-body intermediate Higgs decay into the 'initial' state γZ and the subsequent post-decay Z → l l.Diagrams type (f), with Z boson as mediator, contribute to the 1-loop amplitude in Eq. ( 2.3) as where the intermediate Z boson propagator was written in terms of the polarization vectors and q = p − + p + .Close to the Z-pole, the Narrow-Width Approximation (NWA) can be applied, which mainly replaces the denominator's contribution of this propagator to the concurrences by a delta function, such that only Z boson on-shell (OS) effects remain.In that case, the previous helicity amplitude is splitted as where M s1 s1 H→γZ was given in Eq. ( 4.3) and the SM polarized decay amplitudes for Z → l l are Following [34][35][36], we consider the subsystems γl and γ l, by tracing over l and l, and compute the one-to-other concurrences C γl| l and C γ l|l .Notice that the 'initial' state γZ, described by Eq. (4.4), is maximally entangled (concurrence equals to 1) and the original autodistillation phenomena is not relevant now.However, we arrive to the one-to-other concurrences which are very close to the corresponding one of the γZ state, and we conclude that the Z boson decaying into a lepton-pair does not significantly decrease the initial concurrence.Notoriously, if we impose the 'MaxEnt Principle' [80] to the previous one-to-other concurrences, i.e. demand C γ l|l = 1, the Weinberg angle should satisfy s 2 w = 0.25, which is surprisingly close to the SM value [81].

Numerical Results
In this section, we present the distributions of the entanglement and Bell non-locality measures in the [m l l, cos(θ γl )] plane of the phase-space.We analyze the three lepton families separately since the 1-loop contribution dominates different energy regimes in each case, and also corresponds to separate experimental channels.The numerical results correspond to the SM, i.e. κ l CP = 1 and δ l CP = 0 in Eq. (B.1), including lepton mass effects and using the corresponding form factors for each kind of computation in Eq. (2.4).Then, in the plots, the range for cos(θ γl ) is [0, 1] since the distributions result symmetric under cos → − cos, and the dilepton invariant mass covers the range in Eq. (2.6).
All the concurrences are defined in such a way that have values in the range [0, 1] and the present computation shows that the theoretical minimum and maximum are almost achieved in some particular configurations.For the distributions in Figs.3-5, the same color-scale is used where the purple(red) regions represent values close to 0(1).For each point of the phase-space, the CKW inequality of Eq. (3.12) was verified.Because CP-effects are negligible, the results in this section correspond to the SM computation.Hence C γ l|l and C γl| l are the same under transformation θ γl → π −θ γl , or equivalently cos(θ γl ) → − cos(θ γl ).The same occurs for C γ l and C γl .The general behaviour of these quantifiers were treated analytically in Section 3. In particular, we have simplified formulas of the concurrences in the limit m l ≪ 0.1m H ≤ m l l for both tree level and 1-loop contributions.Of course, the predictions of the 1-loop form factors for each kind of computation in Eq. (2.4) are determined numerically and provide the precise behaviour of these quantifiers in the phase-space.
Concerning the Mermin operator distributions, a green-scale is chosen for the resulting values in the range [1.8, 4].In particular, just small regions satisfy the Bell inequality of Eq. (3.13), compatible with local-realism, and we conclude that these rare Higgs boson decays are promising observables for test tripartite Bell non-locality in a high-energy regime.
In order to implement the common variables for three-body decays and to be independent of the reference frame, the corresponding distributions using the Dalitz plot representation are collected in Appendix C.

Tau-lepton case
The H → γτ τ process is dominated by the 1-loop (tree level) contribution for energies below (above) 30 GeV, as can be seen from the first row of Fig. 2. Around the Z-pole, both contributions interfere constructively but the resonance comes from the 1-loop diagrams type (f).The resulting one-to-other and one-to-one concurrences C l l|γ and C l l in the full computation are shown in the first row of Fig. 3, together with the Mermin operator M 3 in the [m τ τ , cos(θ γτ )] plane.
The one-to-other concurrence C l l|γ achieves values close to the theoretical maximum 1 in the red regions, that is roughly for i) energies m τ τ in the range [40,60]   The results corresponding to H → γV * → γµμ are the same as the tau-lepton case (second row of Fig. 3) since muon mass effects in this two-body intermediate decay are negligible.In particular, we have same entanglement distributions as in the second row of Fig. 3 for the considered invariant mass range of Eq. (2.6).
The distributions within the hybrid computation are shown in the second row of Fig. 4 and differ respect to the full computation (first row), in contrast to the τ -case in which both computations are very similar.In the present case, the hybrid predictions are a sort of transition between the full to the two-body intermediate decay computations.In particular, the red regions are reduced whereas the purple regions increase.
The M 3 distributions of the last column in Fig. 4 are also changed respect to the tau-lepton case.Now the full computation results in lower values and exhibit a small region satisfying Eq. (3.13) around the Z-pole and 0.6 ≤ cos(θ γµ ) ≤ 0.9 (with minimum ∼ 1.9).Notoriously, the hybrid computation exhibits its largest values in that region.Also, the maxima (almost 4) are located in lower-right corner, i.e. m l l ≤ 30 GeV with 0.9 ≤ | cos(θ γµ )|, in both computations.

Electron case
In this channel, the tree-level contribution is negligible since it is suppressed by the electron Yukawa.Then, the hybrid and two-body intermediate computations are practically identical, leading to distributions as in second row of Fig. 3.
The full computation distributions are shown in Fig. 5.The one-to-other concurrence C l l|γ reaches maxima (almost 1) for 0.8 ≲ | cos(θ γe )| with energies below (above) 40 (115) GeV, while the central region, | cos(θ γe )| ≲ 0.1, corresponds to the minimal values.The one-to-other concurrences C γ l|l and C γl| l are no longer homogeneous around 1, remember Eq. (3.9), and exhibit narrow bands in energy.
In particular, values close to 1 are still present at the Z-pole (as expected in the post-decay) and Finally, the Mermin operator distribution of the last plot in Fig. 5 shows values lower than 2 (with minimum ∼ 1.68) for energies ∼ 80, 100 GeV and 0.45 ≤ cos(θ γe ) ≤ 0.85.Now, the theoretical maxima 4 is reached in the lower-right and upper-right corners.

Summary and perspectives
This work aims to elucidate the mechanisms underlying the generation and distribution of entanglement within the framework of fundamental interactions described by the SM and incorporating CP-violation in the lepton Yukawa sector.Based on the concurrence and Bell operator definitions for tripartite systems, we explore entanglement and non-locality properties of H → γl l decays (for l = τ, µ, e).These three-body decays are Yukawa suppressed at leading-order, then electroweak 1-loop corrections are included.They offer a unique opportunity to examine quantum correlations arising at NLO in perturbation theory within the SM.This paper presents novel observables for these three-body Higgs boson decays and extend our understanding of quantum interactions within such systems.
Our goal was to identify regions of the phase-space where the final particles result entangled after the Higgs boson decay and to determine the feasibility of testing Bell inequality under these kinematical configurations.By expanding beyond traditional bipartite systems to the three-body final state, we compute various entanglement measures including one-to-other and one-to-one concurrences, the area of the concurrence triangle, the conditions for genuine entanglement of 3-qubit systems using the concurrence vector and the three-tangle measure.These quantifiers were derived analytically in terms of the generic 1-loop form factors and numerical predictions in the [m l l, cos(θ γl )] plane for three kind of computations were also explored.
We found that the final photon, lepton and antilepton result entangled after the Higgs boson decay, and this also holds by considering the one-to-one and one-to-other subsystems among them.The amount of entanglement depends on the final state kinematical configuration and maximally entangled subsystems appear in certain regions of the phase-space.Concerning the Bell non-locality, both Mermin and Svetlichny operators for 3-qubit systems were computed.We detect predictions incompatible with local realism in the whole phase-space, except for a few particular configurations, suggesting that H → γl l could serve as an ideal laboratory for testing Bell inequality.The development of possible experimental implementations is out of the scope of this work and it is deferred to future study.
Furthermore, we analyzed post-decay entanglement and autodistillation phenomena at dilepton invariant mass close to the Z-pole mass.We found that the qubit-qutrit system of gauge bosons in the H → γZ decay is maximally entangled, with minimal entanglement loss in the 2-qubit subsystems (photon-lepton and photon-antilepton) of the H → γZ → γl l decay.Interestingly, the 'MaxEnt Principle' applied to these subsystems favors s 2 w = 0.25, notoriously close to the SM value.On the other hand, when introducing CP-violating interactions in the lepton Yukawa sector, we observe that CP-effects on these entanglement measures are suppressed by lepton masses, thus these observables are not suitable for studying such kind of new physics.
Future avenues for exploration include extending this analysis to different three-body decays such as the well-studied π 0 → γe − e + or hadrons decays into three fermions.Furthermore, continuing with Higgs boson decays, a natural multipartite extension is to consider the four-fermion channel, constituting a 4-qubit system.
The two independent kinematical variables are the dilepton invariant mass, m l l = √ s, and the angle between photon and lepton, θ γl .The momentum of each particle is and β l = 1 − 4m 2 l /s is the lepton velocity in this frame.The two transverse polarization vectors of the photon, with the usual normalization, are Of course, we can chose the Higgs rest frame by performing the Lorentz transformation In that frame, the photon momentum is and a lower cut E γ cut to the photon energy is imposed in order to avoid IR divergences in the tree level contribution.Then the upper bound on the dilepton invariant mass in Eq. (2.6) is obtained.

B Helicity amplitudes
This Appendix collects the 8 helicity amplitudes and presents the numerical estimations for the decay width.The amplitudes, which were implemented in the numerical results of Section 5, include lepton mass effects and are written in terms of the generic 1-loop form factors.Using the kinematics of Appendix A, the amplitudes of Eqs.(2.2)-(2.3)as function of s and θ γl are: For each helicity amplitude, the first and second lines correspond to tree level and 1-loop contributions, respectively.The form factors a 1,2 and b 1,2 do not depend on m l , then the complete dependence on lepton mass is explicitly shown in this equation.Of course, the SM is recovered for κ l CP = 1 and δ l CP = 0. Some comments about interesting limits are in order.First, when photon has vanishing energy (s = m 2 H ), the 1-loop contribution vanishes for all helicity amplitudes and the tree level yields to IR divergences for {+ + +, + − −, − + ++, − − −} amplitudes (which are avoided by the lower cut E γ cut ).Secondly, when photon is collinear with lepton (θ γl = 0) or with antilepton (θ γl = π), just two helicity amplitudes are non-vanishing for each case: Thirdly, the massless lepton case (m l = 0, β l = 1) results in the state of Eq. (3.5).In particular, it also has IR divergences when photon is collinear with lepton and antilepton, as can be seen from Eq. (B.2).
Furthermore, we focus on the dilepton invariant mass regime of Eq. (2.6) and no additional cuts over the final particles are applied in this work.The resulting decay width for the three flavors, using both full and hybrid computations, are presented in Table 1 (see also Fig. 2).Two relevant energy ranges are considered in this table: a low-mass dilepton subsystem ∈ [0.1m H , 30 GeV], as in [54], and the complete range ∈ [0.1m H , √ s cut ].These results are compatible with previous computations [48][49][50].
The estimation of the expected number of events for LHC Run 2 + Run 3 data, using the full computation, is also included in the last column of Table 1.In this estimation, the NNNLO Higgs boson production cross section at 13 TeV is 48.61 pb, the total decay width is 4.07 MeV and the luminosity of the combined Runs 2 and 3 is 350 fb −1 [81,82].In addition, it is assumed an identification efficiency of 0.7 for each lepton in final state and consider the inclusive decay channel for the tau-lepton.These results scale trivially with the lepton identification efficiency.For the HL-LHC estimation, the expected cross section at 14 TeV is 54.67 pb and the luminosity is 3000 fb −1 , then we can expect 10 times more events.Furthermore, notice that this three-body decay channel is quite clean and the background originates predominantly from non-resonant l lγ production [53,54]

C Main results using Dalitz plots representation
Introducing the Mandelstam variables s = (p − + p + ) 2 , t = (k + p − ) 2 and u = (k + p + ) 2 , we can plot the main results of Section 5 using Dalitz plots representation, in order to be independent of the reference frame.These variables satisfy the relation s + t + u = m 2 H + 2m 2 l and the angle of the photon and lepton is written as Hence Fig. 7 presents the results of the full computation in the [ √ u, √ t] plane, already discussed in Figs.3-5.This format may be more convenient from the experimental point of view in these three-body decays, however we chose the [m l l, cos(θ γl )] plane since is more intuitive for identifying kinematical configurations.Of course the physics behind these phenomena is the same in both cases.

Figure 1 .
Figure 1.Generic diagrams corresponding to H → γl l decays at O(ℏ) in perturbation theory.The gray blobs denote the one-particle-irreducible (1PI) Green functions renormalized in the on-shell scheme.

Figure 2 .
Figure 2. Differential decay width respect to the dilepton invariant mass for l = τ, µ, e in the three kind of computations (solid lines) of Eq. (2.4).Dashed lines correspond to the tree and complete 1-loop computations in Eqs.(2.2)-(2.3).SM parameters values κ l CP = 1 and δ l CP = 0 are assumed.Shaded region correspond to 2m l ≤ m l l ≤ 0.1mH.

Figure 3 .
Figure 3. Entanglement and non-locality quantifiers for H → γτ τ in the [mτ τ , cos(θγτ )] plane, corresponding to the full computation (first row) and two-body intermediate Higgs decay (second row).Dashed lines are contours corresponding to 0.1 and 0.9 for concurrences and to 2 and 3.5 for Mermin operator.

Figure 4 .
Figure 4. Entanglement and non-locality quantifiers for H → γµμ in the [mµμ, cos(θγµ)] plane, corresponding to the full (first row) and hybrid (second row) computations.Dashed lines are contours corresponding to 0.1 and 0.9 for concurrences and to 2 and 3.5 for Mermin operator.

Figure 5 .
Figure 5. Entanglement and non-locality quantifiers for H → γeē in the [meē, cos(θγe)] plane corresponding to the full computation.Dashed lines are contours corresponding to 0.1 and 0.9 for concurrences and to 2 and 3.5 for Mermin operator.

Figure 7 .
Figure 7. Entanglement and non-locality quantifiers in the [ √ u, √ t] plane corresponding to the full computation for tau-lepton (first row), muon (second row) and electron (third row).

Table 1 .
. lepton case √ s min [GeV] √ s max [GeV] Γ full [KeV] Γ hybrid [KeV] N Run 2+3Decay width Γ H→γl l for both full and hybrid computations, in the range of invariant masses ∈ [ √ s min , √ s max ].The estimation of the expected number of events for LHC Run 2 + Run 3 data is also included in the last column (see text for details).