∎

<EMAIL_ADDRESS><EMAIL_ADDRESS>
11institutetext:  Instituto Galego de Física de Altas Enerxías (IGFAE),
Universidade de Santiago de Compostela, E-15782 Galicia, Spain

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

João Baratae2,addr1  Carlos A. Salgadoe3,addr1

(Received: date / Accepted: date)

###### Abstract

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.

###### Keywords:

Jet quenching Quantum Simulation

††journal: Eur. Phys. J. C

##  1 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 Feynman:1981tf
. Since then quantum simulation as seen widespread application in physics,
chemistry and many other areas Georgescu:2013oza ; Schuld_2014 ; Or_s_2019 .

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

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

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

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

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

##  2 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 ${\cal A}$. We will
consider the cases where the propagating parton, a quark, is in the singlet or
fundamental color representations, though we will mostly assume that the
evolution in color space is trivial. Detailed and ample discussions on jet
quenching theory can be found, for example, in Blaizot:2012fh ;
CasalderreySolana:2007zz ; Blaizot:2015lma ; Mehtar-Tani:2013pia  and
references therein.

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

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

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

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

$\left(i\partial_{t}+\frac{\partial^{2}_{\boldsymbol{x}}}{2\omega}+g{\cal
A}^{-}(t,{\boldsymbol{x}})\cdot
T\right)G(t,{\boldsymbol{x}};0,{\boldsymbol{y}})=i\delta(t)\delta({\boldsymbol{x}}-{\boldsymbol{y}})\,,$
(1)

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

${\cal H}(t)=\frac{{\boldsymbol{p}}^{2}}{2\omega}+g{\cal
A}^{-}(t,{\boldsymbol{x}})\cdot T={\cal H}_{K}+{\cal H}_{\cal A}(t)\,,$ (2)

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

##  3 Quantum simulating momentum broadening

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

${\cal U}(t,0)\equiv{\cal T}\exp\left[-i\int_{0}^{t}ds\,{\cal H}(s)\right]\,,$
(3)

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

$|\psi_{t}\rangle={\cal U}(t,0)|\psi_{0}\rangle\,.$ (4)

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

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

  1. 1.

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

  2. 2.

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

  3. 3.

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

  4. 4.

Time evolution – Implementation of Eq. 4;

  5. 5.

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

  6. 6.

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

Figure 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 $\blacksquare$ denotes that the time evolution gates parameters
are to be determined from the field $A$.

###  3.1 Encoding

We begin by discretizing the problem in position space, such that
$\ket{{\boldsymbol{x}}}=\ket{a_{s}{\boldsymbol{n}}}$, with $a_{s}$ the spatial
lattice spacing and ${\boldsymbol{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 ${\boldsymbol{x}}_{\rm
max}=a_{s}(N_{s}-1,N_{s}-1)$. Also, the spatial discretization induces a
lattice discretization in momentum space with
$\ket{{\boldsymbol{p}}}=\ket{a_{d}{\boldsymbol{q}}}$ and
$a_{d}=\frac{2\pi}{a_{s}N_{s}}$ the momentum space lattice spacing, with
${\boldsymbol{q}}=(q_{1},q_{2})$ a two dimensional vector with each component
also taking integer values between $0$ and $N_{s}-1$222In the following
discussion we will consider only positive values for the position and momentum
of the quark, see A..

One can rewrite the Hamiltonian ${\cal H}$ in terms of the dimensionless
Hamiltonian $H={\cal H}a_{s}$ (see A)

$H=\frac{{\boldsymbol{P}}^{2}}{2E}+gA(t,{\boldsymbol{X}})\cdot
T=H_{K}+H_{A}(t)\,,$ (5)

with
$\hat{{\boldsymbol{P}}}\ket{{\boldsymbol{q}}}={\boldsymbol{q}}\ket{{\boldsymbol{q}}}$
and
$\hat{{\boldsymbol{X}}}\ket{{\boldsymbol{n}}}={\boldsymbol{n}}\ket{{\boldsymbol{n}}}$
the dimensionless position and momentum operators. Also
$A(t,{\boldsymbol{n}})\cdot T=a_{s}{\cal A}^{-}(t,a_{s}{\boldsymbol{n}})\cdot
T$ and $E=\frac{N_{s}^{2}\omega a_{s}}{4\pi^{2}}$ 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 NielsenChuang  to
denote the single up spin state $\ket{\uparrow}=\ket{0}=[1,0]^{T}$ in the
computational basis (with the last equality giving the associated vector
representation) and $\ket{\downarrow}=\ket{1}=[0,1]^{T}$. Then, any component
of the vector $\ket{{\boldsymbol{n}}}$ can be represented by a product of many
spins, in a binary basis (see A). The associated momentum state vector
$\ket{{\boldsymbol{q}}}$ is obtained by applying a standard quantum Fourier
Transform (qFT).

###  3.2 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
$\ket{\psi_{0}}$ from the fiducial state $\ket{0}^{\otimes 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
parameter $\hat{q}$ from the quantum simulation output, we wish to avoid
contributions coming from initial state physics. Therefore, we shall mainly
focus on the case where
$\ket{\psi_{0}}=\ket{{\boldsymbol{p}}={\boldsymbol{0}}}$.

However, including a localized initial state distribution might be important
for certain digitizations where one can not prepare the state
$\ket{{\boldsymbol{p}}={\boldsymbol{0}}}$ exactly or if one is simply
interested in studying how different production mechanisms influence the final
state. Several strategies to prepare $\ket{\psi_{0}}$ from an integrable
template distributions can be found in the literature Barata:2020jtq ;
kitaev2009wavefunction ; grover2002creating ; kaye2004quantum . Depending
exactly on what $\ket{\psi_{0}}$ one wants to prepare, in principle, one can
devise a routine which only requires $\mathcal{O}(n_{Q})$ basic quantum gate
operations.

###  3.3 Time evolution

After the initial state $\ket{\psi_{0}}$ has been prepared, we time evolve it
for a time $L$, producing the final state $\ket{\psi_{L}}$. The time evolution
operator in Eq. 3 can be written in terms of the dimensionless Hamiltonian $H$
and medium length $L^{\prime}\equiv L/a_{s}$333In general, one could choose
another length scale to make time dimensionless, leading to the appearance of
a ratio between $a_{s}$ and such scale..

$U(L^{\prime},0)\equiv{\cal
T}\exp\left[-i\int_{0}^{L^{\prime}}dt\,H(t)\right]\,.$ (6)

Directly implementing $U(L^{\prime},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 Poulin_2011 , decomposing $U$ as

$U(L^{\prime},0)\approx\prod_{k_{t}=1}^{N_{t}}\left\\{\exp\left[-iH_{K}\frac{L^{\prime}}{N_{t}}\right]\exp\left[-iH_{A}\left(k_{t}\cdot\frac{L^{\prime}}{N_{t}}\right)\frac{L^{\prime}}{N_{t}}\right]\right\\}\equiv\prod_{k_{t}=1}^{N_{t}}\left\\{U_{K}(\varepsilon_{t})U_{A}(k_{t}\cdot\varepsilon_{t},\varepsilon_{t})\right\\}\,,$
(7)

where we have effectively sliced time into $N_{t}$ steps, each with a length
$\varepsilon_{t}\equiv L^{\prime}/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}\cdot\varepsilon_{t},(k_{t}+1)\cdot\varepsilon_{t})$ the field $A$ is
taken to be constant, leading to the constraint
$\varepsilon_{t}^{-1}\gg||\partial_{t}H_{A}(t)||$; there exist algorithms
Poulin_2011  which circumvent this constraint, as well as other strategies
(see for example Berry_2015 ; Berry_2020 ; Wiebe_2010 ; Berry_2014 ) 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 $\mathcal{O}(\varepsilon_{t}^{2})$.

Let us now consider the $k_{t}^{\rm th}$ 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 $\ket{{\boldsymbol{p}}}\to\ket{{\boldsymbol{x}}}$, time evolve
with $H_{A}$ and transform back to the $\ket{{\boldsymbol{p}}}$ basis, the
generated state being the input to the $k_{t}+1^{\rm th}$ time slice; this
strategy is illustrated in Fig. 2.

Figure 2: Outline of the implementation of the time evolution operator $U$.
Here we detail the $k_{t}^{\rm th}$ time step, as indicated in Eq. 7.

The time evolution operator $U_{K}$ is diagonal in the
$\ket{{\boldsymbol{p}}}$ basis

$\begin{split}U_{K}(\varepsilon_{t})\ket{{\boldsymbol{p}}}=\exp\left(-i\frac{\varepsilon_{t}}{2E}{\boldsymbol{p}}^{2}\right)\ket{{\boldsymbol{p}}}\,,\end{split}$
(8)

thus one only needs to implement a circuit which generates a state dependent
phase. This can be achieved using the algorithm introduced in Zalka_1998 ,
which we detail in B.

After performing the qFT, using standard implementations of the circuit
NielsenChuang , one has to compute the action of $U_{A}$. Although this
operator is diagonal in the $\ket{{\boldsymbol{x}}}$ basis,

$\begin{split}U_{A}(k_{t}\cdot\varepsilon_{t},\varepsilon_{t})\ket{{\boldsymbol{x}}}=\exp(-ig\varepsilon_{t}A(k_{t}\cdot\varepsilon_{t},{\boldsymbol{x}}))\ket{{\boldsymbol{x}}}\,,\end{split}$
(9)

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
B). However, this assumes that one could construct $N_{t}$ oracles which
quantum compute $A(k_{t}\cdot\varepsilon_{t},{\boldsymbol{x}})$ for every
${\boldsymbol{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
$\mathcal{O}(N_{t}\times N_{s}^{2})$ 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 B, that requires solving a system of linear equations with
$N_{s}^{2}$ independent variables, which following the same arguments as above
should be doable in practice, at least for near term small system
applications444We note however that this linear system only has to be solved
once for each $n_{Q}$.. 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 Heeres_2015 555Quantum
strategies to simulate the time evolution of the background field could also
be coupled to our strategy. This could in principle simplify the
implementation of $U_{A}$..

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.

###  3.4 Measurement

Having prepared the state
$\ket{\psi_{L}}=\sum_{\boldsymbol{q}}\psi_{L}^{\boldsymbol{q}}\ket{{\boldsymbol{q}}}$
one could simply measure all the $2n_{Q}$ qubits, obtain the probabilities
$|\psi_{L}^{\boldsymbol{q}}|^{2}$ for every ${\boldsymbol{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 ${\boldsymbol{p}}={\boldsymbol{0}}$. In this case the
coefficients $|\psi_{L}^{\boldsymbol{q}}|^{2}$ are directly related to the
single particle broadening distribution; see 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
$\chi$ (to be detailed below). Then, one averages over all $m$ expectation
values

$\langle\chi\rangle_{\rm
M}=\frac{1}{\sum_{i=1}^{m}p_{A^{(i)}}}\sum_{i=1}^{m}p_{A^{(i)}}\langle\chi\rangle_{\rm
QM}^{(i)}\,,$ (10)

where $p_{{\cal A}^{-}}=p_{A}$, the $i$ superscript denotes a particular field
configuration, running up to $m$, and $\langle.\rangle_{\rm M}$ denotes the
average over field configurations while $\langle.\rangle_{\rm 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 McLerran:1993ka ; McLerran:1993ni  and are
encapsulated in the field-field correlator Kovchegov:1996ty ; Kovchegov:1997pc
; McLerran:1993ka ; McLerran:1993ni . 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 MehtarTani:2006xq ,
the physical origin of ${\cal 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 parameter $\hat{q}$, both in time and in
orthogonal spatial directions Ipp:2020nfu , 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 $\hat{q}$ for that ${\cal A}^{-}$. We add an
ancilla qubit to the circuit and perform the Hadamard test detailed in Fig. 3.

Figure 3: Circuit representation of the measurement strategy.

One first transforms the ancilla, which can be either prepared in the state
$\ket{0}$ or in the superposition $1/\sqrt{2}(\ket{0}+i\ket{1})$, by the
Hadamard gate $H=H^{\dagger}$, and then applies a unitary transformation $V$
on the physical state if the ancilla is in the state $\ket{1}$. Finally the
transformation on the ancilla is reversed and one measures the qubit. We
associate the measured value to a random variable $\chi$ which takes the
values $-1$ if we observe the state $\ket{0}$ and $+1$ if the state $\ket{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 $\ket{0}$ (see D),
then

$\langle\chi\rangle_{\rm
QM}\equiv\bra{\psi_{L}}V+V^{\dagger}\ket{\psi_{L}}=\Re\bra{\psi_{L}}V\ket{\psi_{L}}\,.$
(11)

On the other hand, if the ancilla is prepared in the state
$1/\sqrt{2}(\ket{0}+i\ket{1})$, we have that

$\langle\chi\rangle_{\rm QM}=\Im\bra{\psi_{L}}V\ket{\psi_{L}}\,,$ (12)

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_{\alpha}=\exp(i\alpha{\boldsymbol{P}}^{2})$. Then

$\Re\bra{\psi_{L}}V_{\alpha}\ket{\psi_{L}}=\langle\cos(\alpha{\boldsymbol{P}}^{2})\rangle_{\rm
QM}\,,$ (13)

and

$\Im\bra{\psi_{L}}V_{\alpha}\ket{\psi_{L}}=\langle\sin(\alpha{\boldsymbol{P}}^{2})\rangle_{\rm
QM}\,,$ (14)

from which one extracts $\langle e^{i\alpha{\boldsymbol{P}}^{2}}\rangle_{\rm
QM}$, by definition. We also have that

$\begin{split}\langle e^{i\alpha{\boldsymbol{P}}^{2}}\rangle_{\rm
QM}&=1+\sum_{k=1}^{\infty}\frac{i\alpha^{k}}{k!}\langle\langle
2k\rangle\rangle\,,\end{split}$ (15)

where $\langle\langle
2k\rangle\rangle\equiv\langle{\boldsymbol{P}}^{2k}\rangle_{\rm 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_{d}^{2}\langle\langle 2\rangle\rangle=\hat{q}L$, where we inserted
$a_{d}^{2}$ to get the correct dimensions.

Furthermore, one has the freedom to vary $\alpha$ such that, for small enough
$\alpha$, only linear variations are relevant

$\langle e^{i\alpha{\boldsymbol{P}}^{2}}\rangle_{\rm QM}\approx
1+i\frac{\alpha}{a_{d}^{2}}\hat{q}L\to\langle\sin(\alpha{\boldsymbol{P}}^{2})\rangle_{\rm
QM}\approx\frac{\alpha}{a_{d}^{2}}\hat{q}L\,.$ (16)

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 $\alpha$ corrections, then one has access to the
even moments of the momentum distribution and the respective cumulants. One
can thus imagine varying $\alpha$ and from the observed evolution retrieving
$\langle\langle 2k\rangle\rangle$ 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 $\alpha$ values to be
explored.

If one is only interested in extracting $\hat{q}$ (which is the most relevant
medium parameter for jet quenching), one could consider the unitary
$V=\exp(iF({\boldsymbol{P}}^{2}))$, with
$F({\boldsymbol{P}}^{2})=\arccos({\boldsymbol{P}}^{2})$. Then, for the case
where the ancilla is initially set to $\ket{0}$, we obtain

$\langle X\rangle_{\rm
QM}=\bra{\psi_{L}}\cos(\arccos({\boldsymbol{P}}^{2}))\ket{\psi_{L}}=\langle\langle
2\rangle\rangle\,.$ (17)

In principle, one could implement this protocol following Zalka_1998  (see
also B), provided an efficient arithmetic oracle could be constructed.

##  4 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\cdot
T=A^{a}T^{a}=\frac{1}{2}A^{a}\lambda^{a}$, where $\lambda^{a}$ denotes the
eight Gell-Mann matrices. To deal with this modification, we further split the
time evolution operator to take the form $U=U_{K}\cdot U_{A^{1}}\cdot
U_{A^{2}}\cdots U_{A^{8}}$. Additionally, we must track the color of the quark
as it evolves. To do that, we add a new register with two qubits, which stores
the color state of the quark. In particular we use the following map between
the logical and physical states: $\ket{0,0}\equiv\ket{\rm red}=\ket{R}$,
$\ket{0,1}\equiv\ket{\rm green}=\ket{G}$, $\ket{1,0}\equiv\ket{\rm
blue}=\ket{B}$ and $\ket{1,1}\equiv\ket{W}$, with the latter state not being
physical and therefore absent from any calculation.

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

$\lambda^{1}=\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\to\begin{pmatrix}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{pmatrix}\equiv\tilde{\lambda}^{1}\,,$ (18)

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

$CH=\begin{pmatrix}1/\sqrt{2}&1/\sqrt{2}&0&0\\\ 1/\sqrt{2}&-1/\sqrt{2}&0&0\\\
0&0&1&0\\\ 0&0&0&1\end{pmatrix}\,,$ (19)

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

$e^{-\frac{ig\varepsilon_{t}}{2}A^{1}\otimes\tilde{\lambda}^{1}}=(1\otimes
CH)e^{-\frac{ig\varepsilon_{t}}{2}A^{1}\otimes\tilde{\sigma}^{Z}}(1\otimes
CH)\,.$ (20)

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

$\begin{split}e^{-i\frac{g\varepsilon_{t}}{2}A^{1}\otimes\tilde{\sigma}^{Z}}\ket{{\boldsymbol{x}}}\otimes\ket{c}&=\sum_{n}\frac{(-ig\varepsilon_{t})^{n}}{2^{n}n!}(A^{1}({\boldsymbol{X}})\tilde{\sigma}^{Z})^{n}\ket{{\boldsymbol{x}}}\ket{c}=\ket{{\boldsymbol{x}}}\sum_{n}\frac{(-ig\varepsilon_{t}A^{1}({\boldsymbol{x}}))^{n}}{2^{n}n!}(\tilde{\sigma}^{Z})^{n}\ket{c}\,,\end{split}$
(21)

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

Figure 4: Implementation of the (infinitesimal) time evolution operator
generated by $H_{A^{1}}$.

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
classically777This is not necessary if the qubits storing the color
information are not measured..

##  5 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,
$\langle{\boldsymbol{p}}^{2}\rangle\sim\hat{q}L\equiv Q_{s}^{2}$. We are
interested in length scales $L$ of the order of the nuclear radius of heavy
elements, like Pb or Au, which we take to be $L\sim\mathcal{O}(10\,\rm
fm)=\mathcal{O}(50\,\rm GeV^{-1})$. The value of the jet quenching parameters
$\hat{q}$ varies drastically between different experimental set-ups, due to
the different energy scales being explored. To bridge RHIC, LHC and EIC
experimental conditions, we assume that $\mathcal{O}(0.1\,\rm
GeV^{2}fm^{-1})\leq\hat{q}\leq\mathcal{O}(10\,\rm GeV^{2}fm^{-1})$
Feal:2019xfl ; Liu:2018trl ; Barata:2020jtq . The saturation scale $Q_{s}^{2}$
is then approximately bounded by $Q_{s}^{2}\sim\mathcal{O}(1-100\,{\rm
GeV}^{2})$.

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

$|{\boldsymbol{p}}_{\rm
max.}|\approx\frac{2\pi}{a_{s}}\gg\mathcal{O}(1-10\,{\rm GeV})\,,$ (22)

thus

$a_{s}\ll\mathcal{O}(1-10\,{\rm GeV^{-1}})=\mathcal{O}(0.1-1\,{\rm fm})\,.$
(23)

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

$\mu<a_{d}<Q_{s}\sim\frac{\mu}{Q_{s}}<\frac{1}{N_{s}}<1\,,$ (24)

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

$1<N_{s}<100\iff 0<n_{Q}<7\,.$ (25)

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

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

$1\geq\frac{\lambda}{L}\gg\frac{1}{\mu L}\,.$ (26)

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

$1\leq N_{t}\ll\mu L\implies 1\leq N_{t}\ll\mathcal{O}(100)\,.$ (27)

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

##  6 Conclusions and Outlook

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

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

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

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

###### Acknowledgements.

JB is grateful to Niklas Mueller, Raju Venugopalan and Andrey Tarasov for
helpful discussions on QIS and to Yacine Mehtar-Tani for pointing out quantum
simulating broadening could be relevant. We also are grateful to Nestor
Armesto for useful comments. JB is supported by a fellowship from “la Caixa”
Foundation (ID 100010434). The fellowship code is LCF/BQ/ DI18/11660057. This
work has received financial support from European Union’s Horizon 2020
research and innovation program under the grant agreement No. 82409; from
Xunta de Galicia (Centro singular de investigación de Galicia accreditation
2019-2022); from the European Union ERDF; from the Spanish Research State
Agency by “María de Maeztu” Units of Excellence program MDM-2016-0692 and
project FPA2017-83814-P and from the European Research Council project
ERC-2018-ADG-835105 YoctoLHC

## References

  * (1) R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21 (1982) 467–488.  doi:10.1007/BF02650179.
  * (2) I. M. Georgescu, S. Ashhab, F. Nori, Quantum Simulation, Rev. Mod. Phys. 86 (2014) 153.  arXiv:1308.6253, doi:10.1103/RevModPhys.86.153.
  * (3) M. Schuld, I. Sinayskiy, F. Petruccione, An introduction to quantum machine learning, Contemporary Physics 56 (2) (2014) 172–185.  doi:10.1080/00107514.2014.964942.
URL http://dx.doi.org/10.1080/00107514.2014.964942

  * (4) R. Orús, S. Mugel, E. Lizaso, Quantum computing for finance: Overview and prospects, Reviews in Physics 4 (2019) 100028.  doi:10.1016/j.revip.2019.100028.
URL http://dx.doi.org/10.1016/j.revip.2019.100028

  * (5) J. Carlson, D. J. Dean, M. Hjorth-Jensen, D. Kaplan, J. Preskill, K. Roche, M. J. Savage, M. Troyer, Quantum computing for theoretical nuclear physics, a white paper prepared for the us department of energy, office of science, office of nuclear physics, Tech. rep., USDOE Office of Science (SC)(United States) (2018).
  * (6) I. C. Cloët, M. R. Dietrich, J. Arrington, A. Bazavov, M. Bishof, A. Freese, A. V. Gorshkov, A. Grassellino, K. Hafidi, Z. Jacob, et al., Opportunities for nuclear physics & quantum information science, arXiv preprint arXiv:1903.05453 (2019).
  * (7) K. Matchev, S. Mrenna, P. Shyamsundar, J. Smolinsky, Quantum computing for hep theory and phenomenology, Quantum (2020).
  * (8) K. Bepari, S. Malik, M. Spannowsky, S. Williams, Towards a Quantum Computing Algorithm for Helicity Amplitudes and Parton Showers (10 2020).  arXiv:2010.00046.
  * (9) C. W. Bauer, W. A. de Jong, B. Nachman, D. Provasoli, Quantum Algorithm for High Energy Physics Simulations, Phys. Rev. Lett. 126 (6) (2021) 062001.  arXiv:1904.03196, doi:10.1103/PhysRevLett.126.062001.
  * (10) D. Pires, P. Bargassa, J. a. Seixas, Y. Omar, A Digital Quantum Algorithm for Jet Clustering in High-Energy Physics (1 2021).  arXiv:2101.05618.
  * (11) A. Y. Wei, P. Naik, A. W. Harrow, J. Thaler, Quantum Algorithms for Jet Clustering, Phys. Rev. D 101 (9) (2020) 094015.  arXiv:1908.08949, doi:10.1103/PhysRevD.101.094015.
  * (12) C. W. Bauer, M. Freytsis, B. Nachman, Simulating collider physics on quantum computers using effective field theories (2 2021).  arXiv:2102.05044.
  * (13) W. A. De Jong, M. Metcalf, J. Mulligan, M. Ploskoń, F. Ringer, X. Yao, Quantum simulation of open quantum systems in heavy-ion collisions (10 2020).  arXiv:2010.03571.
  * (14) M. Gyulassy, P. Levai, I. Vitev, Reaction operator approach to multiple elastic scatterings, Phys. Rev. D 66 (2002) 014005.  arXiv:nucl-th/0201078, doi:10.1103/PhysRevD.66.014005.
  * (15) J.-P. Blaizot, F. Dominguez, E. Iancu, Y. Mehtar-Tani, Medium-induced gluon branching, JHEP 01 (2013) 143.  arXiv:1209.4585, doi:10.1007/JHEP01(2013)143.
  * (16) N. Mueller, A. Tarasov, R. Venugopalan, Deeply inelastic scattering structure functions on a hybrid quantum computer, Physical Review D 102 (1) (2020) 016007\.
  * (17) B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, M. Troyer, Hybrid quantum-classical approach to correlated materials, Phys. Rev. X 6 (2016) 031045\.  doi:10.1103/PhysRevX.6.031045.
URL https://link.aps.org/doi/10.1103/PhysRevX.6.031045

  * (18) B. M. Terhal, Quantum error correction for quantum memories, Reviews of Modern Physics 87 (2) (2015) 307–346.  doi:10.1103/revmodphys.87.307.
URL http://dx.doi.org/10.1103/RevModPhys.87.307

  * (19) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, S. Zhu, Theory of trotter error with commutator scaling, Physical Review X 11 (1) (Feb 2021).  doi:10.1103/physrevx.11.011020.
URL http://dx.doi.org/10.1103/PhysRevX.11.011020

  * (20) J. Casalderrey-Solana, C. A. Salgado, Introductory lectures on jet quenching in heavy ion collisions, Acta Phys. Polon. B 38 (2007) 3731–3794.  arXiv:0712.3443.
  * (21) J.-P. Blaizot, Y. Mehtar-Tani, Jet Structure in Heavy Ion Collisions, Int. J. Mod. Phys. E 24 (11) (2015) 1530012.  arXiv:1503.05958, doi:10.1142/S021830131530012X.
  * (22) Y. Mehtar-Tani, J. G. Milhano, K. Tywoniuk, Jet physics in heavy-ion collisions, Int. J. Mod. Phys. A 28 (2013) 1340013.  arXiv:1302.2579, doi:10.1142/S0217751X13400137.
  * (23) J.-P. Blaizot, Y. Mehtar-Tani, The Classical field created in early stages of high energy nucleus-nucleus collisions, Nucl. Phys. A 818 (2009) 97–119.  arXiv:0806.1422, doi:10.1016/j.nuclphysa.2008.11.010.
  * (24) M. A. Nielsen, I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2010.
  * (25) J. Barata, N. Mueller, A. Tarasov, R. Venugopalan, Single-particle digitization strategy for quantum computation of a $\phi^{4}$ scalar field theory (11 2020).  arXiv:2012.00020.
  * (26) A. Kitaev, W. A. Webb, Wavefunction preparation and resampling using a quantum computer (2009).  arXiv:0801.0342.
  * (27) L. Grover, T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, arXiv preprint quant-ph/0208112 (2002).
  * (28) P. Kaye, M. Mosca, Quantum networks for generating arbitrary quantum states, arXiv preprint quant-ph/0407102 (2004).
  * (29) D. Poulin, A. Qarry, R. Somma, F. Verstraete, Quantum simulation of time-dependent hamiltonians and the convenient illusion of hilbert space, Physical Review Letters 106 (17) (Apr 2011).  doi:10.1103/physrevlett.106.170501.
URL http://dx.doi.org/10.1103/PhysRevLett.106.170501

  * (30) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, R. D. Somma, Simulating hamiltonian dynamics with a truncated taylor series, Physical Review Letters 114 (9) (Mar 2015).  doi:10.1103/physrevlett.114.090502.
URL http://dx.doi.org/10.1103/PhysRevLett.114.090502

  * (31) D. W. Berry, A. M. Childs, Y. Su, X. Wang, N. Wiebe, Time-dependent hamiltonian simulation with l1-norm scaling, Quantum 4 (2020) 254.  doi:10.22331/q-2020-04-20-254.
URL http://dx.doi.org/10.22331/q-2020-04-20-254

  * (32) N. Wiebe, D. Berry, P. Høyer, B. C. Sanders, Higher order decompositions of ordered operator exponentials, Journal of Physics A: Mathematical and Theoretical 43 (6) (2010) 065203.  doi:10.1088/1751-8113/43/6/065203.
URL http://dx.doi.org/10.1088/1751-8113/43/6/065203

  * (33) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, R. D. Somma, Exponential improvement in precision for simulating sparse hamiltonians, Proceedings of the forty-sixth annual ACM symposium on Theory of computing (May 2014).  doi:10.1145/2591796.2591854.
URL http://dx.doi.org/10.1145/2591796.2591854

  * (34) C. Zalka, Simulating quantum systems on a quantum computer, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454 (1969) (1998) 313–322.  doi:10.1098/rspa.1998.0162.
URL http://dx.doi.org/10.1098/rspa.1998.0162

  * (35) R. W. Heeres, B. Vlastakis, E. Holland, S. Krastanov, V. V. Albert, L. Frunzio, L. Jiang, R. J. Schoelkopf, Cavity state manipulation using photon-number selective phase gates, Physical Review Letters 115 (13) (Sep 2015).  doi:10.1103/physrevlett.115.137002.
URL http://dx.doi.org/10.1103/PhysRevLett.115.137002

  * (36) L. D. McLerran, R. Venugopalan, Gluon distribution functions for very large nuclei at small transverse momentum, Phys. Rev. D 49 (1994) 3352–3355.  arXiv:hep-ph/9311205, doi:10.1103/PhysRevD.49.3352.
  * (37) L. D. McLerran, R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, Phys. Rev. D 49 (1994) 2233–2241.  arXiv:hep-ph/9309289, doi:10.1103/PhysRevD.49.2233.
  * (38) Y. V. Kovchegov, NonAbelian Weizsacker-Williams field and a two-dimensional effective color charge density for a very large nucleus, Phys. Rev. D 54 (1996) 5463–5469.  arXiv:hep-ph/9605446, doi:10.1103/PhysRevD.54.5463.
  * (39) Y. V. Kovchegov, Quantum structure of the nonAbelian Weizsacker-Williams field for a very large nucleus, Phys. Rev. D 55 (1997) 5445–5455.  arXiv:hep-ph/9701229, doi:10.1103/PhysRevD.55.5445.
  * (40) Y. Mehtar-Tani, Relating the description of gluon production in pA collisions and parton energy loss in AA collisions, Phys. Rev. C 75 (2007) 034908.  arXiv:hep-ph/0606236, doi:10.1103/PhysRevC.75.034908.
  * (41) A. Ipp, D. I. Müller, D. Schuh, Jet momentum broadening in the pre-equilibrium Glasma, Phys. Lett. B 810 (2020) 135810.  arXiv:2009.14206, doi:10.1016/j.physletb.2020.135810.
  * (42) X. Feal, C. A. Salgado, R. A. Vazquez, Jet quenching tests of the QCD Equation of State (11 2019).  arXiv:1911.01309.
  * (43) X. Liu, F. Ringer, W. Vogelsang, F. Yuan, Lepton-jet Correlations in Deep Inelastic Scattering at the Electron-Ion Collider, Phys. Rev. Lett. 122 (19) (2019) 192003.  arXiv:1812.08077, doi:10.1103/PhysRevLett.122.192003.
  * (44) X.-N. Wang, M. Gyulassy, Gluon shadowing and jet quenching in A + A collisions at s**(1/2) = 200-GeV, Phys. Rev. Lett. 68 (1992) 1480–1483.  doi:10.1103/PhysRevLett.68.1480.
  * (45) P. Aurenche, F. Gelis, H. Zaraket, A Simple sum rule for the thermal gluon spectral function and applications, JHEP 05 (2002) 043.  arXiv:hep-ph/0204146, doi:10.1088/1126-6708/2002/05/043.
  * (46) M. Gyulassy, P. Levai, I. Vitev, Reaction operator approach to nonAbelian energy loss, Nucl. Phys. B 594 (2001) 371–419.  arXiv:nucl-th/0006010, doi:10.1016/S0550-3213(00)00652-0.
  * (47) U. A. Wiedemann, Gluon radiation off hard quarks in a nuclear environment: Opacity expansion, Nucl. Phys. B 588 (2000) 303–344.  arXiv:hep-ph/0005129, doi:10.1016/S0550-3213(00)00457-0.
  * (48) V. Vaidya, X. Yao, Transverse momentum broadening of a jet in quark-gluon plasma: an open quantum system EFT, JHEP 10 (2020) 024.  arXiv:2004.11403, doi:10.1007/JHEP10(2020)024.
  * (49) V. Vaidya, Forward scattering in a thermal Plasma (1 2021).  arXiv:2101.02225.
  * (50) J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2 (2018) 79.
  * (51) J. J. Sakurai, J. Napolitano, Modern Quantum Mechanics, Quantum physics, quantum information and quantum computation, Cambridge University Press, 2020\.  doi:10.1017/9781108587280.
  * (52) N. Klco, M. J. Savage, Digitization of scalar fields for quantum computing, Phys. Rev. A 99 (5) (2019) 052335.  arXiv:1808.10378, doi:10.1103/PhysRevA.99.052335.

##  Appendix A Discretization and encoding details

In this appendix we give the details on the discretization of the quantum
mechanical system considered in the main text and the map to the qubits
available in the quantum computer. The discussion in this appendix is quite
standard and can be found in quantum mechanics/computing textbooks
Sakurai:2011zz ; NielsenChuang .

We discretize space using a two dimensional lattice, with lattice spacing
$a_{s}$ and $N_{s}$ lattice sites per dimension, with the spatial cut-off (per
dimension) given by $a_{s}(N_{s}-1)$. We write the position ket
$\ket{{\boldsymbol{x}}}=\ket{a_{s}{\boldsymbol{n}}}$888Notice that
$\ket{{\boldsymbol{x}}}$ has the same mass dimension as
${\boldsymbol{x}}^{-1}$., with ${\boldsymbol{n}}$ a two dimensional integer
vector. Conversely, the momentum induced lattice has a spacing
$a_{d}=\frac{2\pi}{N_{s}a_{s}}$ (with an associated dimensionless integer
vector ${\boldsymbol{q}}$) and the two bases are related by a Fourier
transform

$\begin{split}|{\boldsymbol{p}}\rangle&=\int_{\boldsymbol{x}}e^{-i{\boldsymbol{p}}\cdot{\boldsymbol{x}}}|{\boldsymbol{x}}\rangle\to
a_{s}^{2}\sum_{\boldsymbol{n}}e^{-2\pi
i\frac{{\boldsymbol{q}}\cdot{\boldsymbol{n}}}{N_{s}}}|{\boldsymbol{n}}a_{s}\rangle\,,\end{split}$
(28)
$\begin{split}|{\boldsymbol{x}}\rangle&=\int_{\boldsymbol{p}}e^{i{\boldsymbol{p}}\cdot{\boldsymbol{x}}}|{\boldsymbol{p}}\rangle\to\frac{a_{d}^{2}}{(2\pi)^{2}}\sum_{\boldsymbol{q}}e^{2\pi
i\frac{{\boldsymbol{q}}\cdot{\boldsymbol{n}}}{N_{s}}}|{\boldsymbol{q}}a_{d}\rangle\,,\end{split}$
(29)

where $\int_{\boldsymbol{x}}=\int d^{2}{\boldsymbol{x}}$ and
$\int_{\boldsymbol{p}}=\int(2\pi)^{-2}d^{2}{\boldsymbol{p}}$ and we provide
the discretized version of the Fourier integrals. Using that

$\langle{\boldsymbol{x}}|{\boldsymbol{p}}\rangle=e^{-i{\boldsymbol{p}}\cdot{\boldsymbol{x}}}\,\to\,e^{-2\pi
i\frac{{\boldsymbol{n}}\cdot{\boldsymbol{q}}}{N_{s}}}\,,$ (30)

one can show that

$\langle{\boldsymbol{x}}|{\boldsymbol{y}}\rangle=\delta^{(2)}({\boldsymbol{x}}-{\boldsymbol{y}})=\frac{\delta_{{\boldsymbol{n}},{\boldsymbol{m}}}}{a_{s}^{2}}\,,$
(31)
$\langle{\boldsymbol{p}}|{\boldsymbol{k}}\rangle=(2\pi)^{2}\delta^{(2)}({\boldsymbol{k}}-{\boldsymbol{p}})=(2\pi)^{2}\frac{\delta_{{\boldsymbol{q}}_{\boldsymbol{k}},{\boldsymbol{q}}_{\boldsymbol{p}}}}{a_{d}^{2}}\,,$
(32)

where we used the closure identity

$\sum_{{\boldsymbol{n}}}e^{2\pi
i\frac{{\boldsymbol{n}}\cdot{\boldsymbol{q}}}{N_{s}}}=N_{s}^{2}\delta_{{\boldsymbol{q}},0}\,.$
(33)

We define the dimensionless basis states

$\begin{split}\ket{{\boldsymbol{n}}}=a_{s}\ket{{\boldsymbol{x}}}\,,\quad\ket{{\boldsymbol{q}}}=\frac{a_{d}}{2\pi}\ket{{\boldsymbol{p}}}\,,\end{split}$
(34)

which satisfy
$\langle{\boldsymbol{n}}|{\boldsymbol{m}}\rangle=\delta_{{\boldsymbol{n}},{\boldsymbol{m}}}$,
$\langle{\boldsymbol{q}}_{\boldsymbol{p}}|{\boldsymbol{q}}_{\boldsymbol{k}}\rangle=\delta_{{\boldsymbol{q}}_{\boldsymbol{p}},{\boldsymbol{q}}_{\boldsymbol{k}}}$
and $\langle{\boldsymbol{n}}|{\boldsymbol{q}}\rangle=N_{s}^{-1}\exp(-2\pi
iN_{s}^{-1}{\boldsymbol{n}}\cdot{\boldsymbol{q}})$. The Fourier transforms in
this normalization take the form

$\ket{{\boldsymbol{n}}}=\frac{1}{\sqrt{N_{s}^{2}}}\sum_{{\boldsymbol{q}}}e^{2\pi
i\frac{{\boldsymbol{q}}\cdot{\boldsymbol{n}}}{N_{s}}}\ket{{\boldsymbol{q}}}\,,$
(35)
$\ket{{\boldsymbol{q}}}=\frac{1}{\sqrt{N_{s}^{2}}}\sum_{{\boldsymbol{n}}}e^{-2\pi
i\frac{{\boldsymbol{q}}\cdot{\boldsymbol{n}}}{N_{s}}}\ket{{\boldsymbol{n}}}\,.$
(36)

It is also natural to introduce the operators
${\boldsymbol{P}}={\boldsymbol{p}}/a_{d}$ and
${\boldsymbol{X}}={\boldsymbol{x}}/a_{s}$, satisfying
$\hat{{\boldsymbol{X}}}\ket{{\boldsymbol{n}}}={\boldsymbol{n}}\ket{{\boldsymbol{n}}}$
and
$\hat{{\boldsymbol{P}}}\ket{{\boldsymbol{q}}}={\boldsymbol{q}}\ket{{\boldsymbol{q}}}$.
Inserting this operator definitions into Eq. 2, one can extract the
dimensionless Hamiltonian $H=a_{s}{\cal 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 ${\boldsymbol{n}}=(n_{1},n_{2})$ in
the binary basis, e.g.

$n_{1}=\sum_{i=0}^{2^{n_{Q}}-1}n_{1}^{(i)}2^{i}\,,$ (37)

where $n_{1}^{(i)}\in\\{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_{1}^{(i)}=0$ then we associate a qubit in the state
$\ket{\uparrow}=\ket{0}$ to it; conversely if $n_{1}^{(i)}=1$ we assign
$\ket{\downarrow}=\ket{1}$. Then, for example, the ket state
$\ket{{\boldsymbol{n}}}=\ket{3,3}$ with $n_{Q}=2$ is given by two registers
storing the overall state $\ket{1,1}\otimes\ket{1,1}$. Following Eqs. 35 and
36, the transformation between the momentum and position basis is achieved by
applying a standard quantum Fourier transform (qFT) NielsenChuang .

Finally, in this appendix and in the main text we have restricted ourselves to
considering lattices over positive integer values of ${\boldsymbol{n}}$ and
${\boldsymbol{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 Barata:2020jtq ) 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 operations999See for example Klco:2018zqz  for an example on how to
restrict the qFT to the first Brillouin zone., 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.

##  Appendix B Time evolution details

In this appendix we detail the key steps to implement the time evolution
operators $U_{K}$ and $U_{A}$. For convenience and clarity, and without loss
of generality, we will discuss both cases in one spatial dimension.

The strategy considered to implement $U_{K}$ was first discussed in Zalka_1998
. Starting from a state $\ket{{\boldsymbol{p}}}$ (with ${\boldsymbol{p}}$ now
having a single component) one wants to generate the state
$\exp(-is_{K}{\boldsymbol{p}}^{2})\ket{{\boldsymbol{p}}}$, with
$s_{K}=\varepsilon_{t}/(2E)$ a pure real number which can be easily computed
once all circuit parameters are fixed. This operation can be implemented by i)
adding an ancilla register with $l$ qubits all in state $\ket{0}$ ii) assuming
that a quantum black-box (quantum oracle) can be constructed that given
$\ket{{\boldsymbol{p}}}$ outputs
$\ket{F({\boldsymbol{p}})}=\ket{{\boldsymbol{p}}^{2}}$.

Regarding the first point, the value of $l$ solely depends on the numerical
accuracy one wants to represent ${\boldsymbol{p}}^{2}$ in a binary basis,
roughly $l\geq n_{Q}$. An efficient quantum oracle implementing the above
operation can always be constructed as long as a classical analog exists; this
is the case for the operation at hands.

Given both these conditions are satisfied, we then perform the following set
of operations

$\ket{{\boldsymbol{p}}}\otimes\ket{0}^{\otimes l}\stackrel{{\scriptstyle\rm
a_{1}}}{{\longrightarrow}}\ket{{\boldsymbol{p}}}\otimes\ket{F({\boldsymbol{p}})}\stackrel{{\scriptstyle\rm
a_{2}}}{{\longrightarrow}}\exp(-is_{K}F({\boldsymbol{p}}))\ket{{\boldsymbol{p}}}\otimes\ket{F({\boldsymbol{p}})}\stackrel{{\scriptstyle\rm
a_{3}}}{{\longrightarrow}}\exp(-is_{K}F({\boldsymbol{p}}))\ket{{\boldsymbol{p}}}\otimes\ket{0}^{\otimes
l}\,.$ (38)

Let us detail the above three steps. In a first step –$\rm a_{1}$– one applies
the quantum oracle, with input $\ket{{\boldsymbol{p}}}$ and stores the output
$F({\boldsymbol{p}})$ in the ancilla register. In step $\rm a_{2}$ one
performs a transformation of the form

$\ket{x}\to\exp(-is_{K}x)\ket{x}\,,$ (39)

with $s_{K}$ a real number and $\ket{x}$ denotes the binary decomposition,
with $l$ qubits, of an integer number. This exponentiation operation can
always be performed by applying $l$ single qubit gates
$R_{j}(\varepsilon)={\rm diag}(1,e^{-is_{K}2^{j}})$, taking into account that
$x$ can be decomposed as

$x=\sum_{j=0}^{l}x_{j}2^{j}\,,$ (40)

where $x_{j}\in\\{0,1\\}$. Acting on a single qubit the above operator has
non-zero matrix elements $\bra{0}R_{j}(s_{K})\ket{0}=1$ and
$\bra{1}R_{j}(s_{K})\ket{1}=\exp(-is_{K}2^{j})$; clearly stringing together
$l$ of such operators with increasing values of $j$

$R(s_{K})\equiv R_{0}(s_{K})\otimes R_{1}(s_{K})\otimes\cdots\otimes
R_{l}(s_{K})\,,$ (41)

results in a multi-qubit operator implementing the desired transformation,
i.e. $R(s_{K})\ket{x}=\exp(-is_{K}x)\ket{x}$.

The final step – $\rm a_{3}$ – consists in erasing the ancilla register back
to the state $\ket{0}^{\otimes l}$, which can be achieved by applying the
Hermitian conjugate circuit used in step $\rm 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 $\ket{{\boldsymbol{x}}}$ output
$\ket{A(t,{\boldsymbol{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}\times
N_{s}^{2}$ 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_{\alpha,\beta}\equiv{\rm diag}(\exp(i\alpha),\exp(i\beta))$.
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_{t}^{\rm th}$ time slice that
$\alpha_{k_{t}}=-g\epsilon_{t}A(k_{t}\cdot\varepsilon_{t},{\boldsymbol{0}})$
and $\beta_{k_{t}}=-g\varepsilon A(k_{t}\cdot\varepsilon_{t},\mathbf{1})$,
where the sub-index denotes the time slice and there are only two spatial
lattice points ($\ket{{\boldsymbol{0}}}$ and $\ket{\mathbf{1}}$). If we now
consider $n_{Q}=2$ but still a single spatial dimension, the respective time
evolution operator would be obtained by

$R_{\alpha,\beta}\otimes R_{\sigma,\gamma}=\begin{pmatrix}{\rm
e}^{i(\alpha+\sigma)}&0&0&0\\\ 0&{\rm e}^{i(\alpha+\gamma)}&0&0\\\ 0&0&{\rm
e}^{i(\beta+\sigma)}&0\\\ 0&0&0&{\rm e}^{i(\beta+\gamma)}\\\ \end{pmatrix}\,,$
(42)

for each time slice. By solving the associated system of linear equations, one
can map $\\{\alpha,\beta,\sigma,\gamma\\}$ to $\\{A(\mathbf{x})\\}$, which can
be done offline for any $t$ in a classical computer.

##  Appendix C Relation between $\ket{\psi_{L}}$ and the single particle
momentum distribution

In this appendix we relate $\ket{\psi_{L}}$ to the broadening distribution.

The single particle broadening probability for observing a quark with momentum
${\boldsymbol{k}}$ due to interactions with the medium for a time $L$ is given
by Blaizot:2012fh ; Blaizot:2015lma ; CasalderreySolana:2007zz

${\cal
P}(L,{\boldsymbol{k}})=\frac{1}{N_{c}}\int_{{\boldsymbol{x}},{\boldsymbol{y}}}{\rm
e}^{-i{\boldsymbol{k}}\cdot({\boldsymbol{x}}-{\boldsymbol{y}})}\text{Tr}\langle{\cal
W}(L,{\boldsymbol{x}}){\cal W}^{\dagger}(L,{\boldsymbol{y}})\rangle_{\rm
M}\,,$ (43)

where ${\cal W}(L,{\boldsymbol{x}})$ is a Wilson line operator along the
future pointing light-cone at a transverse position ${\boldsymbol{x}}$, which
can be written in the gauge choice employed in the main text as

${\cal W}(L,{\boldsymbol{x}})={\cal T}\exp\left(ig\int_{0}^{L}dt\,{\cal
A}^{-}(t,{\boldsymbol{x}})\cdot T\right)\,.$ (44)

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 GW ; HTL .

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 ${\cal 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
$\ket{\psi_{0}}=\ket{{\boldsymbol{p}}={\boldsymbol{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

$\ket{{\boldsymbol{0}}}\stackrel{{\scriptstyle\rm
qFT}}{{\longrightarrow}}\frac{1}{\sqrt{N_{s}^{2}}}\sum_{\boldsymbol{x}}\ket{{\boldsymbol{x}}}\stackrel{{\scriptstyle
U_{A}}}{{\longrightarrow}}\frac{1}{\sqrt{N_{s}^{2}}}\sum_{\boldsymbol{x}}U_{A}(L,{\boldsymbol{x}})\ket{{\boldsymbol{x}}}\stackrel{{\scriptstyle\rm
qFT^{\dagger}}}{{\longrightarrow}}\frac{1}{N_{s}^{2}}\sum_{\boldsymbol{q}}\left[\sum_{\boldsymbol{x}}U_{A}(L,{\boldsymbol{x}}){\rm
e}^{2\pi
i\frac{{\boldsymbol{x}}\cdot{\boldsymbol{q}}}{N_{s}}}\right]\ket{{\boldsymbol{q}}}\,.$
(45)

The probability of measuring the state $\ket{{\boldsymbol{k}}}$, ${\cal
P}_{\boldsymbol{k}}$, is simply given by

${\cal
P}_{\boldsymbol{k}}=\frac{1}{(N_{s}^{2})^{2}}\sum_{{\boldsymbol{x}},{\boldsymbol{y}}}{\rm
e}^{2\pi
i\frac{{\boldsymbol{k}}({\boldsymbol{x}}-{\boldsymbol{y}})}{N_{s}}}U_{A}^{\dagger}(L,{\boldsymbol{y}})U_{A}(L,{\boldsymbol{x}})\,.$
(46)

Averaging over all field configurations and noticing that ${\cal
W}({\boldsymbol{x}})=U_{A}^{\dagger}({\boldsymbol{x}})$ we obtain

${\cal
P}_{\boldsymbol{k}}=\frac{1}{(N_{s}^{2})^{2}}\sum_{{\boldsymbol{x}},{\boldsymbol{y}}}{\rm
e}^{2\pi
i\frac{{\boldsymbol{k}}({\boldsymbol{x}}-{\boldsymbol{y}})}{N_{s}}}\langle{\cal
W}(L,{\boldsymbol{y}}){\cal W}^{\dagger}(L,{\boldsymbol{x}})\rangle_{\rm
M}\,,$ (47)

which is just the discretized version of the single particle broadening
distribution ${\cal P}(L,{\boldsymbol{k}})$, as expected (ignoring the color
average, which can be performed as detailed in section 4). Also, since ${\cal
P}$ is a probability $\int_{\boldsymbol{k}}{\cal P}(L,{\boldsymbol{k}})=1$,
which is trivially true in the discrete version.

Notice that in our strategy, one does not need to explicitly provide a
prescription for the field correlators. Of course, these are embedded in the
generated field configurations and are taken into account (non-perturbatively)
by the time evolution operator.

##  Appendix D Measurement details

In this appendix we provide some details on the measurement protocol outlined
in the main text.

Taking the initial ancilla state to be $\ket{0}$, the measurement protocol
performs

$\begin{split}\ket{0}\ket{\psi_{L}}\stackrel{{\scriptstyle
V\,H}}{{\longrightarrow}}\frac{1}{\sqrt{2}}(\ket{0}\ket{\psi_{L}}+\ket{1}V\ket{\psi_{L}})\stackrel{{\scriptstyle
H}}{{\longrightarrow}}\frac{1}{2}\left[(1+V)\ket{0}\ket{\psi_{L}}+(1-V)\ket{1}\ket{\psi_{L}}\right]\,.\end{split}$
(48)

Then the expectation value for the random variable $\chi$ reads

$\begin{split}\langle\chi\rangle_{\rm
QM}=\frac{+1}{4}|\ket{\psi_{L}}+V\ket{\psi_{L}}|^{2}+\frac{(-1)}{4}|\ket{\psi_{L}}-V\ket{\psi_{L}}|^{2}=\frac{1}{2}\langle
V+V^{\dagger}\rangle_{\rm QM}\,,\end{split}$ (49)

which is equivalent to the expression in the main text.

The case where the initial ancilla state is $\ket{\gamma}\equiv
1/\sqrt{2}(\ket{0}+i\ket{1})$, which can be easily generated from the pure
state $\ket{0}$, reads

$\begin{split}\ket{\gamma}\ket{\psi_{L}}&\stackrel{{\scriptstyle
VH}}{{\longrightarrow}}\frac{1}{2}((1+i)\ket{0}\ket{\psi_{L}}+(1-i)\ket{1}V\ket{\psi_{L}})\\\
&\stackrel{{\scriptstyle
H}}{{\longrightarrow}}\frac{1}{\sqrt{8}}[((1+i)+(1-i)V)\ket{0}\ket{\psi_{L}}+((1+i)-(1-i)V)\ket{1}\ket{\psi_{L}}]\,.\end{split}$
(50)

Then the expectation value for $\chi$ reads

$\langle\chi\rangle_{\rm QM}=\frac{i}{2}\langle V^{\dagger}-V\rangle_{\rm
QM}\,,$ (51)

as indicated in the main text.