# Schwinger model at finite temperature and density with $\beta$-VQE

Akio Tomiya<EMAIL_ADDRESS>Department of Information Technology,
International Professional University of Technology, Osaka, 3-3-1, Umeda,
Kita-Ku, Osaka, 530-0001, Japan

###### Abstract

We investigate a quantum gauge theory at finite temperature and density using
a variational algorithm for near-term quantum devices. We adapt $\beta$-VQE to
evaluate thermal and quantum expectation values and study the phase diagram
for massless Schwinger model along with the temperature and density. By
compering the exact variational free energy, we find the variational algorithm
work for $T>0$ and $\mu>0$ for the Schwinger model. No significant volume
dependence of the variational free energy is observed in $\mu/g\in[0,1.4]$. We
calculate the chiral condensate and take the continuum extrapolation. As a
result, we obtain qualitative picture of the phase diagram for massless
Schwinger model.

###### pacs:

††preprint: August 27, 2024

##  I Introduction

Phase diagram of QCD is one of the most attractive subject in particle physics
and nuclear physics [1, 2, 3]. At zero quark chemical potential, lattice QCD
calculations with Markov-chain Monte-Carlo work well and we can access
quantitative information, e.g. pseudo-critical temperature and cumulants [4,
5, 6, 7, 8, 9]. These values are tightly related to the heavy ion collision
experiments and cosmological observations, and it is now more and more
important to investigate QCD phase diagram to understand nature of our
universe.

Phase structure along with finite baryon density cannot be accessed the method
above, due to the infamous sign problem, and a number of proposals are
suggested (see [10, 11, 3, 12] and references therein). The complex Langevin
algorithm has been applied not only matrix models and toy models, but also QCD
in four dimension [13, 14, 15, 16]. Lefschetz thimble is a similar approach,
but it is an exact algorithm [17, 18, 19, 20, 21]. Methods with change of
variable is also investigated [22, 23, 24, 25, 26, 27, 28].

Algorithms for digital quantum computers for lattice QCD are highly demanded
to overcome the sign problem [29]. Quantum algorithms can be regarded as ways
to realize a unitary operation on quantum state based on unitary operations on
qubits. Quantum algorithms for quantum field theory has been applied for the
real time [30, 31] and $\theta$ vacuum [32].

Figure 1:  Density plot of the chiral condensate from $\beta$-VQE along with
$T/g$ – $\mu/g$. The central values in Fig. 11 are interpolated and shown (see
the main text in detail). This density plot can be seen as the phase diagram.

Noisy Intermediate-Scale Quantum Computers (NISQ) [33] have limited capability
than ones with error correction, but still it is useful to investigate quantum
field theory since it enables us to investigate physics at parameters with
sign problem [34]. While it has been investigated purely quantum algorithm
[35], classical-quantum hybrid variational algorithms is considerably suitable
for NISQ devices [36, 37] rather than exact algorithms.

Quantum field theory at finite temperature on digital quantum computers is not
straightforward because the quantum circuit only realizes the unitary
evolution. Conventional thermodynamics uses mixed states, and it cannot be
realized using unitary operations only. Naively, we can combine the unitary
operations and projective measurements to realize the mixed states, but we
loose signal exponentially in the volume. On the other hand, we can eigenstate
thermalization hypothesis [38] or the thermal pure quantum state argument [39]
which can be realized in the unitary operations/quantum circuit but it is not
widely investigated in this context.

A classical quantum hybrid variational algorithm $\beta$-VQE [40] for NISQ
devices enables us to calculate thermal-quantum average without the sign
problem and exponential signal lost. It uses a parametrized probability
distribution to emulate thermal part of the density matrix, and quantum device
or emulator to calculate matrix elements. Moreover, for small system, we can
evaluate the quality of the variational approximation on a classical computer
because it becomes exact if the variational free energy is zero [40]. It has
been applied on the transverse filed Ising model and it agree with the exact
result.

The Schwinger model [41, 42, 43] is a suitable testbed to examine a new
algorithm, which is two dimensional quantum abelian gauge theory with
fermionic matters which has common features to QCD in four dimension. It is
exactly solvable at $m=0$ [41, 44]. For zero temperature in the path integral
formalism, lattice version of the model is investigate to explore new
algorithms with machine learning [45, 46]. For hamiltonian formalism, a tensor
network has been used to investigate finite temperature or finite chemical
potential [47, 48, 49, 50, 51]. Recently it has been employed to examine
feasibility of application of lattice field theory on quantum computers [30,
32, 52, 31, 34].

In this work, we investigate phase structure of the massless Schwinger model
for finite temperature and density toward to investigate QCD phase diagram
using a classical-quantum hybrid algorithm. More precisely, we evaluate
$T$–$\mu$ dependence of the chiral condensate with $\beta$-VQE on a classical
computer. We measure the variational free energy to evaluate quality of
variational thermal states. We evaluate the large volume limit with data from
two large volume and take the continuum limit. As a result, we obtain a phase
diagram in Fig. 1.

This paper is organized as follows. In next section, we introduce the
Schwinger model and its the spin representation. In addition, we review
$\beta$-VQE briefly. Next, we show the results. We discuss the variational
free energy, which quantify the quality of the approach. In next section, we
evaluate the large volume limit and the continuum limit. In next section, we
show the results for the chiral condensate as a function of the coupling
constant and the chemical potential. Finally, we summarize this work.

##  II Model and algorithm

###  II.1 Schwinger model and spin representation

Here we explain one flavor Schwinger model briefly, and please refer [32] for
details. The Lagrangian of the Schwinger model with $\theta=0$ is given by

$\displaystyle\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\mathrm{i}\overline{\psi}\gamma^{\mu}\left(\partial_{\mu}+\mathrm{i}gA_{\mu}\right)\psi-m\overline{\psi}\psi,$
(1)

where $g$ is the dimension-full coupling constant. To write down the
hamiltonian, we have to fix the gauge and we choose the temporal gauge
$A_{0}(x)=0$. The hamiltonian is,

$\displaystyle H$ $\displaystyle=\int
dx\Big{[}-\mathrm{i}\overline{\psi}\gamma^{1}\left(\partial_{1}+\mathrm{i}gA_{1}\right)\psi+m\overline{\psi}\psi+\frac{1}{2}\Pi^{2}\Big{]},$
(2)

where $\Pi\equiv\dot{A}^{1}$. In addition to the Hamiltonian, we need Gauss’
law to make the Hilbert space to keep be physical,

$\displaystyle 0=-\partial_{1}\Pi-g\overline{\psi}\gamma^{0}\psi.$ (3)

We are interested in temperature and chemical potential dependence of thermal
and quantum expectation values of an operator $O$,

$\displaystyle\left\langle
O\right\rangle=\frac{1}{Z}\mathrm{Tr}\,[O\text{e}^{-\frac{1}{T}H}]=\mathrm{Tr}\,[O\rho]$
(4)

where $\rho=\text{e}^{-\frac{1}{T}H}/Z$ and $Z$ is a normalization constant
and $T$ is the temperature. $\mathrm{Tr}\,[\cdots]$ means functional trace
over the Hilbert space.

The Schwinger model enjoys non-zero chiral condensate at the ground state and
also finite temperature $T$ [44],

$\displaystyle\left\langle\overline{\psi}\psi\right\rangle(T)=\Sigma_{0}\text{e}^{I(m_{g}/T)}$
(5)

where $\Sigma_{0}=m_{g}\text{e}^{\gamma}/(2\pi)$ is the chiral condensate at
zero-temperature and $m_{g}=g/\sqrt{\pi}$. The function $I(x)$ is,

$\displaystyle I(x)=\int_{0}^{\infty}\frac{1}{1-e^{x\cosh t}}dt.$ (6)

This has been calculated using the tensor network.

We descritize the system with the staggered formalism,

$\displaystyle H$
$\displaystyle=-\mathrm{i}\sum_{n=1}^{N_{x}-1}w\left[\chi_{n}^{\dagger}e^{\mathrm{i}\phi_{n}}\chi_{n+1}-\text{
h.c. }\right]$
$\displaystyle\;\;\;\;\;\;+\sum_{n=1}^{N_{x}}m(-1)^{n}\chi_{n}^{\dagger}\chi_{n}+J\sum_{n=1}^{N_{x}-1}L_{n}^{2}$
(7)

where $w=1/(2a)$ $J=g^{2}a/2$, and and $N_{x}$ is dimensionless spatial extent
of the staggered fermions.
$L_{n}-L_{n-1}=\chi_{n}^{\dagger}\chi_{n}-({1-(-1)^{n}})/{2}$ is descritized
Gauss’ law and $L_{1}$ chosen by the boundary condition. By imposing the open
boundary condition and apply a gauge transformation, we can eliminate the
gauge field from the Hamiltonian and as the drawback, we obtain a non-local
interaction term.

We can include the number operator $\chi_{n}^{\dagger}\chi_{n}$ with the
chemical potential $\mu$ to the hamiltonian as [53],

$\displaystyle\sum_{n=1}^{N_{x}}(m(-1)^{n})\chi_{n}^{\dagger}\chi_{n}\Rightarrow\sum_{n=1}^{N_{x}}(m(-1)^{n}+\mu)\chi_{n}^{\dagger}\chi_{n}.$
(8)

In order to calculate on digital quantum computers, operators in the
hamiltonian must be written spin operators. The Jordan-Wignar transformation
makes the spin representation of fermions,

$\displaystyle\chi_{n}=\left(\prod_{1\leq
l<n}-\mathrm{i}Z_{l}\right)\frac{X_{n}-\mathrm{i}Y_{n}}{2},$ (9)

where $X_{n}$, $Y_{n}$ and $Z_{n}$ are Pauli matrices for $x,y,z$ direction
acting on $n$-th qubit. This transformation absorbs difference of statistical
property between Pauli matrices and fermion fields.

After the transformation, we get the Hamiltonian for the Schwinger model in
spin operator representation,

$\displaystyle\hat{H}=H/g=\hat{H}_{ZZ}+\hat{H}_{\pm}+\hat{H}_{Z}.$ (10)

where

$\displaystyle\hat{H}_{ZZ}=\frac{g}{8w}\sum_{n=2}^{N_{x}-1}\sum_{1\leq k<l\leq
n}Z_{k}Z_{l},$ (11)
$\displaystyle\hat{H}_{\pm}=\frac{1}{2}\sum_{n=1}^{N_{x}-1}\frac{w}{g}\left[X_{n}X_{n+1}+Y_{n}Y_{n+1}\right],$
(12)
$\displaystyle\hat{H}_{Z}=\frac{1}{2}\sum_{n=1}^{N_{x}}\left(\frac{m}{g}(-1)^{n}+\frac{\mu}{g}\right)Z_{n}$
$\displaystyle\;\;\;\;\;\;-\frac{g}{8w}\sum_{n=1}^{N_{x}-1}(n\bmod
2)\sum_{l=1}^{n}Z_{l}$ (13)

and $\hat{\cdot}$ indicates dimensionless quantity. The exponent of Boltzmann
weight is now $H/T=\hat{H}g/T$, and $g/T$ is dimensionless inverse
temperature, which conventionally referred as $\beta$ but we do not use this
notation to avoid confusion with $1/g^{2}$.

In this study, we do not use the adiabatic state preparation [32], instead, we
use a classical-quantum hybrid algorithm $\beta$-VQE, to evaluate thermal and
quantum expectation values. We focus on $m=0$ in the numerical calculations,
which does not have additive divergence of the chiral condensate.

###  II.2 $\beta$-VQE: Neural network parametrized ansatz

In this work, we adapt $\beta$-VQE to evaluate quantum and thermal expectation
value [40]. Here we briefly review $\beta$-VQE based on the original paper.
VQE stands for variational quantum eigen-solver, which enables us to evaluate
quantum expectation values with a parametrized shallow quantum circuit. The
shallow quantum circuit is parametrized using rotation gates, which rotates a
qubit in SU(2) space. The initial state is taken to classically tractable
product state and a shallow quantum circuit produces entanglement. The
parameters in the circuit is tuned using classical process to minimize a loss
function including quantum expectation value which is evaluated quantum
computer or quantum circuit emulator. We do not pursue exactness since the
variational method gives precise results based on the variational principle
but is not an exact method in practice, and results must be qualitative.

We consider a $N$-qubit system represented by a bit string $x=x_{1}x_{2}\cdots
x_{N}$ where $x_{i}\in\\{0,1\\}$ ($i=1,\cdots,N$). In particular we employ a
product state,

$\displaystyle\left|x\right\rangle=\left|x_{1}\right\rangle\otimes\left|x_{2}\right\rangle\otimes\cdots\otimes\left|x_{N}\right\rangle.$
(14)

We call this as the classical product state, which is easy to prepare.
Variational approaches use a parametrized state
$U_{\theta}\left|x\right\rangle$, where $\theta$ is a set of parameters for a
quantum state and $U_{\theta}$ a unitary evolution like rotation gates in
quantum circuits. In practice, we use SU(4) parametrization for $U_{\theta}$
as in the original work [40].

We parametrize the density matrix as,

$\displaystyle\tilde{\rho}_{\Theta}=\sum_{x}p_{\phi}(x)U_{\theta}|x\rangle\langle
x|U_{\theta}^{\dagger},$ (15)

where $\Theta=\theta{}^{\cup}\phi$ and $p_{\phi}$ is a parametrized classical
distribution. This density matrix satisfies the normalization condition
$\mathrm{Tr}\,[\tilde{\rho}_{\Theta}]=\sum_{x}p_{\phi}(x)=1$. In the
application, we employ the autoregressive model [54], which is a neural
network. As a unit of the autoregressive model, we utilize an autoencoder,
which has single hidden layer of 500 hidden neurons with rectified linear unit
activation function. This is also the same setup from the original work for
$\beta$-VQE [40]. They show that, $\beta$-VQE works well for two dimensional
transverse field Ising model.

Our target system is written as $\rho=\text{e}^{-\frac{1}{T}H}/Z$ and we
minimize following loss function,

$\displaystyle\mathcal{L}(\Theta)=\overline{D}(\tilde{\rho}_{\Theta}||\rho)=\operatorname{Tr}\left(\tilde{\rho}_{\Theta}\ln\left(\frac{\tilde{\rho}_{\Theta}}{\rho}\right)\right),$
(16)

by tuning parameters $\Theta$. This variational argument is developed in the
Gibbs-Delbruck-Moliere variational principle of quantum statistical mechanics
LABEL:huber1968variational We remark that, this loss function can be seen as,
quantum version of (reversed) Kullback Leibler divergence (an application for
lattice field theory with the classical version can be found in [55, 56, 57]
). The Kullback Leibler divergence for classical probability distributions
$p,q$ is given by,

$\displaystyle{D}(p||q)=\int dx\;p(x)\ln\frac{p(x)}{q(x)},$ (17)

which is referred as the relative entropy in physics context. As same as the
classical Kullback Leibler divergence ${D}(p||q)$, this
$\overline{D}(\tilde{\rho}_{\Theta}||\rho)$ is zero if and only if
$\tilde{\rho}_{\Theta}=\rho$ [40]. This means that, we can quantify the
quality of variational solution through $\mathcal{L}(\Theta)$.

On the real quantum devices as intended in the original paper, we cannot
calculate the partition function $Z$, or equivalently the free energy $-\ln
Z$. In that case, we have to use reversed and shifted loss function,

$\displaystyle\mathcal{L}(\Theta)-\ln Z$
$\displaystyle=\operatorname{Tr}(\tilde{\rho}_{\Theta}\ln\tilde{\rho}_{\Theta})+\frac{1}{T}\operatorname{Tr}(\tilde{\rho}_{\Theta}H).$
(18)

We remark that no bias property is still held. Namely, this loss function is
minimized if and only if $\tilde{\rho}_{\Theta}=\rho$. In addition, the
constant shift does not affect on the derivative and training and the right
hand side is calculable because it is free from the normalization constant.
However in this case, we cannot see that the parameter $\Theta$ is optimal or
not. In this pilot study, we calculate
$\mathcal{L}(\Theta)=\hat{D}(\tilde{\rho}_{\Theta}||\rho)$, not the shifted
one, using state vector to examine the algorithm and observe the phase
structure.

Since it can produce training samples by itself, we rely on the self-training
paradigm for the training We sample from the model distribution as, in the
training process,

$\displaystyle{\mathcal{L}}(\Theta)=\mathbb{E}_{x\sim p_{\phi}(x)}\left[\ln
p_{\phi}(x)+\frac{1}{T}\left\langle
x\left|U_{\theta}^{\dagger}HU_{\theta}\right|x\right\rangle\right]+{\rm
const}.$ (19)

Optimization of the loss function ${\mathcal{L}}(\Theta)$ is performed by a
stochastic gradient optimizer (e.g., adam [58]). We have two kinds of
derivative. One is a derivative with respect to parameters in the quantum
circuit,

$\displaystyle\nabla_{\theta}{\mathcal{L}}(\Theta)=\frac{1}{T}\mathbb{E}_{x\sim
p_{\phi}(x)}\left[\nabla_{\theta}\left\langle
x\left|U_{\theta}^{\dagger}HU_{\theta}\right|x\right\rangle\right].$ (20)

This derivative can be evaluated by the shift rule [59, 60].

The other is a derivative with respect to parameters in the classical
distribution,

$\displaystyle\nabla_{\phi}{\mathcal{L}}(\Theta)=\mathbb{E}_{x\sim
p_{\phi}(x)}\left[(f(x)-b)\nabla_{\phi}\ln p_{\phi}(x)\right],$ (21)

where

$\displaystyle f(x)=\ln p_{\phi}(x)+\frac{1}{T}\left\langle
x\left|U_{\theta}^{\dagger}HU_{\theta}\right|x\right\rangle$ (22)

and $b=\mathbb{E}_{x\sim p_{\phi}}[f(x)]$. This term $b$ is not essentially
needed but it helps to reduce the variance of the optimizer [40].

After training, one can sample a batch of input states $\left|x\right\rangle$
and treat them as approximations of the eigenstates of the system. Since the
unitary circuit preserves orthogonality of the input states, the sampled
quantum states span a low energy subspace of the Hamiltonian.

First we stochastically and iteratively find $\Theta^{*}=\underset{\Theta}{\rm
argmin}\mathcal{L}(\Theta)$. By using $\Theta^{*}$ and (15), we can evaluate
the expectation values with sampling $x$ from $p_{\phi}$. We evaluate
spatially averaged chiral condensate [32],

$\displaystyle\left\langle\overline{\psi}\psi\right\rangle_{\Theta^{*}}\big{/}g=\frac{2w}{N_{x}}\sum_{n}\mathrm{Tr}\,[\overline{\psi}\psi(n)\tilde{\rho}_{\Theta^{*}}].$
(23)

Note that, states, in particular for finite density, might not have
translational invariance [61, 62, 63] but we use this operator as an
indicator.

To evaluate the quantum expectation value, we utilize calculations with state
vectors on a classical computer to avoid uncertainty of circuit operations.
Use of classical computer limits the maximum number of sites in practice but
it enables us to evaluate exact free energy. In addition, the replacement is
straightforward except for error evaluation of the variational method [40].

##  III Results

We perform $\beta$-VQE calculations for following parameters:
$N_{x}\in\\{4,6,8,10\\}$ (dimensionless spatial extent of the spin model),
$w/g\in\\{0.5,0.45,0.4,0.35\\}$ (lattice spacing),
$g/T\in\\{0.1,0.5,1.0,2.0,4.0,10.0,20.0\\}$ (inverse temperature),
$\mu/g=\\{0.0,0.2,0.4,0.6,1.0,1.4\\}$ (chemical potential). Our temperature
range is $T/g\in[0.05,10]$. As we mentioned, we store a state vector in this
study to get quantum expectation value in each step in the training, and we
cannot take the volume more than $N_{x}=10$ in practice.

The number of training is 500 for each $(N_{x}$, $w/g$, $g/T$, $\mu/g)$ and we
employ the adam optimizer to minimize the loss function. During the training,
we monitor the variational free energy to quantify the quality of the
variational density matrix. The number of sample for thermal average is taken
to 5000.

In the analysis of the chiral condensate, we take the continuum extrapolation
($g/(2w)=ag\to 0$) for our largest lattice data. As we will comment, we do not
include data from $N_{x}=4,6$ in the analysis for continuum physics to avoid
finite volume effects. Detailed analysis with physically larger lattice using
other method is left for another publication.

The calculations are carried out using an open source package $\beta$-VQE.jl
[64], which is based on Yao.jl [65], a package for quantum calculations, and
Zygote.jl [66], a package for automatic differentiation in neural networks, in
Julia language [67].

###  III.1 Variational free energy in $\beta$-VQE

Here we discus the variational free energy in $\beta$-VQE,
$\mathcal{L}(\Theta)$.

Tab 1 is the results of the variational free energy for our smallest volume
$N_{x}=4$ and largest volume $N_{x}=10$ with the cutoff $w=0.5$ and $0.35$ and
the temperature $g/T=0.1$ to $20.0$. The variational free energy defined in
(18) in the end of the training at epoch 500 are listed. A column Diff is
$(\mathcal{L}/(-\ln Z))\times 100$, deviation from the exact free energy. As
we mentioned above, $\mathcal{L}$ is zero if and only if the variational
approximation is exact. Except for $g/T=20,\mu/g=0$, $N_{x}$ = 4, difference
of loss function to the exact free energy density takes 0 in $O(1)$ %. In our
parameter regime, we do not observe strong dependence on the chemical
potential. Rather, we observe strong dependence on the temperature, as
expected.

We plot deviation of the variational energy from the exact one as a function
epoch in Fig. 2, Fig. 3, Fig. 4 and Fig. 5 for $N_{x}=10$, $g/T=0.1,20$ and
$\mu/g=0.0,1.4$. In each plot, blue, orange, green, red lines correspond to
$w=1/(2a)=0.5,0.45,0.4$ and $0.35$, respectively. One can see that, the
deviation is $O(1)$ % after 500 epoch. This also indicates that, the error
from the variational method mainly depends on the temperature rather than the
chemical potential in our parameter regime.

$\mu/g$ | $g/T$ | $N_{x}$ | $w/g$ | $\mathcal{L}-\ln Z$ | $-\ln Z$ | Diff (%)
---|---|---|---|---|---|---
0.0 | 0.1 | 4 | 0.5 | -27.779 | -27.781 | 0.00804
0.0 | 0.1 | 4 | 0.35 | -27.807 | -27.808 | 0.005
0.0 | 0.1 | 10 | 0.5 | -70.686 | -70.718 | 0.0459
0.0 | 0.1 | 10 | 0.35 | -71.744 | -71.765 | 0.0302
0.0 | 0.5 | 4 | 0.5 | -5.792 | -5.802 | 0.185
0.0 | 0.5 | 4 | 0.35 | -5.885 | -5.891 | 0.105
0.0 | 0.5 | 10 | 0.5 | -17.133 | -17.25 | 0.68
0.0 | 0.5 | 10 | 0.35 | -18.849 | -18.934 | 0.448
0.0 | 10.0 | 4 | 0.5 | -1.748 | -1.75 | 0.161
0.0 | 10.0 | 4 | 0.35 | -1.829 | -1.829 | 0.0184
0.0 | 10.0 | 10 | 0.5 | -8.218 | -8.341 | 1.48
0.0 | 10.0 | 10 | 0.35 | -9.98 | -10.03 | 0.496
0.0 | 20.0 | 4 | 0.5 | -1.492 | -1.739 | 14.2
0.0 | 20.0 | 4 | 0.35 | -1.653 | -1.806 | 8.46
0.0 | 20.0 | 10 | 0.5 | -8.202 | -8.328 | 1.51
0.0 | 20.0 | 10 | 0.35 | -9.955 | -10.006 | 0.509
1.4 | 0.1 | 4 | 0.5 | -28.021 | -28.023 | 0.00697
1.4 | 0.1 | 4 | 0.35 | -27.989 | -27.991 | 0.00755
1.4 | 0.1 | 10 | 0.5 | -70.842 | -70.874 | 0.0453
1.4 | 0.1 | 10 | 0.35 | -71.742 | -71.763 | 0.0291
1.4 | 0.5 | 4 | 0.5 | -6.784 | -6.789 | 0.0609
1.4 | 0.5 | 4 | 0.35 | -6.644 | -6.647 | 0.0327
1.4 | 0.5 | 10 | 0.5 | -17.989 | -18.104 | 0.636
1.4 | 0.5 | 10 | 0.35 | -19.445 | -19.534 | 0.456
1.4 | 10.0 | 4 | 0.5 | -3.708 | -3.71 | 0.0728
1.4 | 10.0 | 4 | 0.35 | -3.63 | -3.669 | 1.07
1.4 | 10.0 | 10 | 0.5 | -10.067 | -10.243 | 1.71
1.4 | 10.0 | 10 | 0.35 | -11.763 | -11.862 | 0.837
1.4 | 20.0 | 4 | 0.5 | -3.673 | -3.681 | 0.218
1.4 | 20.0 | 4 | 0.35 | -3.621 | -3.669 | 1.31
1.4 | 20.0 | 10 | 0.5 | -10.028 | -10.224 | 1.92
1.4 | 20.0 | 10 | 0.35 | -11.699 | -11.862 | 1.37
Table 1: Variational free energy defined in (18) in the end of the training at
epoch 500. Diff is $(\mathcal{L}/(-\ln Z))\times 100$, deviation from the
exact free energy. $\mathcal{L}$ is zero iff the variational approximation is
exact.  Figure 2:  Deviation of the variational energy from the exact one as a
function epoch for $N_{x}=10$, $g/T=0.1$ and $\mu/g=0.0$. Blue, orange, green,
red lines correspond to $w=1/(2a)=0.5,0.45,0.4$ and $0.35$, respectively.
Figure 3:  Same figure of Fig. 2 but for $g/T=20$ and $\mu/g=0.0$.  Figure 4:
Same figure of Fig. 2 but for $g/T=0.1$ and $\mu/g=1.4$.  Figure 5:  Same
figure of Fig. 2 but for $g/T=20$ and $\mu/g=1.4$.

In the later analysis, we do not include systematic error from the variation.

###  III.2 Evaluation of large volume limit, and continuum extrapolation

In the previous work [32], we observed that data from $N_{x}>10$ are necessary
to take stable large volume limit. However in this work, since we use state
vector calculations for each variation step, and we cannot obtain these data
in practice. We just use data from $N_{x}=8,10$ and evaluate the large volume
limit.

We extrapolate the data to $a\to 0$ limit using linear fit ansatz and results
are shown in Fig. 6 and Fig. 7 for $N_{x}=8$ for $\mu/g=0$ and $\mu/g=1.4$.
For $N_{x}=10$, we plot in Fig. 8 and Fig. 9. The tendency of data is nearly
linear, thus, this limit captures scaling to the continuum physics.

Figure 6:  Continuum limit for $N_{x}=8$ and $\mu/g=0.0$. Different color
corresponds to different temperature. The error bar only contains statistical
error.  Figure 7:  Same figure with Fig 6 but for $\mu/g=1.4$.  Figure 8:
Same figure with Fig 6 but for $N_{x}=10$.  Figure 9:  Same figure with Fig 6
but for $N_{x}=10$ and $\mu/g=1.4$.

###  III.3 Phase diagram for $T$ – $\mu$

Here we discuss phase diagram of the Schwinger model through results for the
chiral condensates. We summarize the chiral condensate as a function of
$\beta$ at $a\to 0$ in Fig. 10 and Fig. 11 for $N_{x}=8$ and 10, respectively.
To avoid overlap, we shift the horizontal axis and the vertical axis is
normalized with the value in the continuum at $m=0$ with $g/T=20.0$,
$\mu/g=0.0$ for each $N_{x}$ for our data. Error bar represents statistical
error but error from the variation is not included. The error from the
variation is estimated as $O(1)$ % in the original data, so if we improve the
variation, we expect it will not be changed very much. Our results for
$\mu/g=0$ is qualitatively agree with the exact results [44] in whole
temperature range. One can observe that, the chiral condensate gradually
decreases along with $\mu/g$ for low temperature regime. This has to be
confirmed with larger lattice.

Density plot of the chiral condensate from $\beta$-VQE along with $T/g$ –
$\mu/g$ is shown in Fig. 1. The central values in Fig. 11 are interpolated and
plotted. The results for the chiral condensate are normalized with the value
in the continuum at $m=0$ with $g/T=20.0$, $\mu/g=0.0$ for each $N_{x}$. The
massless Schwinger model does not have distinct phases but we phrase a
parameter regime with the normalized chiral condensate larger than 1/2, as
broken phase otherwise symmetric phase for showing purpose in the plot. One
can see that, this qualitatively agree with the what suggested for QCD [1, 2]
for small chemical potential regime.

Figure 10:  Chiral condensate as a function of $\beta$ at $a\to 0$ for
$N_{x}=8$. To avoid overlap, we shift the horizontal axis. Error bar
represents statistical error (same size with the symbols) but error from the
variation is not included. The results for the chiral condensate are
normalized with the value in the continuum at $m=0$ with $g/T=20.0$,
$\mu/g=0.0$ at $N_{x}=8$.  Figure 11:  Same figure of Fig. 10 but for
$N_{x}=10$.

##  IV Summary and discussion

In this work, we investigate phase diagram on finite temperature and density
for the Schwinger model using a quantum classical hybrid algorithm called
$\beta$-VQE, which is not affected by the sign problem. Quantum expectation
values are evaluated through state vector calculations on a classical computer
which should be replaced by quantum calculations in the future. Thanks to the
state vector calculations, we evaluate the exact free energy $-\ln Z$ and we
confirm that the variational algorithm gives $O(1)\%$ correct results.

Only continuum limit is taken and large volume limit has not been taken in
this work. While, we observe qualitative agreement along with the temperature
to the exact results for $\mu/g=0$ [44] and the deviation is similar to [48].
To establish physics at infinity volume, we should replace the state-vector
calculations to a tensor network or quantum device to overcome the shortage of
the numerical resource. We enphasize that $\beta$-VQE can be used with a
quantum device, and once we get a trustable quantum machine, fault-tolerant
quantum computer, we can easily apply this calculation on it.

We have not investigated crystal structure for $\mu/g>0$ which is predicted in
several studies [61, 62, 63]. By defining position dependent gauge invariant
correlators, and use a large number of sampling, we could see the crystal
structure and dependence on the temperature.

We have used the autoregressive model [54] to realize the probability
distribution $p_{\phi}(x)$ in the density matrix. It could be replaced by
other neural net model like the normalizing flow [56]. Use of QITE [35] is
also interested.

## Acknowledgement

AT thanks to Lei Wang and Jinguo Liu for helping me developing the code.
Numerical computation in this work was carried out at the Yukawa Institute
Computer Facility. The work of AT was supported by JSPS KAKENHI Grant No.
JP20K14479.

## References

  * [1] N. Cabibbo and G. Parisi.  Exponential hadronic spectrum and quark liberation.  Physics Letters B, 59(1):67–69, 1975.
  * [2] Kenji Fukushima and Tetsuo Hatsuda.  The phase diagram of dense QCD.  Rept. Prog. Phys., 74:014001, 2011.
  * [3] Jana N. Guenther.  Overview of the QCD phase diagram: Recent progress from the lattice.  Eur. Phys. J. A, 57(4):136, 2021.
  * [4] A. Bazavov et al.  Fluctuations and Correlations of net baryon number, electric charge, and strangeness: A comparison of lattice QCD results with the hadron resonance gas model.  Phys. Rev. D, 86:034509, 2012.
  * [5] Heng-Tong Ding, Frithjof Karsch, and Swagato Mukherjee.  Thermodynamics of strong-interaction matter from Lattice QCD.  Int. J. Mod. Phys. E, 24(10):1530007, 2015.
  * [6] A. Bazavov et al.  Chiral crossover in QCD at zero and non-zero chemical potentials.  Phys. Lett. B, 795:15–21, 2019.
  * [7] Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Sandor K. Katz, Kalman K. Szabo, Attila Pasztor, Israel Portillo, and Claudia Ratti.  Higher order fluctuations and correlations of conserved charges from lattice QCD.  JHEP, 10:205, 2018.
  * [8] Frithjof Karsch.  Critical behavior and net-charge fluctuations from lattice QCD.  PoS, CORFU2018:163, 2019.
  * [9] Jishnu Goswami, Frithjof Karsch, Christian Schmidt, Swagato Mukherjee, and Peter Petreczky.  Comparing conserved charge fluctuations from lattice QCD to HRG model calculations.  Acta Phys. Polon. Supp., 14:251, 2021.
  * [10] Christof Gattringer, Thomas Kloiber, and Vasily Sazonov.  Solving the sign problems of the massless lattice Schwinger model with a dual formulation.  Nucl. Phys. B, 897:732–748, 2015.
  * [11] Claudia Ratti.  QCD at non-zero density and phenomenology.  PoS, LATTICE2018:004, 2019.
  * [12] Andrei Alexandru, Gokce Basar, Paulo F. Bedaque, and Neill C. Warrington.  Complex paths around the sign problem.  Rev. Mod. Phys., 94(1):015006, 2022.
  * [13] Gert Aarts, Erhard Seiler, Denes Sexty, and Ion-Olimpiu Stamatescu.  Complex Langevin dynamics and zeroes of the fermion determinant.  JHEP, 05:044, 2017.  [Erratum: JHEP 01, 128 (2018)].
  * [14] Erhard Seiler.  Status of Complex Langevin.  EPJ Web Conf., 175:01019, 2018.
  * [15] Casey E. Berger, Lukas Rammelmüller, Andrew C. Loheac, Florian Ehmann, Jens Braun, and Joaquín E. Drut.  Complex Langevin and other approaches to the sign problem in quantum many-body physics.  Phys. Rept., 892:1–54, 2021.
  * [16] M. Scherzer, D. Sexty, and I. O. Stamatescu.  Deconfinement transition line with the complex Langevin equation up to $\mu/T\sim 5$.  Phys. Rev. D, 102(1):014515, 2020.
  * [17] Edward Witten.  Analytic Continuation Of Chern-Simons Theory.  AMS/IP Stud. Adv. Math., 50:347–446, 2011.
  * [18] Yuto Mori, Kouji Kashiwa, and Akira Ohnishi.  Lefschetz thimbles in fermionic effective models with repulsive vector-field.  Phys. Lett. B, 781:688–693, 2018.
  * [19] Christian Schmidt and Felix Ziesché.  Simulating low dimensional QCD with Lefschetz thimbles.  PoS, LATTICE2016:076, 2017.
  * [20] Genki Fujisawa, Jun Nishimura, Katsuta Sakai, and Atis Yosprakob.  Backpropagating Hybrid Monte Carlo algorithm for fast Lefschetz thimble calculations.  JHEP, 04:179, 2022.
  * [21] Masafumi Fukuma, Nobuyuki Matsumoto, and Yusuke Namekawa.  Numerical sign problem and the tempered lefschetz thimble method.  In 21st Hellenic School and Workshops on Elementary Particle Physics and Gravity, 5 2022.
  * [22] Yuto Mori, Kouji Kashiwa, and Akira Ohnishi.  Application of a neural network to the sign problem via the path optimization method.  PTEP, 2018(2):023B04, 2018.
  * [23] Andrei Alexandru, Paulo F. Bedaque, Henry Lamm, and Scott Lawrence.  Deep Learning Beyond Lefschetz Thimbles.  Phys. Rev. D, 96(9):094505, 2017.
  * [24] Kouji Kashiwa, Yuto Mori, and Akira Ohnishi.  Controlling the model sign problem via the path optimization method: Monte Carlo approach to a QCD effective model with Polyakov loop.  Phys. Rev. D, 99(1):014033, 2019.
  * [25] Scott Lawrence and Yukari Yamauchi.  Normalizing Flows and the Real-Time Sign Problem.  Phys. Rev. D, 103(11):114509, 2021.
  * [26] William Detmold, Gurtej Kanwar, Henry Lamm, Michael L. Wagman, and Neill C. Warrington.  Path integral contour deformations for observables in $SU(N)$ gauge theory.  Phys. Rev. D, 103(9):094517, 2021.
  * [27] Gurtej Kanwar and Michael L. Wagman.  Real-time lattice gauge theory actions: Unitarity, convergence, and path integral contour deformations.  Phys. Rev. D, 104(1):014513, 2021.
  * [28] Yusuke Namekawa, Kouji Kashiwa, Akira Ohnishi, and Hayato Takase.  Gauge invariant input to neural network for path optimization method.  Phys. Rev. D, 105(3):034502, 2022.
  * [29] Scott Lawrence.  Sign Problems in Quantum Field Theory: Classical and Quantum Approaches.  PhD thesis, Maryland U., 2020.
  * [30] E. A. Martinez et al.  Real-time dynamics of lattice gauge theories with a few-qubit quantum computer.  Nature, 534:516–519, 2016.
  * [31] Emilie Huffman, Miguel García Vera, and Debasish Banerjee.  Real-time dynamics of Plaquette Models using NISQ Hardware.  9 2021.
  * [32] Bipasha Chakraborty, Masazumi Honda, Taku Izubuchi, Yuta Kikuchi, and Akio Tomiya.  Classically emulated digital quantum simulation of the Schwinger model with a topological term via adiabatic state preparation.  Phys. Rev. D, 105(9):094503, 2022.
  * [33] Kishor Bharti et al.  Noisy intermediate-scale quantum algorithms.  Rev. Mod. Phys., 94(1):015004, 2022.
  * [34] Arata Yamamoto.  Quantum variational approach to lattice gauge theory at nonzero density.  Phys. Rev. D, 104(1):014506, 2021.
  * [35] Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandao, and Garnet Kin Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  Nature Physics, 16(2):205–210, 2020.
  * [36] Jin-Guo Liu, Yi-Hong Zhang, Yuan Wan, and Lei Wang.  Variational quantum eigensolver with fewer qubits.  Physical Review Research, 1(2), sep 2019.
  * [37] Anthony N. Ciavarella and Ivan A. Chernyshev.  Preparation of the SU(3) lattice Yang-Mills vacuum with variational quantum methods.  Phys. Rev. D, 105(7):074504, 2022.
  * [38] Joshua M Deutsch.  Eigenstate thermalization hypothesis.  Reports on Progress in Physics, 81(8):082001, jul 2018.
  * [39] Sho Sugiura and Akira Shimizu.  Thermal Pure Quantum States at Finite Temperature.  Phys. Rev. Lett., 108:240401, 2012.
  * [40] Jin-Guo Liu, Liang Mao, Pan Zhang, and Lei Wang.  Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits.  Machine Learning: Science and Technology, 2(2):025011, 2021.
  * [41] Julian S. Schwinger.  Gauge Invariance and Mass. 2.  Phys. Rev., 128:2425–2429, 1962.
  * [42] Namik K. Pak and P. Senjanovic.  Axial Anomaly and Chiral Symmetry Breaking in Schwinger Model: A Canonical Approach.  Phys. Lett. B, 71:333–336, 1977.
  * [43] K. Johnson.  gamma(5) invariance.  Phys. Lett., 5:253–255, 1963.
  * [44] Ivo Sachs and Andreas Wipf.  Finite temperature Schwinger model.  Helv. Phys. Acta, 65:652–678, 1992.
  * [45] Jacob Finkenrath.  Tackling critical slowing down using global correction steps with equivariant flows: the case of the Schwinger model.  1 2022.
  * [46] Michael S. Albergo, Denis Boyda, Kyle Cranmer, Daniel C. Hackett, Gurtej Kanwar, Sébastien Racanière, Danilo J. Rezende, Fernando Romero-López, Phiala E. Shanahan, and Julian M. Urban.  Flow-based sampling in the lattice schwinger model at criticality.  2 2022.
  * [47] M. C. Bañuls, K. Cichy, K. Jansen, and J. I. Cirac.  The mass spectrum of the Schwinger model with Matrix Product States.  JHEP, 11:158, 2013.
  * [48] Mari Carmen Bañuls, Krzysztof Cichy, Karl Jansen, and Hana Saito.  Chiral condensate in the Schwinger model with Matrix Product Operators.  Phys. Rev. D, 93(9):094512, 2016.
  * [49] Lena Funcke, Karl Jansen, and Stefan Kühn.  Topological vacuum structure of the Schwinger model with matrix product states.  Phys. Rev. D, 101(5):054507, 2020.
  * [50] Nouman Butt, Simon Catterall, Yannick Meurice, Ryo Sakai, and Judah Unmuth-Yockey.  Tensor network formulation of the massless schwinger model with staggered fermions.  Phys. Rev. D, 101(9):094509, 2020.
  * [51] Masafumi Fukuma, Daisuke Kadoh, and Nobuyuki Matsumoto.  Tensor network approach to two-dimensional yang-mills theories.  PTEP, 2021(12):123B03, 2021.
  * [52] Alexander F. Shaw, Pavel Lougovski, Jesse R. Stryker, and Nathan Wiebe.  Quantum Algorithms for Simulating the Lattice Schwinger Model.  Quantum, 4:306, 2020.
  * [53] Mari Carmen Bañuls, Krzysztof Cichy, J. Ignacio Cirac, Karl Jansen, and Stefan Kühn.  Density Induced Phase Transitions in the Schwinger Model: A Study with Matrix Product States.  Phys. Rev. Lett., 118(7):071601, 2017.
  * [54] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.  MADE: masked autoencoder for distribution estimation.  CoRR, abs/1502.03509, 2015.
  * [55] Gurtej Kanwar, Michael S. Albergo, Denis Boyda, Kyle Cranmer, Daniel C. Hackett, Sébastien Racanière, Danilo Jimenez Rezende, and Phiala E. Shanahan.  Equivariant flow-based sampling for lattice gauge theory.  Phys. Rev. Lett., 125(12):121601, 2020.
  * [56] M. S. Albergo, G. Kanwar, and P. E. Shanahan.  Flow-based generative models for Markov chain Monte Carlo in lattice field theory.  Phys. Rev. D, 100(3):034515, 2019.
  * [57] Denis Boyda, Gurtej Kanwar, Sébastien Racanière, Danilo Jimenez Rezende, Michael S. Albergo, Kyle Cranmer, Daniel C. Hackett, and Phiala E. Shanahan.  Sampling using $SU(N)$ gauge equivariant flows.  Phys. Rev. D, 103(7):074504, 2021.
  * [58] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.
  * [59] Gavin E Crooks.  Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition, 2019.
  * [60] Leonardo Banchi and Gavin E. Crooks.  Measuring analytic gradients of general quantum evolution with the stochastic parameter shift rule.  Quantum, 5:386, jan 2021.
  * [61] Max A. Metlitski.  Is Schwinger model at finite density a crystal?  Phys. Rev. D, 75:045004, 2007.
  * [62] Shinji Maedan.  Influence of current mass on the spatially inhomogeneous chiral condensate.  Prog. Theor. Phys., 123:285–302, 2010.
  * [63] R. Narayanan.  Two flavor massless Schwinger model on a torus at a finite chemical potential.  Phys. Rev. D, 86:125008, 2012.
  * [64] BetaVQE.jl.  https://github.com/wangleiphy/BetaVQE.jl.
  * [65] Xiu-Zhe Luo, Jin-Guo Liu, Pan Zhang, and Lei Wang.  Yao.jl: Extensible, Efficient Framework for Quantum Algorithm Design.  Quantum, 4:341, October 2020.
  * [66] Mike Innes, Alan Edelman, Keno Fischer, Chris Rackauckas, Elliot Saba, Viral B Shah, and Will Tebbutt.  A differentiable programming system to bridge machine learning and scientific computing, 2019.
  * [67] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah.  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98, 2017.