Continuous unitary transformations in two-level boson systems

PHYSICAL REVIEW C 72, 064332 (2005) Continuous unitary transformations in two-level boson systems S´ ebastien Dusuel,1,∗Julien Vidal,2,†Jos´ eM.Arias, 3,‡Jorge Dukelsky,4,§and Jos´ e Enrique Garc´ ıa-Ramos5, 1Institut f¨ ur Theoretische Physik, Universit¨ at zu K¨ oln, Z¨ ulpicher Str. 77, D-50937 K¨ oln, Germany 2Laboratoire de Physique Th´ eorique de la Mati` ere Condens´ ee, CNRS UMR 7600, Universit´ e Pierre et Marie Curie, 4 Place Jussieu, F-75252 Paris Cedex 05, France 3Departamento de F´ ısica At´ omica, Molecular y Nuclear, Facultad de F´ ısica, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain 4Instituto de Estructura de la Materia, CSIC, Serrano 123, E-28006 Madrid, Spain 5Departamento de F´ ısica Aplicada, Universidad de Huelva, E-21071 Huelva, Spain (Received 28 September 2005; published 30 December 2005) Two-level boson systems displaying a quantum phase transition from a spherical (symmetric) to a deformed (broken) phase are studied. A formalism to diagonalize Hamiltonians with O(2L+1) symmetry for large number of bosons is worked out. Analytical results beyond the simple mean-field treatment are deduced by using the continuous unitary transformations technique. In this scheme, a 1/N expansion for different observables is proposed and allows one to compute the finite-size scaling exponents at the critical point. Analytical and numerical results are compared and reveal the power of the present approach to compute the finite-size corrections in such a context. DOI: 10.1103/PhysRevC.72.064332 PACS number(s): 21.60.Fw, 21.10.Re, 05.10.Cc, 75.40.Cx I. INTRODUCTION The study of two-level systems has been a topic of interest since the first steps in the development of quantum mechanics. The main advantage of these models is that they can be numerically diagonalized for very large dimensions and, at the same time, they can model realistic quantum many-body systems. Typical examples are the Jaynes-Cummings model of quantum optics [1], the vibron model (VM) of quantum chemistry [2], the two-level pairing model in condensed matter [3] and in nuclear physics [4], the Lipkin-Meshkov-Glick model (LMG) [5–7], and the interacting boson model (IBM) [8] of nuclear structure. Although some of these models describe two-level fermion systems, the model Hamiltonian can always be written in terms of SU(2) pseudospin operators. Subsequently, the spin Hamiltonian can be expressed in terms of bosons using either the finite Schwinger representation or the infinite Holstein-Primakoff representation of the SU(2) algebra. An example is the LMG model, which has recently been newly revived as a model of quantum spins with long-range interactions to investigate the relationship between entanglement and quantum phase transitions (QPTs) [9–15]. In its boson representation, it has also been recently used as a simplified model to describe the Josephson effect between two Bose-Einstein condensates [16]. In this work, we focus on finite two-level boson Hamiltonians, having the common feature that the lower level is always a scalar L=0 boson, hereafter written as sboson. The upper level can have different multipolarities, generically noted as L, whose value defines a particular model. The LMG model ∗Electronic address: [email protected] †Electronic address: [email protected] ‡Electronic address: [email protected] §Electronic address: dukelsk[email protected] Electronic address: [email protected] in the Schwinger representation has a second scalar (L=0) boson for the upper level. A dipolar (L=1) boson leads to the VM and a quadrupolar (L=2) boson corresponds to the IBM. Higher angular momentum bosons can lead new models, such as, for example, a model of octupole vibrations in terms of s and octupolar (L=3) bosons. A schematic representation of the model is shown in Fig. 1. All these two-level boson models are governed by an algebraic structure that is constructed of all bilinear combinations of creation and annihilation boson operators that generate the algebra of U(2L+2). One of the main features of these models is that one can construct dynamical symmetries in which the Hamiltonian can be written in terms of Casimir invariants of a nested chain of subalgebras of U(2L+2). In these particular cases, the problem is analytically solvable providing a powerful tool to check approximate methods and a reference for more detailed calculations. In addition, if the Hamiltonian is written as a linear combination of Casimir invariants of the subalgebras U(2L+1) and O(2L+2) the model is still quantum integrable but then requires solving Bethe-like equations numerically [17]. The exact solution, given by Richardson almost 40 years ago [18], is reduced to a set of Mnonlinear coupled equations, where Mis the number of boson pairs. For two-level boson models it turns out that the numerical diagonalization of the Hamiltonian presented below is more efficient than solving the Richardson equations. The aim of this work is to study the QPT that occurs in the two-level boson system as it evolves from the spherical vibrational U(2L+1) symmetry to the deformed O(2L+2) symmetry, as a function of a control parameter. Although, strictly speaking, QPTs are defined for macroscopic systems, there is a renewed interest in studying structural changes in finite-size systems as the precursors of a QPT in the thermodynamic limit. Traces of these QPTs are readily observed in finite systems and their properties are then correlated with the idealized thermodynamic system [19]. The understanding of 0556-2813/2005/72(6)/064332(17)/$23.00 064332-1 ©2005 The American Physical Society DUSUEL, VIDAL, ARIAS, DUKELSKY, AND GARC´ IA-RAMOS PHYSICAL REVIEW C 72, 064332 (2005) L=0 L x L=0 (LMG), L=1 (VM), L=2 (IBM), . . . FIG. 1. Schematic representation of the two-level boson model studied in this article. the modifications on the characteristics of the QPT induced by finite-size effects is of crucial importance to extend the concept of phase transitions to finite systems. Several techniques have long been used to extrapolate numerical results obtained by large-scale diagonalizations or Monte Carlo calculations to the infinite system. Here, we focus on a somewhat inverse problem, which is the finite-size corrections to the observables in two-level boson models such as the ground-state energy, the gap, the occupation number, and some electromagnetic transition rates. Although the zeroth order in the boson number Nis given by the Hartree mean-field approach for the ground state and the random phase approximation for the excited states, going beyond this order implies the use of more sophisticated techniques. We make use here of the continuous unitary transformations (CUTs) and give the first 1/N corrections to the observables in the whole U(2L+1) to O(2L+2) transition. The structure of the article is the following. In Sec. II we introduce the two-level boson models and the formalism for the numerical diagonalization for a very large number of bosons. Section III describes the mean-field treatment of the two-level boson models. Section IV is devoted to the study of the symmetric phase using CUTs. Analytical expressions for different orders in the 1/N expansion of the ground-state energy, the gap, the expectation value of the number of Lbosons in the ground state, and the transition matrix element between the ground state and the first excited state are obtained. In Sec. V the broken phase is analyzed, and in Sec. VI the study of the critical point is presented from the spherical phase. In this section, we obtain the finite-size scaling exponents for the quantities cited above by analyzing the divergence of their 1/N expansion. In Sec. VII a comparison of the numerical results obtained using the formalism presented in Sec. II with the analytical CUTs results is presented. Section VIII is for summary and conclusions. Technical details concerning flow equations can be found in appendices. II. TWO-LEVEL BOSON MODELS In this section, we present a simple algorithm for diagonalizing boson Hamiltonians with O(2L+1) symmetry for large boson numbers. The formalism is based on previous studies [20,21] and is a generalization of the one presented recently for treating the IBM [22]. We consider the following boson pairing Hamiltonian H=xnL+1−x 4(N−1)(P† LPL+P† sPs−P† LPs−P† sPL),(1) with nL=+L  µ=−L L† µLµ,(2) P† s=(Ps)†=s†s†,(3) P† L=(PL)†=+L  µ=−L (−1)µL† µL† −µ.(4) L† µcreates a boson in the excited Llevel with projection µ, whereas Lµdestroys it. We have introduced above the pair P† ρoperators with ρ=sor L, which is used later on. The total number of bosons N=ns+nLis a conserved quantity. For x=0,H can be cast into a linear combination of the quadratic Casimir operators of O(2L+2) and the corresponding subalgebras, whereas for x=1,H is the linear Casimir operator of the U(2L+1) algebra. Here, xplays the role of a control parameter, driving the system from the O(2L+2) deformed phase to the U(2L+1) spherical phase. The boson pairing Hamiltonian (1) can be studied in an elegant way by means of the noncompact SU(1,1) algebra of boson pair operators. For the subspace of ρbosons, where ρstands generically either for sor Lbosons, the SU(1,1) generators are the raising operator K+ ρ, the lowering operator K− ρ=(K+ ρ)†and the Cartan operator K0 ρdefined as K+ ρ=1 2P† ρ,(5) K− ρ=1 2Pρ,(6) K0 ρ=1 2 µρ† µρµ+1 2=1 2nρ+1 4Dρ,(7) with Dρ=2ρ+1. The three operators {K+ ρ,K− ρ,K0 ρ}satisfy the SU(1,1) commutator algebra K0 ρ,K± ρ=±δρ,ρK± ρ,(8) [K+ ρ,K− ρ]=−2δρ,ρK0 ρ.(9) The Casimir operator is C2 ρ=1 2(K+ ρK− ρ+K− ρK+ ρ)−(K0 ρ)2=−Dρ 4(Dρ 4−1). (10) The complete set of eigenstates of the pairing Hamiltonian (1) can be constructed as a direct product of subspaces associated to sand Lbosons. Each subspace can be written in terms of the rising operators K+ ρacting on the corresponding subspace of unpaired ρbosons |˜ nρ,ν ρ= 1 C˜ nρ ρ,νρ K+˜ nρ ρ|˜ nρ=0,ν ρ,(11) where νs=νL=0=0,1 and νL=0=0,1,2,... The quantity νρis known as the boson seniority for ρbosons and gives the number of bosons of type ρnot coupled in pairs to zero. Note that from now on the label ˜ nmeans number of boson 064332-2 CONTINUOUS UNITARY TRANSFORMATIONS IN TWO- . . . PHYSICAL REVIEW C 72, 064332 (2005) pairs coupled to zero angular momentum. The total number of bosons is 2˜ ns+2˜ nL+νs+νL. The normalization constant in Eq. (11) can be obtained from the action of K− ρand K0 ρon the ρsubspace |˜ nρ=0,ν ρ K− ρ|˜ nρ=0,ν ρ=0,(12) K0 ρ|˜ nρ=0,ν ρ=1 2νρ+1 4Dρ|˜ nρ=0,ν ρ,(13) and the commutation relation [[K− ρ,K+ ρ],K+ ρ]=2K+ ρ,(14) then K− ρ(K+ ρ)˜ nρ|˜ nρ=0,ν ρ =˜ nρ(˜ nρ+Dρ 2+νρ−1)(K+ ρ)˜ nρ−1|˜ nρ=0,ν ρ,(15) and ˜ nρ,ν ρ|˜ nρ,ν ρ= ˜ nρ 2(2˜ nρ+2ρ+2νρ−1) ×˜ nρ−1,ν ρ|˜ nρ−1,ν ρ,(16) and finally ˜ nρ,ν ρ|˜ nρ,ν ρ= ˜ nρ!(2˜ nρ+2ρ+2νρ−1)!! 2˜ nρ(2ρ+2νρ−1)!! =C˜ nρ ρ,νρ.(17) Remember that the label ρstands for sor Lbosons and takes the following numerical values: 0 for sbosons and L for Lbosons. Once the basis for each subspace is obtained [Eqs. (11) and (17)], the complete basis set for the pairing Hamiltonian [Eq. (1)] is easily constructed, |˜ ns,ν s;˜ nL,ν L= 1 C˜ ns s,νsC˜ nL L,νL (K+ s)˜ ns(K+ L)˜ nL ×|˜ ns=0,ν s;˜ nL=0,ν L.(18) We now diagonalize the Hamiltonian given in Eq. (1) in the basis of the states given by Eq. (18). Note that in the following nrefers to boson number operators but ˜ nare number of boson pairs. The relevant matrix elements for the construction of the Hamiltonian matrix are as follows: ˜ ns,ν s;˜ nL,ν L|ns|˜ ns,ν s;˜ nL,ν L=2˜ ns+νs,(19) ˜ ns,ν s;˜ nL,ν L|nL|˜ ns,ν s;˜ nL,ν L=2˜ nL+νL,(20) ˜ ns,ν s;˜ nL,ν L|K+ sK− s|˜ ns,ν s;˜ nL,ν L=˜ ns˜ ns+νs−1 2, (21) ˜ ns,ν s;˜ nL,ν L|K+ LK− L|˜ ns,ν s;˜ nL,ν L=˜ nL˜ nL+L+νL−1 2, (22) (˜ ns−1),ν s;(˜ nL+1),ν L|K+ LK− s|˜ ns,ν s;˜ nL,ν L =1 2˜ ns(˜ nL+1)(2˜ ns+2νs−1)(2˜ nL+2L+2νL+1). (23) It is clear that the Hamiltonian does not mix states with different boson seniority quantum numbers. Thus, the Hamiltonian matrix is block diagonal. In addition, within one block the matrix is tridiagonal, making the diagonalization simple. The different states are obtained as follows: one starts with the boson subspace containing N/2 boson pairs coupled to zero angular momentum, νs=νL=0. The diagonalization of Hin this subspace provides states with angular momentum zero and the first eigenstate is the ground state. Next, one goes to the block with one broken boson pair. This block is composed of two separate blocks because the two bosons can be one sboson and one Lboson (νs=1,ν L=1) or two Lbosons not coupled to zero because the coupling to zero is included in the first block, (νs=0,ν L=2). Notice that two unpaired sbosons are not possible because they are always coupled to zero and consequently they are counted in the first block. For the case of the LMG model in which L=0, the block νs=0,ν L=2 is not allowed for the same reason. The block νs=1,ν L=1 provides states with angular momentum L, the first of them is the first excited state of the system. The block νs=0,ν L=2 contains states with two Lbosons not coupled to zero angular momentum; it contains angular momenta: 2L, 2L−2,...,2. Next, there is another block with two broken boson pairs composed of two subblocks: νs=1,ν L=3 and νs=0,ν L=4. Again, the block νs=0, νL=4 is absent for the LMG model. This construction continues for 3,4,...,N/2 broken boson pairs. Few first low-lying Hamiltonian eigenstates are depicted schematically in Fig. 2. Direct block diagonalization of the Hamiltonian in the basis of the states given by Eq. (18) provides observables as the ground-state energy or the gap and also the wave functions of the states. With the wave function of the ground state the expectation value of the number of Lbosons in the ground state can be easily calculated. One can also calculate the transition probability from the ground state to the first excited state provided with the appropriate operator. The natural transition operator for the pairing Hamiltonian we are considering is as follows: TLµ=L† µs+(−1)µs†L−µ.(24) N 2boson pairs ν = 0; ν = 0 sL ν =1; ν =1 sL ν =0; ν =2 sL ν =1; ν =3 sL ν =0; ν =4 sL 1 broken boson pair 2 broken boson pairs . . . 0 0 0 L L . . . . . . . . . . . . 2L,2L−2,...,2 2L,2L−2,...,2 3L,3L−2,...,0 3L,3L−2,...,0 4L,4L−2,...,0 4L,4L−2,...,0 FIG. 2. Schematic representation of the level sequence obtained by diagonalizing the Hamiltonian for the two-level boson models studied in this article. Numbers above the lines are angular momenta. 064332-3 DUSUEL, VIDAL, ARIAS, DUKELSKY, AND GARC´ IA-RAMOS PHYSICAL REVIEW C 72, 064332 (2005) The action of TLµon the subspace νs=0,ν L=0 that includes the ground state is as follows: TLµ|˜ ns,0; ˜ nL,0 =1 C˜ ns s,0C˜ nL L,0 [˜ ns(K+ s)˜ ns−1(K+ L)˜ nLs†L† µ+˜ nL ×(K+ s)˜ ns(K+ L)˜ nL−1s†L† µ]|˜ ns=0,0; ˜ nL=0,0,(25) where the state |˜ ns=0,ν s=0; ˜ nL=0,ν L=0≡|0) is the boson vacuum. Then, the matrix elements of TLµconnecting the subspaces νs=0,ν L=0 and νs=1,ν L=1 (which includes the first excited state) are as follows: (˜ ns−1),ν s=1; ˜ nL,ν L=1|TLµ|˜ ns,ν s=0; ˜ nL,ν L=0 =2˜ ns(2˜ nL+2L+1) 2L+1,(26) ˜ ns,ν s=1; (˜ nL−1),ν L=1|TLµ|˜ ns,ν s=0; ˜ nL,ν L=0 =2˜ nL(2˜ ns+1) 2L+1.(27) If we write each eigenstate of the Hamiltonian as |i,ν s,ν L= ˜ ns,˜ nL cνs,νL ˜ ns,˜ nL|˜ ns,ν s;˜ nL,ν L,(28) the matrix element of the TLµ operator between the ground state |0≡|0,0,0and the first excited state |1≡|0,1,1 is 1|TLµ|0= ˜ ns,˜ nL2˜ ns(2˜ nL+2L+1) 2L+1c0,0 ˜ ns,˜ nLc1,1 ˜ ns−1,˜ nL +2˜ nL(2˜ ns+1) 2L+1c0,0 ˜ ns,˜ nLc1,1 ˜ ns,˜ nL−1.(29) The formalism presented in this section provides the exact full solution of the problem. A simpler approach to study ground-state properties in the large Nlimit is provided by the mean-field analysis presented in the next section. This limit is a good benchmark to test more elaborate results. III. MEAN-FIELD ANALYSIS The geometrical interpretation of the Hamiltonian [Eq. (1)] can be obtained by introducing a Hartree axial coherent state that allows us to associate to it a geometrical shape in terms of a deformation variable β. For a system with Nbosons, this state is obtained by acting Ntimes with a condensed boson † on the boson vacuum |0) |N,β= 1 √N!(†)N|0),(30) where the basic condensed boson operator is given by the following: †=1 1+β2(s†+βL† 0),(31) 012 β 0.1 0.2 0.3 0.4 0.5 e0 (arbitrary units) x=0.4 x=0.5 x=0.6 FIG. 3. (Color online) Energy surfaces per boson for the Hamiltonian [Eq. (1)] for different values of the control parameter x. which depends on the βshape variable. The energy surface is defined as the expectation value of Hin the intrinsic state E(N,β)=N,β|H|N,β =Nxβ2 1+β2+1−x 41−β2 1+β22.(32) Minimizing the variational energy E(N,β) with respect to βleads to a critical point at xc=0.5. For x>x c(symmetric phase), the ground state is spherical and is obtained for β=0, whereas for x<x c(broken phase) it is deformed since the minimum of the energy per boson e0=E(N,β)/N is obtained for β=√1−2xas can be seen in Fig. 3. At the critical point, it is worth noting that the energy surface is a flat β4surface near β=0 [22,23]. Within this mean-field (variational) approach, one thus gets the following ground-state energy per boson: e0(x⩾xc)=1−x 4,(33) e0(x⩽xc)=x 4 2−3x 1−x.(34) One can also straightforwardly compute the expectation value of nLin the ground state, which is found to vanish in the symmetric phase and equals nL=N1−2x 2(1 −x),(35) in the broken one. However, other properties, as excitation energies or transition probabilities, that imply excited states require to go one step beyond this simple mean-field level. In the following, we use a combination of several methods already detailed for the simple case L=0 in Ref. [24] that allow us to compute the corrections to these mean-field results as well as the gap or the transition rates that require the knowledge of excited states. 064332-4 CONTINUOUS UNITARY TRANSFORMATIONS IN TWO- . . . PHYSICAL REVIEW C 72, 064332 (2005) IV. THE SYMMETRIC PHASE (1/2<x<1) The starting point of our analysis is the elimination of the sboson by means of a Holstein-Primakoff boson expansion [25] of sand Lbosons (for a review in boson expansion techniques see Ref. [26]). Therefore, we introduce a set of bµbosons such that the mapping L† µLν=b† µbν,(36) L† µs=N1/2b† µ(1 −nb/N)1/2=(s†Lµ)†,(37) s†s=N−nb.(38) fulfills the commutation relations at each order in Nin the Taylor expansion of the square root. With these notations, we have the following: nL=nb,(39) P† LPL=P† bPb,(40) P† sPs=(N−1)(N−2nb)+:n2 b:,(41) P† LPs= µ (−1)µL† µsL† −µs, =N µ (−1)µb† µ(1−nb/N)1/2b† −µ(1−nb/N)1/2,(42) where: A: denotes the normal-ordered form of the operator A. The Holstein-Primakoff mapping eliminates the sboson at the cost of introducing infinitely many boson terms. However, each term in the expansion has a definite 1/N order. As shown in the preceding section, for 1/2<x<1, the number of Lboson in the ground state goes to zero in the thermodynamical limit. Thus, to capture the finite Ncorrections, one performs a 1/N expansion of the Hamiltonian considering that nb/N 1. In this phase, the Hamiltonian [Eq. (1)] reads as follows: H=N11−x 4+N01−x 4 ×2(3x−1) 1−xnb−(P† b+Pb)+N−11−x 4 ×:n2 b:+P† bPb−1 2(P† b+Pb)+P† bnb+nbPb +N−21−x 4 ×:n2 b:+P† bPb−3 8(P† b+Pb)+P† bnb+nbPb +O(1/N3).(43) Here, we have restricted this expansion to order (1/N)2but the method we used can, in principle, be applied beyond this limit as shown in Ref. [24] for the LMG model (L=0). Our aim is to diagonalize Horder by order. At the leading order, one obviously recovers the mean-field ground-state energy per boson e0(N)=(1 −x)/4+O(1/N). At order (1/N)0, the Hamiltonian is quadratic and it can thus be easily diagonalized through a Bogoliubov transform, giving rise to the boson random phase approximation formalism presented in Ref. [27] and more recently exploited to describe the properties of the symmetric and broken-symmetry phases in the interacting boson model [21,28]. Higher-order terms cannot be diagonalized by a Bogoliubov transformation, so that one has to resort a more sophisticated method. A. CUTs formalism The CUTs technique has been conjointly proposed by Wegner [29] and Głazek and Wilson [30,31]. For a pedagogical introduction, we refer the reader to Refs. [32,33]. Here we only sketch the main lines of this simple and powerful approach. The idea of the CUTs is to diagonalize the Hamiltonian in a continuous way starting from the original (bare) Hamiltonian H=H(l=0). A flowing Hamiltonian is thus defined by the following: H(l)=U†(l)H(0)U(l),(44) where lis a scaling parameter such that H(l=∞)is diagonal, and U(l) is a unitary transformation, i.e., satisfying U(l)U†(l)=U†(l)U(l)=1. Taking the derivative of Eq. (44) with respect to lyields the differential (flow) equation ∂lH(l)=[η(l),H(l)],(45) where the generator of the unitary transformation η(l)is η(l)=∂lU†(l)U(l)=−U†(l)∂lU(l).(46) CUTs are also a powerful tool to compute the expectation value of any observable . As for the Hamiltonian, we define a flowing operator (l)=U†(l)(0)U(l),(47) which obeys ∂l(l)=[η(l),(l)],(48) with =(l=0). The expectation value of on an eigenstate |ψof His then given by the following: ψ||ψ=ψ|U(l=∞)(l=∞)U†(l=∞)|ψ,(49) where U†(l=∞)|ψis simply the eigenstate of the diagonal Hamiltonian H(l=∞). The key point of this approach is an appropriate choice of the generator η, which, in fact, depends on the problem under consideration. Here, the Hamiltonian Hexpressed in terms of bboson can be schematically written as follows: H(0) =H0(0) +H+ 1(0) +H− 1(0) +H+ 2(0) +H− 2(0),(50) where H− 1,2=(H+ 1,2)†and 0, 1, or 2 subscripts indicate the number of created (+) or annihilated (−) excitations. To perform the CUTs, we choose the so-called quasiparticle conserving generator first proposed by Mielke [34] in the context of finite matrices and generalized to many-body problems by Knetter and Uhrig [35], which reads η(l)=H+ 1(l)−H− 1(l)+H+ 2(l)−H− 2(l).(51) In the symmetric phase (H± 1=0) this choice coincides with the generator proposed by Stein [36]. The flow equations are 064332-5 DUSUEL, VIDAL, ARIAS, DUKELSKY, AND GARC´ IA-RAMOS PHYSICAL REVIEW C 72, 064332 (2005) then simple quadratic functions of the Hamiltonians: ∂lH0(l)=2([H+ 1(l),H− 1(l)] +[H+ 2(l),H− 2(l)]),(52) ∂lH+ 1(l)=[H+ 1(l),H 0(l)] +2[H+ 2(l),H− 1(l)],(53) ∂lH+ 2(l)=[H+ 2(l),H 0(l)].(54) In the limit l=∞, the Hamiltonian conserves the number of bboson so that H± 1,2(∞)=0 and H(∞)=H0(∞). Following the method developed for the LMG model in Ref. [14,24], we convert these equations, which deal with operators, into equations involving coupling constants. This is achieved by expanding Hamiltonians H0and H± 1,2in powers of 1/N (see Sec. IV B). B. Flow equations for the Hamiltonian In the symmetric phase (H± 1=0), we have three elementary operators : nb:,P† b,P bfrom which H0and H± 2are built. More precisely, the 1/N expansion of these Hamiltonians can be written as follows: H0(l)= α,β,δ∈IN h(δ) 0,α,β (l)P† b β:nα b:Pbβ Nα+2β+δ−1,(55) H+ 2(l)= α,β,δ∈IN h(δ) 2,α,β (l)P† bP† b β:nα b:Pbβ Nα+2β+δ.(56) Note that for L=0, one has P† bPb=:n2 b:, so that h(δ) k,α,β = h(δ) k,α+2β,0. One then readily recovers expressions given in Ref. [24] for the case of a scalar bboson. Using this expansion and Eqs. (52)–(54), we can easily derive the flow equations for the couplings h(δ) k,α,β(l),which are given in Appendix A up to order (1/N)2. These flow equations can be solved exactly and, at order (1/N)2, one finally obtains the following: H(∞)=h0,0,0(∞)+h0,1,0(∞)nb+h0,2,0(∞):n2 b: +h0,0,1(∞)P† bPb+h0,3,0(∞):n3 b: +h0,1,1(∞)P† bnbPb+O(1/N3) (57) with h0,α,β (l)= δ∈IN h(δ) 0,α,β (l) Nα+2β+δ−1,(58) and h(0) 0,0,0(∞)=1−x 4,(59) h(1) 0,0,0(∞)=2L+1 21−3x 2+(x)1/2,(60) h(2) 0,0,0(∞)=(2L+1)(1 −x)x ×−(2L+5) +(6L+13)x 16(x)−L+2 4(x)1/2,(61) h(3) 0,0,0(∞)=−(2L+1)(1 −x)x2(2L+1) −(8L2+6L−33)x+(32L2−2L−149)x2−(24L2−38L−179)x3 128(x)5/2 −2L+5+(2L2−3L−17)x−(2L2−5L−20)x2 16(x)2,(62) h(0) 0,1,0(∞)=(x)1/2,(63) h(1) 0,1,0(∞)=(1 −x)x−1+(2L+5)x 4(x)−L+2 2(x)1/2,(64) h(2) 0,1,0(∞)=−(1 −x)x2L+1−(2L2+3L−6)x+(12L2+15L−23)x2−(10L2+5L−32)x3 16(x)5/2 −1+(2L2+5L−1)x−(2L2+3L−4)x2 4(x)2,(65) h(0) 0,2,0(∞)=(1 −x)x2 4(x),(66) h(1) 0,2,0(∞)=−(1 −x)x2 ×1−3x+(12L+29)x2−3(4L+9)x3 32(x)5/2 −xL+1−Lx 4(x)2,(67) h(0) 0,0,1(∞)=x(1 −x)(3x−1) 8(x),(68) h(1) 0,0,1(∞)=−(1 −x)x2(2L+5)1−3x+11x2−9x3 128(x)5/2 −x1+(L−3)x−(L−4)x2 8(x)2,(69) h(0) 0,3,0(∞)=−(1 −x)2x4 8(x)5/2,(70) h(0) 0,1,1(∞)=−(1 −x)2x2(1 −2x+9x2) 64(x)5/2,(71) where we have set (x)=x(2x−1). 064332-6 CONTINUOUS UNITARY TRANSFORMATIONS IN TWO- . . . PHYSICAL REVIEW C 72, 064332 (2005) The ground-state energy per particle is thus given by the following: e0(N)=h(0) 0,0,0(∞)+1 Nh(1) 0,0,0(∞)+1 N2h(2) 0,0,0(∞) +1 N3h(3) 0,0,0(∞)+O(1/N4),(72) whereas the gap reads (N)=h(0) 0,1,0(∞)+1 Nh(1) 0,1,0(∞)+1 N2h(2) 0,1,0(∞) +O(1/N3).(73) Of course, these expressions coincide for L=0 with those given in Refs. [14,24]. For L=2, one recovers the results given in Ref. [37]. The mean-field result (33) is also recovered in the thermodynamical limit. It is important to note that the Hamiltonian H(∞)=H0(∞) is not diagonal in the eigenbasis of nb(except for L=0) even though it always commutes with nb. Consequently, for each number of excitations, Hmust be diagonalized. As can be observed in Eqs. (59)–(71), some divergences appears, for x=xc, in the subleading corrections. We will see in Sec. VI that the structure of this singular 1/N expansion at the critical point allows us to extract nontrivial scaling exponents whose determination is one of the main motivation of this work. C. Flow equations for b† µ We now proceed to derive the flow equation for the operator b† µ(l) from which any other observable can be obtained. Analogously to the treatment of the Hamiltonian flow equations, the first step consists in transforming the flow equation [Eq. (48)] for (l)=b† µ(l) into a set of flow equations for couplings. Therefore, we expand b† µ(l)inpower of 1/N. Generically, one may expect to generate any terms b† µ αP† b β:nγ b:Pbηbµν. Here, we shall restrict our discussion to order 1/N for which one only has eight operators b† µ(l)=A+(l)b† µ+A−(l)˜ bµ+B+(l)b† µnb+B−(l)nb˜ bµ +C+(l)P† b˜ bµ+C−(l)b† µPb+D+(l)b† µP† b +D−(l)Pb˜ bµ+O(1/N2),(74) where we have introduced ˜ bµ=(−1)µb−µand where, as previously, each function has a canonical 1/N expansion given by the number of bosonic operators it is associated with, namely A±(l)=A(0)±(l)+A(1)±(l) N+O(1/N2),(75) B±(l)=B(0)±(l) N+O(1/N2),(76) C±(l)=C(0)±(l) N+O(1/N2),(77) D±(l)=D(0)±(l) N+O(1/N2).(78) The initial condition is of course given by: b† µ(l=0) =b† µ so that one only has a nonvanishing initial coupling which is A(0) +(0) =1. The flow equations are then obtained for these couplings order by order using Eq. (48). The full set of equations is given in Appendix B. As for the couplings defining the running Hamiltonians, these equations can be solved exactly and lead to A(0) s(∞)=1 2(x)1/4,(79) A(1) s(∞)=−(1 −x) 16x2L+3 (x)7/4−2(L+2) (x)5/4,(80) A(0) d(∞)=(x)1/4 2,(81) A(1) d(∞)=(1 −x) 16x2L+3 (x)5/4−2(L+2) (x)3/4,(82) B(0) s(∞)=− 1−x 8x(x)7/4,(83) B(0) d(∞)=1−x 8x(x)5/4,(84) C(0) s(∞)=− 1−x 16x(x)7/4,(85) C(0) d(∞)=1−x 16x(x)5/4,(86) D(0) s(∞)=(1 −x)(3x−1) 32x2(x)7/4,(87) D(0) d(∞)=1−x2 32x2(x)5/4,(88) where we have set (x)=(2x−1)/x, Fs=1 2(F++F−) and Fd=1 2(F+−F−), for each function F=A, B, C, D. The above expansion of b† µ(∞) allows us to compute ||for any operator that can be expressed in terms of b† µand for any eigenstates |and |of the Hamiltonian H(∞). In the following, we consider two different examples to show the power of this approach. D. Expectation value of the occupation number nL Let us first consider the case where =nLand where |=|is the ground state of H. This quantity normalized by the number of bosons can be computed straightforwardly by means of the Hellmann-Feynman theorem which states nL N=∂ ∂y [(1 +y)h0,0,0],(89) where y=x/(1 −x). Because we have the expansion of h0,0,0 up to order (1/N)3, one can easily get nLat this order. Here, instead, we compute it in terms of the flow equation for the operator b† µobtained in the preceding section, using the fact that: nL(∞)=nb(∞)= µ b† µ(∞)bµ(∞),(90) 064332-7 DUSUEL, VIDAL, ARIAS, DUKELSKY, AND GARC´ IA-RAMOS PHYSICAL REVIEW C 72, 064332 (2005) where b† µ(∞) is given by Eq. (74) with final values Eqs. (79)– (88). The ground state of the Hamiltonian (57) being defined as the zero bboson state |0, one has the following: 0|nL(∞)|0 N =1 N µ0|b† µ(∞)bµ(∞)|0, =1 N µ A2 −(∞)0|bµb† µ|0+O(1/N3), =1 N(2L+1)A(0) −(∞)2+2 NA(0) −(∞)A(1) −(∞) +O(1/N3), =(2L+1) N3x−1 4(x)1/2−1 2+(2L+1) N2 x(1 −x)2 16 ×−(2L+3)x (x)2+2(L+2) (x)3/2+O(1/N3),(91) where, as previously, (x)=x(2x−1). It can be easily verified that this expression coincides with Eq. (89). E. Transition probability between the ground state and the first excited state As explained above, the real power of the CUTs method is that it allows to easily compute off-diagonal matrix elements of any operator between any eigenstates of the Hamiltonian provided one knows the expression of the associated running operator. As an example, we focus here on the transition probability T=|1|TL0|0|2between the ground state |0and the first excited state |1. The operator TLµwas defined in Eq. (24). It is important to note that here the ground state has a zero angular momentum, whereas the first excited state has an angular momentum L. To determine the matrix element of TL0of interest, we shall proceed as for the occupation number and consider its 1/N expansion in terms of the bboson: TL0=s†L0+L† 0s =N1/2b† 0+b0−1 2N(b† 0nb+nbb0)+O(1/N2) .(92) To be consistent, given that we only have the expression of b† 0(∞)+b0(∞) at order 1/N, we need to consider b† 0(∞)nb(∞)+nb(∞)b0(∞) at order (1/N)0. Using the expression (74) of the operator b† 0, one easily gets 1|b† 0(∞)+b0(∞)|0=A+(∞)+A−(∞),(93) and 1|b† 0(∞)nb(∞)+nb(∞)b0(∞)|0 =(2L+3)A+(∞)A−(∞)2+A−(∞) ×[(2L+2)A−(∞)2+A+(∞)2].(94) Truncating these expressions at order 1/N and (1/N)0respec- tively and using Eqs. (75), (79) and (80), one finally obtains the following: T N=x (x)1/2+x2 N ×−(2L+1) −4(2L+1)x+(10L+7)x2 4(x)2 +−L+(3L+2)x 2(x)3/2+O(1/N2).(95) As for the expansion of the spectrum, some singularities appears, and we shall see that they also provide the scaling exponents at the critical point. V. THE BROKEN PHASE (0 <x<1/2) As shown by the mean-field analysis, for x<1/2, the order parameter nL/N has a nonvanishing value. This implies that we have to consider a new vacuum for the Holstein-Primakoff expansion. Therefore, we shift the bosonic modes by a term proportional to √N. We thus define the cbosons by the following: b† µ=√Nλ∗ µ+c† µ(96) where the λµ’s are complex numbers that form a (2L+1)- dimensional vector. Of course, the symmetric phase results are recovered when setting λµ=0. Then, using Eqs. (36)–(38) and assuming that nc/N 1, we expand the Hamiltonian that now contains some a term proportional to √Nthat reads as follows: (c†·˜ λ+λ†·˜ c)x+1−x 4[−2(1 −nλ)+P† λ+Pλ] +1−x 2[(c†·λ†)(Pλ−1+nλ)+(λ·c)(P† λ−1+nλ)],(97) where λ† µ=λ∗ µand ˜ λµ=(−1)µλ−µ. There are several choices of the λµs that allow one to get rid of these terms. Here, we have chosen to set λµ=λ0δµ,0with λ2 0=(1 −2x)/[2(1 −x)]. Note that in the thermodynamical limit, we recover the mean-field value [Eq. (35)] nL N= µ|λµ|2=1−2x 2(1 −x).(98) Further, we emphasize that this choice of the λµsisthesame as the one proposed in the mean-field analysis where we have broken the spherical symmetry by populating macroscopically µ=0 boson level only. With this choice, the Hamiltonian reads as follows: H=−Nx 3x−2 4(1 −x)+N0(1 −3x)(1 −2x) 8(1 −x)+x 2nc +5 4(1 −2x)c† 0c0−x 4(P† c+Pc)+3 8(1 −2x)c† 0 2+c2 0 +O(1/√N).(99) Contrary to the symmetric phase, we do restrict our discussion to this order because, as shown later, the existence of gapless modes at this level does not allow us to go beyond this order with the CUTs. 064332-8 CONTINUOUS UNITARY TRANSFORMATIONS IN TWO- . . . PHYSICAL REVIEW C 72, 064332 (2005) A. The spectrum The Hamiltonian (99) can be easily diagonalized via a Bogoliubov transform. Therefore, we introduce the dbosons defined by the following: c† µ=cosh(µ/2)d† µ+sinh(µ/2) ˜ dµ,(100) ˜ cµ=sinh(µ/2)d† µ+cosh(µ/2) ˜ dµ.(101) The angles µare chosen so that Hwritten in terms of the d’s is diagonal. From Eq. (99), it is clear that modes with µ= 0 and µ=0 plays a different role and actually decouple. As can be easily seen, eliminating off-diagonal terms for µ= 0 implies to set µ=0=∞and gives 2Lgapless modes. Because such a transform is singular, one has to use another route. The contribution of terms with µ= 0 in the Hamiltonian reads as follows: Hµ+H−µ=x 2[c† µcµ+c† −µc−µ−(−1)µ(c† µc† −µ+c−µcµ)]. (102) Introducing the position and momentum operator Xµ=c† µ+cµ √2,P µ=ic† µ−cµ √2,(103) one has Hµ+H−µ=−x 2+x 4[(Pµ+(−1)µP−µ)2 +(Xµ−(−1)µX−µ)2].(104) Because [Pµ+(−1)µP−µ,X µ−(−1)µX−µ]=0,H µis written in a diagonal form and its spectrum is indeed found to be gapless and continuous. The correction to the ground-state energy coming from this contribution is thus −Lx/2. Let us now consider the µ=0 part of the Hamiltonian that reads: Hµ=0=5 4−2xc† 0c0+3 8−xc† 0 2+c2 0.(105) For this contribution, the Bogoliubov transform can be simply achieved and one obtains the following: Hµ=0=1 2√1−2x−5 4+2x+√1−2xd † 0d0.(106) The correction to the ground-state energy coming from this contribution is thus given by the d0boson state. As a result, e0 at this order reads as follows: e0(N)=x 4 2−3x 1−x+1 N−2−2(L−2)x+(2L−1)x2 4(1 −x) +1 2√1−2x+O(1/N2).(107) As in the symmetric phase, the leading corrections coincide with the mean-field result (34) and Lonly appears in the subleading terms. Concerning the gap, the above analysis indicates the existence of 2Lgapless modes and a gapped one with excitation energy (N)=√1−2x+O(1/N).(108) As previously, one can simply obtain nL/N by replacing h0,0,0by e0(N) in Eq. (89), and the result is as follows: nL N=1−2x 2(1 −x)−1 2N ×L+x √1−2x+x x−1+O(1/N2).(109) At this stage, one can understand the difficulty to go beyond this order in the presence of gapless modes. Indeed, computing the next-order corrections would imply to keep on making the 1/N expansion around the (broken) vacuum, but such a procedure does not take into account the degeneracy because of the gapless modes. Note that for L=0 where no gapless modes emerge, we have been able to compute these corrections using CUTs [24]. To conclude this subsection, we wish to underline that in the two-level BCS model where gapless modes also exist, Richardson has obtained the finite-size corrections in the broken phase beyond the Bogoliubov order using the 1/N expansion of the exact solution [38], whereas we computed them more recently using CUTs in the symmetric phase [39]. B. Transition probability between the ground state and the first excited state As in the symmetric phase, we now compute the transition T=|1|TL0|0|2, where TL0=s†L0+L† 0s. However, the important difference is that, in the broken phase, one has 2L gapless modes that renders the definition of the first excited states more tricky. In the thermodynamical limit, the ground state thus becomes infinitely degenerate and one actually has to simply consider the expectation value of TL0over the ground state. To avoid any confusion, we call this quantity Tinstead of T. Using the expansion (93) and the shift (96) with the choice of the λµgiven previously, one obtains the following: T=|0|TL0|0|2, =N24λ2 01−λ2 0+O(N), =N21−2x (1 −x)2+O(N).(110) First, it is important to note that Tis proportional to N2in this phase, whereas Tscales as Nin the symmetric phase. Second, in the broken phase, Tvanishes at the critical point, whereas Tdiverges when approaching from the symmetric phase. This result clearly suggests an anomalous scaling behavior at the critical point that we shall now investigate in details. VI. THE CRITICAL POINT In this section, we analyze the behavior of the 1/N expansion of the quantities considered in this study: the ground-state energy, the gap, the expectation value of nLin the ground state, and the transition rate Tbetween the ground state and the first excited state. The common point of all these expansions is that they become singular at the critical point. Following the arguments presented in a recent series 064332-9 DUSUEL, VIDAL, ARIAS, DUKELSKY, AND GARC´ IA-RAMOS PHYSICAL REVIEW C 72, 064332 (2005) APPENDIX C: SOLVING THE FLOW EQUATIONS The flow equations given in the above appendices have to be solved order by order in 1/N.Atorder(1/N)−1, nothing has to be done, so let us turn to order (1/N)0. The equations for h(0) 0,1,0(l) and h(0) 2,0,0(l) are easily solved by noticing that h(0) 0,1,0(l)2−4h(0) 2,0,0(l)2is a constant of the flow. One gets the hyperbolic solutions h(0) 0,1,0(l)=∞ tanh[2∞(l+l0)],(C1) h(0) 2,0,0(l)=−sgn(ε)∞ 2sinh[2∞(l+l0)].(C2) We have denoted ∞the gap of the system at the thermodynamical limit, that is h(0) 0,1,0(∞) [see Eq. (63)]. The quantity l0is such that the initial conditions are fulfilled, namely h(0) 0,1,0(0) =∞/tanh[2∞l0], and we also introduced ε=−h(0) 2,0,0(0)/[2h(0) 0,1,0(0)]. As already explained in Ref. [24], the best way to solve the flow equations is in fact to introduce a new “time scale” that is more adapted to the problem and defined by the following: t=sgn(ε)exp[2∞(l+l0)],(C3) with initial conditions now given at t0=sgn(ε)exp(2∞l0).(C4) After some algebra t0can also be shown to be equal to t0=1 ε(1 +1−ε2).(C5) Eqs. (C1) and (C2) now read as follows: h(0) 0,1,0(t)=∞ t2+1 t2−1,(C6) h(0) 2,0,0(t)=−∞ t t2−1,(C7) which are rational expressions in t. The renormalized values at l→∞are now found by taking the limit t→t∞=sgn(ε)∞. Let us remark that the off-diagonal coupling h(0) 2,0,0(t) goes to zero and behaves like t−1for t→t∞. This will be true for all off-diagonal couplings creating two excitations because the energy cost of such excitations, in the thermodynamic limit and for large t, is nothing but 2∞, so that the couplings must vanish as exp(−2∞l). The last flow equation (for the spectrum) at order (1/N)0 is solved by noticing that 2h(1) 0,0,0(t)−(2L+1)h(0) 0,1,0(t)is a constant of the flow, equal to its initial value, namely −(2L+1)(3x−1)/2. The last task at order (1/N)0is to obtain the solution for the observable. For this, one simply has to insert Eq. (C7) into Eqs. (B1) and (B2) and then replace ∂lwith 2∞t∂tand solve the resulting equation. This yields A(0) s(t)=1 2(t−1)(t0+1) (t+1)(t0−1),(C8) A(0) d(t)=1 2(t+1)(t0−1) (t−1)(t0+1).(C9) The next orders are solved in the same fashion: one inserts the expressions known from the previous orders, replaces ∂l with 2∞t∂t, and solves the equations (which is most simply achieved thanks to a computer algebra program). We refer the interested reader to the details given in Ref. [24], where some more technical details are given. Let us emphasize that the usefulness of the tvariable comes from the fact that there exists only one basic energy scale in the problem, namely the gap. All energy scales are integer multiples of this gap. As previously mentioned, off-diagonal couplings associated to an energy scale 2∞decay as t−1, and in the general case, an off-diagonal coupling whose energy scale is n∞decay as t−n/2. Such a time variable would be useless in a problem where many different energy scales exist. [1] E. T. Jaynes and F. W. Cummings, Proc. IEEE 51, 89 (1963). [2] F. Iachello and R. D. Levine, Algebraic Theory of Molecules (Oxford University Press, Oxford, 1995). [3]H.Suhl,B.T.Matthias,andL.Walker,Phys.Rev.Lett.3, 552 (1959). [4] J. H¨ ogaasen-Feldman, Nucl. Phys. 28, 258 (1961). [5] H. J. Lipkin, N. Meshkov, and A. J. Glick, Nucl. Phys. 62, 188 (1965). [6] N. Meshkov, A. J. Glick, and H. J. Lipkin, Nucl. Phys. 62, 199 (1965). [7] A. J. Glick, H. J. Lipkin, and N. Meshkov, Nucl. Phys. 62, 211 (1965). [8] F. Iachello and A. Arima, The Interacting Boson Model (Cambridge University Press, Cambridge, 1987). [9] J. Vidal, G. Palacios, and R. Mosseri, Phys. Rev. A 69, 022107 (2004). [10] J. Vidal, R. Mosseri, and J. Dukelsky, Phys. Rev. A 69, 054101 (2004). [11] J. Vidal, G. Palacios, and C. Aslangul, Phys. Rev. A 70, 062304 (2004). [12] J. I. Latorre, R. Or´ us, E. Rico, and J. Vidal, Phys. Rev. A 71, 064101 (2005). [13] F. Leyvraz and W. D. Heiss, Phys. Rev. Lett. 95, 050402 (2005). [14] S. Dusuel and J. Vidal, Phys. Rev. Lett. 93, 237204 (2004). [15] R. G. Unanyan, C. Ionescu, and M. Fleischhauer, Phys. Rev. A 72, 022326 (2005). [16] H.-Q. Zhou, J. Links, R. H. McKenzie, and X.-W. Guan, J. Phys. A36, L113 (2003). [17] F. Pan, X. Zhang, and J. P. Draayer, J. Phys. A 35, 7173 (2002). [18] R. W. Richardson, J. Math. Phys. 9, 1327 (1968). [19] F. Iachello and N. V. Zamfir, Phys. Rev. Lett. 92, 212501 (2004). [20] F. Pan and J. P. Draayer, Nucl. Phys. A636, 156 (1998). [21] D. J. Rowe, Nucl. Phys. A745, 47 (2004). [22] J. E. Garc´ ıa-Ramos,J.Dukelsky,andJ.M.Arias,Phys.Rev.C 72, 037301 (2005). 064332-16 CONTINUOUS UNITARY TRANSFORMATIONS IN TWO- . . . PHYSICAL REVIEW C 72, 064332 (2005) [23] J. M. Arias, C. E. Alonso, A. Vitturi, J. E. Garc´ ıa-Ramos, J. Dukelsky, and A. Frank, Phys. Rev. C 68, 041302(R) (2003). [24] S. Dusuel and J. Vidal, Phys. Rev. B 71, 224420 (2005). [25] T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940). [26] A. Klein and E. R. Marshalek, Rev. Mod. Phys. 63, 375 (1991). [27] J. Dukelsky, G. G. Dussel, R. P. J. Perazzo, S. L. Reich, and H. M. Sofia, Nucl. Phys. A425, 93 (1984). [28] D. J. Rowe, Phys. Rev. Lett. 93, 122502 (2004). [29] F. Wegner, Ann. Phys. (Leipzig) 3, 77 (1994). [30] S. D. Głazek and K. G. Wilson, Phys. Rev. D 48, 5863 (1993). [31] S. D. Głazek and K. G. Wilson, Phys. Rev. D 49, 4214 (1994). [32] F. Wegner, Phys. Rep. 348, 77 (2001). [33] B. H. Bartlett, Master Thesis (2003), nucl-th/0305052. [34] A. Mielke, Eur. Phys. J. B 5, 605 (1998). [35] C. Knetter and G. S. Uhrig, Eur. Phys. J. B 13, 209 (2000). [36] J. Stein, J. Phys. G 26, 377 (2000). [37] S. Dusuel, J. Vidal, J. M. Arias, J. Dukelsky, and J. E. Garc´ ıa- Ramos, Phys. Rev. C 72, 011301(R) (2005). [38] R. W. Richardson, J. Math. Phys. 18, 1802 (1977). [39] S. Dusuel and J. Vidal, Phys. Rev. A 71, 060304(R) (2005). [40] D. J. Rowe, P. S. Turner, and G. Rosensteel, Phys. Rev. Lett. 93, 232502 (2004). [41] F. Pan, Y. Zhang, and J. P. Draayer, J. Phys. G 31, 1039 (2005). [42] R. Botet, R. Jullien, and P. Pfeuty, Phys. Rev. Lett. 49, 478 (1982). [43] R. Botet and R. Jullien, Phys. Rev. B 28, 3955 (1983). 064332-17