Quantum simulation of in-medium QCD jets: momentum broadening, gluon production, and entropy growth

Jets provide one of the primary probes of the quark-gluon plasma produced in ultrarelativistic heavy ion collisions and the cold nuclear matter explored in deep inelastic scattering experiments. However, despite important developments in the last years, a description of the real-time evolution of QCD jets inside a medium is still far from being complete. In our previous work, we have explored quantum technologies as a promising alternative theoretical laboratory to simulate jet evolution in QCD matter, to overcome inherent technical difficulties in present calculations. Here, we extend our previous investigation from the single particle $|q\rangle$ to the $|q\rangle+|qg\rangle$ Fock space, taking into account gluon production. Based on the light-front Hamiltonian formalism, we construct a digital quantum circuit that tracks the evolution of a multi-particle jet probe in the presence of a medium described as a stochastic color field. Studying the momentum broadening of the jet state, we observe sizable sub-eikonal effects by comparing to eikonal estimates. We also study the medium-induced modifications to the gluon emission probability, which exhibit small corrections compared to the vacuum splitting function. In addition, we study the time evolution of the von-Neumann entropy associated with the quark component; we find that the exponential of the entropy grows linearly in time for the bare quark but super-linearly when taking into account gluon emission.


I. INTRODUCTION
Heavy-ion collisions at the Relativistic Heavy-Ion Collider and the Large Hadron Collider, produce collimated particle showers originated from highly energetic quarks and gluons, known as jets, that evolve simultaneously with the hot and dense quark gluon plasma.A similar scenario is expected in future deep inelastic scattering experiments, where these sprays of particles traverse the cold nuclear matter target.
Jets can resolve the underlying medium at different energy scales and thus offer an optimal probe to study the structure of QCD matter.Their evolution in these environments is characterized by sizeable modifications to the jets' structure, which are reflected at the level of the final state distributions.Phenomenologically, these result in the broadening of the transverse momentum distributions, due to instantaneous interactions with the medium, and the alteration of the radiation pattern, leading to an excess of energy flowing at large angles; for recent reviews on medium-induced jet modifications see [1][2][3][4][5].To the present day, the theoretical study of all these effects has been mainly constrained to lower orders in perturbation theory [6][7][8], with a limited number of higher order calculations being available [9,10].Compared to their vacuum counterpart, the slower progress seen in this research field is mainly tied to the highly complex multi-particle interference pattern determining parton fragmentation in matter.
More recently, it has been argued that novel advances in quantum information science could be used to leverage our understanding of in medium jet physics [11,12].In particular, future large-scale and fault-tolerant digital quantum computers can potentially offer a platform to efficiently simulate large quantum systems in real time, using the quantum simulation algorithm.These devices capture the dynamics of the target quantum system by evolving a controlled finite-dimensional quantum system whose dynamics can be engineered appropriately.Nonetheless, their practical implementation is still complex and problem-dependent.
In our preceding work [13], we provided a quantum simulation protocol to simulate the real-time evolution of a single hard parton in the presence of a stochastic background.This can be thought of as describing the propagation of the leading parton inside the jet, and it allowed us to establish the basic aspects of this approach.In this work, we extend the strategy to include gluon radiation, thus allowing the jet to have a non-trivial structure.
Our quantum formulation is based on a nonperturbative light-front Hamiltonian approach, the timedependent basis light-front quantization (tBLFQ) [14]; see its various applications in Refs.[15][16][17][18][19][20].This approach allows us to quantize the QCD Hamiltonian and perform real-time simulations at the amplitude level, which is natural and well-suited for applications in quantum computers.In particular, we can exactly track the jet state in time and extract observables by computing the expectation values for appropriate operators.

arXiv:2307.01792v1 [hep-ph] 4 Jul 2023
Quark jet evolution in a colored medium using tBLFQ has already been simulated on classical computers: first in the |q Fock space [18], and later extended to the |q + |qg space [19,20].These studies showed the interplay between coherence and multiple scattering in gluon emission from nonperturbative perspectives.Our preceding work [13] provides a quantum implementation in the |q space and further investigations in jet momentum broadening.This work takes a step forward by extending to the |q + |qg space.We build a digital quantum circuit that can track the evolution of a jet state in the presence of a medium background field.We focus on the momentum broadening of the jet state and its branching pattern, both in vacuum and medium.We also discuss modifications to the single particle entropy growth due to radiation.
This manuscript is organized as follows.In Sec.II, we review the formulation of a jet evolution in a dense medium within the quark and quark-gluon Fock sectors using light-front Hamiltonian formalism; we then build the quantum simulation algorithm of the jet evolution.In Sec.III, we present numerical results of our approaches via quantum simulation of the jet using Qiskit.In Sec.IV, we summarize our current results and discuss the future avenue of this work.

II. METHODOLOGY
In our preceding work, Ref. [13], we quantum simulated the evolution of a single particle through a colored medium.Our method was based on the tBLFQ, a numerical non-perturbative light-front Hamiltonian approach developed to study real-time problems; see Refs.[18,19] for further details and applications using classical methods. 1 In this work, we extend the Fock space to |q +|qg , thus including gluon emission and absorption.On a classical computer, the simulation of this process has been studied in Refs.[19,20], also within tBLFQ.
In the following, we first briefly review the basics of tBLFQ to simulate the in-medium jet evolution process in the |q + |qg Fock space in Sec.II A and then detail how to apply this method using the quantum simulation algorithm in Sec.II B.

A. Jet evolution in the light-front Hamiltonian formalism of tBLFQ
We consider the propagation of a highly energetic massless jet with light-front momentum p = (p + , p − , p), moving close to the light cone along the x + direction. 2his hard probe evolves in the presence of a dense medium, which can be boosted to move along the x − direction.The quark interacts with the medium over a finite distance in light-front time x + = [0, L η ].This process is illustrated in Fig. 1.The dynamics of this system are set by the QCD Lagrangian in the presence of an external field, where is the field strength tensor, D µ ≡ ∂ µ +igC µ the covariant derivative, and C µ = A µ + A µ is the sum of the quantum gauge field A µ and the background gluon field A µ .In this work, we truncate the Fock space of the jet to the leading two sectors, |q and |qg , such that the full quantum state can be written as where ψ q and ψ qg represent the respective Fock amplitudes.The light-front Hamiltonian can be obtained following the canonical light-front quantization formalism [5,18,19,21] via the standard Legendre transformation, in the light-cone gauge of A + = A + = 0, tudinal coordinate, and x = (x 1 , x 2 ) the transverse coordinates.The letters in bold, such as x, denote transverse vectors, while their magnitude is denoted by x ⊥ ≡ |x|.The non-vanishing elements of the metric tensors g µν and gµν are, g Here, P − KE stands for the kinetic energy part where P − KE,g and P − KE,q are the respective kinetic energies for the dynamical gluon and quark.The second term V qg is the interaction between the quark and gluon, The third term V A (x + ) includes the interaction of the background field with the quark and that with the dynamical gluon, As in our preceding work [13], here we also take McLerran-Venugopalan (MV) model [22,23] to describe the field A that accounts for the background medium and we make use of high energy (eikonal) limit, where p + |p| ≡ p ⊥ , p − .Note, nevertheless, the full Hamiltonian method allows us to go beyond the formal eikonal limit of p + = ∞.The color charge density of the medium is assumed to have a Gaussian and local correlation function where we use • • • to denote the average over medium configurations, and µ controls the strength of the medium.The saturation scale is defined as with the fundamental Casimir C F = (N 2 c −1)/(2N c ).The field is then solved from the reduced classical Yang-Mills equation, where the gluon mass m g is introduced to regularize the infrared (IR) divergence in the field [24].The time evolution of the quark jet, as a quantum state, obeys the time-dependent Schrödinger equation.Written in the form of path-ordered exponential, it reads, where T + is the light-front time ordering operator, and |ψ(x + ) the quantum state of the jet at time x + . 3We solve this equation non-perturbatively by decomposing the time-evolution operator as a sequence of small time steps in the light-front time x + , where x + k = k L η /N t is the intermediate time and N t the total number of time steps.

B. Quantum simulation algorithm
The digital quantum simulation algorithm [25][26][27][28][29] typically involves five generic steps: input, encoding, initial state preparation, time evolution, and measurement.Here, we extend the algorithm developed in Ref. [13] for a jet in the |q Fock space to the |q + |qg sectors.

Basis encoding
We choose the eigenstates of the kinetic energy part of the Hamiltonian P − KE as the basis states, as formulated in the Ref. [19] for |q + |qg .This basis choice is convenient in studying the momentum broadening of the jet state.
We start by considering a generically truncated Fock space of the quark jet state.The full Hilbert space of this theory can be formally decomposed as a tensor product over all single particle subspaces [30,31].Each Fock sector can have a finite projection in each one of these subspaces.Let us consider a generic multi-particle Fock sector in the quark jet state, |q . . .q q . . .qg . . .g , in which the number of quarks is one more than that of the antiquarks.The basis state is in the form of |β q...q q...qg...g = |β q ⊗ . . .⊗ |β q ⊗ |β q ⊗ . . .⊗ |β q ⊗ |β g . . .⊗ |β g .Each single particle state carries five quantum numbers where p + is the longitudinal momentum, {p x , p y } the transverse momenta, λ the light-front helicity, and c the color index.For a basis state in the truncated Fock space with up to n + 1 quarks, n anti-quarks and m gluons, in which N q is the total number of single quark basis states, N q for antiquarks, and N g for gluons.Each register |e qi , |e qi , |e gi encodes the occupancy of quarks, antiquarks, and gluons in the i-th single-quark basis state β i q , and they satisfy Particle exchange symmetry should be satisfied accordingly and implemented on the physical state; see discussion and references in [30] for further details.The encoding of the single particle basis |β l can follow the strategy described in our previous work [13], which simply enumerates all the quantum numbers in the phase basis.
Following this construction, let us describe in detail the encoding for the Fock space |q + |qg .We can save on the number of qubits with the following arrangement and simplification.According to the strategy sketched in Eq. ( 13), we need one qubit to encode the occupancy status of the gluon, e.g., |e g = |0 for |q and |e g = |1 for |qg .We extend |e g to multi-qubits to also encode the p + quantum number of the gluon, denoted as |ζ .In the helicity space, we make a simplification by considering only the helicity-non-flip term, i.e., λ q = λ g =↑.Note that the quark is taken to be massless, so only the quarkhelicity-non-flip terms in V qg are non-zero, then the chosen configuration is the dominant contribution when the emitted gluon is soft.We therefore do not need extra quantum registers for the helicity space.The remaining quantum numbers to be encoded are the transverse momenta and the color indices.Therefore, the complete basis encoding for any basis state in the |q and |qg Fock sectors is written as In the following, we recapitulate the encoding scheme for the transverse momentum and extend that for the color sector following our previous work [13], and elaborate on the construction of the |ζ register.

i. The transverse dimension
We formulate the transverse space as a twodimensional square lattice.Both lattices span a size of 2L ⊥ and a number of 2N ⊥ sites per dimension such that the lattice spacing is a ⊥ = L ⊥ /N ⊥ .
We impose periodic boundary conditions on the lattice, such that this position space and the reciprocal momentum space are related by a discrete Fourier Transform, which on the quantum computer can be implemented via a quantum Fourier Transform (qFT ).An arbitrary momentum state vector ii.The color space We consider N c = 2, such that there are N c = 2 color degrees of freedom for the quark and N 2 c −1 = 3 for the gluon.We use one (two) qubit (s) to encode the color index of the quark (gluon).The quantum registers for the color space are explicitly, iii.The |ζ register We compactify x − to a circle of length 2L (i.e., x + to a circle of length L), and impose periodic boundary conditions for bosons and anti-periodic conditions for fermions, such that the longitudinal momentum p + is discretized, where the zero mode for the gluon is excluded. 4he total p + is preserved by the Hamiltonian.Consider a quark jet with a definite P + , then each basis state has the same total longitudinal momentum, i.e., i p + i = P + in which i enumerates the Fock particles.As such, in the |q sector, p + Q ≡ P + ; in the |qg sector, p + q + p + g ≡ P + . 5We introduce the total longitudinal quanta K such that where K is a positive half integer.
We combine the longitudinal encoding with the gluon occupancy using the quantum register |ζ .
The ζ = 0 state on the register encodes the |q state with k In the latter case, the value of ζ relates to the gluon's longitudinal momentum fraction as where we use * to represent the absence of the gluon in the |q sector.
This basis space contains a total number of 2 3 K (2N ⊥ ) 4 basis states, which scales as a power four with the lattice size.On a classical computer, the problem quickly deteriorates when additional Fock sectors are included, requiring at least (2N ⊥ ) 2n resources for n particles.On the quantum circuit, however, we would only need a total of n Q = (7 + 4 log 2 N ⊥ + log 2 K ) qubits for this problem, dramatically reducing the number of resources.Nevertheless, efficient gate simulation is still needed and not necessarily always guaranteed.
Using this qubit encoding scheme encapsulated by Eq. ( 14), one can prepare any initial state |ψ 0 as a superposition of basis states.Though arbitrary state initialization may be difficult, the preparation of many useful choices of initial states is feasible [33,34].Since we are mostly interested in studying the jet evolution in momentum space, we neglect the initial state effects, and take the initial state |ψ 0 to be a zero transverse momentum and even-color single-quark state, unless specified otherwise.

Gate encoding and time evolution
We implement the product formula decomposition by splitting the evolution along the x + direction into N t time steps, each with a duration of δx + = L η /N t , as in Eq. (11).Note that the time-dependence of the Hamiltonian is from the background field A. We slice the medium into N η layers along x + [13, 18-20, 35, 36], such that the time duration for each layer is τ ≡ L η /N η .A schematic representation of the circuit is presented in Fig. 2. We generate the values of the background field A beforehand in a classical computer, as in our preceding work [13].The details on how to calculate the field numerically are given in App. A.
We consider and compare two treatments on the evolution operator.For a Hamiltonian that is constant in time, the evolution operator reduces to an ordinary exponential, which can be evaluated directly.This is the case for our Hamiltonian within each layer of the medium.Taking τ as the size of the time step δx + , such that N t = N η , we therefore have the single-step evolution operator as

i. Direct exponentiation
in which k = 1, 2, . . ., N η .Here, P − = K +V qg +V A denote the matrix form of the corresponding operator evaluated on the basis space. 6We evaluate the Hamiltonian matrix elements in the transverse momentum basis space, such that the K part is diagonal and the V qg part is off-diagonal and sparse.However, the V A term is more complex since it has scattered elements in momentum space.Alternatively, one can perform the equivalent calculation in the transverse position basis, which favors the evaluation of V A but complicates the evaluations of K and V qg .
The major advantage of this treatment is that the evaluation is exact up to numerical accuracy.Nonetheless, the disadvantage of this treatment is the large time complexity of obtaining the Pauli strings of these respective non-diagonal operators; though for small-sized problems, it is feasible to take this treatment.
ii. Alternating exponentiation (mixed-space simulation) When the time step is sufficiently small, the Hamiltonian within a single step can be considered as constant, and one can further "trotterize" the single-step evolution operator as a product of operations with the different components of the Hamiltonian.In this way, one can factorize the singlestep evolution operator into a series of unitary operators from different components in the Hamiltonian and boost the computational efficiency, e.g., Refs.[13,19].We split the single-step evolution operator as the following, in which k = 1, 2, . . ., N t .Note that here δx + ≤ τ , i.e., N t ≥ N η .In practice and in our simulations, the time step δx + is taken sufficiently small (i.e., N t sufficiently large) to ensure the result is convergent 6 Full expressions of the matrix elements are found in Refs.[18,19].
when comparing to smaller δx + s; discussions on the convergence of N t and N η can be found in App.D.
The evolution of the kinetic energy and gluon emission/absorption, as the operation with K + V qg , is performed in the momentum space; whereas the evolution with the medium, that with V A , is performed in the position space.Since the background field A(x + k , x) is diagonal in the transverse position space, this mixed-space evolution approach is the most economical way of evaluating the Pauli terms and could potentially extend to larger lattice and more Fock spaces.Detailed resource cost comparison between the Pauli term evaluations in momentum/position space is included in App. C. The basis transformation between the momentum and position spaces can be performed by applying a quantum Fourier transform qFT to {|n } → {|k }, and its inverse qFT −1 to {|k } → {|n } per each transverse dimension of the quark and the gluon.
To perform the simulation using either of the aforementioned treatments, one needs to implement the quantum gates of an operation in the format of e iHδx + (with M a Hermitian operator).In our previous work [13], the matrix elements of e iHδx + can be obtained exactly using the properties of the exponential of the Pauli vector and then transcribed to unitary gates using the quantum Shannon Decomposition [37].Though this approach works well for simulating the jet in the quark Fock space, it is inconvenient to obtain the exact exponential of the type of Hamiltonian in this work.
Instead, we find the corresponding Pauli terms of H first and then the associated quantum gates, since there is a direct correspondence between the Pauli exponentials and the quantum gates [28].To obtain the Pauli terms, various strategies can be adopted (see App. C for examples and discussions), and we used the sparse matrix projection methods to take advantage of the property of the Hamiltonian matrix.To time evolve our Pauli terms, we use PauliEvolutionGate class provided by Qiskit [38], which automatically maps the Pauli operators to quantum gates.For small problem sizes, we can perform exact operator evolution via matrix exponentiation; for large problem sizes, we can use the Lie-Trotter formula [39] to approximate the exponential of non-commuting operators at first order.For higher-order approximations, we can use the Suzuki-Trotter product formula [40].All these methods are conveniently implemented in various Synthesis classes [28,41] in Qiskit.We studied the performance of these unitary exponential implementations, especially using MatrixExponential, LieTrotter, and SuzukiTrotter, and found that their performances are almost identical with each other at our problem scale.

Measurement
We extract the information about the final quantum state by directly measuring the prepared state.In practice, since we work in small lattice sizes, such an approach is the most efficient.Note however, it is not always necessary to measure the full quantum state.For example, to obtain the induced gluon probability, one can measure the |ζ quantum register on a log 2 (K)-bit classical register alone, greatly reducing the number of measurement shots needed.
While most of the results presented in this work use the shot-based QasmSimulator backend to extract physical observables such as momentum broadening, we also use the StatevectorSimulator backend to capture the exact quantum state, serving as a benchmark.Though not practical on real quantum devices, it allows us to study the information flow in the evolution of the quark jet in medium.Estimating the entropy directly on the quantum circuit is generally difficult [42][43][44], and requires full-fledged fault-tolerant quantum computers in the future.

III. QUANTUM SIMULATION RESULTS
In this section, we study the quantum simulation results for the evolution of a jet in a dense stochastic medium, using the light-front Hamiltonian formalism and quantum simulation method introduced in the preceding sections.Specifically, we focus on the momentum broadening of the jet, the gluon emission, and the entropy growth, for several backgrounds with different medium strengths.We perform the simulations using the ideal QASM simulators from Qiskit.
For the simulations, we take the transverse lattice with N ⊥ = 1 for the |q + |qg system (we will also use larger N ⊥ when examining the |q system), and the total longitudinal momentum quanta K = 3.5.Although these numbers are small, it still provides us with a two-by-two transverse lattice for both the quark and gluon single particle states, allowing investigation on the effects of momentum broadening.With K = 3.5, we are also able to examine the distribution of the longitudinal momentum.We take L ⊥ = 32 GeV −1 = 12.6 fm.The duration of the medium is taken to be L η = 50 GeV −1 = 9.87 fm.We take the layer number to be N η = 4; one can find the discussions on the convergence of N η and the evolution time steps N t in Appendix.D. The IR regulator for the medium is m g = 0.8 GeV.More details on the determination of parameters for a proper lattice and medium can be found in our previous works [13,19].The total number of qubits required in this setup is therefore n Q = 9 according to the encoding scheme in Sec.II B 1.
Since we are mostly using the shot-based quantum simulator, we make sure a sufficient number of counts are used to sample the true probability distribution.Unless stated otherwise, we always use 819200 shots, which proves to be more than enough to take into account of the noise from statistic sampling for a 9-qubit simulation [13].The uncertainties (i.e. standard deviations) provided on our plots are therefore exclusively related to the medium field fluctuations arising from using a stochastic medium in the MV model.With these in mind, we will present our main results in the following.

A. Momentum broadening
Transverse momentum broadening is an important observable to understand the evolution of the jet inside the medium.We examine the square of the transferred momentum ∆ p 2 ⊥ (∆x + ) at various medium strengths of g 2 µ, which simplifies to p 2 ⊥ (∆x + ) when the initial state has a zero transverse momentum.
In the eikonal limit, p 2 ⊥ (x + ) of a single particle is linear in time, and the proportionality constant can be interpreted as the quenching parameter q.We have provided the explicit expression in the chosen basis representation in our previous work, as in Eq.( 19) of Ref. [13]. 7 We will use the eikonal expectation to verify our simulation results in the eikonal limit and examine non-eikonal effects by studying the deviation from them.
With the final jet probability distribution extracted from the quantum simulation, we are able to reconstruct the total transverse momentum of the jet, given as where P |q (P |qg ) is the probability of the state in the quark (quark-gluon) Fock sector.With a zero momentum initial state, p 2 ⊥ indicates the broadening effect exclusively due to the medium.By comparison, p 2 ⊥ = 0 in vacuum due to momentum conservation.The periodic boundary conditions of the lattice are taken into account 7 Here, we write out the analytical expectation at the special case of N ⊥ = 1 in order to compare with the simulation results.The specialty of the phase space at N ⊥ = 1 is that the lattice UV and IR cutoffs estimated in the usual way as λ U V = π/a ⊥ and λ T R = π/L ⊥ would be the same, then the analytical formula for general N ⊥ no longer hold.One should instead, treat the p ⊥ integral as a sum over the full discrete space of p ⊥ .In this way, we get  when summing the momenta of the quark and gluon state; see the prescription in Appendix C of Ref. [19].
We first verify our method by keeping only the medium interaction term in the Hamiltonian, V A = V qA + V gA .This setup corresponds to the process in the eikonal limit of p + = ∞, therefore the expectation value p 2 should agree with the eikonal expectation, e.g., Eq. ( 21).Specifically, we assign the initial state as both a single quark and a single quark-gluon state with total transverse momentum p = 0,8 and we put in the Hamiltonian V qA and V gA separately and in-combined.We present in Fig. 3 the results of the final state p 2 at various saturation scales Q s .The obtained simulation results agree with the expected eikonal analytical results.Similar to the single quark results shown in our previous work [13], the p 2 exhibits increased uncertainty at larger saturation scales, which is related to the larger Gaussian width in constructing the stochastic background fields; see also on the background field in App. A. In addition, the p 2 starts to bend as the saturation scale increases, as a result of the lattice admitting a UV cutoff of π/a ⊥ .
We then perform the simulations with the full Hamiltonian P − = K + V qg + V A .The eikonal approximation is relaxed by letting the jet state have finite energy, p + = 1, 1000 GeV.We assign the initial state as a single quark with p = 0, k + = K, and use both the |q and the |q + |qg Fock space for the simulation.The results of the final state p 2 at various saturation scales Q s is presented in Fig. 4. We have confirmed that the two simulation treatments discussed in Sec.II B 2, the direct and the alternating exponentiation, led to the same results.In the figure, the results in the eikonal analytical limit in the |q and |qg Fock spaces according to Eq. ( 21) are provided in the dashed and dotted lines for comparison.We find that the p + = 1000 GeV result overlaps with the eikonal limit for the |q with uncertainties taken into account;9 this is because both the kinetic energy and the gluon emission contribution are highly suppressed at the near eikonal limit (p + = ∞), and therefore the occupancy in the |qg sector is negligible. 10By contrast, the p + = 1 GeV result in |q + |qg lies between the two eikonal limits, whose deviation from the single quark's eikonal expectation indicates non-eikonal effects due to gluon emission.A simple and intuitive understanding is that the inclusion of the |qg sector enlarges the phase space, and as a result enhances the momentum broadening effect [19].

B. Gluon production
With the quark jet formulated as a superposition of |q and |qg states, it is interesting to study the gluon production through the evolution.In particular, we study the evolution of the probability of the jet in the |qg sector, i.e., P |qg .Furthermore, we examine the distribution of the gluon's longitudinal momentum fraction z g .We obtain the z g distribution by performing projective measurement on the |ζ register.For example, with K = 7/2, we should have 4 different longitudinal modes across the |q and |qg Fock sectors: where the total longitudinal momentum quanta of each mode is always K. The possible longitudinal momentum fractions of the gluon can be read conveniently as z g = k + g /K = {0.29,0.57, 0.86}.To observe the probability distribution throughout the evolution time L η on the quantum simulator, the same simulation is repeated for different x + to extract the corresponding probability.

Vacuum case
To better compare the medium corrections, we first present the simulation of the initial quark jet in the vacuum, which can be achieved by simply turning off the V A term in the Hamiltonian while keeping the K and V qg terms.
In Fig. 5, we show the total probability of the |qg Fock sector, P |qg , as a function of evolution time x + with and without the kinetic energy term, and using two different coupling strengths, g = 1 and g = 10.The setup of having only the gluon emission term V qg (without K) in the Hamiltonian, though not physical, is important to help understand its effect.For the size of the problem being simulated, it is also feasible to diagonalize the The Hamiltonian contains the gluon emission term, with and without the kinetic energy term, as specified in the legends.Each curve is obtained from both quantum simulations and classical diagonalization, and the two are in agreement.
Hamiltonian and obtain the eigenstates, therefore knowing exactly the evolution of a given initial state.The quantum simulation results agree with the diagonalization results, which help verify our quantum simulation algorithm.In both cases of g = 1 and 10, the |qg probability oscillates periodically under just the V qg term, but this behavior is broken with the inclusion of the kinetic term, as expected.The effect of including the kinetic energy term is akin to the existence of an energy-dependent phase for different states, leading to a decoherence effect when summing over states.Similar behavior was seen in the classical simulations done in Ref. [19].
The comparison between the results with g = 1 and g = 10 is also interesting.In the pure V qg case, the oscillation frequency is proportional to g and the amplitude is 1.In the V qg + K case, the amplitude of the Figure 7.The probability of the quark state at different p + configurations (characterized by Fock sector and zg) in the vacuum at selected time instances of x + = 2.5, 10.0 fm.The solid lines are fits to the given functional form up to an overall constant.
decohered oscillation is much larger with the stronger coupling.This is expected by noting that the oscillation amplitude is approximately proportional to the square of the ratio between the averaged V qg and the K terms [19].We note that the oscillation amplitude is small at g = 1 in Fig. 5(a) for the lattice that we are using.In principle, one can increase the lattice size N ⊥ to obtain larger oscillation amplitude; however, the simulation would be more expensive.For the purpose of this work, we will present our results using the larger coupling strength g = 10 in the V qg term unless specified otherwise.In Fig. 6, we present the evolution of the probabilities of different p + states, including the |q sector and the different segments of the |qg sector characterized by the gluon longitudinal momentum fraction z g .The |qg modes with the smallest z g dominate, as having a more rapid initial growth and a larger oscillation magnitude compared to the other two.Since we are simulating the spin non-flipping case, we expect the distribution of the longitudinal momentum fraction roughly proportional to the reduced splitting function , according to the Hamiltonian matrix element.Note that this splitting function matches the leading order q → q + g Altarelli-Parisi splitting function [45][46][47], up to an integration measure which has to be included at the cross-section level.We numerically demonstrate a good agreement between data and the reduced splitting function P q→qg in Fig. 7, up to a x + dependent state normalization constant.

Medium case
The same observables can also be computed for the case of in-medium propagation and directly compared to the vacuum scenario, allowing us to visualize the modifications to the jet fragmentation pattern.In Fig. 8, we present the probability of the |qg Fock states for an initial quark jet going through the colored mediums. 11Specifically, we used two sets of mediums with g 2 µ = 0.1 GeV 3/2 and 0.2 GeV 3/2 .Note g = 1 in the medium term V A whereas g = 10 in the V qg term.Since the number of momentum modes is small, the decoherence among different modes is not sufficient to suppress the oscillation, even at late times.It is therefore Figure 9.
The probability of the quark state at different p + configurations (characterized by Fock sector and zg) in medium with g 2 µ = 0.2 GeV 3/2 at selected time instances of x + = 2.5, 10.0 fm.The solid lines are fits to the given functional form up to an overall constant.
hard to conclude whether the medium induces or suppresses the production of radiation conclusively.At this level, it is then not possible to fully comment on the relation between our numerical results and the Landau-Pomeranchuk-Migdal (LPM), which is known to determine the gluon radiation spectrum in the medium [5].Much larger lattices are necessary to further investigate the effect, which is beyond our current scope of study. 12 In Fig. 9 we compute the splitting function in the medium and compare it to the estimated vacuum splitting function P q→qg ; see also Fig. 7.For the simulated time duration, we observe that the in-medium data points are compatible with the vacuum splitting kernel.Of course, the possible existence of deviations is shadowed by the small number of data points and the unitarity constraint.We observe larger deviations with respect to the quark probability, with the medium leading to a suppression of the single quark sector.Due to probability conservation, this indicates an excess in the gluon production due to the propagation in the medium.This is in agreement with previous studies, where it is observed that the medium can promote the production of a large amount of radiation at larger.However, understanding the origin of this radiation requires making a differential measurement in transverse space, which requires a larger lattice to resolve the distribution. 12In a closely-related classical study using as large as N ⊥ = 16 and Nη = 4, it is found that the gluon probability depends on the medium strength and is enhanced at the presence of the medium [20].

C. Quark entropy
In our simulations we can directly access the final jet state, and therefore easily compute the associated entropy.Here we are particularly interested in computing the entropy of the reduced density matrix of the single quark.This provides a simple and straightforward way to understand the role played by radiative corrections [48][49][50][51] in in-medium jet evolution. 13n what follows, we study the von Neuman (vN) entropy S vN of the quark component of the jet.At leading order in the strong coupling, and using the single momentum mode initial condition we consider in this work, it can be shown that the quark entropy is related to the classical phase space explored by the state [52].Since for a single particle p 2 ⊥ (t) ∝ qt, one has that entropy should grow logarithmically with time.However, once radiation is included, the growth rate should increase.
The vN entropy of the quark component is defined as The reduced quark density matrix is understood as being averaged over medium configurations, i.e., ρ(x + ) .This averaging removes the medium's degrees of freedom.For a jet state in the |q space, the single-event density matrix is given by in which |ψ(x + ) is the state vector.For a jet state in the |q + |qg space, we trace over the gluon degrees of freedom, In practice, at the level of the circuit introduced in Sec.II B 1, this can be achieved by performing projective measurements over the |ζ ⊗ |g registers, or equivalently taking the partial trace of the full density matrix with partial trace in Qiskit [38].
We study the entropy of the jet state formulated in both the |q and the |q +|qg spaces.In the single parton case, the entropy is expected to behave as log 2 (1 + ax + ) according to Ref. [52], in which a is a parameter related to the average transverse momentum square acquired due to the interactions with the medium.In Fig. 10(a), we present the simulation results for increasing lattice sizes of N ⊥ = 1, 2, 4 at fixed medium strength g 2 µ = 0.1 GeV 3/2 and finite energy p + = 1 GeV.We can see a larger entropy growth with the lattice size, which is expected as the phase space becomes larger.We also notice the apparent logarithmic growth for the different parameter sets used as a function of x + .To further examine the dependence, we fit the data points to the expected functional form above, using a as a free fitting parameter for the different N ⊥ .Since a is related to the average momentum transfer experienced by the quark, one expects it to grow linearly with (g 2 µ) 2 , i.e., q.In Fig. 10(b) we show the evolution of the fitting parameter as a function of (g 2 µ) 2 for the different N ⊥ values considered.Indeed, we observe that the evolution for each lattice is reasonably described by linear regression.
When the gluon is included, the entropy can not only grow due to momentum diffusion but also as a consequence of the recoil experienced by the quark due to the gluon production.As a result, one should expect a larger growth of the associated von-Neumann entropy.We consider two mechanisms in describing such growth.
The first possible mechanism is that including the gluon production can lead to a larger effective q, and therefore a larger value for the fitting parameter a.We test this hypothesis by fitting the same functional form log 2 (1 + ax + ) to the quark entropy computed for the two-parton scenario, and found that such a fit can not properly describe the results obtained from the simulation.This suggests that for the quark entropy, the effect of having gluon radiation can not be reduced to having a larger effective value for q.
We then tend to the second possibility that the production of radiation can lead to an accelerated entropy growth reflected in the (anomalous) time exponent.To examine, we consider the functional form log 2 [1+(ax + ) b ], with b > 1 to fit the simulation results for comparison.We present the data points the fits in Fig. 11, including both the case where the jet evolves in the medium and in the vacuum. 14We observe that this functional form can properly capture data points.The fitting parameters for the medium and vacuum points are compatible and show that the gluon production mechanism is dominant over the medium effects, for the parameter sets used.However, we can not extract further dependencies of the fitting parameters due to numerical limitations.As a result, it is not possible at the moment to further pin down the physical meaning of the fitting parameter values obtained.We leave further research on this topic for future work.

IV. CONCLUSION AND OUTLOOK
In this work, we have implemented a digital quantum circuit to quantum simulate the evolution of a QCD jet.We implement the light-front Hamiltonian formalism and perform the real-time quantum simulation in the |q and the |q + |qg Fock spaces.We have studied the total 14 Notice that for the single particle in vacuum, S vN = 0.
(b) Entropy parameter a as a function of (g 2 µ) 2 momentum broadening of the jet, the gluon production, and the von Neumann entropy associated to the quark state.
We find sizable eikonal effects by comparing the total momentum broadening of the jet in the |q + |qg space at finite energy, to that in the |q space at the eikonal limit of p + = ∞.The underlying physics is that the inclusion of gluon radiation, which depends on p + , enlarges the phase space, and the medium interaction also interferes with this process.Furthermore, when studying the energy distribution of the gluon inside the jet, we recover the reduced vacuum splitting function.When the medium is included, we can not observe significant modifications to the vacuum kernel, but note that there is a larger amount of gluon radiation being produced.We leave the interesting study on the closely related QCD LPM effect for future work when simulations on larger lattices are feasible.Finally, we compute the entropy associated with the quark state, in both the |q and the |q + |qg Fock spaces.In the former case, we recover the classical result for momentum diffusion, in which the entropy grows logarithmically in time.In the latter case, the entropy growth accelerates significantly, mainly due to the production of gluon radiation, irrespective of the medium being present.
Extensions of the algorithm presented here to include higher Fock sectors, such as |qgg , are underway.We note that including this sector would allow to perform numerical calculations beyond known analytical results, see e.g.[9,53].Another interesting avenue to be explored regards to the transition of the prepared final state partonic jet into a hadronic state, such as a pion state obtained on the circuit [54].
Table I.Resource cost of evaluating various Pauli terms with both the orthogonal matrix projection (omp) and the sparse matrix projection (smp) methods.We compare the computational time costs (in seconds) needed for the momentum-space (t P ) and the mixed-space strategies (t M ).The respective sparsities (S) of the medium interaction matrix VA in each method are also provided.Numerical benchmark results are performed on the same Ubuntu 22.04.2LTS machine using 1 CPU core with 32.0 GB memory and an Intel i9 processor of 3.50 GHz.Ideally, both the omp and smp methods can be parallelized using multi-core CPUs.We use N η = 4 for the results presented in the main body of this paper.Here, we show the results of p 2 ⊥ at various N η in the eikonal limit of p + = ∞ limit In Fig. 12(a), we show the transverse momenta at increasing Q 2 s and notice that our results agree with analytical results even at N η = 1.

N
We also examine the convergence of the trotterization step size N t .In particular, we compare the in-medium momentum broadening using the two simulation strategies (momentum-space vs mixed-space) in Fig. 12(b) for a finite initial quark energy, i.e., p + = 1 GeV.We can see that any N t could give reasonable results.At the value of N t used in this work N t = 16, the simulation result is within 1% of the expected value.

Figure 1 .
Figure 1.An illustration of the jet (blue line, dressed by helical lines representing the gluon in the |qg state) evolution in the presence of a highly boosted background medium (orange band) described by a classical field A µ (x).

Nq i e qi =
Nq j e qj = n and Ng k e g k = m.
|p = |p x , p y is represented by a lattice coordinate |k = |k x , k y with p = k b ⊥ and b ⊥ = π/L ⊥ .Similarly, we map any position vector |x = |x x , x y to the lattice position vector |n = |n x , n y with x = n a ⊥ .

Figure 2 .
Figure 2. Schematic representation of the digital circuit used to simulate the multi-parton jet in the medium.(a) The top panel is a quantum circuit of the whole simulation process with registers for the quark, the gluon, and the occupancy status.The bottom panel illustrates that the jet is a superposition of all possible quantum states in the phase space, and the medium shown in the yellow band expands through the process.The single timestep evolution block in the circuit corresponds to a time slice of δx + in the whole process, and τ is the duration of each medium layer.(b) Two treatments of the evolution operator in a single timestep: i. direct exponentiation, and ii.alternating exponentiation in a momentum-position-mixed space.See more discussions in the text.

( 21 )
In analogy, the p 2 ⊥ (x + ) for a gluon state replaces C F = (N 2 c − 1)/(2Nc) by C A = Nc in the above equation.For an uncorrelated quark-gluon state, one should replace C F by C F + C A in the above equation.

Figure 3 . 2 ⊥
Figure 3.The dependence of the momentum broadening p 2 ⊥on the saturation scale Qs, for various The Hamiltonian contains only the medium interaction, as specified in the legends.The eikonal analytical results in the dashed and dotted lines are given by Eq. (21).

Figure 4 . 2 ⊥
Figure 4.The dependence of the momentum broadening 2 ⊥ on the saturation scale Qs, at (a) p + = 1 GeV and (b) p + = 1000 GeV.The initial state is a bare quark with p ⊥ = 0 ⊥ .The results of the simulation in the |q Fock space is in the open triangle, and that in the |q + |qg Fock space is in the disk.The eikonal analytical results in the dashed and dotted lines are given by Eq. (21).

P 10 Figure 5 .
Figure 5. Probability of the |qg sector in vacuum as a function of the evolution time x + with coupling strength (a) g = 1 and (b) g = 10 calculated using quantum simulators.The Hamiltonian contains the gluon emission term, with and without the kinetic energy term, as specified in the legends.Each curve is obtained from both quantum simulations and classical diagonalization, and the two are in agreement.

Figure 6 .
Figure 6.Evolution of the probabilities of different p + states, including the |q sector and the different segments of the |qg sector characterized by the gluon longitudinal momentum fraction zg.

g 2 µ=0. 1 2 Figure 8 .
Figure 8. Probability of the |qg component as a function of the evolution time x + with medium strength g 2 µ = 0.1, 0.2 GeV 3/2 .The result in vacuum is in the solid black line for comparison.

Figure 10 .
Figure 10.Entropy growth of the quark state in the |q Fock space with initial energy p + = 1 GeV.(a) Time evolution of the von Neuman entropy SvN at various lattice sizes.(b) Entropy parameter a as a function of (g 2 µ) 2 at various lattice sizes.

Figure 11 .
Figure 11.Entropy of the quark component SvN in the |q + |qg Fock space, in vacuum and in medium, with initial energy p + = 1 GeV.The result in the |q Fock space is plotted for comparison.Medium strength is fixed as g 2 µ = 0.1 GeV 3/2 .The data point at x + = 10 fm is affected by the lattice boundary and therefore not shown here.

Figure 12 .
Figure 12.Convergence studies of the number of layers (Nη) used in the MV model and time steps (Nt) in the trotterizations.
Using the omp method at N ⊥ = 2 takes an extremely long time, so they are not presented.b We used a fraction of the Pauli terms to estimate the total time because the full calculations are very expensive.Their true time costs are expected to be much larger.