Introduction
In the field of many-body calculations, phase transitions are typically characterized by identifying an appropriate order parameter and determining their locations through the changes or scaling behavior of that parameter. [1]. In recent years, with the continuous development of the field of quantum information, the concept of quantum entanglement has also been introduced into the characterization of many-body problems [2, 3]. Entanglement describes non-local quantum correlations within a system, and such behavior has been found to be powerful for distinguishing different phases and phase transitions [4, 5]. As a result, the study of entanglement provides a unique perspective and has become a widely concerned physical observable in current research. Here, we focus on the simplest and most studied case—bipartite entanglement. Among various bipartite entanglement observables, the entanglement spectrum (ES) can encompass more physical information, since it serves as a spectral structure [6]. Several previous studies have used the ES to investigate problems such as topological phases with bulk-edge correspondence [7–14] in different systems: spin lattices [15–23] and fermionic systems [24–31]. These studies further offer a deeper understanding of the inherent symmetries of the model and the emergent symmetry at phase transition points [32, 33]. It indicates that the property of the ES is an extremely important observable for exploring phase of matter.
However, the study of entanglement requires partial trace of the whole density matrix, named the reduced density matrix (RDM) [34, 35], which is not an easy task in actual many-body calculations due to the huge degree of freedom. Specifically, the existing methods for studying the ES are mainly based on the exact diagonalization method (ED), density matrix renormalization group (DMRG), and its two-dimensional generation—projected entangled-pair states (PEPS) [36, 37], or few analytical derivation methods. All these approaches become particularly difficult to extend to dimensions larger than one, due to the inherent algorithmic structures. Even if applied to the higher dimension, the calculations on large system sizes are limited, because of the complexity in the full trace operation of the environment. Thus, the information uncovered by the ES suffers from significant finite-size effects, especially when dealing with long-range problems. As a result, research on the ES is mainly confined to one-dimensional systems, and the development of methods for studying the ES in high-dimensional systems, as well as the exploration of its physical mechanisms, is still ongoing. Fortunately, the quantum Monte Carlo algorithm (QMC) is an unbiased method and has advantages in solving high-dimensional problems. In recent years, there have been several advancements in QMC algorithms for studying bipartite entanglement with high efficiency, such as entanglement entropy and Rényi negativity [33, 38–45]. The advantage of these newly proposed algorithms is that they are in principle able to handle all two- or higher-dimensional interacting systems between spins or fermions by assigning a segment calculation. Hence, the blanks in the ES calculation by QMC need to be filled in.
In this work, we use the determinant quantum Monte Carlo (DQMC) method to study the low-lying levels of the ES in the interacting fermion models as the first time attempt. To optimize the computational method, we analyze and identify all possible conserved quantities to simplify the computation of the ES. As an example, we calculate the ES of the two-dimensional Su–Schrieffer–Heeger (SSH) model under the many-body Hubbard interaction, and identify the associated topological phase transition, where the conclusions are consistent with that from the ED methods in the one-dimensional system. Remarkably, we also apply the algorithm to the case of nearest-neighbor density interaction, and find a different type of topological phase transition due to the formation of the long-range order. In addition, the feature of the ES and the wavefunction analysis indicates the existence of the long-range interaction of the entangled Hamiltonian (EH).
Calculation of the entanglement spectrum
We provide details on how to treat the RDM as a feasible observable. Pioneered by [46], within the DQMC regime ρM is expressed as the averaged RDM among configurations,
$\begin{eqnarray}{\rho }_{{\rm{M}}}=\displaystyle \sum _{s}{P}_{s}{\rho }_{{\rm{M}},s},\end{eqnarray}$
where s labels the auxiliary field, Ps is the probability for each auxiliary field. Since the single auxiliary field configuration serves as the non-interacting fermionic system, ρs now directly depends on the Green function matrix Gs, $\begin{eqnarray}{\rho }_{{\rm{M}},s}={C}_{s,{\rm{M}}}{{\rm{e}}}^{-{{\boldsymbol{c}}}^{\dagger }\mathrm{log}\left({G}_{s,{\rm{M}}}^{-1}-{\mathbb{1}}\right){\boldsymbol{c}}}={C}_{s,{\rm{M}}}{{\rm{e}}}^{-{{\boldsymbol{c}}}^{\dagger }{H}_{s,{\rm{M}}}{\boldsymbol{c}}},\end{eqnarray}$
c = (c1, c2, ⋯ ) denotes the fermionic vector, and ${G}_{s,{\rm{M}}}=\langle {c}_{i}^{\dagger }{c}_{j}\rangle ,\quad i,j\in {\rm{M}}$ is the Green function matrix defined on M. Hs,M is known as the reduced Hamiltonian or the EH for the configuration s [34]. ${C}_{s,{\rm{M}}}=\,\rm{Det}\,({\mathbb{1}}-{G}_{s,{\rm{M}}})$ normalizes the trace of RDM to unity. In the QMC simulation, ρM,s serves as an observable, which should be measured as many times as possible to control the error-bar. However, the computational cost for measuring such observables is heavy, resulting in the inability to the large entangled size. To be specific, the most costly part, which comes from the exponential computation with diagonalization of the matrix ${d}^{{N}_{{\rm{M}}}}\times {d}^{{N}_{{\rm{M}}}}$, is performed on each sample, where d is the degree of freedom defined on the single site, NM is number of site in M.To deal with the high complexity of the exponential calculation, we begin with the discussion on the properties of Gs,M to search for the simplified approach. In the DQMC algorithm, the Green function matrix of each auxiliary field is determined by the configuration of the auxiliary field. However, the single auxiliary field could be arbitrary, and its impact on Gs,M is similar to applying a certain disorder to the matrix element. Therefore, Gs,M and associated Hs,M are expected to be a totally featureless matrix. Nevertheless, a special case for simplification is that Gs,M is already block-diagonalized. The corresponding Hamiltonian includes the model with a spin degree of freedom or bilayer structure without inter-layer hopping terms. At this time, there are no non-diagonal matrix elements between two sectors at specific configurations.
From the perspective of equation (2 ), we notice that Hs,M serves as a generic matrix with basis c and c†, which has the same form with the non-interacting fermionic Hamiltonian with all-to-all hopping amplitude. We find the only existing symmetry, the particle number conservation, which survives owing to its commutation relation with the hopping term. From another perspective, since the Hubbard–Stratonovich transformation essentially reduces the higher symmetry of the original Hamiltonian to U(1), we could only make use of U(1) symmetry as the subgroup of higher symmetry, which corresponds to the particle number conservation. Considering the general EH matrix with spinless fermions, its full spectrum is block-diagonalized to NM + 1 sectors with respect to the conserved total particle number N, where NM is the size of the entangled region. Each square blocked RDM has the dimension expressed as the combinatorial number ${C}_{{N}_{{\rm{M}}}}^{n}$.
In this work, the particle number conservation and the block property with spin index are both considered to simplify the calculation. In the simulation, we divide the RDM into block sectors, then calculate exponentiation and store separately. Specifically, the total Hilbert space of the entangled region is now ${2}^{{N}_{{\rm{s}}}\times {N}_{{\rm{M}}}}$, where Ns is the number of spin flavors. We separate the whole RDM into ${({N}_{{\rm{M}}}+1)}^{{N}_{{\rm{s}}}}$ small blocks labeled {n1, n2, ⋯ , ns} of the particle numbers of each spin flavor, with its dimension ${C}_{{N}_{{\rm{M}}}}^{{n}_{1}}\times {C}_{{N}_{{\rm{M}}}}^{{n}_{2}}\times \cdots \times {C}_{{N}_{{\rm{M}}}}^{{n}_{{\rm{s}}}}$. In the final step, we do the diagonalization for each block to obtain its eigenvalue spectrum and eigenstates of each particle number N and spin index σ. These simplifications guarantee enough of the system size of the entangled region to approach the nature in the thermodynamics limits.
Based on this, we provide a rough estimation of the algorithm complexity to compare with the ED algorithm. The most costly process among both algorithms is the diagonalization process to obtain eigenvalues. However, we emphasize that above blocked diagonalization properties also exist in the original ED methods. The dominating simplification is the combination of the Monte Carlo algorithm. Reminiscent of the ED algorithm, the complexity is equivalent to the exponential function of the total Hilbert space, O(${d}^{{N}_{{\rm{L}}}}$), where NL is the total system size. Thus, the complexity of Monte Carlo algorithm is the exponential function of the entangled region size, O(${d}^{{N}_{{\rm{M}}}}$). More accurately, in addition to the Monte Carlo sampling process, the ES calculation by DQMC should be O(${N}_{{\rm{L}}}^{3}{d}^{{N}_{{\rm{M}}}}$), which reduces the complexity from the exponential order O(${d}^{{N}_{{\rm{L}}}}$) to the algebraic order of the total system size. When NM ≪ NL, the algorithm could demonstrate its advantages. Therefore, it is better to solve the ES with the Monte Carlo approach, where the system has a relative small entangled region and a large degree of freedom from the environment.
At the end of this section, we make a simple comparison between the current methods using QMC and several other commonly used algorithms, such as ED, DMRG and PEPS. For the ED, its biggest advantage is the accuracy of the full spectrum, while it is limited in one-dimension or a few lattice sites. For the DMRG and PEPS, it is natural to obtain the ES due to its tensor representation, thus offering the most convenient approach to extract the entanglement information, only if the ground state wavefunction is constructed. However, the essential issue is finding the appropriate Hamiltonian, where the ground state representation is of a lower bond dimension. It is known that the DMRG is enabled to handle most one-dimensional problems, but for two-dimensions, only certain Hamiltonians are suitable to employ PEPS to study the ES, for example the AKLT model, toric code model [47]. At this point for the QMC algorithm, it could be well applied to most sign-problem-free fermionic models in two or higher-dimensions with various lattice type. Though the computational complexity of the QMC algorithm is higher than the DMRG/PEPS in case of the low bond dimension, the sampling algorithm still demonstrates superiority when the system size is large (scale with O(${N}_{{\rm{L}}}^{3}$). Therefore, the QMC algorithm provides an optional approach for exploring entanglement properties in two-dimensions.
Two-dimensional SSH model
We apply the algorithm to the entanglement study on the SSH model [48]. The one-dimensional SSH model exhibits the topological phase transition from a trivial gapped state to the topological phase with zero-energy modes by tuning the ratio of two staggered coupling constants. Here, we extend the one-dimensional SSH model to two-dimensions, which links several one-dimensional SSH chains with staggered hopping strength along the vertical direction, i.e., SSH-type chains along y-axis. On the one hand, the computation on this two-dimensional model enables us to reveal the superiority of our method compared to the ED and DMRG methods. On the other hand, a total of four edge states contribute to the ground state degeneracy, bringing more complexity. The non-interacting Hamiltonian H0 writes as
$\begin{eqnarray}\begin{array}{rcl}{H}_{0} & = & {H}_{x}+{H}_{y,}\\ {H}_{x} & = & -{t}_{xa}\displaystyle \sum _{j,i=\,\rm{even}\,}^{{L}_{x}}{c}_{(i,j)}^{\dagger }{c}_{(i+1,j)}-{t}_{xb}\displaystyle \sum _{j,i=\,\rm{odd}\,}^{{L}_{x}}{c}_{(i,j)}^{\dagger }{c}_{(i+1,j)}+h.c.,\\ {H}_{y} & = & -{t}_{ya}\displaystyle \sum _{i,j=\,\rm{even}\,}^{{L}_{y}}{c}_{(i,j)}^{\dagger }{c}_{(i,j+1)}-{t}_{yb}\displaystyle \sum _{i,j=\,\rm{odd}\,}^{{L}_{y}}{c}_{(i,j)}^{\dagger }{c}_{(i,j+1)}+h.c..\end{array}\end{eqnarray}$
The sketch map of the two-dimensional SSH model is shown in Figure 1. The alternating hopping strengths are along both x and y directions, denoted as txa, txb, tya, tyb, respectively. Here, under the condition of open boundary conditions (OBC) in the y-direction, the entangled region consists of two boundaries. Here we show the reason of such a model design in the following.
Figure 1. Sketch map of the two-dimensional SSH model. (a) The thick and thin bonds represent the weak hopping strength txa, tya and strong hopping strength txb, tyb, respectively. We choose periodic boundary conditions (PBC) along the x-axis and the OBC along y-axis. The aim of the OBC in the y-direction is to separate the two boundaries, and mimic the physics of the square entangled region. The entangled region M consists of two parts on the boundaries and is colored by red. In particular, for the Hubbard interaction, the entanglement region M has Mx = 4, and for the density interaction, Mx = 6. Note that the four lattice sites are the smallest size along the x direction to guarantee the presence of the zero-energy modes by cutting the region. The total system size is fixed with Lx = 16, Ly = 8. Since the finite-size effect comes from both the entangled region and the total system size, we numerically find that the total system size here is large enough to avoid finite size effect generated by the total system. (b) The lattice sketch map in the three-dimensional perspective. The two boundaries are not directly connected in the region M. |
Thereafter, considering the many-body interactions, we choose two types of them: (1) onsite Hubbard interaction HU; (2) spinless nearest-neighbor interaction HV:
$\begin{eqnarray}\begin{array}{rcl}{H}_{{\rm{U}}} & = & U\displaystyle \sum _{I}\left({n}_{I,\uparrow }-\frac{1}{2}\right)\left({n}_{I,\downarrow }-\frac{1}{2}\right),\\ {H}_{{\rm{V}}} & = & V\displaystyle \sum _{\langle IJ\rangle }\left({n}_{I}-\frac{1}{2}\right)\left({n}_{J}-\frac{1}{2}\right).\end{array}\end{eqnarray}$
To avoid the sign problem, we focus only on the half-filled case. The inverse temperature is β = 100, which is low enough to reach the ground state. The two-dimensional square lattice site is represented by the capital Roman I, J, with its coordinates denoted as (i, j).Moreover, we mention some algorithm details on the simplification of the ES calculation. In fact, the block-diagonalized property is used in the transformation process from Hs to ρM,s. The basis vectors defined on the Fock space are firstly appointed, thus we are able to classify these in different particle number sectors. Based on this, we allocate the memory for each blocked matrix instead of the total density matrix. Then we could simply keep the blocked matrix sequence to calculate its average and do diagonalization for the ES.
Entanglement spectrum under the Hubbard interaction
In this section, we perform the ES algorithm by adding the Hubbard interaction accompanied by the spin degree of freedom. For parameter setting, we fix txa = tya = 0.2, txb =tyb = 0.8, and the edge cutting of the entangled region M is also located on the weak hopping bonds [Figure 1(a)]. The ES results and the ground state degeneracy with respect to the particle number of spin up and down are shown in Figure 2.
Figure 2. The ES versus the total particle number at the Hubbard case. From left to right panels, U = 0, 2, −2, respectively. (a) The degenerate ground states are circled by gray lines. Due to the finite size of the entangled region, the energy levels have some split with each other, and will become completely degenerate at infinite NM. In any case, one could still distinguish the ground states due to the their large gaps with the first excited state. At the presence of the Hubbard interaction, the contributions from the edge state and bulk state are framed by the blue lines, respectively. (b) The ground state degeneracy distribution with respect to the particle number of the two spin indexes. The topological phase transition occurs by adding the Hubbard interaction, manifesting ground state degeneracy changes from 256 to 16, but with different distribution for positive and negative U. |
At the very beginning, we focus on the physical mechanism between the one-dimensional and two-dimensional models. Since the two-dimensional SSH model owns staggered hopping strength at both the x and y-direction, we could also extend the entangled region from a bar in one-dimension to a square shape in two-dimensions, so as to take advantage of the topological properties at both directions. However, due to the high computational complexity, we could only reach a small entangled region as the above analysis. In the following text, one could see such a choice also captures the two-dimensional physics. Nonetheless, in such a case, the number of edges changes from two sites to four sites distributed on the corner of the square, which is illustrated in the change of the ground state degeneracy. As for the interaction term, at the topological phase, the Hubbard-type interaction acts on the edge modes, which further alters the distribution of the ES. This is tenable in both one-dimension and two-dimensions, and the only difference between them is the number of edge states. Based on the above analysis, we expect that the conclusion of the one-dimensional model in [49] could be generalized to the two-dimensional system. In the following text, we would display the numerical simulation to verify the deduction.
We first discuss the ground state degeneracy for benchmarking the code. Notice that the subsystem has a total of eight zero-energy modes (two edge states along both the x and y directions, two spin’s degree of freedoms). Therefore, the associated ES of the non-interacting case has 256-fold (28) ground state degeneracy [figures 2(a1) and (b1)]. When adding HU, for a one-dimensional SSH model, the ground state degeneracy changes from 16-fold (two edge modes and two spin flavors) to 4-fold [49]. Applying to the two-dimensions, the ground state degeneracy changes from 162 to 42. This is verified in both figures 2(b2) and (b3).
For the particle number distributions, we find the positive and negative U cases are different. For example, considering negative U, the pair of fermions with spin up and spin down tends to either locate both or not locate on the same site to minimize the energy. Therefore, it distributes on the diagonal line, as shown in Figure 2(b3), which is in contrast to the positive U case in (b2) with single particles occupied on each edge state.
Next, we focus on the characters of the whole spectrum. We find that the whole spectrum consists of distinguishable information of the edge modes and the bulk, denoted by the blue circles in figures 2(a2) and (a3). Especially at positive U, the contribution of the edge modes own discrete and equal energy gap distribution. It is equivalent to the energy difference for the gain or loss of a particle. For the local Hubbard interaction, the energy gap is proportional to a constant (U), consistent with the results in Figure 2(a2). Moreover, since each edge mode is mutually independent, the joint effect of the four lattice sites on the boundary of M are simply expressed as four replicas. Therefore, the edge state with the maximum energy has four more particles and an equal energy increase than the ground state, as also verified by Figure 2(a2). Except for the discrete portion of the ES, the continuous portion comes from the bulk state, which also corresponds to the higher-energy levels.
The above analysis is also appropriate for the negative U case [Figure 2(a3)]. The edge spectrum is the same as the positive U case, if one transfers the horizontal axis from the total particle number N to the total fermionic spin ns = n↑ − n↓. Unfortunately, the bulk spectrum covers a large proportion of the edge spectrum, so that we could only clearly identify the degenerate ground states.
However, on account of the finite size effect Mx = 4 of the entangled region, the eigenvalues are not perfectly degenerate as in the thermodynamic limit, even if in the non-interacting case (the gray circle in Figure 2(a1)). We could only distinguish the ground state degeneracy by both numerical results and theoretical analysis. In addition, the continuous feature of the high-energy part originates from two reasons: (1) physically, the quasi-particles have relatively strong interactions with each other as a typical feature at high-energy, leading to the breakdown of the one-body approximation. (2) Numerically, the eigenvalues of the high-energy part of the EH are exponentially small in the raw data of the ES, which are greatly influenced by the machine error and sampling error. Therefore, for the ES calculated by the QMC simulation, only the low-energy part could be credible and provide valuable content.
At the end of this section, we would like to emphasize that for a normal QMC algorithm by directly simulating the system with an open boundary SSH model, it is hard to recognize the ground state degeneracy, and hence distinguish the topological phase transition. In any case, the investigation of the ES could provide insight into the state properties of the total system. And, most importantly, our DQMC computation is beyond the capacity of the ED method, which has been used to study the one-dimensional SSH model [49].
Entanglement spectrum under the density interaction
We proceed to employ the spinless density interaction HV to the system, and explore its topological properties. In particular for the QMC algorithm, the Hubbard–Stratonovich transformation decomposes the fermionic quadratic form to the linear form with the auxiliary field s, namely, ${{\rm{e}}}^{-\beta {H}_{{\rm{V}}}}\to {\int }_{s}{C}_{s}{{\rm{e}}}^{\sqrt{V}s({n}_{i}-{n}_{j})}$, where β is the inversed temperature, and Cs is the coefficient dependent on the s. Clearly, the attractive interaction leads to the imaginary auxiliary field value, which further renders the complex weight in QMC. Therefore, to avoid the sign problem for such an interaction, we only focus on the V > 0 region in the following QMC simulation. Since we consider the spinless system, the non-interacting ground state degeneracy is now 16-fold and stretches across various particle number sectors [Figure 3(a1)]. For comparison, we also calculate the uniform hopping case, i.e., txa = txb = tya = tyb = 0.5, where only a single ground state exists [Figure 3(c1)].
Figure 3. The ES versus the total particle number at the existence of HV. (a,b) for the SSH model and (c,d) for uniform hopping. V = 0, 0.5, 1 from the left to the right panels, respectively. The ground states are circled by gray lines. The lowest energies of the bulk contribution fit well with the quadratic function as the function of particle number N for both cases at V = 1. The topological phase transition of SSH model occurs by adding the HV interaction, manifesting ground state degeneracy changes from 16 to 2. |
Firstly, we explore the variation of the ground state degeneracy. Unlike the Hubbard-type interaction defined on the single site, the nearest-neighbor interaction could lead to the phase transition of the whole system. In particular, the repulsive density interaction promotes the charge density wave (CDW) order for a spinless fermion. Such an order introduces the correlation between edge states, and only the occupation numbers that satisfy the CDW order remain on the ground state. For SSH-type hopping in the one-dimensional case, the degeneracy is expected to transfer from 4 to 2, i.e. (01) and (10) for particle numbers of two edge states. The conclusion also applies to the two-dimensional system, once we turn on the interaction [Figure 3(a)]. In contrast, the uniform case has a single ground state at V = 0. We observe that as V increases, the gap of the first excited state rapidly decreases. Although the gap is small enough, the degenerated ground states could only exist in the presence of SSH-type hopping.
Going deep into the topological phase transition, we analyze the degenerated ground state wavefunction of the EH in Figure 4. The average particle density at each site ni is found to be 0.5 and the average value of the site-resolved nearest-neighbor density correlation dIJ = nI(1 − nJ) is close to 0.5 for both degenerate states. Even the density correlations between the disconnect two edge modes ${d}_{(i,1)(i,{L}_{y})}$ are also approximate to 0.5, revealing a global CDW order along both directions. Therefore, the two degenerate states are the linear combination of ∣Φ1⟩ = ∣0, 1, 0, 1, ⋯ ⟩ and ∣Φ2⟩ = ∣1, 0, 1, 0, ⋯ ⟩ in the occupying number representation, demonstrating the 2-fold degeneracy of the ES.
Figure 4. Site-resolved CDW correlation function for entangled region of SSH case. The thin and thick bonds represent the weak and strong hopping strength, respectively. All average values are closed to 0.5, dated from the strong global CDW order. Furthermore, we find d(1,0)(2,0) and d(3,0)(4,0) is larger than d(0,0)(1,0), d(2,0)(3,0) and d(4,0)(5,0), consistent with the alternating hopping strength of original Hamiltonian. Note that both the two entangled wavefunctions are indistinguishable from the studied observables. |
Next, we still focus on the features of the whole spectrum. We notice that at large V = 1, the spectrum of the two cases becomes similar, where figures 3(a3) and (c3) show the large excitation gaps and quadratic form with respect to the particle number of the continuous spectrum. That indicates that all states in the continuous spectrum have a small overlap with ∣Φ1⟩ and ∣Φ2⟩, and the gap is related to V. Moreover, from the perspective of the quasiparticle, the quadratic form can be interpreted as the all-to-all interaction between quasi-particles to the lowest order. It means placing an extra one particle or hole leads to the energy increase proportional to the number of existing quasi-particles. The all-to-all interaction has a close relation with the symmetry breaking and the formation long-range order of the whole system.
In summary, we identify the topological phase transition of the two-dimensional SSH model from 16-fold to 2-fold by adding the density interaction and analyzing its ES characters. The most interesting thing is that the phase transition in the ground state of the EH also reflects the transition in the ground state of the original Hamiltonian. In particular, the variation of the EH degeneracy is related to the type of the interaction, which is further connected with the different topological phase transitions of the total system. This could answer a long-standing controversy in the field [23, 50].
Summary and outlook
We report an efficient way to numerically calculate the ES and its entangled wavefunction of the two-dimensional interacting system based on large-scale QMC algorithms. We analyze its algorithm optimization, and leverage the particle number and spin conservation to simplify the computation. First, we verify the algorithm with the two-dimensional SSH model with the Hubbard interaction. As an application, we investigate the topological phase transition of the SSH model under the density-density interaction, by analyzing both the structure of the ES and the eigenfunctions. Importantly, our scheme could be extended to two- or higher-dimensions with a large system size, which surpasses the computation capacity of other exact or numerical methods. However, the current available size of entangled region is limited to a small subsystem. To reach the larger entangled region, one could use some approximation technique to extract the ground state or low-excited state information of the small blocked matrix, for example the Lanzcos algorithm. Nevertheless, we show the ability of the QMC algorithm to study the ES, which gives researchers more avenues in numerical techniques. Our method thus provides an opportunity to explore the physical mechanism of interacting fermion systems via entanglement properties.