larger systems than classical computers can handle.

More specifically, here we focus on quantum simulation of chiral phase
transition, which is a phenomenon where the chirality of particles is
imbalanced due to temperature or density. The studies presented in this thesis
include quantum simulation of chiral phase transitions and studying chirality
imbalance with quantum algorithms for a low energy model of QCD, the Nambu-
Jona-Lasinio (NJL) model in 1+1 dimensions. These studies have potential
implications for future high-energy experiments and the development of new
theoretical frameworks in QCD.

####  2 Preliminaries of quantum computation

The basic unit of information in a quantum computer is the qubit, which can
exist in a superposition of two states, typically denoted as $|0\rangle$ and
$|1\rangle$ as depicted in fig. 1 shown on a bloch sphere. This means that a
qubit can exist in both states simultaneously, with a probability of being in
each state given by the amplitude of the corresponding state vector.
Furthermore, two or more qubits can be entangled, which means that the state
of one qubit is dependent on the state of the other, even if they are
physically separated. This property of entanglement allows quantum computers
to perform certain tasks exponentially faster than classical computers.

Figure 1: Quantum states $|0\rangle$ and $|1\rangle$ on the bloch sphere. Plot
generated from QISKIT [AAB19b] (IBM).

In the rest of this subsection, we introduce the four main postulates of
quantum mechanics related to Quantum Computing [NC00]. All postulates concern
closed quantum systems (i.e., systems isolated from environments) only.

#####  1 State space postulate.

The set of all quantum states of a quantum system forms a complex vector space
with inner product structure called the state space $\mathcal{H}$. If the
$\mathcal{H}$ is finite dimensional state space, it is isomorphic to some
$\mathbb{C}^{N}$, i.e. $\mathcal{H}\cong\mathbb{C}^{N}$. And one may naively
take $\mathcal{H}=\mathbb{C}^{N}$ without loss of generality. Always we assume
$N=2^{n}$ for some non-negative integer $n$ which is called the number of
quantum bits (qubits). In the Dirac notation, a quantum state
$|\psi\rangle\in\mathbb{C}^{N}$ and its Hermitian conjugate $\langle\psi|$ can
be expressed in terms of

$|\psi\rangle=\left(\begin{array}[]{c}\psi_{0}\\\ \psi_{1}\\\ \vdots\\\
\psi_{N-1}\end{array}\right)\,,\quad\langle\psi|=|\psi\rangle^{\dagger}=\left(\begin{array}[]{llll}{\psi}_{0}^{*}&{\psi}_{1}^{*}&\cdots&{\psi}_{N-1}^{*}\end{array}\right)\,.$
(1)

The inner product is thus

$\langle\psi|\varphi\rangle:=\langle\psi,\varphi\rangle=\sum_{i\in[N]}{\psi}_{i}^{*}\varphi_{i}\,,$
(2)

with $[N]=\\{0,\ldots,N-1\\}$. Let $\\{|i\rangle\\}$ be the standard basis of
$\mathbb{C}^{N}$. The $i$-th entry of $\psi$ can be written as an inner
product $\psi_{i}=\langle i|\psi\rangle$. Then $|\psi\rangle\langle\varphi|$
should be interpreted as an outer product, with $(i,j)$-th matrix element
given by

$\langle i|(|\psi\rangle\langle\varphi|)|j\rangle=\langle
i|\psi\rangle\langle\varphi|j\rangle=\psi_{i}\bar{\varphi}_{j}$ (3)

Two state vectors $|\psi\rangle$ and $c|\psi\rangle$ for some $0\neq
c\in\mathbb{C}$ always denote to the same physical state. In other words, the
complex scalar $c$ has no observable effects. Therefore, for simplicity, we
can always assume that $|\psi\rangle$ is normalized to be a unit vector,
meaning $\langle\psi|\psi\rangle=1$.

Sometimes, for convenience, we may prefer to work with unnormalized states,
denoted by $\psi$ without the ket notation $|\cdot\rangle$. However, when
considering normalized state vectors, the complex number
$c=e^{\mathrm{i}\theta}$, where $\theta\in[0,2\pi)$, is referred to as the
global phase factor.

As an example, one can define in a single qubit system corresponds to a state
space $\mathcal{H}\cong\mathbb{C}^{2}$ that

$|0\rangle=\left(\begin{array}[]{c}1\\\
0\end{array}\right),\quad|1\rangle=\left(\begin{array}[]{c}0\\\
1\end{array}\right)\,.$ (4)

Since the state space of the spin- $1/2$ system is also isomorphic to
$\mathbb{C}^{2}$, this is also called the single spin system, where
$|0\rangle$, $|1\rangle$ are referred to as the spin-up and spin-down state,
respectively. A general state vector in $\mathcal{H}$ takes the form

$|\psi\rangle=a|0\rangle+b|1\rangle=\left(\begin{array}[]{c}a\\\
b\end{array}\right)\,,\quad a,\ b\in\mathbb{C}\,,$ (5)

with the normalization condition $|a|^{2}+|b|^{2}=1$.

#####  2 Quantum operator postulate.

The evolution of a quantum state from
$|\psi\rangle\rightarrow\left|\psi^{\prime}\right\rangle\in$ $\mathbb{C}^{N}$
is always achieved via a unitary operator $U\in\mathbb{C}^{N\times N}$, i.e.,

$\left|\psi^{\prime}\right\rangle=U|\psi\rangle,\quad U^{\dagger}U=I_{N}$ (6)

Here $U^{\dagger}$ is the Hermitian conjugate of a matrix $U$, and $I_{N}$ is
the $N$-dimensional identity matrix. When the dimension is apparent, we may
also simply write $I\equiv I_{N}$. In quantum computation, a unitary matrix is
often referred to as a gate.

For example, for a single qubit, the Pauli matrices are

$\sigma_{x}=X=\left(\begin{array}[]{ll}0&1\\\
1&0\end{array}\right),\quad\sigma_{y}=Y=\left(\begin{array}[]{cc}0&-\mathrm{i}\\\
\mathrm{i}&0\end{array}\right),\quad\sigma_{z}=Z=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right).$ (7)

Together with the two-dimensional identity matrix, they form a basis of all
linear operators on $\mathbb{C}^{2}$ Some other commonly used single qubit
operators include, to name a few:

  * •

Hadamard gate

$H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\\
1&-1\end{array}\right)\,.$ (8)

  * •

Phase gate

$S=\left(\begin{array}[]{ll}1&0\\\ 0&\mathrm{i}\end{array}\right)\,.$ (9)

  * •

T gate:

$T=\left(\begin{array}[]{cc}1&0\\\ 0&e^{\mathrm{i}\pi/4}\end{array}\right)\,.$
(10)

To avoid notation conflicts, we will use the typewriter font such as
$\mathtt{H},\mathtt{X}$ for these single-qubit gates (one common scenario is
to distinguish the Hadamard gate $\mathtt{H}$ from a Hamiltonian $H$ ). An
operator acting on an $n$-qubit quantum state space is called an $n$-qubit
operator.

In the quantum circuit language, time flows from the left to right, i.e., the
input quantum state appears on the left, and the quantum operator appears on
the right, and each “wire” represents a qubit as shown in fig. 2. Note here
$\displaystyle|+\rangle=\frac{|0\rangle+|1\rangle}{\sqrt{2}}$ represents the
eigenstate of Pauli-$X$ matrix.

Figure 2: A few examples of a quantum circuit.

Starting from an initial quantum state $|\psi(0)\rangle$, the quantum state
can evolve in time, which gives a single parameter family of quantum states
denoted by $\\{|\psi(t)\rangle\\}$. These quantum states are related to each
other via a quantum evolution operator $U$ :

$\psi\left(t_{2}\right)=U\left(t_{2},t_{1}\right)\psi\left(t_{1}\right),$ (11)

where $U\left(t_{2},t_{1}\right)$ is unitary for any given $t_{1},t_{2}$. Here
$t_{2}>t_{1}$ refers to quantum evolution forward in ime, $t_{2}<t_{1}$ refers
to quantum evolution backward in time, and $U\left(t_{1},t_{1}\right)=I$ for
any $t_{1}$.

The quantum evolution under a time-independent Hamiltonian $H$ satisfies the
time-independent Schrödinger equation

$\mathrm{i}\partial_{t}|\psi(t)\rangle=H|\psi(t)\rangle$ (12)

Here $H=H^{\dagger}$ is a Hermitian matrix. The corresponding time evolution
operator is

$U\left(t_{2},t_{1}\right)=e^{-\mathrm{i}H\left(t_{2}-t_{1}\right)},\quad\forall
t_{1},t_{2}$ (13)

In particular, $U\left(t_{2},t_{1}\right)=U\left(t_{2}-t_{1},0\right)$. On the
other hand, for any unitary matrix $U$, we can always find a Hermitian matrix
$H$ such hat $U=e^{\mathrm{i}H}$.

#####  3 Quantum measurement postulate.

Next we shall focus our discussion solely on a specific type of quantum
measurements known as projective measurements, without losing generality. It
is worth noting that all quantum measurements, which can be described as
positive operator-valued measures (POVMs), can be equivalently represented
using projective measurements in an enlarged Hilbert space through the Naimark
dilation theorem.

In a finite-dimensional context, a quantum observable can always be associated
with a Hermitian matrix $M$, which possesses the spectral decomposition as
follows:

$M=\sum_{m}\lambda_{m}P_{m}$ (14)

Here $\lambda_{m}\in\mathbb{R}$ are the eigenvalues of $M$, and $P_{m}$ is the
projection operator onto the eigenspace associated with $\lambda_{m}$, i.e.,
$P_{m}^{2}=P_{m}$.

When a quantum state $|\psi\rangle$ is measured by a quantum observable $M$,
the outcome of the measurement is always an eigenvalue $\lambda_{m}$, with
probability

$p_{m}=\left\langle\psi\left|P_{m}\right|\psi\right\rangle$ (15)

After the measurement, the quantum state becomes

$|\psi\rangle\rightarrow\frac{P_{m}|\psi\rangle}{\sqrt{p_{m}}}$ (16)

Note that this is not a unitary process! In order to evaluate the expectation
value of a quantum observable $M$, we first use the resolution of identity:

$\sum_{m}P_{m}=I$ (17)

This implies the normalization condition,

$\sum_{m}p_{m}=\sum_{m}\left\langle\psi\left|P_{m}\right|\psi\right\rangle=\langle\psi|\psi\rangle=1.$
(18)

Together with $p_{m}\geq 0$, we find that $\left\\{p_{m}\right\\}$ is indeed a
probability distribution. The expectation value of the measurement outcome is

$\mathbb{E}_{\psi}(M)=\sum_{m}\lambda_{m}p_{m}=\sum_{m}\lambda_{m}\left\langle\psi\left|P_{m}\right|\psi\right\rangle=\langle\psi|\left(\sum_{m}\lambda_{m}P_{m}\right)|\psi\rangle=\langle\psi|M|\psi\rangle$
(19)

#####  4 Tensor product postulate.

For a quantum state consists of $m$ components with state spaces
$\left\\{\mathcal{H}_{i}\right\\}_{i=0}^{m-1}$, the state space is their
tensor products denoted by $\mathcal{H}=\otimes_{i=0}^{m-1}\mathcal{H}_{i}$.
Let $\left|\psi_{i}\right\rangle$ be a state vector in $\mathcal{H}_{i}$, then

$|\psi\rangle=\left|\psi_{0}\right\rangle\otimes\cdots\otimes\left|\psi_{m-1}\right\rangle$
(20)

in $\mathcal{H}$. However, not all quantum states in $\mathcal{H}$ can be
written in the tensor product form above. Let
$\left\\{\left|e_{j}^{(i)}\right\rangle\right\\}_{j\in\left[N_{i}\right]}$ be
the basis of $\mathcal{H}_{i}$, then a general state vector in $\mathcal{H}$
takes the form

$|\psi\rangle=\sum_{j_{0}\in\left[N_{0}\right],\ldots,j_{m-1}\in\left[N_{m-1}\right]}\psi_{j_{0}\cdots
j_{m-1}}\left|e_{j_{0}}^{(0)}\right\rangle\otimes\cdots\otimes\left|e_{j_{m-1}^{(m-1)}}^{(m)}\right\rangle.$
(21)

Here $\psi_{j_{0}\cdots j_{m-1}}\in\mathbb{C}$ is an entry of a $m$-way
tensor, and the dimension of $\mathcal{H}$ is therefore
$\prod_{i\in[m]}N_{i}$. The state space of $n$-qubits is
$\mathcal{H}=\left(\mathbb{C}^{2}\right)^{\otimes n}\cong\mathbb{C}^{2^{n}}$,
rather than $\mathbb{C}^{2n}$. We also use the notation

$|01\rangle\equiv|0,1\rangle\equiv|0\rangle|1\rangle\equiv|0\rangle\otimes|1\rangle,\quad\left|0^{\otimes
n}\right\rangle=|0\rangle^{\otimes n}$ (22)

Furthermore, $x\in\\{0,1\\}^{n}$ is called a classical bit-string, and
$\left\\{|x\rangle|x\in\\{0,1\\}^{n}\right\\}$ is called the computational
basis of $\mathbb{C}^{2^{n}}$.

As an example, in a two-qubit system corresponds to a state space
$\mathcal{H}=\left(\mathbb{C}^{2}\right)^{\times 2}\cong\mathbb{C}^{4}$, one
can define the standard basis

$|00\rangle=\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\
0\end{array}\right),\quad|01\rangle=\left(\begin{array}[]{c}0\\\ 1\\\ 0\\\
0\end{array}\right)\,,\quad|10\rangle=\left(\begin{array}[]{c}0\\\ 0\\\ 1\\\
0\end{array}\right)\,,\quad|11\rangle=\left(\begin{array}[]{c}0\\\ 0\\\ 0\\\
1\end{array}\right)\,.$ (23)

The Bell state (also called the EPR pair) is defined to be

$|\psi\rangle=\frac{1}{2}\left(|00\rangle+|11\rangle\right)=\frac{1}{2}\left(\begin{array}[]{c}1\\\
0\\\ 0\\\ 1\end{array}\right)\,,$ (24)

which cannot be written as any product state $|a\rangle\otimes|b\rangle$.
There are many important quantum operators on the two-qubit quantum system.
One of them is the CNOT gate, with matrix representation

$\mathrm{CNOT}=\left(\begin{array}[]{llll}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{array}\right).$ (25)

The quantum circuit for the CNOT gate is shown in fig. 3.

Figure 3: The quantum circuit for the CNOT gate

In other words, when acting on the standard basis, we have

$\operatorname{CNOT}\begin{cases}|00\rangle&=|00\rangle\\\
|01\rangle&=|01\rangle\\\ |10\rangle&=|11\rangle\\\
|11\rangle&=|10\rangle\end{cases}$ (26)

This can be compactly written as

$\mathrm{CNOT}|a\rangle|b\rangle=|a\rangle|a\oplus b\rangle\,.$ (27)

Here $a\oplus b=(a+b)\bmod 2$ is the “exclusive or” (XOR) operation.

####  3 Applications in quantum field theory

Noisy Intermediate-Scale Quantum (NISQ) technology has demonstrated its
ability to solve complex problems such as real-time dynamics [CTG19, SKP19,
ZPH17, al13, FFK20, Fey82, Llo96, DMM21b, JLM21], relativistic behaviors,
many-body systems [WSB04, Maj07, JLP12c, ZCR12, ZCR13, BBD13, BDM12, Wie13,
Wie14, JLP14a, GCM15, MWR14, BMT15, ZCR15, MRS15, DM16, ZFR17, MMS16, BAM17,
GCS17, KSP18, MSA18, ZHJ18, ZUZ18, KDM18b, KS19b, GMU19, NAB19, MDF20, JLP19,
LKL19, KS20b, LL18, KSS20b, ABH19, MTV20c, LLY20b, CHI20, BTR18, ZTL20,
ZTL21], ground state estimations [Aru20, MGG20, KTC19, OM16, KMT17, Per14,
Col18], and finite-temperature properties of various systems [BBM20, TD00b,
PW09, RGE12, TOV11a, YA12, LGL22, ZXY21]. While digital quantum simulations of
thermal physical systems have been studied previously, finite-temperature
physics on quantum computers remains a challenging area with limited
understanding [SMT21]. Recently, various algorithms for imaginary time
evolution on quantum computers have been introduced, including the Quantum
Imaginary Time Evolution (QITE) algorithm [MST20, MJE19, YEZ19, BMG19, NKM21,
GZB20, YSP21, YPS20], which uses a unitary operation to simulate imaginary
time evolution. QITE has been applied to calculate finite-temperature
observables such as energy [SMT21, Vil21], magnetization in the Transverse
Field Ising Model (TFIM) [Vil21], and more. These advancements open up new
opportunities for quantum simulation in various fields, including materials
science, condensed matter physics, and quantum chemistry.

Quantum field theory is the branch of theoretical physics that deals with the
study of quantum mechanics applied to fields, such as the electromagnetic
field and the Higgs field. It is the framework used to describe the behavior
of subatomic particles, including quarks and leptons. Quantum field theory
plays a crucial role in the study of high-energy physics and hadron physics.

Quantum computing has the potential to provide significant improvements in the
field of quantum field theory, especially in the study of high-energy physics
and hadron physics. One area where quantum computing can be useful is in
simulating the behavior of subatomic particles, which is a computationally
intensive task that is currently beyond the reach of classical computers. By
leveraging the power of quantum computing, researchers can simulate larger and
more complex systems, enabling more accurate predictions of the behavior of
particles.

Another area where quantum computing can be useful is in solving certain
optimization problems that arise in quantum field theory. For example, the
problem of finding the ground state of a system of interacting particles is an
important problem in quantum field theory. This problem can be reduced to an
optimization problem, which can be solved using quantum annealing, a specific
form of quantum computing that is well-suited for solving optimization
problems.

In theoretical physics, the complex phase structure of strongly interacting
matter, as described by Quantum Chromodynamics, holds great significance
[Fuk08, CdR07, FZL08, XCW15, JLZ11, ZCJ14, SWJ14]. In the infinite quark
limit, the deconfinement phase transition is expected as a result of
spontaneous $\mathbb{Z}_{N}$ symmetry breaking, with the transition occurring
as a function of temperature [RB21, HW00]. Alternatively, in the massless
quark limit, the chiral phase transition has been extensively studied and has
an order parameter [Raj95, BEF10, Baz19, BFG20, Din19, DLM21]. QCD is
dominated by non-perturbative effects at low energies, which makes it
challenging to study the chiral phase transition of strongly interacting
matter. Furthermore, other phases of strongly interacting matter are thought
to exist, such as the Color-Flavor Locked (CFL) phase, which occupy the high
density, low temperature region of the phase diagram.

In conclusion, quantum computing has the potential to revolutionize many
fields of science, including physics, chemistry, and cryptography. In the
field of quantum field theory, quantum computing can provide significant
improvements in the study of high-energy physics and hadron physics by
enabling more accurate predictions of the behavior of subatomic particles and
solving optimization problems. As quantum computing technology continues to
evolve, it is likely that we will see more applications in quantum field
theory and other areas of physics.

#####  1 Discretization of the Continuous Theory

In order to address the non-perturbative nature of lattice Quantum
Chromodynamics, a transformation is employed to promote all SU(3)-valued gauge
fields $A_{\mu}(x)$ to elements in the SU(3) group, denoted as $U_{\mu}(x)$.
This is accomplished through the expression:

$U_{\mu}(x)=\exp\left[\mathrm{iag}_{s}A_{\mu}(x)\right]\,,$ (28)

where $a$ and $g_{s}$ correspond to the lattice spacing and the strong
coupling constant, respectively. Throughout the remainder of this document,
the focus will primarily be on the treatment of $U_{\mu}(x)$.

To facilitate further analysis, a Wick rotation is performed, resulting in the
substitution:

$x_{0}\rightarrow-ix_{4}\,.$ (29)

As a consequence, the metric takes on a Euclidean form:

$x^{2}=x_{\mu}x_{\mu}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}\,,$ (30)

and the weight factor in the path integral becomes real as well,

$\left\langle
0\left|T\left\\{\psi\left(x_{1}\right)\bar{\psi}\left(x_{2}\right)\right\\}\right|0\right\rangle=\frac{1}{Z}\int[d\bar{\psi}][d\psi]\left[dU_{\mu}\right]\psi\left(x_{1}\right)\bar{\psi}\left(x_{2}\right)\exp\left[-S\left[U_{\mu},\bar{\psi},\psi\right]\right]\,.$
(31)

Here the symbol $S$ is still used to denote the Euclidian action.

The fermionic component of Lagrangian Equation eq. 1 can also be expressed in
a straightforward manner as

$S_{F}=a^{3}a_{0}\sum_{n}\bar{\psi}(n)\left(\sum_{\mu=1}^{4}i\gamma_{\mu}\frac{U_{\mu}(n)\psi(n+\hat{\mu})-U_{-\mu}(n)\psi(n-\hat{\mu})}{2a}+m\psi(n)\right)$
(32)

Consequently, the lattice Dirac matrix $D(n\mid m)$ can be written as

$D(n\mid m)_{\alpha i,\beta
j}=\sum_{\mu=1}^{4}i\left(\gamma_{\mu}\right)_{\alpha\beta}\frac{U_{\mu}(n)_{ij}\delta_{n+\hat{\mu},m}-U_{-\mu}(n)_{ij}\delta_{n-\hat{\mu},m}}{2a}+m\delta_{\alpha\beta}\delta_{ij}\delta_{nm}\,,$
(33)

where we explicitly indicate the lattice sites $m$ and $n$, Dirac indices
$\alpha$ and $\beta$, as well as color indices $a$ and $b$. However, this
naive form of the Dirac matrix exhibits the fermion doubling problem [NN81],
characterized by the presence of unphysical poles. Various well-established
methods have been developed to resolve this issue, including Wilson fermions
[Wil74], Kogut-Susskind staggered fermions [KS75], and domain wall fermions
[Kap92].

After addressing the fermion doubling problem, it becomes advantageous to
integrate out the fermionic degrees of freedom in the path integral by
employing Gaussian integration over Grassmann numbers. This allows us to
obtain a path integral solely in terms of gauge fields. For instance, the
partition function can be written as:

$Z=\int\mathcal{D}[U]\operatorname{det}[D]e^{-S_{W}[U]}$ (34)

Here, the matrix $D$ represents the Dirac matrix with a chosen method for
handling the fermion doubling problem. The size of the Dirac matrix scales
linearly with the lattice size, resulting in a computational cost
significantly higher than that of pure Yang-Mills theory.

With the path integral expressed solely in terms of gauge fields, we can now
proceed to numerically evaluate it. Given the high dimensionality of the
integral, it is practical to employ Monte Carlo methods with importance
sampling. In the case of thermal calculations without real-time evolution, the
weight factor $e^{-S}$ is always real and positive, enabling it to be
interpreted as the probability distribution function for the configuration of
gauge fields. Through Markov chain Monte Carlo (MCMC) methods, we sample
$N_{s}$ configurations of gauge links across the lattice according to the
distribution function $e^{-S}$.

Given that $S\left[U_{\mu},\bar{\psi},\psi\right]$ is a real valued number,
the path integral formulation defines a probability that can be used in a
Monte Carlo algorithm. Samples of $U_{\mu},\bar{\psi}$ and $\psi$ are drawn
according to the probability exp
$\left[-S\left[U_{\mu},\bar{\psi},\psi\right]\right]$ and the desired quantity
$\psi\left(x_{1}\right)\bar{\psi}\left(x_{2}\right)$ is measured on these
samples. By taking the average of the values we acquire an approximation of
the true correlation function

$\langle\mathcal{O}\rangle\simeq\frac{1}{n}\sum_{j}[\mathcal{O}]_{j}\,.$ (35)

with the notation $[\mathcal{O}]_{j}$ measured on the $j$-th sample. According
to the strong law of large numbers, the right-hand side of the expression
serves as an unbiased estimator for the true expectation, and its variance
decreases as $\mathcal{O}(1/n)$ as the sample size $n$ grows larger.

#####  2 Staggered fermions

Staggered fermions are a specific type of fermion discretization used in
lattice QCD. They involve a clever arrangement of fermion fields on the
lattice that exploits the taste symmetry to represent multiple fermion
flavors. In staggered fermions, each lattice site is associated with a single
fermion flavor, reducing the computational cost compared to other
discretization schemes. This approach retains chiral symmetry at the level of
the lattice action, making staggered fermions particularly advantageous for
large-scale QCD simulations.

The transformation from a Dirac fermion field to a staggered fermion field is
achieved through the process known as rooting or taking the nth root of the
fermion determinant. By taking the nth root, where n is the number of fermion
flavors, each flavor of the original Dirac fermion field is associated with a
distinct staggered fermion field. This procedure effectively eliminates
unwanted degrees of freedom known as doublers that arise due to the lattice
discretization, resulting in a reduced number of fermion degrees of freedom.
The resulting staggered fermion formulation retains essential properties of
QCD and enables efficient numerical simulations on the lattice grid.

More specifically, the method of staggered fermion fields can be described
using the following steps:

  1. 1.

The transformation from a Dirac fermion field to a staggered fermion field
involves rooting:

$\Psi(x)\rightarrow\sqrt[n]{\det(D)}\Psi_{s}(x)\,.$ (36)

Here, $\Psi(x)$ represents the original Dirac fermion field, $\Psi_{s}(x)$
denotes the staggered fermion field, and $D$ represents the fermion operator
on the lattice.

  2. 2.

The staggered fermion fields can be expressed in terms of staggered spinors:

$\Psi_{s}(x)=\sum_{\alpha}\chi_{\alpha}(x)\eta_{\alpha}\,.$ (37)

In this equation, $\chi_{\alpha}(x)$ represents the staggered spinor
associated with the lattice site $x$ and flavor $\alpha$, while
$\eta_{\alpha}$ denotes the staggered spinor phase factor.

  3. 3.

The staggered fermion action, incorporating the hopping and mass terms, can be
written as:

$S=\sum_{x,y,\alpha}\bar{\chi}_{\alpha}(x)D_{xy}\chi_{\alpha}(y)\,.$ (38)

Here, $\bar{\chi}_{\alpha}(x)$ represents the staggered spinor conjugate, and
$D_{xy}$ represents the staggered fermion hopping matrix on the lattice.

The rooting transformation reduces the number of fermion degrees of freedom,
and the staggered fermion fields are represented using staggered spinors,
enabling the formulation of the staggered fermion action for lattice QCD
simulations.

In summary, staggered fermion fields are primarily employed in lattice QCD
simulations as a discretization scheme for representing fermions on a lattice.
They possess taste symmetry and are computationally efficient. Staggered
fermions provide a way to describe multiple fermion flavors on the lattice and
retain certain symmetries of the continuum QCD theory. While staggered
fermions are widely used in lattice QCD, they are not directly related to the
Jordan-Wigner transformation.

#####  3 Jordan-Wigner transformation

The Jordan-Wigner transformation [JW28] is a powerful mathematical technique
used to establish a connection between lattice fermions and quantum spin
systems. By employing this transformation, fermion creation and annihilation
operators are mapped onto spin operators, enabling the study of lattice QCD
and other quantum field theories through the analogy of condensed matter
systems. This transformation plays a crucial role in bridging different
physical systems, allowing researchers to explore QCD phenomena by leveraging
techniques and insights from condensed matter physics. The Jordan-Wigner
transformation provides a valuable tool for investigating the behavior of
lattice fermions and advancing our understanding of fundamental interactions
in high-energy physics.

In practical terms, the Jordan-Wigner transformation involves expressing
fermion operators in terms of spin operators by introducing string operators
that keep track of the fermion parity. This transformation enables the
formulation of lattice fermion systems as spin models, which can be analyzed
using powerful techniques from condensed matter physics. Through the
application of the Jordan-Wigner transformation, researchers gain a deeper
understanding of lattice QCD dynamics and can explore various phenomena by
drawing analogies to condensed matter systems. This transformative technique
broadens the scope of investigation and facilitates the study of complex
quantum field theories, ultimately contributing to advancements in our
understanding of fundamental particles and their interactions.

In the Jordan-Wigner transformation, one first maps fermion creation and
annihilation operators onto spin operators:

$c_{j}^{\dagger}=\prod_{k<j}\sigma_{k}^{z}\sigma_{j}^{+}\,,\quad
c_{j}=\prod_{k<j}\sigma_{k}^{z}\sigma_{j}^{-}\,,$ (39)

where $c_{j}^{\dagger}$ and $c_{j}$ represent the creation and annihilation
operators for fermions at lattice site $j$, respectively, and $\sigma_{j}^{+}$
and $\sigma_{j}^{-}$ are the spin raising and lowering operators at the
corresponding site.

Next, the fermionic number operator can be expressed in terms of spin
operators:

$n_{j}=\frac{1}{2}\left(1-\sigma_{j}^{z}\right)\,,$ (40)

with $n_{j}$ representing the fermionic number operator at site $j$, and
$\sigma_{j}^{z}$ is the spin operator corresponding to the $z$-component of
spin at that site.

Then the fermionic hopping term can be written in terms of spin operators:

$c_{i}^{\dagger}c_{j}=\frac{1}{4}\left(\sigma_{i}^{x}-i\sigma_{i}^{y}\right)\left(\sigma_{j}^{x}+i\sigma_{j}^{y}\right)\,.$
(41)

Here, $c_{i}^{\dagger}$ and $c_{j}$ are the creation and annihilation
operators for fermions at lattice sites $i$ and $j$, respectively, while
$\sigma_{i}^{x}$, $\sigma_{i}^{y}$, $\sigma_{j}^{x}$, and $\sigma_{j}^{y}$ are
spin operators corresponding to the $x$ and $y$ components of spin at the
respective sites.

These equations illustrate the mapping of fermion operators to spin operators
through the Jordan-Wigner transformation. By expressing fermionic creation,
annihilation, and number operators in terms of spin operators, researchers can
analyze lattice fermion systems using the formalism of spin models,
facilitating the exploration of quantum field theories and condensed matter
physics phenomena.

The Jordan-Wigner transformation establishes a connection between lattice
fermions and quantum spin systems. It is a mathematical mapping that enables
the representation of fermion operators in terms of spin operators. The
Jordan-Wigner transformation is used in various areas of theoretical physics,
particularly in the study of quantum field theories and condensed matter
physics. In the context of lattice QCD, the Jordan-Wigner transformation is
not typically employed directly, as other discretization schemes, such as
staggered fermions, are more commonly utilized.

To discretize a fermion field, a two-step process involving the transformation
into staggered fermion fields and subsequent application of the Jordan-Wigner
transformation can be employed. Initially, the fermion field is transformed
into staggered fermion fields using a rooting procedure. This involves taking
the $n$th root of the fermion determinant, where $n$ corresponds to the number
of fermion flavors. By doing so, each flavor of the original fermion field is
associated with a distinct staggered fermion field, effectively reducing the
number of fermion degrees of freedom.

Once the staggered fermion fields are obtained, the Jordan-Wigner
transformation is applied to each of these fields individually. This
transformation maps the staggered fermion fields to quantum spin systems,
allowing for the representation of the staggered fermion operators in terms of
spin operators. This enables the study of the fermionic lattice system through
the analogy of condensed matter spin models. By applying the Jordan-Wigner
transformation to each staggered fermion field, researchers can explore the
behavior of lattice fermions using the language and techniques of condensed
matter physics, gaining insights into the underlying physics of the system.

####  4 Study QCD chiral phase diagram with Quantum Simulations

#####  1 QCD phase diagram

The study of the QCD phase diagram has experienced tremendous progress with
lattice simulations in the last few years. The phase transition line is
typically determined by extrapolating chiral observables to finite chemical
potential $\mu_{B}$, and finding the temperature at which the chiral
condensate has an inflection point, or the chiral susceptibility has a peak.
The transition temperature as a function of $\mu_{B}$ can be written as
${T_{c}(\mu_{B})}/{T_{c}(\mu_{B}=0)}=1-\kappa_{2}\left({\mu_{B}}/{T_{c}(\mu_{B})}\right)^{2}-\kappa_{4}\left({\mu_{B}}/{T_{c}(\mu_{B})}\right)^{4}+\cdots$.
A high-precision result for the crossover temperature $T_{c}(\mu_{B}=0)$ has
become available [BFG20, Baz19, BFG20]. Similar coefficients for the
extrapolation of the transition temperature to finite strangeness, electric
charge and isospin chemical potentials were obtained in Ref. [Baz19]. No sign
of criticality is observed from lattice QCD simulations up to $\mu_{B}\simeq
300$ MeV [BFG20, Din19]. Future challenges include the extrapolation of the
phase transition line to larger values of chemical potential and more
stringent constraints on the location of the critical point.

Based on lattice QCD calculations, the QGP-hadron gas transition at vanishing
net-baryon density is understood to be a smooth crossover with the transition
temperature $T_{c}=156\pm 1.5$ MeV [AEF06]. Model studies indicate a first-
order phase boundary at large net-baryon density (baryon chemical potential
$\mu_{B}$) [FS13]. If there is a crossover and a first order transition line,
they will be joined at the QCD critical point [SRS98, SRS99, BEK20]. State-of-
the-art lattice calculations further predicted that the chiral crossover
region extends into the finite chemical potential region $\mu_{B}/T\leq 2$
[Baz20], see fig. 4. Precise calculations in the higher $\mu_{B}$ region
become more difficult and experimental measurements are essential to determine
if a QCD critical point exists.

Figure 4: Sketch of the QCD phase diagram, incorporating a conjectured
critical end point and first order transition. The yellow line indicates the
region of the phase diagram where lattice QCD can reliably predict the smooth
crossover region of the hadron-QGP transition, up to $\mu_{B}/T\leq 2$. Figure
adapted from [Apr15a].

#####  2 Quantum Simulations

In order to study the structure of the QCD chiral phase diagram, which divides
the temperature-chemical potential plane into regions based on the physical
properties of strongly interacting matter or “phases”, researchers typically
employ lattice QCD calculations. However, this method has limitations due to
the fermion sign problem. At nonzero baryochemical potential, the QCD action
is no longer necessarily real, making the use of a Monte Carlo evaluation of
thermal expectation values based on $e^{-S}$ impossible [RB21, Phi07]. To
address this issue, simulations on imaginary chemical potentials are conducted
[FP02, DL03, WLC07, DDL07, CD07, FP08, DS09, MWL10], and the expectation value
of a given observable is expanded in powers of $\mu/T$ [AEH02, ADE05, GG08,
Bas08, KKL11]:
$\langle\mathcal{O}\rangle=\sum_{k}\mathcal{O}_{k}(T)(\mu/T)^{k}$, then the
coefficients $\mathcal{O}_{k}(T)$ are evaluated on the lattice [RB21, ADE05,
Phi07]. This approach is limited to a potentially small range of baryochemical
potentials and assumes that the expectation values are analytic in $\mu$
[RB21]. This assumption is unlikely to hold for some observable in the
vicinity of a finite-order phase transition. Despite this obstacle,
researchers can use quantum computers to model the lattice system, eliminating
the need for a Monte Carlo analysis by taking advantage of the statistical
properties of the quantum computer. This approach has been suggested by
Feynman [Fey82] and has been shown in [KK20b] to avoid the sign problem.

In the so-called near-term noisy intermediate quantum (NISQ) era [Pre18],
directly simulating QCD is not possible. Instead, we will study chiral phase
transition via a low energy model of QCD, the Nambu-Jona-Lasinio (NJL) model
in 1+1 dimension. The NJL model has proven to be a practical and convenient
tool to investigate the QCD chiral phase transition, owing to its ability to
provide insight into the mechanism of chiral symmetry breaking and access to
the dynamical mass [Fuk08, CdR07, LDC15, DLX15, CHS15, DCX13, SYX15, MPA09b,
MPA09a, GMC07, BBH12]. This effective model, tailored for low-energy two-
flavor QCD, can also be utilized for analytical calculations at finite
temperature and chemical potential. An even simpler version of the NJL model
is the Gross-Neveu (GN) model [GN74], a renormalizable and asymptotically free
$(1+1)$-dimensional theory that is composed of $N$ fermion species interacting
via four-fermion contact interaction, and it exhibits chiral symmetry breaking
at low energy by means of the generation of a nonvanishing chiral condensate
$\langle\bar{\psi}\psi\rangle$.

In our work [CKM22, CKT22], we aim to examine the chiral condensate of the GN
model at finite temperature and chemical potential [Thi20] and compare our
findings with quantum simulations. In particular, we illustrate how the
quantum algorithm can be used to study the chiral phase transition of the GN
model in 1+1 dimension. In this section, we provide our study that leverages a
4-qubit system with each staggered fermion field located in one qubit to
simulate the $(1+1)$ dimensional NJL Hamiltonian by implementing single-qubit
gates and the CNOT gate via the Jordan-Wigner transformation [JW28]. This
Hamiltonian, when approached through its continuum limit, corresponds to the
NJL model. We will provide detailed quantum simulations in the next
subsections, here let us emphasize the following important points:

  * •

Ultimately, to make contact with the continuum field theory, any such
simulation will have to be performed on a series of increasing lattices, and
the result extrapolated to that where the lattice sites $N\to\infty$ and the
lattice spacing $a\to 0$. Any parameters of the theory present in the
continuum must be suitably matched for this procedure to yield meaningful
results. Standard procedures for how this should be done derived from the
lattice QCD community can be accessed in reputable sources [BBD98, FRG10,
Phi21], and for the NJL model in 3+1 dimensions in [Wal03]. In this chapter,
we work at fixed lattice size $a$ with a small number of lattice sites $N$ and
we leave the detailed investigation of these issues to future work. In this
regard, our work [CKM22, CKT22] is primarily centered around the formulation
and development of the quantum algorithm to illustrate how this can be done
for a system with a small number of qubits.

  * •

In the community, there has been a lot of important studies that use
relatively small number of qubits, e.g. [Cle20, XGX22, BJN21], which allow
further development on the actual quantum algorithms and directly use it in
the current simulations. While its contributions have yet to yield novel
insights into new realms of physics, the concerted efforts in this direction
have nurtured a fertile ground for the exploration of new research
possibilities.

The rest of this subsection is organized as follows: In section 3, we furnish
analytical calculations to explore the finite temperature behavior of the NJL
model with consideration of chemical potentials. Sec. 4 provides an overview
of the discretization of the field theory on a lattice and section 5 shows the
quantum algorithm used, namely the QITE algorithm. Subsequently, we present a
comparison between the results obtained by analytical calculations, exact
diagonalization, and quantum simulation at various chemical potentials in
section 6, where we find a remarkable consistency between them.

#####  3 Investigated model: Nambu-Jona-Lasinio model

This section provides a comprehensive exposition on the Nambu-Jona-Lasinio
(NJL) model, focusing on its phase transition in the context of finite
temperature and chemical potential.

The NJL model is governed by a Lagrangian density in $(1+1)$-dimensional
Minkowski space, which takes the form of the following expression: [NJ61b,
NJ61c]

$\displaystyle\mathcal{L}_{\textrm{N}JL}$
$\displaystyle=\bar{\psi}(i\not{\partial}-m)\psi+g\left[(\bar{\psi}\psi)^{2}+\bar{\psi}i\gamma_{5}\tau_{a}\psi)^{2}\right]\,,$
(42)

where the operator $\not{\partial}=\gamma_{\mu}\partial^{\mu}$, the $\tau_{a}$
are the Pauli matrices in isospin space, $m$ denotes the bare quark mass, and
$g$ is the dimensionless coupling constant. Notably, the gamma matrices
$\gamma_{0}$, $\gamma_{1}$, and $\gamma_{5}$ in $(1+1)$ dimensions are given
by the following expressions:

$\displaystyle\gamma_{0}=Z,\quad\gamma_{1}=-iY,\quad\gamma_{5}=\gamma_{0}\gamma_{1}=-X\,,$
(43)

where, $X$, $Y$, and $Z$ represent the Pauli matrices:

$\displaystyle X=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix},\
Y=\begin{pmatrix}0&-i\\\ i&0\end{pmatrix},\ Z=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}$ (44)

A simplified version of the NJL model is the Gross-Neveu (GN) model [GN74],
which is governed by the following Lagrangian expression:

$\displaystyle\mathcal{L_{\textrm{G}N}}$
$\displaystyle=\bar{\psi}(i\not{\partial}-m)\psi+g(\bar{\psi}\psi)^{2}\,.$
(45)

One can also include the chiral chemical potential $\mu_{5}$ that simulates
the chiral asymmetry between right- and left-handed quarks coupled with the
chirality charge density operator $n_{5}=\bar{\psi}\gamma_{0}\gamma_{5}\psi$
and obtain:

$\displaystyle\mathcal{L}=$
$\displaystyle\bar{\psi}(i\not{\partial}-m)\psi+g(\bar{\psi}\psi)^{2}+\mu\bar{\psi}\gamma_{0}\psi+\mu_{5}\bar{\psi}\gamma_{0}\gamma_{5}\psi\,.$
(46)

We will study the QCD phase diagram at both $\mu_{5}=0$ and $\mu_{5}\neq 0$ in
this section.

###### QCD phase diagram at finite chemical potential

To begin with, we study the dependence on chemical potential and provide an
analytical calculation with zero chiral potential. In the limit where the bare
quark mass is zero, the Lagrangian satisfies a discrete symmetry group known
as $\mathbb{Z}_{2,L}\times\mathbb{Z}_{2,R}\equiv{\pm 1,\pm\gamma_{5}}$, under
which the fields transform by a left-action, i.e. $\psi\mapsto G\cdot\psi$ for
any group element $G\in\mathbb{Z}_{2,L}\times\mathbb{Z}_{2,R}$. In this work,
we investigate the chiral phase transition of the Gross-Neveu (GN) model with
a non-zero chemical potential $\mu$ [Thi20]:

$\displaystyle\mathcal{L}$
$\displaystyle=\bar{\psi}(i\not{\partial}-m)\psi+g(\bar{\psi}\psi)^{2}+\mu\bar{\psi}\gamma_{0}\psi\,,$
(47)

This model, which is a simplified version of the Nambu-Jona-Lasinio (NJL)
model, has been previously utilized for studying the chiral phase transition
at finite values of both temperature $T$ and chemical potential $\mu$. For the
sake of consistency throughout this section, we will refer to the Lagrangian
in eq. 47 as the NJL model.

In order to discern between different phases, we investigate the chiral
condensate $\langle\bar{\psi}\psi\rangle$ in the vacuum, which is widely
recognized as an order parameter [GPY81, NJ61b, MS81, Pol78, FWZ18] in the
chiral limit with zero quark mass. As $\bar{\psi}\psi$ transforms nontrivially
under $\mathbb{Z}_{2,L}\times\mathbb{Z}_{2,R}$, a nonzero value of
$\langle\bar{\psi}\psi\rangle$ requires spontaneous breaking of this symmetry,
thus a transition from non-vanishing to vanishing values of the condensate
indicates a chiral phase transition. Hence, with non-zero quark mass, the
chiral condensate $\langle\bar{\psi}\psi\rangle$ can be considered as a quasi-
order parameter [BIL18].

The transition between the chirally symmetric and broken phases can be
examined in the mean field approximation. This part is widely studied and its
results are also well-known in the literature. Below, we review the important
results in details in order to validate our quantum simulation results shown
in the next section. It is important to emphasize again that the purpose is to
illustrate how specific quantum algorithms can be achieved to simulate the NJL
model in 1+1 dimensions in quantum circuits, and that it is not meant to
provide new physics insights. In the mean field approximation, the quantity
$\bar{\psi}\psi$ can be expressed as
$\bar{\psi}\psi=\langle\bar{\psi}\psi\rangle+\sigma$ [GN74, Wal74], where
$\sigma$ is a small, real scalar field (“fluctuations”), i.e.,
$|\sigma/\langle\bar{\psi}\psi\rangle|\ll 1$, and
$\langle\bar{\psi}\psi\rangle$ is a constant that is independent of
coordinates. $\langle\bar{\psi}\psi\rangle$ depends on both the chemical
potential $\mu$ and the temperature $T$. Thus

$\displaystyle\left(\bar{\psi}\psi\right)^{2}$
$\displaystyle=\left(\langle\bar{\psi}\psi\rangle+\sigma\right)^{2}=\langle\bar{\psi}\psi\rangle^{2}+2\sigma\langle\bar{\psi}\psi\rangle+\mathcal{O}(\sigma^{2})$
$\displaystyle=\langle\bar{\psi}\psi\rangle^{2}+2\left(\bar{\psi}\psi-\langle\bar{\psi}\psi\rangle\right)\langle\bar{\psi}\psi\rangle+\mathcal{O}(\sigma^{2})$
$\displaystyle=-\langle\bar{\psi}\psi\rangle^{2}+2\bar{\psi}\psi\langle\bar{\psi}\psi\rangle+\mathcal{O}(\sigma^{2})\,,$
(48)

and the Lagrangian is expressed as follows:

$\displaystyle\mathcal{L}=\bar{\psi}(i\not{\partial}-m+2g\langle\bar{\psi}\psi\rangle+\mu\gamma_{0})\psi-g\langle\bar{\psi}\psi\rangle^{2}+\mathcal{O}(\sigma^{2})\,.$
(49)

By disregarding the terms of order $\mathcal{O}(\sigma^{2})$, we obtain the
Lagrangian $\mathcal{L}_{\mathrm{Dirac}}$ that describes a free Dirac fermion
(subject to finite chemical potential) possessing a mass
$M=m-2g\langle\bar{\psi}\psi\rangle$, along with a constant potential
$\mathcal{V}=g\langle\bar{\psi}\psi\rangle^{2}=(M-m)^{2}/4g$. This yields an
efficient, linearized Lagrangian that reads as follows:

$\displaystyle\mathcal{L}_{\mathrm{eff}}=\bar{\psi}(i\not{\partial}-M+\mu\gamma_{0})\psi-\frac{(M-m)^{2}}{4g}=\mathcal{L}_{\mathrm{Dirac}}-\mathcal{V}\,,$
(50)

In consequence, the Lagrangian is now an effective one, denoted as
$\mathcal{L}_{\mathrm{eff}}$, with the small fluctuation terms represented by
$\mathcal{O}(\sigma^{2})$, and the mass of the fermion field $M$ is now a
function of the chemical potential $\mu$ and temperature $T$, namely
$M=M(\mu,T)$.

The determination of the chiral condensate $\langle\bar{\psi}\psi\rangle$, or
equivalently the effective mass $M$, can be achieved by minimizing the Grand
Canonical Potential $\Omega(\mu,T;M)=-\frac{T}{L}\log\mathcal{Z}$ under the
condition of thermal equilibrium. It is worth noting that in this context, the
“volume” refers to a one-dimensional space that can be represented by a length
$L$. The partition function $\mathcal{Z}$, expressed in the path integral
formulation, is given by [KG06]

$\mathcal{Z}=\int\mathcal{D}\psi\int\mathcal{D}\bar{\psi}\exp\left[\int_{0}^{\beta}d\tau\int
dx\mathcal{L}_{E}\right]\,.$ (51)

Here Euclidean Lagrangian $\mathcal{L}_{E}$ is derived through a
transformation of variables to imaginary time, with $\tau=it$, while the
coldness $\beta=1/T$. Notably, the effective Lagrangian (section 3) leads to
$\mathcal{L}_{\mathrm{eff},E}=\mathcal{L}_{\mathrm{Dirac},E}-\mathcal{V}$,
which implies…

$\displaystyle\mathcal{Z}$
$\displaystyle=\int\mathcal{D}\psi\int\mathcal{D}\bar{\psi}\exp\left[\int_{0}^{\beta}d\tau\int
dx\,\mathcal{L}_{\mathrm{eff},E}\right]$
$\displaystyle=\int\mathcal{D}\psi\int\mathcal{D}\bar{\psi}\exp\left[\int_{0}^{\beta}d\tau\int
dx\,(\mathcal{L}_{\mathrm{Dirac},E}-\mathcal{V})\right]$
$\displaystyle=e^{-\beta
L\mathcal{V}}\int\mathcal{D}\psi\int\mathcal{D}\bar{\psi}\exp\left[\int_{0}^{\beta}d\tau\int
dx\,\mathcal{L}_{\mathrm{Dirac},E}\right]$
$\displaystyle=e^{-L\mathcal{V}/T}\mathcal{Z}_{\mathrm{Dirac}}$
$\displaystyle=e^{-\frac{L}{T}(\Omega_{\mathrm{Dirac}}+\mathcal{V})}\,.$ (52)

Given by the above derivation, one obtains the Grand Canonical Potential
$\Omega=\Omega_{\mathrm{Dirac}}+\mathcal{V}$, with $\Omega_{\mathrm{Dirac}}$
representing the Grand Canonical Potential of a free Dirac field of mass $M$
in $(1+1)$ dimensions [KG06, Bub05]:

$\displaystyle\Omega_{\mathrm{Dirac}}(\mu,T;M)=-\frac{2}{\pi}\int_{0}^{\infty}dk\left[\omega_{k}+T\log(1+e^{-\beta(\omega_{k}+\mu)})+T\log(1+e^{-\beta(\omega_{k}-\mu)})\right]\,.$
(53)

Here $\omega_{k}=\sqrt{k^{2}+M^{2}}$. Next, by defining the potential
$\mathcal{V}=(M-m)^{2}/4g$, one has

$\displaystyle\Omega(\mu,T;M)=\frac{(M-m)^{2}}{4g}-\frac{2}{\pi}\int_{0}^{\infty}dk\left[\omega_{k}+T\log(1+e^{-\beta(\omega_{k}+\mu)})+T\log(1+e^{-\beta(\omega_{k}-\mu)})\right]\,.$
(54)

This quantity displays a divergent behavior, which requires the implementation
of a regularization procedure. To facilitate a comparative analysis with the
Lattice model detailed in section 4 below, we naturally introduce a hard
cutoff $\Lambda=\pi/a$ for the upper bound in the momentum integration, where
$a$ represents the lattice spacing. With this explicit finite cutoff
introduced, we left the results un-removed by renormalization and dependent on
the lattice spacing. As a result, one can compute the effective mass $M$ for
fixed values of $\mu$ and $T$ by minimizing the Grand Canonical Potential
$\Omega(\mu,T;M)$ given by eq. 54 with respect to $M$. Subsequently, the
chiral condensate can be obtained from
$\langle\bar{\psi}\psi\rangle=(m-M)/2g$. For our simulation, we adopted the
parameter values $m=100$ MeV and $g=a=1$ MeV-1, which are well-suited for our
calculations.

In particular, we analyzed $\mu$ and $T$ over the interval $[0\
\mathrm{MeV},300\ \mathrm{MeV}]$ and present the resulting mass surface in
fig. 5 below.

Figure 5: The mass of the fermion field, expressed as an effective quantity
$M$, is depicted as a three-dimensional surface that varies with the
baryochemical potential $\mu$ and temperature $T$. In the accompanying contour
plot, the contours of constant $M$ are shown to provide a more detailed
representation of the mass landscape.

The findings indicate that when $\mu^{2}+T^{2}$ is sufficiently low, a
dynamically induced mass of approximately $\Delta m\equiv M-m=4$ MeV emerges
in addition to the mass $m=100$ MeV present in the Lagrangian. Importantly,
$M$ approaches $m$ asymptotically, as expected.

It is worth mentioning that historically, the GN model has primarily been
studied due to its property of being asymptotically free [GN74]. Therefore, at
extremely high temperatures or chemical potentials, it is anticipated that a
free field theory will be restored.

###### Study Chirality Imbalance with Quantum Simulations

Another very interesting feature for phase diagram is the significance of a
chiral chemical potential $\mu_{5}$. The existence of $\mu_{5}$ has shown
important consequences at the so-called chiral magnetic effect [KZ07, KMW08,
FKW08] in heavy ion collisions, which has been under intensive experimental
investigation [Abd22a, ZJ22]. In this section, we focus on studying the
effects of phase transition at finite $T$, $\mu$, and $\mu_{5}$.

The examination of the chiral chemical potential $\mu_{5}$ effects on the QCD
chiral phase transition involves studying the QCD phase diagram and the
behavior of the total chirality charge $N_{5}$ under the influence of an
external magnetic field at finite temperature $T$ and baryon chemical
potential $\mu$. To capture the effects of topological charge changing
transitions, a finite chiral chemical potential $\mu_{5}$ is introduced. In
the high temperature/density regime, a chirality charge may arise due to
fermion helicity flipping during the deconfinement phase transition from
hadronic matter to quark-gluon plasma. Experimental observations indicate that
the chirality charge reaches an equilibrium value shortly after a heavy-ion
collision, making it essential to explore the chiral imbalance in the QCD
phase diagrams for an accurate description of heavy-ion collisions.

The Nambu-Jona-Lasinio (NJL) model has been a prominent choice for studying
the chiral magnetic effect and the QCD chiral phase transition [NJ61b, NJ61c,
Fuk08, CdR07, LDC15, DLX15, CHS15, DCX13, SYX15, MPA09b, MPA09a, GMC07,
BBH12]. As an effective model for QCD, it allows for analytical calculations
at finite temperature $T$, chemical potentials $\mu$, and $\mu_{5}$.

###### Effective mass

Next we present theoretical calculations for the vacuum chiral condensate
$\langle\bar{\psi}\psi\rangle$ at various temperatures $T$, chemical
potentials $\mu$, and $\mu_{5}$. The chiral condensate
$\langle\bar{\psi}\psi\rangle$, which is a well-known order parameter [GPY81,
NJ61a, MS81, Pol78, FWZ18] for the chiral phase transition in the chiral limit
($m\sim 0$), has been studied in the mean field approximation. To achieve
this, as first introduced in [GN74, Wal74], we define
$\bar{\psi}\psi=\langle\bar{\psi}\psi\rangle+\sigma$ with a constant term
$\langle\bar{\psi}\psi\rangle$ and a small real scalar field $\sigma$, which
corresponds to fluctuations about the vacuum value. Then, we drop terms that
are $\mathcal{O}(\sigma^{2})$. By doing this, the four-fermion contact
interaction $g(\bar{\psi}\psi)^{2}$ can be expressed as

$\displaystyle g(\bar{\psi}\psi)^{2}=$ $\displaystyle
g\left(\langle\bar{\psi}\psi\rangle+\sigma\right)^{2}=2g\bar{\psi}\psi\langle\bar{\psi}\psi\rangle-g\langle\bar{\psi}\psi\rangle^{2}+\mathcal{O}(\sigma^{2})\,.$
(55)

Moreover, one can define the effective mass $M(\mu,\mu_{5},T)$ through the
following expression:

$\displaystyle
M(\mu,\mu_{5},T)=m-2g\langle\bar{\psi}\psi\rangle(\mu,\mu_{5},T),$ (56)

with the Lagrangian is defined as
$\mathcal{L}=\mathcal{L}_{\mathrm{eff}}+\mathcal{O}(\sigma^{2})$, and

$\displaystyle\mathcal{L}_{\mathrm{eff}}=$
$\displaystyle\bar{\psi}(i\not{\partial}-M+\mu\gamma_{0}+\mu_{5}\gamma_{0}\gamma_{5})\psi-\frac{(M-m)^{2}}{4g}$
$\displaystyle=$ $\displaystyle\mathcal{L}_{\mathrm{Dirac}}-\mathcal{V}\,.$
(57)

In accordance with previous studies [KG06, Bub05], we express the Grand
Canonical Potential $\Omega_{\mathrm{Dirac}}$ of
$\mathcal{L}_{\mathrm{Dirac}}$, which has a mass $M$, as follows:

$\displaystyle\Omega_{\mathrm{Dirac}}(\mu,\mu_{5},T;M)=-\frac{2}{\pi}\sum_{s=\pm
1}\int_{0}^{\infty}\bigg{[}T\ln(1+e^{-\beta(\omega_{k,s}+\mu)})+T\ln(1+e^{-\beta(\omega_{k,s}-\mu)})+\omega_{k,s}\bigg{]}dk\,,$
(58)

Here, the energy spectrum of the free fermions
$\omega_{k,s}=\sqrt{(k+s\mu_{5})^{2}+M^{2}}$ with $s=\pm 1$. To obtain the
grand canonical potential for the NJL model, we add the potential
$\mathcal{V}=(M-m)^{2}/4g$ to the expression.

$\displaystyle\Omega(\mu,\mu_{5},T;M)=\mathcal{V}-\frac{2}{\pi}\sum_{s}\int_{0}^{\infty}\bigg{[}T\ln(1+e^{-\beta(\omega_{k,s}+\mu)})+T\ln(1+e^{-\beta(\omega_{k,s}-\mu)})+\omega_{k,s}\bigg{]}dk\,.$
(59)

To obtain the value of the effective mass, the regularization of the grand
canonical potential $\Omega(\mu,\mu_{5},T;M)$ is necessary to account for its
divergent behavior. For consistency with numerical simulations using a lattice
spacing $a$, a natural hard momentum cutoff $\Lambda=\pi/a$ is adopted for the
integral in eq. 59. Subsequently, the effective mass $M$ is determined at
fixed $\mu,\mu_{5}$, and $T$ through numerical minimization of
$\Omega(\mu,\mu_{5},T;M)$ with respect to $M$, which amounts to solving the
gap equation.

$\displaystyle\frac{\partial\Omega(\mu,\mu_{5},T;M)}{\partial M}=0\,.$ (60)

This leads to the expression for the chiral condensate as
$\langle\bar{\psi}\psi\rangle=(m-M)/2g$ as given in eq. 56.

In the presence of a magnetic field and under an imbalance of right/left-
handed chirality, a finite induced current is generated along the magnetic
field direction. This phenomenon is known as the Chiral Magnetic Effect [KZ07,
KMW08, FKW08]. Specifically, when the number of right-handed quarks ($N_{R}$)
is not equal to that of left-handed quarks ($N_{L}$), a separation of positive
and negative charges takes place along the magnetic field. The CME arises from
the axial anomaly and topological objects in QCD and can produce observable
effects that are useful for studying topological $\mathcal{P}$\- and
$\mathcal{CP}$-odd excitations [Wit79, Ven79, SS98, VP09].

###### Chirality charge density

In the Lagrangian equation eq. 47, we incorporate a novel term featuring the
chiral chemical potential $\mu_{5}$ coupled with the chirality charge density
operator $n_{5}=\bar{\psi}\gamma_{0}\gamma_{5}\psi$. This operator is a
distinguishing attribute of hot and dense QCD matter since it violates the
conservation law owing to the chiral anomaly. Although $\mu_{5}$ is introduced
to examine topological charge fluctuations, it is considered a time-
independent variable representing the chirality imbalance. Similar to the
chiral condensate $\langle\bar{\psi}\psi\rangle$, the chirality charge density
$n_{5}=\langle\bar{\psi}\gamma_{0}\gamma_{5}\psi\rangle$ is also a global
quantity that does not depend on the coordinates. To relate the induced
electric current density with the chirality density, the correlation between
$n_{5}$ and $\mu_{5}$ is crucial, and the chirality charge density $n_{5}$ can
be computed using the method proposed in [FRG10].

$\displaystyle n_{5}=$
$\displaystyle-\frac{\partial\Omega(\mu,\mu_{5},T;M)}{\partial\mu_{5}}\,,$
(61)

The expression for the grand potential $\Omega(\mu,\mu_{5},T;M)$ is provided
in Eqn. eq. 59. In the subsequent section, we will furnish the analytical
computations of $n_{5}$, which will be contrasted with the results obtained
using QITE and exact diagonalization.

#####  4 Discritization of NJL Hamiltonian

In this subsection, we present the constituent parts of the quantum algorithms
employed for the NJL model. Specifically, we begin by outlining the
discretization of the NJL Hamiltonian using a lattice grid, before introducing
the quantum imaginary time evolution algorithm.

In order to investigate the chiral phase transition and chirality imbalance in
the NJL model, we incorporate terms associated with non-zero chemical
potential $\mu$ and chiral chemical potential $\mu_{5}$. These terms simulate
the chiral asymmetry between right- and left-handed quarks coupled with the
chirality charge density operator $n_{5}=\bar{\psi}\gamma_{0}\gamma_{5}\psi$.
Hence, the revised Lagrangian is given by111In our work [CKM22], we explored
the characteristics of the chiral condensate $\langle\bar{\psi}\psi\rangle$
under the conditions of finite and non-zero temperature $T$ and chemical
potential $\mu$, while setting $\mu_{5}=0$ and in [CKT22] we included finite
chiral potential in the study.:

$\displaystyle\mathcal{L}=$
$\displaystyle\bar{\psi}(i\not{\partial}-m)\psi+g(\bar{\psi}\psi)^{2}+\mu\bar{\psi}\gamma_{0}\psi+\mu_{5}\bar{\psi}\gamma_{0}\gamma_{5}\psi\,.$
(62)

Then the corresponding Hamiltonian
$\mathcal{H}=i\psi^{\dagger}\partial_{0}\psi-\mathcal{L}$ can be expressed as
a function of eq. 62.

$\displaystyle\mathcal{H}=$
$\displaystyle\bar{\psi}(i\gamma_{1}{\partial_{1}}+m)\psi-g(\bar{\psi}\psi)^{2}-\mu\bar{\psi}\gamma_{0}\psi-\mu_{5}\bar{\psi}\gamma_{0}\gamma_{5}\psi\,.$
(63)

To avoid ambiguity, we specify that in this section, the term “NJL model” is
employed to refer to the Hamiltonian as stated in Equation eq. 63.

For the discretization of a theory of the Dirac fermion field $\psi(x)$ with
components $\rho(x)$ and $\eta(x)$, we set a staggered fermion field
$\chi_{n}$ on a spatial lattice of spacing $a$. Using an even integer $N$ and
integer values $n=0,1,\cdots N/2$, the sites with even indices $2n$ are
assigned to the staggered fermion field $\chi_{2n}$ with spinor
$\rho(x=n)/\sqrt{a}$, while sites with odd indices $2n+1$ are assigned to the
staggered fermion field $\chi_{2n+1}$ with spinor $\eta(x=n)/\sqrt{a}$. The
Dirac fermion fields $\psi(x)$ are then represented by these staggered fermion
fields [BEF10, BFH14, AFK06, Baz14, ABT20],

$\displaystyle\psi(x=n)=\left(\begin{array}[]{cc}&\rho(x=n)/\sqrt{a}\\\
&\eta(x=n)/\sqrt{a}\end{array}\right)=\left(\begin{array}[]{cc}&\chi_{2n}\\\
&\chi_{2n+1}\end{array}\right)\,.$ (68)

Once one converts the Dirac fermions to staggered fermions, by applying the
Jordan-Wigner transformation [JW28],

$\displaystyle\chi_{n}=\frac{X_{n}-iY_{n}}{2}\prod_{i=0}^{n-1}\left(-iZ_{i}\right)\,,$
(69)

After the application of the Jordan-Wigner transformation to convert Dirac
fermions into staggered fermions, the spin representation of the operators is
obtained for quantum simulation. In eq. 69, the symbol $X_{n}$ represents the
Pauli-$X$ matrix operating on the $n$-th grid on the lattice, and so on.
Consequently, the following equivalence up to a constant term can be readily
confirmed,

$\displaystyle\int
dx\bar{\psi}\psi=\sum_{n=0}^{N-1}(-1)^{n}\chi_{n}^{\dagger}\chi_{n}=\sum_{n=0}^{N-1}(-1)^{n}\frac{Z_{n}}{2}\,,$
(70) $\displaystyle\int
dx\bar{\psi}\gamma_{0}\psi=\sum_{n=0}^{N-1}\chi_{n}^{\dagger}\chi_{n}=\sum_{n=0}^{N-1}\frac{Z_{n}}{2}\,,$
(71) $\displaystyle\int
dx\bar{\psi}i\gamma_{1}\partial_{1}\psi=-\frac{i}{2a}\sum_{n=0}^{N-1}\left[\chi_{n}^{\dagger}\chi_{n+1}-\chi_{n+1}^{\dagger}\chi_{n}\right]=\sum_{n=0}^{N-1}\frac{1}{4a}\left(X_{n}X_{n+1}+Y_{n}Y_{n+1}\right)\,,$
(72) $\displaystyle\int
dx(\bar{\psi}\psi)^{2}=\frac{1}{a}\left[\sum_{n=0}^{N-1}(-1)^{n}\chi_{n}^{\dagger}\chi_{n}\right]^{2}=\sum_{n,m=0}^{N-1}(-1)^{n+m}\frac{Z_{n}Z_{m}}{4a}\,,$
(73) $\displaystyle\int
dx\bar{\psi}\gamma_{0}\gamma_{5}\psi=-\sum_{n=0}^{N/2-1}\left(\chi_{2n}^{\dagger}\chi_{2n+1}+\chi_{2n+1}^{\dagger}\chi_{2n}\right)=\frac{1}{2}\sum_{n=0}^{N/2-1}\left(X_{2n}Y_{2n+1}-Y_{2n}X_{2n+1}\right)\,.$
(74)

where the subscripts denote the index of qubit upon which the single-qubit
gate is applied. Subsequently, the Hamiltonian in spin representation takes
the form of

$\displaystyle H=$ $\displaystyle\int
dx\left[\bar{\psi}(i\gamma_{1}\partial_{1}+m)\psi-g(\bar{\psi}\psi)^{2}-\mu\bar{\psi}\gamma_{0}\psi]\right]$
$\displaystyle=$
$\displaystyle-\frac{i}{2a}\sum_{n=0}^{N-1}\left[\chi_{n}^{\dagger}\chi_{n+1}-\chi_{n+1}^{\dagger}\chi_{n}\right]+m\sum_{n=0}^{N-1}(-1)^{n}\chi_{n}^{\dagger}\chi_{n}-\frac{g}{a}\left[\sum_{n=0}^{N-1}(-1)^{n}\chi_{n}^{\dagger}\chi_{n}\right]^{2}$
$\displaystyle-\mu\sum_{n=0}^{N-1}\chi_{n}^{\dagger}\chi_{n}+\mu_{5}\sum_{n=0}^{N/2-1}\left[\chi_{2n}^{\dagger}\chi_{2n+1}+\chi_{2n+1}^{\dagger}\chi_{2n}\right]\,.$
(75)

Note that hewe we choose a periodic boundary condition, namely
$\chi_{N}\to\chi_{0}$, and thus the $\chi_{N-1}$ field would be coupled with
the $\chi_{0}$ field. It is important to emphasize that in the continuum
limit, the true physics should, of course not depend, on the boundary
condition except for quantities that specifically pertain to the boundary.
However, since we choose a very small number of lattice sites $N=4$, it would
be important to study the dependence on the boundary condition. We leave such
a study in the future where we would also increase the number of lattice sites
in the simulation. Applying the Jordan-Wigner transformation to the
Hamiltonian given in eq. 75 yields the partitioning of $H$ into five distinct
components, which are necessary for the quantum algorithm construction.

$\displaystyle H_{1}$
$\displaystyle=\sum_{n=0}^{N/2-1}\frac{1}{4a}\left(X_{2n}X_{2n+1}+Y_{2n}Y_{2n+1}\right)\,,$
(76) $\displaystyle H_{2}$
$\displaystyle=\sum_{n=1}^{N/2-1}\frac{1}{4a}\left(X_{2n-1}X_{2n}+Y_{2n-1}Y_{2n}\right)+\frac{(-1)^{N/2}}{4a}\left(X_{N-1}X_{0}+Y_{N-1}Y_{0}\right)\prod_{i=1}^{N-2}Z_{i}\,,$
(77) $\displaystyle H_{3}$
$\displaystyle=\frac{m}{2}\sum_{n=0}^{N-1}(-1)^{n}Z_{n}\,,$ (78)
$\displaystyle H_{4}$
$\displaystyle=-\frac{g}{4a}\sum_{n,m=0}^{N-1}(-1)^{n+m}Z_{n}Z_{m}\,,$ (79)
$\displaystyle H_{5}$ $\displaystyle=-\frac{\mu}{2}\sum_{n=0}^{N-1}Z_{n}\,,$
(80) $\displaystyle H_{6}$
$\displaystyle=-\frac{\mu_{5}}{2}\sum_{n=0}^{N/2-1}\left(X_{2n}Y_{2n+1}-Y_{2n}X_{2n+1}\right)\,.$
(81)

Specifically, eq. 77 includes the second term to enforce periodic boundary
conditions. The resulting Hamiltonians outlined in eq. 76-eq. 81 can be
evolved using a quantum simulation through the application of the Suzuki-
Trotter decomposition [Tro59, Suz76], which enables us to examine the
influence of the chiral chemical potential $\mu_{5}$ on the chiral condensate
$\langle\bar{\psi}\psi\rangle$ and the chirality charge density $n_{5}$ of the
$(1+1)$-dimensional NJL model at finite temperature on a quantum simulator.

#####  5 Quantum imaginary time evolution algorithm

In the present subsection, we propose the use of the quantum imaginary time
evolution (QITE) algorithm [MST20] to evaluate the temperature dependence of
the NJL model for different values of the baryochemical potential $\mu$ and
chiral chemical potential $\mu_{5}$. As highlighted in [MST20], the QITE
algorithm possesses several advantages over other methods for quantum thermal
averaging procedures [TD00a, TOV11b, CS16, BK19]. Specifically, it enables the
computation of thermal averages of observables without any ancillary qubits or
deep circuits. In addition, the QITE algorithm is highly resource-efficient,
requiring exponentially less space and time for each iteration than classical
counterparts.

###### Trotterized evolution

In general, the Suzuki-Trotter decomposition [Tro59, Suz76] can be employed to
approximate the (Euclidean) evolution operator $e^{-\beta H}$ for a given
Hamiltonian $H$.

$\displaystyle e^{-\beta H}=\left(e^{-\Delta\beta
H}\right)^{N}+O(\Delta\beta^{2}),$ (82)

In this regard, the chosen imaginary time step $\Delta\beta$ determines the
number of iterations $N=\beta/\Delta\beta$ required to attain the desired
imaginary time $\beta=1/T$ at temperature $T$. Nevertheless, as the evolution
operator $e^{-\Delta\beta H}$ is not unitary, it is not feasible to realize it
as a sequence of unitary quantum gates. Therefore, to achieve the Euclidean
time evolution of a state $\ket{\Psi}$ on a quantum computer, an approximation
of the operator $e^{-\Delta\beta H}$ by some unitary operator is necessary.
Fortunately, the QITE algorithm offers a scheme to tackle this issue.

The QITE algorithm employs the introduction of a Hermitian operator $A$ to
approximate the non-unitary operator $e^{-\Delta\beta H}$ acting on a quantum
state $\ket{\Psi}$ with a unitary operator $e^{-i\Delta\beta A}$. It should be
noted that quantum states are represented by
${\alpha\ket{\Psi}:\alpha\in\mathbb{C}}$ in a Hilbert space, and not by the
vectors themselves, as the normalization and phase of the state vectors are
nonphysical.

$\displaystyle\frac{1}{\sqrt{c(\Delta\beta)}}e^{-\Delta\beta
H}\ket{\Psi}\approx e^{-i\Delta\beta A}\ket{\Psi}\,.$ (83)

Here the normalization is given by
$c(\Delta\beta)=\langle\Psi|e^{-2\Delta\beta H}\ket{\Psi}$.

In the regime where $\Delta\beta$ is very small, one may expand eq. 83 up to
the first nontrivial term, namely up to $\mathcal{O}(\Delta\beta)$, and
truncate the higher-order terms. By doing so, one can approximate the change
of quantum states under the operators $e^{-\Delta\beta H}$ and
$e^{-i\Delta\beta A}$ over a small imaginary time interval $\Delta\beta$ at
imaginary time $\beta$. This can be expressed as:

$\displaystyle\ket{\Delta\Psi_{H}(\beta)}$
$\displaystyle=\frac{1}{\Delta\beta}\left(\frac{1}{\sqrt{c(\Delta\beta)}}e^{-\Delta\beta
H}\ket{\Psi(\beta)}-\ket{\Psi(\beta)}\right)\,,$ (84)
$\displaystyle\ket{\Delta\Psi_{A}(\beta)}$
$\displaystyle=\frac{1}{\Delta\beta}\bigg{(}e^{-i\Delta\beta
A}\ket{\Psi(\beta)}-\ket{\Psi(\beta)}\bigg{)}\,,$ (85)

In accordance with the approach presented in [MST20], the Hermitian operator
$A$ can be obtained by parameterizing it in terms of Pauli matrices as
follows:

$\displaystyle A(\bm{a})=\sum_{\mu}a_{\mu}\hat{\sigma}{\mu},.$ (86)

The various Pauli strings are labeled by the subscript $\mu$, and the
corresponding Pauli string $\hat{\sigma}{\mu}=\prod{l}{\sigma}{\mu_{l},l}$ is
utilized. The evaluation of the Hermitian operator $A$ requires the
minimization of an objective function $F(a)$, where $a{\mu}$ denotes the
coefficients of the Pauli strings, and

$\displaystyle F(a)=$
$\displaystyle||\big{(}\ket{\Delta\Psi_{H}(\beta)}-\ket{\Delta_{A}\Psi(\beta)}\big{)}||^{2}$
(87) $\displaystyle=$
$\displaystyle||\,\ket{\Delta\Psi_{H}(\beta)}\,||^{2}+\sum_{\mu,\nu}a_{\nu}a_{\mu}\bra{\Psi(\beta)}\hat{\sigma}_{\nu}^{\dagger}\hat{\sigma}_{\mu}\ket{\Psi(\beta)}$
$\displaystyle+i\sum_{\mu}\frac{a_{\mu}}{\sqrt{c(\Delta\beta)}}\bra{\Psi(\beta)}\big{(}H\hat{\sigma}_{\mu}-\hat{\sigma}_{\mu}^{\dagger}H\big{)}\ket{\Psi(\beta)}\,.$

The first term $||\,\ket{\Delta\Psi_{H}(\beta)}\,||^{2}$ is independent of
$a_{\mu}$. Therefore, we can differentiate with respect to $a_{\mu}$ and set
the resulting expression equal to zero. This leads to a linear equation of the
form $({\bm{S}}+{\bm{S}}^{T})\,{\bm{a}}={\bm{b}}$, where the matrix ${\bm{S}}$
and vector ${\bm{b}}$ are defined as follows:

$\displaystyle S_{\mu\nu}$
$\displaystyle=\bra{\Psi(\beta)}\hat{\sigma}_{\nu}^{\dagger}\hat{\sigma}_{\mu}\ket{\Psi(\beta)}\,,$
(88) $\displaystyle b_{\mu}$
$\displaystyle=-\frac{i}{\sqrt{c(\Delta\beta)}}\bra{\Psi(\beta)}\big{(}H\hat{\sigma}_{\mu}+\hat{\sigma}_{\mu}^{\dagger}H\big{)}\ket{\Psi(\beta)}\,.$
(89)

Given the equation $({\bm{S}}+{\bm{S}}^{T})\,{\bm{a}}={\bm{b}}$ derived
earlier by setting the derivative of the objective function $F(a)$ to zero, we
can solve for $a_{\mu}$ and use it to evolve an initial quantum state to any
imaginary time $\beta$ by Trotterization. Specifically, we can apply the
Suzuki-Trotter formula to approximate the evolution operator $e^{-\beta H}$
and write it as a product of operators $e^{-i\Delta\beta A}$ and
$e^{-\Delta\beta H}$. By iteratively applying these operators, we can evolve
the initial quantum state to any desired imaginary time with high accuracy,
without the need for ancilla qubits or deep quantum circuits.

$\displaystyle\ket{\Psi(\beta)}$ $\displaystyle=\left(e^{-i\Delta\beta
A}\right)^{N}\ket{\Psi(0)}+\mathcal{O}(\Delta\beta)\,.$ (90)

In fig. 6, we present the implementation of the procedure for evolving a
quantum state under the unitary operator $e^{-i\Delta\beta A}\ $, where
$U_{j}|\Psi(\beta_{j})\rangle=e^{-i\Delta\beta
A}|\Psi(\beta_{j})\rangle=|\Psi(\beta_{j}+\Delta\beta)\rangle$. Specifically,
each block in the figure performs the evolution of the state
$\ket{\Psi(\beta_{j})}$ to $\ket{\Psi(\beta_{j}+\Delta\beta)}$. By repeating
this procedure for $n$ steps, we can obtain the state at the target total
evolution time $\beta$.

Figure 6: The scheme of quantum circuit for generating the thermal state
evolved from the initial state $\ket{\Psi}$ at imaginary time
$\beta=n\Delta\beta$. The dashed-line box contains $n$ blocks, each of which
evolves the state for a small imaginary time interval $\Delta\beta$. Together,
these blocks constitute the QITE algorithm for thermal state preparation.

The QMETTS algorithm [MST20] is employed to calculate the expectation value of
an observable $\hat{O}$ at finite temperature $T=1/\beta$ by a Markov-like
process that starts from a product state. The process involves the imaginary
time evolution of the state using QITE, measurement of $\hat{O}$, and
measurement in a product basis to collapse back to a product state.

Fig. 7 demonstrates the imaginary time evolution of the chiral condensate
$\langle\bar{\psi}\psi\rangle$ at chemical potentials $\mu=0$ MeV (left),
$100$ MeV (middle), and $\mu=150$ MeV (right) using small-time steps
$\Delta\beta=0.001$ MeV-1 (solid line) and $\Delta\beta=0.005$ MeV-1 (dashed
line). The QITE algorithm is applied to an initial equal superposition state
$|++++\rangle$ on a $(1+1)$-dimensional NJL model with coupling constant $g=1$
MeV-1, quark mass $m=100$ MeV, and lattice spacing $a=1$ MeV-1, followed by
exact diagonalization of the NJL Hamiltonian given by eq. 75 for reference.
Here we chose the number of qubits to be 4, i.e. $N=4$ in our simulation, and
thus each staggered fermion field locates on one qubit. Smaller evolution
steps result in more accurate quantum simulations. For the rest of this work,
we use $\Delta\beta=0.001$ MeV-1 for the QITE algorithm and compare the
effective mass $M$ among quantum simulation, exact diagonalization, and
theoretical analysis.

Figure 7: The chiral condensate $\langle\bar{\psi}\psi\rangle$ at chemical
potentials $\mu=0$ MeV, $\mu=100$ MeV, and $\mu=150$ MeV is simulated as a
function of the coldness parameter $\beta=1/T$ (MeV-1) for the
$(1+1)$-dimensional NJL model with a coupling constant of $g=1$ MeV-1, bare
quark mass $m=100$ MeV, and lattice spacing $a=1$ MeV-1. The simulation is
performed using the equal superposition state $|\Psi\rangle=\ket{++++}$ as the
initial state. The results obtained using the QITE algorithm with an imaginary
time step of $\Delta\beta=0.001$ MeV-1 and $\Delta\beta=0.005$ MeV-1 are shown
in solid and dashed lines, respectively. The exact diagonalization results are
presented as points for comparison.

###### Compare to other algorithms

In the ensuing section, we undertake a comparative study of the efficacy of
the QITE algorithm and the Variational Quantum Eigensolver (VQE) algorithms in
evaluating the ground-state energy of the $(1+1)$-dimensional NJL model
described earlier. As highlighted in [MST20], the QITE algorithm is notably
efficient for ground-state energy calculations. Fig. 8 portrays the ground-
state energy with bare mass $m=100$ MeV, chemical potentials $\mu=100$ MeV and
$\mu_{5}=10$ MeV at $g=1$, as a function of the number of operation steps
executed by the algorithms. The VQE algorithms have made considerable strides
in obtaining numerous notable outcomes on NISQ hardware [BVC22, JKG22, CYW21,
ONK22, JKG22], which have recently garnered significant attention. However,
the efficacy of the algorithm is limited by the reliance on an ansatz, as the
Hilbert space segment that the VQE algorithm can investigate is affected by
the specific variational ansatz employed, and the classical component of the
algorithm requires optimization. The QITE algorithm, in contrast, does not
necessitate an ansatz and systematically evolves the prepared state toward the
ground state after each time-step in a controlled manner. The state is
expected to converge to the ground state provided that the initial state
shares some overlap with it, with an error that can be accurately estimated.

As can be observed from the graph in fig. 8, the performance of the QITE and
Variational Quantum Eigensolver (VQE) algorithms are compared by plotting the
ground-state energy with varying operation steps for the $(1+1)$-dimensional
NJL model with bare mass $m=100$ MeV, chemical potentials $\mu=100$ MeV and
$\mu_{5}=10$ MeV at $g=1$. The QITE algorithm, implemented using the QFORTE
quantum algorithms library based on PYTHON, utilizes an imaginary time step of
$\Delta\beta=0.001$. On the other hand, VQE algorithm, with results obtained
using QISKIT of IBMQ, implements various optimizers for comparison with QITE
simulations.

It is worth noting that the QITE algorithm achieves the ground-state energy
with higher accuracy in fewer operation steps than the VQE implementations
shown in the plot. The performance of the VQE algorithm is restricted by its
reliance on a variational ansatz, which limits its ability to explore the
entirety of the Hilbert space. Moreover, the classical optimization component
of the algorithm also affects its efficiency. Conversely, the QITE algorithm
evolves the prepared state closer to the ground state in a controlled manner
without requiring an ansatz. Provided the initial state has some overlap with
the ground state, the state should converge to the ground state with a
quantifiable error.

One can also observe that the error of the VQE implementations begins to level
off at around $1\%$, owing to the finite set of the variational ansatz
scanning distance from the true vacuum. The graph in fig. 8 shows $500$
operation steps for both QITE and VQE algorithms, where the former is
represented by blue points with a curve, and the latter by light-blue, green
and orange curves for various optimizers.

Figure 8: Comparative analysis between the QITE and VQE algorithms utilizing
various optimizers, namely QOBYLA, L_BFGS_B, and SLSQP for $g=1$ case. To
investigate the efficacy of these algorithms in obtaining the ground-state
energy of the NJL Hamiltonian, we employ a 4-qubit setup on $N=2$ lattice
sites with $m=100$ MeV, $\mu=100$ MeV, and $\mu_{5}=10$ MeV.

Notably, in fig. 8 we plot the number of optimization steps for the VQE
algorithm and the number of thermal evolution steps $x=\beta/\Delta\beta$
($\Delta\beta=0.001$) for the QITE algorithms on the $x$-axis. We implement
the quantum circuit of the QITE algorithm in QFORTE, a Python-based quantum
algorithms library. Conversely, we utilize QISKIT of IBMQ for the VQE
algorithm, wherein the maximum optimization steps for all optimizers is set as
500. Our results, as depicted in fig. 8, highlight that the QITE algorithm
consistently achieves a higher degree of accuracy with fewer operation steps,
thus outperforming the VQE algorithm with various optimizers. In contrast, the
error in the VQE implementations seems to level off at around $1\%$,
indicating that the variational ansatz scans a finite set of possibilities
that are at a finite distance from the true vacuum.

In addition to its utility in ground state energy calculations, the QITE
algorithm also exhibits great potential in the simulation of thermal processes
at different temperatures. In this study, we aim to employ the QITE algorithm
introduced earlier to generate the thermal state $|\Psi(\beta/2)\rangle$, and
subsequently compute the thermal average of an observable $\hat{O}$ using the
following relation:

$\displaystyle\langle\hat{O}\rangle_{\beta}=\frac{{\textrm{T}r}(e^{-\beta\hat{H}}\hat{O})}{{\textrm{T}r}(e^{-\beta\hat{H}})}=\frac{\sum_{i\in\mathcal{S}}\bra{i}e^{-\beta\hat{H}/2}\hat{O}e^{-\beta\hat{H}/2}\ket{i}}{\sum_{i\in\mathcal{S}}\bra{i}e^{-\beta\hat{H}}\ket{i}}\,.$
(91)

In this context, $\mathcal{S}$ constitutes a full set of basis states for the
ground state [MST20]. In this study, we will focus on the computation of the
thermal mean of two observables, namely the chiral condensate
$\langle\bar{\psi}\psi\rangle$ and the chirality charge density
$n_{5}=\langle\bar{\psi}\gamma_{0}\gamma_{5}\psi\rangle$, by selecting
$\hat{O}=\bar{\psi}\psi$ and $\hat{O}=\bar{\psi}\gamma_{0}\gamma_{5}\psi$,
respectively.

#####  6 Results

In the present subsection, we present a demonstration of the NJL chiral phase
transition in $(1+1)$ dimensions through the depiction of the quark condensate
$\langle\bar{\psi}\psi\rangle$ as a function of temperature $T$ and chemical
potential $\mu$. The adoption of the $(1+1)$-dimensional model facilitates us
to derive analytical results by solving the gap equation as explained in
section 3, thereby allowing for a comparison with quantum simulations.
Moreover, the constraint imposed by the $(1+1)$ dimensionality enables the
design of a quantum circuit that employs a relatively small number of qubits,
rendering it feasible for implementation on presently available hardware for
further research purposes.

Our computations employ the following set of parameters: lattice spacing $a=1$
MeV-1, bare mass $m=100$ MeV, and coupling constant $g=1$ MeV-1. Our outcomes
have been derived via three distinct avenues:

  1. 1.

Simulation of thermal states through the QITE algorithm;

  2. 2.

Exact diagonalization of the Hamiltonian in spin representation;

  3. 3.

Numerical solution of the gap equation, leading to an analytical calculation.

###### Chiral phase diagram at $\mu_{5}=0$

We first explore the chiral condensate $\langle\bar{\psi}\psi\rangle$ as a
function of temperature $T$ and chemical potential $\mu$ in the context of the
$(1+1)$-dimensional NJL chiral phase transition with $\mu_{5}=0$. To achieve
this, we employ the QITE algorithm with the QForte package, and modify it for
calculating the properties of the NJL Hamiltonian in quantum simulation. The
initial state is set to be the equal superposition state $|+\cdots+\rangle$,
and the algorithm is utilized to compute the quark condensate at different
values of $\beta=1/T$. The simulation results are presented alongside the
theoretical calculations and exact diagonalization for comparison.
Specifically, the diamond data points represent the results obtained from the
QITE simulation, the solid curves indicate the theoretical calculation
described in section 3, and the dashed curves denote the exact diagonalization
results based on the matrix form of the discretized NJL Hamiltonian given by
Eqs.eq. 76 \- eq. 80. Our findings indicate that the quantum simulations align
well with both the exact diagonalization and theoretical analyses.

In fig. 9, we show the temperature dependence of the chiral condensate
$\langle\bar{\psi}\psi\rangle$ for different values of the chemical potential
$\mu$, ranging from 0 to 200 MeV. It is noteworthy that the quantum
simulations are carried out with a lattice spacing of $a=1$ MeV-1, a bare mass
of $m=100$ MeV, and a coupling constant of $g=1$ MeV-1. At chemical potentials
below 100 MeV, the quark condensate increases with temperature. Conversely,
for large chemical potentials, the quark condensate exhibits a contrasting
behavior within a small temperature range before gradually increasing with
temperature, as revealed by the figure. Moreover, the quark condensate at
various finite chemical potentials $\mu$ tends to converge at high
temperatures.

Figure 9: Chiral Condensate $\langle\bar{\psi}\psi\rangle$ as a function of
temperature $T$ at finite chemical potentials $\mu$, ranging from 0 to 200
MeV. Dashed curves represent the results from exact diagonalization, while
solid curves correspond to analytical calculations.

fig. 10 depicts the temperature dependence of the quark condensate
$\langle\bar{\psi}\psi\rangle$ for various ratios of the chemical potential
and temperature, $\mu/T$, ranging from 0 to 8. The simulation is conducted
using the QITE algorithm and the same set of parameters as fig. 9. The diamond
data points represent the quantum simulation results, whereas the solid and
dashed lines represent the analytical calculation and exact diagonalization
results, respectively. At low temperatures, the quark condensate reaches a
minimum value of about $-2$ MeV. As the temperature increases, the quark
condensate gradually increases to approach zero. It is observed from the
curves that, for a larger ratio of $\mu/T$, the quark condensate increases
more rapidly from $-2$ to 0 MeV as the temperature increases.

Figure 10: This plot displays the temperature dependence of the Chiral
Condensate $\langle\bar{\psi}\psi\rangle$ for different ratios of the chemical
potential and temperature, $\mu/T$, ranging from 0 to 8. The diamond points
represent the results obtained from the QITE algorithm, while the solid curves
correspond to the analytical calculation and the dashed curves to the exact
diagonalization.

We present the temperature dependence of the quark condensate
$\langle\bar{\psi}\psi\rangle$ at various chemical potentials in the
$\langle\bar{\psi}\psi\rangle-\mu$ plane in fig. 11, where the temperatures
are chosen as $T=25,\ 50,\ 100,\ 125,\ 250$ MeV. We observe that the absolute
value of the quark condensate is larger at lower temperatures and chemical
potentials. Moreover, as the chemical potential $\mu$ increases, the quark
condensate experiences a rapid increase and approaches zero. Conversely, at
higher temperatures, the quark condensate remains nearly unchanged for
different chemical potentials. As such, the $\langle\bar{\psi}\psi\rangle-\mu$
plot exhibits a more pronounced phase transition at higher temperatures.

Figure 11: In the $\langle\bar{\psi}\psi\rangle-\mu$ plane, we plot the quark
condensate $\langle\bar{\psi}\psi\rangle$ at temperatures $T=25,\ 50,\ 100,\
125,\ 250$ MeV for various chemical potentials. The solid curves represent the
analytical calculation, while the dashed curves depict the results from exact
diagonalization.

In this study, we have successfully developed a quantum simulation for the
chiral phase transition of the 1+1 dimensional NJL model at finite temperature
and chemical potentials using the QITE algorithm. Specifically, we use a
4-qubit quantum circuit, i.e. each staggered fermion field is located in one
qubit. On one aspect, our work [CKM22] serve as a demonstration using a
limited number of lattice sites like the above examples at the frontier of
this field. At the same time, we are aware of efforts in the field of quantum
computing with larger lattice sizes ($>100$ sites) and actively engaging in
endeavors aimed at augmenting the lattice size to approach the continuum limit
and gain deeper insights into the NJL model.

Our results reveal a consistency between the digital quantum simulation and
the standard method for the lattice NJL Hamiltionian with a limited lattice
sites. Even though there is still a long way to go for comparing the quantum
simulation with the continuum limit of the NJL model, this small illustration
makes it promising for further development where the potential of quantum
computing in simulating finite-temperature behaviors for QCD could be
possible, and offers an exciting avenue for exploring finite density effects
in QCD and other field theories. With the promise of scalable quantum computer
technology on the horizon, this study and previous ones showcase the potential
for NISQ quantum computers to tackle physical problems that are challenging or
impossible to solve using classical computing algorithms.

###### Chiral phase diagram at $\mu_{5}\neq 0$

In this subsection, we investigate the thermal behavior of the effective mass
$M$ and chirality charge density
$n_{5}=\langle\bar{\psi}\gamma_{0}\gamma_{5}\psi\rangle$ in the
$(1+1)$-dimensional NJL model as given by eq. 63. As we have highlighted in
previous sections, our approach employs a quantum algorithm to simulate the
thermal properties of physical observables. To demonstrate the accuracy and
reliability of our simulations, we present the results with the same three
approaches as mentioned in the previous subsection. For consistency among the
three methods, we fix the bare mass $m=100$ MeV and lattice spacing $a=1$
MeV-1, and consider the effects of the four-fermion interaction term in the
Lagrangian by testing the coupling constant at $g=1$.

In order to implement quantum circuits, various quantum simulation packages
have been developed and have shown consistent results for quantum simulations.
These software libraries, including PYQUILL [SCZ16] (Rigetti), TEQUILA [al21],
Q# [qsh] (Microsoft), QISKIT [AAB19b] (IBM), QFORTE [SE21], XACC [MLD20], FQE
[RGW21], and CIRQ [Dev21], are coded in PYTHON, and generate expected outputs
of an ideal quantum computer. General quantum simulation packages can be found
in [Bha22], while specific implementations of quantum algorithms can be found
in [ASA21]. In our case, the QITE algorithm is executed using a quantum
circuit constructed with the open-source software package QFORTE [SE21], which
has implemented many useful quantum algorithms.

Figure 12: Effective mass $M$ as a function of temperature $T$ (left panels)
or $\mu_{5}$ (right panels) at coupling constant $g=1$. Each panel adopts a
different value of chemical potential $\mu$.

In this study, we investigate the effects of chiral imbalance on the chiral
condensate at finite temperatures and chemical potentials. We present several
plots of the effective mass $M$ at different chemical potentials $\mu$,
$\mu_{5}$, and temperatures $T$. In fig. 12, we compare the temperature (left
panels) and chiral chemical potentials $\mu_{5}$ (right panels) dependence of
the effective mass $M$ at coupling constant $g=1$ for various baryochemical
potentials $\mu$. The results are obtained using the QITE algorithm, exact
diagonalization, and analytical calculations. Notably, we observe a change in
the pattern of the effective mass between $\mu=100$ MeV and $\mu=150$ MeV
panels, from decreasing to increasing with temperature. Additionally, we find
that the effective mass at smaller $\mu_{5}$ is larger for $\mu\leq 100$ MeV
panels and smaller for $\mu\geq 150$ MeV panels. We also observe that at
$\mu\leq 100$ MeV, $M$ is lower for higher temperature, as expected by
asymptotic freedom, while at higher $\mu$, the effective mass $M$ becomes
smaller at lower temperatures.

Figure 13: Effective mass $M$ as a function of temperature $T$ MeV at coupling
constant $g=1$ at a fixed chiral chemical potential $\mu_{5}$ in each panel.

In fig. 13, we fix the value of chiral chemical potentials $\mu_{5}$ and plot
the $M\mathrm{-}T$ curves for various baryochemical potentials $\mu$ at $g=1$.
We observe a non-trivial phase transition at $100<\mu<150$ MeV in all panels.
The pattern of the $M\mathrm{-}T$ curves looks similar at various $\mu_{5}$,
with the effective mass $M$ changing more rapidly as a function of $T$ at
smaller $\mu_{5}$.

Figure 14: Chirality charge density $n_{5}$ as a function of temperature $T$
(left panels) or $\mu_{5}$ (right panels) at coupling constant $g=1$. Each
panel adopts a different value of chemical potential $\mu$.

In fig. 14, we present the chirality charge density as a function of
temperature $T$ (left panels) or as a function of chiral chemical potential
$\mu_{5}$ (right panels) for coupling constants $g=1$ and $g=5$. We adopt
distinct hues to differentiate various chiral chemical potentials $\mu_{5}=0,\
10,\ 30$ and $50$ MeV, while fixing the chemical potential $\mu$ at $0,\ 50,\
100,\cdots$ or $250$ MeV for each panel. In the absence of any other mechanism
for generating a non-zero
$\mathcal{N}_{5}=\bar{\psi}\gamma_{0}\gamma_{5}\psi=\psi^{\dagger}_{R}\psi_{R}-\psi^{\dagger}_{L}\psi_{L}$,
$n_{5}=0$ at $\mu_{5}=0$ MeV. Thus, $n_{5}=\langle\mathcal{N}_{5}\rangle$ only
exists when $\mu_{5}\neq 0$.

We also observe that, across all chemical potentials $\mu$ and temperatures
$T$, the chirality charge density $n_{5}$ is approximately proportional to the
chiral chemical potential $\mu_{5}$. Specifically, $n_{5}$ starts at zero at
$\mu_{5}=0$ MeV and increases linearly with $\mu_{5}$. Thus, in this model,
the chiral chemical potential serves as a direct measure of the chiral
imbalance present in the plasma. However, the rate of increase of $n_{5}$ with
$\mu_{5}$ decreases at higher chemical potentials $\mu$. Notably, a phase
transition can be discerned between $\mu\leq 100$ MeV and $\mu\geq 150$ MeV.
In the former case, $n_{5}$ decreases with temperature for each $\mu_{5}$,
while in the latter case, $n_{5}$ increases with temperature from $T=50$ MeV
to $T=100$ MeV before subsequently decreasing from $T=100$ MeV to $T=200$ MeV.

Figure 15: Chirality charge density $n_{5}$ as a function of temperature at
$g=1$ with a fixed chiral chemical potential $\mu_{5}$ in each panel.

Similarly to the findings from the effective mass plots in fig. 13, a non-
trivial phase transition emerges between $\mu=100$ MeV and $\mu=150$ MeV. For
$\mu\leq 100$ MeV, the chirality charge density decreases with increasing
temperature, while for $\mu\geq 150$ MeV, it first increases before declining
and ultimately approaching zero as temperature increases.

In fig. 15, we depict the chirality charge density $n_{5}$ as a function of
temperature $T$ for coupling constants $g=1$ and $g=5$. We utilize distinct
hues to distinguish various chemical potentials $\mu\in{0,\ 50,\cdots,250}$
MeV, while fixing the chiral chemical potential $\mu_{5}$ at $0,\ 10,\ 30$ or
$50$ MeV for each panel. As expected from the preceding figure, the chirality
charge density is $0$ for all temperatures and chemical potentials when
$\mu_{5}=0$ MeV. At non-zero $\mu_{5}$ values, we also observe the phase
transition between $\mu=100$ MeV and $\mu=150$ MeV. Specifically, for $\mu\leq
100$ MeV, the chirality charge density commences with some non-zero value at
$T=0$ MeV and diminishes as temperature increases, while for curves with
$\mu\geq 150$ MeV, the chirality charge density starts at $0$ at $T=0$ MeV and
initially increases before decreasing and converging to 0 with increasing
temperature. Furthermore, for curves with greater $\mu_{5}$ values, those with
$\mu\leq 100$ MeV start with higher chirality charge densities.

###  Chapter 6 Conclusion

In this PhD dissertation, we have explored the fascinating world of particle
physics, with a particular focus on Quantum Chromodynamics (QCD) and its
pivotal role in describing the strong force. Our investigation revolved around
two significant aspects of QCD studies, aiming to deepen our understanding of
elementary particles and their behavior.

The first aspect of our research delved into the exciting realm of
perturbative QCD methods for quantum 3D imaging of hadrons. Through QCD
factorization and perturbative expansions, we harnessed these powerful tools
to predict the behavior of high-energy interactions, specifically in inclusive
jet and hadron-in-jet production at colliders. This allowed us to gain
valuable insights into the intricate dynamics of quarks and gluons within
hadrons, shedding light on phenomena like hadron production inside jets and
spin asymmetries. This new concept “polarized jet fragmentation functions”
allow us to probe multi-dimensional structure of the nucleons and hadrons, an
important thrust for the future Electron-Ion Collider.

Simultaneously, we recognized the significance of non-perturbative QCD
studies, particularly in investigating the QCD phase diagram and understanding
the properties of hadrons. Embracing the advancements in quantum computing and
simulators, we explored the potential of enhancing finite-temperature behavior
simulations. This breakthrough opens up avenues to scrutinize extreme
temperatures and densities with unprecedented accuracy and detail, offering
fresh perspectives on the nature of QCD.

Throughout this research, we pursued two approaches to enrich our
comprehension of QCD. Firstly, we focused on investigating the nucleon
structure by analyzing polarized jet fragmentation functions. Secondly, we
harnessed quantum computing techniques to explore chiral phase transitions.
These enabled us to delve into the non-perturbative facets of QCD and discern
new applications for quantum computing in the realm of particle physics.

By undertaking these multifaceted studies, we have contributed to the ongoing
efforts to decipher the intricate interactions between elementary particles.
Our findings provide essential insights into the behavior of quarks and gluons
within hadrons, advancing our knowledge of particle physics. Moreover, the
integration of quantum computing methodologies demonstrates the immense
potential it holds in unraveling the mysteries of QCD and enhancing the
accuracy of simulations at extreme conditions.

As we conclude this thesis, we envision a promising future where continued
research and innovation in QCD studies, combined with quantum computing
advancements, will lead us towards new frontiers of knowledge. The pursuit of
understanding the fundamental building blocks of our universe remains an ever-
evolving and deeply rewarding journey, and we hope that this work serves as a
stepping stone towards further discoveries in the captivating realm of
particle physics.

### References

  * [Aab18] Morad Aaboud et al.  “Measurement of jet fragmentation in 5.02 TeV proton-lead and proton-proton collisions with the ATLAS detector.”  Nucl. Phys., A978:65, 2018.
  * [Aab19a] Morad Aaboud et al.  “Comparison of fragmentation functions for light-quark- and gluon-dominated jets from $pp$ and Pb+Pb collisions in ATLAS.”  Submitted to: Phys. Rev. Lett., 2019.
  * [AAB19b] Gadi Aleksandrowicz, Thomas Alexander, Panagiotis Barkoutsos, Luciano Bello, Yael Ben-Haim, David Bucher, Francisco Jose Cabrera-Hernández, Jorge Carballo-Franquis, Adrian Chen, Chun-Fu Chen, Jerry M. Chow, Antonio D. Córcoles-Gonzales, Abigail J. Cross, Andrew Cross, Juan Cruz-Benito, Chris Culver, Salvador De La Puente González, Enrique De La Torre, Delton Ding, Eugene Dumitrescu, Ivan Duran, Pieter Eendebak, Mark Everitt, Ismael Faro Sertage, Albert Frisch, Andreas Fuhrer, Jay Gambetta, Borja Godoy Gago, Juan Gomez-Mosquera, Donny Greenberg, Ikko Hamamura, Vojtech Havlicek, Joe Hellmers, Łukasz Herok, Hiroshi Horii, Shaohan Hu, Takashi Imamichi, Toshinari Itoko, Ali Javadi-Abhari, Naoki Kanazawa, Anton Karazeev, Kevin Krsulich, Peng Liu, Yang Luh, Yunho Maeng, Manoel Marques, Francisco Jose Martín-Fernández, Douglas T. McClure, David McKay, Srujan Meesala, Antonio Mezzacapo, Nikolaj Moll, Diego Moreda Rodríguez, Giacomo Nannicini, Paul Nation, Pauline Ollitrault, Lee James O’Riordan, Hanhee Paik, Jesús Pérez, Anna Phan, Marco Pistoia, Viktor Prutyanov, Max Reuter, Julia Rice, Abdón Rodríguez Davila, Raymond Harry Putra Rudy, Mingi Ryu, Ninad Sathaye, Chris Schnabel, Eddie Schoute, Kanav Setia, Yunong Shi, Adenilton Silva, Yukio Siraichi, Seyon Sivarajah, John A. Smolin, Mathias Soeken, Hitomi Takahashi, Ivano Tavernelli, Charles Taylor, Pete Taylour, Kenso Trabing, Matthew Treinish, Wes Turner, Desiree Vogt-Lee, Christophe Vuillot, Jonathan A. Wildstrom, Jessica Wilson, Erick Winston, Christopher Wood, Stephen Wood, Stefan Wörner, Ismail Yunus Akhalwaya, and Christa Zoufal.  “Qiskit: An Open-source Framework for Quantum Computing.”, January 2019.
  * [Aad11] Georges Aad et al.  “Measurement of the jet fragmentation function and transverse profile in proton-proton collisions at a center-of-mass energy of 7 TeV with the ATLAS detector.”  Eur. Phys. J., C71:1795, 2011.
  * [Aad12] Georges Aad et al.  “Measurement of $D^{*+/-}$ meson production in jets from pp collisions at sqrt(s) = 7 TeV with the ATLAS detector.”  Phys. Rev., D85:052005, 2012.
  * [Aad14] Georges Aad et al.  “Measurement of inclusive jet charged-particle fragmentation functions in Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV with the ATLAS detector.”  Phys. Lett., B739:320–342, 2014.
  * [Aad15] Georges Aad et al.  “Measurement of the transverse polarization of $\Lambda$ and $\bar{\Lambda}$ hyperons produced in proton-proton collisions at $\sqrt{s}=7$ TeV using the ATLAS detector.”  Phys. Rev., D91(3):032004, 2015.
  * [Aad19] Georges Aad et al.  “Properties of jet fragmentation using charged particles measured with the ATLAS detector in $pp$ collisions at $\sqrt{s}=13$ TeV.”  Phys. Rev. D, 100(5):052011, 2019.
  * [Aad20] Georges Aad et al.  “Measurement of the transverse momentum distribution of Drell–Yan lepton pairs in proton–proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector.”  Eur. Phys. J. C, 80(7):616, 2020.
  * [Aai16] Roel Aaij et al.  “Measurement of the forward Z boson production cross-section in pp collisions at $\sqrt{s}=13$ TeV.”  JHEP, 09:136, 2016.
  * [Aai17] Roel Aaij et al.  “Study of J/$\psi$ Production in Jets.”  Phys. Rev. Lett., 118(19):192001, 2017.
  * [Aai19] Roel Aaij et al.  “Measurement of charged hadron production in $Z$-tagged jets in proton-proton collisions at $\sqrt{s}=8$ TeV.”  2019\.
  * [Aal12] T. Aaltonen et al.  “Transverse momentum cross section of $e^{+}e^{-}$ pairs in the $Z$-boson region from $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV.”  Phys. Rev. D, 86:052010, 2012.
  * [Aba08] V. M. Abazov et al.  “Measurement of the shape of the boson transverse momentum distribution in $p\bar{p}\to Z/\gamma^{*}\to e^{+}e^{-}+X$ events produced at $\sqrt{s}$=1.96-TeV.”  Phys. Rev. Lett., 100:102002, 2008.
  * [Abd22a] Mohamed Abdallah et al.  “Search for the chiral magnetic effect with isobar collisions at $\sqrt{s_{NN}}$=200 GeV by the STAR Collaboration at the BNL Relativistic Heavy Ion Collider.”  Phys. Rev. C, 105(1):014901, 2022.
  * [Abd22b] R. Abdul Khalek et al.  “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report.”  Nucl. Phys. A, 1026:122447, 2022.
  * [Abe09] B. I. Abelev et al.  “Longitudinal Spin Transfer to Lambda and anti-Lambda Hyperons in Polarized Proton-Proton Collisions at s**(1/2) = 200-GeV.”  Phys. Rev., D80:111102, 2009.
  * [Abe12] B. Abelev et al.  “Neutral pion and $\eta$ meson production in proton-proton collisions at $\sqrt{s}=0.9$ TeV and $\sqrt{s}=7$ TeV.”  Phys. Lett. B, 717:162–172, 2012.
  * [ABH19] Andrei Alexandru, Paulo F. Bedaque, Siddhartha Harmalkar, Henry Lamm, Scott Lawrence, and Neill C. Warrington.  “Gluon Field Digitization for Quantum Computers.”  Physical Review D, 100(11):114501, December 2019.  arXiv: 1906.11213.
  * [ABT20] Christopher Aubin, Thomas Blum, Cheng Tu, Maarten Golterman, Chulwoo Jung, and Santiago Peris.  “Light quark vacuum polarization at the physical point and contribution to the muon $g-2$.”  Phys. Rev. D, 101(1):014503, 2020.
  * [Acc16] A. Accardi et al.  “Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all.”  Eur. Phys. J. A, 52(9):268, 2016.
  * [Ach19a] Shreyasi Acharya et al.  “Charged jet cross section and fragmentation in proton-proton collisions at $\sqrt{s}$ = 7 TeV.”  Phys. Rev. D, 99(1):012016, 2019.
  * [Ach19b] Shreyasi Acharya et al.  “Measurement of the production of charm jets tagged with D0 mesons in pp collisions at $\sqrt{\mathrm{s}}=7$ TeV.”  JHEP, 08:133, 2019.
  * [Ach23] P. Achenbach et al.  “The Present and Future of QCD.”  3 2023.
  * [Ada15] A. Adare et al.  “An Upgrade Proposal from the PHENIX Collaboration.”  1 2015.
  * [Ada18a] Jaroslav Adam et al.  “Improved measurement of the longitudinal spin transfer to $\Lambda$ and $\bar{\Lambda}$ hyperons in polarized proton-proton collisions at $\sqrt{s}$ = 200 GeV.”  Phys. Rev., D98(11):112009, 2018.
  * [Ada18b] L. Adamczyk et al.  “Transverse spin-dependent azimuthal correlations of charged pion pairs measured in p↑+p collisions at $\sqrt{s}$ = 500 GeV.”  Phys. Lett., B780:332–339, 2018.
  * [Ada18c] Leszek Adamczyk et al.  “Azimuthal transverse single-spin asymmetries of inclusive jets and charged pions within jets from polarized-proton collisions at $\sqrt{s}=500$ GeV.”  Phys. Rev., D97(3):032004, 2018.
  * [ADE05] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and K. Redlich.  “Thermodynamics of two flavor QCD to sixth order in quark chemical potential.”  Phys. Rev. D, 71:054508, 2005.
  * [AEF06] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo.  “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics.”  Nature, 443:675–678, 2006.
  * [AEH02] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, C. Schmidt, and L. Scorzato.  “The QCD thermal phase transition in the presence of a small chemical potential.”  Phys. Rev. D, 66:074507, 2002.
  * [AFG14] Francois Arleo, Michel Fontannaz, Jean-Philippe Guillet, and Chi Linh Nguyen.  “Probing fragmentation functions from same-side hadron-jet momentum correlations in p-p collisions.”  JHEP, 04:147, 2014.
  * [AFK06] Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo.  “The Equation of state in lattice QCD: With physical quark masses towards the continuum limit.”  JHEP, 01:089, 2006.
  * [AFL19] E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mantysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita.  “The electron–ion collider: assessing the energy dependence of key measurements.”  Rept. Prog. Phys., 82(2):024301, 2019.
  * [Agh18] M. Aghasyan et al.  “Transverse-momentum-dependent Multiplicities of Charged Hadrons in Muon-Deuteron Deep Inelastic Scattering.”  Phys. Rev. D, 97(3):032006, 2018.
  * [Aid19] C. Aidala et al.  “Measurements of $\mu\mu$ pairs from open heavy flavor and Drell-Yan in $p+p$ collisions at $\sqrt{s}=200$ GeV.”  Phys. Rev. D, 99(7):072003, 2019.
  * [Air13] A. Airapetian et al.  “Multiplicities of charged pions and kaons from semi-inclusive deep-inelastic scattering by the proton and the deuteron.”  Phys. Rev. D, 87:074029, 2013.
  * [AKK08] S. Albino, B. A. Kniehl, and G. Kramer.  “AKK Update: Improvements from New Theoretical Input and Experimental Data.”  Nucl. Phys., B803:42–104, 2008.
  * [AKP20] Miguel Arratia, Zhong-Bo Kang, Alexei Prokudin, and Felix Ringer.  “Jet-based measurements of Sivers and Collins asymmetries at the future electron-ion collider.”  Phys. Rev. D, 102(7):074015, 2020.
  * [AKS17] Daniele P. Anderle, Tom Kaufmann, Marco Stratmann, Felix Ringer, and Ivan Vitev.  “Using hadron-in-jet data in a global analysis of $D^{*}$ fragmentation functions.”  Phys. Rev., D96(3):034028, 2017.
  * [al13] R. Islam et. al.  “Emergence and Frustration of Magnetism with Variable-Range Interactions in a Quantum Simulator.”  Science, 340(6132):583–587, 2013.
  * [al21] Jakob S Kottmann et. al.  “TEQUILA: a platform for rapid development of quantum algorithms.”  Quantum Science and Technology, 6(2):024009, mar 2021\.
  * [AP77] Guido Altarelli and G. Parisi.  “Asymptotic Freedom in Parton Language.”  Nucl. Phys. B, 126:298–318, 1977.
  * [Apr15a] Ani Aprahamian et al.  “Reaching for the horizon: The 2015 long range plan for nuclear science.”  10 2015.
  * [Apr15b] Ani Aprahamian et al.  “Reaching for the horizon: The 2015 long range plan for nuclear science.”  10 2015.
  * [Arr21] Miguel Arratia.  Talk given at XXVIII International Workshop on Deep-Inelastic Scattering and Related Subjects, 2021.
  * [Aru20] Frank et. al. Arute.  “Hartree-Fock on a superconducting qubit quantum computer.”  Science, 369(6507):1084, Aug 2020.
  * [ASA21] Abhinav Anand, Philipp Schleich, Sumner Alperin-Lea, Phillip W. K. Jensen, Sukin Sim, Manuel Díaz-Tinoco, Jakob S. Kottmann, Matthias Degroote, Artur F. Izmaylov, and Alán Aspuru-Guzik.  “A Quantum Computing View on Unitary Coupled Cluster Theory.”  9 2021.
  * [ASR20] Miguel Arratia, Youqi Song, Felix Ringer, and Barbara V. Jacak.  “Jets as precision probes in electron-nucleus collisions at the future Electron-Ion Collider.”  Phys. Rev. C, 101(6):065204, 2020.
  * [ASY20] A. Avkhadiev, P. E. Shanahan, and R. D. Young.  “Accelerating Lattice Quantum Field Theory Calculations via Interpolator Optimization Using Noisy Intermediate-Scale Quantum Computing.”  Phys. Rev. Lett., 124:080501, Feb 2020.
  * [BAM17] A. Bermudez, G. Aarts, and M. Müller.  “Quantum Sensors for the Generating Functional of Interacting Quantum Field Theories.”  Physical Review X, 7(4):041012, October 2017.  Publisher: American Physical Society.
  * [Bas08] S. Basak et al.  “QCD equation of state at non-zero chemical potential.”  PoS, LATTICE2008:171, 2008.
  * [Baz14] A. Bazavov et al.  “Equation of state in ( 2+1 )-flavor QCD.”  Phys. Rev. D, 90:094503, 2014.
  * [Baz19] A. Bazavov et al.  “Chiral crossover in QCD at zero and non-zero chemical potentials.”  Phys. Lett. B, 795:15–21, 2019.
  * [Baz20] A. Bazavov et al.  “Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data.”  Phys. Rev. D, 101(7):074502, 2020.
  * [BBB22] Alessandro Bacchetta, Valerio Bertone, Chiara Bissolotti, Giuseppe Bozzi, Matteo Cerutti, Fulvio Piacenza, Marco Radici, and Andrea Signori.  “Unpolarized transverse momentum distributions from a global fit of Drell-Yan and semi-inclusive deep-inelastic scattering data.”  JHEP, 10:127, 2022.
  * [BBC13] Richard D Ball, Valerio Bertone, Stefano Carrazza, Christopher S Deans, Luigi Del Debbio, Stefano Forte, Alberto Guffanti, Nathan P Hartland, José I Latorre, Juan Rojo, et al.  “Parton distributions with LHC data.”  Nuclear Physics B, 867(2):244–289, 2013.
  * [BBD98] Claude W. Bernard, Tom Blum, Carleton E. Detar, Steven A. Gottlieb, Urs M. Heller, James E. Hetrick, Craig McNeile, K. Rummukainen, R. Sugar, and D. Toussaint.  “Continuum limit of lattice QCD with staggered quarks in the quenched approximation: A Critical role for the chiral extrapolation.”  Phys. Rev. Lett., 81:3087–3090, 1998.
  * [BBD13] D. Banerjee, M. Bögli, M. Dalmonte, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller.  “Atomic Quantum Simulation of $\mathbf{U}(N)$ and $\mathrm{SU}(N)$ Non-Abelian Lattice Gauge Theories.”  Physical Review Letters, 110(12):125303, March 2013\.  Publisher: American Physical Society.
  * [BBH08] D. Boer, C. J. Bomhof, D. S. Hwang, and P. J. Mulders.  “Spin asymmetries in jet-hyperon production at LHC.”  Phys. Lett., B659:127–136, 2008.
  * [BBH12] G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, and A. Schäfer.  “QCD quark condensate in external magnetic fields.”  Phys. Rev. D, 86:071502, Oct 2012.
  * [BBL95] Geoffrey T. Bodwin, Eric Braaten, and G. Peter Lepage.  “Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium.”  Phys. Rev. D, 51:1125–1171, 1995.  [Erratum: Phys.Rev.D 55, 5853 (1997)].
  * [BBM05] A. Bacchetta, C. J. Bomhof, P. J. Mulders, and F. Pijlman.  “Single spin asymmetries in hadron-hadron collisions.”  Phys. Rev. D, 72:034030, 2005.
  * [BBM20] Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  “Quantum Algorithms for Quantum Chemistry and Quantum Materials Science.”  Chemical Reviews, 120(22):12685–12717, Oct 2020.
  * [BCH21] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne.  “Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs.”  Eur. Phys. J. C, 81(4):341, 2021.
  * [BDD04] Alessandro Bacchetta, Umberto D’Alesio, Markus Diehl, and C. Andy Miller.  “Single-spin asymmetries: The Trento conventions.”  Phys. Rev. D, 70:117504, 2004.
  * [BDG07] Alessandro Bacchetta, Markus Diehl, Klaus Goeke, Andreas Metz, Piet J. Mulders, and Marc Schlegel.  “Semi-inclusive deep inelastic scattering at small transverse momentum.”  JHEP, 02:093, 2007.
  * [BDH16] Reggie Bain, Lin Dai, Andrew Hornig, Adam K. Leibovich, Yiannis Makris, and Thomas Mehen.  “Analytic and Monte Carlo Studies of Jets with Heavy Mesons and Quarkonia.”  JHEP, 06:121, 2016.
  * [BDL17] Reggie Bain, Lin Dai, Adam Leibovich, Yiannis Makris, and Thomas Mehen.  “NRQCD Confronts LHCb Data on Quarkonium Production within Jets.”  Phys. Rev. Lett., 119(3):032002, 2017.
  * [BDM11] Daniel Boer, Markus Diehl, Richard Milner, Raju Venugopalan, Werner Vogelsang, et al.  “Gluons and the quark sea at high energies: Distributions, polarization, tomography.”  2011\.
  * [BDM12] D. Banerjee, M. Dalmonte, M. Müller, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller.  “Atomic Quantum Simulation of Dynamical Gauge Fields coupled to Fermionic Matter: From String Breaking to Evolution after a Quench.”  Physical Review Letters, 109(17):175302, October 2012\.  arXiv: 1205.6366.
  * [BDP22] Alessandro Bacchetta, Filippo Delcarro, Cristian Pisano, and Marco Radici.  “The 3-dimensional distribution of quarks in momentum space.”  Phys. Lett. B, 827:136961, 2022.
  * [BEF10] Szabolcs Borsanyi, Gergely Endrodi, Zoltan Fodor, Antal Jakovac, Sandor D. Katz, Stefan Krieg, Claudia Ratti, and Kalman K. Szabo.  “The QCD equation of state with dynamical quarks.”  JHEP, 11:077, 2010.
  * [BEK20] Adam Bzdak, Shinichi Esumi, Volker Koch, Jinfeng Liao, Mikhail Stephanov, and Nu Xu.  “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan.”  Phys. Rept., 853:1–87, 2020.
  * [BFG13] Richard D Ball, Stefano Forte, Alberto Guffanti, Emanuele R Nocera, Giovanni Ridolfi, and Juan Rojo.  “Unbiased determination of polarized parton distributions and their uncertainties.”  Nuclear Physics B, 874(1):36–84, 2013.
  * [BFG20] Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Ruben Kara, Sandor D. Katz, Paolo Parotto, Attila Pasztor, Claudia Ratti, and Kalman K. Szabo.  “QCD Crossover at Finite Chemical Potential from Lattice Simulations.”  Phys. Rev. Lett., 125(5):052001, 2020.
  * [BFH14] Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo.  “Full result for the QCD equation of state with 2+1 flavors.”  Phys. Lett. B, 730:99–104, 2014.
  * [BFL00] Christian W. Bauer, Sean Fleming, and Michael E. Luke.  “Summing Sudakov logarithms in $B\to X_{s}\gamma$ in effective field theory.”  Phys.Rev., D63:014006, 2000.
  * [BFP01] Christian W. Bauer, Sean Fleming, Dan Pirjol, and Iain W. Stewart.  “An Effective field theory for collinear and soft gluons: Heavy to light decays.”  Phys.Rev., D63:114020, 2001.
  * [BFP20] Ignacio Borsa, Daniel de Florian, and Iván Pedron.  “Jet Production in Polarized Deep Inelastic Scattering at Next-to-Next-to-Leading Order.”  Phys. Rev. Lett., 125(8):082001, 2020.
  * [BFP22] Ignacio Borsa, Daniel de Florian, and Iván Pedron.  “The full set of polarized deep inelastic scattering structure functions at NNLO accuracy.”  Eur. Phys. J. C, 82(12):1167, 2022.
  * [BGM11] Daniel Boer, Leonard Gamberg, Bernhard Musch, and Alexei Prokudin.  “Bessel-Weighted Asymmetries in Semi Inclusive Deep Inelastic Scattering.”  JHEP, 10:021, 2011.
  * [Bha22] Kishor Bharti et al.  “Noisy intermediate-scale quantum algorithms.”  Rev. Mod. Phys., 94(1):015004, 2022.
  * [BIL18] Florian Burger, Ernst-Michael Ilgenfritz, Maria Paola Lombardo, and Anton Trunin.  “Chiral observables and topology in hot QCD with two families of quarks.”  Phys. Rev. D, 98:094501, Nov 2018.
  * [BJ93] M. Burkardt and R. L. Jaffe.  “Polarized q —¿ Lambda fragmentation functions from e+ e- —¿ Lambda + X.”  Phys. Rev. Lett., 70:2537–2540, 1993.
  * [BJN21] Christian W. Bauer, Wibe A. de Jong, Benjamin Nachman, and Davide Provasoli.  “Quantum Algorithm for High Energy Physics Simulations.”  Phys. Rev. Lett., 126(6):062001, 2021.
  * [BK19] Fernando G. S. L. Brandão and Michael J. Kastoryano.  “Finite Correlation Length Implies Efficient Preparation of Quantum Thermal States.”  Commun. Math. Phys., 365(1):1–16, 2019.
  * [BKL18] Maarten G. A. Buffing, Zhong-Bo Kang, Kyle Lee, and Xiaohui Liu.  “A transverse momentum dependent framework for back-to-back photon+jet production.”  2018\.
  * [BM00] A. Bacchetta and P. J. Mulders.  “Deep inelastic leptoproduction of spin-one hadrons.”  Phys. Rev., D62:114004, 2000.
  * [BMG19] Matthew J. S. Beach, Roger G. Melko, Tarun Grover, and Timothy H. Hsieh.  “Making trotters sprint: A variational imaginary time ansatz for quantum many-body systems.”  Phys. Rev. B, 100:094434, Sep 2019.
  * [BMM16] Reggie Bain, Yiannis Makris, and Thomas Mehen.  “Transverse Momentum Dependent Fragmenting Jet Functions with Applications to Quarkonium Production.”  JHEP, 11:144, 2016.
  * [BMP03] Daniel Boer, P. J. Mulders, and F. Pijlman.  “Universality of T-odd effects in single spin and azimuthal asymmetries.”  Nucl. Phys., B667:201–241, 2003.
  * [BMP04] C.J. Bomhof, P.J. Mulders, and F. Pijlman.  “Gauge link structure in quark-quark correlators in hard processes.”  Phys.Lett., B596:277–286, 2004.
  * [BMP10] Vincenzo Barone, Stefano Melis, and Alexei Prokudin.  “The Boer-Mulders effect in unpolarized SIDIS: An Analysis of the COMPASS and HERMES data on the cos 2 phi asymmetry.”  Phys. Rev. D, 81:114026, 2010.
  * [BMS21] Khadeejah Bepari, Sarah Malik, Michael Spannowsky, and Simon Williams.  “Towards a quantum computing algorithm for helicity amplitudes and parton showers.”  Physical Review D, 103(7), Apr 2021.
  * [BMT15] A. Bazavov, Y. Meurice, S.-W. Tsai, J. Unmuth-Yockey, and Jin Zhang.  “Gauge-invariant implementation of the Abelian-Higgs model on optical lattices.”  Physical Review D, 92(7):076003, October 2015.  Publisher: American Physical Society.
  * [BN11] Thomas Becher and Matthias Neubert.  “Drell-Yan Production at Small $q_{T}$, Transverse Parton Distributions and the Collinear Anomaly.”  Eur. Phys. J. C, 71:1665, 2011.
  * [BNP07] Thomas Becher, Matthias Neubert, and Ben D. Pecjak.  “Factorization and Momentum-Space Resummation in Deep-Inelastic Scattering.”  JHEP, 01:076, 2007.
  * [BNW12] Thomas Becher, Matthias Neubert, and Daniel Wilhelm.  “Electroweak Gauge-Boson Production at Small $q_{T}$: Infrared Safety from the Collinear Anomaly.”  JHEP, 02:124, 2012.
  * [BNW13] Thomas Becher, Matthias Neubert, and Daniel Wilhelm.  “Higgs-Boson Production at Small Transverse Momentum.”  JHEP, 05:110, 2013.
  * [Boe10] Daniel Boer.  “Transverse Lambda polarization at high energy colliders.”  PoS, DIS2010:215, 2010.
  * [Bou23] Renaud Boussarie et al.  “TMD Handbook.”  4 2023.
  * [BPS02] Christian W. Bauer, Dan Pirjol, and Iain W. Stewart.  “Soft collinear factorization in effective field theory.”  Phys.Rev., D65:054022, 2002.
  * [BPS16] Thomas Becher, Benjamin D. Pecjak, and Ding Yu Shao.  “Factorization for the light-jet mass and hemisphere soft function.”  JHEP, 12:018, 2016.
  * [BPV21a] Marcin Bury, Alexei Prokudin, and Alexey Vladimirov.  “Extraction of the Sivers function from SIDIS, Drell-Yan, and $W^{\pm}/Z$ boson production data with TMD evolution.”  JHEP, 05:151, 2021.
  * [BPV21b] Marcin Bury, Alexei Prokudin, and Alexey Vladimirov.  “Extraction of the Sivers Function from SIDIS, Drell-Yan, and $W^{\pm}/Z$ Data at Next-to-Next-to-Next-to Leading Order.”  Phys. Rev. Lett., 126(11):112002, 2021.
  * [BPX18] Radja Boughezal, Frank Petriello, and Hongxi Xing.  “Inclusive jet production as a probe of polarized parton distribution functions at a future EIC.”  Phys. Rev. D, 98(5):054031, 2018.
  * [Bra11] N. Brambilla et al.  “Heavy Quarkonium: Progress, Puzzles, and Opportunities.”  Eur. Phys. J. C, 71:1534, 2011.
  * [BS01] Christian W. Bauer and Iain W. Stewart.  “Invariant operators in collinear effective theory.”  Phys.Lett., B516:134–142, 2001.
  * [BS10] Thomas Becher and Matthew D. Schwartz.  “Direct photon production with effective field theory.”  JHEP, 02:040, 2010.
  * [BTR18] A. Bermudez, E. Tirrito, M. Rizzi, M. Lewenstein, and S. Hands.  “Gross–Neveu–Wilson model and correlated symmetry-protected topological phases.”  Annals Phys., 399:149–180, 2018.
  * [Bub05] Michael Buballa.  “NJL model analysis of quark matter at large density.”  Phys. Rept., 407:205–376, 2005.
  * [Bus96] D. Buskulic et al.  “Measurement of Lambda polarization from Z decays.”  Phys. Lett., B374:319–330, 1996.
  * [BVC22] Stefano Barison, Filippo Vicentini, Ignacio Cirac, and Giuseppe Carleo.  “Variational dynamics as a ground-state problem on a quantum computer.”  4 2022.
  * [CD07] Simone Conradi and Massimo D’Elia.  “Imaginary chemical potentials and the phase of the fermionic determinant.”  Phys. Rev. D, 76:074501, 2007.
  * [CDA19] Ian C. Cloët, Matthew R. Dietrich, John Arrington, Alexei Bazavov, Michael Bishof, Adam Freese, Alexey V. Gorshkov, Anna Grassellino, Kawtar Hafidi, Zubin Jacob, Michael McGuigan, Yannick Meurice, Zein-Eddine Meziani, Peter Mueller, Christine Muschik, James Osborn, Matthew Otten, Peter Petreczky, Tomas Polakovic, Alan Poon, Raphael Pooser, Alessandro Roggero, Mark Saffman, Brent VanDevender, Jiehang Zhang, and Erez Zohar.  “Opportunities for Nuclear Physics & Quantum Information Science.”, 2019.
  * [CdR07] P. Costa, C.A. de Sousa, M.C. Ruivo, and Yu.L. Kalinovsky.  “The QCD critical end point in the SU(3) Nambu–Jona-Lasinio model.”  Physics Letters B, 647(5):431–435, 2007.
  * [CGG18] J. Currie, T. Gehrmann, E. W. N. Glover, A. Huss, J. Niehues, and A. Vogt.  “N3LO corrections to jet production in deep inelastic scattering using the Projection-to-Born method.”  JHEP, 05:209, 2018.
  * [CGH17] James Currie, Thomas Gehrmann, Alexander Huss, and Jan Niehues.  “NNLO QCD corrections to jet production in deep inelastic scattering.”  JHEP, 07:018, 2017.  [Erratum: JHEP 12, 042 (2020)].
  * [CGH19] Chia Cheng Chang, Arjun Gambhir, Travis S. Humble, and Shigetoshi Sota.  “Quantum annealing for systems of polynomial equations.”  Scientific Reports, 9(1):10258, Jul 2019.
  * [CGK20] Justin Cammarota, Leonard Gamberg, Zhong-Bo Kang, Joshua A. Miller, Daniel Pitonyak, Alexei Prokudin, Ted C. Rogers, and Nobuo Sato.  “Origin of single transverse-spin asymmetries in high-energy collisions.”  Phys. Rev. D, 102(5):054002, 2020.
  * [Cha12] Serguei Chatrchyan et al.  “Measurement of jet fragmentation into charged particles in $pp$ and PbPb collisions at $\sqrt{s_{NN}}=2.76$ TeV.”  JHEP, 10:087, 2012.
  * [Cha14] Serguei Chatrchyan et al.  “Measurement of jet fragmentation in PbPb and pp collisions at $\sqrt{s_{NN}}=2.76$ TeV.”  Phys. Rev., C90(2):024908, 2014.
  * [CHI20] Bipasha Chakraborty, Masazumi Honda, Taku Izubuchi, Yuta Kikuchi, and Akio Tomiya.  “Digital Quantum Simulation of the Schwinger Model with Topological Term via Adiabatic State Preparation.”  arXiv:2001.00485 [cond-mat, physics:hep-lat, physics:hep-ph, physics:hep-th, physics:quant-ph], February 2020.  arXiv: 2001.00485.
  * [CHS15] Zhu-Fang Cui, Feng-Yao Hou, Yuan-Mei Shi, Yong-Long Wang, and Hong-Shi Zong.  “Progress in vacuum susceptibilities and their applications to the chiral phase transition of QCD.”  Annals of Physics, 358:172–205, 2015.  School of Physics at Nanjing University.
  * [CJN12] Jui-Yu Chiu, Ambar Jain, Duff Neill, and Ira Z. Rothstein.  “A Formalism for the Systematic Treatment of Rapidity Logarithms in Quantum Field Theory.”  JHEP, 05:084, 2012.
  * [CKM22] Alexander M. Czajka, Zhong-Bo Kang, Henry Ma, and Fanyi Zhao.  “Quantum simulation of chiral phase transitions.”  JHEP, 08:209, 2022.
  * [CKR16] Yang-Ting Chien, Zhong-Bo Kang, Felix Ringer, Ivan Vitev, and Hongxi Xing.  “Jet fragmentation functions in proton-proton collisions using soft-collinear effective theory.”  JHEP, 05:125, 2016.
  * [CKT20] Daniel Callos, Zhong-Bo Kang, and John Terry.  “Extracting the transverse momentum dependent polarizing fragmentation functions.”  Phys. Rev. D, 102(9):096007, 2020.
  * [CKT22] Alexander M. Czajka, Zhong-Bo Kang, Yuxuan Tee, and Fanyi Zhao.  “Studying chirality imbalance with quantum algorithms.”  10 2022.
  * [Cle20] Giuseppe Clemente et al.  “Quantum computation of thermal averages in the presence of a sign problem.”  Phys. Rev. D, 101(7):074510, 2020.
  * [Col93] John C. Collins.  “Fragmentation of transversely polarized quarks probed in transverse momentum distributions.”  Nucl.Phys., B396:161–182, 1993.
  * [Col02] John C. Collins.  “Leading twist single transverse-spin asymmetries: Drell-Yan and deep inelastic scattering.”  Phys.Lett., B536:43–48, 2002.
  * [Col13] John Collins.  Foundations of perturbative QCD, volume 32.  Cambridge University Press, 11 2013.
  * [Col18] J. I. et. al. Colless.  “Computation of Molecular Spectra on a Quantum Processor with an Error-Resilient Algorithm.”  Phys. Rev. X, 8:011021, Feb 2018.
  * [CRW19] Pedro Cal, Felix Ringer, and Wouter J. Waalewijn.  “The jet shape at $\mathrm{NLL}^{\prime}$.”  JHEP, 05:143, 2019.
  * [CS81] John C. Collins and Davison E. Soper.  “Back-To-Back Jets in QCD.”  Nucl.Phys., B193:381, 1981.
  * [CS16] Anirban Narayan Chowdhury and Rolando D. Somma.  “Quantum algorithms for Gibbs sampling and hitting-time estimation.”  2016\.
  * [CSS85] John C. Collins, Davison E. Soper, and George F. Sterman.  “Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production.”  Nucl. Phys. B, 250:199–224, 1985.
  * [CSS89] John C. Collins, Davison E. Soper, and George F. Sterman.  “Factorization of Hard Processes in QCD.”  Adv. Ser. Direct. High Energy Phys., 5:1–91, 1989.
  * [CSS08] Matteo Cacciari, Gavin P. Salam, and Gregory Soyez.  “The anti-$k_{t}$ jet clustering algorithm.”  JHEP, 04:063, 2008.
  * [CTG19] A. Chiesa, F. Tacchino, M. Grossi, P. Santini, I. Tavernelli, D. Gerace, and S. Carretta.  “Quantum hardware simulating four-dimensional inelastic neutron scattering.”  Nature Physics, 15(5):455–459, 2019.
  * [CV17a] Vincent Cheung and Ramona Vogt.  “Polarization of prompt $J/\psi$ and $\Upsilon$(1S) production in the color evaporation model.”  Phys. Rev. D, 96(5):054014, 2017.
  * [CV17b] Vincent Cheung and Ramona Vogt.  “Polarized Heavy Quarkonium Production in the Color Evaporation Model.”  Phys. Rev. D, 95(7):074021, 2017.
  * [CV18] Vincent Cheung and Ramona Vogt.  “Production and polarization of prompt $J/\psi$ in the improved color evaporation model using the $k_{T}$-factorization approach.”  Phys. Rev. D, 98(11):114029, 2018.
  * [CV19] Vincent Cheung and Ramona Vogt.  “Production and polarization of prompt $\varUpsilon$($n$S) in the improved color evaporation model using the $k_{T}$-factorization approach.”  Phys. Rev. D, 99(3):034007, 2019.
  * [CV21] Vincent Cheung and Ramona Vogt.  “Production and polarization of direct J/$\psi$ to O($\alpha$s3) in the improved color evaporation model in collinear factorization.”  Phys. Rev. D, 104(9):094026, 2021.
  * [CYW21] Chenfeng Cao, Yunlong Yu, Zipeng Wu, Nic Shannon, Bei Zeng, and Robert Joynt.  “Mitigating Algorithmic Errors in Quantum Optimization through Energy Extrapolation.”  9 2021.
  * [DCX13] Yi-lun Du, Zhu-fang Cui, Yong-hui Xia, and Hong-shi Zong.  “Discussions on the crossover property within the Nambu–Jona-Lasinio model.”  Phys. Rev. D, 88:114019, Dec 2013.
  * [DDL07] Massimo D’Elia, Francesco Di Renzo, and Maria Paola Lombardo.  “The Strongly interacting quark gluon plasma, and the critical behaviour of QCD at imaginary mu.”  Phys. Rev. D, 76:114509, 2007.
  * [Dev21] Cirq Developers.  “Cirq.”, August 2021.  See full list of authors on Github: https://github .com/quantumlib/Cirq/graphs/contributors.
  * [DGK11] Umberto D’Alesio, Leonard Gamberg, Zhong-Bo Kang, Francesco Murgia, and Cristian Pisano.  “Testing the process dependence of the Sivers function via hadron distributions inside a jet.”  Phys. Lett. B, 704:637–640, 2011.
  * [Din19] H. T. Ding et al.  “Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD.”  Phys. Rev. Lett., 123(6):062002, 2019.
  * [DKL16] Lin Dai, Chul Kim, and Adam K. Leibovich.  “Fragmentation of a Jet with Small Radius.”  Phys. Rev., D94(11):114023, 2016.
  * [DL03] Massimo D’Elia and Maria-Paola Lombardo.  “Finite density QCD via imaginary chemical potential.”  Phys. Rev. D, 67:014505, 2003.
  * [DLM21] H. T. Ding, S. T. Li, Swagato Mukherjee, A. Tomiya, X. D. Wang, and Y. Zhang.  “Correlated Dirac Eigenvalues and Axial Anomaly in Chiral Symmetric QCD.”  Phys. Rev. Lett., 126(8):082001, 2021.
  * [DLX15] Yi-Lun Du, Ya Lu, Shu-Sheng Xu, Zhu-Fang Cui, Chao Shi, and Hong-Shi Zong.  “Susceptibilities and critical exponents within the Nambu–Jona-Lasinio model.”  International Journal of Modern Physics A, 30(34):1550199, 2015.
  * [DM16] M. Dalmonte and S. Montangero.  “Lattice gauge theories simulations in the quantum information era.”  Contemporary Physics, 57(3):388–412, July 2016.  arXiv: 1602.03776.
  * [DMG21] Olivia Di Matteo, Anna McCoy, Peter Gysbers, Takayuki Miyagi, R. M. Woloshyn, and Petr Navrátil.  “Improving Hamiltonian encodings with the Gray code.”  Physical Review A, 103(4), Apr 2021.
  * [DMH18] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski.  “Cloud Quantum Computing of an Atomic Nucleus.”  Phys. Rev. Lett., 120:210501, May 2018.
  * [DMM21a] Wibe A. De Jong, Mekena Metcalf, James Mulligan, Mateusz Płoskoń, Felix Ringer, and Xiaojun Yao.  “Quantum simulation of open quantum systems in heavy-ion collisions.”  Phys. Rev. D, 104(5):051501, 2021.
  * [DMM21b] Wibe A. De Jong, Mekena Metcalf, James Mulligan, Mateusz Płoskoń, Felix Ringer, and Xiaojun Yao.  “Quantum simulation of open quantum systems in heavy-ion collisions.”  Phys. Rev. D, 104(5):051501, 2021.
  * [DMP11] Umberto D’Alesio, Francesco Murgia, and Cristian Pisano.  “Azimuthal asymmetries for hadron distributions inside a jet in hadronic collisions.”  Phys. Rev. D, 83:034021, 2011.
  * [DMP17] Umberto D’Alesio, Francesco Murgia, and Cristian Pisano.  “Testing the universality of the Collins function in pion-jet production at RHIC.”  Phys. Lett. B, 773:300–306, 2017.
  * [DMZ20] Umberto D’Alesio, Francesco Murgia, and Marco Zaccheddu.  “First extraction of the $\Lambda$ polarising fragmentation function from Belle $e^{+}e^{-}$ data.”  3 2020.
  * [Dok77] Yuri L. Dokshitzer.  “Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics.”  Sov. Phys. JETP, 46:641–653, 1977.
  * [DRS21] Zohreh Davoudi, Indrakshi Raychowdhury, and Andrew Shaw.  “Search for efficient formulations for Hamiltonian simulation of non-Abelian lattice gauge theories.”  Phys. Rev. D, 104(7):074505, 2021.
  * [DS05] M. Diehl and S. Sapeta.  “On the analysis of lepton scattering on longitudinally or transversely polarized protons.”  Eur. Phys. J. C, 41:515–533, 2005.
  * [DS09] Massimo D’Elia and Francesco Sanfilippo.  “Thermodynamics of two flavor QCD from imaginary chemical potentials.”  Phys. Rev. D, 80:014502, 2009.
  * [EER21] M. G. Echevarria, I. L. Egusquiza, E. Rico, and G. Schnell.  “Quantum simulation of light-front parton correlators.”  Phys. Rev. D, 104(1):014512, 2021.
  * [EGS22] Markus A. Ebert, Anjie Gao, and Iain W. Stewart.  “Factorization for azimuthal asymmetries in SIDIS at next-to-leading power.”  JHEP, 06:007, 2022.
  * [EIK14] Miguel G. Echevarria, Ahmad Idilbi, Zhong-Bo Kang, and Ivan Vitev.  “QCD Evolution of the Sivers Asymmetry.”  Phys. Rev. D, 89:074013, 2014.
  * [EIS12] Miguel G. Echevarria, Ahmad Idilbi, and Ignazio Scimemi.  “Factorization Theorem For Drell-Yan At Low $q_{T}$ And Transverse Momentum Distributions On-The-Light-Cone.”  JHEP, 07:002, 2012.
  * [EIS13a] Miguel G. Echevarria, Ahmad Idilbi, Andreas SchÃ€fer, and Ignazio Scimemi.  “Model-Independent Evolution of Transverse Momentum Dependent Distribution Functions (TMDs) at NNLL.”  Eur.Phys.J., C73(12):2636, 2013.
  * [EIS13b] Miguel G. Echevarría, Ahmad Idilbi, and Ignazio Scimemi.  “Soft and Collinear Factorization and Transverse Momentum Dependent Parton Distribution Functions.”  Phys. Lett. B, 726:795–801, 2013.
  * [EIS14] Miguel G. Echevarria, Ahmad Idilbi, and Ignazio Scimemi.  “Unified treatment of the QCD evolution of all (un-)polarized transverse momentum dependent functions: Collins function as a study case.”  Phys. Rev. D, 90(1):014003, 2014.