Numerical linked cluster expansions for quantum quenches in one dimensional lattices

TL;DR

This paper applies numerical linked cluster expansions (NLCEs) to analyze one-dimensional lattice systems in thermal equilibrium and after quantum quenches, demonstrating NLCEs' effectiveness in the thermodynamic limit and comparing different expansion methods.

Contribution

It introduces and compares two NLCE methods for quantum lattice systems, highlighting the site-based NLCE's advantages and divergence issues, and demonstrates NLCEs' superiority over exact diagonalization.

Findings

NLCEs provide accurate thermodynamic limit results.

Site-based NLCE works best in thermal equilibrium.

NLCEs outperform exact diagonalization in periodic systems.

Abstract

We discuss the application of numerical linked cluster expansions (NLCEs) to study one dimensional lattice systems in thermal equilibrium and after quantum quenches from thermal equilibrium states. For the former, we calculate observables in the grand canonical ensemble, and for the latter we calculate observables in the diagonal ensemble. When converged, NLCEs provide results in the thermodynamic limit. We use two different NLCEs - a maximally connected expansion introduced in previous works and a site-based expansion. We compare the effectiveness of both NLCEs. The site-based NLCE is found to work best for systems in thermal equilibrium. However, in thermal equilibrium and after quantum quenches, the site-based NLCE can diverge when the maximally connected one converges. We relate this divergence to the exponentially large number of clusters in the site-based NLCE and the behavior of the weights of observables in those clusters. We discuss the effectiveness of resummations to cure the divergence. Our NLCE calculations are compared to exact diagonalization ones in lattices with periodic boundary conditions. NLCEs are found to outperform exact diagonalization in periodic systems for all quantities studied.

Tables

Table 1. Table 1: Total number of clusters and of topological clusters in the site-based NLCE.

$L$	Symmetry distinct clusters	Total number of clusters
1	1	1
2	2	2
3	3	4
4	6	8
5	10	16
6	20	32
7	36	64
8	72	128
9	136	256
10	272	512
11	528	1024
12	1056	2048
13	2080	4096
14	4160	8192
15	8256	16384
16	16512	32768
17	32896	65536

Equations

$$\hat{H}$$

$$\mathcal{O}/N=\sum _c\mathcal{L}(c)\times W_{\mathcal{O}}(c),$$

$$W_{\mathcal{O}}(c)=\mathcal{O}(c)-\sum _{s\subset c}{W}_{\mathcal{O}}(s).$$

$$S_L=\sum _{\{c_L\}}\mathcal{L}(c_L)\times W_{\mathcal{O}}(c_L).$$

$${\mathcal{O}}_l/N=\sum _{L=1}^l{S}_L.$$

$${ϵ}_n^{(-1)}$$

$${ϵ}_n^{(k)}$$

$${\text{Wynn}}_k({\mathcal{O}}_l)={ϵ}_{l-2k}^{2k}.$$

$${\text{Euler}}_k({\mathcal{O}}_l)=\sum _{L=0}^{l-k}{S}_L+(-1{)}^{l-k+1}\sum _{L=0}^k\frac{1}{2^{L+1}}{T}_{k,L}^l,$$

$$\text{where}{T}_{k,L}^l=\left[\sum _{j=0}^L\left(\genfrac{}{}{0ex}{}{L}{j}\right)S_{l-k+L-j+1}\right].$$

$$\hat{U}$$

$${\hat{m}}_k=\frac{1}{N}\sum _{j,j^{\prime }}{e}^{ik(j-j^{\prime })}{\hat{b}}_j^{†}{\hat{b}}_{j^{\prime }}.$$

$${\Delta }_l(\mathcal{O})=\left|\frac{{\mathcal{O}}_{l_{\text{max}}}-{\mathcal{O}}_l}{{\mathcal{O}}_{l_{\text{max}}}}\right|,$$

$${\Delta }_L(\mathcal{O})=\left|\frac{{\mathcal{O}}_{L_{\text{max}}}-{\mathcal{O}}_L}{{\mathcal{O}}_{L_{\text{max}}}}\right|.$$

$${\hat{\rho }}_c^{\text{GE}}=\frac{\exp (-\beta {\hat{H}}_c-\mu {\hat{N}}_c)}{\text{Tr}[\exp (-\beta {\hat{H}}_c-\mu {\hat{N}}_c)]},$$

$${\hat{\rho }}_c^I=\frac{\exp (-{\beta }_I{\hat{H}}_c^I-{\mu }_I{\hat{N}}_c)}{\text{Tr}[\exp (-{\beta }_I{\hat{H}}_c^I-{\mu }_I{\hat{N}}_c)]}.$$

$${\hat{\rho }}_c(\tau )$$

$${\hat{\rho }}_c^{\text{DE}}=\underset{{\tau }^{\prime }\to \infty }{\lim }\frac{1}{{\tau }^{\prime }}{\int }_0^{{\tau }^{\prime }}{\hat{\rho }}_c(\tau )d\tau =\sum _m(⟨m|{\hat{\rho }}_c^I|m⟩)|m⟩⟨m|,$$

$$\text{Tr}[{\hat{\rho }}^{\text{GE}}\hat{H}]$$

$$\text{Tr}[{\hat{\rho }}^{\text{GE}}\hat{N}]$$

$${\delta }_l(\mathcal{O})=\left|\frac{{\mathcal{O}}_{18}^{\text{GE}}-{\mathcal{O}}_l^{\text{DE}}}{{\mathcal{O}}_{18}^{\text{GE}}}\right|$$

$${\delta }_L(\mathcal{O})=\left|\frac{{\mathcal{O}}_{18}^{\text{GE}}-{\mathcal{O}}_L^{\text{DE}}}{{\mathcal{O}}_{18}^{\text{GE}}}\right|,$$

$$W(\mathcal{O},L)=\frac{{\sum }_{\{c_L\}}\mathcal{L}(c_L)\times W_{\mathcal{O}}(c_L)}{{\sum }_{\{c_L\}}\mathcal{L}(c_L)}.$$

Full text

I Introduction

The study of strongly interacting lattice systems is in general very challenging and only special (integrable) models admit exact analytical solutions Korepin et al. (1993); Sutherland (2004); Baxter (2007); Cazalilla et al. (2011). In this context, numerical calculations based on series expansions have proved to be very useful Domb and Green (1972); Guttmann (1989); Oitmaa et al. (2006). An even more challenging question that has attracted much attention recently is the far-from-equilibrium quantum dynamics of isolated many-particle lattice systems D’Alessio et al. (2016); Eisert et al. (2015); Polkovnikov et al. (2011). This has been motivated by recent experiments with ultracold atoms in optical lattices, in which it is possible to create strongly interacting quantum lattice systems with tunable parameters in nearly isolated environments. The far-from-equilibrium dynamics can be generated, e.g., by sudden changes of the depth of the optical lattice Greiner et al. (2002); Will et al. (2010, 2015) or by engineering special initial states Kinoshita et al. (2006); Gring et al. (2012); Trotzky et al. (2012); Langen et al. (2015); Clos et al. (2016); Kaufman et al. (2016). A question that has been addressed in the latter set of experiments is whether, under (nearly) unitary dynamics, observables relax to time independent values that can be described using traditional statistical mechanics.

The same question has been studied theoretically mostly for the quantum dynamics of pure states D’Alessio et al. (2016); Eisert et al. (2015); Polkovnikov et al. (2011). Results from numerical simulations in a variety of nonintegrable lattice models in one and two dimensions have indicated that, under unitary dynamics, few-body observables relax to time-independent values that are described using traditional ensembles of statistical mechanics Rigol et al. (2008); Rigol (2009a); *rigol2009quantum; Eckstein et al. (2009); Rigol and Santos (2010); Bañuls et al. (2011); Khatami et al. (2013); Zangara et al. (2013); Sorg et al. (2014), a phenomenon we call thermalization and which has been the central topic of a recent review D’Alessio et al. (2016). On the other hand, few-body observables in integrable systems relax to time-independent values that are described by the generalized Gibbs ensemble Rigol et al. (2007a); Cazalilla (2006); Cassidy et al. (2011); Calabrese et al. (2011); Ilievski et al. (2015), a phenomenon we call generalized thermalization and which has been a central topic of several recent reviews Vidmar and Rigol (2016); Essler and Fagotti (2016); Cazalilla and Chung (2016); Caux (2016); Ilievski et al. (2016). Thermalization in nonintegrable systems has been understood to be the result of eigenstate thermalization Srednicki (1994); Deutsch (1991); Rigol et al. (2008); Rigol and Srednicki (2012); D’Alessio et al. (2016), while generalized thermalization in integrable systems has been understood to be the result of generalized eigenstate thermalization Cassidy et al. (2011); Vidmar and Rigol (2016).

While integrable models can be studied using exact analytic approaches Vidmar and Rigol (2016); Essler and Fagotti (2016); Cazalilla and Chung (2016); Caux (2016); Ilievski et al. (2016), nonintegrable ones require the use of numerical techniques such as full exact diagonalization Rigol et al. (2008); Rigol (2009a); *rigol2009quantum; Rigol and Santos (2010); Khatami et al. (2013); Zangara et al. (2013); Sorg et al. (2014), for which one can access arbitrarily long times but is limited to small system sizes, or density matrix renormalization group Schollwöck (2005, 2011) and dynamical mean field theory Georges et al. (1996); Aoki et al. (2014) like approaches, for which the limitation is not in the system size but in the accessible times Eckstein et al. (2009); Bañuls et al. (2011); Sorg et al. (2014). Recently, a numerical linked cluster expansion (NLCE) approach was introduced that, when converged, allows one to calculate infinite-time averages of observables in the thermodynamic limit after quenches from initial thermal equilibrium Rigol (2014a) and pure Rigol (2014b) states. This approach was used to probe fundamental differences in quenches to and away from integrability starting from thermal states Rigol (2016), and to study quantum quenches from Néel and tilted ferromagnetic states in the spin-1/2 XXZ chain Wouters et al. (2014); Piroli et al. (2017). Its main limitation is that the expansion may fail to converge when the temperature of the initial state is low Rigol (2014a), as well as for Hamiltonian parameters of interest in quenches from pure states Wouters et al. (2014); Rigol (2014b); Piroli et al. (2017). For systems in thermal equilibrium, it has been shown that NLCEs converge to lower temperatures than traditional high-temperature expansions Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2, and exponentially faster than exact diagonalization calculations Iyer et al. (2015). NLCEs can also outperform other computational approaches in ground state Khatami et al. (2011); Yang and Schmidt (2011); Coester et al. (2015); Ixert et al. (2015) and entanglement entropy Kallin et al. (2013); Stoudenmire et al. (2014); Sherman et al. (2016) calculations.

Here, we compare two complementary NLCEs, one based on maximally connected clusters (used in previous works Rigol (2014a, 2016)) and the second one is a site-based expansion. We compare the effectiveness of both expansions in one dimensional (1D) lattices with nearest and next-nearest neighbor couplings. Furthermore, we compare the NLCE results with exact diagonalization calculations in the grand canonical ensemble in finite systems with periodic boundary conditions. We consider systems in thermal equilibrium as well as infinite-time averages after quenches from thermal states. We show that NLCEs are superior to plain exact diagonalization. When the NLCEs fail to converge, we use resummation techniques to extend their region of convergence.

The presentation is organized as follows. In Sec. II, we introduce the model Hamiltonian considered. The NLCEs used to study thermal equilibrium states and quantum quenches are discussed in Sec. III. There, we also review resummation methods used to accelerate the convergence, which were previously discussed in Ref. Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2 in the context of two-dimensional lattice systems in thermal equilibrium. The results of our calculation in 1D are presented in Sec. IV for thermal states and in Sec. V for quantum quenches. In Sec. VI, we discuss why the site-based NLCE fails to converge after quantum quenches for initial temperatures at which the thermal equilibrium results converge. A summary of our results is presented in Sec. VII.

II Model Hamiltonian

We consider the $t$-$V$-$t^{\prime}$-$V^{\prime}$ model for hard-core bosons in 1D lattices Cazalilla et al. (2011)

[TABLE]

where $\hat{b}^{\dagger}_{i}(\hat{b}_{i})$ is the hard-core boson creation (annihilation) operator at site $i$, and $\hat{n}_{i}=\hat{b}^{\dagger}_{i}\hat{b}_{i}$ is the number operator. The creation-annihilation operators obey bosonic commutation relations: $[\hat{b}_{i},\hat{b}_{j}]=0$, $[\hat{b}^{\dagger}_{i},\hat{b}^{\dagger}_{j}]=0$, $[\hat{b}_{i},\hat{b}_{j}^{\dagger}]=\delta_{ij}$, with the constraints $\hat{b}_{i}^{2}=(\hat{b}_{i}^{\dagger})^{2}=0$, which prevent multiple occupation of the lattice sites. Note that hopping and interaction terms are restricted to nearest $(t,V)$ and next-nearest $(t^{\prime},V^{\prime})$ neighbors sites.

This 1D lattice system can be represented graphically as depicted in Fig. 1. In the absence of next-nearest neighbor hoppings and interactions, the Hamiltonian (1) is integrable (it is the spin-1/2 XXZ chain in the spin language), and can be exactly solved using the Bethe ansatz Cazalilla et al. (2011). Next-nearest neighbor hoppings and interactions make the system nonintegrable.

III Numerical Linked Cluster Expansions

NLCEs are based on the linked cluster theorem, which states that an extensive quantity per site $\mathcal{O}/N$, on a lattice with $N$ sites, can be calculated as the sum of the contributions from all connected clusters that can be embedded on the lattice:

[TABLE]

where $\mathcal{L}(c)$ is the number of ways per site that cluster $c$ can be embedded on the lattice, and $W_{\mathcal{O}}(c)$ is the weight of observable $\mathcal{O}$ in cluster $c$. $W_{\mathcal{O}}(c)$ can be calculated, using the inclusion-exclusion principle, subtracting the weight of the observable in all the connected subclusters of $c$ to the value of the observable in cluster $c$:

[TABLE]

For the smallest cluster, $W_{\mathcal{O}}(c)=\mathcal{O}(c)$.

The observable for a finite cluster with density matrix $\hat{\rho}_{c}$ is calculated as $\mathcal{O}(c)=\text{Tr}[\hat{O}\hat{\rho}_{c}]$. In high-temperature expansions, $\mathcal{O}(c)$ is calculated expanding $\hat{\rho}_{c}$ in powers of the inverse temperature $\beta=1/(k_{B}T)$. (We set $k_{B}=1$ in what follows.) In NLCEs, one calculates $\mathcal{O}(c)$ numerically using full exact diagonalization.

When correlations are short ranged, the weights $W_{\mathcal{O}}(c)$ are expected to decrease rapidly with increasing the cluster size beyond the correlation length. Hence, calculating Eq. (2) to some finite order can be sufficient to estimate thermodynamic limit result to machine precision. One can reduce significantly the number of clusters to be diagonalized by identifying symmetries and topologies that relate clusters with identical expectation values of a given observable of interest. For a pedagogical introduction to NLCEs and their implementation, see Ref. Tang et al. (2013).

Another remarkable feature of NLCEs is that one has freedom to select different building blocks to construct the clusters to be used in the expansion, e.g., in the square lattice one can use bonds, sites, and squares Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2. The order of the NLCE is then determined by the largest clusters having similar characteristics, e.g., in the square lattice it could be the number of bonds, sites, or squares, depending on the building blocks chosen.

For the 1D lattice of interest in this work (see Fig 1), there are three straightforward ways of constructing the clusters. We will use two of them.

(i) NLCE based on maximally connected clusters (NLCE-M): Starting from a single site, a nearest neighbor site is added each time along with all possible bonds it can have with the rest of the cluster (see Fig. 2). This procedure generates clusters that have the maximum number of bonds for a given number of sites. The order of the “maximally connected” NLCE is set by the number of sites of the largest cluster (note that there is only one cluster for each given number of sites) Rigol (2014a, 2016). When only nearest neighbor interactions are present, this is the only NLCE possible for a 1D lattice Rigol (2014b); Wouters et al. (2014); Piroli et al. (2017). Hence, this NLCE is expected to be best suited for weak next-nearest neighbor couplings.

(ii) Site-based NLCE (NLCE-S): Starting from a single site, a nearest or next-nearest neighbor site is added one at a time along with all possible bonds it can have with the rest of the cluster (see Fig. 3). The order of the NLCE is set by the number of sites of the largest clusters (note that this time there are many clusters that can have the same number of sites). This NLCE is expected to outperform the maximally connected NLCE in the presence of strong next-nearest neighbor couplings. (This is apparent if one compares the clusters involved in both NLCEs for vanishing nearest neighbor couplings and nonvanishing next-nearest neighbor couplings.) However, this comes at a price as the number of clusters in the site-based NLCE grows exponentially with the order of the expansion. For a given number of sites $L$, there are $2^{L-1}$ clusters. Identifying which clusters are not related by symmetries (symmetry distinct clusters), allows one to reduce the number of clusters that need to be diagonalized by about a factor 2, as shown in Table 1. The top two clusters with three sites in Fig. 3 are an example of clusters that are related by a symmetry (a reflection about the center of the next-nearest neighbor bond).

(iii) Bond-based NLCE: Starting from a single site, a bond with a nearest or next-nearest neighbor site is added one at a time. This NLCE includes all possible clusters that can be drawn on the lattice. The order of the NLCE is set by the number of bonds of the largest clusters, and the number of clusters that have a given number of bonds $N$ is $\sim 3^{N}$. The clusters in this NLCE are the ones that appear in the traditional high-temperature expansion. However, in the bond-based NLCE a full exact diagonalization calculation is carried out for each cluster, instead of the expansion in powers of $\beta$ carried out in the high-temperature expansion. A drawback of this expansion, when compared to the previous two, is that it is computationally more expensive as the number of clusters increases more rapidly with the number of bonds than for the other two with the number of sites. In addition, its convergence is generally worse than that of the site-based NLCE, as discussed in Ref. Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2 for various two-dimensional lattice geometries. Because of this, we do not consider the bond-based NLCE any further here.

III.1 Resummations

To accelerate the convergence of NLCEs, we use Wynn’s and Euler algorithms Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2.

First, we group together the contributions of all clusters with the same number of sites $L$

[TABLE]

Next, we make explicit that the sum in Eq. (2) can be carried out in clusters with up to $l$ sites, which gives the prediction of the $l^{\text{th}}$ order of the NLCE

[TABLE]

The goal of resummation algorithms is to predict the result for $\mathcal{O}_{l\rightarrow\infty}$ from a finite sequence $\{\mathcal{O}_{l}\}$.

Wynn’s ($\epsilon$) algorithm is generally observed to give the best results for thermal states Rigol et al. (2006); *rigol2007numerical1; *rigol2007numerical2. In this resummation algorithm, $\epsilon_{n}^{(k)}$ is defined as

[TABLE]

Here $k$ is the number of Wynn cycles. Only even cycles are expected to converge to the $l\rightarrow\infty$ result. Note that, every two cycles, the new sequence generated is reduced by two terms. We denote the estimate after $2k$ cycles as

[TABLE]

We also find the Euler transformation to be useful, especially when studying quenches using the site-based NLCE. This resummation accelerates the convergence when the series $\{S_{L}\}$ alternates in sign. When the Euler resummation is used in the last $k$ terms of the series $\{S_{L}\}$, we write

[TABLE]

III.2 Observables and convergence of NLCEs

We focus on a set of observables that includes: (1) The total energy per site $E=\text{Tr}[\hat{H}\hat{\rho}]/N$, (2) the variance of the energy per site $\Delta E^{2}=(\text{Tr}[\hat{H}^{2}\hat{\rho}]-(\text{Tr}[\hat{H}\hat{\rho}])^{2})/N$, (3) the entropy per site $S=-\text{Tr}[\hat{\rho}\ln(\hat{\rho})]/N$, (4) the total interaction energy per site $U=\text{Tr}[\hat{U}\hat{\rho}]/N$, where

[TABLE]

and (5) the occupation of the zero momentum mode, $m_{k=0}=\text{Tr}[\hat{m}_{k=0}\hat{\rho}]$. The momentum distribution $\hat{m}_{k}$ is the discrete Fourier transform of the one-body density matrix

[TABLE]

$m_{k}$ is of particular interest as it can be measured in experiments with ultracold atoms using time-of-flight expansion Bloch et al. (2008). Note that while $U$ is a local observable $m_{k}$ is a nonlocal one. In the expressions above, $\hat{\rho}$ stands for the many-body density matrix.

The convergence of the NLCE to the thermodynamic limit result is probed, for an observable $\mathcal{O}$, by calculating the normalized difference

[TABLE]

where $l_{\text{max}}$ is the highest order of the NLCE that we are able to calculate. When $\Delta_{l}(\mathcal{O})$ reaches machine precision, one can think of the $l^{\text{th}}$ order of the NLCE as converged to the thermodynamic limit result.

The site-based NLCE has a higher computational cost and a higher computational error than the maximally connected one due to the exponentially large number of clusters that need to be considered. Because of its computational cost, the site-based NLCE is carried out up to the $15^{\text{th}}$ order ($\text{NLCE-S}_{15}$), while the maximally connected one is carried out up to the $18^{\text{th}}$ order ($\text{NLCE-M}_{18}$).

We also report results for exact diagonalization calculations in chains with periodic boundary conditions (ED-PBC). They are carried out in lattices with up to 20 sites ($\text{ED-PBC}_{20}$). In the ED-PBC calculations, we only consider even number of sites as those chains exhibit smaller finite-size effects. We compute the equivalent of Eq. (11) as a function of the chain size $L$

[TABLE]

IV Thermal equilibrium:

Grand canonical ensemble

In this section, we discuss the performance of the NLCEs for 1D systems in thermal equilibrium (in the grand canonical ensemble) described by Hamiltonian (1). We benchmark them against full exact diagonalization calculations in chains with periodic boundary conditions.

In the grand canonical ensemble, the density matrix of cluster $c$ is given by the expression

[TABLE]

where $\hat{H}_{c}$ and $\hat{N_{c}}=\sum_{i\in c}\hat{n}_{i}$ are the Hamiltonian and the total particle number operator in cluster $c$, respectively. $\hat{H}_{c}$ is given by Eq. (1) with all the nearest $(t,V)$ and next-nearest $(t^{\prime},V^{\prime})$ neighbor couplings present in the cluster. We set $t=V=1$, and focus on three sets of parameters $t^{\prime}=V^{\prime}=0.5,$ 1, and 2, which help understand the effect of increasing next-nearest neighbor couplings in the NLCE computations. We take $\mu=0$, which results in the systems being at half-filling as Hamiltonian (1) is particle-hole symmetric.

In Figs. 4(a) and 4(b), we plot the energy $E$ and the variance of the energy $\Delta E^{2}$ vs temperature as obtained using NLCEs and exact diagonalization for $t^{\prime}=V^{\prime}=0.5,$ 1, and 2. For each value of $t^{\prime}=V^{\prime}$, the plots are indistinguishable from each other for $T\gtrsim 1$. Below $T\approx 1$, the site-based NLCE exhibits a sharp divergence at a temperature that increases with increasing the value of $t^{\prime}=V^{\prime}$. The maximally connected NLCE and the exact diagonalization results differ from each other only for the energy $E$ at the lowest temperatures.

The normalized differences [Eq. (11)] for $E$ and $\Delta E^{2}$ vs temperature, between the last two orders of the NLCEs and between the two largest chains in the exact diagonalization, are reported in Figs. 4(c) and 4(d) for $t^{\prime}=V^{\prime}=0.5$ and 2. For each value of $t^{\prime}=V^{\prime}$, one can see that the exact diagonalization results are the first to depart from machine precision as the temperature is decreased. The maximally connected and site-based NLCEs results depart from machine precision about the same temperature for $t^{\prime}=V^{\prime}=2$, while for $t^{\prime}=V^{\prime}=0.5$ the site-based NLCE departs from machine precision at higher temperature than the maximally connected one. The results for $t^{\prime}=V^{\prime}=1$ (not shown) are in between.

The results obtained for other observables are qualitatively similar to the ones for the energy and the variance of the energy. In Fig. 5, we report results for the entropy $S$, the total interaction energy $U$, and the occupation of the zero momentum mode $m_{0}$. The main difference between those results and the ones for the energy and the variance of the energy is seen for $m_{0}$. Figure 5(f) shows that the normalized differences of the site-based NLCE are the ones that, when increasing $T$, reach machine precision at the lowest temperature. Namely, the site-based NLCE is the one that converges at the lowest temperatures.

Next one can ask, for a given temperature at which the NLCEs converge to machine precision, how the converged result is obtained with increasing the order of the expansion. In Fig. 6(a)–6(e), we plot the normalized distances for the five observables studied in Figs. 4 and 5 as a function of the order $l$ of the NLCEs and of the chain size $L$ of the exact diagonalization. We report results for $t^{\prime}=V^{\prime}=0.5$ and 2, and for $T=5.96$, for which the NLCEs have converged to machine precision for all observables studied. Figures 6(a)–6(e) show that: (1) the convergence of NLCEs (exact diagonalization) with increasing the order of the expansion (the chain size) is exponential Iyer et al. (2015), (2) for $E$, $\Delta E^{2}$, and $S$, the convergence of the maximally connected and the site-based NLCEs is similar, and they converge more rapidly than the exact diagonalization calculations, (3) for $U$ and $m_{0}$, the site-based NLCE exhibits the fastest convergence, followed by the maximally connected NLCE and exact diagonalization calculations, and (4) in all calculations, increasing the strength of the next-nearest neighbor couplings slows down the convergence.

Qualitative differences between the site-based NLCE, the maximally connected NLCE, and exact diagonalization mostly occur at low temperatures, at which the site-based NLCE diverges. In Fig. 6(f), we plot the normalized differences for the zero momentum mode occupation at $T=0.1$. [$\Delta_{l}(m_{0})$ in the site-based NLCE was computed using the result for the $18^{\text{th}}$ order of the maximally connected NLCE.] The normalized differences decrease for the maximally connected NLCE (with $l$) and the exact diagonalization (with $L$) calculations, but increase nearly exponentially with $l$ for the site-based NLCE. We discuss the origin of this divergence in Sec. VI.

Whenever NLCE calculations fail to converge, one can use resummations to accelerate the convergence. In Fig. 7(a)–7(c), we show $S$, $m_{0}$, and $E$ for $t^{\prime}=V^{\prime}=0.5$, 1 and 2, respectively, after various Wynn resummation cycles (see Sec. III.1) applied to both the site-based and the maximally connected NLCEs. Wynn’s algorithm accelerates the convergence the most for systems in thermal equilibrium. It is remarkable that the resummed results of the two NLCEs agree with each other for all observables and most temperatures reported in Fig. 7. They also agree with the maximally connected NLCE results below temperatures at which the site-based NLCE exhibits a divergence, and are clearly different from the exact diagonalization ones, which suffer from finite-size effects. Having an agreement between resummed results of two different NLCEs suggests that the resummed results are converged and that resummations allow one to extend the regime of applicability of the NLCEs beyond the one provided by the bare NLCE sums in Eq. (2).

V Quantum Quenches

In this section, we study quantum quenches starting from initial thermal equilibrium states. Namely, the initial state is in thermal equilibrium for an initial Hamiltonian $\hat{H}^{I}$, at temperature $T_{I}$ [$\beta_{I}=(T_{I})^{-1}$] and chemical potential $\mu_{I}$. This means that, in each cluster $c$ of the NLCEs, the initial density matrix in the grand canonical ensemble is

[TABLE]

As in Sec. IV, $\hat{H}^{I}_{c}$ and $\hat{N}_{c}$ stand for the Hamiltonian and total particle number operator in cluster $c$, respectively.

The quench consists of changing the initial Hamiltonian $\hat{H}^{I}$ ($\hat{H}^{I}_{c}$ for cluster $c$) into a new time-independent (final) Hamiltonian $\hat{H}$ ($\hat{H}_{c}$ for cluster $c$) and at the same time disconnecting the system from the bath. As a result, the dynamics that follows the quench is unitary, and, as a function of time $\tau$, the density matrix of cluster $c$ can be written as (we set $\hbar=1$)

[TABLE]

where $|m\rangle$ and $E_{m}$ are the energy eigenkets and eigenvalues for $\hat{H}_{c}$, respectively. The expectation value of an observable at time $\tau$ in cluster $c$ is given by $\mathcal{O}_{c}(\tau)=\text{Tr}[\hat{\rho}_{c}(\tau)\hat{O}]$. Here, we are interested in the infinite-time average $\bar{\mathcal{O}}_{c}=\lim_{\tau^{\prime}\rightarrow\infty}(1/\tau^{\prime})\int_{0}^{\tau^{\prime}}O_{c}(\tau)d\tau$, which describes observables after relaxation D’Alessio et al. (2016). One can also write $\bar{\mathcal{O}}_{c}=\text{Tr}[\hat{\rho}^{\text{DE}}_{c}\hat{O}]$, where $\hat{\rho}^{\text{DE}}_{c}$ is the diagonal ensemble Rigol et al. (2008) density matrix in cluster $c$

[TABLE]

assuming that there are no degeneracies in the energy spectrum. This is the density matrix used in the NLCEs to obtain observables after relaxation in the thermodynamic limit following a quantum quench. The entropy in the diagonal ensemble Rezek and Kosloff (2006); Polkovnikov (2011); Santos et al. (2011), also known as the diagonal entropy, is computed as $S^{\text{DE}}_{c}=-\text{Tr}[\hat{\rho}^{\text{DE}}_{c}\ln\hat{\rho}^{\text{DE}}_{c}]$.

For nonintegrable systems, which are the ones of interest in this work, it has been shown that the NLCE results for the diagonal ensemble agree with those of the grand canonical ensemble with the temperature $T$ and chemical potential $\mu$ [see Eq. (13)] selected such that Rigol (2014a, 2016)

[TABLE]

which means that the observables studied thermalize.

In what follows, we restrict our analysis to half-filled systems, which is achieved by setting $\mu_{I}=\mu=0$ (the total number of particles per site is conserved during the quench). Also, the initial Hamiltonian $\hat{H}^{I}$ is always taken to have only nearest neighbor couplings $t_{I}=0.5$ and $V_{I}=1.5$ ($t^{\prime}_{I}=V^{\prime}_{I}=0$). After the quench, the Hamiltonian parameters change to $t=V=1$, and $t^{\prime}=V^{\prime}=0.5$, 1, and 2 (the same parameters considered for systems in thermal equilibrium in Sec. IV).

The temperature of the grand canonical ensemble that describes observables after relaxation following the quench (in short, the “final effective temperature”) $T$ is calculated using the maximally connected NLCE requiring that the relative difference between the total energy per site predicted within the diagonal and grand canonical ensembles for $l=18$ is below $10^{-12}$. In Fig. 8, we plot $T$ vs $T_{I}$ for the quenches of interest in this work. For large values of $T_{I}$, the dependence of $T$ on $T_{I}$ is an almost linear one (see the fits in Fig. 8).

Quantities such as the energy $E$, the variance of the energy $\Delta E^{2}$, and the variance of the total particle number $\Delta N^{2}$ are conserved during the quench, e.g., $\text{Tr}[\hat{\rho}^{I}\hat{H}]=\text{Tr}[\hat{\rho}^{\text{DE}}\hat{H}]$ (note that this is different from $\text{Tr}[\hat{\rho}^{I}\hat{H}^{I}]$). This means that they converge exponentially fast with the order of the NLCEs, as discussed in Sec. IV for various observables in systems in thermal equilibrium and as seen in Refs. Rigol (2014a, 2016) for quantum quenches. Also, since there are no next-nearest neighbor couplings in $\hat{H}_{I}$, the site-based and maximally connected NLCEs provide identical results within machine precision (the clusters with nonzero weight in both expansions are the same).

In Fig. 9(a) and 9(b), we show $E$ and $\Delta E^{2}$ calculated using the maximally connected NLCE in the diagonal ensemble vs $T_{I}$. The energies after the quench are rather close to each other for the values of $t^{\prime}=V^{\prime}$ considered here, but the variance of the energy increases significantly as the value of $t^{\prime}=V^{\prime}$ is increased. Results for the normalized differences between the last two orders of the maximally connected and site-based NLCEs are presented in the main panels of Fig. 9(c) and 9(d) for $E$ and $\Delta E^{2}$, respectively. They can be seen to be within machine precision for the temperatures considered. The insets in Fig. 9(c) and 9(d) also show normalized differences but plotted as a function of the order of the NLCE for a fixed temperature $T_{I}=2$. As expected, the site-based and the maximally connected NLCEs give identical results at low orders. However, as the order is increased, the site-based NLCE saturates at a higher normalized difference (machine precision for this expansion) because it involves exponentially many more clusters than the maximally connected one.

To determine whether observables in the diagonal ensemble approach the grand canonical ensemble predictions as the order of the NLCE is increased, we compute the normalized difference

[TABLE]

The value of this quantity for the energy must be within machine precision in the $18^{\text{th}}$ order of the maximally connected NLCE, as the temperature is determined so that $\delta_{l=18}(E)<10^{-12}$. The inset in Fig. 9(a) shows that this is indeed the case in our calculations. On the other hand, the variance of the energy after the quench does not generally agree with that of the corresponding thermal equilibrium ensemble. This can be seen in the main panel in Fig. 9(b), in which we also plot the grand canonical ensemble predictions for $\Delta E^{2}$, and in the inset in Fig. 9(b), in which we plot $\delta_{l}(\Delta E^{2})$ vs $l$ for $T_{I}=2$. This is the result of $\Delta E^{2}$ being a conserved quantity during the quench. As discussed in Ref. Rigol (2014a), $\Delta E^{2}$ can be used to distinguish the diagonal ensemble from thermal ensembles. We note that an agreement between $\Delta E^{2}$ in the diagonal and thermal ensembles is not required for few-body observables to thermalize. So long as the variance of the energy is extensive (as it is in all the quenches considered here), such that its square root (the width of the energy distribution) is subextensive, and eigenstate thermalization occurs for a given observable, then the observable is guarantied to thermalize in the thermodynamic limit D’Alessio et al. (2016).

For observables that are not conserved during the quench, the overwhelming majority of observables of interest, the convergence of NLCEs is slower for the diagonal ensemble than for the grand canonical one Rigol (2014a). In what follows, we focus on the behavior of the entropy $S$ (a thermodynamic observable), the potential energy $U$ (a local observable), and the occupation of the zero momentum mode $m_{0}$ (a nonlocal observable). Since we restrict our analysis to quenches in which the final Hamiltonian has $t^{\prime}=V^{\prime}\neq 0$, namely, a nonintegrable Hamiltonian, those observables are expected to thermalize D’Alessio et al. (2016). Hence, we can use converged results from the corresponding grand canonical ensemble to study the convergence of the NLCEs for the diagonal ensemble. This means that the normalized difference that will be the focus of the analysis that follows is the one given by Eq. (19). In all calculations of $\delta_{l}(\mathcal{O})$ in the maximally connected and site-based NLCEs, as well as of

[TABLE]

in exact diagonalization, we use $\mathcal{O}_{18}^{\text{GE}}$ obtained from the maximally connected NLCE. We restrict our analysis to initial temperatures such that, for the grand canonical ensemble that describes the systems after the quench, $\Delta_{17}(S)\lesssim 10^{-10}$, $\Delta_{17}(m_{0})\lesssim 10^{-10}$, and $\Delta_{17}(U)\lesssim 10^{-8}$.

In Figs. 10(a)–10(c), we show $S$, $U$ and $m_{0}$ in the diagonal ensemble as obtained using the $18^{\text{th}}$ order of the maximally connected NLCE and exact diagonalization in a periodic chain with 20 sites (the results of the site-based NLCE for the diagonal ensemble fall beyond the scale of these plots). We also report the results obtained for the same observables in the grand canonical ensemble using the $18^{\text{th}}$ order of the maximally connected NLCE. The agreement between the diagonal and grand canonical ensemble results within the maximally connected NLCE is excellent. Differences between those two and exact diagonalization results are only apparent for $S$ and $m_{0}$ at the lowest temperatures, and for $U$ at most temperatures.

In Figs. 10(d)–10(f), we show the normalized differences $\delta_{l}(S)$, $\delta_{l}(U)$ and $\delta_{l}(m_{0})$ for the last order of the NLCE calculations and $\delta_{L}(S)$, $\delta_{L}(U)$ and $\delta_{L}(m_{0})$ for the largest chain in the exact diagonalization calculations. Those differences can be seen to decrease with increasing $T_{I}$ for $S$ and $m_{0}$, as expected as increasing $T_{I}$ should decrease the convergence (NLCE) and finite-size (exact diagonalization) errors. The differences remain fairly constant for $U$ because $U_{18}^{\text{GE}}$ [which is in the denominator of $\delta_{l}(U)$] approaches zero with increasing $T_{I}$. The other fact that is apparent in Fig. 10(d)–10(f) is that the maximally connected NLCE exhibits normalized differences that are between one and two orders of magnitude smaller than those of the exact diagonalization, and that the site-based NLCE exhibits differences that are several orders of magnitude larger at all initial temperatures. The latter is qualitatively different from what was discussed for systems in thermal equilibrium in the context of Fig. 5.

Figures 10(g)–10(i) depict the same normalized differences but as a function of $l$ (NLCEs) and $L$ (exact diagonalization), at a fixed initial temperature $T_{I}=10$. Those plots make apparent that while the maximally connected NLCE and exact diagonalization calculations converge (though more slowly than for systems in thermal equilibrium) with increasing the order of the expansion and the chain sizes, respectively, the site-based NLCE diverges even at this relatively high initial temperature. We discuss the origin of this divergence in the next section.

As for systems in thermal equilibrium, one can use resummations to accelerate the convergence of NLCEs for the diagonal ensemble. In Fig. 11, we show results for resummations of the site-based NLCE for the entropy [Fig. 11(a)] and the occupation of the zero momentum mode [Fig. 11(b)]. The results after resummation can be seen to be close to the thermal equilibrium predictions, and in some instances [such as the $6^{\text{th}}$ cycle of the Wynn algorithm for $m_{0}$ in Fig. 11(b)] they are actually closer to the thermal equilibrium predictions than the exact diagonalization ones. However, the direct sums for the maximally connected NLCE exhibit a better agreement with the thermal equilibrium predictions than the resummed results for the site-based NLCE. In summary, our study indicates that the maximally connected NLCE is the best approach (among the ones considered in this work) for diagonal ensemble calculations after quantum quenches in the thermodynamic limit.

VI Divergence of the site-based NLCE

In the systems in thermal equilibrium discussed in Sec. IV, we have seen that below some temperature the site-based NLCE diverges, while for quantum quenches it appears to diverge for all initial temperatures considered. Here, we provide an understanding of the origin of those divergences, and their differences, by studying the average weight of the clusters with a given number of sites, and comparing it to the number of clusters with that number of sites.

The average cluster weight $W(\mathcal{O},L)$ [see Eq. (2)], of an observable $\mathcal{O}$ in clusters with $L$ sites, is

[TABLE]

When the clusters in the NLCE are not large enough compared to the correlation length associated with a given observable, the weights of the observable in those clusters need not be small even if they decrease with increasing the cluster size. Since the number of clusters in the site-based NLCE grows exponentially fast with the number of sites in the clusters, the combination of non-rapidly-enough decreasing weights with the exponential increase of the number of clusters can lead to a divergence. For the maximally connected NLCE, there is only one cluster with any given number of sites. Hence, no divergence is expected with increasing cluster size no matter the cluster sizes considered.

In Fig. 12, we plot the absolute value of the average weights in the grand canonical ensemble (at various temperatures) for two observables ($S$ and $m_{0}$), as well as the inverse of the number of clusters $\text{nc}^{-1}$, as a function of the cluster sizes $L$. The site-based NLCE can converge only if the average weight decreases more rapidly than $\text{nc}^{-1}$. In Figs. 12(a) and 12(c), one can see that this happens for all $L$ only at the highest temperature shown. At the second highest temperature, the average weights decrease more slowly than $\text{nc}^{-1}$ for small $L$ and then more rapidly for large $L$ (hence the importance of calculating the NLCE up to largest possible cluster sizes). For the third highest temperature, the average weights decrease more slowly than $\text{nc}^{-1}$. For $T=0.1$, the lowest temperature shown, the average weights decrease only for the largest clusters for $S$ and increase for $m_{0}$.

For the maximally connected NLCE [Figs. 12(b) and 12(d)], $\text{nc}^{-1}=1$ and the average weights can be seen to decrease exponentially fast with increasing $L$ for all but the lowest temperature (ensuring the convergence of this NLCE). We note that the average weight of $m_{0}$ in the site-based NLCE decreases much more rapidly than in the maximally connected one at high temperatures. This helps understand why the former NLCE could converge faster than the latter one in Fig. 6(c), and why the site-based NLCE can be used to complement the maximally connected one to study 1D lattice systems in thermal equilibrium.

The behavior of the average weights of observables in the diagonal ensemble after quantum quenches [Fig. 13] is fundamentally different from that observed for systems in thermal equilibrium. Both, in the site-based and maximally connected NLCEs, the rate at which the weights decrease with $L$ does not appear to depend significantly on $T_{I}$, or the effective final $T$ (see Fig. 8 for the relation between the two), as opposed to the behavior seen in Fig. 12. All the effect that those temperatures seem to have is to offset the decay of the weights at large $L$. In addition, while the weights in the maximally connected NLCE [see Figs. 13(b), 13(d), and 13(f)] exhibit a decrease that is consistent with an exponential for large $L$ [except for the lowest temperatures in Fig. 13(f)], the weights in the site-based NLCE exhibit an erratic behavior [see Figs. 13(a), 13(c), and 13(e)] and fail to decrease faster than $\text{nc}^{-1}$. This explains why the site-based NLCE diverges with increasing the order of the expansion for all temperatures and cluster sizes considered. What determines the rate of decay of the weights in the diagonal ensemble after quantum quenches, and why their behavior is fundamentally different from the one in the grand canonical ensemble, is a topic that we will explore in future studies.

VII Summary

In summary, we presented a comprehensive analysis of the performance of maximally connected and site-based NLCEs when applied to the study of hard-core boson systems in 1D lattices with nearest and next-nearest couplings in thermal equilibrium and after quantum quenches. We compared the results of the NLCE calculations with those of exact diagonalization in chains with periodic boundary conditions. Overall, we find that the maximally connected NLCE is the most efficient and provides the most accurate results in most regimes and for most observables studied. Only for some observables, such as the occupation of the zero momentum mode, in systems in thermal equilibrium at high temperatures, did we find that the site-based NLCE exhibits a faster convergence than the maximally connected one. We showed that when the NLCEs fail to converge, resummations can be used to accelerate convergence and extend the regime of applicability of NLCEs to lower temperatures.

We showed that, for observables that are not conserved during the quench, the convergence of the maximally connected NLCE with increasing cluster sizes is slower in the diagonal ensemble than in the grand canonical one. On the other hand, the site-based expansion was found to diverge with increasing cluster size for all initial temperatures and quenches considered. We argued that this is a result of the qualitatively different behavior between the weights of observables in the diagonal and the grand canonical ensembles (in the clusters considered here). Understanding the origin of those differences, as well as what determines the rate of decay of the weights in the diagonal ensemble, is left for future work.

Acknowledgements.

This work was supported by the U.S. Office of Naval Research, Grant No. N00014-14-1-0540. The computations were done at the Institute for CyberScience at Penn State.