A quantum strategy to compute the jet quenching parameter $\hat{q}$

The idea of simulating the dynamics of complex quantum physical systems by using other simpler and controllable quantum systems (which we shall refer to as quantum computers) was first realized by Feynman in the 80’s Feynman:1981tf . Since then quantum simulation as seen widespread application in physics, chemistry and many other areas Georgescu:2013oza ; Schuld_2014 ; Or_s_2019 .

In recent years, a big effort has been made towards exploring to what extent quantum computers might enhance our understanding of High Energy/Nuclear Physics (HEP/NP) phenomena carlson2018quantum ; cloet2019opportunities ; matchev2020quantum . Some of the most recent studies in the application of quantum computing to HEP/NP phenomenology have resulted in the formulation of quantum parton showers Bepari:2020xqi ; Bauer:2019qxa , quantum jet clustering algorithms Pires:2021fka ; Wei:2019rqy , digital simulation of effective field theories Bauer:2021gup and of the propagation of hard probes in a thermal QCD bath deJong:2020tvx .

In this paper, we propose a strategy to use a quantum digital computer to simulate the evolution of a single parton in the presence a QCD background field. In particular, we are interested in the $\alpha_{s}^{0}$ effect, corresponding to the broadening of the parton’s momentum. Although this effect has been extensively studied in jet quenching theory Gyulassy:2002yv ; Blaizot:2012fh , it is only easily computed for isotropic and homogeneous media, where the field fluctuations behave as white noise. More interestingly, at the amplitude level, the associated in-medium propagators are the building blocks of jet quenching formulation for e.g. medium-induced gluon radiation. For this reason, we argue that our algorithmic implementation can be considered as a first step towards a complete simulation of the in-medium parton cascade with quantum color coherence.

We consider an energetic parton that emerges from a hard scattering event and then propagates inside a QCD medium. The net effect of the medium is to alter the initial transverse momentum of the parton. The underlying gauge field is treated stochastically, in line with the usual approach employed in jet quenching theory and phenomenology. As a consequence, our algorithm consists in a hybrid classical-quantum strategy mueller2020deeply ; PhysRevX.6.031045 ; deJong:2020tvx , with the parton evolution in time being tracked, at the amplitude level, by the quantum computer and the gauge fields being provided as an input to the circuit. We will not make an attempt to improve here on these dual description as an eventual future implementation of the quantum computation of the gauge fields would be straightforward to implement in our procedure.

Although for an actual implementation, crucial aspects such as quantum error correction Terhal_2015 , encoding details or Trotter error analysis Childs_2021 have to be considered, we will mostly stay at a more conceptual level and leave such an analysis for future work. In addition, we will try to highlight the connection between our approach and standard jet quenching treatment of momentum broadening, thus making what follows more relevant for an interested reader.

The present manuscript is divided as follows: section 2 briefly reviews the physics of a hard parton propagating in a gauge classical background field, while in section 3 we provide the equivalent Hamiltonian formulation. In the remainder of the section we detail how such a problem can be implemented in a digital quantum computer. Finally, in section 4 we detail how to deal with non-trivial evolution in color space and finally section 6 gives the paper’s main conclusions and outlines possible future research paths. Details on most sections are provided in four appendices.

In this section we review the theoretical description of a highly boosted parton propagating in a classical background gluon field ${\cal A}$. We will consider the cases where the propagating parton, a quark, is in the singlet or fundamental color representations, though we will mostly assume that the evolution in color space is trivial. Detailed and ample discussions on jet quenching theory can be found, for example, in Blaizot:2012fh ; CasalderreySolana:2007zz ; Blaizot:2015lma ; Mehtar-Tani:2013pia and references therein.

We begin by considering a highly energetic quark, being produced inside a QCD medium (whose exact origin is not relevant) from a hard process. Due to the highly boosted kinematics, the quark’s four-momentum $p=(p^{0},{\boldsymbol{p}},p^{z})$ is more conveniently expressed in light-cone coordinates $p\equiv(\omega,{\boldsymbol{p}},p^{-})=((p^{0}+p^{z})/2,{\boldsymbol{p}},p^{0}-p^{z})$, with ${\boldsymbol{p}}$ the transverse momentum, $\omega$ the light-cone energy and $p^{-}$ the minus component of the quark’s momentum. The quark is assumed to be moving in the plus direction, so that $\omega$ is the large momentum component.

It is also convenient to work in the light-cone gauge for the background field ${\cal A}^{\mu}$, taking ${\cal A}^{+}=0$ and fixing the residual gauge freedom such that ${\cal A}^{-}$ is the only non-vanishing component of the background field Blaizot:2008yb . Additionally, because $p^{+}$ is the large component of the quark’s momentum, its interaction with the gluon field is highly localized in $x^{-}$, so that we can simplify the spacetime dependence of the field to be ${\cal A}^{-}(x^{+},{\boldsymbol{x}},x^{-})\equiv{\cal A}(x^{+},{\boldsymbol{x}},0)$, dropping the $x^{-}$ dependence in what follows.

Finally, in the boosted regime the local quark-field spin-flip interactions are energy suppressed and can be ignored. Thus, for each field insertion ${\cal A}^{-}$ and in the strict high energy limit, the large momentum component $\omega$ and transverse component ${\boldsymbol{p}}$ are conserved, with the quark state acquiring an eikonal phase. It is however usual to relax this approximation and allow for a small transverse momentum transfer at each vertex, while light-cone energy is still conserved. Accounting for this sub-eikonal corrections, the quark propagation in the QCD field can be reduced to the study of a two dimensional non-relativistic quantum system Blaizot:2015lma .

To make this discussion more quantitative, we consider the in-medium scalar quark propagator $G(t,{\boldsymbol{x}};0,{\boldsymbol{y}})$ in the transverse plane, between spacetime points $(0,{\boldsymbol{y}})$ and $(t,{\boldsymbol{x}})$ Blaizot:2012fh . This propagator is the Green’s function to the following two dimensional Schrodinger equation

where we have contracted the background gauge field with the respective $SU(3)$ generators $T$ in the adequate representation. The remaining terms are diagonal in color space. This equation explicitly shows that the quark propagation is equivalent to a non-relativistic two dimensional quantum mechanical system, describing a single particle evolving in time with the Hamiltonian Blaizot:2015lma

where $\omega$ plays the role of a mass and light-cone time plays the role of time¹¹1In the boosted regime, light-cone time $x^{+}$ becomes the same as time $x^{0}$ since $x^{+}=(x^{z}+x^{0})/2\sim x^{0}$.. In the strict eikonal limit, where $\frac{{\boldsymbol{p}}^{2}}{\omega}\to 0$, the kinetic term drops out and the evolution leads to the state acquiring a field dependent phase, as mentioned above.

From Eq. 2 one can construct the time evolution operator (with ${\cal T}$ the time ordering operator)

which acts on the infinite dimensional Hilbert space of a single free particle in two spatial dimensions, such that from an initial state $\ket{\psi_{0}}$ at time $t=0$ one obtains the time evolved state $|\psi_{t}\rangle$ via

The Hilbert space is spanned by the position eigenvectors $\ket{{\boldsymbol{x}}}$ or by their Fourier pair $\ket{{\boldsymbol{p}}}$. These two bases are convenient since $\hat{{\boldsymbol{p}}}\ket{{\boldsymbol{p}}}={\boldsymbol{p}}\ket{{\boldsymbol{p}}}$ and $\hat{{\cal A}}^{-a}(t,\hat{{\boldsymbol{x}}})\ket{{\boldsymbol{x}}}={\cal A}^{-a}(t,{\boldsymbol{x}})\ket{{\boldsymbol{x}}}$, where we used the hats to highlight the difference between operators and c-numbers; we also used the fact that the quark-medium interaction is localized in position space (and conversely delocalized in momentum space).

We now detail how to frame single particle momentum broadening in terms of a digital quantum simulation algorithm, implementing Eq. 4. The algorithm, summarized in Fig. 1, can be divided as follows:

1. 1.

   Input – i) Template distribution to be loaded as an initial state $\ket{\psi_{0}}$ ii) A list of $m$ field configurations ${\cal A}^{-}$ with the associated weights $p_{{\cal A}^{-}}$, storing the probability of generating each configuration;

2. 2.

   Encoding – Map between the degrees of freedom of the quantum system and the qubits;

3. 3.

   Initial state preparation – Preparation of $\ket{\psi_{0}}$;

4. 4.

   Time evolution – Implementation of Eq. 4;

5. 5.

   Measurement – Retrieving physical information by measuring the qubits, according to a sensible protocol;

6. 6.

   Output – For each field configuration the algorithm will output the expected value of a random variable $\chi$, which should be then medium-averaged over all $m$ configurations.

In this section we assume that the initial quark probe is in the fundamental $SU(3)$ representation. As a consequence, the $H_{A}$ component of the Hamiltonian now has a non-trivial color structure, i.e. $A\cdot T=A^{a}T^{a}=\frac{1}{2}A^{a}\lambda^{a}$, where $\lambda^{a}$ denotes the eight Gell-Mann matrices. To deal with this modification, we further split the time evolution operator to take the form $U=U_{K}\cdot U_{A^{1}}\cdot U_{A^{2}}\cdots U_{A^{8}}$. Additionally, we must track the color of the quark as it evolves. To do that, we add a new register with two qubits, which stores the color state of the quark. In particular we use the following map between the logical and physical states: $\ket{0,0}\equiv\ket{\rm red}=\ket{R}$, $\ket{0,1}\equiv\ket{\rm green}=\ket{G}$, $\ket{1,0}\equiv\ket{\rm blue}=\ket{B}$ and $\ket{1,1}\equiv\ket{W}$, with the latter state not being physical and therefore absent from any calculation.

We now detail how to implement $H_{A^{1}}$, with the other values of $a$ following analogous implementations. The first Gell-Mann matrix is given by

where in the second step we have embedded $\lambda^{1}$ into the two qubit Hilbert space. The action of $\tilde{\lambda}^{1}$ is thus to color rotate the quark state between the $\ket{R}$ and $\ket{G}$ states, which can lead to a non-trivial evolution in color space. One can diagonalize the above matrix using a control Hadamard gate $CH$

such that we can write $H_{A^{1}}$, in $k_{t}^{\rm th}$ time interval, in terms of a diagonal operator (here we drop all spacetime dependence for readability)

Here we made use of the extended Pauli operator $\tilde{\sigma}^{Z}={\rm diag}(1,-1,0,0)$⁶⁶6To be more precise, this definition takes $\tilde{\sigma}^{Z}$ to be non-unitary, unlike $\sigma^{Z}$. This is done, in order to ensure that only the $\ket{R}$ and $\ket{G}$ states transform non-trivially.. Finally, to compute the exponential of the tensor product we notice that

where $\ket{c}$ denotes the two qubits register storing the state of the quark in color space. From the previous equation it is easy to observe that only $\ket{0,0}$ and $\ket{0,1}$ states result in a phase, the former with a $-i$ prefactor and the latter with a $+i$. Notice however, that due to the application of the diagonalizing gate $1\otimes CH$, the evolution in the physical RGBW basis is off-diagonal. The implementation of Eq. 20 is given in Fig. 4.

Clearly this strategy is only possible as long as the quark is in a small color representation – in the previous example, the color degrees of freedom were treated by adding only two extra qubits and doubling the number of time evolution operators $U_{A^{a}}$, for each $a$.

Another important consequence of including non-trivial color evolution is the fact that the final and initial state are differential in color. Therefore, when preparing the state one has to set colors either according to some initial state prescription or in an equitative way. Consequently, in the measurement protocol the output must be color averaged, which can be performed classically⁷⁷7This is not necessary if the qubits storing the color information are not measured..

In the previous sections we gave a conceptual outline on how to quantum compute the single particle momentum broadening distribution and from it extract meaningful physical information. Here we give a rough estimate on the typical values for the circuit parameters based on estimates for the typical physical scales obtained from jet quenching/saturation physics phenomenology.

Let us first estimate the lattice spacing $a_{s}$ and the number of qubits necessary per dimension $n_{Q}$. For that we recall that, when traversing a dense medium of length $L$, the quark will acquire an average transverse momentum of the order of the saturation scale, $\langle{\boldsymbol{p}}^{2}\rangle\sim\hat{q}L\equiv Q_{s}^{2}$. We are interested in length scales $L$ of the order of the nuclear radius of heavy elements, like Pb or Au, which we take to be $L\sim\mathcal{O}(10\,\rm fm)=\mathcal{O}(50\,\rm GeV^{-1})$. The value of the jet quenching parameters $\hat{q}$ varies drastically between different experimental set-ups, due to the different energy scales being explored. To bridge RHIC, LHC and EIC experimental conditions, we assume that $\mathcal{O}(0.1\,\rm GeV^{2}fm^{-1})\leq\hat{q}\leq\mathcal{O}(10\,\rm GeV^{2}fm^{-1})$ Feal:2019xfl ; Liu:2018trl ; Barata:2020jtq . The saturation scale $Q_{s}^{2}$ is then approximately bounded by $Q_{s}^{2}\sim\mathcal{O}(1-100\,{\rm GeV}^{2})$.

Setting the ultraviolet momentum cutoff induced by the digitization ${\boldsymbol{p}}_{\rm max.}$ to be much larger than the saturation scale $Q_{s}$, we obtain

thus

Conversely, we require that the momentum space discretization is neither too coarse nor too fine. A simple way to ensure this is to impose

with $\mu$ an infrared (model dependent) regulator, related to the medium Debye mass; typically $\mu~{}\sim\mathcal{O}(0.1-1\,\rm GeV)$ Barata:2020jtq ; GW ; HTL and we used the previous estimates to reduce the problem to the ratio between the soft and hard scales at play. Recalling that $N_{s}=2^{n_{Q}}$, we obtain

Thus, one roughly needs $\mathcal{O}(2^{7}=128)$ states per dimension to adequately discretize the problem. In practice this number will have to be larger since the correct energy ratio should be $\mu/|{\boldsymbol{p}}_{\rm max.}|$, which here we took $|{\boldsymbol{p}}_{\rm max.}|=Q_{s}$. This is a rather rough lower bound, and larger values should be considered such that the peak of the broadening distribution is well captured. Even so, one would expect that (roughly) $n_{Q}<20$, which means that $N_{s}<\mathcal{O}(10^{6})$. This allows us to argue that the classical operations detailed in the previous sections needed to implement $U_{A}$ can be performed in a classical computer.

Let us now consider the longitudinal scales entering the problem and estimate the number of time steps $N_{t}$. We recall that in the multiple soft scattering approximation, one usually requires that the mean free path of the quark $\lambda$ is much larger than the typical correlation length in the medium $1/\mu$. This ensures that spatially delocalized scattering centers are not color correlated. On the other hand we also have that in order for a scattering to occur $\lambda\leq L$, leading to

It is also typical to define the opacity of the medium as $\chi_{\rm med.}\equiv L/\lambda$ Gyulassy:2000er ; Wiedemann:2000za , corresponding to the expected number of in-medium scatterings. Therefore, it is natural to identify $\chi_{\rm med.}\sim N_{t}=L^{\prime}/\varepsilon_{t}$. We can then write

The remaining circuit parameter that directly depends on the physics one wants to explore is $m$, the number of field configurations to be generated. As alluded above, the numerical value for $m$ intrinsically depends on the model/prescription for the gauge field and its fluctuations, and therefore it is tied to the its physical origin. As such, and since in our treatment we have avoided discussing details of the background, we leave an estimation of this parameter for future work where a model for ${\cal A}$ is chosen.

In this paper we have outlined a quantum simulation algorithm to extract the jet quenching parameter $\hat{q}$. Our approach is a hybrid one, with the in-medium jet evolution being quantum simulated, while the background field is treated as an external stochastic parameter, given as an input to the algorithm. The connection to standard jet quenching language is immediate, unlike more recent efforts which rely on open quantum system formulation of jet quenching deJong:2020tvx , still not fully developed (see however Vaidya:2020cyi ; Vaidya:2021mjl ).

The overall algorithm requires $2n_{Q}+l$ qubits (assuming one can re-use ancillas) and $\mathcal{O}(N_{t}\times{\rm polylog}N_{s})$ basic gate operations. However, there is an underlying classical cost coming from the $m\times N_{t}\times N_{s}^{2}$ evaluations of the gauge field. This is the major drawback of our strategy, since it is not guaranteed that the classical evaluations of ${\cal A}$ can be performed efficiently. Additionally, there is an overall additional polynomial cost in the measurement section, if one decides to scan several values of $\alpha$. For an actual implementation in a NISQ device preskill2018quantum , these constraints should not be limiting and we hope in the future more efficient algorithms can also be found. Nonetheless, we expect that our method can not outperform current classical approaches.

In future work, we plan to implement our strategy in one spatial dimension (assuming azimuthal symmetry), bench-marking to known results for $\hat{q}$ from jet quenching phenomenology. This would allow a better understanding of the merits of our quantum approach, compared to known classical methods.

Going beyond $\alpha_{s}^{0}$ effects is of course of extreme relevance and the main motivation of our work. Indeed, since broadening is a classical effect, there is little advantage in applying quantum computing techniques to study it. However, a natural but non-trivial next step would be to include parton branching into the evolution operator. This is a purely quantum effect. If one is able to quantum simulate such a process efficiently, then interference contributions, inaccessible to classical Monte Carlo codes, can be exactly taken into account. The major obstacle to overcome is the fact that particle number is no longer conserved, and thus a new formulation of the problem is needed. Nonetheless, since broadening is a key element of in-medium propagation, the present algorithm provides a first step in this direction.