text:
# References

Chiral 3$`\pi `$-exchange NN-potentials: Results for dominant next-to-leading order contributions

N. Kaiser

Physik Department T39, Technische Universität München,
D-85747 Garching, Germany

## Abstract

We calculate in (two-loop) chiral perturbation theory the local NN-potentials generated by the three-pion exchange diagrams with one insertion from the second order chiral effective pion-nucleon Lagrangian proportional to the low-energy constants $`c_{1,2,3,4}`$. The resulting isoscalar central potential vanishes identically. In most cases these $`3\pi `$-exchange potentials are larger than the ones generated by the diagrams involving only leading order vertices due to the large values of $`c_{3,4}`$ (which mainly represent virtual $`\mathrm{\Delta }`$-excitation). A similar feature has been observed for the chiral $`2\pi `$-exchange. We also give suitable (double-integral) representations for the spin-spin and tensor potentials generated by the leading-order diagrams proportional to $`g_A^6`$ involving four nucleon propagators. In these cases the Cutkosky rule cannot be used to calculate the spectral-functions in the infinite nucleon mass limit since the corresponding mass-spectra start with a non-vanishing value at the $`3\pi `$-threshold. Altogether, one finds that chiral $`3\pi `$-exchange leads to small corrections in the region $`r1.4`$ fm where $`1\pi `$\- and chiral $`2\pi `$-exchange alone provide a very good strong NN-force as shown in a recent analysis of the low-energy pp-scattering data-base.

PACS: 12.20.Ds, 12.38.Bx, 12.39.Fe, 13.75.Cs.

To be published in The Physical Review C (2000)

Over the last years effective field theory methods have been successfully applied to the two-nucleon system at low and intermediate energies . The idea of constructing the NN-potential from effective field theory was put forward by Weinberg  and this was taken up by van Kolck and collaborators . Their NN-potential contained one- and two-pion exchange graphs based on time-ordered perturbation theory (regularized by a Gaussian cut-off) and a host of contact interactions. The method of unitary transformations was used by Epelbaoum, Glöckle and Meißner  to construct an energy-independent NN-potential from the effective chiral pion-nucleon Lagrangian. Based on one- and two-pion exchange and a few contact interactions which contribute only to S- and P-waves a good description of the deuteron properties as well as the NN phase-shifts and mixing angles below $`T_{lab}=300`$ MeV was found in that framework .

The issue of the chiral $`3\pi `$-exchange was first addressed by Pupin and Robilotta  who considered one particularly simple diagram with no internal nucleon propagators. In a recent work  we have started to calculate the static NN-potentials generated by the (two-loop) $`3\pi `$-exchange diagrams with all possible interaction vertices taken from the leading-order effective chiral $`\pi N`$-Lagrangian. In ref. we restricted ourselves to the evaluation of four relatively simple classes of diagrams with the common prefactor $`g_A^2/f_\pi ^6`$. It was stressed in ref. that the chiral $`3\pi NN`$-contact vertex employed in ref. depends on the choice of the interpolating pion-field and therefore one has to consider representation invariant classes of diagrams by supplementing graphs involving the chiral $`4\pi `$-vertex. Another obvious consequence of this is that there is no unique definition of the so-called ”pion-nucleon form factor” which is often introduced in phenomenological models. The only meaningful separation of the NN-potential is the one into point-like $`1\pi `$-exchange and $`3\pi `$-exchange. The four classes of $`3\pi `$-exchange diagrams evaluated in ref. gave rise only to isovector spin-spin and tensor potentials with a strength of $`\pm 0.1`$ MeV or less at $`r=1.0`$fm. Compared to the chiral $`2\pi `$-exchange NN-potentials  these are negligibly small corrections.

In a subsequent work  the static NN-potentials generated by the $`3\pi `$-exchange diagrams proportional to $`g_A^4`$ and $`g_A^6`$ (grouped into classes V, VI, VII, VIII and IX) have been worked out. The resulting coordinate space potentials have typical strengths of a few MeV at $`r=1.0`$ fm. A surprising result of ref. was that the total isoscalar central $`3\pi `$-exchange NN-potential vanishes identically. The static $`3\pi `$-exchange NN-potentials have been constructed in refs. from the mass-spectra (i.e. the imaginary parts) of the respective two-loop diagrams in the infinite nucleon mass limit $`M\mathrm{}`$. This method could however not be successfully applied to the spin-spin and tensor potentials generated by the diagrams proportional to $`g_A^6`$ (classes VIII and IX) involving four nucleon propagators. In these cases the integrand of the $`3\pi `$-phase space integral $`_{z^21}𝑑\omega _1𝑑\omega _2\mathrm{}`$ developed a non-integrable singularity of the form $`(1z^2)^{3/2}`$ on the boundary $`z^2=1`$. By considering the analytically solvable example of the crossed three-boson exchange with scalar coupling in ref. the failure of the spectral-function method could be understood from the fact that in the limit $`M\mathrm{}`$ the mass-spectrum starts with a non-vanishing value at the 3-boson threshold $`\mu =3m`$. In coordinate space this feature manifests itself as a large-$`r`$ asymptotics of the form $`e^{3mr}r^2`$.

The purpose of this work is to present suitable double-integral representations for these remaining $`3\pi `$-exchange spin-spin and tensor NN-potentials proportional to $`g_A^6`$ together with numerical results. Furthermore, will work out here the next-to-leading order corrections to the chiral $`3\pi `$-exchange NN-potential. In general, the next-to-leading order contributions are given by all diagrams with one insertion from the second order chiral $`\pi N`$-Lagrangian $`_{\pi N}^{(2)}`$. However, from an analogous calculation of the chiral $`2\pi `$-exchange in ref. one has experienced that the dominant contributions come from the second order chiral $`\pi \pi NN`$-contact vertex proportional to the low-energy constants $`c_{1,2,3,4}`$, whereas the relativistic $`1/M`$-corrections are much smaller. The exception to this is, of course, the spin-orbit NN-potential which is a truly relativistic effect vanishing in the infinite nucleon mass limit $`M\mathrm{}`$. Therefore we will consider here only the two-loop $`3\pi `$-exchange diagrams with exactly one insertion proportional to $`c_{1,2,3,4}`$. The corresponding effective chiral $`\pi N`$-Lagrangian reads

$$_{\pi N}^{(2)}=\overline{N}\left\{2c_1m_\pi ^2(U+U^{})+c_2u_0^2+c_3u_\nu u^\nu +\frac{i}{2}c_4\stackrel{}{\sigma }(\stackrel{}{u}\times \stackrel{}{u})\right\}N,$$
(1)
with the axial vector quantity $`u^\nu =i\{\xi ^{},^\nu \xi \}`$. The unitary matrix $`\xi ^2=U=1+i\stackrel{}{\tau }\stackrel{}{\pi }/f_\pi \stackrel{}{\pi }^2/2f_\pi ^2+\mathrm{}`$ collects the Goldstone pion-fields and $`f_\pi =92.4`$ MeV denotes the weak pion decay constant. The low-energy constants $`c_{1,2,3,4}`$ have been determined in a one-loop chiral perturbation theory analysis of low-energy $`\pi N`$-scattering in the physical region  and inside the Mandelstam triangle . Their central values are $`c_1=0.8`$ GeV-1, $`c_2=3.2`$ GeV-1, $`c_3=4.7`$ GeV-1 and $`c_4=3.4`$ GeV-1. We prefer to use the values obtained in ref. since at the small sub-threshold pion-energies relevant inside the Mandelstam triangle the chiral expansion of the $`\pi N`$-scattering amplitude is expected to work best. As outlined in ref. the large low-energy constants $`c_{2,3,4}`$ receive dominant contributions from virtual $`\mathrm{\Delta }(1232)`$-isobar excitation. In ref. it has been demonstrated that the $`2\pi `$-exchange NN-potentials generated by explicit $`\mathrm{\Delta }(1232)`$-isobar excitation are very similar to those obtained from the $`NN\pi \pi `$-contact interaction  proportional to $`c_{1,2,3,4}`$. In the so-called small scale expansion of ref. the $`\mathrm{\Delta }N`$ mass-splitting $`\mathrm{\Delta }=293`$ MeV is treated as a further small expansion parameter and $`2\pi `$-exchange graphs with $`\mathrm{\Delta }`$-excitation are counted as leading order ones. This is another reason to focus on the $`3\pi `$-exchange potential proportional to $`c_{1,2,3,4}`$. The second order chiral $`\pi \pi NN`$-contact vertex following from eq.(1) reads

$$\frac{i}{f_\pi ^2}\left\{2\delta ^{ab}\left[2c_1m_\pi ^2+c_2q_a^0q_b^0+c_3q_aq_b\right]+c_4ϵ^{abc}\tau ^c\stackrel{}{\sigma }(\stackrel{}{q}_a\times \stackrel{}{q}_b)\right\},$$
(2)
where both pion four-momenta $`q_a`$ and $`q_b`$ are out-going ones. Note that the isoscalar part in eq.(2) is spin-independent while the isovector part is spin-dependent.

Next, let us give some basic definitions in order to fix our notation. In the static limit and considering only irreducible diagrams the on-shell NN T-matrix takes the following form

$`𝒯_{NN}`$ $`=`$ $`V_C(q)+V_S(q)\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2+V_T(q)\stackrel{}{\sigma }_1\stackrel{}{q}\stackrel{}{\sigma }_2\stackrel{}{q}`$ (3)
$`+\left[W_C(q)+W_S(q)\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2+W_T(q)\stackrel{}{\sigma }_1\stackrel{}{q}\stackrel{}{\sigma }_2\stackrel{}{q}\right]\stackrel{}{\tau }_1\stackrel{}{\tau }_2,`$
where $`q=|\stackrel{}{q}|`$ denotes the momentum transfer between the initial and final-state nucleon. The subscripts $`C,S`$ and $`T`$ refer to the central, spin-spin and tensor components, each of which occurs in an isoscalar ($`V_{C,S,T}`$) and an isovector version ($`W_{C,S,T}`$). As indicated, the (real) NN-amplitudes $`V_C(q),\mathrm{},W_T(q)`$ depend only on the momentum transfer $`q`$ in the static limit $`M\mathrm{}`$. We are here interested only in the coordinate-space potentials generated by certain diagrams in which three pions are simultaneously exchanged between both nucleons. For this purpose it is sufficient to calculate the imaginary parts of the NN-amplitudes, Im $`V_{C,S,T}(i\mu )`$ and Im $`W_{C,S,T}(i\mu )`$, analytically continued to time-like momentum transfer $`q=i\mu 0^+`$ with $`\mu 3m_\pi `$. These imaginary parts are then the mass-spectra entering a representation of the local coordinate-space potentials in the form of a continuous superposition of Yukawa-functions,

$`\stackrel{~}{V}_C(r)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2r}}{\displaystyle _{3m_\pi }^{\mathrm{}}}𝑑\mu \mu e^{\mu r}\mathrm{Im}V_C(i\mu ),`$ (4)
$`\stackrel{~}{V}_S(r)`$ $`=`$ $`{\displaystyle \frac{1}{6\pi ^2r}}{\displaystyle _{3m_\pi }^{\mathrm{}}}𝑑\mu \mu e^{\mu r}\left[\mu ^2\mathrm{Im}V_T(i\mu )3\mathrm{Im}V_S(i\mu )\right],`$ (5)
$`\stackrel{~}{V}_T(r)`$ $`=`$ $`{\displaystyle \frac{1}{6\pi ^2r^3}}{\displaystyle _{3m_\pi }^{\mathrm{}}}𝑑\mu \mu e^{\mu r}(3+3\mu r+\mu ^2r^2)\mathrm{Im}V_T(i\mu ).`$ (6)
The isoscalar central, spin-spin and tensor potentials, denoted here by $`\stackrel{~}{V}_{C,S,T}(r)`$, are as usual those ones which are accompanied by the operators $`1`$, $`\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2`$ and $`3\stackrel{}{\sigma }_1\widehat{r}\stackrel{}{\sigma }_2\widehat{r}\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2`$, respectively. For the isovector potentials $`\stackrel{~}{W}_{C,S,T}(r)`$ a completely analogous representation holds. The imaginary parts entering eqs.(4,5,6) are calculated from the pertinent 2-loop $`3\pi `$-exchange diagrams as integrals of the $`\overline{N}N3\pi \overline{N}N`$ transition amplitudes over the Lorentz-invariant $`3\pi `$-phase space making use of the Cutkosky cutting rule as explained in refs.. The sub-leading $`3\pi `$-exchange considered here occurs at fifth order in small momentum (or chiral) expansion. Of course at this order there are also many two-loop $`2\pi `$-exchange diagrams proportional to $`c_{1,2,3,4}`$. Since the corresponding spectral-functions will (in general) contain one-loop sub-divergences (related to the fourth order chiral $`\pi N`$-amplitude) an adequate discussion of these $`2\pi `$-exchange potentials requires good knowledge of the strength of the counter-terms $`\overline{e}_i`$ in $`_{\pi N}^{(4)}`$ . With these included the corresponding $`2\pi `$-exchange potential is well-defined, however, its explicit computation goes far beyond the scope of the present paper.

We start with the diagrams of class X shown in Fig.1. The heavy dot symbolizes the second order chiral $`\pi \pi NN`$-contact vertex eq.(2) proportional to $`c_{1,2,3,4}`$. Diagrams for which the role of both nucleons is interchanged are not shown. These lead to the same contribution to the NN-potential, i.e. a factor of 2. As stressed in ref. diagrams involving the chiral $`3\pi NN`$-vertex or the chiral $`4\pi `$-vertex depend on an arbitrary parameter and therefore one should consider the full class X as one entity. Obviously, the last two pion-pole diagrams in Fig.1 contribute via coupling constant renormalization also to the point-like $`1\pi `$-exchange. This effect is however automatically taken care by working with the physical (on-shell) $`\pi NN`$-coupling constant $`g_{\pi N}`$.

Fig.1: $`3\pi `$-exchange diagrams of class X proportional to $`g_A^2`$. Solid and dashed lines represent nucleons and pions, respectively. The heavy dot symbolized an insertion from second order chiral $`\pi N`$-Lagrangian. The combinatoric factor of these diagrams is 1/2. Diagrams for which the role of both nucleons is interchanged are not shown. They lead to the same NN-potential.

Note that we are considering here only the finite-range Yukawa-parts of the $`1\pi `$-exchange and we disregard all zero-range $`\delta ^3(\stackrel{}{r})`$-terms. Transformed into momentum space the latter become polynomials in $`q^2`$ with possible contributions from higher-derivative operators. From an inspection of the spin- and isospin factors occurring in the diagrams of class X one finds immediately that only non-vanishing isovector spin-spin and tensor NN-amplitudes $`W_{S,T}`$ will be obtained. We find the following imaginary parts from class X (dropping from now on the argument $`i\mu `$),

$`\mathrm{Im}W_S^{(X)}`$ $`=`$ $`{\displaystyle \frac{g_A^2}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{[2C_{13}(w)+{\displaystyle \frac{c_2}{3}}(w^24m_\pi ^2)]\lambda (w)`$ (7)
$`+{\displaystyle \frac{c_4}{3}}(w^24m_\pi ^2)[(w^2m_\pi ^2)^2+\mu ^2(2m_\pi ^2+2w^23\mu ^2)]\},`$
$`\mathrm{Im}W_T^{(X)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}W_S^{(X)}+{\displaystyle \frac{g_A^2}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{[4C_{13}(w)+{\displaystyle \frac{2c_2}{3}}(w^24m_\pi ^2)]`$ (8)
$`\times (\mu ^2+m_\pi ^2w^2)\left[\mu ^2+m_\pi ^2w^2+2\mu ^2(2w^2\mu ^2)(\mu ^2m_\pi ^2)^1\right]`$
$`+{\displaystyle \frac{2c_4}{3}}(w^24m_\pi ^2)\lambda (w)(m_\pi ^2+3\mu ^2)(m_\pi ^2\mu ^2)^1\},`$
where we have introduced the abbreviations $`C_{13}(w)=4c_1m_\pi ^2+c_3(w^22m_\pi ^2)`$ and $`\lambda (w)=w^4+\mu ^4+m_\pi ^42w^2\mu ^22w^2m_\pi ^22\mu ^2m_\pi ^2`$. The variable $`w`$ denotes the invariant mass of a pion-pair and its kinematically allowed range is $`2m_\pi w\mu m_\pi `$. The $`dw`$-integrals in eqs.(7,8) could of course be solved easily in terms of square-root and logarithmic functions. However, we want to avoid the resulting rather lengthy expressions.

Fig.2: $`3\pi `$-exchange diagrams of class XI proportional to $`g_A^2`$. These give a vanishing contribution to the NN-potential.

Next, we consider the diagrams of class XI shown in Fig.2. Each of them gives a vanishing contribution to the NN T-matrix for the following reason. Since the leading order (Weinberg-Tomozawa) $`\pi \pi NN`$-vertex is of isovector nature (i.e. proportional to $`ϵ^{abc}\tau ^c`$) a non-zero isospin-factor is obtained only from the spin-dependent isovector term in eq.(2) proportional to $`c_4`$. When combined with the spin-independent Weinberg-Tomozawa vertex on the other side of the bubble-type sub-diagram one finds immediately that the pertinent 1-loop integral is zero in the heavy baryon formalism (basically because the spin-vector $`\stackrel{}{\sigma }`$ has no time-component).

Fig.3: $`3\pi `$-exchange diagrams of class XII proportional to $`g_A^2`$. For further notation, see Fig. 1.

Next, we consider the diagrams of class XII shown in Fig.3. The isoscalar contribution comes exclusively from the $`c_4`$-term in eq.(2). Altogether one obtains the following imaginary parts of the isoscalar and isovector spin-spin and tensor NN-amplitudes from class XII,

$$\mathrm{Im}V_S^{(XII)}=\frac{2g_A^2c_4}{(4f_\pi )^6\pi ^2\mu ^3}_{2m_\pi }^{\mu m_\pi }𝑑w\sqrt{w^24m_\pi ^2}\lambda (w)\left[7m_\pi ^2\mu ^23w^2+2m_\pi ^2(m_\pi ^2\mu ^2)w^2\right],$$
(9)
$`\mathrm{Im}V_T^{(XII)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}V_S^{(XII)}+{\displaystyle \frac{4g_A^2c_4}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}𝑑w\sqrt{w^24m_\pi ^2}(\mu ^2+m_\pi ^2w^2)`$ (10)
$`\times \left[3w^42w^2(5m_\pi ^2+2\mu ^2)+\mu ^4+2\mu ^2m_\pi ^2+5m_\pi ^4+2m_\pi ^2(\mu ^2m_\pi ^2)^2w^2\right],`$
$`\mathrm{Im}W_S^{(XII)}`$ $`=`$ $`{\displaystyle \frac{g_A^2}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{[2C_{13}(w)+{\displaystyle \frac{c_2}{3}}(w^24m_\pi ^2)]\lambda (w)`$ (11)
$`+{\displaystyle \frac{c_4}{3}}(w^24m_\pi ^2)[(\mu ^2m_\pi ^2)(3\mu ^2+m_\pi ^22w^2)w^4]\},`$
$`\mathrm{Im}W_T^{(XII)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}W_S^{(XII)}+{\displaystyle \frac{g_A^2}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{[(w^2m_\pi ^2)^2\mu ^4]`$ (12)
$`\times [4C_{13}(w)+{\displaystyle \frac{2c_2}{3}}(w^24m_\pi ^2)]{\displaystyle \frac{2c_4}{3}}(w^24m_\pi ^2)\lambda (w)\}.`$
Fig.4: $`3\pi `$-exchange diagrams of class XIII proportional to $`g_A^4`$. For further notation, see Fig. 1.

Next, we turn to the diagrams of class XIII shown in Fig.4. The technique to separate off the iterative part from the first two reducible diagrams has been explained in ref.. The iterative (or reducible) part is defined here (entirely) within perturbation theory as that part which carries in the numerator the large scale enhancement factor $`M`$, the nucleon mass. Such an iterative part does therefore not obey the naive chiral power counting rules (see also section 4.3 in ref. on the so-called iterated $`1\pi `$-exchange). A detailed discussion of possible ambiguities showing up in (non-perturbative) iterations to infinite orders (e.g. via a Schrödinger equation) can be found in ref.. The isovectorial spin-dependent contact vertex proportional to $`c_4`$ (and in fact only this one) produces now also a central NN-amplitude. Interestingly, its isoscalar and isovector components come with a fixed ratio. The corresponding imaginary parts read

$$\mathrm{Im}V_C^{(XIII)}=\frac{3}{4}\mathrm{Im}W_C^{(XIII)}=\frac{12g_A^4c_4}{(4f_\pi )^6\pi ^2\mu }_{2m_\pi }^{\mu m_\pi }𝑑w\sqrt{w^24m_\pi ^2}\lambda (w).$$
(13)
Furthermore, the isoscalar and isovector spin-spin and tensor NN-amplitudes generated by the diagrams of class XIII have the following imaginary parts,

$`\mathrm{Im}V_S^{(XIII)}`$ $`=`$ $`{\displaystyle \frac{g_A^4c_4}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{7\mu ^637\mu ^4m_\pi ^2+45\mu ^2m_\pi ^415m_\pi ^6+9w^6`$ (14)
$`w^4(11\mu ^2+37m_\pi ^2)+w^2(45m_\pi ^4+32\mu ^2m_\pi ^25\mu ^4)+2m_\pi ^2(\mu ^2m_\pi ^2)^3w^2\},`$
$`\mathrm{Im}V_T^{(XIII)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}V_S^{(XIII)}+{\displaystyle \frac{2g_A^4c_4}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}𝑑w\sqrt{w^24m_\pi ^2}`$ (15)
$`\times \{9w^6w^4(19\mu ^2+37m_\pi ^2)+w^2(11\mu ^4+52\mu ^2m_\pi ^2+45m_\pi ^4)`$
$`\mu ^613\mu ^4m_\pi ^235\mu ^2m_\pi ^415m_\pi ^6+2m_\pi ^2(m_\pi ^2\mu ^2)(\mu ^4m_\pi ^4)w^2\},`$
$`\mathrm{Im}W_S^{(XIII)}`$ $`=`$ $`{\displaystyle \frac{g_A^4}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{3C_{13}(w)[4\mu ^2(\mu ^2m_\pi ^2w^2)\lambda (w)]`$ (16)
$`+{\displaystyle \frac{c_2}{2}}(8m_\pi ^25w^2)\lambda (w)+4c_4[w^6w^4(\mu ^2+4m_\pi ^2)`$
$`+w^2(\mu ^42\mu ^2m_\pi ^2+5m_\pi ^4)\mu ^6+4\mu ^4m_\pi ^2\mu ^2m_\pi ^42m_\pi ^6]\},`$
$`\mathrm{Im}W_T^{(XIII)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}W_S^{(XIII)}+{\displaystyle \frac{g_A^4}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{(m_\pi ^2+\mu ^2w^2)`$ (17)
$`\times \left[6C_{13}(w)(3\mu ^2+w^2m_\pi ^2)+c_2(8m_\pi ^25w^2)(m_\pi ^2+\mu ^2w^2)\right]`$
$`+8c_4(w^22m_\pi ^2)[(w^2m_\pi ^2)^2\mu ^4]\}.`$
Fig.5: $`3\pi `$-exchange diagrams of class XIV proportional to $`g_A^4`$. For further notation, see Fig. 1.

Finally, we consider the (irreducible) diagrams of class XIV shown in Fig.5. In this case the isovectorial $`c_4`$-term in eq.(2) does not make a contribution to the isovector NN-amplitudes $`W_{C,S,T}`$. One obtains from class XIV an isoscalar central NN-amplitude which is however exactly canceled by the one coming from class XIII,

$$\mathrm{Im}V_C^{(XIV)}=\mathrm{Im}V_C^{(XIII)}.$$
(18)
Therefore, one has the interesting result that there is no isoscalar central NN-potential from chiral $`3\pi `$-exchange even at next-to-leading order (considering the dominant $`c_{1,2,3,4}`$-terms). For the imaginary parts of the spin-spin and tensor NN-amplitudes generated by the diagrams of class XIV one obtains the following expressions,

$`\mathrm{Im}V_S^{(XIV)}`$ $`=`$ $`{\displaystyle \frac{g_A^4c_4}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{9w^65\mu ^613\mu ^4m_\pi ^2+33\mu ^2m_\pi ^415m_\pi ^6`$ (19)
$`w^4(23\mu ^2+37m_\pi ^2)+w^2(45m_\pi ^4+8\mu ^2m_\pi ^2+19\mu ^4)+2m_\pi ^2(\mu ^2m_\pi ^2)^3w^2\},`$
$`\mathrm{Im}V_T^{(XIV)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}V_S^{(XIV)}+{\displaystyle \frac{2g_A^4c_4}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}𝑑w\sqrt{w^24m_\pi ^2}`$ (20)
$`\times \{9w^6w^4(19\mu ^2+37m_\pi ^2)+w^2(11\mu ^4+52\mu ^2m_\pi ^2+45m_\pi ^4)`$
$`\mu ^613\mu ^4m_\pi ^235\mu ^2m_\pi ^415m_\pi ^6+2m_\pi ^2(m_\pi ^2\mu ^2)(\mu ^4m_\pi ^4)w^2\},`$
$`\mathrm{Im}W_S^{(XIV)}`$ $`=`$ $`{\displaystyle \frac{g_A^4}{(4f_\pi )^6\pi ^2\mu ^3}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}𝑑w\sqrt{w^24m_\pi ^2}\lambda (w)`$ (21)
$`\times \left\{C_{13}(w)+{\displaystyle \frac{c_2}{6}}\left[w^2+2\mu ^26m_\pi ^2+4m_\pi ^2(\mu ^2m_\pi ^2)w^2\right]\right\},`$
$`\mathrm{Im}W_T^{(XIV)}`$ $`=`$ $`{\displaystyle \frac{1}{\mu ^2}}\mathrm{Im}W_S^{(XIII)}+{\displaystyle \frac{g_A^4}{(4f_\pi )^6\pi ^2\mu ^5}}{\displaystyle _{2m_\pi }^{\mu m_\pi }}dw\sqrt{w^24m_\pi ^2}\{2C_{13}(w)`$ (22)
$`\times (m_\pi ^2+\mu ^2w^2)^2+{\displaystyle \frac{c_2}{3}}[w^64w^4(\mu ^2+2m_\pi ^2)+2(m_\pi ^2\mu ^2)^3`$
$`+w^2(5\mu ^4+6\mu ^2m_\pi ^2+9m_\pi ^2)+4m_\pi ^2(\mu ^4m_\pi ^4)(m_\pi ^2\mu ^2)w^2]\}.`$
This completes the presentation of analytical results. In Table 1, we present numerical results for the coordinate-space NN-potentials generated by the $`3\pi `$-exchange graphs of classes X, XII, XIII and XIV for inter-nucleon distances 0.8 fm$`r`$1.4 fm. We use the parameters $`f_\pi =92.4`$ MeV, $`m_\pi =138`$MeV (average pion mass) and the central value of the nucleon axial-vector coupling-constant $`g_A=1.267\pm 0.004`$ . For the second order low-energy constants $`c_{1,2,3,4}`$ we use the values $`c_1=0.8`$ GeV-1, $`c_2=3.2`$ GeV-1, $`c_3=4.7`$ GeV-1, $`c_4=3.4`$ GeV-1 mentioned above. One observes from Table 1 that the next-to-leading order $`3\pi `$-exchange NN-potentials are often larger than the leading order ones (comparing e.g. values at $`r=1.0`$ fm). The reason for this are the large values of the low-energy constants $`c_3`$ and $`c_4`$ representing the strong $`\mathrm{\Delta }`$-excitation effects. A similar feature has also been observed for the chiral $`2\pi `$-exchange NN-potential in ref.. In that case the effects occured however more selectively in the various channels. For example, there was no $`2\pi `$-exchange isoscalar central potential at leading order and at next-to-leading a large isoscalar central potential with a strength of about $`300`$ MeV at $`r=1.0`$ fm was generated mainly by the $`\pi \pi NN`$-contact vertex proportional to $`c_3`$. Note that the low-energy constant $`c_3`$ is related to the nucleon axial polarizability . The repulsive isovector central potential $`\stackrel{~}{W}_C^{(XIII)}(r)`$ from class XIII with a strength of almost $`50`$ MeV at $`r=1.0`$ fm is also fairly large. The largest potential is however the attractive isovector tensor potential $`\stackrel{~}{W}_T^{(XIII)}(r)`$. At distances $`r1.1`$ fm it actually overwhelms the repulsive isovector tensor potential of the (point-like) $`1\pi `$-exchange. (Since the eigen-values of $`\stackrel{}{\tau }_1\stackrel{}{\tau }_2`$ are $`3`$ and $`+1`$ the $`1\pi `$-exchange tensor potential can also act attractively, e.g. in the $`{}_{}{}^{3}S_{1}^{}`$-state relevant for deuteron binding.) It is well-known that phenomenology requires a reduction of the $`1\pi `$-exchange NN-tensor force at intermediate distances. The effect obtained here is presumably too strong and in addition the isovector tensor potential gets reduced by the chiral $`2\pi `$-exchange (see ref.). Because of the finite size (quark-substructure) of the nucleon one can trust the present calculation based on point-like chiral $`\pi N`$-interactions at most for distances $`rr_01`$ fm. Clearly, as soon as the nucleons start to overlap substantially the picture of an undisturbed multi-pion exchange between nucleons breaks down. The limiting distance scale $`r_04/m_\rho 1`$ fm (with $`m_\rho =770`$ MeV the $`\rho `$-meson mass) should be considered as a conservative estimate based on the argument that the nucleon axial radius $`r_A\sqrt{6}/m_\rho `$ receives no substantial contributions from chiral pion-loops . Only by confronting chiral multi-pion exchange potentials with NN-scattering data one is able to determine this limiting distance scale $`r_0`$ more precisely.

By analyzing the elastic proton-proton scattering data-base below $`350`$ MeV laboratory kinetic energy it has been shown recently in ref. that $`1\pi `$\- plus chiral $`2\pi `$-exchange gives a very good strong NN-force at least as far inwards as $`r=1.4`$ fm (i.e. the length scale corresponding to the pion Compton wave-length $`1/m_\pi `$). For such distances ($`r1.4`$ fm) the chiral $`3\pi `$-exchange NN-potential calculated here and previously in refs.) turns out to be a negligibly small correction. By inspection of Table 1 one observes that $`3\pi `$-exchange potentials typically decrease by about a factor of 10 each $`0.3`$ fm step outwards.

Fig.6: $`3\pi `$-exchange diagrams of class VIII proportional to $`g_A^6`$.

Finally, we like to present results for the leading order spin-spin and tensor NN-potentials generated by the $`3\pi `$-exchange diagrams proportional to $`g_A^6`$ (classes VIII and IX shown in Figs.6 and 7) which could not be evaluated with the help of the spectral-function method in ref.. The spin- and isospin factors of the two diagrams of class VIII (shown in Fig.6) are such that the irreducible part of the first graph and the second graph give rise only to a non-vanishing isoscalar central and isovector spin-spin and tensor potential. For the latter the spectral-function method lead to non-integrable boundary singularities in the $`3\pi `$-phase space integral. Therefore an alternative calculational technique is required.

The coordinate-space potential is generally given as the 3-dimensional Fourier transform $`d^3q\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r})\mathrm{}`$ of the 2-loop momentum-space amplitude constructed from the propagators and vertices of a diagram. This requires in total 11 integrations, $`d^3qd^4l_1d^4l_2`$. Usually, one performs first the integrations over the loop-energies $`l_{01}`$ and $`l_{02}`$ using residue calculus and one obtains this way the sum of energy denominators of all time-orderings of the diagram. In the present case this would lead to a complicated rational function of the three (on-shell) pion-energies $`\omega _{1,2,3}`$ which could not be further integrated analytically. The strategy is to perform a Wick-rotation of both loop-energies, $`l_{0j}il_{4j},j\{1,2\}`$, (of course one has to proof that this is indeed allowed) and to evaluate the remaining three 3-dimensional integrations analytically. The latter is possible because the original integrand factorizes in the variables $`\stackrel{}{l}_1`$, $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3=\stackrel{}{q}\stackrel{}{l}_1+\stackrel{}{l}_2`$ (see also the appendix in ref.). The result of this involves among other factors (originating from the vertices of the diagram) a product of three Yukawa-functions $`\mathrm{exp}(r\sqrt{m_\pi ^2+l_{4j}^2})/4\pi r`$. This way one derives the following double-integral representations for the isovector spin-spin and tensor potentials from class VIII,

$`\stackrel{~}{W}_S^{(VIII)}(r)`$ $`=`$ $`{\displaystyle \frac{7g_A^6}{6(16\pi f_\pi ^2)^3r^7}}\{e^{3x}(8+18x+16x^2+6x^3+2x^4x^5)`$ (23)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\zeta _1d\zeta _2}{\pi ^2\zeta _1\zeta _2}}{\displaystyle \frac{^2}{\zeta _1\zeta _2}}[e^{x\mathrm{\Sigma }}(12x^2+12x^1\mathrm{\Sigma }+2\mathrm{\Sigma }^2+8Q`$
$`+x(2\mathrm{\Sigma }Q+6\mathrm{\Pi })+2x^2\mathrm{\Sigma }\mathrm{\Pi }x^4\mathrm{\Pi }^2)]\},`$
$`\stackrel{~}{W}_T^{(VIII)}(r)`$ $`=`$ $`{\displaystyle \frac{7g_A^6}{6(16\pi f_\pi ^2)^3r^7}}\{e^{3x}(1+x)(16+26x+18x^2+6x^3+x^4)`$ (24)
$`+{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\zeta _1d\zeta _2}{\pi ^2\zeta _1\zeta _2}}{\displaystyle \frac{^2}{\zeta _1\zeta _2}}[e^{x\mathrm{\Sigma }}(36x^2+36x^1\mathrm{\Sigma }+10\mathrm{\Sigma }^2+16Q`$
$`+x(10\mathrm{\Sigma }Q+6\mathrm{\Pi })+x^2(4\mathrm{\Sigma }\mathrm{\Pi }+3Q^2)+3x^3Q\mathrm{\Pi }+x^4\mathrm{\Pi }^2)]\},`$
where we have used the dimensionless variable $`x=m_\pi r`$ together with various abbreviations:

$$\alpha =\sqrt{1+\zeta _1^2},\beta =\sqrt{1+\zeta _2^2},\gamma =\sqrt{1+(\zeta _1\zeta _2)^2},$$

$$\mathrm{\Sigma }=\alpha +\beta +\gamma ,Q=\alpha \beta +\beta \gamma +\gamma \alpha ,\mathrm{\Pi }=\alpha \beta \gamma .$$
(25)
The asymptotic behavior $`e^{3m_\pi r}r^2`$ of $`\stackrel{~}{W}_{S,T}^{(VIII)}(r)`$ in eqs.(23,24) indicates that the corresponding mass-spectra start with a non-vanishing value at the $`3\pi `$-threshold $`\mu =3m_\pi `$. For the numerical evaluation of the principal-value integrals in eqs.(23,24) it is advantageous to employ the identity $`_{\mathrm{}}^+\mathrm{}𝑑\zeta f(\zeta )\zeta ^1=_0^+\mathrm{}𝑑\zeta [f(\zeta )f(\zeta )]\zeta ^1`$, where the right hand side is singularity-free. As a check we have derived an analogous double-integral representation for the isoscalar central potential $`\stackrel{~}{V}_C^{(VIII)}(r)`$ and we have found perfect agreement with the numerical results of ref. based on the simpler spectral-function method.

Lastly, we come to the diagrams of class IX shown in Fig.7. For the isoscalar spin-spin and tensor NN-amplitudes the relevant sum of energy denominators is $`(\omega _1^2+\omega _2^2)/(2\omega _1^4\omega _2^4\omega _3^2)`$. Since its denominator factorizes in the three variables $`\stackrel{}{l}_1`$, $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3=\stackrel{}{q}\stackrel{}{l}_1+\stackrel{}{l}_2`$ one can perform the corresponding 9 integrations analytically and one gets the following closed form expressions,

$$\stackrel{~}{V}_S^{(IX)}(r)=\frac{2g_A^6}{(16\pi f_\pi ^2)^3}\frac{e^{3x}}{r^7}(2+2x+x^2)(22x4x^2x^3),$$
(26)
$$\stackrel{~}{V}_T^{(IX)}(r)=\frac{2g_A^6}{(16\pi f_\pi ^2)^3}\frac{e^{3x}}{r^7}(2+2x+x^2)(2+xx^2x^3).$$
(27)
We have corrected here in eq.(28) a small misprint (factor 2) which has occured in eq.(24) of ref.. The same features concerning factorization hold for the isovector central NN-amplitude and one reproduces the result $`\stackrel{~}{W}_C^{(IX)}(r)=(g_A^6/2048\pi ^3f_\pi ^6r^7)e^{3x}(1+x)(4+5x+3x^2)`$ of ref.. For the isovector spin-spin and tensor potentials generated by class IX one repeats (in a slightly modified form) the technique used previously for class VIII and one ends up with the following representations,

$`\stackrel{~}{W}_S^{(IX)}(r)`$ $`=`$ $`{\displaystyle \frac{g_A^6}{3(16\pi f_\pi ^2)^3r^7}}\{e^{3x}(2+2x+x^2)(22x4x^2x^3){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\zeta _1d\zeta _2}{\pi ^2\zeta _1\zeta _2}}`$ (28)
$`\times {\displaystyle \frac{^2}{\zeta _2^2}}[e^{x\mathrm{\Sigma }}(12x^2+12x^1\mathrm{\Sigma }+2\mathrm{\Sigma }^2+8Q+12\beta ^2+x((2\mathrm{\Sigma }+12\beta )Q6\mathrm{\Pi })`$
$`+x^2(2\mathrm{\Sigma }\mathrm{\Pi }+4\beta ^2(\alpha ^2+3\alpha \gamma +\gamma ^2))+4x^3(Q\alpha \gamma )\mathrm{\Pi }+x^4\mathrm{\Pi }^2)]\},`$
$`\stackrel{~}{W}_T^{(IX)}(r)`$ $`=`$ $`{\displaystyle \frac{g_A^6}{3(16\pi f_\pi ^2)^3r^7}}\{e^{3x}(2+2x+x^2)(2+xx^2x^3)`$ (29)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\zeta _1d\zeta _2}{\pi ^2\zeta _1\zeta _2}}{\displaystyle \frac{^2}{\zeta _2^2}}[e^{x\mathrm{\Sigma }}(2\mathrm{\Sigma }^24Q+x(2\mathrm{\Sigma }Q6\mathrm{\Pi })`$
$`+x^2(Q^2+2\alpha ^2\gamma ^2)+x^3(Q+2\alpha \gamma )\mathrm{\Pi }+x^4\mathrm{\Pi }^2)]\}.`$
Again, an analogous (double-integral) representation for the isoscalar central potential $`\stackrel{~}{V}_C^{(IX)}(r)`$ turned out to be in perfect numerical agreement with the result of ref. based on the spectral-function method.

Fig.7: $`3\pi `$-exchange diagrams of class IX proportional to $`g_A^6`$.

Numerical values for the spin-spin and tensor NN-potentials generated by the diagrams of classes VIII and IX are presented in Table 1. With a typical strength of a few MeV (or less) at $`r=1.0`$ fm they are comparable to the other leading order chiral $`3\pi `$-exchange potentials calculated in ref..

In summary, we have completed in this work the calculation of the leading order chiral $`3\pi `$-exchange NN-potentials. Furthermore, we investigated the dominant next-to-leading order contributions. The (analytical) results presented here and in refs. are parameterfree and in a form such that they could be easily implemented in a future empirical analysis of low-energy elastic NN-scattering. This way one would be able to determine the maximal range of validity $`rr_0`$ of chiral multi-pion exchange in the NN-interaction.


format: text

text:
# Monte Carlo determination of the phase diagram of the double-exchange model

## I Introduction.

The main magnetic properties of doped manganites and other materials are described in terms of double exchange interactions. An extensive body of work has been devoted to the understanding of this model, and the related ferromagnetic Kondo lattice, since the pioneering work of de Gennes (for a recent review, see).

While the model is extremely simple to formulate, as it can be written in terms of two dimensionless parameters only (see next section), the elucidation of its phase diagram is yet to be completed. Initial studies suggested that the phase diagram could be understood in terms of a ferromagnetic and an antiferromagnetic phase with an intermediate phase with a canted arrangement of moments. More recent work showed that in the region of the phase diagram where canting was expected, phase separation is more likely. By now, there are extensive studies of the phase diagram of the model in one dimension, two dimensions, and infinite dimensions. Even in these relatively simple cases, the model shows a rich phase diagram, which is not completely elucidated.

The situation is far from understood in the more relevant case of three dimensions. There is a general consensus that in the limit of infinite Hund’s coupling (see next section) the system has a ferromagnetic phase, although phase separation near the Curie temperature at low dopings has also been found. If the Hund’s coupling is not infinite, or if there are antiferromagnetic superexchange interactions present, the situation becomes more complicated. There is evidence that an antiferromagnetic coupling can induce a first order phase transition at relatively large dopings. Emphasis on the presence of first-order transitions at low temperatures has also been done in Refs. . In three dimensions, Berry’s phases in the wave function of the electrons may arise, leading to another source of complexity.

The present work presents a numerical study of the double exchange model in three dimensions, using a method discussed elsewhere. We show that the interplay between antiferromagnetic interactions and a three dimensional structure leads to a rich phase diagram, where, in addition to ferro- and antiferromagnetic phases, other ordered phases with cubic symmetry are possible. The model is discussed in the next section. Then, we present the calculation at zero temperature and the finite temperature results. In the last section, we discuss the main features of our findings, and their experimental implications.

## II The model.

### A The Hamiltonian.

We study the double exchange model. In its simplest form, electrons in a cubic lattice are coupled to localized spins defined at the same lattice sites. The spin of the electrons is constrained to be parallel to that of the localized spins. This restriction leads to a modulation of the electron hopping between lattice sites, which depends on the relative orientation of the core spins. In addition, one can define a direct coupling between spins at different sites. The Hamiltonian is:

$$=J_{\mathrm{A}F}\underset{ij}{}𝑺_i𝑺_j+t\underset{ij}{}(\theta _i\varphi _i|\theta _j\varphi _jc_i^{}c_j+\mathrm{h}.c.).$$
(1)
where $`\theta _i\varphi _i|\theta _j\varphi _j`$ is the overlap between the spinors defined by the polar angle $`\theta `$, and the azimuthal angle $`\varphi `$:

$$\theta _i\varphi _i|\theta _j\varphi _j=\mathrm{cos}\frac{\theta _i}{2}\mathrm{cos}\frac{\theta _j}{2}+\mathrm{sin}\frac{\theta _i}{2}\mathrm{sin}\frac{\theta _j}{2}e^{i(\varphi _i\varphi _j)}.$$
(2)
We describe the spins $`𝑺`$ by classical variables, normalized to one, $`|𝑺|=1`$. The system depends on two dimensionless parameters only, the filling of the electron band, $`x`$, and the ratio $`J_{\mathrm{A}F}/t`$ (notice that our $`J_{\mathrm{A}F}`$ is opposite in sign to the convention of )

The double exchange Hamiltonian, Eq. (1), can be derived from the model which describes a lattice of atoms with strong intra-atomic Hund’s coupling between electrons in different orbitals. The simplest Hamiltonian which includes this effect is that of the ferromagnetic Kondo lattice,

$$_{\mathrm{F}K}=t\underset{i,j,s}{}(c_{i,s}^{}c_{j,s}+\mathrm{h}.c.)J_\mathrm{H}\underset{i,s,s^{}}{}c_{i,s}^{}c_{i,s^{}}𝝈_{ss^{}}𝑺_i,$$
(3)
where $`J_\mathrm{H}>0`$, and the $`\sigma `$’s are Pauli matrices. When $`J_\mathrm{H}t`$, the spin sub-band with an antiparallel orientation to the core spin lies at high energy, and can be projected out. To first order in $`1/J_\mathrm{H}`$, one finds an antiferromagnetic interaction between nearest neighbor core spins, leading to the double exchange Hamiltonian, Eq. (1), with $`|J_{\mathrm{A}F}|=t^2/(4J_\mathrm{H})`$ (note that we have normalized the core spins). This effect adds to an antiferromagnetic superexchange interaction between the Mn $`𝐭_𝐠`$ core spins. Therefore, the most phenomenologically interesting range for $`|J_{\mathrm{A}F}|/t`$ is $`0.008`$$`0.15`$ .

### B Other interactions.

The model discussed above is supposed to be a reasonable starting point to the study of the magnetic properties of doped manganites, La$`{}_{1x}{}^{}A_{x}^{}`$MnO3, where $`A`$ stands for a divalent cation. The electrons of the conduction band arise from the Mn $`𝐞_𝐠`$ orbitals, while the core spin is built up from three electrons in the Mn $`𝐭_𝐠`$ orbitals. In a cubic lattice, the $`𝐞_𝐠`$ orbitals are doubly degenerate. This degeneracy is ignored in Eq. (1). Thus, we cannot study effects associated to orbital ordering. These effects are probably important in explaining the different phases with charge ordering observed in these materials. The model studied here can be considered a first approximation to the study of the magnetic properties, in situations where there is no charge ordering and the lattice symmetry is cubic.

The existence of double degeneracy in the $`𝐞_𝐠`$ orbitals implies the possibility of lattice Jahn-Teller distortions, which have been considered the cause of the breakdown of cubic symmetry at low dopings. In addition, dynamic Jahn-Teller fluctuations can play a role in some properties of doped compounds with cubic symmetry. We will ignore these interactions. Previous studies suggest that the model used here suffices to understand the main features of the magnetic transitions in the manganites, in the doping range $`0.1x0.4`$ where most experiments are done.

The treatment of the core spins as classical variables will overestimate the tendency of the system toward long range order. Previous studies suggest that this effect is small, as results obtained for the physical value $`S=3/2`$ are very close to those found in the limit $`S\mathrm{}`$. In addition, the only transitions at zero temperature found in our work are discontinuous (see next section). Hence, quantum spin fluctuations are bounded, and cannot change qualitatively the results.

## III Results.

### A Phase diagram at zero temperature.

The phase diagram at zero temperature can be analyzed by minimizing the total energy, Eq. (1), as function of the spin configuration, for the different values of $`x`$ and the ratio $`J_{\mathrm{A}F}/t`$. However, previous studies in one, two and infinite dimensions  have shown that the system has a strong tendency toward phase-separation. This enforces us to use the Maxwell construction to correctly find the phase boundaries. The easiest way of performing the Maxwell construction consist in minimizing the grand-canonical energy

$$^{\mathrm{G}C}=\mu \underset{i}{}c_i^{}c_i,$$
(4)
as a function of the spin configuration for the different values of $`\mu `$ and $`J_{\mathrm{A}F}`$. In this way, we obtain the band filling, $`x`$, as a function of the chemical potential $`\mu /t`$ and of $`J_{\mathrm{A}F}`$. At the phase boundaries $`x`$ is a discontinuous (growing) function of $`\mu `$: the jump at the discontinuity is the phase-separation compositional range.

For the above mentioned minimization we will limit ourselves to the possible spin textures observed in Monte Carlo calculations at finite temperatures (see next subsection) and to those which seem plausible on physical grounds. Some of the most stable phases analyzed are described in Table I. We have studied only the range $`x0.5`$, due to the particle-hole symmetry of the Hamiltonian (1).

We find the ferromagnetic (FM) and antiferromagnetic (G-AFM) phases. The A-AFM configuration is ferromagnetic within planes, and antiferromagnetic between neighboring planes. The C-AFM configuration is ferromagnetic along lines, and antiferromagnetic between neighboring lines. These phases have already been discussed in the literature. The twisted phase interpolates between the A-AFM and the C-AFM phases. Island phases have been reported in calculations in one and two dimensional systems. These are structures of spins aligned in one direction with the sense varying periodically. For the values of $`J_{\mathrm{A}F}/t`$ considered here, we only find the $`(\frac{\pi }{2},\frac{\pi }{2},\pi )`$ (that corresponds $`2\times 2\times 1`$ blocks of equal spins) and $`(\frac{\pi }{3},\pi ,\pi )`$ ($`3\times 1\times 1`$ blocks), although we have tried different combinations. In these phases, the electrons are localized, and the electronic density of states is built up of delta functions. Helix phases are ordinary spin density waves. It is interesting to note that helix phases are completely irrelevant in the thermodynamic limit, but they are remarkably stable on finite lattices, and should be considered when analyzing numerical simulations.

We find two new phases, intrinsic to the cubic lattice, labeled flux and skyrmion in Table I. The spin configuration is not coplanar, and the spins are parallel to the diagonals of the unit cube (note that the energy is invariant under a global rotation, so that only the relative angles are relevant). In the skyrmion phase, the spin directions around a unit cube point out from the center, with a hedgehog shape. The induced Berry phase can be thought of as generated by an effective magnetic monopole inside each of the unit cubes which build up the three dimensional lattice. The charge of the monopoles change sign in neighboring cubes, leading to zero total charge. The flux within each cube is isotropic, and equal to $`\pi /3`$ (skyrmion configuration) and $`2\pi /3`$ (flux configuration) per plaquette. In both cases, the lattice unit contains eight sites, and the global symmetry remains cubic.

The electronic dispersion relation for all these phases can be calculated analytically, and it is given in Table I. The corresponding density of states is shown in Fig. 1 . The dispersion relation in the flux and skyrmion phases depend linearly on $`𝒌`$ near $`𝒌=(\pm \frac{\pi }{2},\pm \frac{\pi }{2},\pm \frac{\pi }{2})`$, and resembles the dispersion of Dirac massless particles in the lattice.

The flux and twisted phases show similarities to the phase with spins at right angles found in two dimensions.

In the twisted phase, for $`a=b=1/\sqrt{2}`$, the planes $`xz`$ and $`yz`$ have the same structure as the 2D pattern studied in. The flux phase, on the other hand, has a dispersion relation which is the natural extension to 3D of the one for the 2D flux phase. Note, however, that the flux phase in 2D is coplanar, and all phases in the hopping elements can be rendered real by an appropriate gauge transformation, while this does not happen in the 3D phase studied here.

The stability of the flux and skyrmion phases arise from the canting of the spin orientations, coupled to a shift in the electronic density of states toward the band edges. It is interesting to note that the electronic energy of the flux phase is close to that of the C-AFM one, and the energy of the skyrmion phase is also close to that of the A-AFM phase, for the entire range of electronic concentrations.

The phase diagram at zero temperature is shown in Fig. 2. All transitions are first order. The most interesting result is the stabilization of the skyrmion phase for a range of dopings and values of $`J_{\mathrm{A}F}/t`$ where the model is applicable to real materials. The twisted phase changes continuously from $`a^2=0.45`$ (right) to $`a^2=0.77`$ (left). Notice that although the A-AFM phase is beside the twisted phase, there is a discontinuity in the $`a`$ parameter ($`a=0`$ for A-AFM). We find island phases in a more restricted range of parameters than in two dimensions, possibly due to the competition with the flux and skyrmion phases. We have not found any tridimensional island as $`(\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2})`$. For more negative values of $`J_{\mathrm{A}F}/t`$ we also find a $`(\frac{\pi }{2},\pi ,\pi )`$ island phase starting at $`0.19`$.

\[    Type spin direction $`ϵ(𝒌)/t`$    Ferromagnetic (FM) $`(0,0,1)`$ $`2_{\mu =1}^3\mathrm{cos}k_\mu `$  A-type (A-AFM) $`(0,0,(1)^z)`$ $`2_{\mu =1}^2\mathrm{cos}k_\mu `$  C-type (C-AFM) $`(0,0,(1)^{x+y})`$ $`2\mathrm{cos}k_3`$  G-type (G-AFM) $`(0,0,(1)^{x+y+z})`$ 0  Twisted $`(a(1)^{x+y},b(1)^z,0)`$ $`\pm 2\sqrt{a^2\mathrm{cos}^2k_3+b^2(\mathrm{cos}k_1\pm \mathrm{cos}k_2)^2}`$   $`a^2+b^2=1`$   Flux $`((1)^{y+z},(1)^{x+z},(1)^{x+y})/\sqrt{3}`$ $`\pm 2\sqrt{\frac{1}{3}_{\mu =1}^3\mathrm{cos}^2k_\mu }`$  Skyrmion $`((1)^x,(1)^y,(1)^z)/\sqrt{3}`$ $`\pm 2\sqrt{\frac{2}{3}}\sqrt{_{\mu =1}^3\mathrm{cos}^2k_\mu \pm \sqrt{3_{\mu \nu }\mathrm{cos}^2k_\mu \mathrm{cos}^2k_\nu }}`$  Helix $`(\mathrm{cos}(\omega z),\mathrm{sin}(\omega z),0)`$ $`2_{\mu =1}^2\mathrm{cos}k_\mu 2\mathrm{cos}(\omega /2)\mathrm{cos}(k_3+\omega /2)`$  Island $`(\frac{\pi }{2},\pi ,\pi )`$ $`(0,0,(1)^{[\frac{x}{2}]+y+z})`$ $`1,1`$  Island $`(\frac{\pi }{3},\pi ,\pi )`$ $`(0,0,(1)^{[\frac{x}{3}]+y+z})`$ $`\sqrt{2},0,\sqrt{2}`$  Island $`(\frac{\pi }{2},\frac{\pi }{2},\pi )`$ $`(0,0,(1)^{[\frac{x}{2}]+[\frac{y}{2}]+z})`$ $`2,0,0,2`$

Table I. Spin configurations and electronic dispersion relations of the different phases considered in the text. The notation $`[]`$ stands for the integer part. For a lattice of linear dimension $`L`$, $`\omega `$ and $`k_\mu `$ can be written as $`\frac{2\pi }{L}n`$, with $`n`$ integer (in twisted, flux and skyrmion configurations $`n`$ is even as the unitary cell is a $`2^3`$ cube).

\]

### B Finite temperature results.

We have extended the previous studies to finite temperatures by using the hybrid Monte Carlo method reported elsewhere , which allows to thermalize $`16\times 16\times 16`$ clusters in reasonable computer time. We have studied the model (1), on the cubic lattice of side $`L`$ with periodic boundary conditions, using the grand-canonical ensemble (4). The use of an efficient algorithm has been crucial, because around two hundred points of the $`(T,\mu ,J_{\mathrm{A}F}/t)`$ phase-diagram have been studied.

Before presenting our results it is worth to comment on the strong finite-size effects of this model. Fortunately, one can have a good analytical command on them (at least at zero temperature), by repeating the minimization of subsection III A on a finite lattice (this is straightforward using Table I). In this way one discovers that the phase boundaries change considerably with the lattice size. Even worse, for lattice sizes $`L=4`$ and $`8`$, the helix phase of minimal frequency $`\omega =2\pi /L`$ (see Table I) is more stable than the ferromagnetic phase in important regions of the phase diagram. Also the C-AFM phase is stable in the $`L=8`$, although it is unstable in the $`L\mathrm{}`$ limit. On the other hand, the zero temperature phase diagrams of the $`L=6`$ and $`L=10`$ lattices are rather close to the infinite volume limit. The only caveat of this lattice sizes is that they cannot accommodate island-phase whose unit cell is not commensurate with them (in particular, the $`(\frac{\pi }{2},\frac{\pi }{2},\pi )`$ island phase, that has a $`4\times 4\times 2`$ unit cell). We have therefore chosen to perform most of the numerical simulations on the $`6^3`$ cluster and, only when necessary for elucidation, we have employed as a larger cluster the $`10^3`$. Exception to this rule has been the simulation of the $`8^3`$ in the range $`0.29J_{\mathrm{A}F}/t0.18`$ range, to check for the presence of the $`2\times 2\times 1`$ island phase.

In the simulation, we have used the perfect action formulation , with $`\lambda =0.25`$. The molecular dynamics parameters have been chosen to ensure an acceptance rate greater than an 80%. The selection depends mainly on the temporal lattice size (inverse temperature). Typical values used are $`N_\tau =20`$, $`\mathrm{\Delta }\tau =0.02`$ for $`T/t=1/20`$, and $`N_\tau =50`$, $`\mathrm{\Delta }\tau =0.01`$ for $`T/t=1/45`$. We have systematically controlled the correct thermalization of our data, by comparing the results of hot and cold starts. In some cases, in particular at $`T=t/100`$, we have also compared the results of ordered starts with a slow annealing from high temperatures. The total run-length (up to 2000 trajectories in some points) has always been at least larger than three times the thermalization time. The typical size of our statistical errors can be found in figure 3. The total CPU time devoted to this project has been the equivalent of 10 years of a Pentium III at 800 MHz.

As we have said in the introduction, our scope has been the determination of the magnetic phase diagram of the model. Since one cannot break a symmetry on a finite lattice, it is necessary to use pseudo order-parameters for obtaining the phase diagram . We have defined our pseudo order-parameters in terms of the spin structure factor ($`V=L^3`$ is the total number of spins),

$`G(𝒌)`$ $`=`$ $`|\widehat{𝑺}(𝒌)|^2,`$ (5)
$`\widehat{𝑺}(𝒌)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{V}}}{\displaystyle \underset{j}{}}𝑺_j\mathrm{e}^{\mathrm{i}𝒌𝒓_j}.`$ (6)
In the above equation, $``$ stands for the thermal average. In a paramagnetic phase, $`G(𝒌)`$ will be of order $`1/V`$ for all $`𝒌`$. On the other hand, on an ordered phase it will be of order one for a very small set of $`𝒌`$’s (see below), and of order $`1/V`$ for the rest of wave numbers. If one now calculates the structure factor for the configurations in Table I, it follows that the appropriate pseudo order-parameters are:

$`G_\mathrm{f}`$ $`=`$ $`G(0,0,0),`$ (7)
$`G_\mathrm{s}`$ $`=`$ $`G(\pi ,\pi ,\pi ),`$ (8)
$`G_\mathrm{l}`$ $`=`$ $`G(\pi ,\pi ,0)+G(\pi ,0,\pi )+G(0,\pi ,\pi ),`$ (9)
$`G_\mathrm{p}`$ $`=`$ $`G(\pi ,0,0)+G(0,\pi ,0)+G(0,0,\pi ),`$ (10)
$`G_\mathrm{h}`$ $`=`$ $`G(\frac{2\pi }{L},0,0)+G(0,\frac{2\pi }{L},0)+G(0,0,\frac{2\pi }{L}).`$ (11)
Let us explain them in details. In a ferromagnetic phase, only $`G_\mathrm{f}`$ is of order 1, while for a G-AFM phase, the appropriate pseudo order-parameter would be $`G_\mathrm{s}`$. For an A-AFM phase, the selected order parameter will be $`G_\mathrm{p}`$ (we sum for $`𝒌=(\pi ,0,0),(0,\pi ,0)`$ and $`(0,0,\pi )`$ because we do not know a priori the direction of the AFM lattice planes). A skyrmion phase will also have $`G_\mathrm{p}`$ as pseudo order-parameter, and so the question arise of how to differentiate a skyrmion from an A-AFM phase. The easiest way to do that is to control that for the skyrmion phase the three quantities $`G(\pi ,0,0)`$, $`G(0,\pi ,0)`$ and $`G(0,0,\pi )`$ are simultaneously large (we have also numerically checked the orthogonality of the Fourier transformed spin field at the three $`𝒌`$ values). Also the C-AFM phase and the flux phase have a common pseudo order-parameter, namely $`G_\mathrm{l}`$, but one can differentiate both phases in exactly the same way. A twisted phase will have as pseudo order-parameter both $`G_\mathrm{p}`$ and $`G_\mathrm{l}`$. One can distinguish the twisted phase from a tunneling phenomena between an A-AFM and a C-AFM phases, checking again that both pseudo order-parameters are simultaneously large (we have also checked the orthogonality of the Fourier transformed spin-field at the two $`𝒌`$ for which the structure factor is not negligible). Finally the helix phase is signaled by the $`G_\mathrm{h}`$ order parameter, while the structure factor for the $`(k_1,k_2,k_3)`$ island phase is heavily peaked at $`𝒌=(k_1,k_2,k_3)`$.

The results for the $`G_\alpha `$ pseudo order-parameters and for the equilibrium fermionic density at $`T/t=1/45`$ are shown in Fig. 3, for $`\mu =0`$ (half filling, $`x=0.5`$) and $`\mu =1.5t`$ which includes the interesting region of $`x1/4`$$`1/3`$. The numerical results shown in this plot allow us to identify the FM, G-AFM, twisted and flux phases, and the onset of the skyrmionic phase.

For $`\mu =0`$, the flux-twisted transition is clearly discontinuous. In addition, there is a clear metastable behavior (the duplicated points correspond respectively to two metastable phases). At the right end of the twisted phase, $`G_\mathrm{p}`$ and $`G_\mathrm{l}`$ vanish abruptly and there are clear signs of metastability with a phase where $`G_\mathrm{f}`$ is small but non-vanishing. To check if it is a finite size effect, we have performed some simulations in this region with a larger lattice ($`10^3`$), obtaining that $`G_\mathrm{f}`$ is much smaller (and $`G_\mathrm{p}`$ and $`G_\mathrm{l}`$ vanish). Thus we conclude that there is a first order transition at this temperature between the twisted and the paramagnetic phases. On the contrary, the ferromagnetic-paramagnetic transition is smoother. We recall that in Ref.  it is shown that at $`J_{\mathrm{A}F}=0`$, $`T/t=0.14`$ the transition is second order. Although we find a strengthening of this transition when $`|J_{\mathrm{A}F}|`$ grows, to confirm that it actually becomes of the first-order would require a careful finite-size scaling analysis that is beyond the scope of this work.

For $`\mu =1.5`$ (see Fig. 3) we find a very different landscape. At the left part we find an almost saturated G-AFM region (notice that $`x0`$ and the fermions essentially play no role). In the center, we find a region with $`x0.25`$ where $`G_\mathrm{p}`$ is about 0.05 and the others are much smaller. We interpret this result as the onset of the skyrmionic region (that appears clearly at lower temperatures, see Fig. 5). The start of the ferromagnetic region is much abrupter than at $`\mu =0`$, but as said above, to conclude that it is first order, a study in different lattice sizes is mandatory.

To explore the phase diagram, we have simulated at several values of $`T`$. The phase diagram at half filling is plotted in Fig. 4. The paramagnetic region separates the ferromagnetic from the twisted phase for $`T/t1/45`$ but it is difficult to know the point where it finishes (maybe $`T=0`$). As the ferromagnetic magnetization goes to zero rather smoothly, to obtain a more precise value of the transition point we have used data from larger lattices (mainly $`L=10`$). The dashed lines have been draw approximately at the phase transitions.

Results for the different $`G_\alpha `$ order parameter for $`\mu =1.5t`$ are shown in Fig. 5. Skyrmion correlations are suppressed at relatively low temperatures, but for $`T=t/100`$ we find a clear skyrmionic structure ($`G_\mathrm{p}=0.7`$$`0.8`$, with the appropriate orthogonality of the Fourier components).

We have combined all our results at $`T=t/45`$ in the phase diagram shown in Fig. 6. The results are consistent with the zero temperature calculation for the $`6^3`$ lattice, whose transitions are plotted as continuous lines. The FM, G-AFM and Flux phases (and the twisted one near half filing) are less affected by finite temperature effects than the other phases, that become paramagnetic for temperatures as small as $`T=t/45`$. Notice also that in the $`6^3`$ lattice the island phases are so narrow that they are hardly observed in a simulation. However, we have checked that in a $`8^3`$ lattice at $`\mu =0`$, $`T=t/45`$ the system is in the $`2\times 1\times 1`$ phase for $`J_{\mathrm{A}F}/t(0.29,0.18)`$. We have also simulated a $`8^3`$ lattice in several values of $`(\mu ,J_{\mathrm{A}F})`$ that at $`T=0`$ correspond to the $`2\times 2\times 1`$ region, observing that even at $`T=t/45`$ the system is disordered.

### C The skyrmion phase.

As mentioned above, we find a novel phase, the skyrmion phase, for a realistic range of dopings and antiferromagnetic couplings. We now discuss some physical features of this phase. Similar properties can be expected in the flux phase.

i) The symmetry of the magnetic phase is cubic. This phase is compatible with experiments done in the range $`x1/4`$$`1/3`$ which show phases with cubic structure but no macroscopic magnetization.

ii) The magnetic arrangement is, at first sight, similar to that expected in a canted phase. Measurements of the microscopic magnetic structure will probably be indistinguishable from those expected for a canted phase.

iii) The magnon spectrum will resemble that of other systems with non collinear spins. The dispersion relation will be linear at low energies.

iv) The phase is metallic. However, the chemical potential, at the appropriate dopings, coincides roughly with the position when the density of states has a pronounced minimum (see Fig. 1 ). Hence, it cannot be ruled out that residual disorder will induce insulating behavior. If this is the case, it will be behave in a similar way to the canted insulating phases reported in the literature (see, for instance ).

v) The unusual spin configuration leads to an effective coupling between an external magnetic field and a charge density wave. An applied field tilts the spins. The sign of the change in the angle between neighboring spins depends on the relative orientation of the spins and the magnetic field. A reduction in the angle leads to an increase of the hopping, and the opposite happens if the angle is increased. Hence, electron charge will be displaced from the weak bonds to the strong ones, and a charge density wave will arise. This feature can clearly identify the skyrmion phase.

## IV Conclusions.

We have analyzed the phase diagram of the double exchange model in a cubic lattice, both at zero and at finite temperature. Our results at finite temperature have been obtained using an efficient Hybrid Monte Carlo algorithm , which has been crucial for the exploration of a three dimensional phase diagram. We have found that the competition between the ferromagnetic interaction mediated by the conduction electrons and antiferromagnetic couplings which arise from superexchange effects or from the underlying intra-atomic Hund’s coupling lead to the existence of a variety of phases.

We have expressed our results in terms of the electron hopping, $`t=W/12`$, where $`W2`$ eV is the width of the conduction band, and the effective antiferromagnetic coupling, $`J_{\mathrm{A}F}`$, arising from superexchange interactions between the Mn $`𝐭_𝐠`$ orbitals ($`|J_{\mathrm{A}F}|/t0.005`$$`0.012`$ ) and from the finite value of the intra-atomic Hund coupling, $`J_\mathrm{H}`$. It follows that the most phenomenologically interesting range for $`|J_{\mathrm{A}F}|/t`$ is $`0.008`$$`0.15`$ . Using the previous values, we find that the highest transition temperature ($`x=0.5,J_{\mathrm{A}F}=0`$) is close to $`350`$ K, and therefore most of the phase transitions reported in this work occur in the $`0`$$`300`$K range, in agreement with experiments.

It is somewhat surprising that conventional spin density wave structures are not stable. The three dimensional phase diagram shares some features with solutions of the double exchange model in one, two and infinite dimensions, like phase separation. A number of phases, however, have no counterparts in other dimensions.

We have identified, among others, a phase where the spins are locally arranged in a hedgehog manner, labelled skyrmion phase. The existence of this structure is probably intrinsic to the double exchange mechanism, and cannot be realized in spin systems with short range couplings. The skyrmion phase is stable for dopings $`x0.2`$$`0.3`$ and antiferromagnetic couplings $`|J_{\mathrm{A}F}|/t0.05`$$`0.07`$. These values are realistic for doped manganites in the range where Colossal Magnetoresistance effects are found. While it is probably difficult to characterize this phase experimentally, it will be interesting to verify its existence.

## V Acknowledgements.

We want to warmly thank Victor Laliena for a pleasant and fruitful collaboration. It is also a pleasure to thank B. Alascio, A. Aligia, E. Dagotto and B. Normand, for discussions on this problem.

We acknowledge financial support from grants PB96-0875, AEN97-1680, AEN97-1693, AEN99-0990 (MEC, Spain) and (07N/0045/98) (C. Madrid). V.M.-M. is supported by the Ministerio de Educación y Cultura. The simulations have been carried out in RTNN computers at Zaragoza and Madrid.