Recent Advances in the Calculation of Dynamical Correlation Functions

We review various theoretical methods that have been used in recent years to calculate dynamical correlation functions of many-body systems. Time-dependent correlation functions and their associated frequency spectral densities are the quantities of interest, for they play a central role in both the theoretical and experimental understanding of dynamic properties. The calculation of the relaxation function is rather difficult in most cases of interest, except for a few examples where exact analytic expressions are allowed. For most of systems of interest approximation schemes must be used. The method of recurrence relation has, at its foundation, the solution of Heisenberg equation of motion of an operator in a many-body interacting system. Insights have been gained from theorems that were discovered with that method. For instance, the absence of pure exponential behavior for the relaxation functions of any Hamiltonian system. The method of recurrence relations was used in quantum systems such as dense electron gas, transverse Ising model, Heisenberg model, XY model, Heisenberg model with Dzyaloshinskii-Moriya interactions, as well as classical harmonic oscillator chains. Effects of disorder were considered in some of those systems. In the cases where analytical solutions were not feasible, approximation schemes were used, but are highly model-dependent. Another important approach is the numerically exact diagonalization method. It is used in finite-sized systems, which sometimes provides very reliable information of the dynamics at the infinite-size limit. In this work, we discuss the most relevant applications of the method of recurrence relations and numerical calculations based on exact diagonalizations.


Introduction
Dynamical correlation functions are central to the understanding of time-dependent properties of many-body systems. They appear ubiquitously in the formulation of the fluctuation-dissipation theory, where the response of a system to a weak external perturbation is cast in terms of a time-dependent relaxation function of the unperturbed system [1,2].
In this article, we are concerned with the recent calculations of such correlation functions. We shall cover two lines of approach, namely the method of recurrence relations and the method of exact diagonalization.
The method of recurrence relations was developed in the early 1980s [3,4,5,6,7] following the ideas of the Mori-Zwanzig projection operator formalism [8,9]. Essentially one solves the Heisenberg equation of motion for an operator of an interacting system, from which one obtains dynamic correlation functions, a generalized Langevin equation, memory functions, etc. Review articles found in the literature cover the earlier developments [10,11,12].
On the other hand, exact diagonalization methods have also been used in several areas of physics [13,14,15,16,17]. In this method one numerically determines the eigenvalues and eigenfunctions of a given Hamiltonian of a finite system to find the dynamical correlations of interest. The main drawback is that one is bound by computer limitations and must deal with finite systems. In addition, being a numerical method, it does not provide any new general insight in the form of theorems, etc. Nevertheless, one can obtain surprisingly good results which can be readily extended to the thermodynamic limit. In a way, exact diagonalization complements the method of recurrence relations, especially when solutions become hard to obtain by analytic means. Other approaches can be found in Refs. [18,19,20,21,22,23,24].

Dynamical correlation functions
Consider a system of N elements such as particles, spins, etc., governed by a time-independent Hamiltonian H, in thermal equilibrium with a heat bath at temperature T . For two dynamical variables X and Y of the system, the time-dependent correlation function is given by the average: where Tr[. . . ] denotes a trace over a complete set of states. Here, β = 1/k B T is the inverse temperature, Z ≡ Tr exp(−βH) is the canonical partition function, and X(t) is a time-dependent operator in Heisenberg representation X(t) = exp(iHt/h)X exp(−iHt/h), which satisfies: where [X(t), H] is the quantum commutator. In a classical system, the operators are replaced by classical dynamic variables, the trace by integral over the phase space, and the commutators by Poisson brackets. For a given variable, the time-dependent correlation function C(t) reads: Its Fourier transform S(ω) is called the spectral density, or frequency spectrum: If we use the integral representation of the Dirac δ-function: then we obtain Since the Hamiltonian is time-independent, it follows that C(t) in Eq. (3), has the property < X(0)X(t) > = < X(τ )X(t + τ ) >. If we take τ = −t, then < X(0)X(t) > = < X(−t)X(0) >. Also, it follows that S(ω) is real. Due to the invariance of the trace under cyclic permutations, one finds that S(−ω) = exp(−βhω)S(ω). In the classical limit (h = 0 ) or at infinite temperature β = 1/k B T = 0, it follows that S(ω) is even in ω. In general, the asymmetry in S(ω) is a typical quantum feature, and is referred to as the detailed balance. Dynamical correlation functions appear in the relaxation function R(t) from linear response theory [25,2]: where < · · · > is a canonical average, and X and Y are operators. Time-dependent correlation functions appear in the dynamical structure factor, are related to the inelastic neutron-scattering cross section, where the neutron energy changes upon the scattering process. For a system of interacting spins on a lattice, the dynamic structure factor reads: where S are spin variables and the sum runs over all the lattice sites.
In light scattering experiments, the scattered intensity is given by the differential cross section, proportional to: where the form of operator A is system dependent. It also depends on the the particular frequency of the incoming light.

The method of recurrence relations
The time evolution of a Hermitian operator A(t) is governed by the Heisenberg equation: where: Consider a time-independent and Hermitian Hamiltonian H. From now on we will be using a system of units in whichh = 1. We seek a solution to Eq. (10) for t ≥ 0, thus we set A(t) = 0 for t < 0. In the method of recurrence relations, the formal solution: is cast as an orthogonal expansion in a realized Hilbert space S of d dimensions.
That Hilbert space S is realized by the scalar product: where X, Y ⊂ S, β is the inverse temperature, X(λ) = exp(λH)X exp(−λH), and < · · · > denotes a canonical ensemble average. Thus, the time evolution of A(t) is written as: where {f ν } is a complete set of states in S, while the time-dependence is carried out by the coefficients a ν (t). The dimensionality d of the realized Hilbert space S is still unknown, but it will be determined later. If d turns out to be finite, the solutions are oscillatory functions. However, in most interesting cases d is infinite. The method of recurrence relations imposes constraints on which type of solutions are admissible.
By choosing the basal vector f 0 = A(0) = A, the remaining basis vectors are obtained following the Gram-Schmidt orthogonalization procedure, which is equivalent to the recurrence relation: with f 1 ≡ 0, ∆ 0 ≡ 0. The quantity ∆ ν is defined as the ratio between the norms of consecutive basis vectors: The ∆'s are referred to as the recurrants whereas Eq. (15) is termed the first recurrence relation, or RRI. The time-dependent correlation function C(t) is given by: The basal coefficient a 0 (t) is just the time-dependent correlation function. The time-dependent coefficients a ν (t) obey a second recurrence relation (RRII): whereȧ ν (t) = da ν (t)/dt, and a −1 ≡ 0. It follows from Eq. (14) that the initial choice f 0 = A(0) implies a 0 = 1 and a ν (0) = 0 for ν ≥ 1. Thus the complete time evolution of A(t) is obtained by the two recurrence relations RRI and RRII. One should note that only in very few cases a closed analytic solution to a model can be found. More often, as in many-body problems, approximations are required. A generalized Langevin equation can be derived for A(t) [3,4,26]: where φ is the memory function and F [t] the random force. The random force is given as an expansion in the subspace of S: where the coefficients b ν satisfy the convolution equation: The memory function φ(t) is φ(t) = ∆ 1 b 1 (t). The remaining b ν 's, b 2 , b 3 , are the second memory function, the third memory function, . . . , etc.
Consider now the Lapace transform a ν (z) of a ν (t): Then RRII can be transformed in the following way: These equations can be solved for a 0 (z): resulting in a continued fraction. As can be seen from Eq. (14) and the recurrence relation RRII, that the time-dependence actually depends on the recurrants ∆ ν only. Therefore, the knowledge of all recurrants provides the necessary means to obtain the time correlation function. Moreover, the structure of RRII must be obeyed by time correlation functions. Thus, a pure exponential decay as well as special polynomials can be ruled out as solutions, since their recursion relations are not congruent to RRII. Also, from RRII one obtains (da 0 (t)/dt)| 0 = 0, which precludes a pure time exponential as well as other functions that do not have zero derivative at t = 0. The method of recurrence relations have since been applied to a variety of problems, such as the electron gas [27,28,30,29], harmonic oscillator chains [31,32,33,34,35,36,37,38,39,40], many-particle systems [44,43,41,42], spin chains [45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60], plasmonic Dirac systems [61,62], etc.

The method of exact diagonalization
Given a system governed by a Hamiltonian H, one wishes to numerically determine the time correlation function C(t), defined by: where A(t) = exp(iHt)A exp(−iHt),h = 1, and the brackets denote canonical averages. We consider here self-adjoint operators A and the Hamiltonian H. One numerically diagonalizes H and then uses its eigenvalues E n and eigenvectors |n >, H|n >= E n |n >, to calculate C(t) in Eq. (26), where: with the partition function Z = n exp(−βE n ). Notice that the time correlation function is normalized to unity at t = 0, that is, C(0) = 1.
Another quantity of interest is the moment µ k , also referred to as the frequency moments which can be obtained from the Taylor expansion of C(t) about t = 0: Since C(0) = 1 , it follows that µ 0 = 1. The moments are given by: where L is the Liouville operator, Eq. (11).
From the moments, one can use conversion formulas to obtain the recurrants ∆'s of the method of recurrence relations from the frequency moments [11]. Suppose the moments µ 0 = 1 and µ 2k , k = 1, . . . , K are known. The first K recurrants ∆ ν are determined by the equations: for k = ν, ν + 1, . . . , K and ν = 1, 2, . . . , K, with µ = 0.s For instance, if the first moments µ 0 = 1, µ 2 , µ 4 , . . . , µ 10 are given, the recurrences are obtained from Eq. (31): Conversely, suppose one has the first K known recurrants, ∆ ν , ν = 1, 2, . . . , K, and ∆ −1 = ∆ 0 = 1. Then, the moments µ (0) 2ν = µ 2ν are obtained from the following conversion formula: for ν = k, k − 1, . . . , 1 and k = 1, 2, . . . , K, with µ In case the first recurrants ∆ 1 , ∆ 2 , . . . , ∆ 4 , are known, the moments µ are found to be: Typically, the analytical determination of the recurrants becomes increasingly time consuming. In practice, only a few of them can be obtained to be used in an extrapolation scheme to obtain higher-order recurrants. Several extrapolation schemes have been used. One of the simplest is to set the unknown recurrants to zero, thus truncating the continued fraction for a 0 (z), which leads to a finite number of poles in the complex plane [18]. In other problems, it is most appropriate to introduce a Gaussian termination, that is, a sequence of recurrants ∆ν that grow linearly with its index ν, ∆ ν = ν∆ [45,11]. Other extrapolation schemes are tailored to the problem at hand, especially if the recurrants are not expected to grow indefinitely.

Applications to interacting systems
The dynamics of spin chains has attracted a great deal of attention in recent decades. Exact results for the longitudinal dynamics of the one dimensional XY model have been obtained with the Jordan-Wigner transformation [63]. Later, exact results for the transverse time correlation functions of the XY and the transverse Ising chain were obtained at the high temperature limit by using different methods [64,65,45].
A great deal of progress was achieved in the calculations of the dynamic correlation functions of spin models in one dimension. It was soon recognized that exact solutions using the method of recurrence relations were difficult to obtain, however a notable exception is the classical harmonic oscillator chain where the time correlation functions were obtained exactly [31].
The problem of a mass impurity in the harmonic chain was solved later, and its dynamical correlation functions were found to have the same form as in the quantum electron gas in two dimensions, thus showing that unrelated quantities in these two models displayed the same dynamical behavior, that is, the have dynamic equivalence [66]. It should be mentioned that harmonic oscillator chains have been the subject of a considerable amount of work with the method of recurrence relations [32,33,34,35,36,37,38,39,40].
The method of recurrence relations provides important insights on how to proceed to obtain reliable approximate solutions. The cornerstone quantity in the dynamics is the recurrant, which is the only quantity that ultimately determines the dynamics of the model. Often it is only possible to determine a few of the recurrants analytically. The calculations become too lenghty so that one must stop at a given order. Thus, an extrapolation method must be devised for the higher order recurrants, which hopefully will have the essential ingredients to produce reliable time-dependent correlation functions for longer times as well as spectral densities with the expected behavior near the origin ω = 0 [56,67,68,43].
The dynamics of the transverse Ising model in two dimensions was studied with the method of recurrence relations [69,70,72,71]. The dynamic structure factor of that model compares well to the experimental data of the compound LiTbF 4 [73].
The dynamics of spin ladders has also attracted interest from researchers. The dynamical correlation functions were obtained for a two-leg spin ladder with XY interaction along each leg and interchain Ising couplings in a random magnetic field. More recently, the dynamics of a ladder with Ising couplings in the legs and steps as well as four-spin plaquette interactions in a magnetic field [74] have been also investigated.
The dynamical correlations of the Heisenberg model in one dimension have been have been a subject of great interest in the recent decades [13,14,15,21,75]. The method of recurrence relations has been employed in various works [77,76,78,79,80,81,82]. In spite of the progress made thus far, the long-time dynamics of the Heisenberg spin model is still an open problem. For instance, there is the standing problem on the power law exponent α > 0 ( ∼ t −α ) of the time correlation function as t → ∞. From the work of Fabricius et al. [15], we find that the time correlation functions of the Heisenberg model decay more slowly than that in the XY model, for which the exact solution is known, C(t) ∼ t −1 for large t. Thus we infer that the numerical evidence suggests that α ≥ 1 for the Heisenberg model.
There has been a great deal of work that uses exact diagonalization to study the dynamics of spin systems [13,14]. Earlier works with the Heisenberg model used this technique. Later on, other systems were scrutinized by using exact diagonalization. One of those systems is the Ising model with four-spin interactions in a transverse field. The time correlation function was obtained for one dimensional and infinite temperature [83], where the Gaussian behavior shown in the usual transverse Ising model was ruled out. The effects of disorder on the dynamics of that model were obtained for the cases where the random variables are drawn from bimodal distributions of random couplings and fields [16,84,85]. Dynamical correlation functions were also obtained for the system at finite temperatures, ranging from T = 0 to T = ∞ [86].

Heisenberg model with Dzyaloshinskii-Moriya interactions
The dynamical structure factor for a quantum spin Heisenberg chain with Dzyaloshinskii-Moriya (DM) interactions [87,88] has been investigated by different approaches, such as spin wave theory [89], mean-field [90], and projection operator techniques [91]. The dynamics of the related XY model with DM interactions was also studied by employing Jordan-Wigner fermions [92,93,94].
The dynamical correlation functions of the spin-1/2 Heisenberg model with DM interactions in a transverse magnetic field was studied recently with the method of recurrence relations. The model Hamiltonian for a one-dimensional chain is given by: where J is the Heisenberg coupling, D is the Dzyaloshinskii-Moriya interaction, and B i is a magnetic field perpendicular to the DM axis. The quantities σ x,y,z i are the usual Pauli operators. The effects of a uniform magnetic field B j = B on the dynamics are investigated in the infinite temperature limit [95]. The purpose is to determine the time correlation function C(t) =< σ z j σ z j (t) > and its associated spectral density S(ω). The first four recurrants are determined analytically and an extrapolation scheme is devised to obtain higher order recurrants. Such scheme must take into account what is already known from the solutions of related problems.
One crucial point is to determine whether or not the extrapolated recurrants grow indefinitely. The time correlation function of the longitudinal spin component in the XY chain is known exactly at T = ∞, C(t) = J 2 0 (t) ∼ t −1 asymptotically for large times, where J 0 is the Bessel function of first kind [63]. In this case, the recurrants tend to a constant finite value as ν → ∞. There are numerical indications that the time correlation function of the Heisenberg model decays as a power law [15], which suggests that the extrapolated recurrants grow asymptotically to a finite value. For the Hamiltonian Eq. (35) the extrapolation is the following power-law: where n c is the order of the last exactly calculated recurrant. The limit value ∆ ∞ is obtained by extrapolating the last two recurrants to the origin of 1/ν. The constants b and β guarantee smooth behavior of the recurrants above and below ν = n c . Once the recurrants are obtained, the relaxation function and its spectral density can be readily obtained. For the special case without DM interaction, the result shows good agreement with the known results for the XY and Heisenberg models. The full calculation reveals that the effects of the external field are to produce stronger and more rapid oscillations in the relaxation functions, as well as a suppression of the central peak in the spectral density. In addition a peak centered at a well defined frequency appears, which is attributed to an enhancement of the collective mode of spins precessing about the external field. It should be noted that the method of recurrence relations was also used to study the dynamics of the XY model with DM interaction [59].
The effects of disorder in a transverse magnetic field on the dynamical correlation functions are investigated with the bimodal distribution for B i : The method of recurrence relations is then applied to obtain the dynamical correlation functions for a given realization of disorder [60]. Next, the average over the random fields is performed by using the distribution Eq. (37). This is accomplished by defining the scalar product in the Hilbert space S, so as to include an average over the random variables in addition to the thermal average. Four recurrants ∆ ν are obtained, and an extrapolation is made for the remaining recurrants. In practice only the first dozen are needed to attain convergence. The time correlation function and its associated spectral density were obtained for D = 1 and B A = 0 and B B = 4 in units of the Heisenberg coupling J. When the probability p is very small, a strong central mode appears, as well as a shoulder in the spectral density. As p increases, there is supression of the central mode as well as the shoulder. On the other hand, for large values of the probability p, a nonzero frequency peak appears, resulting from the precession of the spins around the magnetic field adding further suppression of the central peak. This central mode behavior versus collective dynamics, is a known feature of the dynamics of spin systems and they are in some sense universal. However, in the present case the appearence of a shoulder for small p is an interesting novel feature.

Random transverse Ising model
Consider the s = 1/2 spin model in one dimension: where J i and B i are exchange couplings and transverse fields, respectively. These couplings and fields are random variables drawn from distribution functions. The quantities σ α i (α = x, y, z) are Pauli matrices. The model is referred to as the random transverse Ising model (RTIM), and its dynamical correlation function in the infinite temperature limit has been investigated by using the method of recurrence relations [96].
The time correlation C(t) is defined by: where the line indicates that an average over the random variables is performed after the statistical average < · · · >. The time evolution of σ z j (t) in a system governed by the Hamiltonian Eq. (38) is given as an expansion in a Hilbert space S of d dimensions, where d is to be determined later: where f ν are orthogonal vectors spanning S. The time dependence is contained in the coefficients a ν (t).
The inner product in S in the infinite temperature limit is defined in such a way that it encompasses both the thermal average in a realization of disorder and the average over the random variables: where A and B are vectors in S. This definition of scalar product ensures that the form of the recurrence relations in unchanged.
The zeroth basis vector f 0 is chosen as the variable of interest, f 0 = σ x j . Thus, the zeroth-order coefficient a 0 (t) can be identified with the time-dependent correlation function of interest: The remaining basis vectors f ν , ν = 1, 2, ..., d − 1, are obtained from the recurrence relation RRI, Eq. (15). The first vectors are then: etc. The vectors f 4 , f 5 , . . . , f 9 were obtained analytically but not reported because of their length [96]. However, they were used in all of the subsequent calculations. The first three recurrants are the following, Notice that the couplings and fields are site-dependent.
There are two types of disorder considered in Ref. [96], random fields and random spin couplings. Each case is treated separately. In both cases a simple bimodal distribution is used for the random variable. The field B i (or the coupling J i ) can assume two distinct values, with probalities q (p) and 1 − q (1 − p), respectively. The time correlaton function and the spectral density are then obtained numerically. For the pure cases, (p = q = 1), two types of behavior emerge, depending on the relative strength between J and B. For J > B, the dynamics is dominated by a central-mode behavior, whereas for J < B a collective-mode is the prevailing dynamics. In the disordered cases, the dynamics is neither central-mode nor collective-mode type, but something in between those types of dynamics.

Transverse Ising model with next-to-nearest neighbors interactions
Consider the transverse Ising model with an additional axial next-nearest-neighbor interaction (transverse ANNNI model) [17]. The Hamiltonian for a chain with L spins can be written as: where σ α i are the usual spin-1/2 operators, α = x, y, z. Periodic boundary conditions are imposed on this model, namely σ α i+L = σ α i . Consider antiferromagnetic (J 1 > 0) Ising interactions. A competing ferromagnetic interaction is assumed for the next-nearest-interaction (J 2 > 0). The transverse magnetic field (B) induces the quantum fluctuations. In what follows we set J 1 = 1 as the unity of energy.
In the absence of a transverse magnetic field and of thermal fluctuations (T = 0) the ground-state properties of the model are exactly soluble and several phases are present. For J 2 < 0.5 the ground state is ordered ferromagnetically. For J 2 > 0.5, a phase consisting of two up-spins followed by two down-spins is periodically formed. The phase is known as < 2, 2 >-phase or an anti-phase. For J 2 = 0.5, the model has a multiphase point. The ground-state is highly degenerate with many phases of the type < p, q > corresponding to a periodic phase with p-up spins followed by q-down spins, among other spin configurations. The number of degenerate states increases exponentially with the size of the system L. In the case where J 2 = 0 and the magnetic field is switched on, the model becomes the Ising model in a transverse field which was exactly solved by Pfeuty [97]. In this model, due to quantum fluctuations induced by the transverse magnetic field, a second order phase transition occurs at B = 1, which separates a ferromagnetic phase at low magnetic fields from a paramagnetic phase at high magnetic fields. For the full Hamiltonian, Eq. (45), the competing interaction between the ferromagnetic and antiferromagnetic terms induces frustration in the magnetic ordering. This will give rise to a much richer variety of phases when either the transverse magnetic field or the spin-spin interations are varied, such as ferromagnetic or antiferromagnetic phases, disordered or paramagnetic phases, and floating phases [17]. Such variety of phases in the ground-state could carry over their effects into the dynamics at the high temperature limit, like the known transverse Ising model. In this model, a signal of the ground-state transition is manifested in the Gaussian behavior at criticality of the dynamical correlation functions at T = ∞ [45].
The main quantity of interest is the time-dependent correlation function: where σ x j (t) = e iHt σ x j e −iHt and < O > is a canonical average of the operator O. The method of exact diagonalization will be employed to study the dynamics, however, the recurrants of the method of recurrence relations will also be obtained.
The numerical calculations will be performed at the high-temperature limit, T = ∞, hence: One of the properties of C(t) is that it is real and an even function of the time t. Therefore, the Taylor expansion about t = 0 has only even powers of t: where the frequency moments are expressed in terms of the trace over iterated commutators: with L defined such that: where H is the Hamiltonian and A an operator. The correlation function is calculated in the Lehman representation. First, we consider the energies E n and eigenstates |n > of the Hamiltonian, obtained from the eigenvalue equation H|n >= E n |n >. Then, the correlation function takes the form: where the moments µ 2k are given by: The spectral density S(ω) is simply the Fourier transform of C(t): After using Eq. (51), the spectral density can be cast in the form: where ǫ nm ≡ E n − E m . The Dirac δ-function is approximated by a rectangular window of width a and unit area, centered at the zeros of their arguments. The width a, can be adjusted to reduce fluctuations. Another approach could be the use of histograms, such as in Ref. [86]. However, the general shape of the spectral density S(ω) is the same, although the rectangle approximation gives more accurate results. Therefore, both dynamical correlation functions C(t) and S(ω) can be calculated directly via exact diagonalization.
As a case test, Guimarães et al. [17] consider B = 1 and J 2 = 0, the usual transverse Ising model (TIM) with dynamical correlation functions known exactly in the high temperature limit [65,45]. Figure 1 shows their numerical Convergence toward the thermodynamic result increases as the system size grows. However, already for L = 13 the numerical calculations reproduce the Gaussian behavior found by the exact calculation. The corresponding spectral density is shown in Fig. 2 for different chain sizes. The Dirac δ functions are approximated by a rectangle of unit area and width a = 0.1. That is the best value for the width a to reduce the fluctuations due to finite-size effects. Those fluctuations decrease in amplitude as one considers larger system sizes. The frequency-dependent Gaussian of the exact result is already masked by the curve for L = 13. Therefore, the method works just fine with the transverse Ising model, and very likely will do so with the transverse ANNNI model. In the following, consider the representative cases B = 0.5, 1.0, and 2.0. These cases should cover the relevant possibilities for B in the transverse ANNNI model. Consider first B = 0.5. The time correlation function C(t) is shown in Fig. 3 for different next-to-nearest neighbor couplings J 2 . There are pairs of curves for a given J 2 , dashed lines for L = 12 and solid lines for L = 13. Those two lines agree very well with each other for the range of time t displayed. The quantitative agreement between the L = 12 and L = 13 curves is an indication that within the accuracy used, the thermodynamic value has already been obtained. The features shown in Fig. 3 are real and will not change in the thermodynamic limit. They possibly could be traced back to the rich groundstate phase diagram, however, a careful investigation is still necessary to clarify that point. In general, the decay of C(t) with time is slower for larger J 2 . The corresponding spectral density is displayed in Fig. 4, calculated for L = 13. Other than the height near the origin ω = 0, the remaining plots should not change essentially for larger chain sizes, or at the thermodynamic limit. The distinctive feature is the enhancement of the central mode as J 2 increases.
The time correlation function for B = 1 is depicted in Fig. 5 for several J 2 . The curves shown are from L = 12 and L = 13. Note that for J 2 = 0 the calculation reproduces the known Gaussian solution of the TI model. When J 2 = 0.5, oscillations are present in C(t). For J 2 ≥ 1 the curves decay at much slower rate. A careful examination of the figures shows oscillations of relatively small amplitudes. The spectral density S(ω) is shown in Fig. 6, where the calculations were done with L = 13. For J 2 = 0 the Gaussian of the TI model is reproduced. We observe an enhancement of the central model behavior as the values of J 2 are increased.
Finally, consider the case where the transverse field is larger (B = 2) than the Ising coupling. The time correlation function is shown in Fig. 7 for several values of J 2 . The curves were obtained from a chain of size L = 13. For small values of J 2 the correlation function, C(t), shows oscillations typical of collective mode, such as that found in the TI model (J 2 = 0). As J 2 becomes larger, the amplitude of the oscillations decreases. For large enough J 2 , the system displays an enhancement of the central model. Figure 8 depicts the corresponding spectral density S(ω) for the values of J 2 used in the previous figure. For J 2 = 0, the dynamics is dominated by the twopeak structure characteristic of collective mode. As J 2 increases, a reduction of the intensity of the peaks of the collective mode is observed in tandem with a growth of the central peak. For J 2 ≥ 2 the dynamics seems to be dominated entirely by the central mode. The recurrants ∆ ν of the method of recurrence relations are calculated numerically for J 2 = 1 and several values of B. First the moments µ are obtained by using Eq. (54). Next, use the conversion formulas Eq. (31). Table 1 shows some numerical results for the recurrants when B = 1.0, and J 2 = 1.0, obtained for L = 11, 12, and 13. The rightmost column shows the extrapolated value of ∆ ν for L = ∞. As can be seen, with relatively small chain sizes (L ≤ 13), one can infer the thermodynamic value of the lower-order recurrants. Higher order recurrants are still obtained, but with lesser accuracy.
The results for the thermodynamic estimates of the recurrants are shown in Fig. 9 for B = 1.0 and various values for J 2 . For J 2 = 0 the linear behavior that leads to Gaussian behavior is recovered [45]. As J 2 increases, ∆ ν increases at higher rates on the average and becomes rather erratic, therefore, it is difficult to predict a trend based on their behavior. Still the results shown are already the thermodynamic values, and it is very difficult to devise extrapolation schemes for the ∆. Notwithstanding, such an endeavor will not uncover any new physics in regard to the dynamics of the transverse ANNNI model considered here.

Summary and perspectives
The dynamical correlation functions play a crucial role in the fluctuation-dissipation theorem and in the linear response theory. However, the calculation of those quantities is often a very complicated problem in itself. The method of recurrence relations is an exact procedure that allows one to obtain of time correlation functions, spectral densities, and dynamical structure factors. We have shown the main features of the method and the inherent difficulties one might encounter in an attempt to apply to a many-body problem. Another method         Table 1: Recurrants for the transverse ANNNI model, with B = 1.0, J 2 = 1.0 and several chain sizes. The rightmost column is the extrapolation for the thermodynamic limit (L = ∞). that is showing great potential is exact diagonalization, a numerical method which relies mostly on computer capabilities. Nevertheless, the two methods can be used together, one complementing the other, to achieve progress in the calculation of dynamical correlation functions.