Thermodynamic Behavior of Spin-1 Heisenberg Chain: a Comparative Study

In this work, we present a comparative study of thermodynamic quantum equilibrium observables of spin-1 Heisenberg chain. As a theoretical approach, the modified spin-wave theory is chosen. From the numerical side, we consider finite-temperature Lanczos, kernel polynomial, and density matrix renormalization techniques. The results show general consistency of thermodynamic quantities, such as heat capacity and magnetic susceptibility, calculated within the approaches. For the modified spin-wave, the results show some inaccuracy especially at low temperature, failure in capturing the gapped features of spin-1 Heisenberg chain, while the numerical approaches illustrate a good achievement at T→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T\rightarrow 0$$\end{document} within the proper parameters chosen.


Introduction
While the exact diagonalization technique has been used extensively to explore the ground state study of quantum many-body [1][2][3] and also dynamic in spin chains [4], it is of less use at finite temperature. This can be attributed to the exponential growth of Hilbert dimension of systems, which for the spin models with spin-S and lattice size L, it grows as dim = (2S + 1) L . For spin-1/2 the ground state normally can be found for up to L ≲ 30 lattice sites, while the full spectrum accessible for systems only up to L ≲ 16 . However for spin-1 models, accessible systems size scales about to half. Despite exploring the ground states, it is also of interest to study thermodynamics properties as can be accessible experimentally.
The first goal of this paper is to address a comparative study of some approaches to explore accuracy in capturing the thermodynamics properties of spin-1 Heisenberg chain. As a second purpose, we present technical details of the implementation of the methods. Although the particular emphasis of this paper is on numerical techniques, to have a self-contained context, we also revisit the modified spin-wave theory (MSW) [5] where has been invented to overcome the shortcoming of the conventional spin-wave theory (SW) [6][7][8] at finite temperature. For numerical approaches, we consider the finite-temperature Lanczos method (FTLM) [9,10], the kernel polynomial method (KPM) [11], and finite-temperature density renormalization group (DMRG) [12,13] technique.
The FTLM allows us to evaluate thermodynamic observables approximately, without full numerical diagonalization of the Hamiltonian. The key approximations made in the FTLM are (i) the trace estimator of operator Ô and (ii) the Krylov space expansion of the exponential Boltzmann factor exp (− H) by the Lanczos method, where is the inverse temperature. We consider this approach as a benchmark against other approaches in this work. The second approximation in FTLM can be replaced by other approaches which solve the imaginary-time Schrödinger equation iteratively such as Chebyshev polynomials [11]. The main privilege of KPM approach is that the involved Chebyshev polynomials expansion does not suffer from the loss of orthogonality during the Lanczos procedures used in Krylov space methods.
On the other hand, Hilbert space dimension growth with system size is a fundamental problem of Lanczos-based approaches. The density matrix renormalization group (DMRG) approach circumvents this problem by discarding the less important states in a systematic way [14,15]. Likewise to the Lanczos-based methods, finite-temperature DMRG methods also invented [12,13,[16][17][18][19].
As said, we aim to explore the mentioned techniques accuracy to capture correct thermodynamics observable and address proper parameters or possible tricks to remove some unphysical behaviors. The paper is organized as follows. In Sect. 2, we introduce the model and methods. In Sect. 3, we present results and discussions.

Models
We consider archetypical quantum spin systems described by the Heisenberg spin-1 model is the spin vector at site i, J < 0 (> 0) exchange coupling with ferromagnetic (antiferromagnetic) model, and sum is limited to nearest neighbor. The quasione-dimensional quantum magnets like CsNiCl 3 [20][21][22], Ni(C 2 H 8 N 2 ) 2 NO 2 ClO 4 (NENP) [23,24], or LiVGe 2 O 6 [25,26] are well described within it. The second term is singleion anisotropy. In the absence of anisotropy term K = 0 and at low temperature, the classical ground of Eq. (1) depends on exchange coupling is ferromagnetic (colinear Néel) order with the classical energy E cla ∕LS = zJS(−zJS) , where z is the coordinate number. It is worth mentioning that continuous symmetries cannot be spontaneously broken at finite temperature in systems with sufficiently short-range interactions in dimensions d ≤ 2 [27].
For the S = 1 case which is the main focus of this work, the first term with antiferromagnetic exchange interactions stabilizes the Haldane phase [28,29] with a unique disordered ground state and an excitation gap. When the single-ion anisotropy K > 0 is tuned, the system undergoes a phase transition between two gapped, symmetric phases at K ≈ 1 [30 -32]. It has been shown numerically that the Haldane phase is adiabatically connected to the AKLT state [33]. The so-called large-K phase is adiabatically connected to a simple product state.

Modified Spin Wave Method
To have access to the thermodynamic properties, one theoretical approach would be the conventional spin-wave (SW) theory. However, the Dyson-Maleev [34,35] or the Holstein-Primakoff (HP) [36] SW theory works well for the ground state [6][7][8]. To develop a finite-temperature SW theory, Takahashi [5] introduced a constraint in respect to the Mermin-Wagner theorem [27] with the vanishing local magnetization, ⟨S z i = 0⟩ . The modified SW by Takahashi's has shown drawbacks to capture thermodynamics properties, e.g., failing to see a Schottky-like peak of the specific heat and also signaling the first order artificial phase transition to the trivial paramagnetic solution [5], which pertains to meanfield divergence issue [37]. In this work, we specifically follow a slightly simple approach introduced in Refs. [38,39], as it is shown good applicability for a wide range of temperature for S = 1∕2.
For bipartite lattices, one can separate the lattice (with 2N sites) into two sublattices, where sublattice A (B) is defined for the up (down) spins in the Néel order state. Following Takahashi [5] and Refs. [38,39], we introduce a Lagrange constraint into the Hamiltonian, H = ∑ i∈A,B l S z i . The Lagrange multipliers are assumed as i = if i ∈ A and − if i ∈ B , and is redefined as ≡ S(Jz + K)( − 1) for convenience. This Lagrange Hamiltonian fulfills the Mermin-Wagner theorem with a constraint of zero staggered magnetization. This leads to a finite gap in the spin-wave energy spectrum and ends up the number of bosons not diverging with increasing temperature.
Starting from Hamiltonian in Eq. (1) one needs to rotate spins on the B sublattice to ensure an opposite magnetization direction for two sublattices, which can be done by S x → −S x , S y → S y , and S z → −S z . This leads new sets of HP operators, Application of the HP to the H msw = H + H leads to the following effective magnon model where = {a, b} refers the sublattices A and B, respectively. The first term is magnon on-site energies = JzS − , and the second term is magnon hopping between different sublattices. Now, by defining ∕z is structure factor. We see that the states with k and −k are coupled and the Hamiltonina is not number conserving. Therefore, directly diagonalizing Eq. (3) will not give us the eigenvalues. The quadratic part of the above Hamiltonian can be diagonalized via a standard Bogoliubov transformation The quasiparticle operators k and k satisfy bosonic commutation relations provided u 2 k − 2 k = 1 . Then, diagonalized Hamiltonian can be read as follows

The free energy per site is given by
To fulfill the vanishing staggered magnetization constraint we minimize the free energy with respect to , which leads to a self-consistent equation Having energy at finite temperature, thermodynamics observable including the free energy, the internal energy, the entropy, and the specific heat are directly accessible as follows: where n k = 1 e k ∕T −1 is Bose-Einstein function. Then the magnetic susceptibility is given through the z-component correlation func- Using the Holstein-Primakoff transformation, ⟨S z i S z j ⟩ one arrives the following formula for linear spin-wave approximation [40]

Finite-temperature Lanczos Method
From a numerical point, measuring the thermodynamic average of an operator A at finite temperature needs the access of all eigenstates �Ψ i ⟩ and corresponding energies E i , where = 1∕T is the inverse temperature (with k B = 1 ) and N st is the dimension of Hamiltonian H . One possible approach would be the full exact diagonalization(ED) which is restricted to small systems size as dimension of Hamiltonian growth exponentially, e.g., for S = 1 it grows as 3 L , and normally it is possible to have access to all of the eigenstates for L < 12.
To overcome such obstacle, the FTLM [9, 10] method evaluates of the expectation value in Eq. (8) using a general orthonormal basis �n⟩ ( e.g., the Lanczos basis), However, the above equation still involves the summation over the complete set of N st states �n⟩ , which is not doable in practice. The FTLM introduces two approximations. The first one is to replace the summation of n with random sampling r with R times. This approximation is also called typicality or (microcanonical) thermal pure quantum states [9,[41][42][43], the second one is the number of Krylov subspace dimensions M. Then we get where �r⟩ is normalized random initial state and � r j ⟩ are an eigenvector in the M-the Krylov subspace of Hamiltonian H , with corresponding energy r j . The thermodynamic of interest can be read as follows: At high temperatures, a few samplings are enough for high accuracy as the error of physical quantities is proportional [10]. However, at T → 0 to have the correct result, the ground state must be well converged, i.e., � r 0 ⟩ ∼ �Ψ r 0 ⟩ which needs more sampling R. To avoid this, the low temperature Lanczos method (LTLM) has been proposed [44]. The idea of LTLM is very simple, rewrite Eq. (9) in symmetric form now with inserting two set of Lanczos basis, we end up with a double sum

Kernel Polynomial Method
The FTLM relies on a Krylov space expansion for Boltzmann weights exp (− H) and shown to produce accurate approximations with averaged over random sampling r [45][46][47][48]. However, one has to worry about errors in the orthogonality during recursive state generation used in the Lanczos method.
To avoid this, an alternative approximation introduced by authors [11], which is based on the expansion of the density of states; (E) = ∑ i (E − E i ) in terms of Chebyshev polynomials. This property has already been proved by achieving high accuracy results in studying numerical unitary time evolution, see, e.g., [11,[49][50][51]. In the following, we give a concise description of the KPM, while for a more detailed see Ref. [11].
Quantity of interest in a variety of systems can be written in the general form where is introduced to avoid truncation of the approximated delta peaks [11]. After this transformation, the Chebyshev moments n are defined as (16) These traces are estimated using the typicality approach with randomly chosen states �r⟩ [11]. For the numerical method, the number of moments in Eq. (18) is finite and it causes so-called Gibbs oscillations, which leads to an approximated density of the states to have a negative value. To minimize these a convolution of f(x) with a kernel is applied. In this paper, we are using the Jackson Kernel [11] whose coefficients read Then rescaled density of states is expanded in terms of Chebyshev polynomials T n (x): where set {x k } choose as Having a density of state, thermodynamic quantity like canonical partition function and heat capacity can be evaluated as where and can be computed with discrete cosine transform (type III) with standard scientific library.

Finite-temperature Density Matrix Renormalization Group
The original DMRG method [14,15] was conceived for T = 0 . Soon after, finite-temperature variants have been developed, such as a quantum transfer-matrix formulation [16,17], exact representations through purification [12,13], minimally entangled typical thermal states [18], and multiplication of matrix-product operators [19]. In this paper we take the purification route, which enlarges the Hilbert space, introducing auxiliary sites (called "ancillas"). This doubles the size of the lattice and allows one to work directly with the grand canonical density matrix which is expressed as a partial trace over a pure state � ⟩ by tracing out auxiliary degrees of freedom Q that encode the thermal bath.
The state at infinite temperature � = 0⟩ can be easily prepared as the ground state of the Hamiltonian where a(i) indicates the ancilla site attached to the physical site i. This initial state can be evolved as a pure state within the larger space as � ⟩ = e − H∕2 � = 0⟩ , whereby the Hamiltonian H only acts on the physical sites. A major advantage of this approach is that all the standard methods for time-evolving quantum states within DMRG are directly applicable. The internal energy is obtained by taking E = ⟨H⟩ , while the specific heat per site can be directly computed as the energy variance in the physical system c = 2 ∕L The magnetic susceptibility follows from the application of a small magnetic field in z-direction

Results and Discussion
We begin with MSW theory, so let us cite numerical results of the self-consistent equation Eq. (6). The temperature dependence of the Lagrange multiplier is shown in the main panel of Fig. 1.
For K = 0 and at zero temperature the excitation spectrum has a zero gap at k = 0 which results in a long-range antiferromagnetic order [5,52]. While as argued by Haldane [28, (26)

29]
, there is an excitation gap for integer spin and long-range order could not occur. While as can be seen from inset of Fig. 1, at low but finite temperatures, the long-range magnetic order is destroyed by quantum fluctuations. The spin-wave thus behaves as bosons with a finite gap Δ = zJ √ (1 + K∕zJ) 2 2 − 1 , which results in the destruction of the long-range antiferromagnetic order at any finite temperature. The inset of Fig. 1 shows the excitation spectrum with increasing temperature. It can be seen for finite temperature the gapless Goldstone mode vanishes. The effect of anisotropy K on the Lagrange multiplier is shown on the main panel, which increases with decreasing K.
Having , we calculate the thermodynamic quantities and compare them with numerical approaches. So, here we address the parameters used in each numerical technique separately.
For FTLM, the Krylov subspace dimension and the number of sampling are fixed through the paper as M = 50 and R = 100 . Indeed, for all possible accessible system sizes, these two parameters show persistent results, and no numerical instabilities are seen. We consider a chain with length L = 18 on a ring to reduce the edge effect. The results of FTLM are considered as a reference to compare with other approaches accuracy.
For KPM, there are four parameters configuration, as listed in Table 1. Throughout the paper we fix kernel parameter to g n = 1 and truncation parameter = 0.01 , as the influences of these two parameters are already explored for spin-1/2 [53]. Chain size is L = 18 , if it is not stated, and placed on a ring. As discussed in Ref. [53], for the number Fig. 1 The temperature dependence of the Lagrange multiplier in MSW theory of spin-1 antiferromagnetic chain for two single-ion anisotropy parameters. Inset shows, the excitation spectrum k vs k∕ for single-ion anisotropy K = 0 with changing color from cold to hot temperature of points to integrate discrete cosine-transform Ñ , a possible optimized choice is to set it equal to the number of moments N . We only change the number of moment and sampling parameters, and their influences on the accuracy are discussed.
The finite-temperature DMRG, implemented using the C++ iTensor library [54]. System size L = 41 and we approximate the evolution operator by a matrixproduct operator (MPO) [55] with a fourth-order [56] Trotter scheme. For initial "wave-function," we choose maximally entangled state between neighboring sites We adapt the open boundary condition DMRG and its effects are discussed.
As said before, we are interested in the thermodynamic observables. So we consider the heat capacity and the magnetic susceptibility, as are plotted in Figs. 2 and 4, respectively. Different methods output are plotted together, and two single-ion anisotropies results are depicted in separate (a) and (b) panels.
One can see that results of the heat capacity (Fig. 2) match perfectly at high temperatures, namely at T∕J ⪆ 2 , although the result of MSW shows offset less than 5% . All method gives the main (Schottky anomaly) peak about T∕J ≈ 1 with some deviation in the maximum.
Insets zoom on low temperature with the logarithmic scale. One can see that for K = 0.0 case, FTLM and DMRG results are very close, while a deviation at T → 0 is seen for MSW and KPM cases. In the MSW, this shortcoming is caused by the failure of the method to correctly capture the Haldane gapped ground state. This is also apparent in Fig. 2, where single-ion anisotropy already develops a gap in the model.
For the KPM technique, the problem of inaccuracy is explored in Fig. 3. As stated, the parameters of in charge are the number of moment N , the number of sampling R , Table 1 The parameters configuration used in the KPM. The kernel g n and truncation parameters are fixed and chain length L. We noticed that the number of sampling R has no major effect on capturing gap feature at T → 0 , although increasing R smooths the curve. Then, we consider two chain lengths L = 14 , and L = 18 , and let the number of moment changing between N = 100 − 800. In Fig. 3, influence of the number of moment is apparent. For temperature range T∕J ⪅ 0.1 , we have found that the heat capacity gets spurious increase with decreasing temperature. This issue gets remedy by increasing the number of moment, as it can be seen that for a fixed temperature, the heat capacity decreases with increasing the number of moment. The optimized number of moments N is also dependent on the chain length, compare panels (a) and (b) of Fig. 3. Figure 4 shows the magnetic susceptibility versus temperature for two single-ion anisotropies. The results of all methods are plotted together for comparison. The numerical approaches outcomes show a good match in all temperature range, and the gap evolution of the model is well captured by numerical techniques, namely FTLM, KPM, and DMRG. On the other hand, the result of MSW shows an offset in the temperature window presented. Although we could not see any gap evolution evidence within this approach, the sharp increases of the magnetic susceptibility at low temperature missed the gap feature in the model (see Fig. 4(b)). Notwithstanding straight access of correct results at low temperature within the FTLM technique, we noticed some subtlety to access the low temperature results in KPM and DMRG techniques, which explained in the following.
For the KPM, we illustrate effect of the number of moments in Fig. 5. As can be seen for T∕J → 0 , the magnetic susceptibility diverges. This issue gets worse for bigger chain length, as for instance from chain L = 14 to L = 18 , we found an order of increases. As the heat capacity case, we fixed a temperature and monitories the number of moments N effects. What we found, can be clearly seen from the inset panels of Fig. 5. The incriminating number of moments removes the spurious contribution of low temperature energies on the magnetic susceptibility.
Regarding the finite-temperature DMRG technique, to converge a correct result specifically at low temperatures, one needs to approximate infinite systems. We consider two possible tricks, namely periodic and soft (setting the edge spins to 1/2) boundary conditions. We found that periodic boundary condition perform badly, especially for T → 0 , so here we just address the soft boundary condition. Figure 6 shows the magnetic susceptibility computed with DMRG for the open boundary condition (OBC) and soft boundary condition (SBC). One can see for finite chain length in OBC, the magnetic susceptibility diverges at T → 0 and the situation getting worse when the single-ion anisotropy increases and gap becomes profound in the model. Considering the SBC trick, which can be easily implemented in the DMRG method with replacing edges spin-1 with spin-1/2 (see the cartoon scheme depicted in Fig. 6), we converge to correct behavior at low temperature.
To summarize, we conduct a comparative study on the accuracy of theoretical and numerical methods of capturing the thermodynamic observables. As a case model, gapped spin-1 Heisenberg chain is considered. We have seen, although MSW method could capture the general thermodynamic observables trend over the wide temperature range, but it could not detect gapped features of spin-1 Heisenberg model in the thermodynamic limit. Comparing the numerical methods, namely FTLM, KMP, and DMRG approaches, we noticed that by choosing good parameters for KPM or applying a proper trick for DMRG one achieves accurate results. Specifically, at low temperature KPM and DMRG needs subtlety treatment. Fig. 6 The magnetic susceptibility plot computed with DMRG method with open boundary condition (OBC) and soft boundary condition (SBC). Insets show the zoom of low temperature behavior and also a cartoon scheme illustrates spin chain in two boundary conditions