{'fraction_non_alphanumeric': 0.039757875911404596, 'fraction_numerical': 0.01508689870225157, 'mean_word_length': 3.965391621129326, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density Matrices for a Chain of Oscillators

We consider chains with an optical phonon spectrum and study the reduced density matrices which occur in density-matrix renormalization group (DMRG) calculations. Both for one site and for half of the chain, these are found to be exponentials of bosonic operators. Their spectra, which are correspondingly exponential, are determined and discussed. The results for large systems are obtained from the relation to a two-dimensional Gaussian model.


Introduction
The success of the density-matrix renormalization method (DMRG) in treating onedimensional quantum systems [1,2] 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 [3] 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 [4]- [13]. 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 [8,4]. 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 ω 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 and will frequently use the form ω 0 = 1 − k, so that for k = 0 there is no dispersion, while for k → 1 the spectrum becomes acoustic and the system critical. We first consider in Section 2 the density matrix ρ 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 [14]. The eigenfunctions have the character of squeezed states and are used later for numerical calculations. In Section 3, we turn to the density matrix ρ h for half of the system. We treat the case of small and large L explicitely and find that ρ 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 [3] 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.

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 [8]. Here, it can be found analytically.
The ground state of H in (1) has the form where x = (x 1 , x 2 , · · · , x L ). The matrix is determined by the frequencies ω q and the eigenvectors φ q (i) of the normal modes. From the total density matrix one then obtains the reduced one for oscillator l by integrating over all other coordinates with the constants 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 where y 2 = b x 2 l /2 and ω 2 /4 = (a − b)/b. Writing this in terms of Bose operators α , α † one can bring it into diagonal form by an equation-of-motion method. The necessary Bogoljubov transformation is with e θ = (1 + As a result, one finds that ρ 1 has the form is the Hamiltonian of a harmonic oscillator with energy Therefore the eigenvalues of ρ 1 are w n = K e −εn , n ≥ 0, and the spectrum is purely exponential. The constant K follows from the sum rule Tr(ρ 1 ) = n w 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 ρ 1 (x l , x l ) in (5), this quantity depends on the difference (a − b) and thus has no direct relation to ε.
In Fig.1, ε is shown as a function of k, putting ω 0 = (1 − k). For k → 0 it diverges logarithmically. In this limit the influence of the second oscillator vanishes, Ψ(x) becomes a product state and one is left with only one nonzero eigenvalue w 0 = 1. For k → 1, on the other hand, ε goes to zero as √ 1 − k 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 → ∞, 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 ε varies with the position along the chain. The result is shown in Fig.2 for several values of k. One sees that ε is large at the ends. This corresponds to the fact that the influence of the environment is smaller there. At the next site, however, ε 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 ε-values are not very large, though, as seen in the figure.
Due to the form of ρ 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 = ω 0 /γ, . 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.

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 ρ h has already been obtained in Section 2. We therefore proceed immediately to the case L = 4. First, we note that the square root ρ 1/2 h follows directly from Ψ. If the coordinates along the chain are ( Taking into account the form (2) and the symmetries, this leads to where the symmetric (2 × 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 ρ 1/2 h (and thus ρ h ) into diagonal form. Some of the details are given in Appendix A. The final result is that ρ h has also the form (11) where now describes two harmonic oscillators with energies ε 1 and ε 2 . Thus one obtains a simple generalization of the case L = 2. Also the variation of ε 1 with k is very similar to that of ε in Section 2. This is shown in Fig.3, where both quantities are plotted. In particular, one finds that they coincide in the limit k → 0. The ratio ε 2 /ε 1 equals 3 for small k, drops to a minimum of 2.866 for k = 0.34 and then increases continuously, because ε 2 , in contrast to ε 1 , stays finite as k → 1. The shape of the spectrum, which depends on the ratio ε 2 /ε 1 , will be discussed in the next section. At this point one can already conjecture that the structure of ρ 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 [17,3] one first relates ρ 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 ρ h the form (11) with an operator H 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 [18] with where I(k) is the complete elliptic integral of the first kind and k ′ = √ 1 − k 2 . Therefore H describes an infinite set of harmonic oscillators with energies ε j = (2j − 1)ε and is a straightforward extension of the results for small L.
The parameter ε ≡ ε 1 is also shown in Fig.3. For k → 0, it has exactly the same expansion as ε for L = 2 and ε 1 for L = 4. For k → 1, 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 [3]. The only difference is that there the operators β, β † are fermionic (so that β † β = 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 [19]. The consequences for the spectrum of ρ h are discussed below.

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 ρ 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 ε 1 , ε 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 ε 2 ∼ = 3ε 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 ε 2 deviates more from 3ε 1 in this case. For the two largest k's, ε 2 ∼ = 4ε 1 and ε 2 ∼ = 6ε 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 ρ 1 , a small number (5 − 7) is sufficient for not too large k. For example, if k = 0.5 and one chooses the same r states (with some average ε-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ε 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 ε corresponding to L = 30. With m = 7 kept states per block, the error in E 0 /L was about 3 · 10 −7 for k = 0.5. Fig.6 shows the resulting spectra for L = 6 and L = 14, together with the thermodynamic limit according to  , (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 ε 3 ∼ = 5ε 1 , which first appears for L = 6. This shows that, indeed, the number of oscillators in ρ 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 ε 2 = 3ε 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 ε j are also numerically very close to the large-L limit. For example, ε 1 agrees with the exact result 4.0189 up to three decimal places. This can be understood from the short correlation length ξ ∼ 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 ε j = 2εj [3]. 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 ′ 2j−1 . Therefore the degeneracies for the bosonic case are not as large as one might expect at first.  Figure 6: Density-matrix spectrum for k = 0.5 and two sizes L, calculated with DMRG using 7 states and m = 7. Also shown is the analytical result for L → ∞.

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 [15].
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 [5] 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.

Appendices A Four Oscillators
In order to diagonalize ρ 1/2 h in (17) one proceeds as follows. First, new coordinates are introduced by a rotation (x 1 , x 2 ) → (y 1 , y 2 ) with angle ϕ and analogously for the primed quantities. This leads to new quadratic forms in the exponent, with coefficients a ij and b ij . Choosing tan 2ϕ = 2b 12 /(b 11 − b 22 ), the cross-term b 12 becomes zero. One then considers the factors which contain only (y i , y ′ i ). These can be transformed as in Section 2 and one obtains exponentials of harmonic oscillators with energies where In terms of the new coordinates z i one then has Here z i = y i /λ i , µ = a 12 λ 1 λ 2 and the λ i are given by In the final step, one expresses (24) in terms of bosonic operators α i , α † i and considers Heisenberg-like operators ρ 1/2 h α i ρ −1/2 h which are found to be linear combinations of the α i , α † i . Therefore, a transformation as in the analogous fermionic case [16] brings ρ 1/2 h into the form (11), (18) with energies ε j /2. These energies follow from a simple quadratic equation, namely where c i = cosh ν i , s i = sinh ν i and ρ = µ/2 √ ν 1 ν 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 ρ-term in (27) is unimportant, which leads to ε j ∼ = 2ν j and ε 2 ∼ = 3ε 1 . The ratio ε 2 /ε 1 thus has the same value as in the thermodynamic limit. Moreover, ε 1 has the same asymptotic form, ε 1 ∼ = 2 ln(4/k), as for L = 2 and L → ∞. This can be attributed to the short correlation length which suppresses size effects in this limit.

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 [3]. Consider a lattice with variables x (−∞ < x < ∞) at each site, a nearest-neighbour coupling energy 1 2 K (x−x ′ ) 2 and an on-site energy 1 2 ∆ x 2 , all in units of k B T . 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, provided that k = K 2 and ω 2 0 = ∆(∆ + 4K). In this case, T and H have common eigenfunctions and Ψ in (2) gives the maximal eigenvalue for T . This allows one to obtain Ψ and also ρ h from the partition function of a two-dimensional system [3,17]. 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 ≫ ξ, where ξ is the correlation length given by ξ = 2/ 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 [20]. To do this, one uses an elliptic parametrization, so that the two couplings are, for example, with the Jacobi function sn of modulus k. This parameter also determines the on-site energy ∆ and thus the distance to the critical point ∆ = 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 ′ = √ 1 − k 2 . The isotropic case corresponds to u = I(k ′ )/2. (Our notation differs slightly from that in [20]. We have interchanged k ↔ k ′ , written u instead of λθ, used x = √ λφ for the Gaussian variables and have set α = −1.) The derivative (∂ ln T /∂u) then leads again to (1) with ω 0 = (1 − k), which is the reason for choosing this parametrization in H.
As discussed in [3], the density matrix ρ h for half of the system is, for L ≫ ξ and up to a prefactor, where A, B, C, D are the corner transfer matrices for the four infinite quadrants of the twodimensional system. Due to the integrability of the Gaussian model, i.e. the Yang-Baxter equations, these have the exponential form and similarly for B, C, D, with H CTM given by This operator was studied in [18] in connection with the Hamiltonian limit u → 0 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 (2j − 1)π/2I(k ′ ). For ρ 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 ε j = (2j − 1)πI(k ′ )/I(k) and one arrives at the result (19), (20).