# Density Matrices for a Chain of Oscillators

## 1 Introduction

The success of the density-matrix renormalization method (DMRG) in treating one-dimensional quantum systems  is closely related to the properties of the involved density matrices. In the procedure, one determines the eigenvectors of these matrices and uses those with the largest eigenvalues as a truncated basis. To be able to single out a relatively small number, however, the density-matrix spectrum has to decrease rapidly enough. Indeed, it is usually found in the numerical calculations that the eigenvalues decay roughly exponentially.

In a previous publication  it was pointed out that, for non-critical integrable models, the exponential behaviour is ultimately a consequence of the Yang-Baxter equations. For two spin one-half models, the transverse Ising chain and the uniaxial Heisenberg chain, analytical formulae were given and verified in detail in DMRG calculations.

In the present article, we want to extend these considerations to phonons, i.e. to a bosonic problem. So far, comparatively few DMRG studies have been dealing with bosons -. The difference to spin systems is that the full Hilbert space always has infinite dimension. Therefore any numerical treatment has to start with a truncation. One can do this in analogy with the DMRG procedure by selecting local states via the density matrix for a single site . This is still a nontrivial quantity with an infinite number of eigenstates in a full treatment, and it is of interest to find its properties in a solvable case. The same holds, of course, for the more complicated density matrix of a half-chain which is used in the DMRG algorithm.

The system we study is a purely bosonic model, a chain of $`L`$ harmonic oscillators with frequency $`\omega _0`$, coupled together by springs. It has a gap in the phonon spectrum and is a non-critical integrable system as the spin models mentioned above. We write the Hamiltonian

$$H=\underset{i=1}{\overset{L}{}}(\frac{1}{2}\frac{^2}{x_i^2}+\frac{1}{2}\omega _0^2x_i^2)+\underset{i=1}{\overset{L1}{}}\frac{1}{2}k(x_{i+1}x_i)^2$$
(1)   
and will frequently use the form $`\omega _0=\mathrm{\hspace{0.17em}1}k`$, so that for $`k=0`$ there is no dispersion, while for $`k1`$ the spectrum becomes acoustic and the system critical.

We first consider in Section 2 the density matrix $`\rho _1`$ for one oscillator and show that it can always be written as the exponential of the Hamiltonian of a (new) harmonic oscillator. The spectrum therefore is purely exponential, with a decay rate depending on $`k`$ and (weakly) on the chosen site. This generalizes the known result for the case $`L=2`$ . The eigenfunctions have the character of squeezed states and are used later for numerical calculations. In Section 3, we turn to the density matrix $`\rho _h`$ for half of the system. We treat the case of small and large $`L`$ explicitely and find that $`\rho _h`$ has the same exponential form, with the number of oscillators in the exponent determined by the size of the system. The result in the thermodynamic limit is derived by relating the chain to a massive two-dimensional Gaussian model and its corner transfer matrices (CTMs). It is very similar to that for the spin chains in  which lead to fermionic operators instead of bosonic ones. In particular, the spectrum without the degeneracies is purely exponential. Its form for different values of $`k`$ and different sizes $`L`$ is discussed in more detail in Section 4, including numerical results obtained by truncation and by DMRG calculations. These also illustrate, to what extent the degeneracies are reproduced in an approximate treatment. The concluding Section 5 contains a summary and additional remarks. Some details concerning the case $`L=4`$ and the Gaussian model are given in two appendices.

## 2 Density Matrix for One Oscillator

In this section we consider the case where one oscillator is singled out and all other form the environment. The corresponding density matrix (determined numerically) was used previously in the study of an electron-phonon system . Here, it can be found analytically.

The ground state of $`H`$ in (1) has the form

$$\mathrm{\Psi }(𝒙)=C\mathrm{exp}(\frac{1}{2}\underset{ij}{}A_{ij}x_ix_j)$$
(2)   
where $`𝒙=(x_1,x_2,\mathrm{},x_L)`$. The matrix

$$A_{ij}=\underset{q}{}\omega _q\varphi _q(i)\varphi _q(j)$$
(3)   
is determined by the frequencies $`\omega _q`$ and the eigenvectors $`\varphi _q(i)`$ of the normal modes. From the total density matrix

$$\rho (𝒙\text{,}𝒙^{\mathbf{}})=\mathrm{\Psi }(𝒙)\mathrm{\Psi }(𝒙^{\mathbf{}})$$
(4)   
one then obtains the reduced one for oscillator $`l`$ by integrating over all other coordinates $`x_i=x_i^{}`$. This leads to

$$\rho _1(x_l,x_l^{})=C_1\mathrm{exp}(\frac{1}{2}(ab)x_l^2)\mathrm{exp}(\frac{b}{4}(x_lx_l^{})^2)\mathrm{exp}(\frac{1}{2}(ab)x_{l}^{}{}_{}{}^{2})$$
(5)   
with the constants

$`a`$ $`=`$ $`A_{ll}`$ (6)   
$`b`$ $`=`$ $`{\displaystyle \underset{i,jl}{}}A_{li}[A^{(l)}]_{ij}^1A_{jl}`$ (7)   
where $`A^{(l)}`$ is the matrix obtained from $`A`$ by deleting the $`l`$-th row and column. The second exponential in (5) can be transformed into a differential operator, giving

$$\rho _1=C_2\mathrm{exp}(\frac{1}{4}\omega ^2y^2)\mathrm{exp}(\frac{1}{2}\frac{^2}{y^2})\mathrm{exp}(\frac{1}{4}\omega ^2y^2)$$
(8)   
where $`y^2=bx_l^2/2`$ and $`\omega ^2/4=(ab)/b`$. Writing this in terms of Bose operators $`\alpha `$ , $`\alpha ^{}`$ one can bring it into diagonal form by an equation-of-motion method. The necessary Bogoljubov transformation is

$$\beta =\mathrm{cosh}\theta \alpha +\mathrm{sinh}\theta \alpha ^{}$$
(9)   
with

$$e^\theta =(1+\frac{\omega ^2}{4})^{1/4}$$
(10)   
As a result, one finds that $`\rho _1`$ has the form

$$\rho _1=K\mathrm{exp}()$$
(11)   
where

$$=\epsilon \beta ^{}\beta $$
(12)   
is the Hamiltonian of a harmonic oscillator with energy

$$\epsilon =\mathrm{\hspace{0.33em}2}\mathrm{sinh}^1(\frac{\omega }{2})=\mathrm{\hspace{0.33em}2}\mathrm{sinh}^1\left(\sqrt{a/b1}\right)$$
(13)   
Therefore the eigenvalues of $`\rho _1`$ are $`w_n=Ke^{\epsilon n},n0`$, and the spectrum is purely exponential. The constant $`K`$ follows from the sum rule $`\mathrm{Tr}(\rho _1)=_nw_n=1`$.

This result is completely general. The details of the oscillating system and the position of the chosen oscillator enter only via the ratio $`a/b`$. The same constants $`a`$ and $`b`$ also determine the probability to find a certain elongation $`x_l`$. However, as seen from $`\rho _1(x_l,x_l)`$ in (5), this quantity depends on the difference ($`ab`$) and thus has no direct relation to $`\epsilon `$.

In the simplest case of two oscillators ($`L=\mathrm{\hspace{0.17em}2}`$) one finds explicitely

$$\epsilon =\mathrm{\hspace{0.25em}2}\mathrm{sinh}^1\left(\sqrt{4\omega _1\omega _2/(\omega _1\omega _2)^2}\right)$$
(14)   
or, equivalently,

$$\epsilon =\mathrm{ln}(\mathrm{coth}^2(\frac{\eta }{2}))$$
(15)   
where $`\omega _1=\omega _0`$ , $`\omega _2=\sqrt{\omega _0^2+2k}`$ are the two eigenfrequencies and $`e^{2\eta }=\omega _2/\omega _1`$. This is the result obtained in  in a different way.

In Fig.1, $`\epsilon `$ is shown as a function of $`k`$, putting $`\omega _0=(1k)`$. For $`k0`$ it diverges logarithmically. In this limit the influence of the second oscillator vanishes, $`\mathrm{\Psi }(𝒙)`$ becomes a product state and one is left with only one nonzero eigenvalue $`w_0=\mathrm{\hspace{0.17em}1}`$. For $`k1`$, on the other hand, $`\epsilon `$ goes to zero as $`\sqrt{1k}`$ and the eigenvalues $`w_n`$ decrease only very slowly, which reflects the strong coupling. These features are encountered also in all other cases. For $`L=3`$ one can still give explicit analytical formulae, but for larger $`L`$ the problem has to be treated numerically. In the figure, two additional cases, $`L=10`$ and $`L=100`$, are shown. The limit $`L\mathrm{}`$, which is approached exponentially in $`L`$ with a correlation length increasing with $`k`$, is indistinguishable from $`L=100`$ on the given scale.

One can also investigate how $`\epsilon `$ varies with the position along the chain. The result is shown in Fig.2 for several values of k. One sees that $`\epsilon `$ is large at the ends. This corresponds to the fact that the influence of the environment is smaller there. At the next site, however, $`\epsilon `$ drops and then approaches the bulk value from below as one moves into the interior. The approach becomes slower as $`k`$ increases. The overall differences in the $`\epsilon `$-values are not very large, though, as seen in the figure.

Due to the form of $`\rho _1`$, its eigenstates are standard oscillator functions of a coordinate $`z`$ which differs from $`x_l`$ by a scale factor. Compared with the eigenfunctions of the uncoupled oscillator $`l`$, they are squeezed states whose spatial extent is reduced by a factor $`q=\sqrt{\omega _0/\gamma }`$, where $`\gamma =\sqrt{a(ab)}`$. For small $`k`$, $`q`$ approaches one and the two sets of functions coincide. With increasing $`k`$, the amount of squeezing increases, and it is then advantageous to choose the squeezed states as a local basis. This was done in the numerical calculations which will be presented in Section 4.

## 3 Density Matrix for a Half-Chain

We now turn to the central quantity in the usual DMRG calculations, the reduced density matrix for half of the system. It enters each time the system is enlarged in the infinite-size algorithm. We will determine its spectrum in the two limits of small and large $`L`$.

For $`L=2`$, one-half of the system is just one oscillator and $`\rho _h`$ has already been obtained in Section 2. We therefore proceed immediately to the case $`L=4`$. First, we note that the square root $`\rho _h^{1/2}`$ follows directly from $`\mathrm{\Psi }`$. If the coordinates along the chain are $`(x_2,x_1,x_1^{},x_2^{})`$, one has

$$\rho _h^{1/2}(x_1,x_2;x_1^{},x_2^{})=\mathrm{\Psi }(x_2,x_1,x_1^{},x_2^{})$$
(16)   
Taking into account the form (2) and the symmetries, this leads to

$$\rho _h^{1/2}=C\mathrm{exp}\left\{\frac{1}{2}\underset{ij}{}a_{ij}(x_ix_j+x_i^{}x_j^{})\underset{ij}{}b_{ij}x_ix_j^{}\right\}$$
(17)   
where the symmetric $`(2\times 2)`$ matrices $`a_{ij}`$ and $`b_{ij}`$ follow from the matrix $`A`$ of Section 2. Altogether one has six different coefficients which couple the variables as shown in the following figure.

The cross-couplings, shown as dashed lines, can be eliminated by introducing new coodinates $`y_i,y_i^{}`$. After that, a sequence of transformations similar to those in Section 2 brings $`\rho _h^{1/2}`$ (and thus $`\rho _h`$) into diagonal form. Some of the details are given in Appendix A. The final result is that $`\rho _h`$ has also the form (11) where now

$$=\underset{j=1}{\overset{2}{}}\epsilon _j\beta _j^{}\beta _j$$
(18)   
describes two harmonic oscillators with energies $`\epsilon _1`$ and $`\epsilon _2`$. Thus one obtains a simple generalization of the case $`L=2`$. Also the variation of $`\epsilon _1`$ with $`k`$ is very similar to that of $`\epsilon `$ in Section 2. This is shown in Fig.3, where both quantities are plotted. In particular, one finds that they coincide in the limit $`k0`$. The ratio $`\epsilon _2/\epsilon _1`$ equals $`3`$ for small k, drops to a minimum of $`2.866`$ for $`k=0.34`$ and then increases continuously, because $`\epsilon _2`$, in contrast to $`\epsilon _1`$, stays finite as $`k1`$. The shape of the spectrum, which depends on the ratio $`\epsilon _2/\epsilon _1`$, will be discussed in the next section.

At this point one can already conjecture that the structure of $`\rho _h`$ remains the same also for larger $`L`$. A direct computation as above does not seen feasible, though. In the limit of large $`L`$, however, a different approach is possible. As in  one first relates $`\rho _h`$ to the partition function of a two-dimensional classical system, which is a massive Gaussian model in our case, in the form of an infinite strip of width $`L`$ with a perpendicular cut. This connection is discussed in more detail in Appendix B. One then expresses the partition function as the product of four corner transfer matrices. In the case where $`L`$ is much larger than the correlation length, one can use the thermodynamic limit of these CTMs and finds for $`\rho _h`$ the form (11) with an operator $``$ which is very similar to $`H`$ in (1). The coefficients, however, are multiplied by additional site-dependent factors which increase linearly along the chain and reflect the corner geometry. Up to a prefactor it is the operator given in (32) in Appendix B and its diagonalization amounts to finding the normal modes of the corresponding vibrational problem. From the results in  one obtains

$$=\underset{j1}{}(2j1)\epsilon \beta _j^{}\beta _j$$
(19)   
with

$$\epsilon =\pi \frac{I(k^{})}{I(k)}$$
(20)   
where $`I(k)`$ is the complete elliptic integral of the first kind and $`k^{}=\sqrt{1k^2}`$. Therefore $``$ describes an infinite set of harmonic oscillators with energies $`\epsilon _j=(2j1)\epsilon `$ and is a straightforward extension of the results for small $`L`$.

The parameter $`\epsilon \epsilon _1`$ is also shown in Fig.3. For $`k0`$, it has exactly the same expansion as $`\epsilon `$ for $`L=2`$ and $`\epsilon _1`$ for $`L=4`$. For $`k1`$, it vanishes only logarithmically, i.e. more slowly than the quantities for finite $`L`$.

One should note that the results (19),(20) are formally the same as for the transverse Ising chain in the disordered phase . The only difference is that there the operators $`\beta ,\beta ^{}`$ are fermionic (so that $`\beta ^{}\beta =0,1`$), whereas here they are bosonic. Such similarities can also be observed in the row transfer matrices of the Gaussian and the Ising model, if one uses the corresponding parametrizations . The consequences for the spectrum of $`\rho _h`$ are discussed below.

## 4 Spectra and Numerics

In the following, we show the density-matrix spectra for half-chains and discuss some numerical aspects. In the figures, the eigenvalues $`w_n`$ of $`\rho _h`$ are ordered according to magnitude and plotted on a semilogarithmic scale.

Figure 4 shows spectra for $`L=4`$ and several values of $`k`$. These results were obtained by calculating the two energies $`\epsilon _1,\epsilon _2`$ numerically from the formulae in Appendix A. Apart from the rapid decrease, one notes a clear ladder structure for the smallest three $`k`$’s. It results from the relation $`\epsilon _23\epsilon _1`$ which leads to the approximate degeneracies $`(1,1,1,2,2,2,3)`$ for the first seven levels. The steps for $`k=0.3`$ are less perfect, since $`\epsilon _2`$ deviates more from $`3\epsilon _1`$ in this case. For the two largest $`k`$’s, $`\epsilon _24\epsilon _1`$ and $`\epsilon _26\epsilon _1`$, so that the first step appears at these levels and the spectra look more stretched out.

It is interesting to see, how these results are recovered in a numerical treatment using a truncated Hilbert space. If one works with the eigenstates of $`\rho _1`$, a small number $`(57)`$ is sufficient for not too large $`k`$. For example, if $`k=0.5`$ and one chooses the same $`r`$ states (with some average $`\epsilon `$-value) for all four sites, the error in the ground-state energy $`E_0/L`$ is of the order $`10^r`$. The spectra which then result are shown in Fig.5 for three values of $`r`$. The first $`w_n`$ are always quite accurate, but there are characteristic differences for the following ones, which are connected with the number of steps, i.e. with the degeneracies. One can see that if $`r`$ states are kept, the pattern is correct for the first $`r`$ levels (counted from the top). At the next level, the state with energy $`r\epsilon _1`$ is missing and the corresponding step is absent. Thus there is a certain correspondence between the states in the local basis and in the density matrix. For smaller $`w_n`$, however, the situation is less clear, and the spectrum finally becomes irregular. The tails of the approximate spectra always lie below the exact one.

In order to obtain results also for $`L>4`$, we have carried out DMRG calculations, using $`7`$ states at each site, with an $`\epsilon `$ corresponding to $`L=30`$. With $`m=7`$ kept states per block, the error in $`E_0/L`$ was about $`310^7`$ for $`k=0.5`$. Fig.6 shows the resulting spectra for $`L=6`$ and $`L=14`$, together with the thermodynamic limit according to (19),(20). One notes, that the spectra for the two $`L`$’s are similar, though not identical. Compared to $`L=4`$, the degeneracies have changed to $`(1,1,1,2,2,3,4)`$. The latter two result from a third energy $`\epsilon _35\epsilon _1`$, which first appears for $`L=6`$. This shows that, indeed, the number of oscillators in $`\rho _h`$ is equal to the size $`L/2`$ of the half-chain. One also sees that for $`L=14`$ the first two steps have become perfect, so that $`\epsilon _2=3\epsilon _1`$ as for the infinite system. Up to some small deviations, this also holds for the next two steps. Only for the remaining levels $`8`$ and $`9`$, the degeneracies are not correct. This is the same effect as found above for $`L=4`$.

For $`L=14`$, the $`\epsilon _j`$ are also numerically very close to the large-$`L`$ limit. For example, $`\epsilon _1`$ agrees with the exact result $`4.0189`$ up to three decimal places. This can be understood from the short correlation length $`\xi 3`$ for $`k=0.5`$ which makes size-effects small. Finally, we want to mention that, in the thermodynamic limit, the multiplicities are just one-half of those found in the fermionic case for the ordered phase where $`\epsilon _j=2\epsilon j`$ . This is because the number $`P_j`$ of partitions without repetition is the same as that of the odd integers with repetition, $`P_j=P_{2j1}^{}`$. Therefore the degeneracies for the bosonic case are not as large as one might expect at first.

## 5 Conclusion

We have investigated a bosonic system, where the ground-state density matrices can be determined explicitely in various cases. It turned out that they are exponentials of oscillator Hamiltonians, so that all results are quite transparent. The spectra have exponential character and the eigenfunctions are oscillator states. For the single-site density matrix, these states are related to those of the chain oscillators by squeezing. For the half-chain density matrix, they are connected with certain normal modes concentrated near the middle of the system. The thermodynamic limit was obtained in the same way as for the integrable spin chains treated previously, and the spectra are very similar to those found there. By counting the degeneracies, one would arrive at formulae as given in .

Taking all this together, the chain treated here may serve as a standard example where one can see the features of the density matrices in detail. In this context, it would still be interesting to determine the half-chain density matrix for arbitrary sizes, in particular at the critical point, where the vibrational spectrum becomes acoustic. This case has already been studied by DMRG  but, as for the critical spin models, the density-matrix spectra have yet to be explained. Another question is, whether the model of coupled oscillators, for which the ground state is known explicitely, could be used to study density matrices in higher dimensions. For the DMRG method, it would be quite important to know if the spectral properties change in this case.

Acknowledgements

We thank M. Kaulke, E. Jeckelmann, G. Babudjan, A. Pelster and H. Kleinert for discussions. M.C. Chung acknowledges the support of Deutscher Akademischer Austauschdienst (DAAD).

## Appendices

## Appendix A Four Oscillators

In order to diagonalize $`\rho _h^{1/2}`$ in (17) one proceeds as follows. First, new coordinates are introduced by a rotation $`(x_1,x_2)(y_1,y_2)`$ with angle $`\phi `$ and analogously for the primed quantities. This leads to new quadratic forms in the exponent, with coefficients $`\widehat{a}_{ij}`$ and $`\widehat{b}_{ij}`$. Choosing $`\mathrm{tan}2\phi =2b_{12}/(b_{11}b_{22})`$, the cross-term $`\widehat{b}_{12}`$ becomes zero. One then considers the factors

$$\mathrm{exp}(1/2\widehat{a}_{ii}y_i^2)\mathrm{exp}(\widehat{b}_{ii}y_iy_i^{})\mathrm{exp}(1/2\widehat{a}_{ii}y_{i}^{}{}_{}{}^{2}),i=1,2$$
(21)   
which contain only $`(y_i,y_{}^{}{}_{i}{}^{})`$. These can be transformed as in Section 2 and one obtains exponentials of harmonic oscillators with energies

$$\nu _i=\mathrm{\hspace{0.33em}2}\mathrm{sinh}^1(\mathrm{\Omega }_i/2)$$
(22)   
where

$$\mathrm{\Omega }_i/2=\sqrt{(\widehat{a}_{ii}+\widehat{b}_{ii})/(2\widehat{b}_{ii})}$$
(23)   
In terms of the new coordinates $`z_i`$ one then has

$$\rho _h^{1/2}=C\mathrm{exp}(\mu z_1z_2)\mathrm{exp}(\underset{i}{}(\frac{1}{2}\frac{^2}{z_i^2}+\frac{1}{2}\nu _i^2z_i^2))\mathrm{exp}(\mu z_1z_2)$$
(24)   
Here $`z_i=y_i/\lambda _i,\mu =\widehat{a}_{12}\lambda _1\lambda _2`$ and the $`\lambda _i`$ are given by

$$\lambda _i=\left(\frac{\nu _i}{\widehat{b}_{ii}\mathrm{\Omega }_i}\right)^{1/2}\left(1+\frac{\mathrm{\Omega }_i^2}{4}\right)^{1/4}$$
(25)   
In the final step, one expresses (24) in terms of bosonic operators $`\alpha _i,\alpha _i^{}`$ and considers Heisenberg-like operators $`\rho _h^{1/2}\alpha _i\rho _h^{1/2}`$ which are found to be linear combinations of the $`\alpha _i,\alpha _i^{}`$. Therefore, a transformation as in the analogous fermionic case 

$$\beta _j=\underset{i}{}(g_{ji}\alpha _i+h_{ji}\alpha _i^{})$$
(26)   
brings $`\rho _h^{1/2}`$ into the form (11),(18) with energies $`\epsilon _j/2`$. These energies follow from a simple quadratic equation, namely

$$\mathrm{cosh}\frac{\epsilon _j}{2}=\frac{1}{2}(c_1+c_2)\pm \sqrt{\frac{1}{4}(c_1c_2)^2+4\rho ^2s_1s_2}$$
(27)   
where $`c_i=\mathrm{cosh}\nu _i,s_i=\mathrm{sinh}\nu _i\text{and}\rho =\mu /2\sqrt{\nu _1\nu _2}.`$

These quantities have to be evaluated, starting from the initial constants $`a_{ij}`$ and $`b_{ij}`$, which are simple analytic expressions involving the four eigenfrequencies of the chain. It turns out that, for small values of $`k`$, the $`\rho `$-term in (27) is unimportant, which leads to $`\epsilon _j2\nu _j`$ and $`\epsilon _23\epsilon _1`$. The ratio $`\epsilon _2/\epsilon _1`$ thus has the same value as in the thermodynamic limit. Moreover, $`\epsilon _1`$ has the same asymptotic form, $`\epsilon _12\mathrm{ln}(4/k)`$, as for $`L=2`$ and $`L\mathrm{}`$. This can be attributed to the short correlation length which suppresses size effects in this limit.

## Appendix B Relation to the Gaussian Model

The Hamiltonian $`H`$ in (1) has a close relation to the transfer matrix of a two-dimensional Gaussian model (GM). The connection is the same as between the transverse Ising chain and the two-dimensional Ising model . Consider a lattice with variables $`x(\mathrm{}<x<\mathrm{})`$ at each site, a nearest-neighbour coupling energy $`\frac{1}{2}K(xx^{})^2`$ and an on-site energy $`\frac{1}{2}\mathrm{\Delta }x^2`$, all in units of $`k_BT`$. If the lattice is oriented diagonally, the appropriate transfer matrix $`T`$ involves the piece shown in the figure below.

One can then verify by a simple direct calculation (using two interpenetrating lattices) that, with periodic boundary conditions,

$$[H,T]=\mathrm{\hspace{0.25em}0}$$
(28)   
provided that $`k=K^2`$ and $`\omega _0^2=\mathrm{\Delta }(\mathrm{\Delta }+4K)`$. In this case, $`T`$ and $`H`$ have common eigenfunctions and $`\mathrm{\Psi }`$ in (2) gives the maximal eigenvalue for $`T`$. This allows one to obtain $`\mathrm{\Psi }`$ and also $`\rho _h`$ from the partition function of a two-dimensional system . If the GM has open boundaries, one has to modify $`H`$ at the end, so as to preserve (28). However, for a system with $`L\xi `$, where $`\xi `$ is the correlation length given by $`\xi =2/\mathrm{ln}(1/k)`$, this effect is not important and can be neglected.

An alternative approach is to consider a Gaussian model with anisotropic couplings for periodic boundary conditions, to show that the $`T^{}`$s for different anisotropies commute and to realize that a proper derivative leads to $`H`$ . To do this, one uses an elliptic parametrization, so that the two couplings are, for example,

$$K_1=i/\mathrm{sn}(iu,k);K_2=ik\mathrm{sn}(iu,k)$$
(29)   
with the Jacobi function $`sn`$ of modulus $`k`$. This parameter also determines the on-site energy $`\mathrm{\Delta }`$ and thus the distance to the critical point $`\mathrm{\Delta }=\mathrm{\hspace{0.17em}0}`$, as well as the correlation length. The parameter $`u`$, on the other hand, specifies the ratio $`K_1/K_2`$. It varies between $`0`$ and $`I(k^{})`$, where $`I`$ is the complete elliptic integral of the first kind and $`k^{}=\sqrt{1k^2}`$. The isotropic case corresponds to $`u=I(k^{})/2`$. (Our notation differs slightly from that in . We have interchanged $`kk^{}`$, written $`u`$ instead of $`\lambda \theta `$, used $`x=\sqrt{\lambda }\varphi `$ for the Gaussian variables and have set $`\alpha =1`$.) The derivative $`(\mathrm{ln}T/u)`$ then leads again to (1) with $`\omega _0=(1k)`$, which is the reason for choosing this parametrization in $`H`$.

As discussed in , the density matrix $`\rho _h`$ for half of the system is, for $`L\xi `$ and up to a prefactor,

$$\rho _h=ABCD$$
(30)   
where $`A,B,C,D`$ are the corner transfer matrices for the four infinite quadrants of the two-dimensional system. Due to the integrability of the Gaussian model, i.e. the Yang-Baxter equations, these have the exponential form

$$A=e^{u_{\mathrm{CTM}}}$$
(31)   
and similarly for $`B,C,D`$, with $`_{\mathrm{CTM}}`$ given by

$$_{CTM}=\underset{n1}{}\left\{\frac{1}{2}(2n1)\frac{^2}{x_n^2}+\frac{1}{2}(2n1)(1k)^2x_n^2+\frac{1}{2}2nk(x_{n+1}x_n)^2\right\}$$
(32)   
This operator was studied in  in connection with the Hamiltonian limit $`u0`$ of $`A`$, where one can determine its form simply by inspection. It is associated with a corner of Ramond type, i.e. without a site at the tip. In terms of vibrations, it describes a system of coupled oscillators, where the spring constants and inverse masses increase along the chain. It can be diagonalized with the help of Carlitz polynomials and then becomes the sum of harmonic oscillators with eigenvalues $`(2j1)\pi /2I(k^{})`$. For $`\rho _h`$ one needs $`ABCD`$, or $`A^4`$ if one has an isotropic model. In either case this gives a factor $`2I(k^{})`$, so that the energy becomes $`\epsilon _j=(2j1)\pi I(k^{})/I(k)`$ and one arrives at the result (19),(20).