Thermal bosons in 3d optical lattices via tensor networks

Ultracold atoms in optical lattices are one of the most promising experimental setups to simulate strongly correlated systems. However, efficient numerical algorithms able to benchmark experiments at low-temperatures in interesting 3d lattices are lacking. To this aim, here we introduce an efficient tensor network algorithm to accurately simulate thermal states of local Hamiltonians in any infinite lattice, and in any dimension. We apply the method to simulate thermal bosons in optical lattices. In particular, we study the physics of the (soft-core and hard-core) Bose–Hubbard model on the infinite pyrochlore and cubic lattices with unprecedented accuracy. Our technique is therefore an ideal tool to benchmark realistic and interesting optical-lattice experiments.


The TgPEPS method
Traditionally, TN methods for zero-temperature properties target ground states of local Hamiltonians by, e.g., variational optimization 41,72 and imaginary-time evolution (ITE) [73][74][75] . However, at finite-temperature we target the thermal density matrix (TDM) of Hamiltonian H, i.e., ρ = e −βH , β = 1/T being the inverse temperature. To approximate this state, one typically evolves in imaginary-time for a time β/2 both the bra and ket degrees of freedom starting from the infinite-temperature state, i.e., ρ = e −βH/2 · I · e −βH/241, 63,64,66,69 . In TN language, the TDM of a system on a lattice with coordination number z can be described by a Projected Entangled-Pair Operator (PEPO) 36 . The typical tensor describing this PEPO is of the type T i,j,α 1 ,··· ,α z , with i, j = 1, . . . , p and α k = 1, . . . , D , p being the dimension of the local Hilbert space and D the bond dimension of the bond indices, Scientific Reports | (2020) 10:19051 | https://doi.org/10.1038/s41598-020-75548-x www.nature.com/scientificreports/ which in turn controls the amount of correlations (classical and quantum) that can be handled by the ansatz, see Fig. 1(b) for the example of the cubic lattice. Let us consider, without loss of generality, the case of a Hamiltonian with nearest-neighbour interactions, H = �i,j� h ij . The TDM can be sliced into m ITE steps using the Suzuki-Trotter decomposition, i.e, where δ ≪ 1 and m · δ = β . The thermal density matrix of the system at inverse temperature β is then obtained after m successive applications of the gates g ≡ e −δh ij /2 on the corresponding links of the lattice, as shown in Fig. 1(b,c) for the cubic lattice.
After applying a gate on the PEPO, the bond dimension of the index connecting the local sites grows from D to p 2 D . To truncate it back to its original size one can use a variety of methods, including the so-called "simple" update (SU) 61,76,77 , or the more involved "fast-full" 78 and "full" updates (FFU, FU) 52,55,57 . These algorithms differ in the way the handle the correlations around the link to be truncated. While the SU handles these correlations by a mean-field-like approach and truncates the index directly with a Singular Value Decomposition (SVD), FFU and FU approximate the full effect of the environment via powerful techniques such as Tensor Renormalization Group (TRG) 79,80 and Corner Transfer Matrix Renormalization Group (CTMRG) 52,77,[81][82][83] . The price to pay, though, is that FFU and FU are computationally much more expensive than SU. Besides, application of the FU and FFU algorithms to lattices with higher dimensionality and connectivity is highly non-trivial. However, for thermal systems with a lot of connectivity (such as higher-dimensional systems), the mean-field-like approximation of the environment in the SU is actually good, in turn making the SU a quite accurate option in these situations. And this are good news, because the computational cost of the FFU and FU algorithms is really high for thermal high-dimensional systems. In this tensor network diagram, shapes are tensors, lines are indices, and connected lines are contracted common indices. Red and blue indices at every tensor correspond to the local bra and ket degrees of freedom, respectively. (c) Action of the Suzuki-Trotter gate g on both ket and bra indices of two nearest-neighbouring sites, for the PEPO tensors of the cubic lattice. Small red tensors correspond to the matrices of singular values obtained from the simple update 76 . (d) Thermal expectation value of a local two-body operator O in the cubic lattice, with a mean-field approximation of the environment (see main text for more details).
Using the "Choi's" isomorphism 69,70 , the density operator ρ(β/2) can be replaced by its vectorized form |ρ(β/2)� in the product space of the ket and bra. In this vectorized representation, the partition function then becomes an inner product of |ρ(β/2)� and its vector conjugate, namely and the un-normalized expectation value of local operators reads The tensor network representation of Eq. (4) for a two-body operator is demonstrated in Fig. 1(d) using the mean-field approximation to the environment which, as we argued before, works well in the considered regimes (notice also that, unlike at zero temperature, here we do not target variational approximations of ground-state energies). In the end, our approach allows us to push the simulation of thermal 3d models in any geometry to very large dimensions with a very cheap computational cost of O(p 2 D z ).

Numerical results
To prove the validity of our approach, we apply the TgPEPS technique to the Bose-Hubbard (BH) model 18 in 3d, for which we study the low-temperature phase diagram in the cubic and pyrochlore lattices, with maximum occupation number n oc = 2 (soft-core) and n oc = 1 (hard-core) respectively. The pyrochlore lattice is one of the most challenging structures from the point of view of TN simulations. Besides, the hard-and soft-core cases have different physical dimensions, i.e., p = 2 and p = 3 , respectively resulting in different computational complexity. Our choice of these lattices and occupation numbers was to seriously challenge our algorithm and assess the stability and accuracy of the simulations for different physical systems and lattice geometries. The generic Hamiltonian of the Bose-Hubbard model is given by where a † (a) are the bosonic creation (annihilation) operators, n = a † a is the particle number operator, t is the hopping rate between nearest-neighbour sites, U is the on-site repulsive interaction, and µ the chemical potential. At zero-temperature T = 0 , in the extreme regime where t/U ≪ 1 the repulsive interaction is very strong and only one particle per site is allowed, resulting in a Mott insulating (MI) phase. Complementarily, in the t/U ≫ 1 regime the bosons are highly delocalized, and the system is in a coherent superfluid (SF) phase. One therefore expects a quantum phase transition (QPT) between these two regimes. Let us first consider the T = 0 properties for the soft-core ( n oc = 2 ) BH model on the cubic lattice, which we can compute accurately using our technique from Ref. 61 . Figure 2(a) shows the T = 0 SF-MI transition at t/U = 0.01 . This is captured by the superfluid order parameter ρ s = |�a i �| 2 which is nonzero in the SF phase and zero in the MI phase. The QPT takes place at (µ/U) c ≈ 1.119 . The discontinuity in the second derivative of the ground state energy per site (inset of Fig. 2(a)) confirms that this transition is second-order. We have further mapped out the full T = 0 phase diagram in the t/U-µ/U plane up to n oc = 2 in Fig. 3(a). The figure depicts the contour-plot of the superfluid density, which is zero in the two Mott lobes (dark regions).
It is already known that there exist two different types of transitions between the Mott insulator and superfluid regions of the Bose-Hubbard phase diagram. The trivial phase transition occurs when the MI-SF phase boundary is crossed at fixed t/U and a critical phase transition at fixed integer density when the phase boundary www.nature.com/scientificreports/ is crossed at fixed µ/U such as the tip of the Mott lobes. While the former transition can be described by the physics of a weakly interacting Bose gas and trivial (mean-field) universality class, the latter transition which occurs at relatively large t/U is described by highly fluctuating delocalized bosons with emerging particle-hole symmetry and relativistic dispersion 43,44,86 that belongs to the (d+1)-XY universality class 18,43 . Our accurate analysis shows that the tip of the first lobe (blue star in Fig. 3(a)) is located at µ/U ≈ 0.393 , in excellent agreement with previous studies 43,44,86 . We find that the QPT occurs at (t/U) c ≈ 0.0331 , which is well detected by ρ s as well as by the one-site Von-Neumann entanglement entropy shown in Fig. 2(b). On the 3d cubic lattice, the QPT at the tip of the lobe is in the four-dimensional XY universality class 43,86 which for d > 3 is equivalent to mean-field universality class hence, is trivial. Next, we study the thermal properties of the BH model. In experiments with ultracold gases in optical lattices, bosonic atoms are cooled down to the nanoKelvin regime by efficient techniques such as laser cooling 87,88 and the desired quantum states are engineered by inducing quantum correlations between the atoms by microwave pulses. Increasing the temperature of the system destroys quantum correlations due to the extra kinetic energy imposed on the atoms by thermal fluctuations, and eventually the system will end up in a normal gas (NG) phase at large T 42 . A thermal phase transition (TPT) is therefore expected between the underlying quantum state and the NG phase, depending both on T as well as the couplings of the Hamiltonian.
We further use TgPEPS to provide a good insight into the stability of the Mott and superfluid phases at finite-T. For this, we first dive deeply into the MI and SF regions of the zero-T phase diagram computed previously, and then ramp up the temperature from zero up to large T. To this end, we fixed t/U = 0.01 and computed the TDM of the BH model in the cubic lattice for µ/U = 0.5 and µ/U = 1.0 , which respectively correspond to the MI and SF phases. Figure 2(c,d) show the energy ε and specific heat C = ∂ε/∂T of the BH model on the cubic lattice versus T in the MI and SF phases. The ideal MI phase only exist at T = 0 . However, at very low-temperature www.nature.com/scientificreports/ some Mott-like features still persist in the system but with a small finite compressibility. The boundary between the MI and NG phases is therefore a thermal crossover not a thermal phase transition. Our analysis shows that the MI-NG crossover occurs at T ′ ≈ 0.38 as can be captured by the peak of the specific heat in Fig. 2(c). On the other hand, the superfluid phase can exist even at finite-temperature and. The SF-NG phase boundary is therefore a true thermal phase transition which for the example couplings of our choice is located at T c ≈ 0.87 (see the peak of C v in Fig. 2(d)). Our TN thermometry analysis thus, indicates that both MI and SF states are stable only at very low temperatures which evaporate to a normal gas for T/t 1.
Additionally, we have mapped the finite-T phase diagram of the model in the t/U-µ/U plane. Figure 3(b-d) illustrate the superfluid density ρ s for T = 0.05, 0.1, 0.2 , revealing how the thermal fluctuations shrink the SF region. The distinction between the MI and NG phases (dark regions) is not visible in the ρ s plot because ρ s = 0 for both. However, one can distinguish them by observing the particle density ρ 0 = �a † i a i � , which is shown in Fig. 4(a): it holds integer values in the MI phase, and non-integer values in the NG phase.
Finally, we used TgPEPS to study the hard-core BH model on the pyrochlore lattice at finite temperature. This lattice is well-known for being responsible for interesting effects in frustrated quantum magnets and Kitaev materials. Figure 4(b,c) shows the superfluid density and the particle density of this model for various temperatures. The phase boundaries between SF and NG and between MI and NG can clearly be identified from the curves in both plots.
Let us note that in principle there is no limit on the choice of occupation number, n oc , in our TN simulation of the BH model. However, one should note that by increasing the n oc , the computational cost of the simulation is increased as well. Lets us further stress that there is no soft-core cut-off n oc in our algorithms for simulating the physics of each lobe and simply by fixing the proper physical dimension to p = n oc + 1 (including the unoccupied sites) for each lobe, one obtains the correct physics for all lobs with occupation n ≤ n oc .

Conclusions and outlook
In this paper we introduced TgPEPS, an efficient TN algorithm for finite-temperature simulation of strongly correlated systems in the thermodynamic limit for any lattice and in any spatial dimension. The method follows the ideas that we introduced in Ref. 61 , extending them to the finite-temperature case, and allowing for very efficient and accurate simulations of thermal systems with large connectivity. We benchmarked the method by computing the zero-and finite-temperature of the 3d Bose-Hubbard model in the pyrochlore and cubic lattices, with unprecedented accuracy. We believe that the TgPEPS algorithm can serve an essential tool for both benchmarking experiments with ultracold atoms in optical lattices with complex geometric structures. While the TgPEPS algorithm presented in this study has been developed for infinite systems, the optical lattice experiments Nevertheless, the TgPEPS algorithm can be adopted to finite systems with minor modifications. Let us further point out that while the method has been proven to be very successful for simulating local lattice Hamiltonians, some of the longer-range correlations which go beyond the nearest-neighbours may not be fully captured due to the simple mean-field-like environments that have been used for calculating the expectation values and correlators. Besides, after a detailed investigation of the TgPEPS algorithms for simulating some highly entangled systems such as the Kitaev spin liquids, we found out that the algorithms have difficulties in convergence at very low-temperature regime below T < 10 −3 due to the proximity of thermal states to the T = 0 ground state 89,90 .
Last but not least, for the lattice geometries that can be adopted to renormalization techniques, such as TRG and CTMRG, the thermal states obtained from TgPEPS method can be used to approximate the full environment from which one can calculate variational energies and expectation values. One can therefore increase the accuracy of the simulations by capturing longer-range correlations in the system. Finally, let us point out that the TgPEPS algorithm can also be applied to simulate fermionic atoms on optical lattices 91 , which we will address in the future extensions of this work. All in all, we think that our method will become a very helpful tool in the discovery of new exotic phases of quantum matter.

Appendix
Structure-matrix for the cubic and pyrochlore lattices. The structure-matrix (SM) codifies the connectivity information of the tensor network (TN) corresponding to a given lattice geometry. More specifically, the connectivity information of two neighbouring tensors along their shared edges are stored in the columns of the SM. The explicit SMs that we used in our calculations for TN simulation of the cubic (see Fig. 5) and     www.nature.com/scientificreports/ pyrochlore lattices (see Fig. 6) are provided in Table 1 and 2, respectively for translationally invariant unit-cells of 8-sites. For example, according to the second column of Table 1, the edge E 3 connects the bond matrix 3 and the dimensions 4 and 2 of tensors T 1 and T 4 , respectively. Thanks to this information, the algorithm can automatically identify the links and tensors where two-body gates are applied, and implement a simple update. This is done by looping over the columns of the SM and systematically updating the iPEPS tensors along their corresponding edges, which can now be done automatically and regardless of the underlying lattice. Last but not least, the non-zero elements of the SM at each row start from 2 which is due to the fact that the first two dimensions (0, 1) of the tensors T i in our notation corresponds to the physical bonds of ket and bra, respectively, and play no role in the connectivity of the underlying TN. The TgPEPS algorithm for any infinite lattice. The thermal density matrix of the system is evaluated by iteratively applying g on every shared link of the two neighbouring tensors T i , T j and updating the tensors along the corresponding links. In this scheme, the update changes only the tensors along the link where a given gate is acting. Therefore, one can update lower-rank sub-tensors related to them and substantially reduce the computational cost of the algorithm 78 , thus allowing to achieve larger bond dimension D.
Let us briefly revisit how the SU proceeds for the sub-tensors, in the context of TgPEPS. Given a tensor network and its corresponding structure matrix, the SU consists of the following iterative main steps: 1. Do for all edges E k , k ∈ [1, N Edge ] (columns of SM matrix) (a) Find tensors T i , T j and their corresponding dimensions connected along edge E k . (b) Absorb bond matrices m to all virtual legs m = k of T i , T j tensors. (c) Group all virtual legs m = k to form P l , P r PEPO tensors. (d) QR/LQ decompose P l , P r to obtain Q 1 ,R and L, Q 2 sub-tensors, respectively 78 . (e) Contract the ITE gate U i,j , with R, L along the legs of both ket and bra together with k to form tensor. (f) Obtain R , L , ˜ k tensors by applying an SVD to and truncating the tensors by keeping the D largest singular values. (g) Glue back the R , L , sub-tensors to Q 1 , Q 2 , respectively, to form updated PEPO tensors P ′ l , P ′ r . (h) Reshape back the P ′ l , P ′ r to the original rank-(z + 2) tensors T ′ i , T ′ j . (i) Remove bond matrices m from virtual legs m = k to obtain the updated PEPO tensors T i and T j . Figure 7-(a) shows all these steps graphically for the case of a 2d honeycomb lattice (the generalization to 3d is straightforward). By repeating the above ITE iterations m/2 times where m.δ = β , the system will cool down to the desired temperature T = 1/β . Once the TDM of the system at desired T is obtained, we can calculate the expectation values of local operators as shown graphically for one-and two-body operators in Fig. 7(b,c). In contrast to the gPEPS algorithm 61 where the optimization is performed for ground state tensors with one physical leg, in the TgPEPS the optimization is performed over thermal PEPO tensors with two physical degrees of freedom. (b) One-site and (c) two-site expectation values, as computed with a mean-field environment, in the TgPEPS scheme. All these diagrams are for a 2d honeycomb lattice, but they can be straifgtforwardly generalized to any 3d structure, such as the ones considered in the main text.
Scientific Reports | (2020) 10:19051 | https://doi.org/10.1038/s41598-020-75548-x www.nature.com/scientificreports/ In order to have an efficient and universal algorithm applicable to any infinite lattice, the following remarks are in order: (i) In steps (b), (c), (g) and (h) one can locate the lambda matrices corresponding to each leg of a tensor from rows of the SM. One can therefore design clever functions for absorbing (removing) matrices to (from) each tensor legs as well as for grouping (un-grouping) the non-updating tensor legs by using the information stored in each row of the SM. (ii) In our SU optimization, we perform the ITE iteration for δ = 10 −3 . (ii) Furthermore, one can increase the stability of the SU algorithm by applying the gauge-fixing 61 .
Let us further note that the computational cost of the SU scales as O(p 2 D z ) , and evidently depends on the coordination number of the underlying lattice. Henceforth, the maximum achievable bond dimension D is lattice dependent and is larger for structures with less coordination number, though structures with large z usually need low D because of entanglement monogamy.