{'fraction_non_alphanumeric': 0.03897088743387711, 'fraction_numerical': 0.01783005955683794, 'mean_word_length': 4.420250626566416, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, '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}}
Paramagnetic State of Degenerate Double Exchange Model: a Non-Local CPA Approach

The off-diagonal disorder caused by random spin orientations in the paramagnetic (PM) state of the double exchange (DE) model is described by using the coherent-potential-approximation (CPA), which is combined with the variational mean-field approach for the Curie temperature (T_C). Our CPA approach is essentially non-local and based on the perturbation theory expansion for the T-matrix with respect to the fluctuations of hoppings from the ''mean values'' specified by matrix elements of the self energy, so that in the first order it becomes identical to the DE theory by de Gennes. The second-order effects, considered in the present work, can be viewed as an extension of this theory. They are not negligible, and lead to a substantial reduction of T_C in the one-orbital case. Even more dramatic changes are expected in the case of orbital degeneracy, when each site of the cubic lattice is represented by two e_g orbitals, which also specify the form of interatomic transfer integrals. Particularly, the existence of two Van Hove singularities in the spectrum of degenerate DE model may lead to the branching of CPA solutions, when the Green function and the self energy become double-valued functions in certain region of the complex plane. Such a behavior can be interpreted as an intrinsic inhomogeneity of the PM state, which consists of two phases characterized by different electronic densities. The phase separation occurs below certain transition temperature, T_P, and naturally explains the appearance of several magnetic transition points, which are frequently seen in manganites.


I. INTRODUCTION
The nature of paramagnetic (PM) state in perovskite manganese oxides (the manganites) is one of the fundamental questions, the answer to which is directly related with understanding of the phenomenon of colossal magnetoresistance.
There is no doubts that any theoretical model for manganites should include (at least, as one of the main ingredients) the double exchange (DE) physics, which enforces the atomic Hund's rule and penalizes the hoppings of polarized e g electrons to the sites with the opposite direction of the localized t 2g spins. [1][2][3] If the spins are treated classically, the corresponding Hamiltonian is given, in the local coordinate frame specified by the directions e i =(cos φ i sin θ i , sin φ i sin θ i , cos θ i ) of the spin magnetic moments, by 4 where t ij are the bare transfer integrals between sites i and j (which can be either scalars or 2×2 matrices depending on how the orbital degeneracy of the e g states is treated in the model), and ξ ij =cos θ i 2 cos θ j 2 +sin θ i 2 sin θ j 2 e −i(φ i −φ j ) describes their modulations caused by the deviation from the ferromagnetic (FM) alignment in the pair i-j.
The form of Hamiltonian (1) implies that the spin magnetic moments are saturated (due to the strong Hund's rule coupling) and the spin disorder corresponding to the PM state is in fact an orientational spin disorder. Despite an apparent simplicity of the DE Hamiltonian (1), the description of this orientational spin disorder is not an easy task and so far there were only few theoretical model, which were based on rather severe (and presumably unsatisfactory) approximations.
The first one was proposed by de Gennes more than forty years ago. 3 In his theory, all ξ ij are replaced by an averaged value ξ, so that the spin disorder enters the model only as a renormalization (narrowing) of the e g -bandwidth. The effect is not particularly strong and the fully disordered PM state corresponds to ξ= 2 3 . A generalization of this theory to the case of quantum spins was given by Kubo and Ohata. 5 The same idea was exploited recently in a number of theories aiming to study the behavior of orbital degrees of freedom at elevated temperature, 6 but based on the same kind of simplifications for the spin disorder and its effects on the kinetic energy.
Another direction, which features more recent activity, is the single-site dynamical meanfield theory (DMFT, see Ref. 7 for a review) for the FM Kondo lattice model. 8,9 The model itself can be viewed as a prototype of the DE Hamiltonian (1) before projecting out the minority-spin states in the local coordinate frame. 4 If the localized t 2g spins are treated classically (that is typically the case), this method is similar to the disordered local moment approach proposed by Gyorffy et al., 10 and based on the coherent-potential-approximation (CPA) for the electronic structure of the disordered state. 11 Now it is almost generally accepted that both approaches are inadequate as they fail to explain not even all, but a certain number of observations in manganites of a principal character such as the absolute value and the doping dependence of the Curie temperature (T C ), 12 the insulating behavior above T C , 13 and a rich magnetic phase diagram along the temperature axis, which typically show a number of magnetic phase transitions 14 and the phase coexistence [15][16][17] in certain temperature interval. Therefore, it is clear that the theory must be revised. Then, there are two possible ways to precede. One is to modernize the model itself by including additional ingredients such as the Jahn-Teller distortion, 18 the Coulomb correlations, 6,19 and the disorder effects caused by the chemical substitution. 20 Another possibility is to stick to the basic concept of the DE physics and try to formulate a more advanced theory of the spin disorder described by the Hamiltonian (1), which would go beyond the simple scaling theory by de Gennes 3 as well as the single-site approximation inherent to DMFT. 8,9 An attention to the second direction was drawn recently by Varma. 21 The main challenge to the theoretical description of the orientational spin disorder in the DE systems comes from the fact that it enters the Hamiltonian (1) as an off-diagonal disorder of interatomic transfer integrals, which presents a serious and not well investigated problem. In the present work we try to investigate some possibilities along this line by employing a non-local CPA approach.
What do we expect?
1. It was realized very recently that many aspects of seemingly complicated lowtemperature behavior of the doped manganites (the rich magnetic phase diagram, optical properties, etc.) can be understood from the viewpoint of DE physics, if the latter is considered in the combination with the realistic electronic structure for the itinerant e g electrons and takes into account the strong dependence of this electronic structure on the magnetic ordering (see, e.g., Ref. 22 and references therein). If this scenario is correct and can be extended to the high-temperature regime (that is still a big question), there should be something peculiar in the electronic structure of the disordered PM state, which can be linked to the unique properties of manganites. In order to gain insight into this problem, let us start with the DE picture by de Gennes, and consider it in the combination with the correct form of the transfer integrals between two e g orbitals on the cubic lattice. 23 Then, the properties of the PM state should be directly related with details of the electronic structure of the FM state, which is shown in Fig. 1 and connected with the PM electronic structure by the scaling transformation. This electronic structure is indeed very peculiar because of two Van Hove singularities at the (π, π, 0) and (0, π, 0) points of the Brillouin zone, which are responsible for two kinks of density of states at ±1. It is also interesting to note that the first singularity appears near the Fermi surface when the hole concentration is close to 0.3, i.e. in the most interesting regime from the viewpoint of colossal magnetoresistance. 24 Such a behavior was discussed by Dzero, Gor'kov and Kresin. 25 They also argued that this singularity can contribute to the T 3/2 dependence of the specific heat in the FM state.
If so, what is the possible role of these singularities in the case of the spin disorder? Note that apart from the single-site approximation, all recent DMFT calculations 8,9 employed a model semi-circular density of states, and therefore could not address this problem.
2. There are many anticipations in the literature that there is some hidden degree of freedom which controls the properties of perovskite manganites. A typical example is the picture of orbital disorder proposed in Ref. 19 in connection with the anomalous behavior of the optical conductivity in the FM state of perovskite manganites. According to this picture, the large on-site Coulomb interaction gives rise to the orbital polarization at each site of the system. The local orbital polarizations remain even in the cubic FM phase, but without the long-range ordering. In the present work we will show that, in principle, by considering non-local effects in the framework of pure DE model, one may have an alternative scenario, when there is a certain degree of freedom which does control the properties of manganites. However, contrary to the orbital polarization, this parameter is essentially non-local and attached to the bond of the DE system rather than to the site.
3. There were many debated about the phase separation in perovskite manganites, 26,27 and according to some scenarios this effect plays an important role also at elevated temperatures, being actually the main trigger behind phenomenon of the colossal magnetoresistance. 20 The problem was intensively studied numerically, using the Monte Carlo techniques. 20,27 If this is indeed the case, what does it mean on the language of analytical solutions of the DE model (or its refinements)? Presumably, the only possibility to have two (and more) phases at the same time is to admit that the self energy (and the Green function) of the model can be a multi-valued function in certain region of the complex plane. Such a behavior of non-linear CPA equations was considered as one of the main troublemakers in the past, 28,29 but may have some physical explanation in the light of newly proposed ideas of phase separation.
The rest of the paper is organized as follows. In Sec. II we briefly review the main ideas of the variational mean-field approach. In Sec. III we describe general ideas of the non-local CPA to the problem of orientational spin disorder in the DE systems. In Sec. IV we consider the CPA solution for the PM phase of the one-orbital model and evaluate the Curie temperature. We will argue that two seemingly different approaches to the problem of spin disorder in the DE model, one of which was proposed by de Gennes 3 and the other one is based on the DMFT 8,9 , have a common basis and can be regarded as CPA-type approaches, but supplemented with different types of approximations. In Sec. V we consider more realistic two-orbital case for the e g electrons and argue that it is qualitatively different from the behavior of one-orbital model. Particularly, the CPA self-energy becomes the double-valued function in certain region of the complex plane, that can be related with an intrinsic inhomogeneity of the PM state of the DE model. In Sec. VI we summarize the main results of our work, discuss possible connections with the experimental behavior of perovskite manganites as well as possible extensions.

II. CALCULATION OF THERMAL AVERAGES
In order to proceed with the finite temperature description of the DE systems we adopt the variational mean-field approach. 3,5,30 Namely, we assume that the thermal (or orientational) average of any physical quantity is given in terms of the single spin orientation distribution function, which depends only on the angle between the local spin and an effective molecular field λ, while all correlations between different spins are neglected. 31 In the case of PM-FM transition, the effective field can be chosen as λ=(0, 0, λ).
The first complication comes from the fact that the DE Hamiltonian (1) is formulated in the local coordinate frame, in which the spin quantization axes at different sites are specified by the different direction {e i }. Therefore, we should clarify the meaning of orientational averaging in the local coordinate frame. 32 From our point of view, it is logical that in order to calculate the thermal averages associated with an arbitrary chosen site 0, the global coordinate frame should be specified by the direction e 0 =(cos φ 0 sin θ 0 , sin φ 0 sin θ 0 , cos θ 0 ), so that at each instant the spin moment at the site 0 is aligned along the z-direction. The averaging over all possible directions e 0 of the global coordinate frame in the molecular field λ can be performed as the second step.
Then, corresponding distribution function at the site 0 is given by Eq.
(2). The distribution functions at other sites can be constructed as follows. The transformation to the coordinate frame associated with the site 0 is given by the matrix: In the new coordinates, the i-th moment has the direction e , and the effective field becomes Obviously, while λ ′ and e ′ i depend on θ 0 and φ 0 , the new distribution function p i (e ′ i )∝exp(λ ′ · e ′ i ) -does not, and up to this stage the transformation to the new coordinate frame was only the change of the notations. In order to obtain the total distribution function P i (e ′ i ) for i =0, formulated in the local coordinates of the site 0 and taking into account the motion of e 0 in the molecular field λ, p i (e ′ i ) should be averaged over e 0 with the weight p 0 (e 0 ): The normalization constant ν is obtained from the condition: The form of Eq. (3) implies that the directions of magnetic moments are not correlated (in the spirit of the mean-field approach) and the averaging over e 0 can be performed independently for different sites of the system. For the analysis the PM state and the magnetic transition temperature, it is sufficient to consider the small-λ limit. Then, Eq. (3) becomes: The thermal average of the function F (e ′ i ), formulated in the coordinate frame of the site 0 and taking into account the motion of both e 0 and e i , is given by The spin entropy can be computed in terms of the molecular field λ as: 3 In the second order of λ this yields (both for i=0 and i =0): Then, the free energy of the DE model is given by 3,30 where E D (T, λ) is the electron free energy (or the double exchange energy): calculated in terms of the (orientationally averaged) integrated density of states n(z, λ).
] −1 is the Fermi-Dirac function (µ being the chemical potential). The best approximation for the molecular field λ is that which minimizes the free energy (6). Assuming that the transition to the FM state is continuous (of the second order), 33 T C can be found from equation: In practice, the derivative ∂ 2 E D (T, λ)/∂λ 2 near λ=0 can be calculated using the variational properties of n(z, λ) in CPA and the Lloyd formula. 34,35

III. NON-LOCAL CPA FOR THE DOUBLE EXCHANGE MODEL
In this section we discuss some general aspects of the non-local CPA to the problem of orientational spin disorder in the DE model specified by the Hamiltonian (1). We attempt to describe the disordered system in an average by introducing an effective energy-dependent Hamiltonian, where Σ ij is the non-local part of the self energy, which is restricted by the nearest neighbors; and Σ ii is the local (site-diagonal) part. The non-local formulation of CPA is essential because in the low-temperature limit H ij should be replaced by the conventional kinetic term in which Σ ij plays a role of the bare transfer integral t ij . Therefore, Σ ij cannot be omitted. On the other hand, Σ ii will be needed in order to formulate a closed system of CPA equations. We require the effective Hamiltonian (9) to preserve the cubic symmetry of the system and be translationally invariant. The first requirement has a different form depending on the degeneracy of the problem and the symmetry properties of basis orbitals, and will be considered separately for the one-orbital and degenerate DE models. In any case, using the symmetry properties, all matrix elements of the self energy {Σ ii , Σ ij } on the cubic lattice can be expressed through {Σ 00 , Σ 01 } for one of the dimers (for example, 0-1 in Fig. 2). Then, the Hamiltonian (9) can be Fourier transformed to the reciprocal space, , and the first equation for the orientationally averaged Green function can be written as where the integration goes over the first Brillouin zone with the volume Ω BZ . Unfortunately, the non-local form of the Hamiltonian (9) in the combination with the translational invariance do not necessarily guaranty the fulfillment of causality principle, that is the single-particle Green function should be analytic in the upper half of complex energy plane and satisfy a certain number of physical requirements. 28,29,[36][37][38] This is still a largely unresolved problem, despite numerous efforts over decades. We do not have a general solution to it either. What we try to do here is simply to investigate the behavior of this particular model and try to answer the question whether it is physical or not. Note also that neither dynamical cluster approximation 37 nor cellular DMFT method 38 , for which the causality can be rigorously proven, can be easily applied to the problem of off-diagonal disorder.
The calculations near Van Hove singularities in the case of degenerate DE model requires very accurate integration in the reciprocal space. In the present work we used the mesh consisting of 374660 nonequivalent q-points and corresponding to 258×258×258 divisions of the reciprocal lattice vectors.
In order to formulate the CPA equations we consider only site-diagonal and nearestneighbor elements of G ij . Again, using the symmetry properties they can be expressed through {G 00 , G 01 }. In addition, there is a simple relation between G 00 and G 01 for given Σ 00 and Σ 01 : which follows from the definition of the Green function (10) and the Hamiltonian (9). 39 Using this identity, some matrix elements of the Green function can be easily excluded from CPA equations. In order to obtain the closed system of CPA equations which connects {Σ 00 , Σ 01 } with {G 00 , G 01 } we construct the T -matrix: 34 and require the average of scattering due to the fluctuations ∆ H= H− H to vanish on every site and every bond of the system, i.e.: 11,34 T 00 (z) = T 01 (z) = 0.
The hat-symbols in Eq. (12) means that all the quantities are infinite matrices in the real space and the matrix multiplications imply also summation over the intermediate sites. Note that our approach is different from the so-called cluster-CPA 28,29,36,40 because the matrix operations in Eq. (12) are not confined within a finite cluster (the dimer, in our case). We believe that our approach is more logical and more consistent with the requirements of the cubic symmetry and the translational invariance of the system, because all dimers are equivalent and should equally contribute to the averaged T -matrix. This equivalence is artificially broken in the cluster-CPA approach, which takes into account the contributions of only those atoms which are confined within the cluster. However, our approach also causes some additional difficulties, because the non-local fluctuations ∆ H tend to couple an infinite number of sites in Eq. (12). Therefore, for the practical purposes we restrict ourselves by the perturbation theory expansion up to the second order of ∆ H: As we will show, the first term in this expansion corresponds to the approximation considered by de Gennes, 3 and the next term is the first correction to this approximation. In order to evaluate the matrix elements T 00 (z) and T 01 (z) in the approximation given by Eq. (14) it is necessary to consider the interactions confined within the twelve-atom cluster shown in Fig. 2 (obviously, an additional term in the perturbation theory expansion for the T -matrix would require a bigger cluster). All such contributions are listed in Table I.

IV. ONE-ORBITAL DOUBLE EXCHANGE MODEL
In the one-orbital case, the effective DE Hamiltonian takes the following form, in the reciprocal space: where Σ 00 (z) and Σ 01 (z) are C-numbers, c γ =cos q γ , and all energies throughout in this section are in units of the effective transfer integral t 0 , which is related with the e g -bandwidth W in the FM state as t 0 =W/12. Elements of the Green function, G 00 (z) and G 01 (z), are obtained from Eq. (10). In subsequent derivations we will retain both G 00 (z) and G 01 (z), but only for the sake of convenience of the notations, because formally G 01 (z) can be expressed through G 00 (z) using identity (11). The self-consistent CPA equations are obtained from the conditions T 00 (z)=T 01 (z)=0. In the local coordinate frame associated with the site 0, e ′ 0 =(0, 0, 1) and all contributions to T 00 (z) and T 01 (z) shown in Table I should be averaged over the directions of magnetic moments of remaining sites of the cluster with probability functions given by Eq. (4). This is a tedious, but rather straightforward procedure. In order to calculate the thermal averages, it is useful to remember identities listed in Ref. 41. We drop here all details and, just for reader's convenience, list in Appendix A the averaged values for all contributions shown in Table I without the derivation, so that every step can be easily checked. We also introduce short notations for the self energies: Σ 00 ≡σ 0 and Σ 01 ≡2/3+σ 1 , and for the Green function: G 00 ≡g 0 and G 01 ≡g 1 . Then, the CPA equations T 00 (z)=0 and T 01 (z)=0 can be written in the form (for α=0 and 1, respectively): These equations should be solved self-consistently in combination with the definition (10) for the Green function. Different elements of G ij (z) and Σ ij (z) obtained in such a manner for the PM state (λ=0) are shown in Fig. 3. We note the following: 1. both G ij (z) and Σ ij (z) are analytic in the upper half of the complex plane; 2. ImG 00 (z)≤0 in the upper half plane; 3. the (numerically obtained) integrated density of states lies in the interval 0≤n(µ)≤1 and takes all values within this interval as the function of chemical potential µ, meaning that our system is well defined for all physical values of the electronic density.
We believe that for our purposes, the fulfillment of these three causality principles is quite sufficient and some additional requirements do not necessary apply here. 42 Note also that in the one-orbital case there is only one CPA solution in the complex energy plane.
In first order expansion for the T -matrix with respect to ∆H (thereafter all quantities corresponding to such approximation will be denoted by tilde), we naturally reproduce parameters of the DE model by Gennes: 3Φ 0 =Φ 1 =Ψ 0 =0 andΨ 1 =2/45. The corresponding elements of the Green function are also shown in Fig. 3. Although the second-order approach gives rise to the large matrix elements of the self energy, the elements of the Green function obtained in the first and second order with respect to ∆H are surprisingly close (apart from a broadening in the second-order approach, caused by the imaginary part of the self energy), meaning that there is a good deal of cancellations of different contributions to G 00 and G 01 .
However, it is not true for the Curie temperatures. While in the first order T C is solely determined by G 00 , in the second-order it explicitly depends on both G 00 and matrix elements of the self energy, which make significant difference from the canonical behavior. 3 Indeed, T C can be obtained from Eq. (8). In order to evaluate the DE energy, we start with the PM solution (λ=0) and include all contributions of the first order of λ 2 as a perturbation. Employing variational properties of the averaged integrated density of states, 34,35 n(z, λ) can be found using the Lloyd formula: where Ψ 0 , Ψ 1 , g 0 , and g 1 correspond to the PM state. 43 Moreover, g 1 (z) can be expressed through g 0 (z) using identity (11). Then, the DE energy takes the form: Taking into account the explicit expression for the entropy term, Eq. (5), we arrive at the following equation: In the first order with respect to ∆H (corresponding to the choiceσ 0 =σ 1 =Ψ 0 =0 and Ψ 1 =2/45) and after replacing f T (z−µ) by Θ(z−µ),D(T ) can be expressed through the DE energy of the PM state:D=− 1 15Ẽ D (PM). Taking into account thatẼ D (PM)= 2 3 E D (FM), we arrive at the well known expression obtained by de Gennes: The results of calculations are shown in Fig. 4. 44 A more accurate treatment of ∆H in the expression for the T -matrix significantly reduces T C (up to 20% in the second-order approach). The values T C ≃0.20t 0 and 0.17t 0 obtained correspondingly at n=0.5 and n=0.25 are substantially reduced in comparison with the (local) DMFT approach. 45 We also note a good agreement with the results by Alonso et al., 30 who used similar variational mean-filed approach supplemented with the moments-method for computing the averaged density of states in the one-orbital DE model. Using the value W ≈4eV obtained in bands structure calculations for the FM state, 22 T C can be roughly estimated as T C ≤800K. The upper bound, corresponding to n=0.5, exceeds the experimental data by factor two. T C can be further reduced by taking into consideration the antiferromagnetic (AFM) superexchange (SE) interactions between the localized spins 3,30,46,47 and spatial spin correlations. For example, according to recent Monte Carlo simulations, the latter can reduce T C up to 0.12t 0 ≈460K at n=0.5. 48

V. DEGENERATE DOUBLE EXCHANGE MODEL FOR THE e g ELECTRONS
In this section we consider a more realistic example of the double exchange model involving two e g orbitals, which have the following order: |1 ≡|x 2 −y 2 and |2 ≡|3z 2 −r 2 . Then all quantities, such as the transfer integrals t ij , the averaged Green function G ij (z), and the self energy Σ ij (z) become 2×2 matrices in the basis of these orbitals (thereafter, the bold symbols will be reserved for such matrix notations). The DE model is formulated in the same way as in the one-orbital case. Namely, modulations of the transfer integrals are described by Eq. (1) with the complex multipliers ξ ij . What important is the peculiar form of the t ij matrices on the cubic lattice, given in terms of Slater-Koster integrals of the ddσ type. 23 For example, for the 0-1 and 0-2 bonds parallel to the z-axis (see Fig. 2) these matrices have the form: In the other words, the hoppings along the z-direction are allowed only between the 3z 2 −r 2 orbitals. Throughout in this section, the absolute values of the parameter (ddσ)=W/6≈0.7eV will be used as the energy unit. The remaining matrix elements in the xy-plane can be obtained by 90 • rotations of the 0-1 bond around the x-and y-axes. Corresponding transformations of the e g orbitals are given by: and respectively. Then, it is easy to obtain the well known expressions: 23 for the hoppings parallel to the x-axis, and for the hoppings parallel to the y-axis.
We remind these symmetry properties as an introduction to the analysis of the self energy in the case of the spin disorder, which obeys very similar symmetry constraints. Namely, since the cubic symmetry is not destroyed by the disorder, the x 2 −y 2 and 3z 2 −r 2 orbitals belong to the same representation of the point symmetry group. Therefore, the site-diagonal part of the self energy (as well as of the Green function) will be both diagonal and degenerate with respect to the orbital indices, i.e.: where Σ 11 00 =Σ 22 00 . On the other hand, the matrix elements associated with the bond 0-1 should transform to themselves according to the tetragonal (C 4v ) symmetry. Since the x 2 −y 2 and 3z 2 −r 2 orbitals belong to different representations of the C 4v group (a 1g and b 2g , respectively), the corresponding non-local part of the self energy has the form: Note that Σ 11 01 (z) is not necessarily zero. The identity t 11 01 =0, which holds for the transfer integrals, reflects the hidden symmetry of the ordered FM state. However, there is no reason to expect that the same identity will be preserved in the case of the spin disorder. Moreover, as we will show below, the condition Σ 11 01 (z) =0 is indispensable in order to formulate the closed system of CPA equations.
Thus, in the case of the orbital degeneracy, there are three independent matrix elements of the self energy: Σ 11 00 (z)=Σ 22 00 (z), Σ 11 01 (z) and Σ 22 01 (z). In comparison with the one-orbital case we gain an additional non-local parameter, which may control the properties of disordered DE systems, even in the case of the strictly imposed cubic symmetry.
The matrix elements of the self energy in the xy-plane can be obtained using the transformations (20) and (21), which yield: .
Then, the effective Hamiltonian takes the following form, in the reciprocal space: The orientationally averaged Green function can be then obtained from Eq. (10). Similar to the self energy, all matrix elements of the Green function can be expressed through G  (22) and (23). In addition, they satisfy the matrix equation (11), from which one can easily exclude one of the matrix elements.
The self-consistent CPA equations are obtained from the conditions (13). We drop here all details and present only the final result (some intermediate expressions for thermal averages of different contributions listed in Table I 14), one can obtain the following expressions for the coefficients Φ α and Ψ α : and Ψ 2 (σ, g) = 2 45 − 4 45 σ 0 g 0 + 1 90 17σ 1 + 7σ 2 + 2 3 A. CPA solution for the paramagnetic phase In this section we consider CPA solutions for the PM phase of the degenerate DE model, and argue that the situation is qualitatively different from the one-orbital model, even in the case of cubic symmetry. Let us first consider the limit |z|→∞. In this case, the matrix elements of the Green function have the following asymptotic behavior: 49 g 0 (z)→ 1 z , g 1 (z)→ A z 2 , and g 2 (z)→ B z 2 ; and the CPA solution in the second order of 1 z can be easily obtained analytically from Eqs. (24)-(26) as σ 0 (z)→ 1 6z , σ 1 (z)→− 2B 3z 2 , and σ 2 (z)→− 7B 54z 2 . It is a single-valued solution and in this sense the situation is similar to the one-orbital case. The analysis is supported by numerical calculations for large but finite |z|, and represents a typical behavior when Im(z)≥0.75 (Fig. 5). The Van Hove singularities are smeared due to the large imaginary part of z and the self energy (Fig. 6).
However, when we start to approach the real axis, the situation changes dramatically. In the first calculations of such type we fix Im(z) and solve CPA equations by moving along the real axis and each time starting with the self-consistent self energy obtained for the previous value of Re(z). Then, in certain region of the complex plane we obtain two different solutions, depending on whether we move in the positive or negative direction of Re(z). A typical hysteresis loop corresponding to Im(z)=0.7 is shown in Figs. 5 and 6.
The behavior is related with the existence of two Van Hove singularities at the (0, π, 0) and (π, π, 0) points of the Brillouin zone, which are responsible for some sort of instability in the system. The singularities become increasingly important when z approaches the real axis. The positions of these singularities depend on matrix elements of the self energy, and given by Re(Σ 11 00 −3Σ 11 01 +Σ 22 01 ) and Re(Σ 11 00 +Σ 11 01 −3Σ 22 01 ), respectively. Therefore, by choosing different starting points for the self energy, the singularities can be shifted either 'to the right' or 'to the left' (Fig. 6). Since the CPA equations are non-linear, this may stabilize two different solutions.
A more complete picture can be obtained from Fig. 7, where we plot Re(Σ 22 01 −Σ 11 01 ) in the complex plane. 50 Depending on the location in the complex plane, the CPA equations (15) are again converged to either the same or two different solutions. The double-valued behavior of Re(Σ 22 01 −Σ 11 01 ) typically occurs within the shaded area. Note that this area is the result of numerical calculations which depend on the choice of the starting point. We do not exclude here the possibility that our result may be incomplete and that with a better choice of the starting conditions the double-valued area may be enlarged.
Our analysis is limited by Im(z)≈0.5. Further attempts to approach the real axis were conjugated with serious difficulties: the topology of solutions becomes increasingly complicated and may include several additional branches, some of which are presumably unphysical. At the present stage we do not have a clear strategy of how to deal with this problem. Of course, formally the situation can be regarded as the violation of causality principles. 28,29 Nevertheless, we would like to believe that the double-valued behavior of the Green function and the self energy obtained for Im(z)≥0.5 does have a physical meaning and is not an artifact of the model analysis. That is because of the following reasons: 1. for both solutions, the Green function satisfies the inequality ImG 11 00 ≤0; 2. the system is well defined in the whole interval of densities 0≤n(µ)≤1 (corresponding to the doping range 0≤x≤1 or the position of the chemical potential µ≤1) if and only if to take into account the superposition of the two solutions. Conversely, by considering only one of the solutions (and disregarding the other one as unphysical), the density will exhibit a finite jump (a discontinuity) at certain value of µ, and the system will be undefined within the discontinuity range. Presumably, the discontinuity of the electronic density would present even more unphysical behavior than the fact of the existence of two CPA solutions.
3. In some sense the appearance of two CPA solutions fits well into the logic of our work.
As it was discussed before, the non-local part of the self energy in the degenerate case is characterized by two different matrix elements and, in comparison with the one-orbital case, acquires an additional degree of freedom (which may be related with some non-local order parameter attached to the bond of the system). At the same time, the positions of Van Hove singularities, which control the behavior of non-linear CPA equations, depend on these matrix elements. Therefore, the appearance of several CPA solutions seems to be natural.
Note that Im(z)≈0.5 corresponds to the position of the first Matsubara pole for k B T ≈0.159. Taking into account the realistic value of the parameter (ddσ)≃0.7eV, 22 it roughly corresponds to T ≈1200K, which can be regarded as the lowest estimate for the temperature for which our analysis is strictly justified.
Below we discuss possible physical consequences of the existence of two CPA solutions in the PM state. Our scenario is based on the following observations (see Fig. 7): 1. The existence of the branch-point (B in Fig. 7), which forms two physical branches of CPA solutions in certain area of the complex plane. The requirement implies that there is a continuous path around the branch-point, which connects the points located on two different branches.
2. Both branches of the (multi-valued) Green function and the self energy are analytic (perhaps except the branch-point itself and the branch-edges). The requirement allows us to use the standard theorems of the contour integration: for example, the contour integral around the branch-point does not depend on the form of the contour, etc.
Strictly speaking, the second requirement is a postulate which is solely based on results of numerical calculations and at the present stage we do not have a general proof for it. Assume, it is correct. Then, the physical interpretation of the multi-valued behavior becomes rather straightforward and two CPA solutions can be linked to two PM phases (with different densities) corresponding to the same chemical potential µ. The situation has many things in common with the phenomenon of inhomogeneous phase separation, which was intensively discussed for manganites. 26,27 The new aspect in our case is that both phases are paramagnetic. The position of the branch-point itself can be related with the temperature, below which the PM state becomes intrinsicly inhomogeneous.

B. energy integration and appearance of two paramagnetic phases
In this section we discuss some aspects of the energy integration in the complex plane, related with the existence of two CPA-branches. Let us consider the integral: where X(µ) is a physical quantity, which can be the density of e g electrons, the double exchange energy or the change of either of them, and X(z) has the same topology in the complex plane as the self energy shown in Fig. 7. Then, the behavior of the integral (30) will depend on the position of the chemical potential µ with respect to the double-valued area. Generally, we should consider three possibilities (see Fig. 8 for notations): 1. µ<Υ 1 . In this case there is only one CPA solution and the integral (30) can be evaluated by using the standard techniques (see, e.g., Ref. 51 and references therein).
2. Υ 1 ≤µ≤Υ 2 . In this case the integral (30) takes two values for each value of the chemical potential µ: X 1 (µ), if the integrand is solely confined within one branch; and X 2 (µ)=X 1 (µ)+∆X(µ), if it is extended to the second branch. The discontinuity ∆X(µ) is given by the contour integral C 1 around the branch-point. It is important here that we do not try to define X(z) as a single-valued function by introducing the branch-cuts, which by itself is largely arbitrary procedure. 28 Instead, we treat both branches on an equal footing, which inevitably leads to the multi-valued behavior of X(µ). This integral (30) can be replaced by the contour integral C 2 spreading in the single-valued area plus residues calculated at Matsubara poles z n =µ+iπk B T (2n+1): In order to calculate X 1 and X 2 , the Matsubara poles should lie on the first and second branches of X(z), respectively. 52 3. µ>Υ 2 . In this case the integration along the real axis shall be combined with the discontinuity ∆X given by the contour integral C 1 . The integration can be replaced by Eq. (31). In this case, all Matsubara poles lie on the single branch. Therefore, X(µ) is a single-valued function.
As an illustration, we show in Fig. 9 the behavior of the averaged density as the function of µ. For T =0.24, corresponding to Im(z 0 )≃0.754, which is slightly above the branch-point, the first Matsubara pole falls beyond the double-valued area and n(µ) shows a normal behavior when for each value of µ there is only one value of the electronic density. For smaller T , n(µ) may take two different values for the same µ. This fact can be interpreted as the coexistence of two different phases. Such a behavior typically occurs in the interval −1.4≤µ≤−0.8, which may also depend on the temperature. As it was mentioned before, any single phase fails to define the electronic density in the whole interval 0≤n(µ)≤1 because of the discontinuity of n(µ). For example, for T =0.23 one of the phases is not defined in the interval 0.34≤n≤0.38, corresponding to the discontinuity at µ=−1.25, and the other one -for 0.46≤n≤0.66 corresponding to µ=−0.80. The problem can be resolved only by considering the combination of these two phases for which n is defined everywhere in the interval 0≤n(µ)≤1.
T P ≈0.23 (about 1800K), roughly corresponding to the position of the branch-point, can be regarded as the transition temperature to the two-phase state (a point of PM phase separation).

C. phase coexistence
In this section we briefly consider the problem of phase coexistence using a semiquantitative theory of non-interactive pseudoalloy. Namely, we assume that the free energy of the mixed paramagnetic state is given by where E (1) D and E (2) D are the DE energies of two PM phases (with lower and higher densities of the e g electrons, respectively); y is the "alloy concentration"; and S mix (y) is the configurational mixing entropy: Then, the equilibrium concentration, which minimizes the free energy (32) is given by D can be calculated using the definition (7) and the formula (31) for the contour integration. Since the contour C 2 is confined within the single-valued area, ∆E D is given by the difference of residues at a limited number of Matsubara poles: Some details of these calculations can be found in Appendix C. The discontinuity of the electronic density ∆n=n (2) −n (1) can be calculated in a similar way. Results of these calculations are shown in Fig. 10. Particularly, the equilibrium alloy concentration y is close to 0.5, meaning that the difference ∆E D is small and the main contribution to the free energy comes from the entropy term. The electronic density n =(1−y)n (1) +yn (2) averaged simultaneously over the spin orientations and the alloy concentrations shows a discontinuity at the edges of the two-phase region, meaning that the system is not defined. This behavior is unphysical and caused by the non-interactive approach to the problem of phase coexistence. The averaged densities range, for which two PM phases may coexist is typically 0.3-0.7 and depends on the temperature.

D. Curie temperature
The Curie temperature can be obtained from Eq. (18). In the case of degenerate DE model, the function D(T ) is given by: 53 Results of these calculations are shown in Fig. 11. For µ<−1, the magnetic transition temperature appears to be lower than the point of PM phase separation (T P ). Taking into account that T P ≃0.23 and using results of Fig. 10, µ<−1 roughly corresponds to the densities n <0.5. In this region, T C should be calculated independently for two different phases. The calculations can be done using Eq. (31). Not surprisingly that different phases are characterized by different T C 's. Therefore, for the degenerate DE model we expect the existence of two magnetic transition points. With the cooling down of the sample, the transition to the FM state takes place first in one of the phases, characterized by lower density (the hole-rich phase). Then, within the interval T (1) C <T < T (2) C , the FM phase continues to coexist with the PM phase, persisting in the hole-deficient part of the sample. Taking into account that (ddσ)≈0.7eV, the difference of two transition temperatures, which depend on µ, can be evaluated as 0<T (2) C the system exhibits the proper FM order established in both phases.
Formally, the opposite scenario when the phase separation occurs below the magnetic transition temperature is also possible, and according to Fig. 11 may take place when µ>−1 ( n ≥0.5). However, the quantitative description of this situation is beyond the small-λ limit, considered in the present work.

VI. CONCLUDING REMARKS
We have applied a non-local CPA approach to the problem of orientational spin disorder in the double exchange systems, which was supplemented by a mean-field theory for the analysis of magnetic transition temperature. Our CPA approach is based on the perturbation theory expansion for the T -matrix with respect to fluctuations of hoppings from the mean value of the DE Hamiltonian specified by the matrix elements of the (non-local) self energy, so that in the lowest (first) order it automatically recovers the main results of the DE theory by de Gennes (the bandwidth narrowing in the PM state, expression for the Curie temperature, etc.). 3 Our main focus was on the correction of this theory by higher-order effects with respect to the fluctuations, which were included up to the second order and treated in the real space.
In the one-orbital case, it led to a substantial reduction of T C (up to 20%). Nevertheless, the obtained value of T C was largely overestimated in comparison with results of Monte Carlo calculations, 48 due to limitations inherent to the mean-field approach. Therefore, a sensible description of spatial spin correlations, beyond the mean-field approximation, presents a very important direction for the improvement of our model.
It appeared, however, that even on the level of mean-field theories the situation is far from being fully investigated, especially if one takes into account the effects of orbital degeneracy and details of realistic electronic structure for the e g electrons in the cubic perovskite structure. Particularly, two Van Hove singularities at the (π, π, 0) and (0, π, 0) points of the Brillouin zone may contribute to the properties of DE systems not only in the ordered FM state, 25 but also in the case of the spin disorder, above the magnetic transition temperature. The singularities lead to the branching of CPA solutions for the PM state, so that for certain values of the chemical potential the system consists of two PM phases with two different densities. In such a situation, the PM state becomes intrinsicly inhomogeneous, that also determines details of the phase transition to the FM state with the decrease of the temperature. The magnetic transition is characterized by the existence of two transition points, and first takes place in one of the phases which constitutes only a fraction of the sample. The FM transition in the second phase takes place at yet another temperature, which can be significantly lower than the first one. Both magnetic transitions are continuous (of the second order). However, they occur separately in two different phases which already exist in the PM region, i.e. above the first magnetic transition point. 54 How is this scenario consistent with the experimental behavior of perovskite manganites? It is true that there is no clear experimental evidence supporting the existence of an individual temperature of the PM phase separation, T P , and the intrinsic inhomogeneity of the PM phase (unless it is caused by external factors such as the chemical and structural inhomogeneities and the grain boundaries). 55 In this sense our result can be regarded as the prediction. On the other hand, the phase coexistence below T C is rather common, and was observed in a number of experiments. [15][16][17] In addition, our result naturally explains the appearance of several magnetic transition points in perovskite manganites. Yet, a clear difference is that experimentally one (or sometimes both) magnetic transitions are antiferromagnetic (typically either to the A-or CE-type AFM state), and the situation when two consecutive transitions go to the FM state is not realized in practice. 24 Presumably, the difference is caused by the limitation of our analysis by the PM and FM states, while in reality for n <0.7 (i.e. when several magnetic transition points are expected) the formation of the (A-, CE-, and C-type) AFM structures seems to be more natural. From this point of view, results of our work present mainly an academical interest: we pointed out at the principal possibility of the existence of several magnetic transition temperatures for the DE systems, however the type of the ordered state considered at the low temperature was not sufficiently general. 56 Finally, we would like to discuss briefly some possible extensions of our model.

1.
As it was already mentioned, any realistic description of the phase diagram of perovskite manganites would be incomplete without the (A-, CE-, and C-type) AFM structures. The minimal model which captures the behavior of doped manganites at the low temperature is the double exchange combined with the isotropic AFM SE interaction between the localized spins (J S ). So, the magnetic phase diagram at T =0 can be understood in terms of the anisotropy of interatomic DE interactions, caused by the anisotropy of magnetic ordering and operating in the background of isotropic AFM SE interactions. 22 In this picture, the SE interactions were needed to shift the reference point simultaneously in all bonds in the direction of the AFM coupling, while the variety of the phase diagram itself is described by the anisotropy of the DE interactions. However, similar shift is not applicable for the magnetic transition temperature, or may not be the major effect of the SE interactions. For example, the shift of the Curie temperature by J S , T C →T C −2|J S |, 47 will also affect the mesh of Matsubara poles in Eq. (31). It will require the careful analysis of the topology of CPA solutions near the real axis, which will present the main obstacle for such calculations.
2. Our main results regarding the nontrivial topology of the CPA solutions for the PM state have been based on the second order of the perturbation theory expansion for the T -matrix. Since the effect was so dramatic, it naturally rises the question about importance of the higher-order terms. It seems to be a very important problem for the future analysis.
3. The main effect discussed in this work is based on the existence of Van Hove singularities in the density of states of degenerate DE model. However, the exact position of these singularities, or even the fact of their existence in realistic compounds depend on many other factors, such as the Mn(3d)-O(2p) hybridization, the cation and structural disorder, the purity of sample, etc. All these factors may significantly alter conclusions of our work, if one try to apply them to realistic compounds.
ACKNOWLEDGMENTS I thank F. Aryasetiawan for discussion of the energy integration around the branch-point in the complex plane; Y. Tomioka -for discussion of the experimental behavior of manganites, and A. I. Liechtenstein -for discussion of causality problems in non-local CPA. The present work is partly supported by New Energy and Industrial Technology Development Organization (NEDO).  Table I) become: ∆H 00 G 01 ∆H 11 = σ 2 0 g 1 , ∆H 01 G 10 ∆H 01 = σ 2 1 + Due to degeneracy, it is sufficient to consider the contributions to only one site-diagonal element of the averaged T -matrix, say T 11 00 (z). Results of the orientational averaging for T 22 00 (z) will be identical. Then, for different contributions listed in Table I we have: For the bond 0-1, results of the orientational averaging of the diagonal 11 and 22 matrix elements will be different. The contributions to T 11 01 (z) are given by: [∆H 01 G 10 ∆H 01 ] 11 = σ 2 1 g 1 , and while similar contributions to T

APPENDIX C: CHANGE OF THE INTEGRATED DENSITY OF STATES IN THE DOUBLE-VALUED REGION
In this appendix we discuss some practical aspects of calculations of the difference n (2) (z)−n (1) (z) between two CPA solutions in the PM states. According to Ducastelle, 34 n(z) is given by the following expression: where Tr is the trace over site and orbital indices, N is the number of atomic sites, and the hat-symbols here stand for the matrices in the orbital and atomic coordinates space. Let us begin with the first term. In the second order of G (2) (z)− G (1) (z) we have: which can be further transformed using the definition of the Green function (10) as The inverse matrix {...} −1 can be calculated in the same way as the Green function (10). Then, if R 0 (z) is the site-diagonal element of {...} −1 , and R 1 (z) and R 2 (z) are the siteoff-diagonal ones corresponding to the x 2 −y 2 and 3z 2 −r 2 states for the bond 0-1, we can write where ∆σ α (z)=σ (2) α (z)−σ (1) α (z). In the second term of Eq. (C1) (the so-called vortex correction) we expand ln up to the second order of [ H− H(z)]. This expansion is necessary to preserve the variational properties of our CPA formalism, 34 which is based on the same approximation (14) for the T -matrix. Thus, we have In order to calculate the thermal average of this expression, we note that which immediately follows from the CPA equations (13) under condition (14). Thus, and corresponding contribution to the integrated density of states is given by TABLES TABLE I. Different contributions to the site-diagonal (T 00 ) and site-off-diagonal elements (T 01 ) of the T -matrix in the second order with respect to the fluctuations ∆H ij =H ij −H ij in the real space. The column 'comment' is used to explain the position of intermediate sites used in the summation in Fig. 2. contribution comment element T 00 : ∆H 00 ∆H 00 G 00 ∆H 00 , (4,8), (5,9), (6,10) (π,π,0) (0,0,0) (π,π,0) (0,π,0) (π,π,π) (0,π,0) Energy FIG. 1. Tight-binding density of states and the dispersion of two e g bands in the ferromagnetic state (in units of ddσ transfer integral). The dotted line shows the positions of the Fermi level as a function of hole-concentration x, which is related with the integrated density of states as x=1−n. Note the existence of two Van Hove singularities at (π, π, 0) and (0, π, 0), responsible for the kinks of density of states at ±1. The first singularity is located near the Fermi level when x≃0.3 (shown by arrow).    (b): An equivalent expression in terms of the contour integral C 2 spreading in the single-valued area of the complex plane and residues calculated at Matsubara poles. The latter contributions are different for two different branches, that is equivalent to the discontinuity term in the scheme (a). D , the equilibrium alloy concentration y, and the density of e g electrons averaged simultaneously over the spin orientations and the alloy concentrations n =(1−y)n (1) +yn (2) as a function of chemical potential µ for three different temperatures (in units of ddσ-integral). The superscripts (1) and (2) stand for the phases correspondingly with lower and higher electronic densities. (2) FIG. 11. Main transition temperatures for the degenerate double exchange model (in units of ddσ-integral). T P is the transition temperature to the two-phase paramagnetic state. The shaded area shows the approximate range of the chemical potentials (µ) when the paramagnetic state becomes intrinsicly inhomogeneous. T (1) C and T (2) C are the Curie temperatures for two different phases (characterized by lower and higher densities of the e g electrons, respectively).