Accelerated quantum circuit Monte-Carlo simulation for heavy quark thermalization

Thermalization of heavy quarks in the quark-gluon plasma (QGP) is one of the most promising phenomena for understanding the strong interaction. The energy loss and momentum broadening at low momentum can be well described by a stochastic process with drag and diffusion terms. Recent advances in quantum computing, in particular quantum amplitude estimation (QAE), promise to provide a quadratic speed-up in simulating stochastic processes. We introduce and formalize an accelerated quantum circuit Monte-Carlo (aQCMC) framework to simulate heavy quark thermalization. With simplified drag and diffusion coefficients connected by Einstein's relation, we simulate the thermalization of a heavy quark in isotropic and anisotropic mediums using an ideal quantum simulator and compare that to thermal expectations. With Grover-like QAE, we calculate physical observables with quadratically fewer resources, which is a boost over the classical MC simulation that usually requires a large sampling number at the same estimation accuracy.


I. INTRODUCTION
Thermalization is one of the most important common features of a non-equilibrium system.An open system that undergoes quantum decoherence by rapidly exchanging information with the environment usually tends to thermalize conventionally and classically.Heavy quark thermalization in the background of quark-gluon plasma (QGP) produced in relativistic heavy-ion collisions (HICs) is such an open system that heavy quarks have distinguished separation of scales compared to the soft QGP medium.With the shear viscosity over entropy density ratio η/s characterizing the speed of hydrodynamization and temperature T as the only scale, the thermalization time of the QGP is characterized by a scale of τ h ≃ 4πη/(T s) [1].Compared to the soft QGP, the relaxation of the heavy quark is prolonged by its heavy mass τ R ≃ m HQ τ h /T .With a typical charm quark mass of ≃1.5 GeV, and the QGP medium temperature of 300-500 GeV in HICs, heavy quark undergoes a time of thermalization caused and dominated by a thermal environment.This is more extreme for the bottom quark with ≃ 4.5 GeV mass, that the thermalization process is not even finished in a 10 fm of the QGP phase.Eventually, the hadronized heavy flavors measured in the detector are not thermalized and the characterization of the heavy quark spectra tells us the medium property of the QGP.
The heavy masses of heavy quarks not only delay the thermalization in the QGP medium but also make the heavy quarks less relativistic compared to the almost massless partons in the QGP.This leads to a wellestablished thermalization description for heavy quarks based on a stochastic process with low-momentum ran-dom kicks from the medium [2][3][4][5][6][7][8][9][10][11].The thermalization in this description is controlled by two competing effects, the energy loss from a drag term and a diffusion from a stochastic term.The energy loss tends to reduce the momentum of a heavy quark while the diffusion tends to broaden the momentum distribution.The competing contributions eventually thermalize the heavy quark to a certain distribution controlled by a fluctuation-dissipation theorem, known as Einstin's relation.In a non-relativistic or static limit, the thermal distribution is given by a classical Maxwell-Boltzmann distribution.For more discussions on heavy quark thermalization in HIC phenomenology, see reviews [12][13][14].Notably, this stochastic process is so generic that it is not limited to the description of a heavy quark thermalization but is broadly utilized in many research topics, such as the Black-Scholes model in quantitative finance, which in part inspired our work.
Quantum computing technology, using laws of quantum mechanics, has already been extensively applied in many areas of nuclear physics [15][16][17][18][19][20][21][22][23][24][25][26][27][28], where the strength of quantum computing is usually exploited from its exponential state space, local Hamiltonian simulation, and near-term variational algorithms.Recently, novel gatebased quantum finance strategy [29][30][31][32] with the quantum amplitude estimation (QAE) [33] exhibits a promising quadratic speed-up over the classical Monte-Carlo (MC) method.In much the same spirit as Grover's algorithm [34,35], the QAE allows efficient estimation of the amplitude of the designated quantum state.The main contribution of this work is the first application of an accelerated quantum circuit Monte-Carlo (aQCMC) strategy using the QAE techniques for heavy-quark thermalization.Different events are simulated as quantum state evolution with sufficient quantum shots, and the physical observables are efficiently extracted with amplitude estimation.With the constant improvements in the QAE algorithms [36,37], the aQCMC may expect to become a more standard approach, especially in future large-scale quantum simulations.
This manuscript is organized as follows.In Sec.II, we review the heavy-quark thermalization formulated as a stochastic differential equation and its standard classical simulation strategy with the MC method.In Sec.III, we discuss the aQCMC strategy utilized in this work to speed up the computation.In Sec.IV, we present our simulation results in isotropic and anisotropic mediums using Qiskit.In Sec.V, we summarize and discuss future avenues of this work.

A. Stochastic description of heavy quark thermalization
The heavy-quark thermalization can be characterized by a stochastic differential equation (SDE) known as the Langevin equation [2], evolving in full positionmomentum phase space ⃗ x, ⃗ p with evolution time t where the random force that sampled as a Wiener process dW ∼ N (0, dt) has correlation ⟨dW i dW j ⟩ = δ ij dt.
The drag coefficient A(⃗ x, ⃗ p, t) and the diffusion coefficient σ ij (⃗ x, ⃗ p, t) in HICs may be calculated from either quantum chromodynamics (QCD) [38][39][40][41][42][43][44][45][46], or QCD-like theories [47][48][49][50][51][52][53][54] with a heavy quark interacting with the medium.Applying Ito's lemma, the Langevin equation Eq. ( 1) can be reformulated as a Kolmogorov-forward equation, known as the Fokker-Planck equation, presenting the time evolution of the heavy quark non-equilibrium distribution f (⃗ x, ⃗ p, t) as with the diffusion coefficient There is no general solution to the Fokker-Planck equation Eq. ( 2) and the evolution would depend on the initial condition.However, the solution to the Fokker-Planck equation would be an attractor towards the thermal limit.This transition from various ordered initial conditions to a unique chaotic limit is the thermalization of heavy quarks within a medium.These transport coefficients are generally medium profile dependent, but in a thermal and homogeneous medium, we may drop the spatial ⃗ x and time t dependencies.The perturbative QCD calculation suggests the drag coefficient A(⃗ p) to be almost a constant at low momentum p < ∼ 2M [44].With an approximately constant drag coefficient, one may simplify the Langevin equation in the nonrelativistic limit at a small momentum p, which may be further rescaled by the heavy quark mass M .Keeping diagonal terms only in the diffusion term, these simplifications lead to a dimensionless Langevin equation In the above equation we have used dimensionless momentum q i = p i /M , time d t = Adt, and anisotropic stochastic terms d Wi ∼ N (0, 2T d t/(M χ 2 i )) with proper Einstein's relation A = σ 2 ii χ 2 i /(2M T ).See App.A for details of the derivations, and also Refs [44,[53][54][55][56][57] for recent discussions on hard probes in anisotropic medium.Notice that the heavy quark relaxation time τ R ≃ 1/A, the value of d t represents the speed of energy loss and thermalization.Thus, a realistic simulation would favor the d t to be as small as possible, and a value of d t ≃ 1/N t takes about N t steps to thermalize (thermalization will also be delayed by a large momentum, for instance, a heavy quark jet).Another relevant scale is the temperature over heavy quark mass ratio T /M in the variance σ2 i d t = 2T d t/(M χ 2 i ).The dimensionless Fokker-Planck equation corresponding to Eq. (3) reads The thermal distribution in terms of these dimensionless quantities reads This stochastic process is usually simulated with the MC methods, by sampling the Wiener process for each time step.The trajectory of a heavy quark contributes to an event and a collection of these events provides a time series of the heavy quark distribution towards thermalization.On a modern digital computer, this MC simulation is straightforward: one starts with whatever heavy quark initial distribution f (⃗ q, t0 ), and samples the heavy quark initial momentum (q t0 x , q t0 y , q t0 z ) accordingly.Similarly, the values of the stochastic variables (dW t x , dW t y , dW t z ) can be uncorrelatedly sampled with a set of independent normal distributions {N (0, 2T d t/(M χ 2 i ))} with i = x, y, z for each time in a diagonalized form.The increment of the momentum follows the Langevin equation Eq. ( 3) and the momentum at the next step can be calculated with the forward-Euler method as Iterating the above algorithm for large enough N t steps from t0 to t0 +N t d t gives a time series of one heavy quark momentum and repeating for a total of N event events produces an emergent phenomenon of heavy quark thermalization, which leads to a thermal distribution with N t d t ≫ 1.Then, for any physical quantity F (⃗ q) at time t, its expectation value would be The MC simulation on a modern computer is straightforward but often requires large computational resources for reasonable precision.By encoding the stochastic process on the quantum circuit and accelerating with the QAE algorithms, one may reduce the inherent problem complexity faced in classical simulations, and obtain a quadratic quantum speed-up compared to the classical method to the same precision.

III. QUANTUM STRATEGY
In this section, we formulate the quantum strategy, the quantum circuit Monte-Carlo (QCMC) to simulate the heavy quark thermalization in a stochastic description.For the QCMC simulation, we encode the particle's momenta q i in each direction as a quantum state.With a generic n-qubit quantum register, one has in principle N = 2 n possible modes for the heavy quark momenta q.By restricting the momentum q ∈ [−q max , q max ), we discretize q into N values with δq = 2q max /N .Then, we further shift the physical momentum q to the positive momentum q by a constant q max so that q = q + q max ∈ [0, 2q max ) and impose a periodic boundary condition, i.e. q = q mod (2q max ).The use of non-negative dimensionless momenta q makes a straightforward binary mapping onto the corresponding quantum states, which can be extended to all three spatial dimensions x, y, and z.
To thermalize with approximately N t steps, reasonable values of the coefficients for the simulation scale as d t ≃ 1/N t with N t > 1.The variance in the stochastic term is chosen to be σ2 i ) according to scales of heavy quark mass M and temperature T in HICs.Since generic quantum multiplication and divisions are complicated [58], we pick d t = 1/2 d with positive integer d in practical simulations, which can be realized on the quantum circuit by shifting the quantum state with d qubits using a sequence of CX gates.
In general, for each of the i = x, y, z directions, we prepare quantum register S i to encode the particle's momentum q i and quantum register W i to encode the diffusion term d Wi .Each register is represented by a set of qubits.The numbers of qubits n S , n W in registers S i , W i are not necessarily the same.The increment at each time step t = nd t contributed from the drag term −q i d t and the diffusion term −d Wi are implemented as unitary quantum operators U A n i following Eq.( 6) so that Here, U A n i = U A is time-independent with a constant drag coefficient A, though it is not required.Now we introduce the quantum gates used in the circuit: 1. Distribution loading gates (U L ) are responsible for loading the initial momentum distribution on the quantum register S for the system.In principle, one can start with either a single momentum or any momentum distribution for the heavy quark and evolve it on the circuit.Here, we initialize an arbitrary single momentum each time using X gates.
3. Quantum evolution gates (U A ) are the main building blocks of the QCMC, where we follow Eq. ( 10) to construct the evolution gates.Specifically, we implement and utilize the quantum adders and multipliers (see App. B for a brief review) to build the stochastic Langevin evolution at each time step.One additional constant quantum adder is included to remedy the momenta from q to nonnegative q per each step.Notably, these quantum arithmetic gates correspond directly to the classical arithmetic operations, though one still needs to manually manipulate these operations at the quantum-register level for today's quantum computers.Since quantum Fourier transforms are innate to most arithmetic operations, it may be more efficient to use Fourier basis as the encoding basis to abbreviate consecutive operations.
In principle, one could simulate the MC process on the quantum circuits as efficiently as on a classical computer.Nevertheless, since at each time step, the Wiener process dW needs to be uncorrelated and the quantum arithmetic operations are on the register level, the quantum circuit would require additional sets of registers S and W for each time iteration, making the total qubit number scales as O((2N t − 1)n) assuming n Q = n W = n.To circumvent this tower-like quantum circuit, one may include reset gates to economically reuse the quantum registers repeatedly for different time steps, as in Fig. 1(a), leading to only O(3n) qubits.
The quantum circuit Monte-Carlo (QCMC) method can be accelerated by taking advantage of the quantum amplitude estimation [33] (QAE), a generalized version of Grover's search algorithm [34,35].Suppose an operator A F acts on n + 1 qubits, where , and a ∈ [0, 1] is the desired expectation of interests.Specifically, a = ⟨ψ| n F |ψ⟩ n is the expectation of any physical quantity F in terms of heavy quark momentum, whose distribution is represented by the quantum state |ψ⟩ n .Using Grover operator Q = A F S 0 A † F S ψ1 with S x signflipping operator on the state x, the QAE allows for highprobability estimation of a in N q queries of A F with error ϵ = O(1/N q ).This represents a quadratic speed-up over classical MC [33].Intuitively, one can consider Grover's operator Q ∼ R ψ * R ψ * 0 as sequence of two reflection operators R ψ * , R ψ * 0 .This operator successively reflects about the "bad" state |ψ * 0 ⟩ n+1 and the "mean" state |ψ * ⟩ n+1 , such that the amplitude of the "good" state |ψ * 1 ⟩ n+1 is amplified.See App.C for more technical details.In principle, one may use the standard quantum phases estimation (QPE) with extra auxiliary qubits [33,59] to retrieve the amplitude where the estimation success rate is quickly boosted close to unity.
In practice, the QAE approach is usually difficult for two reasons: Firstly, universal oracle implementation for the expectation function F is nontrivial; sec- ondly, the QPE, the key to extract amplitude, requires expensive auxiliary qubits and substantial multi-qubit gates [33].Fortunately, operators U F involving piecewise linear functions can be approximated via Taylor expansion and implemented using controlled R Y gates [30,60], so we are capable of investigating momentum and absolute momentum expectation of the particle, i.e., F (q) = q and F (q) = |q|.Alternative loading methods to reduce the circuit complexity that one may consider include quantum generative adversarial networks [61] and approximate quantum compiling [62].
On the other hand, the complexity of the QPE can be circumvented using novel QPE-free algorithms [36,37,[63][64][65][66][67], which are mostly based on selected Grover iterations Q k A F to estimate the quantum amplitude efficiently, and the same quadratic speed-up can be obtained [65].In particular, we focus on the Iterative QAE (IQAE) algorithm in our simulation result, which proves most economical in estimation accuracy and confidence level [36] for our simulation resources.Alternatively, one Quantum simulation using the aQCMC with Iterative QAE [36] (in shaded area) compared to QCMC with direct measurement (in lines) for the earlier time steps.Absolute error, IQAE, N q,max sim, DM sim, aQCMC FIG. 4. Quantum advantage using the aQCMC approach with Iterative QAE [36] to estimate physical observable with quadratically less resources.
may also consider using variational QAE with constantdepth circuits [68].Nonetheless, it is crucial to point out that by having Grover operators in the QAE we cannot use the non-unitary reset gates directly, and consequently, we regress to the tower-like quantum circuit in Fig. 1(b) when the QAE is involved.

IV. SIMULATION RESULTS
With the theory of heavy quark thermalization and its quantum strategy described above, we present the numerical simulation results of 1D and 2D heavy quark thermalization, including both isotropic and anisotropic mediums.While one can use real quantum computers to simulate the heavy quark thermalization, these  5. Probability density distribution for quantum simulation with an anisotropic medium at late time t = 3, compared with analytical thermal equilibriums distribution feq(q) in Eq. ( 5).
devices are still limited today by their short coherence time.Therefore, for the purpose of our investigation, we use the QASM simulator provided by Qiskit to mimic ideal quantum devices.Regarding the physical scales we choose in Eq. ( 3), the heavy quark with a mass M = 1.5 GeV in a typical plasma temperature in heavyion collisions T ≃ 300 -500 MeV gives a range of the variance σ2 i d t = 2T d t/(M χ 2 i ) around 2d t/(5χ 2 i ) -2d t/(3χ 2 i ).For simplicity, we may just consider σ2 i d t ≃ d t/(2χ 2 i ).We first study the heavy quark thermalization using the QCMC approach in a one-dimensional medium, following the quantum circuit in Fig. 1(a).In general, we can simulate the stochastic process with any system size and time step; however, in practice, we are limited by the classical simulation resources.Here, for practical purposes, we consider a small system of n S = n W = 4 qubits and q ∈ [−q max , −q max ) = [−2, 2) with momentum resolution δq = 0.25 and time interval d t = 1/2 d = 0.5.Together with registers to store the intermediate quantum states, a total of 12 qubits and 8192 shots are used.We present the numerical simulation results in Fig. 2. The upper panel of Fig. 2 shows a collection of heavy quark momentum q trajectories as an attractor toward thermal expectation.Simulation with two different parameters χ = 1, and χ = 2 are used.Similarly, the lower panel of Fig. 2 shows a collection of heavy quark momentum absolute value |q| trajectories.In addition, we also present the large time momentum distribution at t = 3 from the simulations with χ = 1, 2 compared to thermal 1-D distributions f eq (q) in Eq. ( 5) with σ2 = 1/(2χ 2 ) = 1/2, 1/8.The comparison is shown in Fig. 5.The simulations agree with the expectations with small discrepancies caused by both momentum and time lattice effects.We also observe that the thermalization with a larger anisotropic parameter χ leads to thermal distribution with a narrower col-lection of momentum, characterized by a smaller value of variance σ2 .
To speed up the QCMC, we implement the QAE to directly extract the expectation ⟨q⟩ and ⟨|q|⟩, following the aQCMC circuit in Fig. 1(b).In particular, we focus on Iterative QAE (IQAE) for its economical complexity in the calculation [36], and present the simulation results1 for the early time steps in Fig. 3. Without using the reset gates, the qubit number grows linearly with the number of time steps, so we simulate to a maximum step of 3, taking up a total of 20 qubits, with an additional ancilla qubit for QAE.We can see that both the aQCMC results agree with the QCMC results within the uncertainty band obtained from 10 simulation batches.Specifically, we used IQAE with an estimation accuracy ϵ = 0.01, a confidence level (1−α) = 95%, and N shots = 30 shots per iteration [36].In total, we used N q ≈ 1000±240 shots per each data point for aQCMC2 , which agrees well with the QCMC results using a large number shots (81920 shots).
To show the quantum advantage using the aQCMC approach with the IQAE, we present the estimation accuracy versus the number of oracle queries N q in Fig. 4. The estimation accuracy is evaluated as the absolute error of the measured momentum to analytical expectation ϵ = | ⟨q⟩ − ⟨q⟩ analytical |, where for convenience, we start with q 0 = 0 and thus ⟨q⟩ analytical = 0. We estimate the absolute error in both direct measurement (DM) and the IQAE using the same total number of shots at t = 0.5.The DM is equivalent to the classical strategy of statistical measurement while the IQAE is its quantum counterpart.We observe the estimation accuracy ϵ follow exactly the O(1/N q ), offering a quadratic boost to classical strategy with O(1/ N q ).Notably, our results are within the theoretical upper bound N q,max = (50/ϵ) log[(2/α) log(π/(4ϵ))] [36].The same quadratic quantum speedup is also found for the Maximum Likelihood QAE [37].
The quantum simulation strategy introduced in this work, both the QCMC and the aQCMC, can also be applied to two-dimensional heavy quark systems with both isotropic and anisotropic mediums.Here, we consider an isotropic medium with χ x = χ y = 1 thus σ2 xx = σ2 yy = 1/2, and an anisotropic medium with χ x = 1, χ y = 2 thus σ2 xx = 1/2, σ2 yy = 1/8.Due to the large amount of total qubits required for a complete twodimensional circuit, the simultaneous calculation of the x, y directions in one circuit is not practical with the current hardware.Instead, since the heavy quark dynamics in x, y directions are decoupled with diagonal coefficients, we can simulate x and y directions separately and pair the events randomly to perform the calculation in two dimensions.With uniformly distributed heavy quarks as the initial condition, we present the heavy quark thermalization over time in Fig. 6.We can see that different thermalization patterns reflect accordingly for the different medium properties, isotropic and anisotropic.In both cases, the collections of heavy quarks reach the thermal distributions provided by Eq. ( 5).
In the anisotropic medium, one may further evaluate the buildup of the elliptic flow v 2 , which characterizes the anisotropization due to the medium profile, f (q, cos(ϕ), t) cos(2ϕ)dϕ f (q, cos(ϕ), t)dϕ where the v 2 in the thermal equilibrium is a ratio of modified Bessel functions I 1 (x) and I 0 (x).In Fig. 7, we calculate v 2 using the simulation result.Despite the discrepancy between the simulated v 2 at a late stage compared to the analytical thermal limit due to insufficient lattice, we observe a gradual buildup of the v 2 for heavy quarks in an anisotropic medium that approaches the limit.

V. CONCLUSION AND OUTLOOK
In this work, for the first time, we present a quantum strategy for stochastically simulated heavy quark thermalization with the QCMC and the aQCMC algorithms on the circuit.Specifically, we simulate the heavy quark thermalization with both the QCMC with a longer evolution time step and the aQCMC with a shorter evolution time step but boosted with amplitude estimation.With these algorithms, we study heavy quark thermalization in both one-dimensional and two-dimensional mediums, as well as isotropic and anisotropic mediums.We show their thermalization patterns and late-time behaviors compared to the analytical expectations.We also calculate the buildup of the elliptic flow for heavy quarks in an anisotropic medium.Remarkably, with Grover-like quantum amplitude estimation, we can estimate physical observables with quadratically a smaller number of simulation shots, compared to classical Monte Carlo methods, which prepares for quantum advantage in future faulttolerant simulation.
Notably, the quantum strategies utilized in the work, the QCMC and aQCMC algorithms, are generic for simulating stochastic processes in even broader contents in physics, where the aQCMC has the potential to speed up the calculation in comparison to classical methods.For instance, with proper modeling of the quark coalescence, this framework can be extended to study quarkonium dissociation and recombination in a medium.Furthermore, real-time dynamics of heavy quark/quarkonium simultaneous production involving many heavy quark/antiquark pairs will benefit from the quantum boost that significantly reduces the number of samples in such stochastic simulation.The classical Wiener process may be replaced by a quantum random walk, which eventually gives quantum statistics instead of classical statistics, and a stochastically quantum thermalization of a heavy quark, or spin-chain system might be achieved on a quantum circuit.We leave these for future studies.
acting on the state and an ancilla qubit, To estimate this a, Grover's operator is defined as )) [36].
The amplitude can also be estimated by applying Q k A F operations, without the need for ancilla qubits.Let a = sin 2 (θ a ) and we observe that  specific choice of k varies in each algorithm [36,37,[63][64][65], and we provide a general circuit in Fig. 9(b) for those QPE-free approaches.
FIG. 1.The schematic quantum circuits of the QCMC and the aQCMC involving Nt time steps in direction i constructed using S and W registers (a) with reset gates, and (b) without reset gates but with the QAE.Quantum gates UL, UW , and UA are used to represent the initial distribution loading, stochastic diffusion, and quantum evolution, respectively.

FIG. 2 .
FIG.2.Quantum simulation using the QCMC method in a one-dimensional system.Panels (a) and (b) show the expectations of ⟨q⟩ and ⟨|q|⟩ for differential initial momenta (marked in colored lines).
FIG. 3.Quantum simulation using the aQCMC with Iterative QAE[36] (in shaded area) compared to QCMC with direct measurement (in lines) for the earlier time steps.