Dynamical Correlation Functions for Linear Spin Systems

Dynamical spin correlation functions in ( q,ω ) -space are calculated numerically for linear quantum-mechanical Heisenberg systems containing up to ten spins with s = 12 and s = 1 . We consider ferro- and antiferromag-nets including single-site and exchange anisotropies. We compare our results with various other theoretical treatments and draw conclusions for the dynamics of inﬁnite chains. In the case of the s = 1 planar Heisenberg ferromagnet direct comparison is made with inelastic neutron scattering cross sections of CsNiF 3 . The theoretical and experimental results are in good agreement.


Introduction
The one-dimensional magnetic chain has been the object of various theoretical treatments in the past (Bethe 1931, Lieb et al 1961, Des Cloizeaux and Pearson 1962, Fisher 1964, Bonner and Fisher 1964. Nowadays one-dimensional spin models are of more than only mathematical interest, since in recent years a number of magnetic compounds have been found which, under certain conditions, can be well described by the linear chain Heisenberg Hamiltonian with suitable anisotropy parameters. The one-dimensional properties are due to the arrangement of the magnetic ions parallel to some crystal axis, with intrachain exchange interaction being very much stronger than interchain magnetic coupling. A few representative substances are: (CH 3 ) 4 NMnCl 3 (TMMC) with s = 5 2 , CuCl 2 · 2N(C 5 D 5 ) (CPC) with s = 1 2 , and CsNiF 3 with s = 1. The static and dynamical properties of many such materials have been investigated extensively by various kinds of experiments. A comprehensive review where references to original work can be found was recently published by Steiner et al (1976). Here we are interested in the dynamics, i.e. in time-dependent quantities, especially in the dynamical two-spin correlation functions S α (i, t)S α (j, 0) . The latter are directly related to experimentally observable quantities like, for example, neutron scattering cross sections and NMR spin-lattice relaxation times. Actually, inelastic neutron scattering provides the most detailed information about two-spin correlation functions since the differential cross section involves the space-time Fourier transforms of the functions S α (i, t)S α (j, 0) . These experimental data have stimulated extensive theoretical work. Apart from some exact analytical results about spin chains (Bethe 1931, Lieb et al 1961, Des Cloizeaux and Pearson 1962, Fisher 1964) one can distinguish the following theoretical approaches to the dynamics of linear spin systems: (i) for classical systems, the use of moments (Lovesey and Meserve 1972, Lovesey 1974, Loveluck and Lovesey 1975; computer simulation (Blume et al 1975); (ii) for quantum systems, Green function techniques (Manson and Sjolander 1975); numerical diagonalisation of the Hamiltonian for finite chains followed by numerical evaluation of the appropriate correlation functions in the form of histograms Richards 1969, Richards andCarboni 1972).
This work uses the fourth method. Properly speaking it is an extension of the work of Richards and Carboni (1972) to various anisotropic systems and to spin quantum numbers s > 1 2 . Our motivation to do this work originates mainly from the recently available neutron scattering data for CsNiF 3 (Steiner et al 1975), which represents an anisotropic s = 1 Heisenberg ferromagnet. The experimental results of Steiner et al (1975) exhibit some quite particular features, different from those of isotropic Heisenberg systems. These facts have already been accounted for in the framework of classical (Loveluck and Lovesey 1975) and semiclassical (Villain 1974) calculations.
Our fully quantum-mechanical treatment allows for an interpretation from a different viewpoint and yields good agreement with experiment. A short account of our work has been presented previously (Müller and Beck 1977).
In Sec. 2 we introduce the Heisenberg Hamiltonian with various anisotropy parameters specifying the different models, as well as appropriate quantum numbers allowing us to simplify the computational procedure and to characterise the eigenstates. We define the dynamical spin correlation functions and their representation appropriate to a finite chain. Finally we discuss some selection rules for the matrix elements appearing in the correlation functions. The s = 1 2 Heisenberg antiferromagnet (Sec. 3.1) has already been extensively investigated with the same method (Richards and Carboni 1972). An exact theoretical result obtained by sum rules (Hohenberg and Brinkman 1974) is verified by extrapolation. Concerning the isotropic Heisenberg ferromagnet (Sec. 3.2) the dominant role of the spin-wave states and the spin-wave bound states for the correlation functions at low temperatures is described. Section 3.3 deals with the s = 1 Heisenberg ferromagnet including a local anisotropy appropriate to CsNiF 3 . We characterise those eigenstates which are most important for the correlation functions at low temperatures. On the basis of these few states we can understand the characteristic shape of the correlation functions. A quantitative comparison between experiment and our calculations is given for the peak position and the linewidth of the correlation function. Also the temperature dependence of the intensity is investigated. Section 3.4 gives an idea of the behaviour of the s = 1 2 Heisenberg ferromagnet with exchange anisotropy. Section 3.5 finally deals with the XY ferromagnet. Our numerical results are compared with those of exact analytical calculations for the correlation functions at T = 0 (Niemeijer 1967, McCoy et al 1971.

Notation and definitions
We treat the Hamiltonian aS z (l)S z (l + 1) + b S x (l)S x (l + 1) + S y (l)S y (l + 1) + c Depending on the sign of J the Hamiltonian (2.1) describes a ferromagnet or an antiferromagnet. In order to diagonalise H in the Hilbert space with dimension (2s + 1) N we make use of the following symmetries: (i) H is invariant with respect to lattice translations. Therefore each eigenfunction can be labelled by a quantum number k = 2πn/N (wavenumber) with n = 0, 1, . . . , N − 1. A translation by I lattice sites thus operates on the eigenfunctions as T l |k = e −ikl |k .
(2.2) (ii) (ii) Furthermore H is invariant under rotations of all spins around the z-axis. This allows us to introduce S T z , the z-component of the total spin, as a second quantum number.
(iii) In the first of the above-mentioned models the total spin of magnitude [S T (S T +1)] 1/2 commutes with H due to full isotropy in spin space. Thus S T is an additional quantum number.
By the use of these symmetries the original (2s + 1) N -dimensional matrix associated with H can be reduced into considerably smaller blocks. Thus, with N = 10 for s = 1 2 and N = 6 for s = 1 as maximum numbers of particles treated in this work, the largest subspace for which the diagonalisation had to be done numerically has dimension 26.
Once the eigenvalues E λ and eigenfunctions |λ are known, the Fourier transform of the timedependent spin correlation function can be determined by Z is the partition function and Instead of G xx and G yy it is more convenient to evaluate G +− and G −+ using For finite systems these functions are best represented, for fixed q, as histograms in frequency space. The number of non-vanishing terms in the double sum of equation (2.3) is greatly reduced by the following selection rules: (2.6c)

Results
3.1. s = 1 2 Heisenberg antiferromagnet This case has already been investigated with the same method by Carboni (1969, 1972). Their histograms, calculated for a chain of eight particles, served as a test for our computer program. At high temperatures the spectra G αα (q, ω) are essentially Gaussians centred at ω0. At low T , on the other hand, relatively sharp spin-wave peaks appear, which are produced by matrix elements between the ground state with wavenumber k 0 and the lowest excited states with wavenumbers k 0 + q. These spin-wave states were identified by Des Cloizeaux and Pearson (1962). In the thermodynamic limit their energies are shown to lie above the ground state energy by the amount In a finite system the dispersion relation is different from (3.1) as discussed by Des Cloizeaux and Pearson (1962). In contrast to the isotropic ferromagnet (see Sec. 3.2), even at T = 0 these lowest excited states do not completely make up the correlation function G αα (q, ω). Hohenberg and Brinkman (1974) investigated the spectral moments of G αα (q, ω) at T = 0 and small q, which can, as usual, be expressed in terms of static correlation functions. They showed the quantityω Thus, the moment ratio (3.3) lies somewhat higher than the threshold value determined by (3.1). We have evaluated the quantity (3.3) for various values of N and the corresponding smallest non-zero wavenumber q min = 2π/N . With increasing N the values of w;a(q) indeed approach the theoretical value (3.4) as is shown in figure 1. 3.2. s = 1 2 Heisenberg ferromagnet Here, in contrast to the antiferromagnet, the correlation function G αα (q, ω) at T = 0 consists for all q only of a single sharp spin-wave peak at the energy SW (q) = 2Js(1 − cos q). (3.5) The only non-zero matrix elements contributing to G αα (q, ω) at T = 0 are those involving one of the degenerate ground states (k = 0, −N s ≤ S T z ≤ N s) and the corresponding member of the group of degenerate spin-wave states with k = q and the appropriate quantum number S T z . At finite T there are additional contributions to G αα (q, ω). We find that in our small systems the most important ones at low T arise from matrix elements between a spin-wave state with wavenumber k and a spin-wave bound state with wavenumber k + q and, more generally, between spin-wave bound states with wavenumbers differing by q. These bound states can be characterised by the number r of reversed spins, givenby r = N s − S T . For s = 1 2 their energies are approximately given by (Jain et al 1975): From (3.5) and (3.6) one can easily derive the following relation, which is approximately valid for small q: i.e. the binding energy goes to zero as k → 0. Therefore, several of the above-mentioned matrix elements, namely those involving a ground state (k = 0) and a spin-wave state (k = q), a spinwave state (k = q) and a spin-wave bound state (k = 2q) as well as a pair of spin-wave bound states (k = nq) and (k = (n + 1)q), all contribute approximately at the same frequency SW (q). The relevance of these statements for an infinite system is not quite clear, however, since in that case there is, for all k, a continuum of states above a threshold energy which is lower than the one-spin-wave energy (3.5).
For systems with spin quantum numbers s > 1 2 the dispersion relation (3.5) is still valid. Again there exist spin-wave bound states. We have evaluated G αα (q, ω) for systems with N = 6 spins of quantum number s = 1 and have found qualitatively the same frequency and temperature behaviour as for s = 1 2 systems.

s = 1 Heisenberg ferromagnet
The model described in this section is still isotropic in the exchange term of the Hamiltonian (2.1) but there is a local anisotropy (c > 0), which makes it favourable for a spin to lie in the XY plane. (For s = 1 2 such an anisotropy term is a c-number). According to Steiner et al (1975Steiner et al ( , 1976 this model is appropriate for CsNiF 3 , which behaves as a linear magnetic system with an easy plane in spin space. Now G zz (out-of-plane correlation function: OPC) and G xx = G yy (in-plane correlation function: IPC) are no longer identical. Experimentally one finds the following facts: (i) The IPC has a fairly broad peak the width of which is comparable with that of isotropic systems, concerning magnitude and temperature dependence.
(ii) The OPC has a narrow peak (its width is essentially determined by experimental resolution) the intensity of which decreases rapidly with rising T . Loveluck and Lovesey (1975) investigated the dynamical properties of classical easy-plane chains by evaluating the correlation functions in a continued-fraction representation, using exact expressions for the static correlation functions. Villain (1974) gave an approximate quantum-mechanical treatment for one-dimensional easy-plane ferro-and antiferromagnets, by making use of a semipolar representation of the spin operators. In both approaches one finds the following dispersion relation for spin-wave excitations ω 2 (q) = 4J 2 s 2 (1 − cos q)(1 − cos q + c/J). (3.8) Moreover the correlation functions calculated by both of the above methods agree well with the experimental findings (i) and (ii). We calculated eigenvalues and eigenfunctions as well as IPC and OPC for a planar Heisenberg ferromagnetic s = 1 chain of six particles with a ratio c/J = 0.21, appropriate for CsNiF 3 (Steiner et al 1975(Steiner et al , 1976. For classifying the eigenstates of the anisotropic model it is useful to start from the eigenstates of the isotropic system and to study the change when the anisotropy is switched on. As described in Sec. 3.2, the low-lying eigenstates for an isotropic Heisenberg ferromagnet (which are important for G αα at small T ) can be grouped as follows: (i) 2N s + 1 degenerate ground states with k = 0 and S T = N s; (ii) for each k = 0 there are 2N s − 1 degenerate spin-wave states with S T = N s − 1; (iii) spin-wave bound states of r reversed spins with S T = N s − r and degeneracy 2S T + 1.
The anisotropy of the planar model partially removes the degeneracies (only the degeneracy with respect to the sign of S T z persists). Figure 2 shows states of the above three categories for the isotropic and the planar s = 1 Heisenberg ferromagnet. The correlation functions G zz (OPC) and G xx (IPC) for a planar system with N = 6 are shown in figure 3. Although histograms for such small systems can only give a qualitative picture for the behaviour of a large system, we can make the following statements: (i) the IPC is relatively broad and only weakly T -dependent; (ii) the OPC is fairly sharp; its intensity decreases with increasing T , while the width increases.  In-plane (Gxx) and out-of-plane (Gzz) correlation function for the planar Heisenberg ferromagnetic s = 1 chain of six particles. The value c = 0.212J for the anisotropy parameter is appropriate for CsNiF3 and q is close to qz = 0.35π used in neutron scattering (Steiner et al 1975). The three temperatures correspond to those of the references: (a) T = 0.208J, (b) T = 0.343J, (c) T = 0.5J. This is in qualitative agreement, both with the experimental data and the previous theoretical work. Furthermore the peak positions for both functions lie very close to the energies given by (3.8).
The different behaviour of G zz and G xx is easy to understand on the basis of figure 2 and the selection rules (2.6b,c). For G zz : S z (q) connects states with equal S T z . These states, however, are affected by the anisotropy in a similar way, i.e. the energy difference between states with the same value of S T z , belonging to different 'multiplets' are almost identical. This leads to a narrow peak in G zz . For G xx : S + (q) and S − (q) have matrix elements between states the quantum numbers S T z of which differ by ±1, i.e. between states which are shifted differently by the anisotropy. Therefore the peak of G xx is broader. Figure 4. T -dependence of the peak position ωp of IPC. Experimental data (•) are by Steiner et al (1975) for q = 0.35π, the numerical results for systems of five (broken curve) and six (full curve) particles were evaluated at the corresponding wavenumber q = (2/N )π and transformed to the value q = 0.35π according to (3.8).  Steiner et al (1975) for q = 0.35π, corrected for finite resolution. The numerical results for systems of five (broken curve) and six (full curve) particles were evaluated at the corresponding wavenumbers q = (2/N )π. Figures 4 and 5 give a quantitative comparison between experiment and our theoretical results for the temperature dependence of peak position ω p and linewidth δ of the IPC, defined as (3.9) The experimental width in figure 5 is corrected for finite resolution. In figure 6, finally, we show the temperature dependence of the intensity I defined by for the isotropic and the planar case. The decrease of the OPC intensity with rising T is in qualitative agreement with experimental data (Steiner et al 1975). Figure 6. T -dependence of the intensity I at q = 2π/N (N = 6) for the isotropic s = 1 Heisenberg ferromagnet (full curve) as well as for IPC (broken curve) and OPC (chain curve) of the planar s = 1 Heisenberg ferromagnet appropriate for CsNiF3.
At first sight the behaviour of the IPC and OPC intensities shown in figure 6 seems surprising: the system is expected to order in the easy plane at T = 0. Thus the in-plane fluctuations should be critically enhanced when T → 0. Table 1 shows that this is indeed the case for q = 0, whereas for q = π/3 the OPC fluctuations dominate. Such a behaviour is also obtained when the correlation functions are approximated by an Ornstein-Zernicke form. Moreover it is also characteristic of another planar model, the XY chain. Table 1. Temperature dependence of IPC (Ix) and OPC (Iz) intensity for q = 0 and q = π/3. T /J I z (0) I x (0) I z (π/3) I x (π/3)

Heisenberg model with exchange anisotropy
Here we study the Heisenberg Hamiltonian with parameters 0 < a < 1, b = 1, c = 0 (the case a =0 is treated separately). The effect of an exchange anisotropy (a < 1) on the eigenvalues of a ferromagnetic chain is essentially the same as that of a local anisotropy (c > 0), described in detail in Sec. 3.3 (for s = 1 2 systems actually only the first kind of anisotropy exists). It is therefore obvious that the correlation functions G zz and G xx = G yy behave in the same way as OPC and IPC for the planar model. In figure 7 we present the linewidth δ at q = 2π/N (N = 8) for an isotropic (a = 1) and an anisotropic (a = 0. 75) chain. Figure 7. T -dependence of the linewidth δ at q = 2π/N (N = 8) for the s = 1 2 Heisenberg ferromagnet (full curve) corresponds to Gαα for the isotropic system (a = 1), (chain curve) and (broken curve) correspond to Gzz and Gxx, respectively for the anisotropic system (a = 0 75). Lieb et al (1961) have given a mapping of the XY chain onto a system of non-interacting fermions, the energies of which are given by Λ(k) = J cos k.

s = 1 2 XY ferromagnet
(3.11) There are 2N possible states in the band; in the ground state all negative energy levels are occupied. The number N F of fermions present in a given eigenstate is related to the z-component of the total spin by (3.12) Due to the cyclic boundary conditions for the spin chain, the fermion Hamiltonian has a different form in the subspace of even N F and odd N F , respectively. The evaluation of the dynamic spin correlation functions is still fairly complicated, since the transformation from spin to fermion operators is non-linear. Niemeijer (1967) calculated the correlation functions G zz (q, ω) for all temperatures, whereas McCoy et al (1971) gave results for G xx = G yy at T = 0. We have calculated these correlation functions for systems of N = 8 spins in order to determine to what extent it is possible to draw conclusions concerning the behaviour of infinite chains from the results for such small systems. G zz (q, ω) and G xx (q, ω) are shown in figure 8 for various values of q at T = 0.
The results for G zz can be understood in the following way: The ground state has N F = 4 fermions occupying the negative energy levels, and therefore S T z = 0. For G zz matrix elements between the ground state and other states with S T z = 0 contribute. These states also belong to the subspace with N F = 4 representing particle-hole, two particle-two hole etc excitations above the ground state. According to our results only one particle-one hole excitations contribute.
From these facts it is straightforward to find the lower and upper bounds for the spectrum of G zz (k, ω) for N → ∞. The possible energies of a particle-hole state with wavenumber k can be written as where k = (2π/N )n, α = m p /n. For N → ∞, n → ∞ α lies between 0 and 1. Equation (3.18) implies Thus the frequency spectrum of G zz (k, ω) lies between J sin k and 2J sin(k/2). The same result can be found from the analytic form of G zz (k, ω) given by Niemeijer (1967). G xx is more complicated even at T = 0: there are small contributions from excitations with higher numbers of particles and holes.

Summary
The exact eigenfunctions and energy eigenvalues of linear Heisenberg systems containing up to ten particles were obtained by numerical diagonalisation of the Hamiltonian. Using these results the dynamical two-spin correlation functions in (q, ω)-space were calculated as functions of ω in the form of histograms at various q and T . We have demonstrated that conclusions can be drawn for the dynamics of infinite chains from our results of small systems. Direct comparison was made with available neutron scattering data as well as with former exact or approximate, classical or quantum mechanical theoretical treatments.
Our main results are: (i) s = 1 2 Heisenberg antiferromagnet: The spin-wave states identified by Des Cloizeaux and Pearson (1962) play a dominant part in the correlation functions at low T , leading to a relatively sharp peak at the spin-wave frequency (3.1). However, even at T = 0, there are contributions at higher frequencies, as is also verified by the exact moment calculation of Hohenberg and Brinkman (1974).
(ii) s = 1 2 Heisenberg ferromagnet: At T = 0 only spin-wave excitations are involved in the correlation functions. At finite but low T and small q the correlation functions are determined by excitations of spin-waves and spin-wave bound states, all contributing near the spin-wave frequency (3.5). A similar behaviour is displayed by s = 1 systems.
(iii) s = 1 planar Heisenberg ferromagnet: Correlation functions were calculated for systems of five and six particles. Comparison with inelastic neutron scattering cross sections measured by Steiner et al (1975) gives good agreement in: (a) the qualitative shape of both IPC and OPC; the T -dependence of peak position and intensity of IPC and OPC; the T -dependence of the line width of IPC. The difference in shape of the OPC and IPC can be understood on the basis of the effect of anisotropy on the low-lying eigenstates.
(iv) Heisenberg model with exchange anisotropy: From finite-chain calculations it is evident that an anisotropy in the exchange interaction (a < 1) of the Heisenberg Hamiltonian (2.1) has a similar influence on the correlation functions as a local anisotropy (c > 0) comparable in magnitude.
(v) s = 1 2 XY ferromagnet: The spin excitations occurring in the correlation functions can be identified with particle-hole excitations in a non-interacting-fermion representation (Lieb et al 1961). From the calculations for systems of eight particles we have drawn conclusions for infinite systems, which are in good agreement with exact theoretical results (Niemeijer 1967).
The present calculations are restricted to correlation functions in (q, ω)-space, being directly related to neutron scattering data. By the same method one can evaluate autocorrelation and paircorrelation functions, which enter the expression for the NMR spin-lattice relaxation time. Indeed many such experiments on quasi-one-dimensional magnetic systems have been performed. They have partly been explained by several (mainly classical) theoretical approaches. A later publication will be devoted to this aspect of magnetic linear chains, starting from a quantum-mechanical point of view.