Thermal DMRG for highly frustrated quantum spin chains: a user perspective

TL;DR

This paper evaluates the effectiveness of thermal DMRG in calculating magnetic properties of frustrated quantum spin chains at finite temperatures, highlighting its accuracy and limitations in different systems.

Contribution

It demonstrates the applicability of thermal DMRG to frustrated spin chains and analyzes its temperature-dependent accuracy in these complex systems.

Findings

Thermal DMRG accurately approximates magnetic observables in unfrustrated chains.

In frustrated sawtooth chains, accuracy is limited to higher temperatures.

Frustration affects the low-energy spectrum, impacting DMRG performance.

Abstract

Thermal DMRG is investigated with emphasis of employability in molecular magnetism studies. To this end magnetic observables at finite temperature are evaluated for two one-dimensional quantum spin systems: a Heisenberg chain with nearest-neighbor antiferromagnetic interaction and a frustrated sawtooth (delta) chain. It is found that thermal DMRG indeed accurately approximates magnetic observables for the chain as well as for the sawtooth chain, but in the latter case only for sufficiently high temperatures. We speculate that the reason is due to the peculiar structure of the low-energy spectrum of the sawtooth chain induced by frustration.

Equations

$$\underset{\sim }{\mathrm{\rho \rho }}(\beta ,B)$$

$$\tau =-i\beta \hslash /2,\beta =\frac{1}{k_{B}T}.$$

$$|{\Psi }_{\beta }⟩\equiv e^{-\beta \underset{\sim }{\mathrm{HH}}(B)/2}|{\Psi }_0⟩.$$

$$⟨\underset{\sim }{\mathrm{OO}}{⟩}_{\beta }=\frac{⟨{\Psi }_{\beta }|\underset{\sim }{\mathrm{OO}}|{\Psi }_{\beta }⟩}{⟨{\Psi }_{\beta }|{\Psi }_{\beta }⟩},\underset{\sim }{\mathrm{\rho \rho }}(\beta )=\frac{Z(0)}{Z(\beta )}{\text{Tr}}_{\text{Q}}|{\Psi }_{\beta }⟩⟨{\Psi }_{\beta }|,$$

$$\stackrel{\text{mbox}1}{\sim }=Z(0)\underset{\sim }{\mathrm{\rho \rho }}(0),Z(\beta )=Z(0)⟨{\Psi }_{\beta }|{\Psi }_{\beta }⟩.$$

$$\underset{\sim }{\mathrm{HH}}$$

$$\chi(T)=\dfrac{\partial M(T,B)}{\partial B}\bigg{|}_{B=0}\approx\dfrac{M(T,0.1)-M(T,0)}{0.1-0}=\dfrac{M(T,0.1)}{0.1}\ .$$

Full text

††journal: J. Magn. Magn. Mater.

1 Introduction

The evaluation of magnetic properties starting from model Hamiltonians allows one to obtain important information on magnetic materials such as magnetic molecules or low-dimensional magnets. A successful modeling of magnetic observables, e.g. magnetization, heat capacity or cross sections for inelastic neutron scattering, then leads to a concrete understanding of magnetic exchange patterns as well as anisotropic contributions to the Hamiltonian [1]. Nevertheless, a preferable complete and numerically exact diagonalization of the Hamiltonian is often spoiled by the vast dimension of the underlying Hilbert space that grows exponentially with system size. Even if group theoretical tools are employed to decompose the Hilbert space into smaller orthogonal subspaces according to the available symmetries and their irreducible representations, the problem of exponential growth cannot be resolved [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

Therefore, approximate methods provide a very valuable means to assess magnetic properties of larger spin systems. Some of these methods aim at only low-lying states, which is sufficient to understand low-temperature spectroscopic data, others address thermodynamic properties over a larger temperature range. Among the approximate methods that were successfully used over the past years are high-temperature series expansions (HTE) [14, 15, 16, 17, 18, 19, 20], density-matrix renormalization group theory (DMRG & DDMRG) [21, 22, 23, 24, 25, 26, 27], quantum Monte Carlo (QMC) [28, 29, 30, 31, 32, 33], as well as Lanczos and in particular finite-temperature Lanczos methods (FTLM) [34, 35, 36, 37, 38, 39, 40, 41, 42].

In the present article we would like to introduce another approximate method to the molecular magnetism community. It extends the applicability of DMRG to non-zero temperatures [43, 44, 45, 46, 47, 48, 49] and allows to treat very large spin systems, in particular larger than FTLM allows, without being hampered by the negative sign problem of QMC. The method – thermal DMRG (ThDMRG) – rests on an approximate time evolution of matrix product states along imaginary time from a state representing infinite temperature down to the temperature of interest. We will provide one example – a spin chain – where this approximation yields very accurate results and another one – a sawtooth chain– where the approximation fails for low temperatures.

The article starts with a short reminder of the thermal DMRG in Sec. 2, presents our results in Sec. 3, before it summarizes our main points in Sec. 4.

2 Thermal DMRG

The idea of finite-temperature, i.e. thermal DMRG (ThDMRG) is published in several references [43, 44, 45, 46, 47, 48, 49]. Here we would like to provide a very short description of the idea.

A density matrix at temperature $T$ and external field $B$ can be thought of as the time-evolved density matrix starting from a density matrix at infinite temperature down to the desired temperature $T$ by means of imaginary time evolution with the appropriate Hamiltonian $\!\hbox{}\mathop{\vphantom{H}\smash{H}}\limits_{\sim}\hbox{}\!(B)$ of the system

[TABLE]

The inverse temperature $\beta$ plays the role of the imaginary time. The density matrix at infinite temperature is proportional to the unit operator.

This approach is chosen since one can then employ time-dependent DMRG for the approximate time evolution of a quantum state describing the large spin system. Technically, the time evolution of the density matrix is replaced by the time evolution of a state vector which is first constructed by purifying the density matrix. Purification turns a density matrix into a ket state belonging to a larger Hilbert space in a way, that a trace over the additional (auxiliary) degrees of freedom returns the desired density matrix [50, 51, 52].

In particular one performs the following steps:

Set up the initial state $|\Psi_{0}\rangle$ that corresponds to infinite temperature. This state is a matrix product state in a Hilbert space consisting of the Hilbert space ${\mathcal{H}}_{P}$ of the physical system one wants to model as well as the Hilbert space ${\mathcal{H}}_{Q}$ of the unphysical system of auxiliary degrees of freedom called ancillas, compare Fig. 1. There is a certain freedom in constructing $|\Psi_{0}\rangle$ [53, 51]. In our implementation we employ the entangler Hamiltonian [51] to generate $|\Psi_{0}\rangle$.

• 2.

Evolve the state in imaginary time by means of time-dependent DMRG [54, 55, 56, 57]

[TABLE]

This is the numerically demanding part since the basis, which was optimized for the state at $\beta=0$, is not well-enough suited for later times and therefore must be adapted by targeting evolved states at certain times. For this purpose we project the time evolution operator onto a Krylov subspace [58]. Other time evolution algorithms such as Runge-Kutta or Chebyshev are equally feasible. A Suzuki-Trotter approach works best for nearest-neighbor interactions. Since additional long-range interactions are introduced by purifying the state, this method may not be optimal here. Grouping spins together to recover nearest-neighbor interaction is possible but quickly becomes tedious.

• 3.

Finally, calculate the thermal expectation value of a physical observable (not acting on the auxiliary degrees of freedom)

[TABLE]

with

[TABLE]

$\text{Tr}_{\text{Q}}$ denotes the partial trace in ${\mathcal{H}}_{Q}$.

• 4.

The accuracy of this method depends on a few parameters, among which the number $m$ of kept basis states plays a major role. We use $m$ values of up to 60, and work with fixed $m$. Alternatively, one could employ the discarded weight of DMRG as a measure of quality; this would lead to varying $m$ values in the course of calculations to keep the discarded weight constant. An important part of basis improvement are the so-called sweeps. We run one full sweep in-between time steps, which is usually enough to obtain a sufficiently adapted basis. One may choose to perform more sweeps to increase accuracy [59]. The time step $dt$ for the numerical time evolution is chosen as $dt=0.1$ (in natural units). Smaller $dt$ did not yield better results; on the contrary, one has to trade off between more accurate shorter time steps and the total time duration one can achieve.

In this outlined way one evaluates observables for decreasing temperatures (increasing $\beta$). The calculation has to be repeated for every field value $B$ one is interested in.

3 Numerical examples

The following examples are evaluated for the respective Heisenberg Hamiltonian

[TABLE]

where $\!\hbox{}\mathop{\vphantom{\vec{s}}\smash{\vec{s}}}\limits_{\sim}\hbox{}\!$ denotes the quantum spin operators. In the following examples the single-spin quantum number is always $s=1/2$. For the linear chain with antiferromagnetic nearest-neigbor exchange interaction $J_{ij}=J$ for adjacent spins on the chain; for the sawtooth chain the spins are connected by $J_{1}$ and $J_{2}$ as depicted in Fig. 3.

3.1 Spin chain

As our first example we would like to discuss the evaluation of magnetic observables for a linear spin chain of $N=20$ spins $s=1/2$ with antiferromagnetic nearest-neighbor interaction. This example, which was already investigated by Feiguin and White [47], serves to show that the method is indeed capable of reproducing exact results with rather high accuracy. Figure 2 shows the specific heat capacity (top & center) and the specific zero-field susceptibility (bottom) plotted vs temperature. The values $m=10$ (green), $m=20$ (blue), $m=40$ (red), and $m=60$ (orange) have been used in the approximation; black solid curves represent the respective observables obtained from exact diagonalization [12, 60]. The zero-field specific heat is evaluated using the variance of the Hamiltonian, whereas the susceptibility is evaluated as difference quotient, i.e.,

[TABLE]

Interestingly, this yields better results than calculating the susceptibility as a variance of the magnetization. As can be seen in the panels of Fig. 2 the results improve as $m$ increases. For the susceptibility the approximation lies almost on top of the ED curve already for $m=20$ and $m=40$, whereas the low-temperature shoulder of the specific heat needs at least $m=60$ to be reasonably well approximated, compare Fig. 2 (center).

3.2 Sawtooth chain

The situation changes drastically when moving to a sawtooth (delta) chain system with ferromagnetic interaction $J_{1}$ between adjacent apical and basal spins and antiferromagnetic interaction $J_{2}$ between adjacent basal spins, compare Fig. 3. We selected this system on purpose since it exhibits a quantum phase transition for $\alpha=|J_{2}/J_{1}|=1/2$ [61, 62, 33]. At the transition point there is a flat one-magnon excitation band that leads to extraordinary features in the low-energy spectrum [61]. In particular, 50 % of the total available entropy is contained in the massively degenerate ground state. Moreover, even for finite sawtooth chains the exication gap is extremely small and the excitations above the massively degenerate ground-state manifold are also highly degenerate. As a result, low-temperature physics is unconventional [61]. In the vicinity of the critical point the degeneracy splits up and leads to an additional well-pronounced low-temperature maximum of the heat capacity. The position of the maximum corresponds to an additional low-temperature scale $T_{\text{low}}$. One would expect that the imaginary time evolution of the DMRG matrix product state runs into problems when this temperature scale is reached from above and thus the state $|\Psi_{\beta}\rangle$ enters a part of Hilbert space with an enormous spectral density and entropy contribution.

Figure 4 displays the results for a sawtooth chain with $N=32$ and open boundary conditions (OBC) as black solid curves for $\alpha=0.45$ (top), $\alpha=0.50$ (center) and $\alpha=0.55$ (bottom); $m=40$ was used for the calculations. This result is compared to FTLM approximations which are virtually exact as discussed elsewhere [42]. In the case of periodic boundary conditions (PBC) translational symmetry was employed by means of the freely available program spinpack [63] (Therefore, much smaller $R$ suffice.). The blue curves display FTLM results that correspond to the situations investigated by ThDMRG, the red and green curves show FTLM results of closely related cases of $N=32$ (PBC) as well as $N=31$ (OBC, no dangling bond, see Fig. 3), respectively. As one can visually verify, the approximation by means of thermal DMRG becomes untrustworthy when the temperature approaches the low-energy scale $T_{\text{low}}$, i.e., the realm of the low-temperature maximum. The reason is that the low-energy part of the spectrum and thus of Hilbert space can (effectively as a whole) be easily represented by a matrix product state (of smaller $m$) as long as the state $|\Psi_{\beta}\rangle$ resides at higher temperatures. But at temperatures of the order of $T_{\text{low}}$ and below, the matrix product state would need to resolve the spectral decomposition at this scale, which becomes virtually impossible in view of the vast dimension that is exponential in system size. This phenomenon worsens the problem of time-dependent DMRG to become unreliable after a characteristic runaway time [64, 65]. Independent of this time scale, $T_{\text{low}}$ introduces a system specific scale where ThDMRG breaks down. An increase of the block dimension from $m=40$ to $m=60$ improves the situation somewhat, compare center of Fig. 4, but only for a limited temperature range. In order to maintain a good accuracy an exponential increase of block dimension would be necessary with time [64, 65].

Despite the failure of thermal DMRG at very low temperatures one must acknowledge that it successfully approximates the thermal behavior above $T_{\text{low}}$. In particular, the Schottky peak around $T\sim|J_{1}|\sim|J_{2}|$ is very well reproduced.

In order to go a bit into details, we would like to discuss the behavior of the entanglement spectrum with imaginary time. This serves as a measure of how good the matrix product state can approximate the target state, which could for instance be the ground state or $|\Psi_{\beta}\rangle$. A spectrum that does not fall off rapidly, signals that a subsystem, i.e. block $A$, is strongly entangled with the environment, i.e. block $B$, and thus – for the chosen value of $m$ – the matrix product state is not a good approximation of the target state.

Figure 5 shows on the l.h.s. the entanglement spectra of the reduced density matrices of block $A$ (half the system size) for a chain and on the r.h.s. for a sawtooth (delta) chain with $\alpha=0.5$, both of length $L=32$. One notices that for both systems the spectra fall off slower and slower with increasing imaginary time, i.e. inverse temperature. This means that the entanglement with the environment grows and the matrix product state of fixed block dimension $m=40$ is less and less capable to represent the thermal state at the respective $\beta$. This effect is somewhat more pronounced for the sawtooth chain, but not as dramatically as the poor approximation of the low-temperature observables would suggest. This means that the entanglement spectrum alone is not a sufficient quantifier of the achieved quality of calculation.

4 Summary

In this article the thermal DMRG was applied to two archetypical spin systems – the linear antiferromagnetic chain with nearest-neighbor exchange interaction as well as the frustrated sawtooth (delta) chain close to the quantum critical point. Our intention was an investigation of the achievable accuracy in view of a later use of the method in molecular magnetism. The general conclusion seems to be that non-trivial low-temperature features are hard to resolve and approximate properly, since such features stem from intricate low-energy spectral properties, which would require a good representation by the matrix product state at this temperature. The latter becomes increasingly unlikely since both, the intricate low-energy spectrum as well as the unavoidable runaway time of DMRG impede a proper representation of the thermal state of the system in terms of a product state of low-dimensional matrices.

Nevertheless, the accuracy for thermal magnetic observables is very good above those scales. And, most importantly, the major advantage to be able to treat very large and possibly frustrated spin systems cannot be accomplished by any other approximate method.

Acknowledgment

This work was supported by the Deutsche Forschungsgemeinschaft DFG (314331397 (SCHN 615/23-1); 355031190 (FOR 2692); 397300368 (SCHN 615/25-1)). The authors thank Andreas Klümper and Salvatore Manmana for helpful discussions.