Block product density matrix embedding theory for strongly correlated spin systems

TL;DR

This paper enhances block product DMET with spin-state optimization to better study strongly correlated spin systems, applying it to models like the Kitaev-Heisenberg and analyzing energy and correlation functions.

Contribution

It introduces an improved variational Ansatz for BPDMET using spin-state optimization and applies it to complex spin models, providing insights into excited states and spectral functions.

Findings

Improved accuracy in energy calculations for spin models

Effective analysis of correlation functions

Insights into excited states via tangent space diagonalization

Abstract

Density matrix embedding theory (DMET) is a relatively new technique for the calculation of strongly correlated systems. Recently, block product DMET (BPDMET) was introduced for the study of spin systems such as the antiferromagnetic $J_1 - J_2$ model on the square lattice. In this paper, we extend the variational Ansatz of BPDMET using spin-state optimization, yielding improved results. We apply the same techniques to the Kitaev-Heisenberg model on the honeycomb lattice, comparing the results when using several types of clusters. Energy profiles and correlation functions are investigated. A diagonalization in the tangent space of the variational approach yields information on the excited states and the corresponding spectral functions.

Figures

Equations

$$|\Psi ⟩=\sum _{i=1}^A\sum _{j=1}^B{\Psi }_{ij}|{\alpha }_i⟩|{\beta }_j⟩=\sum _k^{\min (A,B)}{\lambda }_k|{\tilde{\alpha }}_k⟩|\ast ⟩{\tilde{\beta }}_k.$$

$$\hat{H}=\sum _{m,n}\sum _{\mu ,\nu }{J}_{\mu \nu }^{mn}{\hat{S}}_m^{\mu }{\hat{S}}_n^{\nu }.$$

$$|\Psi ⟩=\sum _i{a}_i|{\alpha }_i⟩|{\mathrm{BPS}}_i⟩,$$

$$|{\mathrm{BPS}}_i⟩=\prod _{C\in \mathrm{bath}\mathrm{clusters}}\sum _{\beta }{b}_{C\beta }^i{|\beta ⟩}_C.$$

$$|\Psi ⟩=\sum _i{a}_i|{\alpha }_i⟩\left(\prod _{C\in \mathrm{bath}cl.}\sum _{\beta }{b}_{C\beta }^i{|\beta ⟩}_C\right),$$

$$\sum _{\beta }{b}_{C\beta }^{i\ast }{b}_{C\beta }^i=1,$$

$$\sum _i{a}_i^{\ast }{a}_i=1$$

$$\begin{array}{rl}|\Psi ⟩=& \sum _{i(\ne i_0)}{a}_i|{\alpha }_i⟩\prod _C\sum _{\beta }{b}_{C\beta }^i{|\beta ⟩}_C+\\ & \sum _{\beta }{a}_{i_0}{b}_{B_0\beta }^{i_0}|{\alpha }_{i_0}⟩{|\beta ⟩}_{B_0}\prod _{C(\ne B_0)}\sum _{{\beta }^{\prime }}{b}_{C{\beta }^{\prime }}^{i_0}{|{\beta }^{\prime }⟩}_C.\end{array}$$

$$|{\varphi }_{{\alpha }_i}⟩$$

$$|{\varphi }_{\beta }⟩$$

$$\mathcal{L}=\text{expectationvalue}\hat{H}\Psi -\lambda \text{innerproduct}\Psi \Psi ,$$

$$|\Psi ⟩=\sum _i{a}_i\sum _m{U}_{im}|m⟩\prod _{C\in \mathrm{bath}cl.}\sum _{\beta }{b}_{C\beta }^i{|\beta ⟩}_C.$$

$$\text{expectationvalue}{\hat{S}}_{\mathrm{tot}}^z=\sum _C\text{expectationvalue}{\hat{S}}_{\mathrm{to}{\mathrm{t}}_{\mathrm{C}}}^z{\Psi }_C=N_C\text{expectationvalue}{\hat{S}}_{\mathrm{to}{\mathrm{t}}_{\mathrm{C}}}^z{\Psi }_C,$$

$$\text{expectationvalue}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j=\frac{1}{2}\text{expectationvalue}\ast {\hat{\mathbf{S}}}_j\cdot {\hat{\mathbf{S}}}_i{\Psi }_{C_i}+\frac{1}{2}\text{expectationvalue}\ast {\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j{\Psi }_{C_j},$$

$$\text{expectationvalue}\hat{H}$$

$$\text{expectationvalue}{\hat{\mathbf{S}}}_{\mathrm{tot}}^2$$

$$\sum _{\beta }{b}_{C\beta }^{i\ast }{b}_{C\beta }^i$$

$$\sum _j{U}_{ij}^{\ast }{U}_{ij}$$

$$\sum _i{a}_i^{\ast }{a}_i$$

$$A(\omega ,\hat{X})=-\frac{1}{\pi }\text{imaginary}[\text{expectationvalue}{\hat{X}}^{†}\frac{1}{\omega -(\hat{H}-E_0)+i\eta }\hat{X}{\varphi }_0],$$

$$A(\omega ,\hat{X})\approx -\frac{1}{\pi }\text{imaginary}\left[\sum _n\frac{|\text{matrixelement}{\varphi }_n\hat{X}{\varphi }_0{|}^2}{\omega -(E_n-E_0)+i\eta }\right]$$

$$(\omega -(\hat{H}-E_0)+i\eta )|{\varphi }_1⟩=\hat{X}|{\varphi }_0⟩.$$

$$A_{i_1{i}_2{i}_3{i}_4\cdots }^m=U_{i_{1}m}{\delta }_{i_1{i}_2}{\delta }_{i_1{i}_3}{\delta }_{i_1{i}_4}\cdots$$

$$\hat{H}=J_1\sum _{⟨i,j⟩}{\hat{S}}_i\cdot {\hat{S}}_j+J_2\sum _{⟨⟨i,k⟩⟩}{\hat{S}}_i\cdot {\hat{S}}_k$$

$$M_N^2(\mathbf{Q})=\frac{1}{N^2}\sum _{ij}⟨{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j⟩e^{i\mathbf{Q}({\mathbf{R}}_i-{\mathbf{R}}_j)},$$

$$\begin{array}{rl}\hat{H}=-J_1& \sum _{x-\mathrm{links}}{\hat{S}}_i^x{\hat{S}}_j^x-J_1\sum _{y-\mathrm{links}}{\hat{S}}_i^y{\hat{S}}_j^y\\ -J_1& \sum _{z-\mathrm{links}}{\hat{S}}_i^z{\hat{S}}_j^z+J_2\sum _{⟨ij⟩}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j.\end{array}$$

$$\begin{array}{r}A(\omega ,\hat{X})=-\frac{1}{\pi d}\text{imaginary}\left[\sum _{i=gs}\sum _{j\ne gs}\frac{|\text{matrixelement}{\varphi }_j\hat{X}{\varphi }_i{|}^2}{\omega -\Delta E_j+i\eta }\right],\end{array}$$

$$|\Psi ⟩=\sqrt{1-{ϵ}^2}|{\Psi }_0⟩+ϵ|{\Psi }_{\mathrm{error}}⟩$$

$$\begin{array}{rl}||\Psi ⟩-|{\Psi }_0⟩|{|}_2& =\sqrt{{\left(\sqrt{1-{ϵ}^2}-1\right)}^2+{ϵ}^2}\\ & =ϵ+\mathcal{O}({ϵ}^2)\end{array}$$

$$\begin{array}{r}\text{expectationvalue}\hat{Q}\Psi -\text{expectationvalue}\hat{Q}{\Psi }_0=\\ \left(\text{matrixelement}{\Psi }_{\mathrm{error}}\hat{Q}{\Psi }_0+\text{matrixelement}{\Psi }_0\hat{Q}{\Psi }_{\mathrm{error}}\right)ϵ\\ +\left(\text{expectationvalue}\hat{Q}{\Psi }_{\mathrm{error}}-\text{expectationvalue}\hat{Q}{\Psi }_0\right){ϵ}^2+\mathcal{O}({ϵ}^3)\end{array}$$

Full text

I Introduction

When studying quantum many-body systems, exactly solving the system becomes unfeasible for large system sizes. The exponential scaling of the Hilbert space dimension with the number of particles inhibits an exact simulation of larger systems, and approximate methods have to be used.

A particular type of quantum many-body systems is the quantum spin-lattice system. For this system, interacting spins are localized on the lattice points of a lattice. In this paper, we extend the block product matrix embedding theory (BPDMET),Fan and Jie (2015) a method that was recently introduced to study quantum spin-lattices. We study the validity of the model by applying it to the $J_{1}-J_{2}$ model on the square lattice with Heisenberg interaction and the Kitaev-Heisenberg model. These are two spin lattice systems of particular interest.

The first system has been a long time subject of research. This partly because of its fundamental interest in its simplicity, but also for its use in Fe-based superconductors and other materials. For example, high-$T_{c}$ superconductivity in iron pnictide (or oxypnictides) has been discovered with LaOFeAs being the first.Si and Abrahams (2008) The Fe atoms form a square lattice in these iron pnictides and they exhibit NN and NNN superexchange interactions that can be described with this $J_{1}-J_{2}$ model (however with $S=1$ or $2$). When the crystal is doped, an effective model with a $t-J_{1}-J_{2}$ Hamiltonian can be suggested, introducing a kinetic component.Si and Abrahams (2008) This should give rise to superconductivity. Also the properties of $\mathrm{Li}_{2}\mathrm{VOSiO}_{4}$ have been investigated through the $J_{1}-J_{2}$ model.Melzi et al. (2001)

The Kitaev-Heisenberg model was first introduced for the theoretical examination of iridium oxides of the form $A_{2}\mathrm{IrO_{3}}$, with $A=\mathrm{Li,\,Na}$.Chaloupka et al. (2010) Herein, $\mathrm{Ir^{4+}}$ ions form honeycomb-like lattice planes and have an effective spin one-half. The interaction between the different effective spins is anisotropic. Experimental evidence has shown that the proposed Kitaev-Heisenberg model is a successful model for the iridium oxides, however some extensions of the model have been introduced for a better description.Kimchi and You (2011); Singh et al. (2012); Chaloupka et al. (2013); Reuther et al. (2014); Rau et al. (2014); Chaloupka and Khaliullin (2015) Particular interest for the Kitaev-Heisenberg model has arisen since the model is able to have a Kitaev spin liquidKitaev (2006) as ground state within a finite parameter region.

To solve these systems, exact diagonalization is unfeasible for larger systems and approximate methods have to be used. Some examples of approximate methods are series expansion,Gelfand et al. (1989); Singh et al. (1999); Kotov et al. (1999); Oitmaa (2015) large-N expansion,Read and Sachdev (1991) density matrix renormalization group,White and Chernyshev (2007); Yan et al. (2011); Jiang et al. (2012) projected entangled pair states,Murg et al. (2009) and coupled cluster methods.Darradi et al. (2008)

An approximate solution of quantum many body systems can also be obtained through embedding theories. Here, the system is divided into two parts: an impurity, cluster or fragment (which is the subsystem of interest) and an environment. Using this division, embedding theories are able to solve the problem approximately.Sun and Chan (2016)

The total Hilbert space of the system is now a direct product of the Hilbert spaces of the impurity and the environment (with dimension $A$ and $B$ respectively). A basis for this Hilbert space is given by $\{\ket{\alpha_{i}}\otimes\ket{\beta_{j}}\}$, where $\ket{\alpha_{i}}$ are states of the impurity and $\ket{\beta_{j}}$ are states restricted to the environment. Every state $\ket{\Psi}$ in the total Hilbert space can be written as

[TABLE]

The latter result is obtained by using a singular value decomposition of the matrix $\Psi_{ij}$. This is known as the Schmidt decomposition of a state. Here, $\ket{\widetilde{\alpha}_{k}}$ are states of the impurity and $\ket*{\widetilde{\beta}_{k}}$ are states of the environment. The summation is restricted to the minimum of the impurity and the environment dimension. Since the dimension of the environment is typically much larger than the dimension of the impurity, the summation is limited by the impurity. It is thus clear that only $A$ states in the environment are needed for the construction of the wave function. If only one of the singular values $\lambda_{k}$ is nonzero, the state $\ket{\Psi}$ can be factorized and impurity and environment are not entangled. However, if several singular values are nonzero, $\ket{\Psi}$ is called entangled.Bennett et al. (1996); Horodecki et al. (2009); Eisert et al. (2010); Eisert (2013)

Embedding theories capitalize on this division of the system. By replacing the environment by an approximate model, one tries to calculate the properties of the impurity accurately and cost-effectively. The simplest option is to approximate the environment in such a way that there is no entanglement with the impurity. This is a good approximation when the Schmidt singular values have one dominant non-zero value. For many systems this is sufficient. However, for systems with strong static correlation between the impurity and environment one has to go beyond the mean-field approximation. Static correlation refers to systems where a product state is not sufficient for a qualitative description of the system, but a superposition of multiple product states is needed, i.e. whenever substantial entanglement is present between the impurity and environment, implying several important Schmidt values $\lambda_{k}$ in Eq. (1).

One of the more powerful and popular embedding theories is dynamical mean-field theory (DMFT)Metzner and Vollhardt (1989); Georges and Krauth (1992); Georges et al. (1996); Zgid and Chan (2011). It maps the system to an impurity and a non-interacting bath in a self-consistent way, using the single-particle Green’s function.

Density Matrix Embedding Theory (DMET) is a new embedding theory first proposed by Knizia and ChanKnizia and Chan (2012) for the Hubbard model. It was later extended to full quantum chemical Hamiltonians.Knizia and Chan (2013) For ground state energies, DMET is a computationally cheaper alternative to DMFT with similar accuracy. The self-consistency for DMET is based on the density matrix, instead of the frequency dependent Green’s function in DMFT. Information about excited states in DMET can still be explored.Booth and Chan (2015); Wouters et al. (2016)

In DMET the entanglement between impurity and environment is explicitly kept and the wave function is of the form given by Eq. (1). Finding the Schmidt basis for the environment $\{\ket*{\widetilde{\beta}_{k}}\}$ can be done if the exact wave function $\ket{\Psi}$ is known. This is, however, not an option since finding the exact wave function is equivalent to solving the many-body problem. DMET solves the lack of a priori knowledge of $\ket{\Psi}$ by embedding the impurity in an approximate bath. Solving this combined impurity and bath system is called the embedded problem. To find this bath space, one can use different techniques. A Fock space of bath orbitals which is obtained from a low-level particle-number conserving mean-field wave function is used in the original Refs. Knizia and Chan, 2012, 2013 and is illustrated extensively in Ref. Wouters et al., 2016. Other methods are also possible: single-particle states from Hartree-Fock-Bogoliubov theoryZheng and Chan (2016); LeBlanc et al. (2015) and antisymmetrized geminal power (AGP) wave functionsTsuchimochi et al. (2015) have also been used. Extensions of DMET to coupled interacting fermion-boson systems through coherent state wave functions for phonons have been described as well.Sandhoefer and Kin-Lic Chan (2016)

Adapting DMET for spin lattices has given rise to the so-called cluster density matrix embedding theory (CDMET), as introduced by Fan et al.Fan and Jie (2015) In this paper, however, we opt for the name block product DMET (BPDMET) in order to avoid possible confusion with fermionic DMET when the impurity consists of a cluster of degrees of freedom. In BPDMET, bath states in a spin lattice system are represented by block product states, which is emphasised by our alternative name. Recently, this method has been further extended by implementing BPDMET with the hierarchical mean-field approach.Qin et al. (2016) In the present work, the original BPDMET is used and the ansatz is further extended with so-called spin-state superpositions in the impurity, yielding improved results.

As will be shown, the BPDMET ansatzFan and Jie (2015) can easily be written as a particular case of a more general tensor network state (TNS). The concept of tensor network states is an increasingly important technique for the description of highly correlated systems. It can be viewed as an extension of the density matrix renormalization group (DMRG) as introduced by White.White (1992, 1993) DMRG was shown to be highly accurate and useful for one-dimensional lattice systems, but also reasonably small two-dimensional lattices can be studied up to high accuracy.Stoudenmire and White (2012) The DMRG algorithm has been later rewritten as an optimization of a matrix product state (MPS) ansatz.Östlund and Rommer (1995); Rommer and Östlund (1997)

The concept of TNS is quite general and includes projected entangled pair states (PEPS)Verstraete and Cirac (2004); Murg et al. (2009) which is the extension of MPS into two and higher dimensions, as well as tree TNS (TTNS)Murg et al. (2010, 2015); Nakatani and Chan (2013); Li et al. (2012a) and complete-graph TNS (CGTNS).Marti et al. (2010); Marti and Reiher (2011)

In section II, the concept of BPDMET is introduced by means of a variational ansatz. The optimization procedure is explained as well as the calculation of the energy and other properties in the BPDMET framework. We also extend the variational ansatz using spin-state superposition in the impurity. Finally, we point out a possible way of extracting information on spectral properties. The link with tensor network states is also made. In section III, the method is applied to the 2D Heisenberg model on a square lattice and to the Kitaev-Heisenberg model on a honeycomb lattice. Results for the energies, correlation functions and the location of quantum phase transitions are studied. The concept of diagonalization in tangent space and the resulting spectral function is applied to both lattice systems. Summary and conclusions are provided in section IV.

II Method

II.1 Block Product DMET

The BPDMET method as introduced by Fan et al.Fan and Jie (2015) can be used for the approximate solution of spin lattice systems. These systems can be used to model magnetic properties of materials. A spin lattice system comprises a number of lattice sites $N$. Each site has a spin degree of freedom interacting with other spins. Only spin-spin interactions are investigated in this paper. For these systems, the total Hamiltonian can be written in its most general form as

[TABLE]

Here, $m$ and $n$ denote lattice indices and $\mu$ and $\nu$ denote the spatial components $x$, $y$ and $z$, i.e. the different measuring directions for the spin operator. Magnetic terms have been excluded but can be introduced in a straightforward manner.

An interesting variational ansatz for the system was provided by Fan et al.Fan and Jie (2015) After division of the spin lattice into an impurity and environment, they propose a replacement of the exact environmental states $\ket*{\widetilde{\beta}_{k}}$ in Eq. (1) by a set of block product states $\ket{\mathrm{BPS}_{k}}$. With this approximation, the wave function of the impurity model becomes:

[TABLE]

where $i$ labels the different states in the Hilbert space of the impurity. To define these block product states, the spin lattice system is divided into different clusters. One of the clusters is the impurity. The other clusters are called the bath clusters. An exemplary division can be seen in Fig. 1. With this division of the lattice system, the block product states are defined as follows:

[TABLE]

Here $\{\ket{\beta}_{C}\}$ is a complete set of states within the Hilbert space restricted to the bath cluster $C$. For example, if each bath cluster contains 3 sites with spin-$\frac{1}{2}$, we have $\{\ket{\downarrow\downarrow\downarrow},\ket{\downarrow\downarrow\uparrow},\ket{\downarrow\uparrow\downarrow},..,\ket{\uparrow\uparrow\uparrow}\}$ in the natural basis (eigenstates of $\hat{S}^{z}$). These block-product states have to be optimized so that Eq. (3) optimally represents the exact ground state wave function. The approximation made is that the correlations within the bath clusters are fully taken into account, while the correlations between the bath clusters are only taken into account via mediation of the impurity.

The dimension of the complete Hilbert space is given by $2^{N}$ with $N$ the number of spins, and hence grows exponentially with the number of spins. When the size of all clusters is chosen equal, the number of degrees of freedom of the BPDMET ansatz is given by $2^{N_{s}}(2^{N_{s}}-1)(N_{C}-1)+2^{N_{s}}$. $N_{s}$ is the number of spins in a cluster and $N_{C}$ is the number of clusters. When the cluster size $N_{s}$ is kept fixed, the degrees of freedom scale linearly with the number of spins (or the number of clusters $N_{C}$). This linear scaling clearly is the major advantage of the BPDMET ansatz.

II.2 Optimizing the wave function

Thanks to the linear scaling, one can target the state $\ket{\Psi}$ (Eq. (4)) variationally. Optimization of the block product states and finding the ground state of the impurity model proceeds in an iterative way. At each step of the iteration, a bath cluster is chosen and its state corresponding to a certain impurity state is optimized. This way, a large number of coefficients of the variational wave function is kept fixed, and only a restricted number of coefficients is optimized in each step. The variational wave function within the impurity model can be written as

[TABLE]

with $C$ being the different bath clusters. For every possible wave function of this form, we can take

[TABLE]

by absorbing appropriate factors in the $a_{i}$’s. Even more, when the wave function is normalized,

[TABLE]

will also be satisfied.

The coefficients of the variational wave function are obtained with a restricted optimization. All coefficients are fixed, except for the $a_{i}$’s and $b^{i_{0}}_{B_{0}\beta}$, i.e. the $b$-coefficients corresponding with the impurity state $i_{0}$ and a chosen bath cluster $B_{0}$. By looping over the different $i_{0}$’s and bath clusters $B_{0}$ we optimize the DMET wave function iteratively.

We now rewrite the wave function given by Eq. (5) as

[TABLE]

Since optimization happens over $a_{i}$ with $i\neq i_{0}$ and $a_{i_{0}}b^{i_{0}}_{B_{0}\beta^{\prime}}$, every iteration is equivalent to a diagonalization in the low-dimensional subspace spanned by $\{\ket{\phi_{\alpha}},\ket{\phi_{\beta}}\}$, with

[TABLE]

To find the optimal solution with every iteration, the following Lagrangian is minimized within this restricted Hilbert space:

[TABLE]

yielding a linear eigenvalue problem.

The Lagrangian multiplier $\lambda$ is the variational energy of the wave function. Within every iteration, the solution corresponding to the smallest $\lambda$ is chosen. It is clear that the solution of the previous iteration can still be chosen within the freedom of the parameters in the current iteration. Because of this, the minimal $\lambda$-value obtained in the current iteration is at least as small as the $\lambda$ value of the previous iteration. Since $\lambda$ decreases with every two consecutive iterations, we converge to a minimal $\lambda$ value, although it is not guaranteed to be the global minimum of the energy. The complexity of a major iteration (i.e. an iteration over all $i_{0}$-values and over all bath clusters $B_{0}$) is of the order $\mathcal{O}(N_{C}^{2})$ with $N_{C}$ the number of clusters. The number of major iterations needed up to convergence may increase when increasing the number of clusters. However, the scaling of the problem when enlarging the number of spins is of course more favorable than the exponential scaling of the exact diagonalization, as long as the size of the clusters does not change.

When changing the size of the clusters, the algorithm scales exponentially due to the exponential scaling of the restricted Hilbert space chosen every minor iteration step. Even more, the number of impurity states also grows exponentially, making the number of minor iteration steps in every major iteration step also blow up. The scaling of every major iteration step is approximately $\mathcal{O}(2^{4N_{s}})$. Making the clusters bigger results quickly in high computational times.

Details of the calculations are investigated in more depth in the Supplemental Material.sup

II.3 Spin-state superposition in the impurity

When looking at the Schmidt decomposition of the exact wave function in a small impurity and a larger environment (Eq. (1)), the set of impurity states $\{\ket{\alpha_{i}}\}$ can be any basis of the Hilbert space restricted to the impurity. Corresponding environmental states will always be found. When approximating the environmental states by block product states, this freedom of choice is lost. The choice of the impurity states influences the corner of the Hilbert space the BPDMET algorithm optimizes in. In the original BPDMET method,Fan and Jie (2015) the states of the impurity $\{\ket{\alpha_{i}}\}$ are given by the natural basis for a spin system, e.g. $\{\ket{\downarrow\cdots\downarrow\downarrow},\ket{\downarrow\cdots\downarrow\uparrow},\ket{\downarrow\cdots\uparrow\downarrow},\cdots,\ket{\uparrow\cdots\uparrow\uparrow}\}$ for $S=\frac{1}{2}$. This choice is arbitrary and influences the obtained results. In this section we extend the BPDMET wave function enabling it to find an optimal set of orthonormal impurity states. Orthonormality is imposed in order to keep the simplifications in the calculation of $\matrixelement{\frac{\partial\Psi}{\partial z*}}{\hat{H}}{\Psi}$, as described in the Supplemental Material.sup

The adapted BPDMET ansatz is again given by Eq. (3). However, the impurity states are now superpositions of the natural basis states: $\ket{\alpha_{i}}=\sum\limits_{m}U_{im}\ket{m}$ with $\ket{m}=\{\ket{\downarrow\cdots\downarrow\downarrow},\ket{\downarrow\cdots\downarrow\uparrow},\cdots\}$ and $U$ a unitary matrix. The BPDMET ansatz is thus given by:

[TABLE]

Now the unitary matrix has to be optimized as well. We simply extend the original algorithm: the block-product states and the unitary matrix are optimized successively until convergence is obtained.

The unitary matrix is optimized by minimizing the variational energy $\lambda$ through successive Jacobi-rotations.Poelmans et al. (2015)

The optimization scheme for BPDMET with spin-state superposition now looks like:

Initialization of of the impurity states $(U)$ and the BPS $(b^{i}_{C\beta})$ 2. 2.

Optimization of the BPDMET wave function:

1. (a)

Optimization of the BPS: Loop over bath clusters $B_{0}$ and impurity states $i_{0}$ and solve Eq. (9) keeping appropriate parameters fixed until convergence 2. (b)

Optimization of the impurity states through successive Jacobi rotations. 3. (c)

Restart from step 2a until convergence.

As a convergence criterion both the variational energy $\lambda$ and the BPDMET energy $E$ can be used. The BPDMET energy (Eq. 13) is an alternative way to calculate the energy of the system within the DMET framework and will be introduced in the next section. In this paper we choose the BPDMET energy, as it converges somewhat more slowly than the variational energy. The faster convergence of the variational energy is clear since it is quadratically dependent on the error of the wave function, while generally, the expectation value of a property (which Eq. (13) represents) has a linear dependency as shown in the appendix.

II.4 Expectation values

In this section, the calculation of expectation values within the BPDMET framework is discussed. The method of calculation is equivalent to the method used in DMET as presented by Wouters et al.Wouters et al. (2016) In BPDMET, we divide the lattice into different clusters and choose one cluster as the impurity cluster. In this paper the division happens in such way that all clusters are equivalent with respect to the lattice symmetry. All clusters can be transformed into each other by using a translation or rotation for which the lattice is invariant. By picking one cluster $C$ as impurity and calculating its corresponding BPDMET wave function $\ket{\Psi_{C}}$, we immediately know all the BPDMET wave functions corresponding to the other choices of the impurity. This is a great advantage as far as computational time is concerned.

Wouters et al.Wouters et al. (2016) noted the fundamental difference between local and nonlocal operators. Local operators act within one impurity while nonlocal operators do not. Just like in the original DMET framework, expectation values for local operators are quite straightforward, while expectation values for nonlocal operators require some inventiveness. When a local operator $\hat{A}$ only acts upon cluster $C$, its expectation value can be calculated by: $\expectationvalue*{\hat{A}}=\expectationvalue{\hat{A}}{\Psi_{C}}$, where $\ket{\Psi_{C}}$ is the calculated BPDMET wave function with cluster $C$ chosen as impurity. Note however that also operators consisting of summations of local operators impose no problem. For example, the expectation value of the total spin in the $z$-direction is given by:

[TABLE]

where $\hat{S}_{\mathrm{tot_{C}}}^{z}$ is the total spin in the $z$-direction restricted to sites belonging to cluster $C$. Since all $\ket{\Psi_{C}}$ are equivalent, the summation over the different clusters is omitted in the last step, where $N_{C}$ is the number of clusters in the system.

For nonlocal operators the original DMET framework suggests splitting these operators appropriately.Wouters et al. (2016) The expectation values of interest for the spin lattice systems studied in this paper are given by a summation of scalar products of spin operators. The expectation value can thus be written as the sum of the expectation values of the different terms $\expectationvalue{\hat{\mathbf{S}}_{i}\cdot\hat{\mathbf{S}}_{j}}$. When both $i$ and $j$ are sites within one cluster, this expectation value is an expectation value of a local operator. However, when this is not the case, this expectation value is an expectation value of a nonlocal operator and will be calculated as:

[TABLE]

where $\ket{\Psi_{C_{i}}}$ and $\ket*{\Psi_{C_{j}}}$ are the BPDMET solutions with the impurity chosen to be the cluster of site $i$ or site $j$, respectively. Since the solutions of the BPDMET for different impurity clusters are equivalent, the calculation of expectation values of operators that respect the lattice symmetry can be simplified. Examples of these are the Hamiltonian $\hat{H}$ in Eq. (2) and the squared total spin $\hat{\mathbf{S}}^{2}_{\mathrm{tot}}$. These are given by:

[TABLE]

By calculating properties in this way, we take into account that BPDMET describes the impurity more accurately than the bath. However, this method is not variational in nature, so energies obtained through Eq. (13) are not upper bounds to the exact energy, and squared total spins can be slightly negative. From now on, this non-variational energy will be called the BPDMET energy and denoted by $E$, while the variational energy is given by the Lagrangian multiplier $\lambda$ in Eq. (9).

II.5 Tangent space and excitations

The BPDMET algorithm can also be extended for the calculation of approximate spectral functions. For the calculation of the spectral function, the Hamiltonian restricted to the tangent space of the BPDMET ansatz (Eq. 10) is diagonalized. Shifting of weights between the different $b^{i}_{C\beta}$ and $U_{im}$ parameters is possible without changing the actual wave function. It is thus clear that the parametrization is redundant. By introducing a set of restrictions, the redundancy can be lifted and normalization of the wave function can be imposed:

[TABLE]

The tangent space is constructed by taking these restrictions into account and differentiating with respect to the nonredundant parameters. Diagonalization of the Hamiltonian restricted to the tangent space now amounts to solving a generalized eigenvalue problem. Eigenvectors with eigenvalues close to zero of the overlap matrix are projected out.

The spectral function is given by

[TABLE]

where $\hat{X}$ is a perturbation operator connecting the ground state with the excited states. By restricting to the tangent space, this can be rewritten as

[TABLE]

where $\phi_{n}$ and $E_{n}$ are the eigenvectors of the generalized eigenvalues problem and their corresponding energies. For $E_{n}$ and $E_{0}$, both the variational energy $\lambda$ and the BPDMET energy $E$ in Eq. (13) can be used. It will be shown that the latter choice yields inferior results.

In this paper, the spectral function is calculated by searching all the excitations within the tangent space. Another option is through solving the linear response equation given by:

[TABLE]

Both methods can be done through sparse iterative solvers. However, since the BPDMET ansatz has a rather low number of parameters, the tangent space has a low dimension, making explicit solving feasible.

II.6 Connection with Tensor Networks

When looking at the BPDMET ansatz, it is clear that it can be represented by the tensor network depicted in Fig. 2. The central tensor corresponds to the impurity cluster and the remaining tensors to the bath clusters. It should be noted, that these bath tensors can be chosen differently from each other, as can be readily seen from Eq. (10). The spin degrees of freedom are combined into one physical index per cluster. The impurity tensor has virtual indexes connected to every bath tensor. It is therefor a very high rank tensor (one virtual index per bath cluster). However, the impurity tensor is heavily restricted, making the tensor network manageable. The BPDMET high rank impurity tensor can be represented as

[TABLE]

with $U$ a unitary matrix. The traditional technique to make TNS manageable is the truncation of the virtual dimension (e.g. as in DMRG and PEPS). In BPDMET, the TNS is made manageable by imposing restrictions on the impurity tensor. There are no truncations in the virtual dimension. Clarifying the link with tensor networks can facilitate the construction of different ansatzes for the bath states.

III Results

In this section, results of the BPDMET method are discussed for two different models. In section III.1, the 2D Heisenberg model on the square lattice with nearest-neighbor (NN) and next-nearest-neighbor (NNN) interaction is studied. The BPDMET results obtained by Fan et al. Fan and Jie (2015) are reproduced and compared with the results when spin-state superposition is introduced (see section II.3). In addition, different order parameters are calculated. In section III.2, the BPDMET algorithm is applied to the Kitaev-Heisenberg model on a Honeycomb latticeChaloupka et al. (2010) with 24 spins and compared with exact results.

Spectral functions obtained through diagonalization in the tangent space as discussed in section II.5 are also presented.

Note that when dividing the complete lattice into two clusters, one impurity and one bath cluster, that the BPDMET ansatz (Eq. (10)) yields a redundant parametrization of the complete Hilbert space. In this situation, BPDMET should coincide with exact diagonalization (ED). Comparison with ED for small systems with only one bath cluster allows us to check the correctness of the algorithms and implementations.

III.1 NN and NNN interaction on the square lattice

The first system under consideration is the square lattice with NN and NNN Heisenberg interactions. The Hamiltonian of this system is given by

[TABLE]

with $\langle i,j\rangle$ denoting NN sites, $J_{1}$ the NN interaction strength, $\langle\langle i,k\rangle\rangle$ denoting NNN sites and $J_{2}$ the NNN interaction strength. We restrict ourselves to anti ferromagnetic (AF) interactions ($J_{1},J_{2}\geqslant 0$). This is the same system as investigated by Fan et al.Fan and Jie (2015) Three distinct phases are present in the infinite lattice. At low NNN interaction, the ground state is in a Néel phase, this is a long-range-ordered phase. At $J_{2}/J_{1}\approx 0.4$ a phase transition from this Néel phase happens to a disordered quantum paramagnetic phase. When tuning the system to stronger NNN interactions, the system undergoes a transition to another long range ordered phase at $J_{2}/J_{1}\approx 0.6$. This is the collinear phase. The nature of the intermediate paramagnetic phase is still undecided. Multiple interpretations for this phase have been proposed, such as spin liquids and valence bond states like the columnar and staggered dimer valence bond crystals and the plaquette resonating valence bond (PRVB) state.Chandra and Doucot (1988); Gelfand et al. (1989); Read and Sachdev (1989); Takano et al. (2003); Darradi et al. (2008); Murg et al. (2009); Jiang et al. (2012); Mezzacapo (2012); Yu and Kao (2012); Li et al. (2012b); Wang et al. (2013); Hu et al. (2013); Doretto (2014); Gong et al. (2014) In Ref. Fan and Jie, 2015, it is shown that BPDMET calculations suggests no rotational symmetry breaking at the intermediate phase, and evidence is found in favor of the PRVB.

The BPDMET energies (see Eq. (13)) are calculated for a $8\times 8$ square lattice with periodic boundary conditions, i.e. an $8\times 8$ square lattice on a torus. The lattice is divided into 16 equal $2\times 2$ clusters of which one is chosen as impurity. Ground-state energies are calculated as explained in section II.4. In Fig. 3, the converged energy values for BPDMET with random initialization are given. At high and low $J_{2}/J_{1}$ convergence of the BPDMET algorithm happens quite consistently to the same energy values. When $J_{2}/J_{1}\in[0.6,0.8]$, multiple energy values are found and the algorithm converges to a variety of local minima. For purposes of reproducibility, we will make use of sweeps (Fig. 4 and Fig. 5). A sweep starts in a region where convergence is consistent to the same minimum, and sweeps through the parameter region using previous converged results as initialization. Sweeps can be done from low to high $J_{2}/J_{1}$ or vice versa. The BPDMET energy is used as selection criterion for the optimal solution. It should be stressed that using the random initializations (40 runs per parameter value) we never found a lower energy than the minimum energies found through sweeps from the left or right. This suggests that the sweep finds more optimal solutions that are hard to find with random initialization. This poses some justification for the use of sweeps.

Introduction of spin-state optimization in BPDMET yields significant changes in the calculated energy, as can been seen in Fig. 4(a). In the two ordered phases, minor changes in energy per spin are observed. In the intermediate paramagnetic phase, larger changes are visible. It is good to note that the BPDMET energy is not variational so a lowering in energy is not necessarily a net improvement of the energy. However, a study has been done on the 40-spin lattice (Fig. 5) and compared with exact results obtained in Ref. Richter and Schulenburg, 2010. The boundary conditions are chosen in the same way as in Ref. Richter and Schulenburg, 2010. When introducing the spin-state optimization for the 40-spin system, there is a substantial improvement in energy observed in the intermediate paramagnetic phase. This makes us confident that the substantial change in BPDMET energy through introduction of spin-state optimization is also a net improvement for the 64-spin lattice.

When looking at the variational energy $\lambda$ (Fig. 4(a)), only a small correction occurs with the introduction of spin-state optimization. At $J_{2}=0$, $\lambda$ changes from $-0.5939$ per spin to $-0.5940$ per spin, while the non-variational BPDMET energy $E$ changes from $-0.6657$ per spin to $-0.6678$ per spin. The change in the BPDMET energy is clearly much larger than the change in the variational energy, and this was found to occur for all $J_{2}/J_{1}$ values. The improvement of the results is largely contained in the impurity spins. The energy per spin at $J_{2}=0$ ($E=-0.667840$) corresponds reasonably well with results obtained by quantum Monte Carlo (QMC) for the $8\times 8$ lattice ($E=-0.673487$).Sandvik (1997)

Fan et al.Fan and Jie (2015) calculated the energy per spin at $J_{2}=0$ as $E=2\expectationvalue{\hat{\mathbf{S}}_{1}\cdot\hat{\mathbf{S}}_{2}}{\psi}$, making use of rotational symmetry present in the Néel phase. This is twice the bond energy between two NN spins in the impurity cluster. Using this expression for the energy per spin, $E=-0.66917$ and $E=-0.67045$ without and with spin-state optimization are obtained. Correspondence with the QMC solution clearly improves using this alternative energy expression.

The Néel and the collinear order parameters are calculated in order to identify the different phase transitions. These order parameters are given by

[TABLE]

where $\mathbf{Q}$ is given by $(\pi,\pi)$ for the Néel parameter and by $(\pi,0)$ and $(0,\pi)$ for the collinear parameter in $x$- and $y$-direction. $R_{i}$ is the position vector of the spin site $i$ and $N$ is the total number of spins. For an infinite lattice rotational symmetry breaking will occur. However, the exact solution for the finite lattice (for instance a $4\times 4$ lattice) yields a collinear order parameter that is equal in the $x$- and $y$-direction. We notice that BPDMET finds a rotation symmetry broken solution, even in finite lattices. In Fig. 4(b) the order parameters are shown. The largest collinear parameter is plotted here.

When introducing spin-state optimization, small changes are noticeable in the order parameters. The Néel order parameter at $J_{2}=0$ ($M^{2}_{N}=0.19825$ and $M^{2}_{N}=0.19670$ without and with spin-state optimization) corresponds well with results obtained through QMC ($M^{2}_{N}=0.17784$).Sandvik and Evertz (2010); Sandvik (1997) Phase transitions at $J_{2}/J_{1}\approx 0.4$ and $J_{2}/J_{1}\approx 0.62$ can be observed, in correspondence with previous studies. The location of these phase transitions does not change with the introduction of spin-state optimization. A strong Néel order is observed at low $J_{2}/J_{1}$ while a strong collinear order is observed at high $J_{2}/J_{1}$. In the intermediate region, both order parameters stay rather small, but do not completely vanish due to finite size effects.

III.2 The Kitaev-Heisenberg model

The applicability of BPDMET is not limited to square lattices, other lattices are also equally feasible. In this section we consider the Kitaev-Heisenberg modelChaloupka et al. (2010); Jiang et al. (2011); Reuther et al. (2011); Price and Perkins (2012); Chaloupka et al. (2013); Oitmaa (2015) on the honeycomb lattice. We study a 24-spin lattice with periodic boundary conditions (see Fig. 6). For this system exact diagonalization of the system is still feasible and the BPDMET method can be benchmarked. This model is a mixture of the Kitaev modelKitaev (2006) and the Heisenberg model on the honeycomb lattice and spin interactions are given by the Hamiltonian:

[TABLE]

$J_{1}$ represents the strength of the anisotropic Kitaev interaction while $J_{2}$ is the strength of the Heisenberg part. Every spin is connected through exactly one $x$-, $y$- and $z$-link to a neighboring spin. The different links are shown in Fig. 6.

The Kitaev and Heisenberg interaction are parametrized as $J_{1}=2\alpha$ and $J_{2}=1-\alpha$. In the interval $\alpha\in[0,1]$, three phases occur.Chaloupka et al. (2010) At low $\alpha$ a Néel AF phase is observed, at intermediate $\alpha$ a stripy AF phase, and at high $\alpha$ a quantum spin liquid is observed. The phase transitions are expected at $\alpha\approx 0.4$ and $\alpha\approx 0.8$.Chaloupka et al. (2010); Jiang et al. (2011); Reuther et al. (2011); Price and Perkins (2012); Chaloupka et al. (2013); Oitmaa (2015) At the intermediate point, $\alpha=0.5$, the system is exactly solvable through the use of a rotated basis.Chaloupka et al. (2010) This rotated basis is defined by dividing the lattice into 4 sublattices and defining different rotated spin operators $\hat{\tilde{\mathbf{S}}}$ on these sublattices as can be seen in Fig. 6. At this intermediate point, the Hamiltonian is reduced to the ferromagnetic Heisenberg model in the rotated basis, and exact solutions are given by states with maximal total spin in the rotated basis (e.g. $\ket{\uparrow\uparrow\dots\uparrow}$). The two trivial states with maximal total spin can be clearly represented by the BPDMET ansatz. The basis transformation suggested in Ref. Chaloupka et al., 2010 leaves the corner of the Hilbert space described by the BPDMET ansatz unchanged.Gunst (2016) In the original basis it is therefore also possible to find exact solutions at this intermediate point through BPDMET.

BPDMET allows freedom in the choice of the impurity and bath clusters. Three different types of clusters that can cover the entire lattice (see Fig. 7) have been examined.

When randomly initializing the BPDMET algorithm, convergence happens quite consistently in the two AF phases ($\alpha<0.8$); however, in the spin liquid phase, the calculations converge to a wide variety of local minima as can be seen in Fig. 8 with the calculation of $\hat{\tilde{\mathbf{S}}}^{2}_{\mathrm{tot}}$ for the S-shaped cluster. Similar results were obtained for the star shaped and hexagonal clusters. It was found that results obtained with the S-shaped cluster are slightly inferior to the other two types of clusters in describing the phase transitions. We will therefore only present these random initialization results for the S-shaped cluster and in the remainder of the paper only the star shaped and hexagonal clusters are discussed.

Sweeps, either from the left or from the right, are again used to ensure reproducibility of the results.

Comparison of the BPDMET results for the 24-spin system with the exact results shows a difference between the star-shaped and hexagonal clusters. Both clusters describe the two AF phases well. When using the star-shaped cluster a phase transition occurs at $\alpha\approx 0.8$ (Fig. 9). Although the phase transition to the spin liquid is observed, the spin liquid itself is poorly described through BPDMET as can be inferred from the calculated BPDMET energies (Fig. 9(a)), spin properties (Fig. 9(b)) and correlation functions (Fig. 9(c)). The phase transition to the spin liquid gets more pronounced when introducing spin-state optimization, but the spin liquid itself is still not represented adequately.

In general, the BPDMET ansatz can capture a larger corner of the Hilbert space with increasing impurity size, so it is expected that the hexagonal cluster performs better than the 4-spin clusters. When using the hexagonal cluster, the phase transition towards a spin liquid is detected and the obtained properties of the spin liquid are in good correspondence with the exact diagonalization at $\alpha=1$ (Fig. 10). However, the phase transition happens at $\alpha\approx 0.92$, which is not the right value. Calculations with spin-state optimization have not been performed for this type of cluster since they were already quite intensive without the optimization.

BPDMET allows to investigate larger systems than the 24 spin lattice. When expanding the honeycomb lattice by one extra layer, a 54 spin honeycomb lattice is obtained which is however not coverable with the S-shaped and star shaped clusters. Also, a rotated basis respecting periodic boundary conditions as proposed in Ref. Chaloupka et al., 2010, cannot be found. With another extra layer, we get a 96 spin honeycomb lattice that is coverable with the three types of clusters (Fig. 7), and is also consistent with the concept of the rotated basis. No shift in the location of the phases is found when extending to 96 spins.

III.3 Spectral functions

In section II.5, a method is introduced for the BPDMET to find the excitations and the spectral functions through the tangent space. For the calculation of the spectral function Eq. (19) is used, with $\eta=0.01$. In Eq. (19), two possible energy values for $E_{n}$ can be used, i.e. the variational energy $\lambda$ and the BPDMET energy as given in Eq. (13).

For the Kitaev-Heisenberg model at $\alpha=0.5$, the exact solution can be found through BPDMET. We expect thus that the spectral function will be reproduced quite well. Since the ground state in $\alpha=0.5$ is degenerate, the equation for the spectral function is changed to

[TABLE]

where $d$ is the dimension of the degenerate ground state space, $\Delta E_{j}$ is the energy difference between the state $\ket{\phi_{j}}$ and the energy of the degenerate ground states, and the index $i$ sums over all ground states while $j$ sums over the excited states. This expression is independent of the chosen basis in the degenerate ground state space.

To calculate the spectral function through BPDMET we simply perform the calculations multiple times with random initializations and average out over the obtained spectral functions. Residues between the ground states are not taken into account in the spectral function, as can be seen in Eq. (25). To make sure that the tangent space method does not take these residues into account, only eigenvectors with an eigenvalue significantly different from the ground state energy are used in the summation in Eq. (19). When diagonalizing the Hamiltonian in the tangent space, a few eigenvalues very close to the ground state energy are obtained (within a 0.03 margin), which are separated from the other excitations by a clear gap (the next eigenvalues are found at $\approx 0.25$). These eigenstates are left out of the spectral function since they originate from the degenerate ground-state solutions within the tangent space.

We calculate the spectral function with $\hat{S}^{+}_{i}$ used as perturbation $\hat{X}$, i.e. the raising operator on spin $i$. In Fig. 11(a), it is observed that the tangent space method finds practically the exact spectrum when using the variational energy. Some small differences in peak location and amplitude are visible.

When using the non-variational DMET energy in Eq. (19), all correspondence with the exact solution is lost (see Fig. 11(b)). It is therefore clear that the variational energy should be used in the calculation of the spectral function in contrary to the calculation of ground state properties. There, the non-variational BPDMET energy has proven to be superior.

In the previous example, BPDMET is able to find the exact ground-state wave function. When this is not the case, the calculation of the spectral function through tangent space diagonalization can be inadequate. This can be seen in Fig. 12, where the spectral function for the same Kitaev-Heisenberg model is calculated at $\alpha=0.25$. BPDMET is, in this case, unable to find the exact ground state. For this system, correspondence of the BPDMET spectral function with exact results is lost.

IV Conclusion

The BPDMET method is used to study the spin-1/2 anti ferromagnetic Heisenberg model on the square lattice with 64 spins as originally done in Ref. Fan and Jie, 2015 and the spin-1/2 Kitaev Heisenberg modelChaloupka et al. (2010) on the honeycomb lattice with different sizes. A systematic approach for the calculation of properties within the BPDMET framework is introduced. Spin-state superposition in the impurity has been added to BPDMET yielding improved results for the energy profile and properties. The calculation of excited states and the spectrum through diagonalization in the tangent space of the BPDMET ground state is investigated, but provides somewhat unsatisfactory results. In any case, it has been shown that for the calculation of the spectral function the variational energy $\lambda$ should be used and not the BPDMET energy.

For the AF Heisenberg model on the square lattice, order parameters are calculated and the right phases are detected in the system. For the Kitaev Heisenberg model, different types of clusters are used. The results for a 24-spin honeycomb lattice are compared with the exact solution. It is clear that the results are dependent on the cluster shape and that some clusters are more appropriate for certain phases than others. However, none of the cluster shapes are able to represent the spin liquid regime. Only with hexagonal clusters, the spin liquid phase is detected, but the phase transition happens at a wrong value of the coupling parameters. BPDMET enables one to investigate larger systems unreachable with exact diagonalization. The ground state for the 96-spin honeycomb lattice is calculated and it is found that the position of the phase transitions does not change when enlarging the system.

Acknowledgements.

K.G. and S.D.B. acknowledge support from the Research Foundation Flanders (FWO Vlaanderen).

Appendix A Convergence

For a variational optimization, we show that the error on the energy scales quadratically with the error on the wave function. For expectation values of general properties, the error scales linearly with the wave function error. In order to see this, we decompose our approximate wave function into a component parallel to the exact ground state and a perpendicular error-component, i.e.Wouters and Van Neck (2014)

[TABLE]

with $\innerproduct{\Psi_{0}}{\Psi_{\mathrm{error}}}=0$ and $\epsilon$ a measure for the error. The error on the wave function is

[TABLE]

and thus linear in $\epsilon$ for small errors. We get

[TABLE]

for the error on the expectation value of general properties. We see that this error is linear with $\epsilon$. Since $\ket{\Psi_{0}}$ is an eigenstate of $\hat{H}$ and $\innerproduct{\Psi_{\mathrm{error}}}{\Psi_{0}}=0$, the linear term vanishes for the error on the energy and the leading error term is thus quadratic in $\epsilon$. In BPDMET, an optimization is performed in a restricted Hilbert space during each minor iteration step. During this optimization, the optimal solution within the restricted space is chosen by diagonalizing the effective Hamiltonian in this restricted Hilbert space (as discussed in section II.2). The updated wave function is therefore an eigenstate (corresponding to the minimal eigenvalue) of the effective Hamiltonian (not necessarily of the full Hamiltonian).

In each minor iteration, the wave function $\ket{\Psi}$ is updated to the ground state $\ket{\Psi_{0}}$ of the effective Hamiltonian. The change during each minor iteration of the variational energy $\lambda$ is given by Eq. (28) with $\hat{Q}=\hat{H}_{\mathrm{eff}}$ and for the same reason as explained above, will have a quadratic leading term. The effective non-variational BPDMET energy operator (i.e. the BPDMET energy operator mapped to the restricted Hilbert space), will have a linear leading term.

When the BPDMET algorithm is close to a local minimum and close to convergence, the correction parameter $\epsilon$ will be small. The corrections on the variational energy $\lambda$ will be quadratic in $\epsilon$, while the corrections on the BPDMET energy will be linear in $\epsilon$. Using the BPDMET energy as convergence criterion is therefore more stringent.