Quantum Simulation of Chiral Phase Transitions

The Nambu-Jona-Lasinio (NJL) model has been widely studied for investigating the chiral phase structure of strongly interacting matter. The study of the thermodynamics of field theories within the framework of Lattice Field Theory is limited by the sign problem, which prevents Monte Carlo evaluation of the functional integral at a finite chemical potential. Using the quantum imaginary time evolution (QITE) algorithm, we construct a quantum simulation for the $(1+1)$ dimensional NJL model at finite temperature and finite chemical potential. We observe consistency among digital quantum simulation, exact diagonalization, and analytical solution, indicating further applications of quantum computing in simulating QCD thermodynamics.

As described in Quantum ChromoDynamics (QCD), the strongly interacting matter is believed to have a complex phase structure, which is of great importance in theoretical physics [76][77][78][79][80][81][82]. In the infinite quark limit, one has the deconfinement phase transition [83] as a function of temperature as a result of the spontaneous Z N symmetry breaking [84], while in the massless quark limit, one finds an order parameter for the chiral phase transition [85], which has been extensively studied in [86][87][88][89][90] 1 . Because QCD is dominated by non-perturbative effects at low energies, it is difficult to study the chiral phase transition of strongly interacting matter.
To investigate the structure of the QCD chiral phase diagram 2 , calculations are typically carried out in lattice QCD. However, the utility of this traditional approach is limited by the fermion sign problem. At nonzero baryochemical potential µ, the (Euclidean) QCD action S is no longer constrained to be real, so that the quantity e S can no longer be used as a weight for a Monte Carlo evaluation of thermal expectation values [83,91]. The traditional approach to this problem involves simulations on imaginary chemical potentials [92][93][94][95][96][97][98][99] and expansion of the expectation value O of a given observable O in powers of µ/T [100][101][102][103][104], namely O = k O k (T )(µ/T ) k , and evaluating the coefficients O k (T ) on the lattice [83,91,101]. While as pointed out in [83], noisy results emerge in the higher order coefficients of the Taylor series from the direct simulation. In addition to the fact that the domain of validity of such an approach is limited to a (potentially small) range of baryochemical potentials, it assumes that for fixed T the expectation values O are analytic in µ, an assumption that should fail for at least some observable in the vicinity of a finite-order phase transition. The sign problem is thus a significant obstacle to understanding the phase structure of nuclear matter. However, this is not a defect of the underlying theory (QCD), but of the attempt to simulate quantum statistics (via the functional integral) with a classical algorithm (Markov Chain Monte Carlo). By modeling the lattice system on a quantum computer, one can leverage the statistical properties of the computer itself to eliminate the need for a Monte Carlo analysis. The fact that quantum simulations can avoid the sign problem has also been pointed out in [105] 3 .
To study the QCD chiral phase transition, the Nambu-Jona-Lasinio (NJL) model [106,107] has been used as a convenient and practical tool due to the easy access to the dynamical mass and the mechanism of chiral symmetry breaking [76,77,[108][109][110][111][112][113][114][115][116]. The NJL model, an effective model for low-energy two-flavor QCD, is also amenable to analytical calculations at finite temperature and chemical potential. A simplified version of the NJL model, the Gross-Neveu (GN) model [117], is a renormalizable and asymptotically free (1 + 1)dimensional theory comprised of N fermion species interacting via four-fermion contact interaction, and it is known to exhibit chiral symmetry breaking at low energy through the generation of a nonvanishing chiral condensate ψ ψ . We will analyze the chiral condensate of the GN model at finite temperature and chemical potential [118] and compare the results to quantum simulations.
In this work, we use a 4-qubit system to simulate the (1 + 1) dimensional NJL Hamiltonian using single-qubit gates and the CNOT gate, via the Jordan-Wigner transformation [119]. In Sec. 2, we provide analytical calculations for the finite temperature behavior of the NJL model including chemical potentials. A description of the discretization of the field theory on a lattice and the quantum algorithm we use, the QITE algorithm, is given in Sec. 3. Then we present a comparison between the results obtained by analytical calculation, exact diagonalization, and quantum simulation at various chemical potentials in Sec. 4, where a strong consistency is found. We summarize our results in Sec. 5, and discuss 2 That is: the division of the temperature-chemical potential plane into regions based on the qualitative physical properties of strongly interacting matter, "phases", with these phases coexisting at boundaries. 3 Evidently, the idea of using quantum systems to simulate other quantum systems was first suggested by Feynman [6].
further directions of using quantum computing for studying finite-temperature properties of QCD.

Nambu-Jona-Lasinio model in 1 + 1 dimensions
In this section, we provide an overview for the Nambu-Jona-Lasinio model and its phase transition at finite temperature and finite chemical potential. The Lagrangian density of the Nambu-Jona-Lasinio (NJL) model in (1 + 1)-dimensional Minkowski space is [106,107] where / ∂ = γ µ ∂ µ , τ a are the Pauli matrices in the isospin space, m is the bare quark mass and g is the dimensionless coupling constant. In (1 + 1) dimensions, the gamma matrices γ 0 , γ 1 and γ 5 are given by, where σ X , σ Y and σ Z are Pauli matrices. A simplified variant of the NJL model is the Gross-Neveu (GN) model [117] with the following Lagrangian In the chiral limit m = 0, this Lagrangian is symmetric under left-action by the discrete symmetry group Z 2,L × Z 2,R ≡ {±1, ±γ 5 } on the fields, i.e. ψ → G · ψ for any group element G ∈ Z 2,L × Z 2,R .
In this work, we study the chiral phase transition of GN type model with non-zero chemical potential µ [118], which has been applied for studying the chiral phase transition at finite chemical potential µ and temperature T . Throughout this paper, we refer to the Lagrangian in Eq. (2.4) as we mention "NJL model".
To distinguish between different phases, we study the chiral condensate ψ ψ in the vacuum, which is known to be an order parameter [106,[120][121][122][123] in the chiral limit with zero quark mass. Since the quantityψψ transforms nontrivially under Z 2,L ×Z 2,R , a nonzero value of ψ ψ requires spontaneous breaking of this symmetry. Thus, the transition from nonvanishing to vanishing values of the condensate signals a chiral phase transition. Thus, with non-zero quark mass, the chiral condensate ψ ψ can be regarded as a quasi-order parameter [124].
The transition between the chirally symmetric and broken phases can be analyzed in the mean field approximation. In this case, the quantityψψ can be written asψψ = ψ ψ + σ [117,125], where σ is a real scalar field ("fluctuations"), assumed to be small, i.e. |σ/ ψ ψ | 1, and ψ ψ is just a constant. Independent of the coordinates, that is. Such a coordinate-independent chiral condensate ψ ψ is sometimes referred to as global chiral condensate [126,127], to be distinguished from the local chiral condensate ψ (x)ψ(x) which does depend on the coordinate. The global chiral condensate ψ ψ will depend on the chemical potential µ and the temperature T . Then, the Lagrangian becomes Neglecting O(σ 2 ) terms, this is the Lagrangian L Dirac of a free Dirac fermion (at finite chemical potential) with mass M = m − 2g ψ ψ with a constant potential V = g ψ ψ 2 = (M − m) 2 /4g. Thus, we have an effective, linearized, Lagrangian is the effective mass of the fermion field. The value of the chiral condensate ψ ψ , or equivalently the effective mass M , can then be determined from the condition of thermal equilibrium, i.e. that the Grand Canonical Potential Ω(µ, T ; M ) = − T L ln Z is minimized 4 . The partition function Z is given, in the path integral formulation, by [128] where the Euclidean Lagrangian L E is obtained from making the change of variables to imaginary time τ = it, and β = 1/T is the coldness. In particular, from our effective From this, we see that the Grand Canonical Potential Ω is given by Ω = Ω Dirac + V, where Ω Dirac is the Grand Canonical Potential of a free Dirac field of mass M (in one spatial dimension), which is given by [128,129] (2.10) 4 Here, the volume is the "volume" of 1 dimensional space, so it is just a length L.
This quantity diverges, so it must be regularized. However, for the sake of comparison with the Lattice model given in Sec. 3 below, we have a natural momentum cutoff Λ = π/a, where a is the lattice spacing. By imposing this hard momentum cutoff in the integral of Eq. (2.10), we may numerically compute the effective mass M at fixed values of µ, T by minimizing Ω with respect to M : ∂Ω/∂M = 0, which leads to the so-called gap equation. The chiral condensate can then be calculated from ψ ψ = (m − M )/2g. The results obtained this way will be referred to as the analytical calculation in the following sections. We performed this analysis for µ, T in the range [0 MeV, 300 MeV] using the model parameters m = 100 MeV, a = 1 MeV −1 and g = 1 that are relevant for our simulation. The resulting mass surface is plotted in Fig. 1 below. As can be seen, for sufficiently low µ 2 + T 2 there Historically the GN model was studied specifically because it is asymptotically free [117]. Thus, at asymptotically high temperatures/chemical potentials we expect to recover a free field theory.

Quantum algorithms for the NJL model
In this section, we provide the ingredients of the quantum algorithms for the NJL model, in particular, we first show the discretization of the NJL Hamiltonian on a lattice grid, then we introduce the quantum imaginary time evolution algorithm.

NJL Hamiltonian at the lattice fermion model
Starting from the Lagrangian density given in Eq. (2.4), one can obtain the NJL Hamiltonian density as follows Given a Dirac fermion field ψ(x) with components ρ(x) and η(x), to discretize the theory we place the fermion field on a spatial lattice, with spacing a, and set a staggered fermion field χ n on each lattice site [130]. Note that staggered fermions are used to avoid a subtlety known as the fermion doubling problem, where naïve discretization of the fermion field leads to a theory with too many degrees of freedom [131].
With even N and integers n = 0, 1, · · · N/2 − 1, sites 2n are assigned to the stagger fermion field χ 2n / √ a with spinor ρ(x = n). And similarly, sites 2n + 1 are assigned to the stagger fermion field χ 2n+1 / √ a with spinor η(x = n). Then the position x is discretized into N/2 position points x 0 , · · · , x N/2−1 with spacing x n − x n−1 = δx = a. And the Dirac fermion field ψ(x) at position x = n on the lattice is represented by staggered fermion fields [86,[132][133][134][135], Also note that the partial derivative of the field ∂ 1 ψ(x) is defined by the forward-backward derivative convention [136] in lattice QCD and given by To show the discretization of the fermion field, we specify ψ n to represent ψ(x = n), then one can easily find the following equivalence Keep in mind that N is an even number and thus N/2 is always an integer in the above formula. Note that we have taken a periodic boundary condition, namely χ N → χ 0 , and the χ N −1 field would be coupled with the χ 0 field. For illustration, we have explicitly written this boundary term out in Eq. (3.6). Therefore, the Hamiltonian in spin representation is given by Once one converts the Dirac fermions to staggered fermions, by applying the Jordan-Wigner transformation [119], one obtains the operators in spin representations for quantum simulation. In Eq. (3.9), σ X,n represents a Pauli-X matrix acting on the n-th grid on the lattice, etc. Then one is able to verify the following equivalence up to a constant term, where the subscripts indicate the index of qubit where the single-qubit gate is acting on.
Here again we have explicitly written out the boundary term in Eqs. (3.10) and (3.13). Note that even though Eq. (3.13) has a slightly different form, it is equivalent to the expression given in [130]. In other words, by applying the Jordan-Wigner transformation from Eq. (3.9) on Eq. (3.8), we obtain the Hamiltonian H = 5 j=1 H j decomposed into the following 5 pieces for constructing the quantum algorithm, 14)

16)
With these Hamiltonians, one is able to apply Suzuki-Trotter decomposition [137,138] and evolve the Hamiltonian using a quantum simulation.

Quantum imaginary time evolution algorithm
To perform a quantum simulation for calculating the thermal properties of the NJL Hamiltonian in Eq. (3.1) at finite temperatures, we use the quantum imaginary time evolution (QITE) algorithm developed in [67]. Quantum simulation of a Hamiltonian H traditionally uses the Suzuki-Trotter decomposition to approximate the evolution e −iHt |Ψ by a series of local unitary operations which are implementable on a quantum computer. In QITE, the imaginary time substitution β = it is first made, followed by a series of decomposition and approximation steps. QITE is especially useful as a subroutine within the quantum minimally entangled typical thermal states (QMETTS) algorithm [67], which is used to calculate thermal averages of observables. The benefit of using imaginary time evolution to calculate thermal averages is that, compared to other quantum algorithms [139? , 140] for computing thermal averages, QITE does not require deep circuits or any ancilla qubits. We will use QITE and QMETTS to calculate the chiral condensate ψ ψ for the purpose of studying chiral phase transition.
The QITE algorithm aims to decompose the imaginary time evolution operator e −βH . Taking the Hamiltonian as given in Sec. 3.1, Trotterization of the imaginary time evolution operator yields where ∆β is the chosen imaginary time step size and n = β/∆β is the number of iterations needed to reach imaginary time β.
In contrast with real-time quantum simulation, each intermediate evolution e −∆βH is non-unitary, and cannot be directly implemented in terms of quantum gates. The crux of QITE lies in approximating the non-unitary evolution with a unitary operation so that (with proper normalization c(∆β) = Ψ|e −2∆βH |Ψ ) 1 c(∆β) e −∆βH |Ψ ≈ e −i∆βA |Ψ , (3.20) where A is a Hermitian operator parameterized with a linear combination of Pauli operators, where a µ are coefficients to be determined below andσ µ is short-hand notation for the Pauli string given above. The subscript µ in a µ andσ µ is a composite index µ = (µ 1 , µ 2 , · · · , µ nq ) with µ l ∈ {I, X, Y, Z} and thus runs over all possible 4 nq subsets of Pauli strings. On the other hand, σ (l) µ l represents a Pauli-µ l operator acting on the l-th qubit in the circuit. In this work, we choose the number of qubits n q = 4. When ∆β is very small, one can expand both sides of Eq. (3.20)  To evaluate the Hermitian operator A, we need to minimize the objective function F (a) to obtain the values of a µ . Here F (a) is given by where F (0) = || |∆Ψ H (β) || 2 corresponds to F (a) at a = 0, thus this term is independent of the values of a µ . To minimize the objective function F (a), one is required to take a derivative of a and solve the linear equation below Therefore, by defining the elements of matrix S and vector b as one obtains the parameters a µ by solving the linear equation With a µ obtained, the unitary evolution e −i∆βA is implementable in terms of quantum gates. We illustrate this implementation in Fig. 2, where U j |Ψ(β j ) = e −i∆βA |Ψ(β j ) = |Ψ(β j + ∆β) . Namely, each block evolves state |Ψ(β j ) to |Ψ(β j + ∆β) , and after n steps, one obtains the state at the target total evolution time β. Using the thermal state |Ψ(β/2) generated by the QITE algorithm, we then find the expectation value of an observableÔ at finite temperature T = 1/β, which can be written as Tr(e −βĤ ) . where S is the complete set of 2 nq -dimensional Hilbert space. In other words, S = {|0000 , |0001 , |0010 , · · · , |1111 } with 16 orthogonal basis for n q = 4 qubits used in our quantum simulation below. Below, we will compute the expectation value of the global chiral condensate ψ ψ . In doing so, we choose the operatorÔ in Eq. (3.31) witĥ O = dxψ(x)ψ(x)/ dx and applying Eq. (3.11). In other words, we take the average of ψ(x)ψ(x) on all the lattice sites. In Fig. 3, we present the imaginary time evolution of the chiral condensate ψ ψ at chemical potentials µ = 0 MeV (left), 100 MeV (middle) and µ = 150 MeV (right) with small-time step ∆β = 0.001 MeV −1 (solid line) and ∆β = 0.005 MeV −1 (dashed line). Our quantum algorithm contains 4 qubits. The quantum simulation is performed on a (1 + 1)-dimensional NJL model with the dimensionless coupling constant g = 1 and the quark mass m = 100 MeV on a lattice grid with spacing a = 1 MeV −1 . The exact diagonalization carried out with discretization of the NJL Hamiltonian as given in Eq. (3.8) is shown in diamond points for reference. As expected, one obtains more accurate quantum simulations with smaller evolution steps. In the following sections, we will apply small-time step ∆β = 0.001 MeV −1 for the QITE algorithm and compare the effective mass M among quantum simulation, exact diagonalization, and theoretical analysis.

Results
In this section, we demonstrate the (1 + 1)-dimensional NJL chiral phase transition by plotting the quark condensate ψ ψ as a function of temperature T and chemical potential µ. By working on a (1 + 1)-dimensional model, from the theoretical point of view, we are able to derive an analytical result by solving the gap equation as illustrated in Sec. 2 as a comparison to the quantum simulations. At the same time, restriction on (1 + 1) dimension enables us to design a quantum circuit with a small amount of qubits, which leads to feasibility of implementing such circuit on currently available hardwares for further research.
The parameters that will be used in our work are set as: lattice spacing a = 1 MeV −1 , bare mass m = 100 MeV and the dimensionless coupling constant g = 1. Our results are obtained from the following three directions: 1. Simulation of the thermal states with the QITE algorithm; 2. Exact diagonalization of the Hamiltonian in spin representation; 3. Analytical calculation by numerically solving the gap equation.
To obtain the quantum simulation results, one can install a quantum simulation package such as PYQUILL [141], TEQUILA [142], QFORTE [143], QISKIT [144], XACC [145], CIRQ [146], and so on to apply the quantum circuits. These quantum simulation packages give similar results for quantum simulations. They are python software libraries and are used for providing expected outputs of an ideal quantum computer. Algorithm availability for some packages is discussed in [147] and a more complete list of general quantum simulation packages can be found in [148]. In this work, we work within the framework of the open-source software package QFORTE [143], which provides an efficient state simulator and quantum algorithms library. We modify the implementation of QITE within QFORTE to calculate properties of the NJL Hamiltonian in quantum simulation.
In Fig. 4, we plot the temperature dependence of the chiral condensate ψ ψ at chemical potentials set at µ = 0, 25, · · · , 200 MeV. The quantum simulation results are presented in comparison with theoretical calculations and exact diagonalization. The diamond data points are from QITE simulation, while the solid curves are from the theory calculation as explained in Sec. 2. On the other hand, the dashed curves are from the exact diagonalization which is performed by applying the matrix form of the discretized NJL Hamiltonian given in Eqs. (3.14) -(3.18). We find that our quantum simulations are in good agreement with both exact diagonalization and theoretical analyses. For chemical potentials below 100 MeV, the quark condensate increases as the temperature rises. While for large chemical potentials, the quark condensate shows opposite behavior at a small temperature range then slowly increases with temperature as indicated by the figure. And the quark condensate at various finite chemical potentials µ tends to converge at large temperatures.
In Fig. 5, we plot the quark condensate ψ ψ as a function of temperature for different ratios of the chemical potential and temperature, µ/T , from 0 to 8 with the same set of parameters chosen for Fig. 4. Here the quantum simulations using the QITE algorithm are plotted in diamond points. For comparison, we also provide the analytical calculation and exact diagonalization results as shown in solid lines and dashed lines respectively. At low temperatures, the quark condensate is about −2 MeV. As the temperature increases, the quark condensate increases to around 0 MeV. As indicated by the curves, at a large ratio µ/T , the quark condensate increases from −2 to 0 MeV more rapidly.
Similarly, we plot the quark condensate ψ ψ for finite chemical potentials at various temperatures in Fig. 6, for T = 25, 50, 100, 125, 250 and 500 MeV in the ψ ψ − µ plane.  At low temperatures and chemical potentials, the quark condensate has a larger absolute value. As the chemical potential µ increases, the quark condensate rapidly increases to around 0 MeV. While for high temperatures, the quark condensate barely changes at different chemical potentials. Therefore, one observes a more obvious phase transition at larger temperatures in the ψ ψ − µ plot.

Conclusion
In this work, we have constructed a quantum simulation for the chiral phase transition of the 1+1 dimensional NJL model at finite temperature and chemical potentials with the QITE algorithm, where we applied a 4-qubit quantum circuit, and simulated the NJL Hamiltonian with quantum gates. We observe a consistency among digital quantum simulation, exact diagonalization and analytical solution, which indicates rich applications of quantum computing in simulating finite-temperature behaviors for QCD in the future.
The demonstrated insensitivity of the quantum simulation's efficacy to the value of the chemical potential (and likely other external parameters) presents an exciting avenue for exploring finite density effects in QCD and other field theories. This aspect of quark physics remains largely unexplored (in comparison to other areas) due to technical constraints (nonperturbative dynamics, the failure of Monte Carlo due to the sign problem, etc.) preventing the use of traditional computational methods (perturbation theory, Lattice Gauge Theory, Semi-Classical Methods, etc.). With scalable quantum computer technology on the horizon, performing Lattice QCD calculations on quantum computers is not only increasingly possible, but practical as well.
This work, as well as previous studies, demonstrates that NISQ quantum computers are capable of producing consistent and correct answers to physical problems, some of which cannot be efficiently or effectively solved with classical computing algorithms. This signals bright prospects for future applications to nonperturbative QCD and beyond.