Quantum variational approach to lattice gauge theory at nonzero density

The simulation of dense fermionic matters is a long-standing problem in lattice gauge theory. One hopeful solution would be the use of quantum computers. In this paper, digital quantum simulation is designed for lattice gauge theory at nonzero density. The quantum variational algorithm is adopted to obtain the ground state at nonzero density. A benchmark test is performed in the lattice Schwinger model.


I. INTRODUCTION
Lattice quantum chromodynamics (QCD) at nonzero baryon density suffers from the sign problem. Many challenging approaches were proposed to overcome the sign problem, but the general solution has not yet been discovered. The sign problem might be beyond the limit of classical supercomputers. In recent years, quantum technology is very rapidly developing. Quantum computers open up new possibilities in lattice QCD [1][2][3]. Although the lattice QCD simulation is difficult for the present noisy-intermediate-scale-quantum (NISQ) devices, it might be feasible in future. Toward the final goal, we should start discussing how to formulate the simulation of dense QCD in quantum computers.
In classical computers, lattice gauge theory is usually formulated by the path integral with the Euclidean action. This is nothing but thermodynamics with the grand canonical partition function. Physical observables are calculated as functions of temperature and chemical potential. On the other hand, in gate-based quantum computers, qubits are state vectors and quantum gates are unitary operations, so the operator formalism with the Hamiltonian is natural. In the Hamiltonian formalism, it is easier to consider energy eigenstates than thermal average. When we are interested in the physics at zero temperature, we need only the ground state. As for density, it is easier to fix particle number density than chemical potential. Since net electric charge is conserved in gauge theory, we can take the basis where the total number of charged particles is fixed. Therefore, we can study the physics at zero temperature and nonzero density by two steps: we prepare the ground state with the total particle number fixed, and then calculate physical observables for the prepared state. This is the standard procedure in condensed matter physics and nuclear physics, where the Hamiltonian formalism is familiar. Although the Hamiltonian formalism has been considered in the long history of lattice gauge theory [4], we should reconsider it from more concrete point of view.
In this paper, we design the computational strategy and algorithm for the quantum simulation of lattice gauge theory at nonzero density. We give some demonstrations in a simple theory, the lattice Schwinger model. For the demonstrations, we used the simulator in IBM Quantum, which is a classical algorithm to simulate quantum computation, not a real quantum computer. The results are free from any quantum device noise. When system size is small, the Hamiltonian can be diagonalized and all the eigenstates can be obtained by a classical computer. We also performed the exact diagonalization for comparison.

II. SCHWINGER MODEL
The lattice Schwinger model is the lattice discretization of quantum electrodynamics in 1 + 1 dimensions. The model is frequently used in the researches of quantum computation, both in theory [5][6][7][8][9][10] and experiment [11][12][13][14][15]. The theoretical formulation and the device implementation are very well explained in the literature. We here summarize necessary equations.
The fundamental building blocks are the spinless fermion χ n , the electric field L n , and the spatial gauge field U n . The Hamiltonian on the N -site lattice is where J and w are c-number parameters. Open boundary condition is assumed. In the staggered fermion formulation, the fermion χ n is assigned as a particle on even lattice sites and as an antiparticle on odd lattice sites. The fermion number density operator is given by The electric fields and the number density are constrained by the Gauss law The total particle number of the fermions is given by We fixed one boundary condition, L 0 = 0, for simplicity. The other boundary condition is automatic, L N = Q, when the total particle number is given. For quantum computation, L n and U n are eliminated by the Gauss law (5) and the field redefinition of χ n , respectively, and χ n is transformed to the Pauli matrices by the Jordan-Wigner transformation [16]. Then we arXiv:2104.10669v2 [hep-lat] 12 Jul 2021 Each term of the Hamiltonian is realized by one Pauli gate The "nonzero density" means that the total particle number is nonzero. From the commutation relation [Q, H] = 0, the total particle number is conserved. In other words, the Hamiltonian is block-diagonal and the Hilbert space is decomposed into subspaces q = N/2, · · · , 1, 0, −1, · · · , −N/2. The states in different particle number sectors cannot mix with each other. The expectation value O of a physical observable can be written as a function of q if O does not create or annihilate particles. Changing the particle number sector by hand, we can study the density-dependence of the physical observable.

III. QUANTUM ADIABATIC ALGORITHM
The quantum adiabatic algorithm is the computational scheme to obtain the ground state of arbitrary Hamiltonian [17,18]. It is based on the adiabatic theorem; if the Hamiltonian is gradually changed from one to another and if the initial state is the ground state, the system keeps staying the ground state. The quantum adiabatic algorithm was applied to many models, including the lattice Schwinger model with q = 0 [8].
In the lattice Schwinger model, the ground state of the electric field Hamiltonian H L is the pure state depicted in Fig. 1. The initial state Ψ 0 (q) is set to the ground state of H L . The remaining part, the fermion Hamiltonian H χ , is treated as the perturbation term. The algorithm is written by In the limit of T → ∞, the evolution is adiabatic and the ground state of H is obtained. The evolution operator is approximated by the Suzuki-Trotter decomposition, We can show the commutation relation [Q, U (s)] = 0. The commutation relation guarantees that the obtained state Ψ(q) has the same total particle number as the initial state Ψ 0 (q).
There is a technical but important issue. In general, the decomposition is not unique. For example, another choice is U (s) = n e −i δt XnXn+1 · · · . This choice is however bad because it does not satisfy the commutation relation. We must choose the decomposition which conserves the total particle number, otherwise the obtained state Ψ(q) is not the eigenstate of Q.
We demonstrate the quantum adiabatic algorithm in the lattice Schwinger model at q = 0. We fix the parameters J = w = 1 and use the dimensionless unit scaled by them. The lattice size is N = 8, so the dimension of the Hilbert space is 2 N = 256. The ground-state energy is shown in Fig. 2. The energy decreases as T increases, and eventually in T > 10, becomes insensitive to T (and δt). This is exactly the expected behavior in this algorithm. The energy in T > 10 is indeed consistent with the exact ground-state energy obtained by matrix diagonalization. The calculation seems successful, but it has a potential problem. The operation U (s) needs to be manipulated many times; e.g., S = 50 for δt = 0.2 and S = 100 for δt = 0.1 in Fig. 2. In real NISQ devices, that enhances decoherence and makes quantum error out of control. A noise-robust algorithm is required for the real quantum simulation.

IV. QUANTUM VARIATIONAL ALGORITHM
We adopt the quantum variational algorithm to suppress the noise. The prototype of the quantum variational algorithm was proposed in quantum chemistry [19], and then applied to other fields, in particular, to the lattice Schwinger model [6,14]. The algorithm consists of the following quantum and classical processes: (a) One generates a certain quantum state and calculate its energy in a quantum computer. (b) Based on the calculated energy, one classically determines the feedback to improve the quantum state. The processes (a) and (b) are repeated until the energy converges to the minimum.
The quality of the variational calculation strongly depends on the variational ansatz, i.e., how to parameterize quantum states. We take the variational ansatz [20,21]  There are 2S variational parameters {γ s , β s }(s = 1, · · · , S). Because V (s) is the same form as U (s) in Eq. (14), it satisfies the commutation relation [Q, V (s)] = 0 and conserves the total particle number. It is trivial that the results are improved by increasing the integer S because the number of variational parameters increases. In this work, we fix the number of variational parameters by further simplification. Let us consider two parameterizations: the one-parameter (α) ansatz γ s = α, β s = α s S (18) and the two-parameter (α, β) ansatz They are inspired by the quantum adiabatic algorithm. In the limit of S → ∞, they can reproduce Eqs. (11) to (14), so the convergence to the exact ground state is a priori ensured.
We tested the variational calculation with the same simulation condition as in the adiabatic calculation in Sec. III. The dependence of the ground-state energy on the integer S is shown in Fig. 3. While the one-parameter ansatz gradually converges as a function of S, the two-parameter ansatz rapidly converges and reproduces the exact value even at S = 2. The two-parameter ansatz also works well in other particle number sectors, as shown in Fig. 4. The value of S is much smaller than that of the adiabatic calculation in Sec. III. The calculation will be less influenced by the device noise in real quantum computers.
In the analysis of nonzero density, we often want to know the density-dependence of physical observables. This is now possible because we obtained the state vector of the ground state in each particle number sector. As an example, we cal- culated the volume-averaged chiral condensate The second term is the offset to make the chiral condensate in the perturbative vacuum, Ψ 0 (q = 0) in Fig. 1, zero. The chiral condensate is plotted as a function of the particle number in Fig. 5. We see that the chiral condensate matches the exact value with good accuracy. Note that the purpose of this calculation is just to show that we can calculate any observable other than energy. Physically speaking, the chiral condensate in the massless lattice Schwinger model is lattice artifact. (The particle-antiparticle symmetry is realized by the even-odd site symmetry in the staggered fermion formulation, but it is explicitly broken by the boundary condition. This explicit symmetry breaking induces the nonzero chiral condensate.) We can do physical analysis if we use physical theory. The results almost perfectly agree with the exact values in this small system. In realistic theory with large Hilbert space, such perfect agreement is not expected. We can only obtain an approximate ground state. To improve the quality of the approximation, we need to increase S or the number of variational parameters. The decoherence will become nonnegligible if S is large and the computational cost to search the minimum will blow up if the number of variational parameters is large; this is the trade-off between the quantum noise and the classical computational cost.

V. TOWARD DENSE QCD
In this work, we proposed the possibility of quantum computers as a tool to study lattice gauge theory at nonzero den- sity. The lattice Schwinger model was used because it is the best we can simulate with the current resource of quantum computers. The basic strategy is quite general and applicable to other lattice gauge theories, e.g., in two-or threedimensions or with non-Abelian gauge group. We might be able to challenge dense QCD in a far future.
This approach has two differences from the conventional lattice gauge theory. One is the computational device: quantum computers or classical computers. Since this approach does not rely on the Monte Carlo method, it does not encounter the sign problem. (The data obtained by quantum computers have statistical error, but its origin is essentially different.) In general, deterministic calculation has the problem that computational time exponentially increases as system size increases, instead of the sign problem. The variational calculation must have the same problem if it is executed in classical computers. Even though quantum computation speeds up one part of the calculation, it is non-trivial whether the total computational cost is reasonable or not. We need further investigation to estimate it. The other is the theoretical formalism: the Hamiltonian with fixed particle number or the Euclidean path integral with fixed chemical potential. Almost all the physics can be discussed in equivalent manners, but there are a few exceptions. An apparent difference is whether the axis of the phase diagram is number density or chemical potential. Along the number density axis, we cannot meet the phenomenon below the critical chemical potential, such as the Silver Blaze phenomenon. Finally, we should comment on another proposal of lattice gauge theory at nonzero density [22], which is based on the Euclidean path integral. We are not ready to say which is better at the moment, so we should try both possibilities.