{'fraction_non_alphanumeric': 0.033222984157011116, 'fraction_numerical': 0.01622132891936628, 'mean_word_length': 4.26531374501992, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bose-Einstein condensation and the magnetically ordered state of TlCuCl3

The dimerized S=1/2 spins of the Cu ions in TlCuCl3 are ordered antiferromagnetically in the presence of a field larger than about 54 kOe in the zero-temperature limit. Within the mean-field approximation all thermal effects are frozen out below 6 K. Nevertheless, experiments show significant changes of the critical field and the magnetization below this temperature, which reflect the presence of low-energetic dimer-spin excitations. We calculate the dimer-spin correlation functions within a self-consistent random-phase approximation, using as input the effective exchange coupling parameters obtained from the measured excitation spectra. The calculated critical field and magnetization curves exhibit the main features of those measured experimentally, but differ in important respects from the predictions of simplified boson models.


I. INTRODUCTION
The concept of Bose-Einstein condensation dates back more than 80 years to the prediction of Einstein, based on Bose's work on the statistics of photons, that a gas of non-interacting massive bosons would condense below a certain critical temperature T c . The condensation implies that below T c a non-zero fraction of the total number of particles occupies the lowest single-particle quantum state. For dilute atomic gases this phenomenon was realized experimentally in 1995 for trapped clouds of alkali atoms (see e.g. Ref. 1).
For a uniform gas of density n the transition temperature is given by kT c ≈ 3.31 2 n 2/3 /m, where m is the particle mass. For particles trapped in a harmonic oscillator potential (trap frequencies ω x , ω y and ω z ) one has kT c ≈ 0.94 (N ω x ω y ω z ) 1/3 , where N is the total number of particles. In the latter case the particle mass enters through the trap frequencies, equal to the square root of the force constants in the three directions divided by the particle mass. When a trapped gas is dilute in the sense that the atom-atom scattering length is much less than the interatomic distance, the observed transition temperatures agree well with theoretical expectation for a noninteracting gas. For less dilute gases interaction effects give rise to an observable small shift of T c proportional to the scattering length.
The condensation of massive bosons into a single quantum state is intimately connected to the conservation of particle number. For massless bosons such as phonons or magnons the particle number is not fixed but depends on temperature, and there is therefore no Bose-Einstein condensation in the traditional sense of the term. However, there has been a wide use of model Hamiltonians for magnetic systems that have features in common with those of interacting, massive bosons. The aim of the present work is to consider one such specific system, that of the dimerized Cu 2+ spins in TlCuCl 3 , and compare predictions of such models with calculations that are based on (approximate) solutions of the many-body problem of interacting spins. A recent review of experimental and theoretical developments concerning the magnetic ordering of TlCuCl 3 and related compounds has been given by Giamarchi et al. 2 Magnetization measurements 3,4 and inelastic neutronscattering experiments 5,6 demonstrate clearly that nearest-neighboring pairs of S = 1 2 spins of the Cu 2+ ions in TlCuCl 3 are dimerized leading to an S = 0 ground state and an S = 1 excited triplet around 5.2-5.7 meV above the singlet. Due to the exchange interactions between the dimers the excitations become strongly dispersive, and, in the zero-temperature limit, the minimum energy of the degenerate singlet-triplet mode, at (001), is only about 0.7 meV. When a field is applied, the energy of one of the three normal modes is reduced and goes to zero at a critical field of about 54 kOe at zero temperature. The Cu spins of isostructural KCuCl 3 are similarly dimerized, but the interdimer interactions are relatively weaker and the critical field is about 230 kOe at T = 0 in this system. 7 The phase transition shown by TlCuCl 3 at the field where the excitation energy vanishes, has been analyzed by Nikuni et al. 8 They assumed the dimer system to be described by an effective Haniltonian of the form where the bosonic operators a † and a denote "magnon" creation and annihilation operators and the positive constant v 0 denotes the strength of the repulsive magnonmagnon interaction, assumed to be a delta function in real space. The quantity µ plays the role of a chemical potential, assumed proportional to the difference between the applied magnetic field and the critical field. Using v 0 and m as fitting parameters, Nikuni et al. were able to give a reasonable account of the temperature dependence of the critical field and the magnetization along the applied field in the ordered antiferromagnetic state. An extended version of their theory based on a more realistic dispersion of the magnetic excitations was presented by Misguich and Oshikawa. 9 However, as we shall see in detail in Sec. IV below, these simplified boson models suffer from inconsistencies that originate in their neglect of the highest level of the triplet. Another, more general problem with boson models is that double occupancy of a local site should be prohibited. In this connection we mention the work of Sirker et al., 10 who used a bondoperator approach to map the spin system onto a model of interacting bosons by introducing an infinite on-site repulsion between local triplet excitations.
In the following treatment of the dimerized spin system in TlCuCl 3 we adopt a different point of view and start from the Hamiltonian for the spin system itself, using as input the effective exchange coupling that has been derived from measured excitation spectra. Our approach is a generalization of the zero-temperature theory by Matsumoto et al., 11 who used the random-phase approximation (RPA) to calculate magnetization curves and excitation spectra for TlCuCl 3 . They also considered the case when the phase transition is induced by the application of a hydrostatic pressure. 12,13 Here we only address the case of a field-induced transition, but the theory of Matsumoto et al. is extended to include both the effects of quantum fluctuations at zero temperature and the effects of thermal fluctuations. The self-consistent version of the RPA, which is the one applied here, is faced with the similar problem of double occupancy as the boson modelling. However, here this problem is found to have a natural solution by a consideration of the higher order modifications of the Green functions. The self-consistent RPA theory for the paramagnetic phase of the dimer-spin system is presented in Sec. II, which, in Sec. III, is followed by an analysis of the antiferromagnetic phase. A closer examination and discussion of the results obtained are referred to the last Sec. IV.

II. EXCITATIONS IN THE PARAMAGNETIC PHASE
A. The self-consistent RPA theory The TlCuCl 3 crystal is monoclinic (space group P 2 1 /c) and the lattice parameters are a = 3.9815Å, b = 14.144 A, c = 8.8904Å and β = 96.32 • at room temperature. 14 The crystal is constructed from layers with configuration Cu 2 Cl 6 stacked on top of each other so as to form two chains of Cu ions parallel to the a axis. The chains are separated by Tl ions and pass through the center and corners of the b-c plane in the unit cell. There are four Cu ions or two dimer pairs per unit cell. The dimer pair in the unit cell belonging to the chain through a corner is located at site 1: (x, y, z) and site 2: (x,ȳ,z), and the pair belonging to the other chain is placed at site 3: (x,ȳ + 1 2 , z + 1 2 ) and site 4: (x, y + 1 2 ,z + 1 2 ). Here x = 0.2338, y = 0.0486, and z = −0.0175, and the numbering of the sites from 1 to 4 defines the four different Cu-sublattices.
The Hamiltonian is assumed to be where s i is the spin-variable of the Cu ion at the ith site. The most important exchange parameter is ∆ = −J(i 1 i 2 ), where (i 1 i 2 ) are the nearest-neighbor Cu pairs (the 1-2 or the 3-4 ions in the unit cell). The Fourier transform of the Heisenberg exchange interactions between spins on sublattice α and β is defined in terms of the remaining coupling parameters where R i is the position of the ith dimer, and the prime indicates that the dominating interaction ∆ is excluded from the sum, (ij) = (i 1 i 2 ). When the interactions between the dimers are neglected, the Hamiltonian may be diagonalized exactly in terms of independent products of single-dimer eigenstates. The total spin of the ith dimer is S i = s i1 + s i2 , where s i1 and s i2 denote the Cu spins belonging to, respectively, the sublattices 1 and 2, or 3 and 4. The total spin defines the basis |SS z , and when the field is along the z axis, the eigenstates of the non-interacting dimer are state |3 = |1 − 1 at the energy ∆ + h, state |2 = |10 at the energy ∆, state |1 = |1 + 1 at the energy ∆ − h, state |0 = |00 at zero energy, with h = gµ B |H| and ∆ positive. When h < ∆ the ground state is the non-magnetic singlet |00 . In the present section we focus on this condition, and we shall assume that the system stays paramagnetic also in the presence of the interdimer interaction J αβ (q). The original Hamiltonian (2) may be rewritten in terms of two dimer-spin variables, the sum S i = s i1 + s i2 and the dif-ferenceS i = s i1 − s i2 . The sum operator only has nonzero matrix elements between the three excited S = 1 states. We consider the case where the populations of these levels are small (at sufficiently low temperatures in the disordered phase), in which case the dynamical effects due to S i are negligible. When S i is neglected, the Hamiltonian involves only a single effective q-dependent The presence of two equivalent dimers per unit cell yields two values for the effective interaction for each value of the wave vector within the first Brillouin zone. Alternatively, one may use an extended zone scheme with an effective basis of one dimer per unit cell, in which case only the upper sign applies. In order to study the spin dynamics of this Hamiltonian we introduce the standard basis operators 15,16,17 for the jth dimer a j µν = (|µ ν|) j , µ, ν = 0, 1, 2, 3 .
In the present case of a dimer system with stationary bonds, these operators serve the same purpose but are of more general use than the "bond operators" applied by Matsumoto et al. 11 In terms of the standard basis operators the components ofS j becomē and the Hamiltonian may be written when only the Fourier transform J(ij) of the effective interaction, Eq. (4), is included. Next we define a 6 × 6 matrix of Green functions 18 A single bracket · · · denotes the thermal expectation value, and a i (t) is a vector operator of site i at time t with components With the short-hand double-bracket notationḠ(ij, ω) = a i ; a † j , the equations of motion for the Green functions are (see for instance Ref. 17) where the new higher-order Green functions introduced by the second term are determined by the Hamiltonian (7) by the use of the commutator relation By using an RPA decoupling of the higher-order Green functions, a i µν a j µ ′ ν ′ ≈ a i µν a j µ ′ ν ′ + a i µν a j µ ′ ν ′ , one finds that the equations of motion are reduced to a closed set of equations, which are solvable by a Fourier transformation. Further, the 6 × 6 matrix equations decouple into 3 sets of 2 × 2 matrix equations, and one of these is The Green functions depend on the Fourier variables, G µν = G µν (q, ω), and we have introduced the following parameters with n µ being the average population of the µth dimer level n µ = a i µµ , and the difference n µν = n µ − n ν . (14) Inversion of the first matrix in Eq. (12) results in The poles of the Green functions determine the excitation energies E ± q , which are given by The (2,5) part ofḠ(q, ω) is given by the same expression as the (1,6) part in Eq. (15), except that ω is replaced by −ω. The result for the (3,4) part is that obtained from Eq. (15) for h = 0, corresponding to the replacement of n 01 and n 03 by n 02 , and it leads to poles at the energies ±E z q , where We introduce the following matrix of equal-time correlation functions where N is the number of dimers. According to the fluctuation-dissipation theorem (see e.g. Ref. 17) whereḠ ′′ denotes the imaginary part of the matrix Green function and β = 1/kT . By definition, the average value of, for instance, the 11-component is Hence, by calculating the q averages of the correlation functions, the Green functions may be used for determining the populations of the dimer levels. From the diagonal part of the A-matrix, we get straightforwardly and where The three equations determine the four population numbers, when they are supplemented by the exact condition that The RPA decoupling is valid in the mean-field (MF) approximation, where the thermal averages are determined by the Hamiltonian H MF for the noninteracting system, i.e. within the approximation a µν ≈ Tr[a µν exp(−βH MF )] = δ µν n MF µ . In the case of strong dispersion, this approximation certainly underestimates the populations of the excited levels. Anticipating that the correlation effects predicted by the RPA theory are reasonably trustworthy, the self-consistent equations above should lead to a more accurate determination of the population numbers than that offered by the MF approximation. Unfortunately, the present theory also predicts that the off-diagonal components of the A matrix are non-zero contradicting that, for instance, A 16 = a i 01 a i 03 should vanish identically. Phrased differently, this inconsistency implies that the RPA result for the occupation numbers is not unique but depends on the particular choice of correlation functions used in the calculation. A non-zero value of a i 01 a i 03 is the equivalent of a double occupancy of bosons at a single site. Instead of introducing an arbitrary repulsive potential, we are here going to consider possible improvements of the RPA-decoupling procedure applied above.
The present system is in many ways similar to the singlet-doublet system encountered in praseodymium metal. Improvements of the RPA for this system have been derived in Ref. 19, and this theory has been applied to the calculation of the linewidths and the energy renormalization of the excitations in two other dimer systems Cs 2 Cr 2 Br 9 and KCuCl 3 for the case of zero field. 20,21 It is straightforward to generalize the singlet-doublet theory so as to account for the presence of a third level. Here, we are going to extend the theory to the case where the field is non-zero, and, in the next section, to consider the modifications produced by an ordered moment. The presence of the interdimer interactions J(ij) implies that a i 00 does not commute with the Hamiltonian. Hence, the assumed ground state, the product state of |0 i , is not an eigenstate of the interacting system. The situation compares with the simple Heisenberg antiferromagnet, where the mean-field Néel state is not the true ground state. In this case the RPA theory predicts a zero-temperature reduction of the antiferromagnetic moment from its saturated Néel state value. Equivalently, the RPA results above imply that n 0 is smaller than its saturation value 1 at zero temperature. The single-dimer population numbers are subject to quantum as well as thermal fluctuations.
One of the terms neglected in the RPA equation (12) involves the Green function (a i 00 − a i 11 − n 01 )a j 01 ; a † k . Since a i 00 and a i 11 are not true constants of motion, this Green function may modify the RPA result. The equations of motion for the Green functions neglected in Eq. (12) have been analyzed in Ref. 19. The consequences are that the RPA parameters are being replaced by effective ones and that the excitations become damped. One of the effective parameters J xy (q) replaces J(q) in J 1 and J 3 in all of the equations above and is 19 where η xy = 2 n 01 + n 03 The most important modification of J(q) is the constant shift introduced by a xy , but first we want to discuss the other term where χ xy (k, ω) is a generalized susceptibility. 19 The imaginary part of b xy (ω) determines the damping effects, which are, however, small at low temperatures. The renormalization of J xy (q) produced by the real part of b xy (ω) is somewhat smaller than the constant shift due to a xy , but it is not entirely negligible. Because of its moderate importance we have simplified the expression for the b xy (ω) term as follows We include only the real part and neglect modifications produced by the field. The quantity ǫ 0 is due to the finite lifetime of the excitations and we take it to be a constant, ǫ 0 ≈ 0.1 meV at zero temperature and field. The linewidths are going to increase rapidly when the temperature becomes comparable to ∆/k. We have accounted for this effect in a rough manner by scaling the result B q at zero temperature and field by the populationdependent factor in front, where n 0 µ is the value of a population number at zero field and temperature. A closer examination indicates that this simple scaling accounts for the increase of the linewidths in a reasonable way. The scale factor is unimportant for the analysis of the low temperature properties, but the power of 3 used in this expression ensures that the scale factor times η xy vanishes in the high-temperature limit.
Finally, the most important renormalization effect, the constant term (1 + η xy )a xy in J xy (q) in Eq. (25), is determined implicitly by The renormalized value of A 16 is equal to this sum over q times n 01 n 03 , and since the sum now vanishes, the condition a i 01 a i 03 = 0 is satisfied. The parameters J 1 and J 3 in Eqs. (12)-(16) are replaced by, respectively, n 01 J xy (q) and n 03 J xy (q), and, similarly, J 2 in Eq. (17) is replaced by n 02 J z (q), where the effective exchange coupling is Here η z = 1/n 2 02 − 1 and b z (ω) ≈ (n 02 /n 0 02 ) 3 B q . The constant term (1 + η z )a z is determined by the condition that a i 02 a i 02 = 0. Besides the modifications of the exchange couplings, the energy splitting ∆ in E 1 and E 3 , or in E 2 , is replaced by, respectively, 19 The renormalization of the energy-level separations has not much influence on the final calculations, and the small additional modifications derived in Ref. 19 might have been neglected. However, the leading order effects of these extra terms are actually included above by assuming 1 + η xy to be a factor n 0 + (n 1 + n 3 )/2 smaller than derived in Ref. 19, and, similarly, 1 + η z has been divided by n 0 + n 2 .

B. Comparison with experiments
Cavadini et al. 5 and Oosawa et al. 6 have measured the dispersion of the magnetic excitations in TlCuCl 3 at zero field in the zero-temperature limit (1. Here a 0 is the value of a xy or a z at zero temperature and field. The quantity B q is roughly proportional to J(q) and its averaged value with respect to q is zero. The RPA equations above have been solved numerically by an iterative procedure. The calculations benefit from the fact that all q summations may be parameterized in terms of J eff (q) [we neglect the minor difference between J(q) + η 0 xy B q and J(q) in Eq. (28) determining B q ]. Hence all summations may be expressed as integrals with respect to J eff (q) times a corresponding "density of states" calculated once and for all from Eq. (32). The value of ∆ eff used in the calculations is 5.671 meV, (almost) equal to the value of 5.68 meV derived by Oosawa et al. 6 Besides this parameter and those defining J eff (q) we have assumed that g = 2.06, which is the generally accepted value for g in the case where the field is applied along the b axis. 3,4,22 Finally, we have added the mean field from the parallel component to the ap- It is worthwhile to notice that although the population numbers of the excited levels are predicted to be nonzero at T = 0, Eq. (21)-(24), the quantum fluctuations do not give rise to any difference between n 1 and n 3 , i.e. m z = 0 at zero temperature as long as the system stays paramagnetic. This is consistent with the condition that i S iz = i (a i 11 − a i 33 ) commutes with the Hamiltonian. Using the model defined here we have calculated the susceptibility as a function of temperature. The result is compared with experiments in Fig. 1. The calculated critical field at which the paramagnetic phase becomes unstable is compared with experiments in Fig. 2. The induced magnetic moment m z at various values of the field has been calculated as a function of T . These results are shown in Fig. 3. This figure also includes results obtained in the ordered phase, which is considered in the following section.

III. EXCITATIONS IN THE ANTIFERROMAGNETIC PHASE
The paramagnetic phase becomes unstable when the energy of the lowest excitation vanishes. The lowest energy mode is the one with the energy E − q at q = Q = (001), and below the transition the expectation value of the dimer-spin variableS i becomes non-zero. The ordering is antiferromagnetic in the sense that S i have opposite signs on the two chain sublattices. The antiferromagnetic ordering may be transformed to the uniform one by an interchange of the two spins in the definition ofS i for the dimers belonging to, for instance, the 3-4 sublattices. The only effect of this transformation is that the interchain coupling between the 1-2 and 3-4 dimers changes sign, i.e. the ± in front of the second term of J(q) in Eq. (4) is being replaced by ∓, and J(Q) and J(0) are being interchanged within the extended zone scheme. That the ordering is antiferromagnetic instead of being uniform does not introduce any further complications. S i is perpendicular to the z-direction of the field, but its direction within the x-y plane is arbitrary (as long as any anisotropy is neglected). For convenience we shall make the choice that the ordered moment is along the x axis, and we define When S x is non-zero, the MF Hamiltonian for the "non-interacting" ith dimer becomes in terms of the standard basis operators of the paramagnetic system. We shall continue to label the eigenstates by |µ , where µ = 0, 1, 2, 3, and the ground state of this MF Hamiltonian may then be written The two angles θ and α minimize 0|H i MF |0 , or may be determined by demanding the off-diagonal terms of the MF Hamiltonian to vanish, and we get tan α = h cos 2θ ∆ cos 2 θ − h sin 2 θ sin 2α (38) and sin 2θ = 2J(Q) cos α ∆ cos 2α cos 2θ + 2 sin 2 θ sin 2 α S x .
(39) By calculating the state vectors of the excited MF levels to order θ 2 , we find that the order parameter is By h c we denote the critical field at which these equations have a non-zero solution for S x in the limit of θ → 0, and this field is found to be determined by If we replace J(Q) by J xy (Q) and ∆ by ∆ 1 this condition is the same as that derived from the requirement that the energy of the lowest paramagnetic excitation E − Q , within the self-consistent RPA, should vanish at the transition. The results above are more general but coincide with those derived in the zero-temperature mean-field theory of Matsumoto et al. 11 [Their angle φ corresponds to our π 4 − α. In our notation h = h 0 + J F (0) S z , whereas in their notation h in Eq. (38) should read h − J(Q) S z corresponding to the replacement of h 0 by h with their implicit assumption that the ferromagnetic interaction J F (0) is equal to −J(Q)].
When the MF Hamiltonian of the antiferromagnetic phase has been diagonalized we may proceed as in the paramagnetic case for calculating the correlation functions. The positions of the four different levels and the matrix elements ofS i may be calculated analytically if only terms to leading order in θ 2 are included. We are not going to present these results, since in the final calculations we chose the more accurate approach of diagonalizing the MF Hamiltonian numerically. In terms of the standard basis operators of the final MF Hamiltonian,S jx is now [m 3 (a j 30 + a j 03 ) − m 1 (a j 10 + a j 01 )]/ √ 2, where m 1 and m 3 are different from 1 and from each other. Here we neglect the extra complication that the matrix elements between the excited states, as for instance 1|S jx |3 ∝ θ, become non-zero. These excited state contributions to the correlation functions get multiplied by n 13 , hence they are unimportant not only when the ordered moment is small, but in most of the regime where the RPA modifications of the MF behavior are of importance.
The final Hamiltonian may be written in the same way as in the disordered case, Eq. (7). The positions of the three excited levels are being shifted and J(ij) is being multiplied by different factors depending on which operator product is considered, and, finally, the remaining off-diagonal products, a i 01 a j 01 , a i 01 a j 30 etc., now appear. The new off-diagonal contributions are all of order θ 2 and affect the diagonal correlation functions only to order θ 4 . Hence, to leading order these extra contributions may be neglected. In this case the matrix equations once again decouple into 3 sets of 2 × 2 equations, which may be solved analytically. The result for the population numbers is the same as that given by Eqs. (21)- (24), except that ∆ and ∆ ± h are being replaced by the energies of the three corresponding MF levels and that J 1 , J 2 , and J 3 are being multiplied by matrix-element factors, which are slightly different from 1 and from each other. The equivalence implies that the modifications of the RPA correlation functions may be calculated as in the paramagnetic case, and, for instance, a xy is still determined by Eq. (29) except that the expressions for E ± q are being modified. We used this approximation, valid in the limit of θ 2 being small, for calculating the magnetization curves. The results were close to those shown in Fig. 3, however, in the final calculations we included the higherorder modifications.
For a given set of population numbers the MF Hamiltonian was diagonalized numerically, determining all possible matrix elements and the four energy levels. The knowledge of the population numbers and the matrix elements is used for determining the two expectation values S x and S z , and for constructing the total Hamiltonian expressed in terms of the standard basis operators. When the interactions between the excited states are neglected, the equations of motion lead to a 4×4 set of matrix equations for the xy part and a 2 × 2 set for the longitudinal Green functions. The two sets of equations were inverted analytically (utilizing Mathematica for handling the set of 4 × 4 equations). In this way we derived an explicit expression for the correlation function matrixĀ(q). The averaged values of the diagonal components were used for calculating the population numbers as in the paramag- netic case. The renormalization parameter a xy and a z are determined by, respectively, A 16 = 0 and A 34 = 0. The condition A 16 = 0 also implies that A 25 = 0, but not necessarily that the new off-diagonal components vanish. We have neglected the possibility that the renormalization of the additional off-diagonal exchange terms might be different, since the new terms are small whenever the renormalization effects are important. Except for the matrixelement modification of the exchange terms we use the same approximate expression, Eq. (28), for the renormalization parameter b xy (ω), and similarly for b z (ω). When the renormalization parameters have been determined we may calculate the renormalized value of J(Q) in the MF Hamiltonian (36), which is being replaced by J xy (Q) determined from Eq. (26). Similarly ∆ ± h and ∆ in Eq.
(36) are replaced by, respectively, ∆ 1 ±h and ∆ 2 given by Eq. (31). The MF Hamiltonian is then consistent with the renormalized RPA expressions, and the whole procedure has been carried out in a self-consistent manner, so that the population numbers assumed as a start are the same as those derived. This theory was used for calculating the properties of the dimer system in the ordered phase. The temperature dependence of the parallel moment in a constant applied field is shown in Fig. 3. In Fig. 4 we show the calculated energies for the three modes at (001) as functions of field at 1.5 K compared with the neutron-scattering results of Rüegg et al. 25,26 The experimental results indicate that ∆ for the longitudinal mode, the energy of which is nearly unaffected by the field in the paramagnetic phase, is slightly larger than for the two other modes, and we have accounted for this effect by adding 0.03 meV to ∆ 2 in Eq. (31). This modification does not affect the xy- polarized modes, and as discussed by Matsumoto et al., 11 the lowest-energy xy-polarized mode becomes the Goldstone mode in the ordered phase, the energy of which depends linearly on |q − Q| and is zero at the ordering wave vector. 26 The spin-resonance experiments of Glazkov et al. 22 indicate that the Goldstone mode develops an energy gap for fields larger than the critical one corresponding to the presence of a small anisotropy within the x-y plane of the same order of magnitude as the anisotropy considered above. The two anisotropy terms are unimportant for the renormalization effects and are neglected elsewhere in the present calculations. The final figure, Fig. 5, shows the field dependencies of the squared primary order parameter m 2 xy and of the parallel magnetization m z in the zero-temperature limit. The comparison with experiments shows that the self-consistent RPA accounts reasonably well for the field dependence of m z , whereas the order parameter is calculated to increase rather faster with field than observed. The corresponding MF model, on the other hand, underestimates the value of m z , but predicts an order parameter which is close to the one observed.

IV. DISCUSSION
The present dimer system is unique because it clearly exhibits the importance of quantum fluctuations. The system is close to a quantum critical point, and a zero- temperature phase transition may be achieved either by the application of a modest magnetic field of 54 kOe 3 or a hydrostatic pressure of 1.1 kbar. 12,13 Here we have considered the case where the transition is approached by applying a magnetic field at ambient pressure. The field removes the degeneracy of the S = 1 triplet states of the dimers, and the collective singlet-triplet excitations separate into a longitudinal, z-polarized wave and two transverse modes. The collective transverse modes are linear combinations of propagating modes due to transitions between the ground state and the lowest and the highest excited states of the single dimers. In the paramagnetic phase, the two transverse modes are subject to the same renormalization effects, because the rigid energy shift of the excitations, E + q − E − q = 2h, does not influence the quantum fluctuations and n 13 = 0 at T = 0. At non-zero temperatures, the thermal fluctuations imply that n 13 is non-zero, but the renormalized exchange interaction J xy (q) is still the same for the upper and lower transverse modes in the paramagnetic phase (within the present approximation scheme).
The RPA renormalization parameters are determined to be rather substantial in the zero-temperature limit at zero field, because the system is close to the critical point. The occupation number of the dimer ground state is calculated to be n 0 0 = 0.935 and a 0 = 0.393 meV. The constant reduction of the exchange interaction is (1 + η 0 xy )a 0 = 0.472 meV, which is about 17% of the maximum value of the effective exchange interaction J eff (Q) = 2.8 meV. The effective singlet-triplet splitting is about ∆ + a 0 = 5.26 meV in the temperature range of the maximum in the susceptibility. As indicated by the comparison in Fig. 1, this is in good agreement with that derived from the experimental data. 3,4 Hence, the theory is able to account for the difference found ex-perimentally between the (effective) energy gap of 5.68 meV, derived from the excitation spectrum in the T = 0 limit, 6 and the smaller value of the gap determined from the susceptibility measurements. At zero temperature the renormalization parameters are independent of the field as long as it stays smaller than the critical one. At non-zero temperature the lower branch E − q is much more easily populated than the corresponding MF level. For comparison, the bulk magnetization at 53 kOe predicted by the corresponding MF model is about a factor of 50 times smaller at 6 K. At a constant non-zero temperature E − q decreases, and the thermal population of the E − qbranch increases, when the field approaches the critical one. This also implies that the renormalization parameter a xy increases with increasing field. When the field becomes larger than the critical one, a xy is reduced as the field is further increased. This reduction is so large that it is able to stabilize the ordered moment also at fields smaller than the critical one, i.e. the transition is so strongly modified that it becomes a first order one. This behavior of the renormalization parameters at the phase transition is illustrated in Fig. 6. The unphysical enhancement of the renormalization effects near a phase transition is a quite general feature of the self-consistent version of RPA. In the close neighborhood of the critical point the renormalization effects show a pronounced sensitivity to small modifications of the model, and the rather good agreement between theory and experiments obtained here for the critical field and for the parallel magnetization, see Figs. 2 and 3, is somewhat fortuitous. Even the minor change introduced by using the exchange parameters derived by Oosawa et al. 6 rather than those defined by Eq. (32) leads to a relative increase of a xy by about 10% and to corresponding changes of the magnetization curves and the critical field (see Fig. 2).
Focusing our attention on the zero-temperature limit, then n 13 is zero when the field is smaller than the critical one, but becomes non-zero in the ordered phase. As long as the ordered moment is small, the contributions due to n 13 may be neglected and the equations determining S x and S z , Eqs. (40) and (41), predict when terms of the order of θ 4 and α 2 θ 2 are omitted. The equation is only weakly influenced by the renormalization effects since n 01 +n 03 ≈ 2 at T = 0. Nevertheless, the experimental low-temperature results are far from obeying this relationship, which circumstance makes it difficult to reproduce the field dependencies of the two magnetization components simultaneously, as illustrated by  Fig. 4, should be multiplied by ∆/h ≃ 2.6 leading to a slope of m 2 xy with respect to field, which is nearly twice the one derived by the present theory, i.e. a factor of 3-4 larger than the experimental one (see Fig. 5). In the paramagnetic phase, the lowest excited state of a single dimer is the |1 + 1 state with S z = 1, and Nikuni et al. are assuming that this state is the one determining the wave functions of the lowest lying mode of collective excitations, and hence the one which defines the condensate in the ordered phase. This corresponds to assuming α = π/4 in our Eq. (37). However, as also stressed by Matsumoto et al., 11 it is crucial to include the presence of the |1 − 1 level in order to get a consistent description of the excitations and of the condensate. In the paramagnetic phase, the matrix element ofS x between the ground state |00 and the lowest excited state |1 + 1 is numerically the same as its matrix element between |00 and |1 − 1 . The same applies toS y , and this means that the collective transverse excitations transmitted via these two operators are mixed S z = ±1 excitations. Being proportional to J(q), the degree of mixing depends on the wave vector and is at its maximum at the ordering wave vector. In correspondence to this, the MF ground state in the ordered phase, Eq. (37), involves |1 − 1 as well as |1 + 1 . The two states are of equal importance in the limit of zero field, and the relative weight of the two states is shifted from 1 in the presence of a field as described by the angle α ≃ h/∆. The two components depend differently on α, e.g. S z = 0 whereas S x has its maximum at α = 0, and the factor h/∆ in Eq. (43) is a simple consequence of this difference.
The present self-consistent RPA theory accounts reasonably well for the paramagnetic properties of the dimer system. Within the MF model, the bulk susceptibility vanishes exponentially in the zero-temperature limit, whereas the present RPA model predicts a power law m z /H ∝ T φ with φ = 1.8 at H = 53 kOe. This is consistent with experiments and the RPA theory also predicts the right critical field for the phase transition. In contrast to the boson model of Eq. (1), the present theory does not rely on any free parameters. Note, however, that if we use the effective exchange parameters, which Oosawa et al. determined from their measurements of the dimer excitation spectrum, 6 the critical field is derived to increase slightly faster with temperature than observed, as shown in Fig. 2. This minor discrepancy was neutralized by a small adjustment of the density of states of J eff (q). Actually, it would have been a surprise if the self-consistent RPA theory had been able to predict the right critical field without any adjustments. In all circumstances, it is clear that the theory needs to be corrected due to critical fluctuations, since the self-consistent RPA predicts the phase transition to the antiferromagnetic phase to be of first order in contradiction with experiment.
The dimer system has a number of unusual magnetic properties. The most outstanding one is that the system is driven into the phase of antiferromagnetic order by the application of a uniform field. Another unusual property shown by the ordered phase is that the bulk magnetization increases, when the temperature is lowered at a constant field, as this happens in spite of the fact that the "degrees of freedom" are being reduced because of the accompanying enhancement of the antiferromagnetic order parameter. This behavior is not in accordance with the MF model, whereas the self-consistent RPA theory accounts, at least, qualitatively for this observation. The self-consistent theory of the ordered phase is complicated, and we have been forced to neglect a number of effects. One of the complications, which has not been mentioned above, is that the matrix elements of S i between the ground state and the excited states are no longer zero in the ordered state implying additional modifications of all normal modes of the system (the two classes of modes may mix because the Cu sites lack inversion symmetry). We may also add that the non-zero values of the diagonal elements ofS x in the ordered phase effectively give rise to additional contributions to ∆ 1 and ∆ 3 . Fortunately, these extra complications should be unimportant within the regime of low temperatures and small order parameter, where the theory is applied. The approximations made for b xy (ω) in Eq. (28) are only acceptable, when the renormalization effects due to this term are small. This term is handled in a more rigorous way by a diagrammatic high-density expansion. 28 To first order in 1/z (where z is the number of interacting neighbors), the diagrammatic theory may be formulated in a self-consistent, "effective medium" fashion, which is the equivalent of the present, self-consistently improved version of RPA. 17,29 This 1/z theory has been applied to the Ising systems HoF 3 and LiHoF 4 , 30,31 and in these cases the phase transitions are predicted to remain of second order. Therefore, we expect that the use of the 1/z expansion theory, to first order in 1/z, is going to reproduce the self-consistent RPA results derived here for the paramagnetic phase, and should lead to an improved description of the ordered phase.
The present RPA theory establishes a classification of the renormalization effects, which affect the properties of this quantum-critical dimer system. The most important one is the constant reduction of the exchange interaction by (1+η xy )a xy . This term is equivalent to that produced by an on-site repulsive interaction J xy (ii) and is here found to be determined directly from the higher-order modifications of the RPA Green functions. Our analysis of the spin model also predicts the presence of other renormalization effects, such as the ω-dependent correction b xy (ω) to J xy (q) and the increase of the effective splitting between the single-dimer energy levels by a xy or a z . The q-independent reduction of the exchange interaction has its parallel in the phenomenological repulsive interaction, v 0 in Eq. (1), in the boson model of Nikuni et al., 8 whereas the two other renormalization effects have no counterpart in their theory. A more problematic sim-plification made by Nikuni et al. is their assumption that the low temperature properties of the system are dominated by one type of bosons, whereas, in reality, the system contains three different kinds, where those corresponding to S z = 1 and S z = −1 are mixed. The degree of mixing depends on wave vector and on field, and is important for the characterization of the bosons in the condensate.
One basic difficulty in the many-body theory of localized spin systems is that the operators describing the dynamics of the single spins are not bosonic but more complicated operators, as indicated by Eq. (11). This complication is responsible for the need to renormalize the simple RPA theory. Although the present self-consistent theory includes the leading-order renormalization effects, the comparison between theory and experiments within the ordered phase of TlCuCl 3 is not satisfactory. The diagrammatic 1/z theory, 17,28,29 represents a more systematic approach and should be able to give a more acceptable description of the ordered phase. We expect, however, that the pronounced experimental violation of the MF/RPA relation Eq. (43), between the bulk magnetization and the ordered antiferromagnetic moment, will remain a challenge to future theory.