1 Introduction↩︎

Quantum simulation of quantum field theories at finite-temperature regime is a challenging endeavor. This complexity is mainly due to the nature of thermal state as a mixed state, which complicates its preparation on quantum circuit. A possible approach to this problem is usage of a “typical pure state”, the so-called Thermal Pure Quantum (TPQ) state [1], [2]. In thermal equilibrium system, the TPQ state is defined as one whose expectation value approximates the thermal ensemble average of local observables at finite volume. It is proven that the expectation value calculated with the TPQ state converges to the thermal ensemble average in the thermodynamic limit [1]. Although this theoretical framework was initially introduced as a novel formulation of statistical mechanics, it also has the potential for use in the quantum simulation of gauge theories. For example, its applicability to \(Z_2\) gauge theory has been discussed in Ref. [3].

We use a specific realization of the TPQ state termed as the canonical TPQ state [2]. This state is defined by applying the imaginary time evolution operator to random state. Implementing imaginary time evolution on quantum circuit is also challenging since it is expressed by non-unitarity operators. Several algorithms have been proposed, mainly to address the ground state energy problem, such as variational methods [4], probabilistic methods [5], power methods [6], and the Quantum Imaginary Time Evolution (QITE) algorithm [7]. Here, we employ the QITE algorithm combining the random-state generation algorithm to generate the TPQ state.

In this work, we simulate the temperature-dependence of the chiral condensate of the Schwinger model using the TPQ state and the QITE algorithm. We perform the massless case to confirm the validity of our method by comparing our results with the analytic result. Furthermore, we explore the massive and non-zero \(\theta\) regime where the sign problem appears in the conventional Monte Carlo method.

2 The Schwinger model↩︎

Let us first briefly overview the Schwinger model. The Schwinger model is a 1+1 dimensional U(1) gauge theory. The qubit description of the Schwinger model Hamiltonian [8] is given by \[H=\frac{1}{4 a} \sum_{n=1}^{N-1}\left[X_{n} X_{n+1}+Y_{n} Y_{n+1}\right]+\frac{m}{2} \sum_{n=1}^{N}(-1)^{n} Z_{n}+\frac{a g^{2}}{2} \sum_{n=1}^{N-1}\left[\sum_{i=1}^{n-1} \frac{Z_{i}+(-1)^{i}}{2}+\frac{\theta}{2\pi}\right]^{2},\label{eq:hamiltonian95spin}\tag{1}\] where \(N,a\) and \(g\) denote the number of lattice site, the lattice spacing, and the coupling constant, respectively. \(X_i, Y_i, Z_i\) are the Pauli matrices. Here, we solve the Gauss’s law constraint and impose the open boundary condition to remove the degree of freedom for the U(\(1\)) gauge field. Note that the non-local terms are induced due to the Gauss’s law constraint in the last term in Eq. 1 . The topological \(\theta\)-term is shown as a constant shift of the electric-field (\(Z_i\)), and thus no additional difficulty emerges by the insertion in this formula.

In this work, we investigate the (discrete) chiral symmetry as a function of temperature. The order parameter is given by the chiral condensate, \(\langle\bar{\psi}\psi\rangle\), where \(\psi\) is a fermionic field operator. The quantity can be also expressed by the Pauli matrices using the Jordan-Wigner transformation as follows: \[\langle\bar{\psi} \psi\rangle=\frac{1}{N} \sum_{i=1}^{N} (-1)^{i}\frac{1+\left\langle Z_{i}\right\rangle}{2a}. \label{eq:chiral95condensate}\tag{2}\] In the massless case, its temperature-dependence has been derived analytically based on the Lagrangian formalism  [9] as \[\langle\bar{\psi} \psi\rangle=-\frac{m_{\gamma}}{2 \pi} e^{\gamma} e^{2 I\left(\beta m_{\gamma}\right)} \quad,\quad I(x)=\int_{0}^{\infty} \frac{1}{1-e^{x \cosh (t)}} d t ,\label{eq:exact95SW}\tag{3}\] where \(\gamma\) is the Euler constant and \(m_{\gamma}\) is the mass gap \(m_\gamma =g / \sqrt{\pi}\). Note that at zero temperature, it takes the value \(\langle\bar{\psi} \psi\rangle \simeq-0.1599\) [10]. On the other hand, it is difficult to solve this model analytically in massive regime. Furthermore, when dealing with non-zero \(\theta\) regime, the conventional lattice Monte Carlo method is not doable due to the sign problem.

3 Calculation strategy↩︎

To express a thermal state, we prepare a TPQ state [1] on a quantum circuit. The TPQ state, \(|\psi \rangle^{\mathrm{TPQ}}\), is defined as a state that satisfies for low-degree polynomials of local operators, \(^{\forall}\hat{O}\), and \(\forall\epsilon>0\), \[\mathrm{P}\left(\left|\frac{\langle \psi | \hat{O}| \psi \rangle^{\mathrm{TPQ}}}{\langle \psi | \psi \rangle^{\mathrm{TPQ}} }-\langle\hat{O}\rangle^{\mathrm{ens}}\right| \geq \epsilon\right) \leq \eta_{\epsilon}(N),\] at finite size system \(N\). Here \(\langle\hat{O}\rangle^{\mathrm{ens}}\) presents the ensemble average and \(\eta_{\epsilon}(N)\) is a function that vanishes at the thermodynamic limit (\(N\rightarrow \infty\)). We use a specific representation of the TPQ state termed as the canonical TPQ state [2], \[|\beta, N\rangle \equiv e^{-\frac{\beta}{2} \hat{H}}\left|\psi_{R}\right\rangle,\label{eq:cTPQ}\tag{4}\] where \(\beta\) denotes the inverse temperature, and the state \(|\psi_{R}\rangle\) is random state. Calculating the expectation value of the local operator \(\langle \hat{\mathcal{O}} \rangle^{\mathrm{TPQ}}_{\beta,N} \mathrel{\vcenter{:}}= \langle \beta , N | \hat{\mathcal{O}} | \beta, N \rangle /\langle \beta,N |\beta,N \rangle\), then taking an average for initial random states, \(\overline{\langle \hat{\mathcal{O}} \rangle_{\beta, N}^{\mathrm{TPQ}} }\), approaches the thermal ensemble average in the thermodynamic limit, \[\left(\overline{\langle\hat{\mathcal{O}}\rangle_{\beta, N}^{\mathrm{TPQ}}}-\langle\hat{\mathcal{O}}\rangle_{\beta, N}^{\mathrm{ens}}\right) \xrightarrow{N\rightarrow\infty} 0. \label{eq:prop-TPQ}\tag{5}\]

To prepare the canonical TPQ state, Eq. 4 , on a quantum circuit, we first prepare the random state \(|\psi_{R}\rangle\). Here, we use an algorithm proposed in Ref. [11], which incorporates a staggered arrangement of local unitary operations.

The problem is how to perform the imaginary time evolution \(e^{-\frac{\beta}{2} \hat{H}}\). Here, we utilize the QITE algorithm [7], which is a quantum-classical hybrid algorithm designed to perform the imaginary time evolution on a quantum circuit. In an infinitesimal time step (\(\Delta \beta\)) of the Trotter decomposition of a given Hamiltonian \(\hat{H} = \sum_{l=1}^{L} \hat{h}[l]\), the imaginary time evolution, which is expressed by a non-unitary operator, can be approximated by a unitary operator, \(e^{-i \Delta \beta \hat{A}[l]}\). The unitary operator can be expanded by the following set of unitary operators, \[e^{-i \Delta \beta \hat{A}[l]}= e^{-i \Delta \beta \sum_{I} a_{I}[l] \hat{\sigma}_{I}}.\label{eq:qite95unitary}\tag{6}\] Here, \(\hat{\sigma}_{I} \in\{I, X, Y, Z\}^{\otimes D}\) denotes the Pauli string, which acts on \(D\)-qubit domain. Since the index \(I\) runs from \(1\) to \(4^D\), the coefficients \(\{ a_{I}[l] \}\) form a real vector of size \(4^D\). On the other hand, we suppose that a part of the Hamiltonian in each Trotter step, \(\hat{h}[l]\), is a \(k\)-local operator. The domain size \(D\) should be larger than \(k\) for each Trotter step. In the case of a local Hamiltonian system, the algorithm works well even in \(D\ll N\) [7]. In our Hamiltonian 1 , it includes all-to-all qubit interactions, so that we take \(D=k(=N)\) in our calculation.

Using the unitary operators, we obtain a state \(|\Phi\rangle_{\text{unitary}} = e^{-i \Delta \beta \sum_{I} a_{I}[l] \hat{\sigma}_{I}}|\psi\rangle\). On the other hand, a desired state is \(|\Phi\rangle_{\text{non-unitary}} = C e^{-\Delta \beta \hat{h}[l]}|\psi\rangle\), where \(C\) is a normalization factor \(C=\langle \psi |e^{-2\Delta\beta \hat{h}[l]}| \psi\rangle\). Here, we initialize the state \(|\psi\rangle\) with a random state for our purpose of preparing the canonical TPQ state. To find an optimized value of \(a_I[l]\), we minimize the norm difference between these two states, \(\| \;|\Phi\rangle_{\text{unitary}}-|\Phi\rangle_{\text{non-unitary}} \;\|,\) using classical computers.

In the process of fixing the value of \(a_I[l]\), we calculate (\(4^D \times 4^D\)) matrix elements containing the expectation values of the Pauli strings. Ideally, the calculation of the expectation values is performed by quantum computers in the QITE algorithm. The bottleneck of this algorithm is the calculation of these numerous expectation values, in particular, for our non-local Hamiltonian case which requires a large value of \(D\). To address this computational obstacle, we developed a method to handle multiple Hamiltonian terms with different localities in a single Trotter step. Furthermore, by utilizing the symmetry of Pauli strings, we reduce the computational cost for both quantum and classical calculations. These improvements allow us to compute with large system sizes.

To obtain the thermodynamic quantity from the calculation data given by the QITE algorithm, we have to take double limits; namely the \(\Delta\beta \rightarrow 0\) limit at finite \(N\) to remove the Trotter errors and then \(N\rightarrow \infty\). After taking the first extrapolation, we can see that the result of QITE will be consistent with that of TPQ by classical calculation. After the second extrapolation, we eventually obtain the value of thermodynamic quantity.

4 Simulation results↩︎

4.1 Calculation methods and simulation parameters↩︎

Now, we calculate the chiral condensate, Eq. 2 , as a function of temperature using three independent methods: the exact diagonalization method, the “TPQ method”, and the “QITE method”. In the second method, we classically prepare the TPQ state exactly as defined in Eq. 4 and calculate the expectation value of the chiral condensate. Also in the third method, we generate the TPQ state, but here we perform the imaginary time evolution using the QITE algorithm. This is a quantum-classical hybrid algorithm, but here we simulate the quantum part classically using a state vector simulator.

The lattice spacing and the coupling constant are set to \(a=0.80\) and \(g=1.00\), respectively. The lattice size \(N\) is \(4\)–\(12\). The mass parameters in this work are massless and \(m/g=0.15\). In the massive case, we also study the non-zero \(\theta\) regime. In both the TPQ and QITE methods, we utilize the same set of \(100\) initial random states and take the average over them. Furthermore, in the QITE method, we take the Trotter step in the range of \(\Delta\beta =0.05\)–\(0.25\).

4.2 Results of the massless case↩︎

First, we study the massless case as a feasibility test of our calculation strategy. In this case, the analytical result is known as shown in Eq. 3 . Figure 1 depicts the comparison between the results of the TPQ method (blue circle symbol) and the exact diagonalization results (red curve) at finite \(N\).

The results of the TPQ method do not agree with the exact diagonalization results for each \(N\), but the discrepancy becomes smaller in larger \(N\). We extrapolate these data at each \(\beta\) at \(N=4,6,8, 10\), and \(12\) with the quadratic function and estimate the value at \(N\rightarrow\infty\).

Figure 2 illustrates the \(\beta\)-dependence of the chiral condensate at the thermodynamic limit. Here, we plot the analytical result obtained by the Lagrangian formalism (black solid curve), the extrapolated value of the exact diagonalization results (red triangle symbol), and the extrapolated value of TPQ method (blue circle symbol). The statistical error of the extrapolated values comes from the fitting error of the extrapolation. These three results agree with each other within the statistical error as expected from Eq. 5 .

Furthermore, we can see that at \(\beta=10\) the result converges to the value at zero temperature \(\langle\bar{\psi} \psi\rangle \simeq-0.1599\) [10]. This indicates that \(\beta = 10\) corresponds to a sufficiently low-temperature.

Next, we examine the QITE method. Here, we fix the lattice size as \(N=6\). We take \(\Delta\beta= 0.05, 0.10, 0.125, 0.20\), and \(0.25\) and investigate the effects of the Trotter error. In Figure 3, we plot the absolute difference between results by the QITE and TPQ methods as a function of the Trotter step size (\(\Delta \beta\)). It can be fitted well using the quadratic function of \(\Delta \beta\) and the extrapolated value at \(\Delta\beta\rightarrow 0\) is consistent with zero within the statistical error. This result indicates that the Trotter error is under control and the QITE method reproduces the results of the TPQ method. Combined with the consistency between the results of the TPQ method and the analytical results in the thermodynamic limit, our strategy based on TPQ and QITE looks promising to study the local observables at finite-temperature even for the Schwinger model Hamiltonian, which includes non-local terms.

4.3 Results of the massive case↩︎

We now consider the massive and finite \(\theta\) regimes. Figure 4 shows the value of the chiral condensate computed by the QITE method in the \(\beta\) – \(\theta\) plane. Here, we take \(m/g=0.15\) and \(0 \leq \theta \leq \pi\) with the \(\theta/2\pi =0.1\) interval.

We can see that the absolute value of the chiral condensate is large only in small \(\theta\) and large \(\beta\) regimes. Outside this region, it gradually approaches zero, indicating that the (discrete) chiral symmetry is being restored. The next task is the double extrapolation of \(\Delta \beta \rightarrow 0\) and \(N \rightarrow \infty\), which is an ongoing work.

5 Summary and discussion↩︎

In this study, we simulated the Schwinger model at finite temperatures using the TPQ state approach and the QITE algorithm. As a feasibility test, we first considered the massless Schwinger model and showed that the (classical) TPQ method reproduces the result given by the exact diagonalization in the thermodynamic limit. Furthermore, these results are consistent with the analytical result in the continuum limit based on the Lagrangian formalism. We next confirmed that the results of the QITE method converge to those of the (classical) TPQ method in the Trotter error \(\Delta \beta \rightarrow 0\) limit. Thus, we found that the Trotter error for the discrete imaginary time evolution is under control in our calculation setup. We conclude that our method based on TPQ and QITE works well even for the Schwinger model Hamiltonian, which includes non-local terms.

We also simulated the massive Schwinger model with non-zero \(\theta\)-term. We explored the \(\theta\) and \(\beta\)-dependence of the chiral condensate for a finite system. The chiral condensate exhibits non-zero values in the region of small \(\theta\) and high \(\beta\), and gradually approaches zero outside this region. Taking both \(\Delta \beta \rightarrow 0\) and \(N\rightarrow \infty\) limits in the massive and non-zero \(\theta\) case will be reported in the near future.

We would like to thank K. Nitadori for the useful discussions. The numerical simulations are supported by the Yukawa-21 in YITP, Kyoto University and HOKUSAI supercomputer at RIKEN (Project ID No. Q22577). The work of J. W. P is supported in part by the JSPS Grant-in-Aid for Research Fellow Number22J14732 and the JST SPRING, Grant Number JPMJSP2108. The work of E. I. and S. Y. is supported by JST Grant Number JPMJPF2221 and Program for Promoting Researches on the Supercomputer “Fugaku”, Joint Institute for Computational Fundamental Science (JICFuS), Grant Number JPMXP1020230411.The work of E. I. is supported by JSPS KAKENHI with Grant Number 23H05439, JST PRESTO Grant Number JPMJPR2113, and JSPS Grant-in-Aid for Transformative Research Areas (A) JP21H05190. The work of S. Y. is supported by JSPS Grant-in-Aid for Scientific Research (A) JP21H03455, and the COE research grant in computational science from Hyogo Prefecture and Kobe City through Foundation for Computational Science.

References↩︎