Magnetocaloric effect in quantum spin-s chains

We compute the entropy of antiferromagnetic quantum spin-s chains in an external magnetic field using exact diagonalization and Quantum Monte Carlo simulations. The magnetocaloric effect, i.e., temperature variations during adiabatic field changes, can be derived from the isentropes. First, we focus on the example of the spin-s=1 chain and show that one can cool by closing the Haldane gap with a magnetic field. We then move to quantum spin-s chains and demonstrate linear scaling with $s$ close to the saturation field. In passing, we propose a new method to compute many low-lying excited states using the Lanczos recursion.


Introduction
The magnetocaloric effect, i.e., a temperature change induced by an adiabatic change of an external magnetic field was discovered in iron by Warburg in 1881 [1]. This effect also has a long history in cooling applications. For example, adiabatic demagnetization of paramagnetic salts was the first method to reach temperatures below 1K [2]. Among the many experimental but also theoretical investigations to follow, we would like to mention in particular the computation of the field variation of the entropy in a quantum spin-1/2 XXZ chain [3]. This investigation is remarkable in so far as it is one of the very first numerical computations for an interacting quantum many-body system. The relation of the entropy of quantum spin chains in a magnetic field to magnetic cooling was also noted [4] and even one exact computation was performed, albeit for a chain with a strong easy-axis anisotropy [5].
In this paper we will investigate magnetocaloric properties of quantum spin-s Heisenberg chains given by the Hamiltonians J > 0 is the antiferromagnetic exchange constant, h an external magnetic field, L the length of the chain, and S i are quantum spin-s operators at site i. We will use periodic boundary conditions, i.e., S L+1 = S 1 . In section 2 we will first summarize the methods which we are going to use, namely exact diagonalization (ED) and Quantum Monte Carlo (QMC) simulations using the Stochastic Series Expansion (SSE) framework. Then we will apply these methods to the s = 1 chain in section 3. The spin s = 1 Heisenberg chain is famous for the presence of the so-called 'Haldane' gap at h = 0 [27][28][29][30]. In particular, we will illustrate with the s = 1 chain that one can exploit the presence of a spin gap for cooling by adiabatic magnetization at low magnetic fields. Next, we discuss the scaling with spin quantum number s in section 4. In particular, we will show that there is a quantum scaling regime at low temperatures close to the saturation field where one has linear scaling with s. Finally, in section 5 we will conclude with a summary.

Exact diagonalization
Thermodynamic quantities can be computed using spectral representations. For example, the entropy can be written as [10] where Z = n e −En/T is the partition function and E n are the eigenvalues of the Hamiltonian. In order to evaluate spectral representations like (2), one needs to diagonalize the Hamiltonian. If this is done numerically exactly, this is called Exact Diagonalization (ED). First, one should perform a symmetry reduction. We have used translations and S z -conservation as well as reflections and spin inversion where appropriate. One can then use a library routine to perform a full diagonalization. Such an approach is very much in the spirit of classic work [3,31,32]. The main technical differences of our computations and [3,31,32] are: (i) we have exploited SU (2)-symmetry to reconstruct the S z = 1 sector from the S z ≥ 2 sectors and the spin-inversion resolved S z = 0 sector, (ii) we have used improved implementations of the diagonalization routines [33,34], (iii) and we have substantially more powerful computers at our disposal. The combination of these factors enables us to compute full spectra for bigger systems than in the 1960s.
If we are interested only in low-energy properties, we can use iterative diagonalization algorithms like the Lanczos method [35][36][37].
The basic Lanczos algorithm [35][36][37] for a Hermitian matrix H proceeds as follows. For a given normalized start vector v 1 , one defines a sequence of normalized Lanczos vectors v j via the recurrence relations with the initial conditions β 1 = 0 and v 0 = 0. The real coefficients α i and β i define the so-called Lanczos matrices The crucial point of the Lanczos algorithm is that the eigenvalues of T n yield a good approximation of the extremal eigenvalues of H already for n much smaller than the total Hilbert space dimension. The Lanczos algorithm suffers from one practical problem: the recursion relations (3) are supposed to guarantee mutual orthogonality of the vectors v i , i.e., v i · v j = 0 for i = j. However, this fails to be correct when the computations are carried out numerically, leading to undesired spurious states, called ghost states [36,37]. Several strategies have been proposed [36,37] in order to deal with the ghost problem during the computation of excited eigenvalues.
Here we suggest another method to compute excited states with controlled accuracy and correct multiplicity. We start by performing a fixed number n of Lanczos iterations with a given start vector v 1 . Next we compute the orthogonal transformation U i,j , 1 ≤ i, j ≤ n diagonalizing the Lanczos matrix (4). Then we fix a number m ≪ n and repeat the Lanczos procedure with the same start vector v 1 . During this second Lanczos pass we construct m vectors The second Lanczos pass is needed in order to avoid storing the n ≫ m Lanczos vectors v j . The vectors u i given by (5) yield approximations to the eigenvectors of H. However, due to the ghost problem, there are spurious vectors which need to be eliminated. We perform this in two steps. First, using the observation that the vectors u i are supposed to be mutually orthogonal, we can eliminate those u i which have big projections on the u j with j < i. Second, we reorthogonalize the remaining vectors yieldingm ≤ m orthonormal vectors˜ u i . Finally, we project the matrix H onto this subspace via A full diagonalization of them×m matrix H i,j and application of the resulting basis transformation to the vectors˜ u i yieldsm orthogonal vectors w i which approximate the eigenvectors of H. The eigenvalues E i associated to the vectors w i and their accuracy can be estimated according to According to our experience, there are neither any missing nor any spurious eigenvalues among the converged ones. It should be mentioned that there is no guarantee that the Lanczos procedure yields the complete spectrum for a given start vector v 1 , in particular if one has not performed a complete symmetry decomposition. Nevertheless, numerical noise seems to prevent this from happening, at least for sufficiently generic start vectors v 1 .
For our procedure one needs to choose the number of Lanczos iterations n and the dimension of the subspace m ≤ n in advance and the whole procedure needs to be repeated if it turns out that less than the desired number of converged eigenvalues have been obtained. In practice, however, one obtains results of a comparable quality for a class of problems where the choice of n and m has been adjusted for one representative case. In the examples to be reported below, we have been able to obtain about 100 eigenvalues in a given symmetry sector with a relative accuracy of 10 −8 or better using n = 2000 Lanczos iterations and m = 500 vectors u i .

Quantum Monte Carlo
Since the system sizes accessible by either full diagonalization or Lanczos diagonalization are limited, it is desirable to have other methods at our disposal. For non-frustrated spin models like the spin-s Heisenberg chain (1), one can in principle use Quantum Monte Carlo (QMC) simulations. There is, however, one problem: we are particularly interested in the entropy which is usually obtained with large statistical errors from integrating a Monte-Carlo result for the specific heat. In this subsection we will summarize how one can circumvent this problem.
In order to obtain a QMC estimate of the temperature dependence of the entropy, we employed an extended ensemble, broad-histogram method [38][39][40] within the Stochastic Series Expansion (SSE) framework [41][42][43][44]. Based on this approach, thermodynamic quantities can be obtained over a broad range of temperatures from a single QMC simulation that provides estimates of the expansion coefficients g(n) of the system's partition function Z in a high-temperature series expansion: Here, β = 1/T is the inverse temperature. In terms of the expansion coefficients g(n), the internal energy is obtained using and the free energy from Finally, the entropy can be calculated using S = (E − F )/T . In the QMC simulation, estimates of the first Λ coefficients g(n), n = 0, ..., Λ are obtained by performing a random walk in the expansion order n, such as to sample efficiently all expansion orders from n = 0 (where g(0) = (2s + 1) L for a spin-s chain of L sites is known exactly) to n = Λ. This is accomplished using an extension of the Wang-Landau flat-histogram sampling algorithm [45,46] to the quantum case, as detailed previously in Refs. [38][39][40]. Knowledge of the first Λ coefficients g(n) allows for the calculation of thermodynamic quantities from β = 0 down to a temperature T m = 1/β m for a system of L sites, where Λ scales proportional to β m × L [38]. In particular, for the spin s = 1 Heisenberg chain, we were able to treat systems with up to L = 40 sites, accessing temperatures down to T m ≈ 0.03, which required an already large value of Λ = 5000. The data shown in this contribution (and also those in [14]) were obtained based on the original version of the algorithm [38], while recently an improved sampling strategy was proposed, based on optimizing the broad histogram Quantum Monte Carlo ensemble [40].

A test case
We test and compare the aforementioned numerical methods using the example of the s = 1 Heisenberg chain (1) with a finite magnetization. Fig. 1 shows the entropy per site S/L as a function of temperature T for a fixed magnetic field h = 3 J. For L = 14, we have performed a full diagonalization. The corresponding ED curve can therefore be considered as the exact result for L = 14.
A full determination of the spectrum is clearly out of reach already for L = 20. In this case, we therefore had to perform a severe truncation to low energies. Thus, the L = 20 ED curve is only a low-temperature approximation. In this paper, the precise temperature range is fixed as follows: Let E be the highest energy until which the spectrum is definitely complete. Then we restrict to temperatures T ≤ E/5. This choice ensures that missing states are suppressed by a Boltzmann factor ≤ exp(−5) ≈ 0.006738. This may seem a small number, but since we are discarding many states, extending truncated data to higher temperatures would yield artifacts which are clearly visible on the figure(s).
QMC (denoted by SSE in Fig. 1) is able to treat the larger system size L = 40. This method, however, is best suited for high temperatures. Indeed, at T ≈ 0.05 J one can see deviations caused by statistical errors in the SSE result of Fig. 1 which prevented us from going to lower temperatures. Finite-size effects can be observed in Fig. 1 in the L = 14 ED curve for T 0.4 J and in the L = 20 ED result. Otherwise it is reassuring that we observe good overall agreement between all three methods.
In passing we note that, for the value of the magnetic field used in Fig. 1, the low-energy physics of the s = 1 chain is described by a Luttinger liquid (see Ref. [47] and references therein). It is well known that the specific heat C of a Luttinger liquid is linear in T [48]. Due to the relation C = T (∂ S/∂T ) and because of S(T = 0) = 0, the entropy of a Luttinger liquid is identical to its specific heat and in particular also linear in T . Indeed, Fig. 1 is consistent with a linear behavior S ∝ T at sufficiently large L and low T .

Spin s = 1 Heisenberg chain
In this section we will discuss the entropy of the s = 1 Heisenberg chain in more detail. First, we show results for the value of the entropy per site S/L as a function of h and T in Fig. 2. The ED curves have been obtained by computing the entropy S on a mesh in the h-T -plane and determining the constant entropy curves, i.e., the isentropes from this data. This leads to some discretization artifacts at the cusps in the the S/L = 0.05 curves (the ED curves at the lowest temperature) in Fig. 2. Conversely, the QMC had problems to resolve the large low-temperature entropy around the saturation field h sat. = 4 J. Accordingly, the L = 40 SSE data points with S/L = 0.05 are missing at h/J = 3.7 and 4. From our QMC simulations we can just conclude that for an L = 40 chain these two points are in the region T < 0.067 J and T < 0.085 J for h = 3.7 J and 4 J, respectively. The ED curves with S/L = 0.05 in Fig. 2 have marked finite-size wiggles. Also at higher temperatures (entropies), finite-size effects can be observed. There are, however, two regions where finite-size effects are evidently small. One is the low-temperature regime for fields h smaller than the Haldane gap [29,30] ∆ ≈ 0.4105 J . (11) Also in the gapped high-field regime above the saturation field h 4 J one observes essentially no finite-size effects at any temperature. This renders ED the method of choice at high fields h h sat. , in particular at low temperatures, while QMC is preferable otherwise because bigger system sizes L are accessible. One can read off from Fig. 2 pronounced temperature changes in two regimes. Firstly, one can obtain cooling by adiabatic demagnetization from a high magnetic field h > h sat. as one lowers the field to the saturation field. Secondly, at low temperatures, one can also cool by adiabatic magnetization from h = 0 to h = ∆. This demonstrates that one can use the gap-closing at a generic field-induced quantum phase transition for cooling purposes. In the regime 0.4105 J ≈ ∆ < h < h sat. = 4 J, the spectrum is gapless (compare Fig. 1 and the related discussion in section 2.3). Accordingly, one observes only small temperature changes induced by adiabatic (de)magnetization in this field range.
We now focus on the low-field regime. Firstly, we consider the behavior of the entropy at h = 0 which is shown in Fig. 3. One observes that finite-size effects are small for L ≥ 14: the L = 14 curve is almost within the error margins of the L = 40 QMC curve even at low temperatures. This can be attributed to a finite correlation length ξ which is evidently sufficiently smaller than 14 for all temperatures at h = 0. Indeed, the correlation length of the s = 1 Heisenberg chain is known to be ξ ≈ 6 at h = 0 and T = 0 [29,30].
Because of the presence of a gap ∆ at h = 0, we expect the entropy to be exponentially activated as a function of temperature. Accordingly, the low-temperature asymptotic behavior at If we fix the gap to the known value ∆ = 0.4105 J [29,30], the only free parameter in the lowtemperature asymptotic behavior (12) is the prefactor A. A fit of the L = 40 QMC results yields A = 0.577(1). The dashed line in Fig. 3 shows that the formula (12) describes the low-temperature regime quite well with the given parameters. Secondly, we consider a magnetic field exactly equal to the Haldane gap h = 0.4105 J. Fig. 4 shows a doubly-logarithmic plot of the entropy S as a function of temperature T for this value of the magnetic field. Since we are now sitting at a field-induced quantum phase transition, we expect a large number of low-lying states and relatedly an infinite correlation length ξ = ∞. On a finite system, these low-lying states will be pushed to higher energies. Indeed, we observe pronounced finite-size effects at low temperatures in Fig. 4.
At a quantum phase transition in one dimension which preserves a U (1)-symmetry, the entropy S is expected [4,6,7,10] to vary asymptotically as a square root of temperature: As before, the amplitude B is the only free parameter. A fit to the L = 40 QMC data with T ≤ 0.07 J yields B = 0.1323 (6). The dashed line in Fig. 4 shows that the L = 40 QMC data is indeed consistent with the asymptotic square root (13) in the range 0.03 T /J 0.07. However, because of finite-size effects, in this case we really need to go to L = 40 to recover the asymptotic behavior. Furthermore, we need to restrict to T 0.07 J since one can see from Fig. 4 that a pure power law is no longer a good description for T 0.07 J.
Comparison of (12) and (13) shows that one can reach exponentially small temperatures T f at h = ∆ by adiabatic magnetization of a spin-1 Heisenberg chain from an initial temperature T i at h = 0. For an initial temperature T i 0.1 J, the final temperature T f can be quantitatively estimated from the following combination of (12) and (13) with the parameters ∆, A, and B given above.

Scaling with spin quantum number s
We will now turn to the high-field region and consider in particular scaling with the spin quantum number s close to the saturation field.
It is straightforward to compute the excitation energy ε of a single flipped spin propagating above a ferromagnetically polarized background (see, e.g., [49]): The minimum of this one-magnon dispersion is located at k = π. The saturation field h sat. of the spin-s Heisenberg chain (1) is given by the one-magnon instability [50]. Accordingly, it can be determined by inserting the condition ε(π) = 0 into (15): The fact that the single-particle energy (15) scales linearly with s suggests to scale all energies close to the saturation field and at low temperatures linearly with s. In order to test this scenario, we have first computed the isentropes for a fixed chain length L = 20 and s ≤ 2. The s = 1/2 curves are obtained by a full diagonalization of the Hamiltonian [10] while we have performed additional ED computations for s = 1, 3/2, and 2 using the truncation procedure described in section 2. Fig. 5 shows the resulting isentropes for S/L = 0.05, 0.01, 0.15, 0.2, and 0.25 close to the saturation field (16). We observe that at low temperatures, scaling both h and T by s leads to a nice collapse of the isentropes around and in particular above the saturation field.
Further details can be read off from the temperature scans at a fixed magnetic field shown in Fig. 6. In this figure, we have chosen to show the specific heat C as a representative for energyrelated quantities 2 and the magnetic susceptibility χ as a representative for magnetic quantities. The top row of Fig. 6 corresponds to a magnetic field h = 1.1 h sat. . Here we observe a nice collapse with s for T /(J s) 0.15. The middle row of Fig. 6 corresponds to a magnetic field exactly equal to the saturation field. Here we observe a good scaling collapse for T /(J s) 0.04. Finally, the bottom row of Fig. 6 corresponds to a magnetic field h = 0.9 h sat. . In this case, one expects the scaling region to be pushed to even lower temperatures. However, for h < h sat. the finite size L = 20 leads to artifacts in the low-temperature behavior. Therefore it is difficult to make definite statements for this case.
To summarize this section, we have provided evidence of linear scaling with s for thermodynamic quantities of spin-s Heisenberg chains close to the saturation field. Note that this is very different from the quadratic scaling, e.g. of T with s 2 [23] or with s (s + 1) [51,52], which is needed to approach the classical limit.

Conclusions and outlook
In this paper we have illustrated the numerical computation of thermodynamic quantities and in particular the entropy for spin-s Heisenberg chains. From a technical point of view, we have described in section 2.1 how to perform a reliable computation of a large number of low-lying states using the Lanczos method and in section 2.2 how to compute the entropy directly by a QMC simulation.
In section 3 we have then focused on the s = 1 Heisenberg chain and shown in particular that one can cool with an adiabatic magnetization process during which the Haldane gap ∆ is closed. Many previous investigations (e.g., [9,10,14,18,19,23]) have focused on the saturation field. The reason for this is that the saturation field is a field-induced quantum phase transition at a known value of the magnetic field which also gives rise to technical simplifications. However, the scenario of quantum phase transitions [6,7] is universal and not restricted to the saturation field. Indeed, Fig. 4 is consistent with the same universal square-root behavior of the entropy S at h = ∆ in the spin-1 Heisenberg chain as observed previously for s = 1/2 chains exactly at the saturation field [4,10]. This may also be important from an experimental point of view since a possible spin gap may be accessible by laboratory magnetic fields even if the saturation field is out of reach. In fact, cooling by adiabatic magnetization when closing a spin gap has presumably been indirectly observed in pulse-field magnetization experiments on SrCu 2 (BO 3 ) 2 [53].
In section 4 we then moved to the saturation field and investigated scaling with the spin quantum number s. In contradistinction to the classical scaling regime where temperature should scale quadratically with s [23,51,52], there is a quantum scaling regime around the saturation field where one observes a collapse using the linearly scaled parameters h/s and T /s. Also in higher dimensions the one-magnon dispersion typically scales with s. By the same arguments as in section 4, we therefore expect linear scaling with s at any continuous transition to saturation. This expectation could be tested numerically with the methods of the present paper.