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

Jet quenching, the modification of the properties of a QCD jet when the parton cascade takes place inside a medium, is an intrinsically quantum process, where color coherence effects play an essential role. Despite a very significant progress in the last years, the simulation of a full quantum medium induced cascade remains inaccessible to classical Monte Carlo parton showers. In this situation, alternative formulations are worth being tried and the fast developments in quantum computing provide a very promising direction. The goal of this paper is to introduce a strategy to quantum simulate single particle momentum broadening, the simplest building block of jet quenching. Momentum broadening is the modification of the quark or gluon transverse momentum due interactions with the underlying medium, modeled as a QCD background field. At the lowest order in $\alpha_s$ that we consider here, momentum broadening does not involve parton splittings and particle number is conserved, greatly simplifying the quantum algorithmic implementation. This quantity is, however, very relevant for the phenomenology of RHIC, LHC or the future EIC.


Introduction
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 [1]. Since then quantum simulation as seen widespread application in physics, chemistry and many other areas [2][3][4].
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 [5][6][7]. 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 [8,9], quantum jet clustering algorithms [10,11], digital simulation of effective field theories [12] and of the propagation of hard probes in a thermal QCD bath [13].
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 α 0 s effect, corresponding to the broadening of the parton's momentum. Although this effect has been extensively studied in jet quenching theory [14,15], 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 [13,16,17], 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 [18], encoding details or Trotter error analysis [19] 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 nontrivial 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.

Hard parton propagating in a background field
In this section we review the theoretical description of a highly boosted parton propagating in a classical background gluon field 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 [15,[20][21][22] 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 , p, p z ) is more conveniently expressed in light-cone coordinates p ≡ (ω, p, p − ) = ((p 0 + p z )/2, p, p 0 − p z ), with p the transverse momentum, ω 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 ω is the large momentum component.
It is also convenient to work in the light-cone gauge for the background field A µ , taking A + = 0 and fixing the residual gauge freedom such that A − is the only non-vanishing component of the background field [23]. 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 Finally, in the boosted regime the local quark-field spin-flip interactions are energy suppressed and can be ignored. Thus, for each field insertion A − and in the strict high energy limit, the large momentum component ω and transverse component 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 [21].
To make this discussion more quantitative, we consider the in-medium scalar quark propagator G(t, x; 0, y) in the transverse plane, between spacetime points (0, y) and (t, x) [15]. 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 [21] where ω plays the role of a mass and light-cone time plays the role of time 1 . In the strict eikonal limit, where p 2 ω → 0, the kinetic term drops out and the evolution leads to the state acquiring a field dependent phase, as mentioned above.

Quantum simulating momentum broadening
From Eq. 2 one can construct the time evolution operator (with 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 |ψ 0 at time t = 0 one obtains the time evolved state |ψ t via The Hilbert space is spanned by the position eigenvectors |x or by their Fourier pair |p . These two bases are convenient sincep |p = p |p andÂ −a (t,x) |x = A −a (t, x) |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. Input -i) Template distribution to be loaded as an initial state |ψ 0 ii) A list of m field configurations A − with the associated weights p A − , storing the probability of generating each configuration; 2. Encoding -Map between the degrees of freedom of the quantum system and the qubits; 3. Initial state preparation -Preparation of |ψ 0 ; 4. Time evolution -Implementation of Eq. 4; 5. Measurement -Retrieving physical information by measuring the qubits, according to a sensible protocol; 6. Output -For each field configuration the algorithm will output the expected value of a random variable χ, which should be then medium-averaged over all m configurations. Fig. 1 Overview of the quantum circuit detailed in the main text. Single lines denote quantum channels while double lines denote classical ones. Above each line we detail the state being store in the circuit (see main text for notation). The denotes that the time evolution gates parameters are to be determined from the field A.

Encoding
We begin by discretizing the problem in position space, such that |x = |a s n , with a s the spatial lattice spacing and n = (n 1 , n 2 ) a two component dimensionless transverse vector, where each component can take integer values between 0 and N s − 1, with N s the number of lattice sites per dimension. The spatial cutoff is given by x max = a s (N s − 1, N s − 1). Also, the spatial discretization induces a lattice discretization in momentum space with |p = |a d q and a d = 2π a s N s the momentum space lattice spacing, with q = (q 1 , q 2 ) a two dimensional vector with each component also taking integer values between 0 and N s − 1 2 .
One can rewrite the Hamiltonian H in terms of the dimensionless Hamiltonian H = Ha s (see Appendix A) withP |q = q |q andX |n = n |n the dimensionless position and momentum operators. Also is the dimensionless energy factor. In what follows, position and momentum vectors are assumed to be given in this dimensionless basis.
With this discretization, the problem can be mapped to the qubits available in a quantum computer. For each spatial dimension, we use a register with n Q qubits (each qubit being equivalent to a 1/2-spin), such that we can generate 2 n Q = N s states. We use the QIS convention [24] to denote the single up spin state |↑ = |0 = [1, 0] T in the computational basis (with the last equality giving the associated vector representation) and |↓ = |1 = [0, 1] T . Then, any component of the vector |n can be represented by a product of many spins, in a binary basis (see Appendix A). The associated momentum state vector |q is obtained by applying a standard quantum Fourier Transform (qFT).

Initial state preparation
Given the above encoding, the first step in the algorithm consists in loading a desired template distribution by constructing the initial state |ψ 0 from the fiducial state |0 ⊗2n Q . The template is meant to represent the relevant physics of the hard scattering which generates the initial parton.
In this manuscript, since we are interested in extracting the jet quenching parameterq from the quantum simulation output, we wish to avoid contributions coming from initial state physics. Therefore, we shall mainly focus on the case where |ψ 0 = |p = 0 .
However, including a localized initial state distribution might be important for certain digitizations where one can not prepare the state |p = 0 exactly or if one is simply interested in studying how different production mechanisms influence the final state. Several strategies to prepare |ψ 0 from an integrable template distributions can be found in the literature [25][26][27][28]. Depending exactly on what |ψ 0 one wants to prepare, in principle, one can devise a routine which only requires O(n Q ) basic quantum gate operations.

Time evolution
After the initial state |ψ 0 has been prepared, we time evolve it for a time L, producing the final state |ψ L . The time evolution operator in Eq. 3 can be written in terms of the dimensionless Hamiltonian H and medium length L ≡ L/a s 3 .
Directly implementing U (L , 0) in terms of a quantum circuit is in general impossible. Rather, one decomposes the full evolution into a sequence of short time evolution steps. Here we do this by considering the simplest product formula [29], decomposing U as where we have effectively sliced time into N t steps, each with a length ε t ≡ L /N t . In each time step, the evolution operator is split into a short evolution according to H K , followed by an evolution in time with H A . Notice that during the time interval (k t · ε t , (k t + 1) · ε t ) the field A is taken to be constant, leading to the constraint ε −1 t ||∂ t H A (t)||; there exist algorithms [29] which circumvent this constraint, as well as other strategies (see for example [30][31][32][33]) to quantum simulate time dependent Hamiltonians with expected higher precision. Although the way one chooses to implement U is of critical importance to determine the efficiency and accuracy of the quantum circuit, since we are aiming to restrict our discussion to a more conceptual level, we limit our analysis to the simple product formula considered above, which has a Trotter error O(ε 2 t ). Let us now consider the k th t time slice of the evolution. As mentioned above, H K has a trivial action in the momentum basis, while H A can be simply written in the position basis; this justifies the decomposition of H taken in Eq. 7. Since these bases are trivially related by a qFT, one can simply first time evolve with H K , perform the transformation |p → |x , time evolve with H A and transform back to the |p basis, the generated state being the input to the k t + 1 th time slice; this strategy is illustrated in Fig. 2. The time evolution operator U K is diagonal in the |p basis thus one only needs to implement a circuit which generates a state dependent phase. This can be achieved using the algorithm introduced in [34], which we detail in Appendix B.
After performing the qFT, using standard implementations of the circuit [24], one has to compute the action of U A . Although this operator is diagonal in the |x basis, the value of the phase depends on the local value of the field A. Notice that here we assume that the quark is a color singlet; see section 4 for details on how to deal with non-trivial color evolution. In principle, one could again use a strategy similar to the one used to implement U K (see Appendix B). However, this assumes that one could construct N t oracles which quantum compute A(k t · ε t , x) for every x in each time slice. Since in general one does not have a closed form expression or a simple numerical routine to compute the field values, such an approach might not be possible. A more feasible approach would consist on first computing the field values for all positions and times. This would require evaluating the field O(N t × N 2 s ) times, which would defeat the purposes of the present strategy since it requires exponentially many classical evaluations of A. Nonetheless, we notice that in practice a small number of qubits n Q is needed to have a sufficiently good discretization (see section 5), and thus the actual number of field evaluations needed could in practice be performed by a classical computer.
Once one has evaluated all the relevant field values, they are stored in a classical memory (double lines in Figs. 1 and 2) which are loaded onto the circuit as parameters to the basic gates implementing Eq. 9. We illustrate this procedure in Appendix B, that requires solving a system of linear equations with N 2 s independent variables, which following the same arguments as above should be doable in practice, at least for near term small system applications 4 . Clearly the implementation of the operator U A would greatly benefit from native implementations of quantum diagonal gates, where each entry exponentiates a circuit input [35] 5 .
After performing this operation and transforming back to the momentum basis, this block is iterated until k t = N t , where the time evolution section of the algorithm terminates.

Measurement
Having prepared the state |ψ L = q ψ q L |q one could simply measure all the 2n Q qubits, obtain the probabilities |ψ q L | 2 for every q and reconstruct the underlying probability distribution. However, such a strategy requires a exponentially large number of measurements. This constraint is a direct consequence of the quantum nature of the simulation, absent from classical simulations where information can be easily retrieved.
In this section, we assume that the initial condition of the state was that of a quark with p = 0. In this case the coefficients |ψ q L | 2 are directly related to the single particle broadening distribution; see Appendix C. This statement is only true after having averaged over all field configurations, the so called medium average, which in our strategy is performed at the end of the algorithm. For each of the m field configurations one runs the algorithm the necessary number of times to extract the expectation value of some classical variable χ (to be detailed below). Then, one averages over all m expectation values where p A − = p A , the i superscript denotes a particular field configuration, running up to m, and . M denotes the average over field configurations while . QM denotes the (quantum mechanical) expectation value. The numerical value for m depends on field fluctuations. In jet quenching, these are typically Gaussianly distributed, following the prescription of the Mclerran-Venugopalan (MV) model [36,37] and are encapsulated in the field-field correlator [36][37][38][39]. We note however that in our approach, one is not constrained to assume the MV prescription, nor does one need to explicitly construct any field correlator. In addition, due to the formal similarities between jet quenching and saturation physics [40], the physical origin of A − , either generated from hot and dense quark gluon plasma, the initial glasma or from cold nuclear matter, is not constrained. This means that our approach should be able to explore the evolution of the jet quenching parameterq, both in time and in orthogonal spatial directions [41], for different medium models. The only practical constraint is that the larger the background field fluctuations become, the larger m must be, despite this only leading to a linear increase in cost for running the full algorithm.
We then focus the remaining discussion on the case of a fixed field configuration and how to extract q for that A − . We add an ancilla qubit to the circuit and perform the Hadamard test detailed in Fig. 3. One first transforms the ancilla, which can be either prepared in the state |0 or in the superposition 1/ √ 2(|0 + i |1 ), by the Hadamard gate H = H † , and then applies a unitary transformation V on the physical state if the ancilla is in the state |1 . Finally the transformation on the ancilla is reversed and one measures the qubit. We associate the measured value to a random variable χ which takes the values −1 if we observe the state |0 and +1 if the state |1 is generated. This strategy is not the only possible one, but it is particularly simple and inexpensive in terms of extra ancillas and number of gate operations.
One can show that if the ancilla is in the initial state |0 (see Appendix D), then On the other hand, if the ancilla is prepared in the state 1/ √ 2(|0 + i |1 ), we have that which when combined give access to both the real and imaginary parts of the expectation value of the unitary operator V . Let us consider first the case where V = V α = exp(iαP 2 ). Then and from which one extracts e iαP 2 QM , by definition. We also have that where 2k ≡ P 2k QM corresponds to the expectation value of the 2k power of the momentum operator. Eq. 15 can be viewed as the (even) moment generating function, and it is easily related to the cumulants of the underlying broadening distribution. Also, in the case where initial state effects are absent, a 2 d 2 =qL, where we inserted a 2 d to get the correct dimensions. Furthermore, one has the freedom to vary α such that, for small enough α, only linear variations are relevant Notice that the left hand side corresponds to a quantity readily extracted from the quantum computer, while the right hand side is written in terms of the physical jet quenching parameter. If one includes higher order α corrections, then one has access to the even moments of the momentum distribution and the respective cumulants. One can thus imagine varying α and from the observed evolution retrieving 2k moments via a numerical fit. Of course, such a strategy, on top of the additional polynomial cost in m, would increase the cost of running the algorithm by the number of α values to be explored.
If one is only interested in extractingq (which is the most relevant medium parameter for jet quenching), one could consider the unitary V = exp(iF (P 2 )), with F (P 2 ) = arccos(P 2 ). Then, for the case where the ancilla is initially set to |0 , we obtain X QM = ψ L | cos(arccos(P 2 )) |ψ L = 2 .
In principle, one could implement this protocol following [34] (see also Appendix B), provided an efficient arithmetic oracle could be constructed.

Treating color evolution
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 · T = A a T a = 1 2 A a λ a , where λ 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 ·U A 1 ·U A 2 · · · 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: |0, 0 ≡ |red = |R , |0, 1 ≡ |green = |G , |1, 0 ≡ |blue = |B and |1, 1 ≡ |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 λ 1 into the two qubit Hilbert space. The action ofλ 1 is thus to color rotate the quark state between the |R and |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 th t 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σ Z = diag(1, −1, 0, 0) 6 . Finally, to compute the exponential of the tensor product we notice that where |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 |0, 0 and |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⊗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 7 .

Numerical estimates for the circuit parameters
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, p 2 ∼qL ≡ Q 2 s . 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 ∼ O(10 fm) = O(50 GeV −1 ). The value of the jet quenching parametersq 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 O(0.1 GeV 2 fm −1 ) ≤q ≤ O(10 GeV 2 fm −1 ) [25,42,43]. The saturation scale Q 2 s is then approximately bounded by Q 2 Setting the ultraviolet momentum cutoff induced by the digitization p max. to be much larger than the saturation scale Q s , we obtain 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 µ an infrared (model dependent) regulator, related to the medium Debye mass; typically µ ∼ O(0.1 − 1 GeV) [25,44,45] 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 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 µ/|p max. |, which here we took |p 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 < 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 λ is much larger than the typical correlation length in the medium 1/µ. 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 λ ≤ L, leading to It is also typical to define the opacity of the medium as χ med. ≡ L/λ [46,47], corresponding to the expected number of in-medium scatterings. Therefore, it is natural to identify χ med. ∼ N t = L /ε 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 A is chosen.

Conclusions and Outlook
In this paper we have outlined a quantum simulation algorithm to extract the jet quenching parameter 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 [13], still not fully developed (see however [48,49]).
The overall algorithm requires 2n Q + l qubits (assuming one can re-use ancillas) and O(N t × polylogN s ) basic gate operations. However, there is an underlying classical cost coming from the m × N t × N 2 s evaluations of the gauge field. This is the major drawback of our strategy, since it is not guaranteed that the classical evaluations of 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 α. For an actual implementation in a NISQ device [50], 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 forq from jet quenching phenomenology. This would allow a better understanding of the merits of our quantum approach, compared to known classical methods.
Going beyond α 0 s 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. where x = d 2 x and p = (2π) −2 d 2 p and we provide the discretized version of the Fourier integrals. Using that where we used the closure identity which satisfy n|m = δ n,m , q p |q k = δ q p ,q k and n|q = N −1 s exp(−2πiN −1 s n · q). The Fourier transforms in this normalization take the form It is also natural to introduce the operators P = p/a d and X = x/a s , satisfyingX |n = n |n andP |q = q |q . Inserting this operator definitions into Eq. 2, one can extract the dimensionless Hamiltonian H = a s H, given in Eq. 5. The map to the 1/2-spin registers in the quantum computer is achieved by decomposing each component of the vector n = (n 1 , n 2 ) in the binary basis, e.g.
where n (i) 1 ∈ {0, 1} and we assume that there are n Q qubits available, such that 2 n Q = N s is total number of possible states. If n (i) 1 = 0 then we associate a qubit in the state |↑ = |0 to it; conversely if n (i) 1 = 1 we assign |↓ = |1 . Then, for example, the ket state |n = |3, 3 with n Q = 2 is given by two registers storing the overall state |1, 1 ⊗ |1, 1 . Following Eqs. A.8 and A.9, the transformation between the momentum and position basis is achieved by applying a standard quantum Fourier transform (qFT) [24].
Finally, in this appendix and in the main text we have restricted ourselves to considering lattices over positive integer values of n and q. In an actual implementation, one would have to consider signed values, since, in general, there is no condition that physically constrains the system to positive values. In principle, signed values can be dealt with by, for example, including an extra qubit that stores the sign of the state (similar to the encoding used in [25]) or using a two's complement encoding. This caveat requires one to modify circuits we detail in the main text to accommodate for the new encodings. In general, one should be able to do this without incurring in an exponential number of extra qubits or basic gate operations 9 , nor does it lead to any new conceptual challenge that must be addressed. As such, we do not further discuss this issue in the paper and leave it for a future work where we tackle a detailed implementation of the algorithm. where x j ∈ {0, 1}. Acting on a single qubit the above operator has non-zero matrix elements 0| R j (s K ) |0 = 1 and 1| R j (s K ) |1 = exp(−is K 2 j ); clearly stringing together l of such operators with increasing values of j results in a multi-qubit operator implementing the desired transformation, i.e. R(s K ) |x = exp(−is K x) |x . The final step -a 3 -consists in erasing the ancilla register back to the state |0 ⊗l , which can be achieved by applying the Hermitian conjugate circuit used in step a 1 .
The implementation of the operator U A could be done following exactly the same strategy as just described. However, as mentioned in the main text, this would require having a way to construct efficient quantum oracles that, for each time t, given |x output |A(t, x) . We expect that for most cases, this will be difficult to do.
As an alternative we consider that one is handed a list of N t × N 2 s values, describing the field values at all the relevant spacetime points. Then one can implement U A by stringing together 2n Q single qubit gates R α,β ≡ diag(exp(iα), exp(iβ)). Such gates can be written as the product of the exponential of Pauli gates and the R j gate. In one spatial dimension and for n Q = 1 one simply has for the k th t time slice that α k t = −g t A(k t · ε t , 0) and β k t = −gεA(k t · ε t , 1), where the sub-index denotes the time slice and there are only two spatial lattice points (|0 and |1 ). If we now consider n Q = 2 but still a single spatial dimension, the respective time evolution operator would be obtained by for each time slice. By solving the associated system of linear equations, one can map {α, β, σ, γ} to {A(x)}, which can be done offline for any t in a classical computer.
Appendix C: Relation between |ψ L and the single particle momentum distribution In this appendix we relate |ψ L to the broadening distribution. The single particle broadening probability for observing a quark with momentum k due to interactions with the medium for a time L is given by [15,20,21]  The above medium average is usually performed by detailing the non-trivial correlators of the background field, in jet quenching typically the MV/Gaussian prescription. Using this further assumption, one can then write the broadening distribution in terms of a so called dipole cross-section, which is typically constrained to recover the Coulomb form at short distances and to have a model dependent form in the infrared [44,45]. It is not difficult to check that, in the strict eikonal limit, where H = H A , the circuit detailed in the main text mirrors the P distribution. For clarity, we ignore the details in the implementation of the time evolution operator; additionally, we assume that the initial state is that of a quark with zero transverse momentum |ψ 0 = |p = 0 .
In this scenario the circuit simplifies significantly since all but an initial and a final qFT cancel out. Then the system state transforms as The probability of measuring the state |k , P k , is simply given by