Searching Quantum Entanglement in pp → ZZ process

Quantum entanglement [1], a cornerstone of quantum mechanics, represents a fascinating phenomenon where two or more particles become intricately correlated, such that the change in one particle's quantum state instantaneously influences the state of another termed as "spooky action at a distance" by Einstein et. al in their paper [2]. This phenomenon might be the most outstanding characteristics of quantum mechanics that truly separates quantum mechanics from classically deterministic theory. This could be proven by testing the violation of Bell inequalities [3] using the theory of quantum mechanics, which is unfeasible within any theory that advocates local realism.

Studies about quantum entanglement of multi-particle system at the highest energy frontiers of particle physics has recently gained enormous interest among both theorists and experimentalists. There have been breakthroughs in the study of quantum entanglement and quantum information science at the high energy frontier [4−11]. For example, several works show that violation of Bell Inequality could be measured in heavy quark system, i.e. $ t\bar{t} $, $ b\bar{b} $ at the Large Hadron Collider (LHC) [12−25]. Considering the heavy mass and very short life time, the $ t\bar{t} $ system can be an ideal platform to perform quantum tests at the LHC. Both ATLAS and CMS collaborations [26−28] have measured entanglement among top quark pairs with high sensitivity. Up to now, search for QE and test Bell nonlocality between other pair of system have continuously gained growing interests at collider experiment.

Quantum state tomography [29−31], determining the density matrix from an ensemble of measurements, has emerged as a cornerstone technique for reconstructing the quantum state of a system, providing invaluable insights into quantum correlations, coherence, and entanglement [32−35]. While initially developed within the realm of quantum optics and low-dimensional systems, advancements in experimental and theoretical physics have extended its applicability to high-energy particle physics, particularly in exploring the quantum properties of massive gauge particles produced at particle colliders [36−52]. The advent of high-energy colliders such as LHC and the proposed muon collider opens a promising avenue for probing fundamental aspects of quantum mechanics in previously inaccessible regimes.

The study of massive spin-1 particles, such as the electroweak gauge bosons WW and ZZ, is particularly compelling in this context. These particles play a critical role in the Standard Model of particle physics, mediating weak interactions and participating in processes that probe electroweak symmetry breaking [53−56]. At a muon collider [57, 58], where precise beam properties and high luminosity enable clean experimental conditions, the production of such particles offers an unparalleled opportunity to investigate their quantum properties. Specifically, quantum state tomography of WW and ZZ boson pairs provides a direct method to study their spin correlations, polarization states, and entanglement, enriching our understanding of the quantum nature of the Standard Model.

In this work, we explore quantum properties of ZZ system such as quantum correlation, concurrence and checked separability of this state. This is realized through quantum state tomography of two massive spin-1 particles produced in a 13 TeV proton-proton collision, whose spin density matrix is determined from a simulated events of $ pp\to ZZ\to4\ell $. Quantum observables like concurrence are measured to examine whether the ZZ pairs are entangled or separable at such an extremely relativistic environment.

Spin density matrix is mathematical structure that describes the quantum state of a given system, including both pure and mixed states and encapsulating all the observable related to quantum information such as concurrence [59], quantum discord [60], quantum magic [61] and entanglement entropy [62]. In quantum mechanics, state vector $ \left|{\Psi}\right \rangle $ is dominantly used to describe the quantum mechanical state of the corresponding system. For pure state, the spin density matrix, $ \rho = \left|{\Psi}\right \rangle\left\langle{\Psi}\right| $ has the equivalent role as $ \left|{\Psi}\right \rangle $, while for a mixed state, the spin density matrix is defined as the convex sum

where $ p_i $ is the classical probability of the i-th pure state $ \left|{\psi_i}\right \rangle $ and satisfies $ \sum p_i=1 $. Density matrix is a positive semi-definite operator, which means $ \left\langle{i}\right|\rho\left|{i}\right \rangle $, where $ \left|{i}\right \rangle $ is the base state in the complex Hilbert space $ {\cal{H}} $. The total sum probabilities $ p_i $ in Eq. 1 results from the requirement that $ \rho $ must have unit trace.

To reconstruct the spin density matrix from experimental or a simulated pseudo data a parametrization of it must be applied. As best know in the spin context of quantum mechanics, the Pauli spin matrices $ \sigma_i $ is usually chosen to parametrize the spin density matrix of s = 1/2 system, which constitutes qubit. For example, the following form of spin density matrix for a two-qubit system is parametrized with $ \sigma_i $:

where $ A_i $ and $ B_i $ are components of the polarization vector of each qubit, and $ C_{ij} $ is the entry of the correlation matrix. For any a two-qutrit system, composed of two spin-1 particles, the general spin density matrix can be represented using traceless Gell-Mann matrices [63]:

where $ T^i $ are the $ 3\times 3 $ Gell-Mann matrices with $ i=1,2,\dots,8 $, and $ I_3 $ is the 3-dimensional identity matrix. The coefficients $ A_i $, $ B_i $ are spin polarization parameters or the components of the polarization for each qutrit, and $ C_{ij} $ are the spin correlation parameters. These parameters can be obtained by projecting the spin density matrix on the desired subspace basis as:

Analytically, the obtained coefficients in Eq. 4 are Lorentz invariant which depend only on the energy E, momentum or velocity and the scattering angle in the CM frame. Once these coefficients are known, it is quite straightforward to determine the observables quantifying the entanglement in the massive gauge boson system.

Quantifying entanglement of the state of a quantum system is generally challenging as the complexity of the problem increases with the system dimensionality. For a pure state described by a single vector in the Hilbert space, or equivalently by a density matrix, this can be solved by their Schmidt decomposition [1]. As for the mixed states, a general approach is to compute the so-called entanglement witness, a quantities that give sufficient conditions to establish the presence of entanglement in the system. Such an observable quantity is concurrence, which already mentioned in the above section. This is reliable entanglement measure for bipartite system, consisting of two particles [64, 65].

Usually for a pure state described by a joint-state vector, the corresponding concurrence is defined as

where $ \rho_r $ is the reduced density matrix of the subsystem obtained by tracing over the degree of freedom of either subsystem. However, concurrence of any mixed state described Eq. 1 is given by means of the concurrence of the pure states as

where infimum is taken over all the possible decompositions of $ \rho $ into pure states. A vanishing value of Eq. 5 implies that the pure state of the bipartite system is separable. However, due to the complexity in the evolution of mixed states, evaluating its concurrence, Eq. 6, to a specific value is quite challenging. Therefore, finding a lower bound of concurrence for mixed states, rather than obtaining the exact value of, also unequivocally indicates the presence of entanglement. This lower bound is analytically computable [66, 67]:

or equivalently, the above lower bound can be written for two-qutrit system of spin-1 massive bosons

By using the explicit form of the spin density matrix parametrized as Eq. 2, the lower bound of concurrence $ \mathscr{C}_2 $ can be written in terms of the coefficients given in Eq. 4:

The state is entangled if the above quantity is positive, while that lower bound is negative or equal to zero the concurrence test is inconclusive.

Due to its short lifetime, Z boson can not be detected directly inside the detector. The most abundant decay channel is hadronic decay, which accounts for up to 70% of the total decay. The leptonic decay channel is convenient to manipulate both in the simulation and in the actual experiments. The information about entanglement of the states are passed to the decay products namely, charged leptons. The angular direction of each of these leptons is correlated with the direction of the spin of their parent Z boson in such a way that the normalized differential cross-section of the process $ ZZ\to \ell^+_1\ell^-_1 \ell^+_2\ell^-_2 $ may be written as [68]

in which the solid angle measure $ {\rm{d}}\Omega^{\pm} = \sin\theta^{\pm}{\rm{d}}\theta^{\pm}{\rm{d}}\phi^{\pm} $ are written in terms of the spherical coordinates for the momenta of the final charged leptons in the respective rest frame of the decaying particles. The density matrix $ \rho_{V_1V_2} $ encodes the production information of the joint ZZ system, while the decay density matrix $ \Gamma_i $ passes the relevant entanglement information about the ZZ state to the final decay products, so that we can recover the spin density matrix from these angular distribution. The explicit form of the decay density matrix $ \Gamma_i $ can be obtained using polarized decay amplitudes [68]

where $ \theta $ and $ \phi $ are the polar angles of the momentum of the negative charged lepton in the rest frame of mother particle, which is Z boson. The factor $ \eta_{\ell} $ is function of electroweak mixing angle [69] and is equal to 0.15 when using the latest measurement of mixing angle [70].

The functions given in Eq. A2 and together with the matrix $ A_i^j $ A3 in Appendix can be used to extract the polarization and correlation coefficients of the density matrix given in Eq. 3:

where $ \tilde{p} $ is the corresponding Wigner P symbols for Z boson decay, and given as

The unit vectors $ \hat{\bf{n}}(\theta, \phi) $ in Eq. 12 is defined to point the momentum of the decay leptons in the rest frame. Because the two Z bosons are indistinguishable, a symmetry factor should be implied by exchanging the labels $ i $ and $ j $, so that the above coefficients are determined as:

Eq. 14 provide the means to reconstruct the correlation functions fo the density matrix from the distribution of the lepton momenta and thus allow to derive the expectation value of the lower bound of concurrence $ \mathscr{C}_2 $ from the data.

A sensitivity test for the entanglement content is done by conducting pseudo-experiments of the $ pp\to ZZ $ process at the LHC with a center-of-mass energy of $ \sqrt{s}=13 $ TeV. Unlike massive gauge boson pair production in the Higgs decay, in which one of the gauge boson is produced off-shell, here both ZZ bosons are produced on-shell. The Feynman diagram of the $ pp\to ZZ $ process, in which $ Z\to\ell^+\ell^- $ $ (\ell = e, \mu ) $, where $ \ell $ can be either electron or muon, is shown in Fig. 1. This is the signal process in our study. Simulation is achieved through publicly available Monte-Carlo software $ \mathtt{MadGraph5\_aMC@NLO} $ [71]. During the simulation the relevant spin-correlation and Breit-Wigner effects are include through $ \mathtt{MadSpin} $ program [72, 73] that comes with $ \mathtt{MadGraph5\_aMC@NLO} $.

Figure 1.  The Feynman diagram of the $ pp\to ZZ $ process.

The number of generated events and corresponding cross section sections are given in the Table 1. The relevant cross sections are computed using $ \mathtt{MadGraph5\_aMC@NLO} $ at the Leading Order (LO) level. The center-of-mass energy (CM) is set to $ \sqrt{s} = 13 $ TeV for three different collider luminosities corresponding to the LHC Run 2+3 and HL-LHC luminosities, respectively. When calculating the cross section for both the signal and background process, we use the default Parton Distribution Functions (PDFs) set $ \mathtt{NNPDF} $23 [74] and set the factorization and normalization scale to the physical mass of the Z boson.

Process	Cross section (fb)	Events(Lumi=137 fb−1)	Events(Lumi=300 fb−1)	Events(Lumi=3 ab−1)
$ p\ p\to ZZ $(signal)	$ 42.29 $	5793	12687	126870
ZZZ	24.16	3309	7248	72480
WWZ	29.76	4077	8928	89280

Table 1.  The signal and background processes considered in this work in the first column, total cross section in the second column and the event numbers given for a total luminosity.

As for the possible background process, we consider the following three process as background to our signal:

● $ pp\to WWZ $

● $ pp\to WZZ $

● $ pp\to ZZZ $

Fig. 2 depicts the event distributions for both the signal and total backgrounds as function of invariant mass of final four leptons. As is shown, the background events does not matter in the signal region. Based on this, we can safely ignore any contamination coming from background events. The $ p p \to Z Z $ signal process exhibits strong background suppression due to the reconstruction of the ZZ pair via isolated leptons, resulting in minimal background contamination. The dominant background, WWZ, involves two $ W $ bosons decaying to charged lepton and corresponding neutrino, leading to a broad four-lepton invariant mass ($ M_{4\ell} $) distribution. This distinct kinematic feature allows clear separation from the narrow resonance of the ZZ signal. The final lepton coming from Z boson pairs can be identified as two electrons and two muons, four electrons or four muons. For two electrons and two muons, it is easy to distinguish their parent Z boson. In the case where the ZZ system decays into four same-flavor leptons, we distinguish the signal by combinations: we iterate through all possible lepton pairings and select the combination in which the invariant mass of the reconstructed Z boson is closest to the on-shell mass. This strategy is effective because both Z bosons in the signal process are on-shell. In steps of the numerical calculation, we extract the polarization and correlation coefficients of the spin density matrix from each single event using Wigner P symbols. Running this procedure over all events gives an average value and the standard deviation of our observables. We produced one million events in total and divide them into 1000 pseudo experiments, each one includes events that matches to the total reachable events in an actual experiments. The statistical uncertainty is determined with these pseudo experiments by repeating over all of the pseudo experiment. The parameters of the spin density matrix, $ A_i $, $ B_i $ and $ C_{ij} $ are shown in a matrix plot depicted in Fig. 3. It is clear that there are large portion of non-zero coefficients indicating the correlations between two Z boson. The concurrence observable that quantifies entanglement is given as the function of $ \cos{\theta} $ which is depicted in Fig. 4. For the Run 2+3 data, the gray band, the computed lower bound of the concurrence is relatively larger than one, which is unphysical, while it is close to unity when $ \theta\to 0 $ or $ \pi $. However, the HL-LHC data, the colored one, gives rather physical results when the scattering aligns with beam direction. Because of the high event numbers, the uncertainty of the concurrence for this run is also quite small. Quantum entanglement test is of highly statistical meaning, which means the results computed for the observables at hand, both the mean and deviation, are dependent on the event numbers at the actual experiment. To confirm this statistical nature of the entanglement, a Gaussian distribution of the lower bound of the concurrence, $ C_{\rm{MB}}^2 $, is obtained corresponding to LHC Run 2+3 and HL-LHC luminosities. The results are shown in the Fig. 5. It can be seen clearly that, $ C_{\rm{MB}}^2 $ distribution converges to the mean value depicted by the red-dashed line for the HL-LHC data, while it diverges for the Run 2+3 data.

Figure 2.  (color online) The event distribution of the signal and background process.

Figure 3.  (color online) The correlation matrix element extracted from the simulated data of ZZ production. The leftmost column contains the $ A_i $ parameters of Z boson and the bottom row contains the parameters $ B_i $ and the center middle table plot contains $ C_{ij} $ coefficients

Figure 4.  (color online) Concurrence distribution as function of scattering angle cosine $ \cos{\theta} $. The thick central lines in are the central value and colored bands represent statistical deviations obtained from events corresponding to LHC Run 2+3 and HL-LHC luminosities.

Figure 5.  (color online) $ C_{\rm{MB}}^2 $ distribution computed from LHC Run 2+3 and HL-LHC luminosity data. Red dashed line is the mean value obtained from HL-LHC luminosity.

Based on the distributed results shown in Fig. 5, The expectation value of $ C_{\rm{MB}}^2 $ observable converges to some point when the luminosity of the collider increases. As a test, we compute a luminosity spectrum for $ C_{\rm{MB}}^2 $ and the results are shown on the left panel of the Fig. 6: the points are the expectation value and the error bar is the uncertainty. Obviously, $ C_{\rm{MB}}^2 $ converges approximately to 0.375 as the luminosity increases. Meanwhile, the corresponding significance can reach to 3.75 $ \sigma $ for the HL-LHC luminosity, which is shown on the right panel of the Fig. 6. This is well-above the 3$ \sigma $ confidence level.

Figure 6.  (color online) The expectation value (left panel) of $ C_{\rm{MB}}^2 $ as function of collider luminosity, and the significance distribution (right panel).

Quantum entanglement is one of the most distinctive features of the quantum mechanics. It is widely being examined and tested by particle physics community especially the high energy physics frontiers such as at the LHC and future colliders. In this work, we apply quantum state tomography to ZZ system to inspect the presence of the quantum entanglement in the $ pp\to ZZ\to 4\ell $ process. The theoretical framework in which how the spin density matrix is computed or parametrized for the ZZ system is given as well as a through explanation of the quantum state tomography is also provided.

Importantly, a fast Monte-Carlo simulation of the $ pp\to ZZ $ is done using $ \mathtt{MadGraph5\_aMC@NLO} $ software so that a large number of pseudo experiments can be used to determine whether quantum entanglement exist or not in the ZZ state. To achieve this, firstly, Gell-Mann spin matrices are used to parameterize the spin density matrix. Secondly, one million events corresponding to LHC Run 2+3 and HL-LHC luminosity are produced for $ pp\to ZZ $ process for the LHC Run 2+3 and HL-LHC luminosity, based on which the total events are sliced into 1000 pseudo experiments. Then using the theoretical framework, all the coefficients of the spin density matrix is determined by the angular distribution of the final lepton pairs coming form Z boson decay. Entanglement quantifier, the lower bound of the concurrence, is measured using these pseudo experiments. The significance of this quantity can reach to 2$ \sigma $ with combined LHC Run 2+3 data and 3.75 $ \sigma $ for HL-LHC data, making it possible to measure quantum entanglement between two massive gauge boson pairs in an actual LHC experiment.

The functions $ \mathfrak{q}^{n}_\pm $ and $ \mathfrak{p} $ are defined as follows according to [75, 76]:

where $ \theta^\pm $ and $ \phi^\pm $ are the polar and azimuthal angles of the decay lepton $ l^\pm $ in the Z boson rest frame, respectively. The function $ \mathfrak{p} $ is defined as:

The matrix $ \mathfrak{a}_m^n $is defined as:

The $ 3\times3 $ Gell-Mann matrices are defined as: