text:
# Maximally–localized Wannier Functions in Antiferromagnetic MnO within the FLAPW Formalism

## I INTRODUCTION

The mean–field one–particle description of the electronic struture of periodic crystalline solids is usually based on extended Bloch functions (BFs). Within the Born–von Karman periodic boundary conditions, the cyclic translational subgroup, which commutes with the effective one–electron Hamiltonian, contains $`N`$ translations by corresponding direct lattice vectors $`𝐑`$. This Abelian subgroup has $`N`$ one–dimensional irreducible representations, which are labeled with wavevectors $`𝐤`$ within the first Brillouin zone (BZ). Therefore, extended Bloch states $`\psi _{m𝐤}(𝐫)`$ in this description are classified with two quantum numbers, the band index $`m`$ and the crystal–momentum $`𝐤`$, and are obtained by diagonalization of the effective one–electron Hamiltonian. An alternative description can be derived in terms of localized Wannier functions Wannier37  (WFs) $`w_n(𝐫𝐑)`$, which are defined in real space via a unitary transformation performed on the Bloch functions. They are also labeled with two quantum numbers, the orbital index $`n`$ and the direct lattice vector $`𝐑`$ indicating the unit cell they belong to. In contrast to BFs, WFs are useful, for example, in visualizing chemical bonds or in describing the dielectric properties of non–metallic materials. They can be considered as the generalization to solids of the concept of “localized molecular orbitals” for finite systems Boys96 . However, one major problem in practical calculations within this representation is the non–uniqueness of WFs, related to the phase arbitrariness of the BFs, and to the arbitrary unitary transformations which can be performed on the BFs at any given $`𝐤`$–point in the BZ. As demonstrated recently by Marzari and Vanderbilt Marzari97 , this non–uniqueness can be resolved in principle by imposing the further condition of maximal localization. However, some residual arbitrariness remains, related to the choice of the localization criterion.

The late transition–metal (TM) monoxides MnO, FeO, CoO, and NiO challenge the theory of electronic states since several decades. These highly correlated materials feature a Mott insulator character, and the conventional local–spin–density (LSD) scheme Terakura84  gets into difficulties with the localized TM $`d`$–orbitals, predicting in particular incorrectly their spectral weight and their energy relative to O $`2p`$ states. All these oxides have the rocksalt structure in their paramagnetic phase, while below the Neel temperature $`T_\text{N}`$, a type II antiferromagnetic (AFM) ordering occurs Roth58 . The experimentally well–documented compound MnO, which is considered to be in the intermediate charge–transfer/Mott–Hubbard regime, is a particularly suitable case study, since the available theoretical schemes apply best to this material. Indeed, because of the exchange stabilization of its half–filled $`d`$ shell, even LSD predicts its insulating character, although the corresponding energy gap and magnetic moment are much smaller than experimental data. Several electronic structure calculations have been performed for MnO, using traditional band–structure schemes, like unrestricted Hartree–Fock (UHF) Massidda99  and local–spin density (LSD), as well as more innovative approaches taking into account at various levels the large value of the on–site Coulomb repulsion for the metal 3$`d`$ states, like the SIC Svane90  and the LSD + $`U`$ Anisimov91  methods, and the more recent model $`GW`$ scheme Massidda95 ; Massidda97 . Among the quantities calculated within the latter approach, special attention has been devoted to the quasi–particle spectrum Massidda95 , and recently to the zone–center optic phonon frequencies and the Born effective charge tensor Massidda99 . Comparison with experimental data of the various physical quantities computed in the different schemes, has demonstrated a monotonic trend with the separation energy between occupied and empty $`d`$ bands, which is too small within LSD, too high within UHF, and has about the correct value within model $`GW`$. This energy separation is strongly related to the on–site interaction $`U`$. The analysis of all these results, has been performed in terms of the extended Bloch states.

On the other hand, no ab–initio investigation exists of the electronic states in AFM MnO in terms of localized Wannier functions. Keeping in mind the inherent limitations of the mean–field one–particle approaches, it would be however instructive to materialize within the WFs description how superexchange (which is responsible for AFM ordering) manifests itself in a first–principles one–electron picture. Furthermore, the maximal localization method of Marzari and Vanderbilt 111A generalization of the method of Ref. Marzari97,  for the case of entangled energy bands has been given by I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001). has been applied in the past to several periodic systems Marzari97 ; Marzari98 ; Souza00 , but not yet to a low (trigonal) symmetry solid like AFM MnO, nor to a compound with partially filled cation $`d`$ shell.

Motivated by the above reasons, we present in this work an ab–initio calculation of the maximally–localized Wannier functions of AFM MnO corresponding to the uppermost occupied bands, using the all–electron, full–potential linearized–augmented–plane–wave (FLAPW) method Jansen84 ; Massidda93 . Because the trends in the electronic properties from UHF to LSD are already known for this system, most of the computations have been performed for convenience within the latter scheme. To investigate the effects related to the on–site Coulomb interaction, we have also used the LSD + $`U`$ method, implemented with the same FLAPW technical ingredients Shick99  used in the LSD computations.

The manuscript is organized as follows. In Sec. II we give technical details of the FLAPW implementation, and structural informations relevant to AFM MnO. The necessary ingredients of our WFs calculations are then introduced. As we follow closely the formalism of Ref. Marzari97, , we will not describe in detail the method itself. A general result regarding maximal localization of $`s`$$`p`$ WFs is deferred to an Appendix. In Sec. III, we present and discuss the results of our maximally–localized WFs calculations. Finally, in Sec. IV we draw our conclusions.

## II METHOD

### A. FLAPW calculations

All quantities presented and discussed in this paper have been computed using the semirelativistic FLAPW method Jansen84 ; Massidda93 , expanded with local orbitals Singh91  (lo) where appropriate. Inclusion of lo’s in addition to the normal FLAPW basis enforces mutual state orthogonality and increases variational freedom. This allows to treat the semicore Mn $`3s`$, $`3p`$ states together with the valence states, and helps in dealing with the linearization of Mn $`3d`$ and O $`2s`$, $`2p`$ states. Core states are calculated fully relativistically and self–consistently in the crystal potential. For the LSD calculations, the Hedin and Lundqvist exchange–correlation functional has been adopted. In the calculations performed within the LSD + $`U`$ method, the values of the Hubbard and exchange constants, $`U=6.9`$ eV and $`J=0.86`$ eV, have been taken from Ref. Anisimov91, . The atomic–sphere radii for Mn and O are chosen to be 2.0 and 1.8 a.u., respectively, and the FLAPW basis size is set to include all plane waves with energy up to 16.0 Ry. Four special points inside the irreducible wedge of the BZ were used for evaluating the charge density during self–consistency cycles.

### B. Structural details

AFM ordering occurs in MnO below $`T_\text{N}=117`$ K along any one of the $`111`$ directions. The transition is accompanied by a small crystallographic distortion that transforms the cubic structure into a rhombohedral one Landolt84 : the $`90^{}`$ angle between the lattice vectors increases only by $`0.624^{}`$. We have shown Massidda99  earlier that the effect of this distortion on the calculated zone–center optic phonon frequencies and Born effective charge tensor is negligible compared to the anisotropy induced solely by the magnetic order. Therefore, throughout the present study we use the ionic positions of the perfect rocksalt chemical cell, with the experimental lattice constant value $`a=4.435`$ Å. In the AFM phase, the magnetic moments of MnO are arranged into ferromagnetic sheets which are parallel to (111), while the direction of magnetization in neighboring planes is reversed. The magnetic and crystalline structure of MnO is shown in Fig. 1. The rhombohedral magnetic cell, whose volume is twice that of the paramagnetic rocksalt one, corresponds to the space group D$`{}_{3\text{d}}{}^{5}(R\overline{3}m)`$. Several choices for the primitive translation vectors are possible. However, in order to minimize discrete mesh effects in $`𝐤`$–space integrations, and also for taking advantage of symmetry properties (see below), we consider here the following primitive translation vectors (in the cubic coordinate system): $`𝐭_1=a(1,\frac{1}{2},\frac{1}{2})`$, $`𝐭_2=a(\frac{1}{2},\mathrm{\hspace{0.25em}1},\frac{1}{2})`$, and $`𝐭_3=a(\frac{1}{2},\frac{1}{2},\mathrm{\hspace{0.25em}1})`$. The two equivalent anions O1 and O2, whose site symmetry is C$`_{3\text{v}}`$, are located at positions $`\pm a(\frac{1}{2},\frac{1}{2},\frac{1}{2})`$, while the two non–equivalent cations Mn1 (spin–up) and Mn2, whose site symmetry is D$`_{3\text{d}}`$, are at $`(0,\mathrm{\hspace{0.25em}0},\mathrm{\hspace{0.25em}0})`$ and $`a(1,\mathrm{\hspace{0.25em}1},\mathrm{\hspace{0.25em}1})`$. With our choice of origin on Mn1, the two Oxygen sites are therefore related by spatial inversion.

### C. Maximally localized Wannier functions

We consider the general case of a composite group of energy bands labeled by their band index $`m=1,2,\mathrm{},f`$, and connected among themselves by degeneracies along high symmetry lines, but isolated from all other bands at lower and higher energy. The phase arbitrariness (gauge dependence) of each $`|\psi _{m𝐤}`$ with respect to the set of allowed $`𝐤`$ values propagates to the corresponding Wannier functions $`|w_{n𝐑}`$, labeled by $`n=1,2,\mathrm{},f`$ and lattice vectors $`𝐑`$. This resulting non–uniqueness manifests itself trough an arbitrary unitary matrix $`U_{mn}^{(𝐤)}`$ appearing in the most general transformation between BFs and WFs, which is given by

$$|w_{n𝐑}=\frac{V}{(2\pi )^3}_{BZ}𝑑𝐤e^{i𝐤𝐑}\underset{m}{}U_{mn}^{(𝐤)}|\psi _{m𝐤},$$
(1)
where $`V`$ is the real–space primitive cell volume. A direct consequence of this arbitrariness is that linear combinations of the WFs in the set $`\{|w_{n𝐑}\}`$ also form a suitable basis.

The strategy of Ref. Marzari97,  is to pick out from this arbitrary choice of WFs, the particular set which is maximally localized according to some criterion. Once this criterion has been chosen, and the composite group of bands specified, the search for the set of maximally localized WFs becomes a problem of functional minimization in the space of the matrices $`U_{mn}^{(𝐤)}`$. The selected functional, which measures the sum of the quadratic spreads of the WFs probability density, is given by

$$\mathrm{\Omega }=\underset{n}{}[r^2_n𝐫_n^2],$$
(2)
where $`𝐫_n=w_{n\mathrm{𝟎}}|𝐫|w_{n\mathrm{𝟎}}=iV/(2\pi )^3_{BZ}𝑑𝐤u_{n,𝐤}|\mathbf{}_𝐤|u_{n,𝐤}`$ and $`r^2_n=w_{n\mathrm{𝟎}}|r^2|w_{n\mathrm{𝟎}}_n`$. In these expressions, $`u_{n,𝐤}(𝐫)=e^{i𝐤𝐫}_mU_{mn}^{(𝐤)}\psi _{m𝐤}(𝐫)`$ is a periodic function which can be obtained from the Bloch functions of the composite group of bands. Practically, one calculates the BFs on equispaced Monkhorst–Pack Monkhorst76  meshes of $`𝐤`$ points in the unit cell $`\overline{\text{BZ}}`$ (whose volume is equal to the conventional BZ one), built on the reciprocal lattice vectors. The grids have been offset in order to include $`\mathrm{\Gamma }`$. With our chosen cell geometry of MnO, the reciprocal lattice vectors are: $`𝐠_1=\pi /a(3,1,1)`$, $`𝐠_2=\pi /a(1,\mathrm{\hspace{0.25em}3},1)`$, and $`𝐠_3=\pi /a(1,1,\mathrm{\hspace{0.25em}3})`$. We then express the matrix elements of the gradient $`_𝐤`$ appearing in the localization functional in terms of finite differences. As shown in Ref. Marzari97, , the only information needed is the overlap matrix $`M_{nn^{}}^{(𝐤,𝐛)}=u_{n𝐤}|u_{n^{},𝐤+𝐛}`$, where $`𝐛`$ are vectors connecting each mesh point to its nearest neighbors. From these quantities, the gradient of the spread functional $`\mathrm{\Omega }`$ with respect to an infinitesimal unitary transformation $`\delta U_{mn}^{(𝐤)}`$ of the BFs can also be evaluated. Once this gradient is computed, the minimization can take place via a steepest descent or conjugate–gradient algorithm.

It is important to stress that if the minimum of $`\mathrm{\Omega }`$ is flat (as it will be shown to be the case for some of the MnO WFs), its precise location may be hindered by the various numerical approximations involved in the calculations. In this work, we have especially taken care of the following features, which clearly become irrelevant in the limit of a dense mesh ($`N\mathrm{}`$, $`b0`$): (i) the regular mesh $`\{𝐤\}_{\overline{BZ}}`$ has been chosen in order to have the lattice symmetry, i.e. it transforms into itself by application of the point group operations $`R`$ of the crystal; (ii) the shell of $`𝐛`$ vectors used in the finite–difference formula for $`\mathbf{}_𝐤`$ has been constructed on the 12 vectors of type $`\{\mathrm{3\hspace{0.25em}3\hspace{0.25em}2}\}`$ and $`\{1\overline{1}\mathrm{\hspace{0.25em}0}\}`$ (in internal units). They are the sides midpoints of cubes which can be built on particular points of the rhombohedral mesh, and this set of 12 points is therefore compatible with the highest possible symmetry.

Finally, and following the procedure described in Eqs. (62)–(64) of Ref. Marzari97, , we start the minimization procedure by constructing a set of “trial functions” in the unit cell, which are an initial guess of the final WFs. We use Gaussians centered on atomic sites and modulated by an appropriate combination of spherical harmonics, with a rms width value such that the gaussian is negligible outside the corresponding atomic sphere. A unitary rotation among the initial BFs is then performed in order to maximize their projection on the trial functions.

It has been shown Krueger72  that Wannier functions can be chosen with well–defined symmetry properties, and forming (generally reducible) representations of the crystal point group. In order to investigate the nature and possible reductions of these representations in the case of MnO, and also to be in the position to use group–theory methods for further analysis, we need to calculate the representation matrix elements $`w_{n^{}\mathrm{𝟎}}|P_R|w_{n\mathrm{𝟎}}`$, where $`P_R`$ is the transformation operator corresponding to the operation $`R`$ of the point group, and the WFs are in the central cell. For simplicity, we restrict the formalism to symmorphic space groups (which is actually the case for D$`{}_{3\text{d}}{}^{}{}_{}{}^{5}`$). Using Eq. (1) with the $`𝐤`$ points discretization and $`𝐑=\mathrm{𝟎}`$, we have for the rotated Wannier function

$$P_R|w_{n\mathrm{𝟎}}=\frac{1}{N}\underset{𝐤\overline{BZ}}{}\underset{m}{}U_{mn}^{(𝐤)}P_R|\psi _{m𝐤},$$
(3)
and, with the “periodic gauge” condition $`\psi _{mR𝐤}(𝐫)=\psi _{mR𝐤+𝐆_𝐤}(𝐫)`$

$`w_{n^{}\mathrm{𝟎}}|P_R|w_{n\mathrm{𝟎}}_{NV}`$ $`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{\stackrel{𝐤,𝐤^{}}{\overline{BZ}}}{}}{\displaystyle \underset{m,m^{}}{}}U_{m^{}n^{}}^{(𝐤^{})}U_{mn}^{(𝐤)}\psi _{m^{}𝐤^{}}|P_R|\psi _{m𝐤}_{NV}`$ (4)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤\overline{BZ}}{}}{\displaystyle \underset{m,m^{}}{}}U_{m^{}n^{}}^{(R𝐤+𝐆_𝐤)}U_{mn}^{(𝐤)}\psi _{m^{}R𝐤+𝐆_𝐤}|P_R|\psi _{m𝐤}_V.`$
The last braket is trivial only in the case of non–degenerate BFs. In the general case, on the other hand, we have

$`P_R\psi _{m𝐤}(𝐫)`$ $`=`$ $`\psi _{m𝐤}(R^1𝐫)`$ (5)
$`=`$ $`{\displaystyle \underset{m^{}}{}}D_{m^{}m}(R,𝐤)\psi _{m^{}R𝐤}(𝐫),`$
where $`D_{m^{}m}(R,𝐤)`$ is a unitary transformation associated with the symmetry operation $`R`$, and the summation is over states degenerate with $`\psi _{m^{}R𝐤}(𝐫)`$. If this latter state is non–degenerate, $`D_{m^{}m}(R,𝐤)=\delta _{m^{}m}`$, and Eq. (5) becomes the usual formula. We note that the above formulas allow to decompose $`|w_{n\mathrm{𝟎}}`$ into a linear combination of basis functions transforming according to the relevant irreducible representations of the point group.

## III RESULTS

For our Wannier–function description of the electron states in AFM MnO, we have considered as a single composite group the topmost thirteen valence bands (O $`2s`$, O $`2p`$, Mn $`3d`$). In this way, the corresponding Wannier functions display the highest localization, and therefore can be considered as the “elementary building blocks” of the occupied electron states. As it is convenient to interpret some of our results in terms of combinations of atomic orbitals, we give in Table 1 the relevant O$`_\text{h}`$ representations corresponding to $`s`$, $`p`$, and $`d`$ orbitals and their reduction into irreducible representations of D$`_{3\text{d}}`$ and C$`_{3\text{v}}`$, together with their corresponding basis. Cartesian coordinates in the parent cubic reference system are written using $`x`$, $`y`$ and $`z`$ symbols, while $`X`$, $`Y`$ and $`Z`$ correspond to the trigonal reference system ($`Z`$ axis along the cubic  direction). Throughout this work, we use the Bethe notation for the irreducible representations Koster63 .

We have first studied the convergence of the spread functional $`\mathrm{\Omega }`$ in terms of the mesh density, using a sampling of the Brillouin zone $`\overline{\text{BZ}}`$ defined by the Monkhorst–Pack meshes $`\nu \times \nu \times \nu `$ with $`\nu =2`$, 4, 6, and 8. Analysis has been done on the various terms of the localization functional $`\mathrm{\Omega }`$ discussed in Ref. Marzari97, . We found that a satisfactory level of convergence is achieved using the 6$`\times `$6$`\times `$6 mesh, which represents a good compromise between accuracy and computer burden. Also, such a mesh is large enough to prevent difficulties related to finite grids and occurence of periodic WFs replica. The 6$`\times `$6$`\times `$6 sampling has been retained throughout this work.

The minimization of the localization functional, which determines the unitary matrix $`U_{mn}^{(𝐤)}`$, has been performed using a mixed strategy: a simple fixed–step steepest–descent procedure is used during the first iterations, followed by several iterations with a conjugate–gradient procedure which is reset to steepest descent every 100 iterations. Because of the delicate convergence of the Mn WFs, and in order to be on the safe side, we used a very large number (40000$``$80000) of iterations.

Several choices of trial orbitals have been tested to initialize efficiently the minimization procedure. In particular (see Section II), we have considered the following sets of initial trial orbitals: (i) Gaussians times real combinations of spherical harmonics $`s`$, $`p_x`$, $`p_y`$, and $`p_z`$ on the O1 and O2 Oxygen sites, and $`d_{z^2}`$, $`d_{x^2y^2}`$, $`d_{xy}`$, $`d_{yz}`$, and $`d_{zx}`$ on the Manganese Mn1 site. These real harmonics, reported in Table 1, form the basis of reducible representations of C$`_{3\text{v}}`$ (O site) and of D$`_{3\text{d}}`$ (Mn1 site), and are oriented along cubic axes; (ii) same centers as above, but with Gaussians modulated by spherical harmonics combinations $`s`$, $`p_X`$, $`p_Y`$, $`p_Z`$, and $`d_{Z^2}`$, $`d_{X^2Y^2}`$, $`d_{XY}`$, $`d_{YZ}`$, and $`d_{ZX}`$, i.e. basis functions of the irreducible representations of C$`_{3\text{v}}`$ and D$`_{3\text{d}}`$; (iii) combinations of the two above possibilities. All these sets of initial trial orbitals eventually lead, at convergence, to the same WFs, which are real, in agreement with the discussion in Ref. Marzari97, . However, the set (i) corresponds to a much faster convergence than the other choices. The reason is directly related to the shape and orientation of the final O $`2s`$/$`2p`$ WFs, and this point will be discussed further below.

We give in Table 2 the converged Wannier centers and spreads corresponding to the four Oxygen WFs at the O1 site for the spin–up chanel.

$`\overline{r}_x`$, $`\overline{r}_y`$, and $`\overline{r}_z`$ are the components of the Wannier center position with respect to the corresponding atomic site. We give also the total spread $`\mathrm{\Omega }`$ of these four WFs, and the sum of the corresponding centers whose interpretation will be discussed below. WFs related to the O2 site have the same $`\mathrm{\Omega }`$ and opposite $`\overline{𝐫}`$. The spreads corresponding to the five LSD WFs (9–13) at the Mn1 site are given in the first line of Table 3.

By symmetry, these WFs are centered on the Mn1 site (origin). Because the two Oxygen sites are equivalent, corresponding results for the spin–down chanel can be simply obtained by reversing Mn1 and Mn2. The total spread (sum of the 13 individual spreads) is 9.332 Å2. We note also that the sum of centers has all cartesian components equal to zero. Indeed, this quantity represents the electronic polarization (modulo a lattice vector) corresponding to valence states, and has to be zero as the system is centrosymmetric. In order to analyze in detail Table 2, it is useful to construct the matrix elements of the representation of the symmetry operations of the D$`_{3\text{d}}`$ point group on the basis of the 13 Wannier functions. This is done by using Eq. (4). First, we have verified by constructing the multiplication table that these 12 matrices form indeed a representation of D$`_{3\text{d}}`$. For non–trivial operations of the point group D$`_{3\text{d}}`$, they display the following structure: (i) there is an 8$`\times `$8 block corresponding to O $`2s`$/$`2p`$ Wannier functions. Each line of this block consists of zeros, except for one element which is equal to 1. Because the origin is on the Mn1 site, the operations of D$`_{3\text{d}}`$ are simply transforming all WFs on both O sites into each other; (ii) there is a 5$`\times `$5 block corresponding to Mn1 $`3d`$ Wannier functions. This block does not display a well–defined structure, and in general all its elements are non–zero; (iii) the two off–diagonal blocks (connecting O $`2s`$/$`2p`$ and Mn1 $`3d`$ WFs) have all their elements equal to zero, at an accuracy better than $`10^7`$. Therefore, the 13$`\times `$13 representation is block–diagonal, and in particular the 5$`\times `$5 Mn1 $`3d`$ block makes a (reducible) representation by itself.

We discuss first the O $`2s`$/$`2p`$ Wannier functions. The two non–equivalent WFs corresponding to site O1 are shown in Fig. 2, with a viewpoint which is about the same as in Fig. 1. The equivalent WFs corresponding to site O2 can simply be obtained by inversion symmetry with respect to the origin. In the absence of discrimination between Mn1 and Mn2 sites, the Oxygen site would have cubic (O$`_\text{h}`$) symmetry, and the 4 converged WFs originating from $`s`$ and $`p`$ states on the same O1 site would be $`sp^3`$ hybrids. This is a general result, which is demonstrated in Appendix A: from localized atom–centered $`s`$$`p`$ orbitals, the maximal localization algorithm always favors $`sp^3`$–like Wannier functions, provided that the symmetry is cubic or higher. This situation occurs, e.g., in a hypothetical ferromagnetic (FM) MnO crystal, which will be briefly described later. In AFM MnO, Mn1 and Mn2 sites are non–equivalent, so that the symmetry of the Oxygen site is reduced to C$`_{3\text{v}}`$. However, reminiscence of the $`sp^3`$ hybrids persists: disregarding for the moment the $`d`$ admixture, we see in panel (a) of Fig. 2 that the WF 1, with $`\mathrm{\Gamma }_1`$ symmetry, displays clearly the shape of a non–symmetric $`p_z`$ orbital, with its center shifted from the atomic site along the  direction ($`p_z`$$`s`$ mixing), as also indicated in Table 2. The three other WFs are equivalent to each other under the operations of C$`_{3\text{v}}`$ (they have the same spread), and their centers are located much closer to the cubic axes than to the diagonals. This can be explained by considering the interaction with the $`3d`$ orbitals of Mn2: all these WFs display a substantial bonding $`3d`$ contribution from the three neighboring unoccupied Mn2 sites. As the first–neighbor shell of Mn sites is an octahedron centered on each O, the above interaction favors the orientation of the equivalent WFs 2, 3, and 4 along cubic axes. We display in panel (b) of Fig. 2 the WF 2 whose center is close to the  axis. The sum of WFs centers given in Table 2 is a measure of the trigonal distortion of the $`sp^3`$–like orbitals, since it would be vanishing in a perfect cubic environment. This configuration of the converged WFs justifies also our choice of trial orbitals: $`p_x`$, $`p_y`$, and $`p_z`$ orbitals are already good guesses of final WFs 2, 3, and 4. During the iterations, most of the effort is devoted to building up WF 1 (the $`p_z`$$`s`$ hybrid oriented along ). We note also in Fig. 2 small features on the neighboring Mn1 sites. They correspond to a weak bonding contribution from unoccupied states on these sites.

We consider next the five Mn1 $`3d`$ Wannier functions. They are exactly centered at the Mn1 site and their spreads do not display any well–defined symmetry. When convergence of O $`2s`$/$`2p`$ Wannier functions is practically achieved, further iterations of the localization process have the only effect of modifying slightly (esentially by amounts comparable to the numerical noise) the individual spreads of the five Mn1 $`3d`$ WFs, while their total spread remains unchanged. These five WFs are not uniquely defined by the criterion of maximal localization. In fact, any set of WFs obtained by applying a unitary transformation to the original set, will in general have different individual spreads, but the same value of the total spread. Since it is advisable that the final WFs display explicitely as much symmetry as possible, we exploit this non–uniqueness and perform a symmetry reduction of the 5–fold representation of the Mn1 $`3d`$ WFs, which decomposes into $`\mathrm{\Gamma }_1^++2\mathrm{\Gamma }_3^+`$ according to D$`_{3\text{d}}`$ symmetry.

We note that the Hermitian operator $`r^2`$, which appears in the localization functional Eq. (2), has the full rotational symmetry, and therefore its matrix representation $`r^2_{n^{}n}=w_{n^{}\mathrm{𝟎}}|r^2|w_{n\mathrm{𝟎}}`$ belongs to the $`\mathrm{\Gamma }_1^+`$ representation of D$`_{3\text{d}}`$. These matrix elements can be calculated through a simple generalization of Eq. (23) of Ref. Marzari97, , and are given by

$`w_{n^{}\mathrm{𝟎}}|r^2|w_{n\mathrm{𝟎}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b[2\delta _{n^{}n}M_{n^{}n}^{(𝐤,𝐛)}`$ (6)
$`M_{nn^{}}^{(𝐤,𝐛)}].`$
We choose the five new WFs as the linear combinations of the original set which diagonalize the $`r^2`$ matrix, Eq. (6). We get in this way an extra unitary transformation which allows us to update the matrices $`U^{(𝐤)}`$. The resulting WFs, which display the expected ($`\mathrm{\Gamma }_1^++2\mathrm{\Gamma }_3^+`$) symmetry, are indicated with the header “s–LSD” in Table 3, and are shown in Fig. 3. We note that an arbitrariness (up to a two–by–two unitary transformation) remains in the definition of the partners in the two $`\mathrm{\Gamma }_3^+`$ representations. The figure shows that the WFs 9–13 are essentially atomic orbitals modified by the D$`_{3\text{d}}`$ crystal field.

The above picture depicting the whole set of WFs in AFM MnO has been confirmed by similar calculations performed for a hypothetical FM MnO crystal. In this case, the symmetry is cubic (space group O$`{}_{\text{h}}{}^{}{}_{}{}^{5}`$), and the majority–spin chanel contains 18 valence states (O $`2s`$, O $`2p`$, and Mn $`3d`$), while the minority–spin chanel has 8 valence states (O $`2s`$, O $`2p`$). As in AFM MnO, the Mn $`3d`$ WFs are found to be essentially atomic orbitals. By contrast, the O $`2s/2p`$ WFs in both chanels are here undistorted $`sp^3`$ hybrids (see Appendix A). In the minority–spin chanel, both lobes of each hybrid display a substantial bonding $`3d`$ contribution from the three neighboring unoccupied Mn sites, while these features are absent in the majority–spin chanel (occupied Mn sites). The $`sp^3`$–like Oxygen WFs, with center along the  direction, are shown in Fig. 4 for both spin chanels.

The covalent interaction in AFM MnO between O and Mn2 $`3d`$ orbitals, clearly singled out in Fig. 2 for the case of the spin–up chanel (we remind that corresponding results for spin–down are obtained by reversing Mn1 and Mn2), is the fingerprint of superexchange in a one–electron picture. Superexchange mechanism results in fact from indirect many–body interactions, involving metal atoms (here Mn1 and Mn2) with an intervening Oxygen. The shared covalency of nearest–neighbor pairs of magnetic ions leads to an antiferromagnetic alignment of their moments. As a matter of fact, supposing, e.g. a spin–up WF on the Mn1 site, the neighboring Mn WF has the choice between two spin orientations, corresponding to a FM or an AFM ordering. The FM choice will be energetically less favourable, due to the required orthogonalization to the Mn1 WF, while in the AFM alignment, the spin part of the wavefunction accounts for this necessary orthogonalization. The above textbook explanation of AFM ordering can be discussed in light of our calculations, showing for the first time a real–space representation of the WFs for an AFM system. The maximally–localized WFs bring out the relevant physical information, namely the covalent interaction with empty Mn $`3d`$ states. We also notice that the FM alignment, on the other hand, is favoured by kinetic energy. This can be seen in Figs. 2 and 4 from the larger number of covalent interactions (6 instead of 3 in AFM case) associated with the O $`sp^3`$–like WFs, which leads eventually to a broader band–width.

It is useful at this point to contrast this rather simple one–electron picture of AFM MnO, obtained in real space using the localized Wannier–function description, with the more entangled picture based on the energy band–structure of the extended Bloch–function description. For this purpose, we show in the left panel of Fig. 5 the spin–up band structure of MnO along the $`X\mathrm{\Gamma }Z`$ directions, calculated within the LSD scheme.

Starting from the lowest energies, we find first a $`\mathrm{\Gamma }_1^+`$ and a $`\mathrm{\Gamma }_2^{}`$ states at $`𝐤=\mathrm{𝟎}`$. They arise from $`\mathrm{\Gamma }_1`$$`2s`$ orbitals on the two O sites (related by inversion symmetry), forming symmetric ($`\mathrm{\Gamma }_1^+`$) and antisymmetric ($`\mathrm{\Gamma }_2^{}`$) combinations, respectively. The primary origin of the next composite group of 6 bands is O $`2p`$ triplets on the two O sites, which are split into a $`\mathrm{\Gamma }_1`$ singlet and a $`\mathrm{\Gamma }_3`$ doublet by trigonal C$`_{3\text{v}}`$ site symmetry. They lead to one set of symmetric combinations ($`\mathrm{\Gamma }_1^+`$, $`\mathrm{\Gamma }_3^+`$), and one set of antisymmetric combinations ($`\mathrm{\Gamma }_2^{}`$, $`\mathrm{\Gamma }_3^{}`$) between the two sites. Next, the occupied $`t_{2g}`$ orbitals of the Mn1 atom are split at $`\mathrm{\Gamma }`$ into a doublet plus a singlet ($`\mathrm{\Gamma }_3^+`$, $`\mathrm{\Gamma }_1^+`$), while the $`e_g`$ orbitals remain a doublet ($`\mathrm{\Gamma }_3^+`$). These Mn1 $`3d`$ orbitals give rise to the uppermost occupied composite group of five bands. Clearly, states of the Mn1 $`3d`$ and O $`2p`$ complexes which belong to the same $`\mathrm{\Gamma }_1^+`$ (or $`\mathrm{\Gamma }_3^+`$) representation interact with each other, so that eventually both sets $`\mathrm{\Gamma }_1^+`$(Mn1 $`d`$), $`\mathrm{\Gamma }_1^+`$(O $`p`$) and $`\mathrm{\Gamma }_3^+`$(Mn1 $`d`$), $`\mathrm{\Gamma }_3^+`$(O $`p`$) consist of bonding and antibonding partners of Mn1 $`3d`$/O $`2p`$ hybrids. A mirror picture (relative to Fermi energy) of the above description applies to the conduction bands which originate from the Mn2 $`3d`$ orbitals, if we exclude the lowest conduction band at $`\mathrm{\Gamma }`$, which is free–electron–like, has mostly Mn $`4s`$/O $`3s`$ character, and is not relevant to the present discussion. The above description would leave completely empty the spin–up $`3d`$ states of the Mn2 atom. However, a covalent interaction (clearly depicted in Fig. 2 within the Wannier–function description) between O $`2p`$ and the unoccupied Mn2 $`3d`$ states takes place, and some amount of spin–up charge appears at the Mn2 sites, which would be totally empty in a completely ionic picture. For instance, the highest antibonding $`\mathrm{\Gamma }_3^+`$ valence state at $`\mathrm{\Gamma }`$ has 60% $`d`$–occupancy in the Mn1 sphere, 7% $`p`$–occupancy in each of the two O spheres, and 14% $`d`$–occupancy in the Mn2 sphere.

With the above results in mind, it is instructive to investigate WFs within other one–electron schemes, leading to different degrees of wavefunction localization, covalent interaction, and magnetic moments. Data corresponding to the LSD + $`U`$ scheme are presented in Tables 2 and 3 for the spin–up chanel. Values for the Mn1 $`3d`$ WFs are those obtained after diagonalization of $`r^2_{n^{}n}`$. All individual spreads, as well as the total spread (whose value is 8.083 Å2), are smaller than in the LSD case: the O $`2s`$/$`2p`$ and Mn $`3d`$ Wannier functions are more localized by 9% and 22%, respectively, and this is consistent with the narrower valence band width. This property appears also in Fig. 6, where we show the LSD + $`U`$ result for the WF 2 derived from O $`2s`$/$`2p`$ states.

This LSD + $`U`$ orbital is similar to the corresponding one calculated within LSD, and displayed in Fig. 2b. The most noticeable difference appears in the weaker covalent interaction between O and the Mn2 $`3d`$ orbitals. This is clearly visible in the right panel of Fig. 6, where we display the $`d`$–like differential charge density ($`\rho _{\text{LSD}+U}\rho _{\text{LSD}}`$) associated with this particular Oxygen WF. The five Mn1 $`3d`$ WFs, obtained within LSD + $`U`$, are very similar to the LSD ones, and are essentially atomic orbitals modified by the crystal field.

The more localized character of all WFs within the LSD + $`U`$ scheme, together with the weaker covalent interaction between O and Mn2 $`3d`$ orbitals, is the fingerprint in the Wannier–function description of the larger on–site interaction $`U`$ of cation $`3d`$ states. It is worth to realize that this simple trend is more easily detectable in terms of Wannier functions than Bloch functions. To this end, we show in the right panel of Fig. 5 the LSD + $`U`$ band structure of MnO. It differs in several respects from the LSD one. The energy gap (2.0 eV), though still smaller than the experimental one ($`3.84.2`$ eV), is closer to it than the LSD value (0.8 eV). The energy separation between occupied and empty $`d`$ states, which is due to the on–site interaction $`U`$ between cation 3$`d`$ states, has also increased. Furthermore, the three Mn1 $`3d`$ bands corresponding to the lower $`\mathrm{\Gamma }_3^+`$ and $`\mathrm{\Gamma }_1^+`$ states ($`t_{2g}`$) do not form any more a composite group of bands by themselves, but rather mix with the O $`2p`$ complex. This results in a narrowing of the valence band width (6.3 eV, against 7.2 eV in LSD).

The distinct features of the WFs computed within LSD and LSD + $`U`$ can be used for discussing the behavior of the Born dynamical charge tensor $`𝖹^{}`$ within these different many–body approximations. In fact, as $`𝖹^{}`$ can be defined in terms of the relative displacement of WFs centers Vanderbilt93 , it carries along a relevant information. In particular, a $`p`$$`d`$ interaction between valence and conduction states in mixed ionic–covalent materials Posternak97  leads to anomalous $`𝖹^{}`$, while atomic–like occupied WFs correspond to a nominal value. Therefore in MnO, one expects smaller $`|𝖹^{}|`$ (i.e. closer to 2, its nominal value) within the LSD + $`U`$ scheme, to which correspond more localized WFs, than within LSD. We performed the calculation of the  component of this tensor, using the WF description, and obtained within LSD the value $`|Z_{}^{}|=2.36`$, which is consistent with the one obtained previously Massidda99  using the Berry phase approach. LSD + $`U`$, on the other hand, gives as expected a smaller value $`|Z_{}^{}|=2.26`$. These results relate the behavior of two important physical quantities (the superexchange energy and $`𝖹^{}`$) to the same information contained in the AFM MnO WFs, namely the magnitude of their amplitudes on unoccupied Mn2 sites in the different approximations.

## IV CONCLUSION

We have used the maximum localization criterion for constructing Wannier functions of AFM MnO within both LSD and LSD + $`U`$ schemes. This is the first application of the method to a magnetic low–symmetry system, with a partially filled $`d`$ band. The LSD + $`U`$ approach has been considered here, as it improves the LSD description of the ground–state properties of strongly correlated materials. The Wannier–function description, which can be viewed as a generalization to extended systems of the localized molecular–orbitals description of molecules, provides in the case of MnO a simple picture complementary to the Bloch–function one. We obtain two distinct groups of WFs: (i) the five Mn $`3d`$ WFs centered on the Mn1 site which are atomic orbitals modified by the crystal field; (ii) the four trigonally distorted $`sp^3`$ WFs associated with each of the O1 and O2 sites, which are not centered on these sites. The latter display a substantial covalent bonding with $`3d`$ states on Mn2 sites, which is consistent with the AFM ordering of this material. Concerning the Mn WFs, they are not uniquely defined by the criterion of maximal localization. We have exploited this feature to obtain symmetrized maximally–localized WFs. We impose this further condition by diagonalization of the $`r^2`$ operator.

###### Acknowledgements.

It is a pleasure to acknowledge stimulating discussions with D. Vanderbilt. We would like to thank N. Bernstein for providing us with his graphics program “dan”. This work was supported by the Swiss National Science Foundation (Grant No. 20–59121.99).

## Appendix A MAXIMALLY LOCALIZED WANNIER FUNCTIONS FROM ONE S–LIKE AND THREE P–LIKE ORBITALS FOR CUBIC SITE SYMMETRY.

We consider the 4–dimensional Hilbert space $`_4`$ spanned by 4 orthonormal localized basis functions $`|\varphi _1`$, $`|\varphi _2`$, $`|\varphi _3`$, and $`|\varphi _4`$, centered on a given site (taken here as the origin) with cubic point–group symmetry (T, T$`_\text{h}`$, O, T$`_\text{d}`$, or O$`_\text{h}`$). $`|\varphi _1`$ is the basis function of a representation $`\mathrm{\Gamma }_1`$ (for T, O, and T$`_\text{d}`$) or $`\mathrm{\Gamma }_1^+`$ (for T$`_\text{h}`$ and O$`_\text{h}`$). $`|\varphi _2`$, $`|\varphi _3`$, and $`|\varphi _4`$ are the basis functions of a threefold representation $`\mathrm{\Gamma }_4`$ (for T and O), $`\mathrm{\Gamma }_5`$ (T$`_\text{d}`$), $`\mathrm{\Gamma }_4^+`$ (T$`_\text{h}`$), or $`\mathrm{\Gamma }_4^{}`$ (O$`_\text{h}`$). More precisely, $`|\varphi _2=|\varphi _x`$ $``$ $`x`$–row of the threefold representation, and $`|\varphi _3=C_3^1|\varphi _x`$, $`|\varphi _4=C_3|\varphi _x`$, where $`C_3`$ is the rotation by $`2\pi /3`$ around the $`[\mathrm{1\hspace{0.33em}1\hspace{0.33em}1}]`$ axis. We consider real $`|\varphi _i`$ functions, as the maximally–localized WFs turn out to be real, apart from an arbitrary overall phase factor Marzari97 . Finally, we suppose that the $`|\varphi _i`$ $`(i=1,\mathrm{},4)`$ are sufficiently localized so that $`\varphi _i|𝐫|\varphi _j`$ and $`\varphi _i|r^2|\varphi _j`$ exist and are finite.

We demonstrate below that minimizing within $`_4`$ the localization functional, Eq. (2), leads to $`sp^3`$–like Wannier functions, as opposed to the original atomic–like $`|\varphi _i`$ functions. The solution is unique and corresponds to one of the two tetrahedra defined by the cubic point–symmetry group if the latter is T or T$`_\text{d}`$. For the higher symmetry groups T$`_\text{h}`$, O, and O$`_\text{h}`$, maximal localization leads to an infinity of $`sp^3`$–like solutions constructed from $`|\varphi _1`$ and unitary transformed three–fold functions $`|\varphi _2`$, $`|\varphi _3`$, $`|\varphi _4`$ according to any of the $`3\times 3`$ orthogonal matrices $`𝐃^{(1)}(\alpha ,\beta ,\gamma )`$ representing a rotation of the cartesian co–ordinates of a vector.

As the operator $`r^2`$ belongs to the identity representation, it is diagonal in the $`|\varphi _i`$ basis, and we have

$$\varphi _i|r^2|\varphi _i=\left(\begin{array}{cccc}\mathrm{\Omega }_s& \mathrm{\Omega }_p& \mathrm{\Omega }_p& \mathrm{\Omega }_p\end{array}\right),$$
(7)
where cubic symmetry has been taken into account. Similarly, writing $`𝐫=x𝐢+y𝐣+z𝐤`$ and using projection operators, we find after some algebra

$$\varphi _i|𝐫|\varphi _j=\left(\begin{array}{cccc}0& \lambda 𝐢& \lambda 𝐣& \lambda 𝐤\\ \lambda 𝐢& 0& \delta 𝐤& \delta 𝐣\\ \lambda 𝐣& \delta 𝐤& 0& \delta 𝐢\\ \lambda 𝐤& \delta 𝐣& \delta 𝐢& 0\end{array}\right),$$
(8)
with $`\lambda =\varphi _s|\xi |\varphi _\xi `$, where $`\xi =x,y,z`$, and $`\delta =\varphi _x|y|\varphi _z=\varphi _y|z|\varphi _x=\varphi _z|x|\varphi _y`$. The phase arbitrariness of the orbitals allows us to set $`\lambda 0`$. We note that symmetry requires $`\delta =0`$ in the case of T$`_\text{h}`$, O, or O$`_\text{h}`$. The total spread of the original atomic–like $`|\varphi _i`$ functions is

$$\mathrm{\Omega }_\varphi =\underset{i}{}[\varphi _i|r^2|\varphi _i|\varphi _i|𝐫|\varphi _i|^2]=\mathrm{\Omega }_s+3\mathrm{\Omega }_p.$$
(9)
We consider now the four orbitals $`|\psi _i`$ obtained from the $`|\varphi _j`$ by a general orthogonal transformation $`O`$, given explicitely by

$$|\psi _i=\underset{j}{}O_{ij}|\varphi _j(i,j=1,\mathrm{},4).$$
(10)
There are 6 independent matrix elements, and 10 orthonormality conditions for columns or rows. Using Eqs. (7), (10) and the orthonormality conditions, we obtain

$`{\displaystyle \underset{i}{}}\psi _i|r^2|\psi _i`$ $`=`$ $`{\displaystyle \underset{i}{}}[O_{i1}^2\mathrm{\Omega }_s+(O_{i2}^2+O_{i3}^2+O_{i4}^2)\mathrm{\Omega }_p]`$ (11)
$`=`$ $`{\displaystyle \underset{i}{}}[O_{i1}^2\mathrm{\Omega }_s+\mathrm{\Omega }_p(1O_{i1}^2)]`$
$`=`$ $`\mathrm{\Omega }_s+3\mathrm{\Omega }_p,`$
which verifies the trace invariance property of Eq. (7) under orthogonal transformations, and gives for the total spread $`\mathrm{\Omega }_\psi `$ of the functions defined by (10) the inequality

$`\mathrm{\Omega }_\psi `$ $`=`$ $`{\displaystyle \underset{i}{}}\psi _i|r^2|\psi _i{\displaystyle \underset{i}{}}|\psi _i|𝐫|\psi _i|^2`$ (12)
$`=`$ $`\mathrm{\Omega }_s+3\mathrm{\Omega }_p{\displaystyle \underset{i}{}}|\psi _i|𝐫|\psi _i|^2`$
$``$ $`\mathrm{\Omega }_s+3\mathrm{\Omega }_p=\mathrm{\Omega }_\varphi .`$
Therefore, the atomic–like $`|\varphi _i`$ functions correspond to the maximum value of the spread, and the problem of the total spread minimization reduces to minimizing the second term $`_i|\psi _i|𝐫|\psi _i|^2`$ of the localization functional, which can be evaluated using Eqs. (8) and (10).

We consider next the two particular orthogonal transformations $`\stackrel{~}{O}`$ and $`\stackrel{~}{O^{}}`$ of the $`|\varphi _i`$ corresponding to the two possible choices for the four $`sp^3`$ hybrid orbitals $`|\stackrel{~}{\psi }_i`$ and $`|\stackrel{~}{\psi }_i^{}`$ (written together below $`|\stackrel{~}{\psi }_i^{()}`$), and given explicitely by

$$\{\begin{array}{c}|\stackrel{~}{\psi }_{111}^{()}=\frac{1}{2}(|\varphi _1\pm |\varphi _2\pm |\varphi _3\pm |\varphi _4)\\ |\stackrel{~}{\psi }_{1\overline{1}\overline{1}}^{()}=\frac{1}{2}(|\varphi _1\pm |\varphi _2|\varphi _3|\varphi _4)\\ |\stackrel{~}{\psi }_{\overline{1}1\overline{1}}^{()}=\frac{1}{2}(|\varphi _1|\varphi _2\pm |\varphi _3|\varphi _4)\\ |\stackrel{~}{\psi }_{\overline{1}\overline{1}1}^{()}=\frac{1}{2}(|\varphi _1|\varphi _2|\varphi _3\pm |\varphi _4)\end{array},$$
(13)
where $`i=1,\mathrm{},4=(111),\mathrm{},(\overline{1}\overline{1}1)`$. Using the $`\stackrel{~}{O}`$ transformation and Eq. (8), we get for the $`\stackrel{~}{\psi }_i|𝐫|\stackrel{~}{\psi }_j`$ matrix elements

$$\frac{1}{2}\left(\begin{array}{cccc}(\lambda +\delta )(𝐢+𝐣+𝐤)& (\lambda \delta )𝐢& (\lambda \delta )𝐣& (\lambda \delta )𝐤\\ (\lambda \delta )𝐢& (\lambda +\delta )(𝐢𝐣𝐤)& (\lambda \delta )𝐤& (\lambda \delta )𝐣\\ (\lambda \delta )𝐣& (\lambda \delta )𝐤& (\lambda +\delta )(𝐢+𝐣𝐤)& (\lambda \delta )𝐢\\ (\lambda \delta )𝐤& (\lambda \delta )𝐣& (\lambda \delta )𝐢& (\lambda +\delta )(𝐢𝐣+𝐤)\end{array}\right).$$
(14)
We note that taking the $`\stackrel{~}{O^{}}`$ transformation corresponds to changing $`\lambda `$ into $`\lambda `$ in Eq. (14) for the $`\stackrel{~}{\psi }_i^{}|𝐫|\stackrel{~}{\psi }_j^{}`$ matrix elements. Using this latter equation, together with (12), we obtain the total spreads of the $`|\stackrel{~}{\psi }_i`$ and $`|\stackrel{~}{\psi }_i^{}`$ WFs

$`\mathrm{\Omega }_{\stackrel{~}{\psi }}`$ $`=`$ $`{\displaystyle \underset{i}{}}[\stackrel{~}{\psi }_i|r^2|\stackrel{~}{\psi }_i|\stackrel{~}{\psi }_i|𝐫|\stackrel{~}{\psi }_i|^2]`$ (15)
$`=`$ $`\mathrm{\Omega }_s+3\mathrm{\Omega }_p3(\delta +\lambda )^2`$
$`=`$ $`\mathrm{\Omega }_\varphi 3(\delta +\lambda )^2\mathrm{\Omega }_\varphi ,`$
and

$$\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}=\mathrm{\Omega }_\varphi 3(\delta \lambda )^2\mathrm{\Omega }_\varphi .$$
(16)
Two cases have to be considered according as $`\delta `$ is zero (point groups T$`_\text{h}`$, O, O$`_\text{h}`$) or not (point groups T and T$`_\text{d}`$).

If $`\delta =0`$, the two spreads $`\mathrm{\Omega }_{\stackrel{~}{\psi }}`$ and $`\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}`$ are both equal to $`\mathrm{\Omega }_\varphi 3\lambda ^2`$. We demonstrate first that this latter value is the minimum of the spread, Eq. (12). To this end, using Eqs. (8), (10), and the orthonormality conditions $`_jO_{ij}^2=1`$ and $`_iO_{i1}^2=1`$, we are lead to minimizing

$`\mathrm{\Omega }_\psi `$ $`=`$ $`\mathrm{\Omega }_\varphi 4\lambda ^2{\displaystyle \underset{i}{}}O_{i1}^2(O_{i2}^2+O_{i3}^2+O_{i4}^2)`$ (17)
$`=`$ $`\mathrm{\Omega }_\varphi 4\lambda ^2+4\lambda ^2{\displaystyle \underset{i}{}}O_{i1}^4,`$
with the constraint

$$C=\underset{i}{}O_{i1}^21=0.$$
(18)
The conditions for an extremum of $`\mathrm{\Omega }_\psi `$, introducing the Lagrangian multiplier $`\mathrm{\Lambda }`$ are

$$\frac{\mathrm{\Omega }_\psi }{O_{i1}}+\mathrm{\Lambda }\frac{C}{O_{i1}}=16\lambda ^2O_{i1}^3+2\mathrm{\Lambda }O_{i1}=0,$$
(19)
from which we get

$$O_{11}^2=O_{21}^2=O_{31}^2=O_{41}^2=\frac{1}{4},$$
(20)
and the result

$$\mathrm{min}\mathrm{\Omega }_\psi =\mathrm{\Omega }_\varphi 3\lambda ^2=\mathrm{\Omega }_{\stackrel{~}{\psi }}\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}\mathrm{\Omega }_\varphi .$$
(21)
We demonstrate next that in addition to $`\stackrel{~}{O}`$ and $`\stackrel{~}{O^{}}`$, there is in fact an infinity of orthogonal transformations producing $`|\psi _i`$ orbitals with this same minimum spread value. For this purpose, we introduce the operator $`R(\alpha ,\beta ,\gamma )`$ corresponding to a geometrical rotation specified by its Euler angles, and its three–dimensional unitary representation $`D^{(1)}(\alpha ,\beta ,\gamma )_{mm^{}}`$, whose basis functions are the spherical harmonics $`Y_{1m}`$, with $`m=1,\mathrm{\hspace{0.33em}0},1`$. Consistently with our choice of the basis functions $`|\varphi _2`$, $`|\varphi _3`$, and $`|\varphi _4`$, we consider instead the equivalent, real orthogonal representation $`D^{(1)}(\alpha ,\beta ,\gamma )_{\xi \xi ^{}}`$ with $`\xi ,\xi ^{}=x,y,z`$, which is the usual transformation of the cartesian co–ordinates of a vector. We then build from this matrix the following orthogonal transformation acting in the 4–dimensional Hilbert space $`_4`$

$$𝓓(\alpha ,\beta ,\gamma )=\left(\begin{array}{cc}1& \begin{array}{cccc}0& \mathrm{\hspace{0.33em}\hspace{0.33em}0}& \mathrm{\hspace{0.33em}\hspace{0.33em}0}& \end{array}\\ & \\ \begin{array}{c}0\\ 0\\ 0\end{array}& D^{(1)}(\alpha ,\beta ,\gamma )_{\xi \xi ^{}}\end{array}\right).$$
(22)
We write $`|\widehat{\psi }_i=_k\stackrel{~}{O}_{ik}_j𝒟(\alpha ,\beta ,\gamma )_{kj}|\varphi _j`$, and get after some algebra, using Eqs. (8) and (10)

$$\mathrm{\Omega }_{\widehat{\psi }}\mathrm{\Omega }_{\stackrel{~}{\psi }}=\mathrm{\Omega }_\varphi 3\lambda ^2\alpha ,\beta ,\gamma .$$
(23)
It is important to remember that the $`|\varphi _i`$ are basis functions of irreducible representations of a cubic group, and not of the full rotation group. Therefore, applying the transformation (22) on the original basis functions $`|\varphi _i`$ does not result in rotated $`\varphi _i(𝐫)`$ functions, i.e., $`\varphi _i(R^1(\alpha ,\beta ,\gamma )𝐫)`$, but in linear combinations of the original basis functions, defining the same Hilbert space $`_4`$, and carrying the same chemical information as the original functions.

When $`\delta 0`$ (case of T and T$`_\text{d}`$), the two transformations $`\stackrel{~}{O}`$ and $`\stackrel{~}{O^{}}`$ produce $`|\stackrel{~}{\psi }_i`$ and $`|\stackrel{~}{\psi }_i^{}`$ orbitals with different total spreads

$$\begin{array}{ccc}\mathrm{\Omega }_{\stackrel{~}{\psi }}<\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}\hfill & \text{ if }\delta >0\hfill & (\lambda \delta 0)\hfill \\ \mathrm{\Omega }_{\stackrel{~}{\psi }}>\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}\hfill & \text{ if }\delta <0\hfill & (\lambda \delta 0)\hfill \end{array}\}.$$
(24)
We prove finally that $`\mathrm{\Omega }_{\stackrel{~}{\psi }}`$ for $`\delta >0`$ (or $`\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}`$ for $`\delta <0`$) is a local minimum of the total spread. To this end, we consider an infinitesimal orthogonal transformation of the $`|\stackrel{~}{\psi }_i`$ (or $`|\stackrel{~}{\psi }_i^{}`$)

$$O(ϵ)=I+ϵA+\frac{1}{2}ϵ^2A^2+𝒪(ϵ^3),$$
(25)
where $`A`$ is an antisymmetrical operator defined by 6 real parameters. Using Eq. (14) and applying (25) to the sets of $`sp^3`$ hybrid orbitals $`|\stackrel{~}{\psi }_i`$ and $`|\stackrel{~}{\psi }_i^{}`$, we obtain for the infinitesimal variations of the spreads, up to second order terms in $`ϵ`$, positive semi–definite quadratic forms in the parameters of the $`A`$ operator. Hence the correction to $`\mathrm{\Omega }_{\stackrel{~}{\psi }}`$ (or $`\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}`$) induced by the infinitesimal transformation, are non–negative. The spread $`\mathrm{\Omega }_{\stackrel{~}{\psi }}`$ ($`\mathrm{\Omega }_{\stackrel{~}{\psi }^{}}`$) therefore corresponds for $`\delta >0`$ ($`\delta <0`$) to a local minimum of the localization functional.


format: text

text:
# Spin-dependent forces of quarks in baryon

## 1 Introduction

The spin structure of baryons presents a still unsolved problem,both on partonic and quark model level. For excited baryon spectrum the apparent small spin-orbit splitting of some baryonic states is a topic of vivid discussions -. Some baryonic states, like Roper resonance $`N(1440)`$ or $`\mathrm{\Lambda }(1405)`$ are not yet explained in the traditional framework of relativistic quark model (RQM) .More detailed information about the spin structure of baryons comes from the polarization experiments on electroproduction of excited resonances,  which effectively measure the convolution of the baryon wave function, and is very sensitive to its structure.

Meanwhile the theoretical knowledge of the quark spin forces in the baryon is limited to the perturbative expressions calculated decades ago , and nonperturbative spin-orbit Thomas term, written in the framework of RQM .

In applying these results to light baryons the notion of constituent quark masses is introduced in RQM,these appear in spin-dependent (SD) forces and play role of fitting parameters.

It is the purpose of the present paper to derive SD forces in a baryon in a most straightforward way from the QCD Lagrangian with the nonperturbative vacuum described by vacuum field correlators . Limiting to the lowest (Gaussian) field correlators one can express all terms of SD forces through the scalar functions $`D`$ and $`D_1`$ representing this Gaussian correlator . High accuracy of such procedure is supported by recent lattice data  and contribution of higher correlators can be estimated to be of the order of few percent .

The functions $`D,D_1`$ are themselves measured on the lattice  and also found in analytic approaches .

An essential element of the present approach is that it is not connected to the heavy mass expansion, and can be applied also to light quarks in a baryon. In this case an effective Hamiltonian is constructed from first principles, which contains einbein (auxiliary) fields. It was shown previously, that stationary point of these einbein fields yields exactly the constituent quark masses which can be expressed unambigiously through the only parameter of this approach – the string tension $`\sigma `$.

The decisive check of this procedure is the calculation of baryon magnetic moments, since they are inversely proportional to the quark constituent masses. That was done in  and results agreed with all known experimental data within $`10\%`$.

The SD forces derived below in the paper are computed as a series in field correlators (cumulants) with growing powers of fields.

The lowest (Gaussian) term yields SD forces inversely proportional to the square of constituent quark masses. Having in mind the high accuracy of the Gaussian approximation  and baryon magnetic moments  one should expect that SD forces found below have the accuracy of the order of 10%.

Analogous expressions for heavy quarkonia  and light mesons  have been reported earlier and estimated for realistic meson system respectively in  and .

The plan of the paper is as follows. In chapter 2 a general expression for the 3q baryon Green’s function is introduced, and the Fock-Feynman-Schwinger (world-line) formalism is used to reveal the dependence on gauge fields with spin operators explicitly written. Averaging over those with the help of the Field Correlator Method (FCM) one finally obtains an expression for the Green’s function written in terms of field correlators.

In chapter 3 a special case of heavy quark masses is considered and all SD forces are obtained in closed form, expressed in terms of correlator functions $`D`$ and $`D_1`$.

In chapter 4 the perturbative contribution to SD forces is written down. In chapter 5 a general case is considered when current quark masses can also be vanishingly small, and SD forces are written again in terms of integrals over functions $`D`$ and $`D_1`$ with constituent (dynamical) masses entering in the denominator.

Chapter 6 is devoted to the discussion of the relativistic structure of SD forces in the excited baryon spectrum. Comparison to other results in the literature is also made and possible extension of the method is suggested. Main points of the paper are summarized in Conclusions.

## 2 $`3q`$ Green’s function with spin insertions

Following  we consider the $`3q`$ Green’s function, which can be written as

$$G_{3q}(x,y)=tr_L\left[\mathrm{\Gamma }_{out}\underset{i=1}{\overset{3}{}}(m_i\widehat{D})D\mu _iDz^{(i)}e^{K_i}W_3\mathrm{exp}g\sigma F\mathrm{\Gamma }_{in}\right]$$
(1)
where $`tr_L`$ is the trace over Dirac matrix indices, $`\mathrm{\Gamma }_{out}(\mathrm{\Gamma }_{in})`$ are final (initial) state operators creating given $`J^{PC}`$ assignment to the $`3q`$ state and we have also denoted as in

$$K_i=_0^T𝑑t\left[\frac{m_i^2}{2\mu _i(t)}+\frac{\mu _i(t)}{2}(\dot{𝐳}^2(t)+1)\right],$$
(2)
$$W_3\mathrm{exp}(g\sigma F)=tr_Y\mathrm{exp}[\underset{n=0}{\overset{\mathrm{}}{}}\frac{(ig)^n}{n!}F(1)\mathrm{}F(n)d\rho (1)..d\rho (n)].$$
(3)
Here and in what follows $`F(1)`$ is always implied to be gauge-transported to one point $`x_0`$,namely $`F(1)F(z(1),x_0)=\mathrm{\Phi }(x_0,z(1))F(z(1))\mathrm{\Phi }(z(1),x_0)`$ where $`\mathrm{\Phi }(x,y)`$ is defined in the Appendix, and finally, $`d\rho (n)=_{i=1}^3d\rho ^{(i)}(n)`$, with

$$d\rho ^{(i)}(n)ds_{\mu _n\nu _n}^{(i)}(u^{(n)})+\frac{1}{i}\sigma _{\mu _n\nu _n}^{(i)}\frac{dt_n}{2\mu _i(t_n)}.$$
(4)
The integration in (3) extends over all 3 lobes of the minimal area surface $`(S_1+S_2+S_3)`$ inside the quark trajectories $`z^{(i)}(t)`$ and the string-junction trajectory $`z^{(Y)}(t)`$. We have also denoted $`\sigma _{\mu \nu }=\frac{1}{4i}(\gamma _\mu \gamma _\nu \gamma _\nu \gamma _\mu )`$ and everywhere Euclidean space-time is used (till the last moment when the resulting Hamiltonian is obtained in Minkowski space-time) with $`\gamma `$-matrices

$$\gamma _4=\gamma _0\beta ;\gamma _i=i\beta \alpha _i,\gamma _\mu \gamma _\nu +\gamma _\nu \gamma _\mu =2\delta _{\mu \nu }.$$

Note also that notation $`tr_Y`$ means

$$tr_YP\frac{1}{6}e_{abc}e_{a^{}b^{}c^{}}P_{abc/a^{}b^{}c^{}}.$$
(5)
In the combination $`F(k)d\rho (k)`$ in (3) one can write

$$F_{\mu \nu }\sigma _{\mu \nu }^{(i)}=\left(\begin{array}{cc}𝝈^{(i)}𝐁,\hfill & 𝝈^{(i)}𝐄\hfill \\ 𝝈^{(i)}𝐄,\hfill & 𝝈^{(i)}𝐁\hfill \end{array}\right),i=1,2,3.$$
(6)
As it was shown in  the spin-independent part of (3) which obtains neglecting the $`\mathrm{\Sigma }`$ term in (4), yields at large quark separations, $`|𝐳^{(i)}𝐳^{(Y)}|T_g,i=1,2,3,`$ the familiar area-law asymptotics

$$W_3=\mathrm{exp}[\sigma (S_1+S_2+S_3)]$$
(7)
implying linear confinement for each quark. In what follows we shall use the general expression (1) to derive the spin-dependent part of interaction both for heavy quarks (expansion in inverse powers of mass) and for light quarks.

## 3 Spin-dependent interaction in $`1/m`$ expansion

To illustrate the method we shall start with the derivation of spin-dependent (SD) forces via $`1/m`$ expansion. Defining the SD potential as $`V_{\mathrm{SD}}`$, one can write $`G_{3q}e^{TV_{\mathrm{SD}}}1TV_{\mathrm{SD}}`$, and for $`V_{\mathrm{SD}}`$ the following general form will be obtained below, similar (but not identical) to the corresponding form for heavy quarkonia ,

$$V_{\mathrm{SD}}(𝐑^{(1)},𝐑^{(2)},𝐑^{(3)})=\underset{i=1}{\overset{3}{}}\frac{𝝈^{(i)}𝐋^{(i)}}{2m_i^2}\left(\frac{1}{R^{(i)}}\frac{dV_1}{dR^{(i)}}+\frac{1}{2R^{(i)}}\frac{d\epsilon }{dR^{(i)}}\right)$$

$$+\frac{1}{N_c1}\underset{i<j}{}\frac{(𝝈^{(i)}𝐋^{(j)}+𝝈^{(j)}𝐋^{(i)})}{2m_im_j}\frac{1}{R^{(j)}}\frac{dV_2(R^{(i)},R^{(j)})}{dR^{(j)}}+$$

$$+\underset{i<j}{}\left[\frac{(𝝈^{(i)}𝝈^{(j)})V_4(R_{ij})}{12m_im_j(N_c1)}+\frac{3(𝝈^{(i)}𝐧)(𝝈^{(j)}𝐧)(𝝈^{(i)}𝝈^{(j)})}{12m_im_j(N_c1)}V_3(R_{ij})\right]+V_5,$$
(8)
where $`𝐧=\frac{𝐑_{ij}}{R_{ij}}`$,$`𝐑_{ij}=𝐑_i𝐑_j`$. We assume that current quark masses are large, $`m_i\sqrt{\sigma },i=1,2,3,`$ and hence also $`\mu _i`$ are large, since the latter are defined through $`m_i`$ and $`\sigma `$ in the stationary point analysis  and always satisfy $`\mu _i=m_i+O(1/m_i)`$. Hence for simplicity we keep in the following $`\mu _i=m_i\sqrt{\sigma }`$ and expand in inverse powers of $`1/m_i`$.

As it was observed in  the SD terms of the lowest order $`(1/m_i^2,\frac{1}{m_im_j})`$ come from 3 different sources:

A) Diagonal terms in (6) are kept together with diagonal terms in $`\mathrm{\Lambda }_i(m_i+\frac{1}{2}\gamma _\mu \dot{z}_\mu )`$, yielding one power of $`1/m_i`$. An additional power of $`\frac{1}{m_i}`$ or $`\frac{1}{m_j}`$ then comes from the expansion of $`W_3`$. This yields spin-orbit terms $`V_1^{},V_2^{}`$ and $`V_5`$.

B) The off-diagonal terms are kept both in (6) and in $`\mathrm{\Lambda }_i`$. This gives spin-orbit potential $`\frac{d\epsilon }{dR}`$.

C) Diagonal terms from two matrices (6) with $`ij`$ are retained. This yields spin-spin potentials $`V_3`$ and $`V_4`$. We now calculate the SD contributions from A) – C) point by point.

A) From (3), (4) one gets for $`i=1`$

$$tr_Y[(1+\frac{g}{2m_1}\sigma _k^{(1)}_0^TB_k(z^{(1)},t^{(1)})𝑑t_1)W_3]1TV_{\mathrm{SD}}^{(1)}$$
(9)
Using the relation $`igF_{\mu \nu }W=\frac{\delta }{\delta \sigma _{\mu \nu }(z)}W`$ which obtains easily with nonabelian Stokes representation for $`W`$, one has

$$tr_YF_{\lambda \sigma }(x,z_0)W(C)=tr_Y\{ig𝑑s_{\mu \nu }(z)F_{\mu \nu }(z,z_0)F_{\lambda \sigma }(x,z_0)W_3(C)\}$$
(10)
one can rewrite the l.h.s. of (9) as

$$tr_YW+\frac{ig^2}{2m_1}\sigma _k^{(1)}_0^T𝑑t_1tr_YB_k(z^{(1)},t^{(1)})_{S(C^{})}𝑑s_{\mu \nu }(u)F_{\mu \nu }(u,x_0)W_3(C^{}).$$
(11)
In (10) and (11) the common reference point $`z_0`$ is chosen to make both expressions gauge invariant; as it will be seen, this point will not appear in the final equations.

In (11) the contour $`C^{}`$ is deformed due to orbital momentum of quarks as compared to the zeroth-order contour $`C_0`$ consisting of straight lines. This is essential since otherwise the vacuum average $`B_kW_3(C_0)`$ vanishes since it is odd with respect to reflection $`z_iz_i`$ $`ik`$. Therefore all nonzero contribution in (11) is due to deflection of the quark path in $`C^{}`$ from the straight line in $`C_0`$.

At this point we shall describe the quark trajectory $`z_\mu ^{(i)}(t)`$ and the corresponding string piece $`W_\mu ^{(i)}`$ from the quark position to the string junction (which we for simplicity take at the origin).

$$w_\mu ^{(i)}(t,\beta )=z_\mu ^{(i)}(t)\beta ,1\beta 0$$
(12)
$$ds_{ik}^{(i)}=d\beta ^{(i)}dte_{ikm}\frac{\beta L_m^{(i)}}{im_i}$$
(13)
where $`L_m^{(i)}`$ is the (Minkowskian) angular momentum of the $`i`$-th quark

$$L_s^{(i)}=im_ie_{skm}R_k^{(i)}\dot{z}_m^{(i)},𝐑^{(i)}=𝐳^{(i)}𝐳^{(Y)}=𝐳^{(i)}.$$
(14)
Similarly $`d\sigma _{k4}^{(i)}=R_k^{(i)}d\beta ^{(i)}du_4`$, and one arrives at the result

$$B_k(z^{(i)},t_i)W_3(C^{})=ig𝑑\beta ^{(i)}𝑑u_4\frac{\beta ^{(i)}L_n^{(i)}}{im_i}B_kB_n(u_4,\beta )W+$$

$$+ig𝑑\beta ^{(i)}𝑑u_4R_l^{(i)}(u_4)B_kE_l(u_4,\beta ^{(i)})W.$$
(15)
Denoting

$$B_k(z^{(1)},t_1)E_i(𝐮,u_4)W_3e_{kin}(u_nz_n^{(1)})\frac{\mathrm{\Lambda }_0}{u_4}$$
(16)
one obtains

$$𝝈^{(i)}𝐋^{(i)}\left(\frac{1}{R}\frac{dV_1}{dR}\right)^{(i)}=g^2_0^1\beta d\beta _0^Tdu_4\frac{\sigma _k^{(i)}L_n^{(i)}}{W_3}B_k(z^{(i)},t_i)B_n(\beta z^i,u_4)$$

$$\frac{𝝈^{(i)}𝐋^{(i)}}{W_3}_0^T𝑑u_4(u_4t_1)𝑑\beta \frac{\mathrm{\Lambda }_0}{u_4}$$
(17)
where we have used the relation $`u_n(u_4)z_n^{(1)}(t_1)\dot{z}_n^{(1)}(u_4)(u_4t_1)`$.

Until now we have not used the Gaussian dominance of the vacuum, i.e. the fact that $`W_3`$ is saturated by the lowest cumulant $`FF`$, which is found to be an accurate approximation . Using it one can write $`FFW_3FFW_3`$, and introduce scalar functions $`D,D_1`$ for tensor $`FF`$, as it was done in .

Referring the reader to the Appendix for the corresponding relations, one finally obtains

$$\left(\frac{1}{R}\frac{dV_1}{dR}\right)^{(i)}=_0^R\frac{d\lambda }{R}(1\frac{\lambda }{R})_{\mathrm{}}^{\mathrm{}}𝑑\nu [D(\lambda ,\nu )+D_1(\lambda ,\nu )+\lambda ^2\frac{D_1}{\lambda ^2}]$$

$$_{\mathrm{}}^{\mathrm{}}\nu ^2𝑑\nu _0^R\frac{d\lambda }{R}\frac{D_1}{\lambda ^2}.$$
(18)
Till now we have taken into account the interaction of the spin of the $`i`$-th quark with the surface $`S_i`$, which yields the term $`(\frac{1}{R}V_1^{})^{(i)}`$ multiplied with $`(m_i^2)^1`$. At this point we consider the interaction of the $`i`$-th quark spin with the (deformed) surface $`S_j`$, which will give the term $`V_2^{}`$ in (8). For this one needs to consider a vacuum average of two $`F`$’s from two different surfaces $`S_i`$ and $`S_j`$.

In general we have for two $`F`$’s transported to the same point $`x`$ $`(\alpha ,\mathrm{}\eta `$ are fundamental color indices)

$$F(u,x)_{\alpha \xi }F(v,x)_{\gamma \eta }=\frac{trFF}{N_c^21}(\delta _{\alpha \eta }\delta _{\xi \gamma }\frac{1}{N_c}\delta _{\alpha \xi }\delta _{\gamma \eta }).$$
(19)
Taking into account (5) and the relation

$$tr_Y\mathrm{\Phi }_{\alpha \alpha ^{}}(x,y)\mathrm{\Phi }_{\beta \beta ^{}}(x,y)\mathrm{\Phi }_{\gamma \gamma ^{}}(x,y)1$$
(20)
one obtains

$$tr_YF_{\alpha \alpha ^{}}(u,x)F_{\beta \beta ^{}}(v,x)=\frac{trF(u,x)F(v,x)}{N_c(N_c1)}$$
(21)
where we have also accounted for different orientation of plaquettes in $`S_i`$ and $`S_j`$.

Now proceeding as in (15) one has

$$B_k^{(i)}W_3(C_j^{})=g_0^1\beta ^{(j)}𝑑\beta ^{(j)}𝑑u_4^{(j)}\frac{L_n^{(j)}}{m_j}B_k^{(i)}B_n^{(j)}(u_4^{(j)},\beta ^{(j)}R^{(j)})W_3+$$

$$+ig𝑑\beta ^{(j)}𝑑u_4^{(j)}R_l^{(j)}(u_4^{(j)})B_k^{(i)}E_l^{(j)}W_3.$$
(22)
At this point one can use Gaussian dominance and relations (21) to obtain finally

$$\left(\frac{1}{R}V_2^{}(R)\right)_j=_0^1𝑑\beta ^{(j)}\beta ^{(j)}_{\mathrm{}}^{\mathrm{}}𝑑\nu [D(r^{(ij)},\nu )+D_1+((r^{(ij)})^2+\nu ^2)\frac{D_1}{(r^{(ij)})^2}].$$
(23)
In a similar way one obtains from the first term on the r.h.s. of (22),with the use of the last term on the r.h.s. of (A1).

$$V_5=\underset{i>j}{}\frac{(𝝈^{(i)}𝐫^{(ij)})(𝐋^{(j)}𝐫^{(ij)})}{2m_im_j(N_c1)}\beta ^{(j)}𝑑\beta ^{(j)}𝑑\nu \frac{D_1(r^{(ij)},\nu )}{(r^{(ij)})^2}.$$
(24)
Here we have defined $`𝐫^{(ij)}=𝐑^{(i)}\beta ^{(j)}𝐑^{(j)}`$, and one should take into account, that $`D,D_1`$ depend on their arguments as

$$D(r,\nu )=D(\sqrt{r^2+\nu ^2}).$$

This concludes derivation of terms with the procedure A) and one goes over to the next point.

B) Following (9) one can write for the corresponding term of a given quark ($`i`$)

$$V_{SD}^{(\epsilon )}T=\frac{g}{(2m_i)^2}\left(\begin{array}{cc}m_i+\mu _i\hfill & 𝝈^{(i)}𝐩^{(i)}\hfill \\ 𝝈^{(i)}𝐩^{(i)}\hfill & m_i\mu _i\hfill \end{array}\right)\left(\begin{array}{cc}0\hfill & 𝝈^{(i)}𝐄\hfill \\ 𝝈^{(i)}𝐄\hfill & 0\hfill \end{array}\right)W_3(C).$$
(25)
Now one can use relation

$$E_k(z^{(i)},t^{(i)})W_3(C_0)=\frac{\delta W_3(C_0)}{ig\delta \sigma _{k4}(z^{(i)},t^{(i)})}=\frac{1}{ig}\frac{\epsilon }{R_k^{(i)}}W_3(C_0)$$
(26)
where the following notation was introduced for the spin-independent potential $`\epsilon (R^{(1)},R^{(2)},R^{(3)})`$.

$$W_3(C_0)=\mathrm{exp}[\epsilon (R^{(1)},R^{(2)},R^{(3)})T].$$
(27)
Note that in (25-27) one can keep in $`W_3(C)`$ the unperturbed (straight-line) contours for quark trajectories since the prefactor in (25) is already $`O(1/m^2)`$. Keeping in mind relation

$$\sigma _k^{(i)}E_k^{(i)}\sigma _l^{(i)}p_l^{(i)}=E_k^{(i)}p_k^{(i)}+ie_{kln}\sigma _n^{(i)}E_k^{(i)}p_l^{(i)}$$
(28)
one recovers the second term on the r.h.s. of (8).

C) Here one considers spin-spin interaction and the corresponding term looks like

$$1V_{\mathrm{SD}}T=tr_Y\underset{i>j}{}(1+\frac{g}{2m_i}𝝈^{(i)}𝐁^{(i)}(z^{(i)},t_i)dt_i)\times $$

$$\times (1+\frac{g}{2m_j}𝝈^{(j)}𝐁^{(j)}(z^{(j)},t_j)dt_j)W_3$$
(29)
Identifying spin-spin terms in (29) and using (21) and relations for $`BBW_3`$ in Appendix one arrives at

$$V_{\mathrm{SD}}^{(\sigma \sigma )}=\underset{i<j}{}_{\mathrm{}}^{\mathrm{}}\frac{d\nu \sigma _k^{(i)}\sigma _k^{}^{(j)}}{4m_im_j(N_c1)}[\delta _{kk^{}}(D+D_1+(𝐮)^2\frac{D_1}{(𝐮)^2})$$

$$u_ku_k^{}\frac{D_1}{(u)^2}]$$
(30)
where notations are used

$$𝐮=𝐑^{(i)}𝐑^{(j)},\nu =t_it_j.$$
(31)
Rewriting (30) as

$$V_{\mathrm{SD}}^{(\sigma \sigma )}=\underset{i<j}{}\frac{𝝈^{(i)}𝝈^{(j)}V_4(u)+S_{ij}V_3(u)}{12m_im_j(N_c1)}$$
(32)
one has

$$V_4(u)=_{\mathrm{}}^{\mathrm{}}𝑑\nu (3D(u,\nu )+3D_1(u,\nu )+2u^2\frac{D_1}{u^2})$$
(33)
$$V_3(u)=_{\mathrm{}}^{\mathrm{}}𝑑\nu u^2\frac{D_1(u,\nu )}{u^2},$$
(34)
$$S_{ij}=3(𝝈^{(i)}𝐧)(𝝈^{(j)}𝐧)𝝈^{(i)}𝝈^{(j)},𝐧=\frac{𝐮}{|𝐮|}.$$
(35)
This concludes definition of all NP spin-dependent terms in (8) to the order $`O(1/m^2)`$ and in the approximation when only lowest, $`FF`$, correlator is retained in the Wilson loop.

Now comparing our expressions for $`V_1^{},V_2^{},\epsilon ^{},V_3,V_4`$ with the corresponding ones for heavy $`Q\overline{Q}`$ case, given in , one can see that they coincide exactly, the only difference being that one should sum up over all 3 quarks for $`V^{}`$ and $`\epsilon ^{}`$, and take a double sum, $`i<j`$, for $`V_2^{},V_3,V_4`$. In addition there is a term $`V_5`$ which is of 3 body character and vanishes in two-body situation since in that situation $`𝐋^{(j)}𝐫^{(ij)}=𝐋^{(j)}𝐫0.`$

## 4 Perturbative spin-dependent forces

We are now in position to consider also perturbative contributions to the SD potentials, which were calculated in . The easiest way for us is to remember, that to the lowest order, $`O(\alpha _s)`$, all perturbative terms are pair-wise interaction of quarks, and they can be reconstructed from the expressions obtained above, Eqs. (18), (23), (24),(32)-(34), using $`O(\alpha _s)`$ contribution to $`D_1(x)`$, while $`D(x)`$ does not have contributions at this order,

$$D(x)=D^{NP}(x),D_1(x)=\frac{16\alpha _s}{3\pi x^4}+D_1^{NP}(x).$$
(36)
It is rewarding to realize that $`D_1`$ does not enter into $`V_1^{}`$ (terms containing $`D_1`$ in (18) cancel exactly), so that perturbative contribution occurs only in the nondiagonal, $`ij`$, terms in (8) and in $`\epsilon ^{}`$.

In the $`Q\overline{Q}`$ case analogous calculations have been done and compared to standard ones in .

We start with the $`\frac{d\epsilon }{dR}`$ term and rewrite it in the original form (26), not assuming that $`\epsilon `$ depends on $`R^{(i)}`$ only, but also on $`𝐑^{(i)}𝐑^{(j)}`$, as it is for the Coulomb term to be added to $`\epsilon `$. In this way one obtains from (25), (26),(28)

$$V_{\mathrm{SD}}^{(\epsilon )}=\frac{e_{kln}p_l^{(i)}\sigma _n^{(i)}}{4m_i^2}\frac{\epsilon }{R_k^{(i)}}$$
(37)
and for $`\epsilon \epsilon +V_{coul},V_{coul}=\frac{2\alpha _s}{3}_{i>j}\frac{1}{|R^{(i)}R^{(j)}|}`$ one obtains

$$V_{\mathrm{SD}}^{(\epsilon ,pert)}=\frac{2\alpha _s}{3}\underset{i>j}{}\left[\frac{(𝐑^{(ij)}\times 𝐩^{(i)})𝝈^{(i)}}{4m_i^2(R^{(ij)})^3}+\frac{(𝐑^{(ji)}\times 𝐩^{(j)})𝝈^{(j)}}{4m_j^2(R^{(ij)})^3}\right].$$
(38)
This expression coincides with the corresponding in . Consider now nondiagonal spin-orbit term, equivalent to $`V_2^{}`$. Instead of using replacement (36) and doing integrations in (23), we start from more general expression (9) to derive

$$V_{\mathrm{SD}}^{(2,pert)}=\frac{2\alpha _s}{3(N_c1)}\underset{i>j}{}\frac{[𝝈^{(i)}(𝐑^{(ij)}\times 𝐩^{(j)})+𝝈^{(j)}(𝐑^{(ji)}\times 𝐩^{(i)})]}{m_im_j(R^{(ij)})^3}$$
(39)
which is again in agreement with .

Next we consider the spin-spin interaction. Here it is straightforward to replace in $`V_4`$ (33) $`D_1`$ as in (36) $`D_1D_1^{(pert)}=\frac{16\alpha _s}{3\pi x^4}`$ to obtain

$$V_4^{(pert)}(r)=_{\mathrm{}}^{\mathrm{}}𝑑\nu (3D_1^{(pert)}(r,\nu )+2r^2\frac{D_1^{(pert)}}{r^2})=$$

$$=\frac{8\alpha _s}{3r}\left(\frac{^2}{r^2}+\frac{2}{r}\frac{}{r}\right)\frac{1}{r}=\frac{32\pi \alpha _s}{3}\delta ^{(3)}(𝐫).$$
(40)
In a similar way one obtains for $`V_3^{(pert)}`$,

$$V_3^{(pert)}(r)=\frac{4\alpha _s}{r^3}.$$
(41)
One can also persuade oneself that $`V_5`$ has no pertubative counterpart since there $`𝐋^{(j)}𝐋^{(ij)}`$ and it is orthogonal to $`𝐫^{(ij)}`$ and therefore the total perturbative SD contribution to the order $`O(\alpha _s)`$ can be written as

$$V_{\mathrm{SD}}^{(pert)}=V_{\mathrm{SD}}^{(\epsilon ,pert)}+V_{\mathrm{SD}}^{(2,pert)}+\underset{i<j}{}\frac{𝝈^{(i)}𝝈^{(j)}V_4^{(pert)}(R^{(ij)})+S_{ij}V_3^{(pert)}(R^{(ij)})}{12m_im_j(N_c1)}$$
(42)
where explicit form of 4 terms on the r.h.s. of (42) is given in (38), (39), (40) and (41).

Our results for $`V_{\mathrm{SD}}^{(2,pert)},V_3^{(pert)},V_4^{(pert)}`$ coincide with the corresponding expressions in , however our $`V_{\mathrm{SD}}^{(\epsilon ,pert)}`$ is 2 times smaller than the corresponding term in . In the next sections we shall argue that for light quark this term gets indeed twice as big, since there one should replace $`m_i\mu _i`$, and for $`V_{SD}^{(\epsilon )}`$ the coefficient appears to be 2 times larger.

## 5 Spin-dependent forces for light quarks

In sections 3, 4 the SD forces have been obtained as an expansion in $`1/m_i,1/m_j`$ taking all 3 quark current masses large, $`m_i\sqrt{\sigma },i=1,2,3.`$

It was noticed before  however, that general expressions (1)-(4) for Green’s functions written in FSR, with the einbein function $`\mu _i(t)`$ introduced as in , allow to obtain expressions for SD forces also for light quarks without $`1/m`$ expansion, and the corresponding terms for the meson case have been written before . Below we demonstrate in this section that the same procedure works also for the $`3q`$ case with light current quarks masses as well.

We start again with general form (3) and instead of expansion in $`1/m`$ (or $`1/\mu `$ which is equivalent for heavy quark) shall do the only approximation: keeping in the sum in the exponent (3) the lowest (Gaussian) cumulant $`F(1)F(2)`$. This approximation was recently supported by lattice data for Casimir scaling , while higher cumulants provide (for Wilson loop) less than 2% .

We start with the spin-spin interaction, which is easily obtained keeping in (3) the bilocal term and in (4) only the $`\sigma `$\- dependent term. In (3) one should take into account that $`F(1)`$ and $`F(2)`$ belong to different lopes $`S_1,S_2`$ of the $`S_{123}`$ surface, hence $`n=1`$ for each of them; moreover one uses (19) and (21), the latter with opposite sign, since orientation of $`S(1,3)`$ and $`S(2,3)`$ is the same in our case. As a result one obtains

$$V_{\mathrm{SD}}^{(\sigma \sigma )}=\underset{i>j}{}_{\mathrm{}}^{\mathrm{}}\frac{d(t_it_j)}{4\mu _i\mu _j(N_c1)}\sigma _{\mu _i\nu _i}^{(i)}\sigma _{\mu _j\nu _j}^{(j)}F_{\mu _i\nu _i}(i)F_{\mu _j\nu _j}(j).$$
(43)
The combination $`\sigma ^{(i)}\sigma ^{(j)}`$ in (43) is a product of two $`4\times 4`$ matrices, which can be split into the product of Pauli spin matrices $`\sigma _i`$ and chiral $`2\times 2`$ matrices $`\widehat{1}`$ and $`\widehat{\rho }_1\left(\begin{array}{cc}0\hfill & 1\hfill \\ 1\hfill & 0\hfill \end{array}\right).`$ Thus one can rewrite (43) as

$$V_{\mathrm{SD}}^{(\sigma \sigma )}=\underset{i>j}{}_{\mathrm{}}^{\mathrm{}}\frac{g^2d(t_it_j)}{4\mu _i\mu _j(N_c1)}\sigma _m^{(i)}\sigma _n^{(j)}[trB_m(i)B_n(j)(\widehat{1}\times \widehat{1})+$$

$$+trE_m(i)E_n(j)(\widehat{\rho }_1\times \widehat{\rho }_1)+trB_m(i)E_n(j)(\widehat{1}\times \widehat{\rho }_1)+$$

$$trE_m(i)B_n(j)(\widehat{\rho }_1\times \widehat{1})]V_{\mathrm{SD}}^{(\sigma \sigma )}(BB)+V_{\mathrm{SD}}^{(\sigma \sigma )}(EE)+$$

$$V_{\mathrm{SD}}^{(\sigma \sigma )}(BE)+V_{\mathrm{SD}}^{(\sigma \sigma )}(EB).$$
(44)
Using formulas from Appendix for correlators of $`B,E`$ one has

$$V_{\mathrm{SD}}^{(\sigma \sigma )}(B_1B)=\underset{i>j}{}\frac{𝝈^{(i)}𝝈^{(j)}V_4(u)+S_{ij}V_3(u)}{12\mu _i\mu _j(N_c1)}(\widehat{1}\times \widehat{1})$$
(45)
where $`𝐮𝐑^{(i)}𝐑^{(j)}`$.

One can see that (45) coincides with (32-35) with substitution $`m_i,m_j\mu _i,\mu _j`$. For $`EE`$ term one obtains

$$V_{\mathrm{SD}}^{(\sigma \sigma )}(EE)=\underset{i>j}{}\frac{𝝈^{(i)}𝝈^{(j)}\stackrel{~}{V}_4(u)+S_{ij}\stackrel{~}{V}_3(u)}{12\mu _i\mu _j(N_C1)}(\widehat{\rho }_1\times \widehat{\rho }_1)$$
(46)
where we have defined

$$\stackrel{~}{V}_4(u)=_{\mathrm{}}^{\mathrm{}}𝑑\nu (3D(u,\nu )+3D_1(u,\nu )+(3\nu ^2+u^2)\frac{D_1(u,\nu )}{\nu ^2})$$
(47)
$$\stackrel{~}{V}_3(u)=_{\mathrm{}}^{\mathrm{}}𝑑\nu u^2\frac{D_1(u,\nu )}{u^2}=V_3(u).$$
(48)
Finally for the last two terms in (44) one has

$$V_{\mathrm{SD}}^{(\sigma \sigma )}(BE)=V_{\mathrm{SD}}^{(\sigma \sigma )}(EB)=\underset{i>j}{}\frac{(𝝈^{(i)}\times 𝝈^{(j)})\frac{1}{2i}(\frac{𝐩^{(i)}}{\mu _i}+\frac{𝐩^{(j)}}{\mu _j})}{4\mu _i\mu _j(N_c1)}_{\mathrm{}}^{\mathrm{}}\frac{D_1(u,\nu )}{u^2}\nu ^2𝑑\nu .$$
(49)
This concludes calculation of spin-spin interaction.

We turn now to the calculation of spin-orbit terms. The corresponding expression in (3) can be written as

$$W_3\mathrm{exp}(g\sigma F)_{so}=\mathrm{exp}\{\underset{i=1}{\overset{3}{}}i\frac{dt_i}{2\mu _i}\sigma _{\mu \nu }^{(i)}ds_{\rho \sigma }^{(i)}(u)D_{\mu \nu ,\rho \sigma }(z,u)+$$

$$+\frac{i}{N_c1}\underset{ij}{}\frac{dt_i}{2\mu _i}ds_{\rho \sigma }^{(j)}(u)\sigma _{\mu \nu }^{(i)}D_{\mu \nu ,\rho \sigma }(z,u)\}$$
(50)
In (50) we have defined as in Appendix

$$D_{\mu \nu ,\rho \sigma }(z,u)=\frac{g^2}{N_c}trF_{\mu \nu }(z(t_i))F_{\rho \sigma }(u)=D(h)(\delta _{\mu \rho }\delta _{\nu \sigma }\delta _{\mu \sigma }\delta _{\nu \rho })$$

$$+\frac{1}{2}[_\mu h_\rho \delta _{\nu \sigma }+perm.]D_1(h)$$

and $`h_\mu =z_\mu (t_i)u_\mu `$.

The Dirac structure of the exponent in (50) is a sum, which can be written with notations from (44) as

$$V^{(so)}=\underset{i=1}{\overset{3}{}}V_i^{(so,diag)}(\widehat{1}_i\times \widehat{1}_{jk})+V_i^{(so,nondiag)}(\widehat{\rho }_1^{(i)}\times \widehat{1}_{jk}),ij,k.$$
(51)
In the $`4\times 4`$ matrix $`\sigma _{\mu \nu }^{(i)}`$ we first consider the diagonal part. Repeating all the steps leading to (18) and (23), (24) one has the same expressions with the replacement $`m_i\mu _i,i=1,2,3,`$ namely

$$V_{\mathrm{SD}}^{(sodiag)}=\left\{\underset{i=1}{\overset{3}{}}\frac{𝝈^{(i)}𝐋^{(i)}}{2\mu _i^2}\frac{1}{R^{(i)}}\frac{dV_1}{dR^{(i)}}+\frac{1}{N_c1}\underset{i<j}{}\frac{𝝈^{(i)}𝐋^{(i)}+𝝈^{(j)}𝐋^{(i)}}{2\mu _i\mu _j}\frac{1}{R^{(j)}}\frac{dV_2}{dR^{(j)}}\right\}.$$
(52)
Let us now consider the nondiagonal part in (51), which can be written as

$$W_3\mathrm{exp}(g\sigma F)_{(so,nondiag)}=\mathrm{exp}\{i\underset{i=1}{\overset{3}{}}\frac{dt_i}{2\mu _i}\sigma _k^{(i)}(D_{k4,l4}^{(ii)}(z,u)ds_{l4}^{(i)}(u)$$

$$+D_{k4,mn}(zu)ds_{nm}^{(i)}(u))+i_{ij}\frac{dt_i}{2\mu _i}\sigma _k^{(i)}(D_{k4,l4}^{(ij)}(z,u)ds_{l4}^{(j)}(u)+D_{k4nm}^{(ij)}ds_{nm}^{(j)}(u))_{ij}$$
(53)
Now taking into account (13), (14) one has writing (53) in the form

$$W_3\mathrm{exp}(g\sigma F)_{(so,nondiag)}=\mathrm{exp}[T\{\underset{i=1}{\overset{3}{}}V_{so,nondiag)}^{(ii)}+\underset{i<j=1}{\overset{3}{}}V_{(so,nondiag)}^{(ij)}]\}$$
(54)
one finds from (53), replacing $`D_{\mu \nu ,\rho \sigma }`$ from (A1) and using (13), (14)

$$V_{(so,nondiag)}^{(ii)}=\mathrm{\Delta }_{EE}^{(ii)}+\stackrel{~}{\mathrm{\Delta }}_{EE}^{(ii)}+\mathrm{\Delta }_{EB}^{(ii)}$$
(55)
and similarly for $`V_{(so,nondiag)}^{(ij)}`$. For terms on the r.h.s. of (55) one obtains

$$\mathrm{\Delta }_{EE}^{(ii)}=i\frac{𝝈𝐑^{(i)}}{2\mu _iR^{(i)}}\mathrm{\Lambda }^{(ii)},$$

$$\mathrm{\Lambda }^{(ii)}=_0^{R^{(i)}}𝑑\nu 𝑑u(D(\nu ,u)+D_1+\nu ^2\frac{D_1}{\nu ^2})$$

$$\stackrel{~}{\mathrm{\Delta }}_{EE}^{(ii)}=i(𝝈𝐮)(𝐑^{(i)}𝐮)\frac{D_1}{u^2}(\nu ,u)𝑑\nu 𝑑\beta R^{(i)},𝐮=𝐑^{(i)}\beta ;$$

$$\mathrm{\Delta }_{EB}^{(ii)}=\frac{i}{2\mu _i^3}_{\mathrm{}}^{\mathrm{}}\nu ^2𝑑\nu \frac{D_1(\nu _1,h)}{\nu ^2}_0^1\beta ^{(i)}𝑑\beta ^{(i)}(𝐋^{(i)}\times 𝝈^{(i)})𝐩^{(i)}.$$
(56)
Similarly for $`\mathrm{\Delta }^{(ij)}`$ one obtains from the second term on the r.h.s. of (53)

$$\mathrm{\Delta }_{EE}^{(ii)}=i\frac{𝝈𝐑^{(j)}}{2\mu _iR^{(j)}}\frac{\mathrm{\Lambda }^{(ij)}}{(N_c1)},$$

$$\mathrm{\Lambda }^{(ij)}=R^{(i)}_0^1𝑑\beta ^{(j)}_{\mathrm{}}^{\mathrm{}}𝑑\nu (D(\nu ,r^{(ij)})+D_1+\nu ^2\frac{D_1}{\nu ^2})$$

$$\stackrel{~}{\mathrm{\Delta }}_{EE}^{(ii)}=i\frac{(𝝈𝐫^{(ij)})(𝐑^{(j)}𝐫^{ij)})}{2\mu _i(N_c1)}𝑑\nu _0^1𝑑\beta ^{(j)}\frac{D_1(\nu ,r^{(ij)})}{(r^{(ij)})^2}$$

$$\mathrm{\Delta }_{EB}^{(ij)}=\frac{i}{2\mu _i\mu _j}_{\mathrm{}}^{\mathrm{}}\nu ^2𝑑\nu \frac{D_1(\nu ,r^{(ij)})}{\nu ^2}_0^1\beta ^{(j)}𝑑\beta ^{(j)}(𝐋^{(j)}\times 𝝈^{(j)})(\frac{𝐩^{(i)}}{\mu _i}+\beta ^{(j)}\frac{𝐩^{(j)}}{\mu _j}).$$
(57)
Here we have defined $`𝐫^{(ij)}=𝐑^{(i)}\beta ^{(j)}𝐑^{(j)}`$.

## 6 Discussion

.

Let us now discuss the results obtained in the paper. For the heavy-quark case the nonperturbative dependent potential is given in (8), and the perturbative part in (42), so that the total SD potential is

$$V_{\mathrm{SD}}^{(total)}=V_{\mathrm{SD}}^{(nonpert)}+V_{\mathrm{SD}}^{(pert)}.$$
(58)
Perturbative part agrees with that obtained long ago in  and repeated in many subsequent papers. However in  the term $`V^{(\epsilon ,pert)}`$ is taken twice as big as (38) (or corresponding term in ). A possible modification for light quarks, which can produce this increase is discussed later.

The nonperturbative part $`V_{\mathrm{SD}}^{(nonpert)}`$ (8) consists of six terms, which were never fully written before.

Only asymptotics at large distances of the first term in (8) has been written before in  and later in , one can find it from (8)

$$V^{(nonpert)}(R^{(i)}\mathrm{})=\frac{𝝈^{(i)}𝐋^{(i)}\sigma }{4m_i^2R^{(i)}}.$$
(59)
In  were instead postulated the pairwise nonperturbative spin-orbit forces, which contradict expressions derived in this paper and in . All other terms, proportional to $`V_2^{},V_3,V_4`$ and $`V_5`$ have never been written for the 3q case, while for the $`q\overline{q}`$ case the corresponding terms (except for $`V_5`$) have been written in , . The term $`V_5`$ which has no counterpart in the $`q\overline{q}`$ case, is completely new, and its physical implication is still unclear.

We now turn to be light quark case. Here the total SD ”potential” is in general a sum of product of $`(4\times 4)`$ matrices, which can be written as

$$\widehat{V}_{\mathrm{SD}}^{(lightquarks)}=\underset{i,j}{}(V_{diag}^{(ij)}\widehat{1}_i\times \widehat{1}_j+\widehat{V}_{nondiag.}^{(ij)})$$
(60)
where $`\widehat{V}_{nondiag}^{(ij)}`$ contains terms like $`\widehat{1}_i\times \widehat{\rho }_{1j},\widehat{\rho }_{1i}\times \widehat{1}_j,\widehat{\rho }_{1i}\times \widehat{\rho }_{1j}`$, and $`\rho _1=\left(\begin{array}{cc}0\hfill & \widehat{1}\hfill \\ \widehat{1}\hfill & 0\hfill \end{array}\right)`$, where each entry in $`\rho _1`$, is $`2\times 2`$ unit matrix.

Now for $`V_{diag}^{(ij)}`$ one has

$$V_{diag}^{(ij)}=V_{\mathrm{SD}}^{(\sigma \sigma )}(BB)+V_{\mathrm{SD}}^{(so,diag)}$$
(61)
where the first term on r.h.s. of (61) is given in (45) and the second in (52). One can see in these expressions for $`V_{diag}^{(ij)}`$ the same terms as in (8) with exchange $`m_i\mu _i`$ except for the spin-orbit term proportional to $`\frac{d\epsilon }{dR}`$. Before discussing two different strategies for obtaining this last term, let us look at the general structure of (60). It has the described above matrix form and depends on einbein fields $`\mu _i,i=1,2,3`$. The latter have been defined previously in  as scalars, $`2\mu _i=\frac{dz_i(\tau )}{d\tau }`$, and are assumed to be found from the stationary point equation in the path-integral form of the meson Green’s function, or from the stationary point of the Hamiltonian. Now the spin-independent part of Hamiltonian is a unit matrix and hence can produce scalar stationary values for $`\mu _i`$. The situation changes however if one tries to incorporate also the SD part of Hamiltonian in the stationary point equation for $`\mu `$, since it would require $`\mu `$ to have a matrix form similar to that of $`\widehat{V}_{\mathrm{SD}}`$.

This is possible in the generalized form of the FFFR, which is now under investigation, but in the present form the only possible way of treatment the SD part of Hamiltonian is to consider it as a perturbation. For light quarks it is not an expansion in $`1/\mu _i`$, and the whole expression (60) is obtained with the only and numerically good approximation – keeping the bilocal (Gaussian) correlator, neglecting all higher ones.

As it was shown in meson case , this perturbation procedure works well even for lowest mesons, where SD corrections produce up to around 15% of the total mass (a similar situation holds true in earlier quark model calculations with fixed and prescribed constituent masses, see e.g. ). For heavier meson and baryon states the masses $`\mu _i`$ grow rapidly with quantum numbers  and validity of perturbative treatment of SD terms becomes even better established. In what follows we describe a perturbative procedure of treating the nondiagonal terms.

To this end we must remember (as in point B of derivation in section 4) that nondiagonal terms are also present in the preexponential factor $`(m\widehat{D})`$ in (1). Consider the largest nondiagonal term $`\mathrm{\Delta }_{EE}^{(ii)}`$ in (56), and take for simplicity its asymptotic form

$$V_{EE}i\underset{i=1}{\overset{3}{}}\frac{𝜶^{(i)}𝐧^{(i)}\sigma }{2\mu _i},𝐧^{(i)}=𝐑^{(i)}/R^{(i)}.$$
(62)
One has (omitting index $`(i)`$ for simplicity)

$$(m\widehat{D})\mathrm{exp}(V_{EE}T)\left(\begin{array}{cc}m+\mu ,\hfill & 𝝈𝐩\hfill \\ 𝝈𝐩,\hfill & m=\mu \hfill \end{array}\right)(1+i\frac{\sigma }{2\mu }𝜶𝐧+\mathrm{}T).$$
(63)
Comparing with the leading term, given by the upper left corner, one normalizes $`(m\widehat{D})`$ by extracting the factor $`(m+\mu )`$ and thus obtains

$$Eq.(63)=(1V_{\mathrm{SD}}^{(\epsilon )}T+\mathrm{}),V_{\mathrm{SD}}^{(\epsilon )}=\underset{i=1}{\overset{3}{}}\frac{\sigma 𝝈^{(i)}𝐋^{(i)}}{2\mu _i(m_i+\mu _i)R^{(i)}}.$$
(64)
In the heavy quark limit, $`\mu _im_i`$ and (61) coincides with the term proportional to $`\frac{1}{R}\frac{d\epsilon }{dR}`$ in (8). For light quarks , when $`\mu _im_i`$ however one has twice as large coefficient in(61) which coincides with the heavy-quark expansion with the light quark expression. Thus our total expression for SD potential treated as perturbation is a $`2\times 2`$ matrix

$$\stackrel{~}{V}_{\mathrm{SD}}^{(lightquarks)}=\underset{i,j}{}V_{diag}^{(ij)}+V_{\mathrm{SD}}^{(\epsilon )}$$
(65)
where $`V_{diag}^{(ij)}`$ is given in (61), (45), (52) and $`V_{\mathrm{SD}}^{(\epsilon )}`$ is given in (64), with general form obtained by replacing $`\sigma \frac{\epsilon }{R}`$. Now one can see from these expressions that we have a full correspondence between terms in (8) and in (65), where each term in (65) is obtained from the corresponding one in (8) by replacement $`m_i\mu _i`$, except for the term with $`\frac{d\epsilon }{dR}`$, where one replaces $`2m_i^2\mu _i(\mu _i+m_i)`$.

## 7 Concluding remarks

We have obtained all perturbative and nonperturbative spin-dependent terms in the 3q system in the approximation when lowest (bilocal) field correlator is retained in Wilson loop. The analogous procedure for mesons in  yielded SD potentials satisfying Gromes relation , with correct asymptotics at large distances of Thomas precession type.For the 3q system we also get this asymptotics for spin-orbit terms in the form of a sum of one-body Thomas terms, in agreement with earlier results in . All other nonperturbative terms and exact nonasymptotic form of Thomas terms are new. The signs of perturbative and nonperturbative spin-orbit terms are different and one may expect some cancellation, which should be checked in exact calculations of baryon spectra with spin splittings. All nonperturbative SD terms in (8) except for $`V_5`$ have the structure similar to that of the $`Q\overline{Q}`$ case,considered in ,except that spin-orbit terms are of one-body rather than the two-body character.The new term $`V_5`$ (24) does not have a $`Q\overline{Q}`$ analog,and after averaging over coordinates has a structure similar to two-body spin-orbit force.

The large $`N_c`$ structure of SD interaction can be clearly seen from explicit expressions and may be represented as leading $`(O(N_c^0))`$ terms of one-body spin-orbit interaction, when both fields in the field correlator are on the same sheet of the 3-sheet surface, and suppressed $`(O(N_c^1))`$ terms of spin-spin interactions and spin-orbit from two different sheets. Hence the 3q dynamics in the large $`N_c`$ limit reduces to the uncorrelated motion of $`N_c`$ quarks around a common center (string junction), which can be taken as infinitely heavy.

The general structure of the SD potential (8) at large $`N_c`$ is in agreement with the classification done in  where the unsuppressed at large $`N_c`$ terms are one-body spin-orbit potentials,while two-body spin-dependent terms are $`1/N_c`$ suppressed.In addition in  appear also spin-flavour terms which can be associated with with pion and kaon exchange forces. The latter were not considered in the present paper,but can be easily included in the same formalism,using the new chiral Lagrangian derived in .It is shown there that in the $`qstring\overline{q}`$ system pions are emitted by quarks with the known amplitude,so that the pion-exchange force can be predicted unambiguosly and added to the those obtained in the present work. This would complete the overall picture of SD forces in baryon. The author is grateful to J.Goity for useful discussion,remarks and suggestions.

This work was supported by DOE contract DE-ACOS-84ER 40150 under which SURA operates the Thomas Jefferson National Accelerator Facility.

Appendix A

Field correlators

¿From general definitions of $`g^2F_{\mu \nu }F_{\lambda \sigma }`$ through $`D,D_1`$ in  one gets

$$\frac{g^2}{N_c}trF_{\mu \nu }(x)\mathrm{\Phi }(x,y)F_{\rho \sigma }(y)\mathrm{\Phi }(y,x)=(\delta _{\mu \rho }\delta _{\nu \sigma }\delta _{\mu \sigma }\delta _{\nu \rho })D(z)+$$

$$+\frac{1}{2}[_\mu z_\rho \delta _{\nu \rho }+perm.]D_1(z)$$
(A.1)
$$\frac{g^2}{N_c}trB_i(x)\mathrm{\Phi }(x,y)B_j(y)\mathrm{\Phi }(y,x)=(\delta _{ij}(D(z)+D_1(z)+𝐳^2\frac{D_1}{z^2}]$$

$$z_iz_j\frac{D_1}{z^2}$$
(A.2)
$$\frac{g^2}{N_c}trE_i(x)\mathrm{\Phi }(x,y)E_j(y)\mathrm{\Phi }(y,x)=(\delta _{ij}(D(z)+D_1(z)+z_4^2\frac{D_1}{z^2}]+$$

$$+z_iz_j\frac{D_1}{z^2}$$
(A.3)
$$\frac{g^2}{N_c}trB_i(x)\mathrm{\Phi }(x,y)E_j(y)\mathrm{\Phi }(y,x)=e_{ijk}z_4z_k\frac{D_1}{z^2}$$
(A.4)
where we have defined

$$z_\mu =x_\mu y_\mu ,\mu =1,2,3,4,\mathrm{\Phi }(x,y)=P\mathrm{exp}ig_y^xA_\mu (u)𝑑u_\mu $$


text:
# Heisenberg model in a random field: phase diagram and tricritical behavior

## I Introduction

The random field Ising model (RFIM) has been investigated extensively both theoretically and experimentally in last years. -  Some questions such as the lower critical dimension and the existence of a static phase transition have been already solved from the theoretical point of view, but others questions as the existence of the tricritical point are still open. However, it is now understood that RFIM describes the essential physics of a rich class of experimentally accessible disordered systems. These include structural phase transitions in random alloys, commensurate charge-density-wave systems with impurity pinning, binary fluid mixtures in random porous media, and the melting of intercalates im layered compounds such as TiS2. The behavior of the RFIM is governed by the competing tendencies of the spins to align ferromagnetically under impetus of the exchange or to follow local fields and so be uncorrelated.

The RFIM has been studied with different methods. The renormalization-group theory (RG) of phase transition and critical phenomena represent one of most powerful tool of theoretical physics in last decades. Until then (early 1970s) all calculations of critical exponents were made either by exactly or by series expansions. Therefore, employing the renormalization-group arguments (which can be applied for weak fields), the proof of the equivalence between an Ising ferromagnet in a random field and a dilute antiferromagnet in a uniform field has been demonstrated. These question have been very explored (see III and references therein).

In the RG scheme, such as, the effective field RG (EFRG),  and the Mean-field RG (MFRG) ,  methods have been successfully employed in spin systems, where the results obtained are in accordance with other more effective approach (series expansion, Monte Carlo simulation among others). Those approaches (MFRG and EFRG) are based on comparison of two clusters of different sizes, each of them simulating infinite systems. On the other hand, the mean field approximation (MFA) has been extensively applied in the study of RFIM. The MFA provides very interesting results about RFIM. For example, a bimodal distribution for random field gives rise to a tricritical point and consequently first order phase transitions, at very low temperatures for certain values of field strength, but does not occur in the case of a gaussian distribution.- In addition, it is possible to establish a complete mapping between the parameters of the Ising ferromagnet in a random field and the dilute Ising antiferromagnet in a uniform field.

By employing the same strategy of EFRG, the effective field theory (EFT) has been used to study cooperative phenomena and phase transitions in various systems and giving useful qualitative and quantitative insights for the critical behavior of a wide variety of classical and quantum lattice spin systems.- The approach of EFT uses, as a starting point, the rigorous Callen-Suzuki spin identities (See Refs. III and III) and the effects of the surrounding spins on each cluster are taken into account by using a convenient differential operator technique introduced by Honmura and Kaneyoshi. In this procedure, all relevant self-spin correlations are taken exactly into account in the EFT equation of state. Therefore, the EFT approach is superior to the standard MFA.

In this work, we study the behavior of phase diagram (in $`Th`$ plane) for Heisenberg classical model of spin-1/2 on a square and simple cubic lattice by employing a bimodal probability distribution for random fields in a cluster with two spins by using the result recently obtained on diluted systems. The outline of the remainder of this paper is as following: in section 1, the formalism and calculations are developed and the results and conclusions are presented in section 2.

## II Formalism and Calculations

The system to be studied is the $`n`$-vectorial model in a random field described by the Hamiltonian

$$\beta =K\underset{(i,j)}{}𝐒_𝐢𝐒_𝐣+\underset{i}{}𝐡_𝐢𝐒_𝐢,$$
(1)
where the summation is performed over all pairs of the nearest-neighboring sites $`(i,j)`$. The quantities $`𝐒_𝐢`$ are isotropically interacting $`n`$ -dimensional classical spins of magnitude $`\sqrt{n}`$ localized on site $`i`$ and the Cartesian components of $`𝐒_𝐢`$ obeys the normalization condition,

$$\underset{u}{\overset{n}{}}(S_i^n)^2=n.$$

$`K`$ \[$`J/k_BT`$, $`k_B`$ is the Boltzmann constant and T the temperature\] is the exchange interaction between the spins, $`h_i`$ \[$`h\mu _BH/k_BT`$ , where $`\mu `$ is the Bohr magneton and $`H`$ is the random magnetic field\] is the reduced random magnetic field at site $`i`$ with probability distribution

$$P(h_i)=\frac{1}{2}[\delta (h_i+h)+\delta (h_ih)].$$
(2)
Hamiltonian (1) reduces to the $`S=\frac{1}{2}`$ Ising, planar (XY), Heisenberg and spherical models for $`n=1,2,3`$ and $`\mathrm{}`$, respectively.

In this work, we follow the EFT procedure (see Ref. III) to study the critical properties of the Hamiltonian described by Eq.(1) by employing the axial approximation. Since the Hamiltonian for a cluster with two spins can be written as

$$H=K𝐒_1𝐒_2+a_1S{}_{1}{}^{1}+a_2S{}_{2}{}^{1},$$
(3)
where $`a_l=h{}_{l}{}^{}+K\mathrm{\Sigma }_{jl}^{z1}S_j^1`$, ($`l=1,2`$) and $`z`$ is the lattice coordination number.
The average magnetization per-spin $`m=\frac{1}{2}\left(S_1^1+S_2^1\right)`$ for cluster with two spins is given by (See Refs. III and III)

$$m=\underset{k1,2}{\overset{z1}{}}\left(\alpha _x+S_k^1\beta _x\right)\underset{l1,2}{\overset{z1}{}}\left(\alpha _y+S_l^1\beta _y\right)g_n(X,Y)|_{X=h_1,Y=h_2}.$$
(4)
$`\alpha _\nu =\mathrm{cosh}(KD_\nu )`$, $`\beta _\nu =\mathrm{sinh}(KD_\nu )`$, $`(\nu =x,y),`$ $`X(Y)=x(y)+h_1(h_2)`$. $`D_\nu `$ ($`\frac{}{\nu }`$) is differential operator which satisfy the mathematical relation

$$sinh(aD_x+bD_y)g_n(X,Y)|_{X=h_1,Y=h_2}=g_n(a+h_1,b+h_2),$$

where $`g_n(X,Y)`$ is given by

$$g_n(X,Y)=\frac{\mathrm{sinh}(X+Y)}{\mathrm{cosh}(X+Y)+\mathrm{exp}(2K)M_n(K)\mathrm{cosh}(XY)},$$
(5)
and

$$M_n(K)=\frac{I_{n/2}(nK)I_{n/21}(nK)}{I_{n/2}(nK)+I_{n/21}(nK)}.$$

Here $`I_n(X)`$ is a modified Bessel function of the first kind. Eq. (4) is exact and will be applied here as a basis of our formalism, since it yields the cluster magnetization $`m`$ and the corresponding multi-spin correlation functions associated with various sites for the cluster under consideration. Here we apply the EFT approximation in both sides of Eq. (4), i.e., the thermal and random average (denoted by $`\mathrm{}_c`$), along with the decoupling procedure which ignores all high-order spin correlations, namely $`S_i^1S_j^1\mathrm{}S_n^1_cS_i^1_cS_j^1_c\mathrm{}S_n^1_c`$ , with $`ij\mathrm{}n.`$ Based on this approximation and replacing each boundary configurational spin average by the symmetry breaking mean-field parameters $`b_i`$ for all $`i`$, one set up the equation of state for $`\overline{m}=m_c`$. By using the properties of differential operator and assuming translational invariance, we obtain for the square (sq) and simple cubic lattice (sc), respectivelly,

$`\overline{m}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{6}{}}}A_{k,sq}(K,n,h)\overline{m}^k,`$ (6)
$`\overline{m}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{10}{}}}A_{k,sc}(K,n,h)\overline{m}^k.`$ (7)
The coefficients $`A_{k,(sq,sc)}(K,n,h)`$ \[$`k=`$ even\] are zero to satisfy the time reversal symmetry of the Ising model, as well the properties of differential operator technique.

## III Results and Conclusions

In this section we apply the conditions for determining of the second-order transition line and the tricritical point (TCP) in the $`n`$ vector model in a random field. We will focus on the case $`n=3`$, which corresponds to the Heisenberg classical model of spin-$`1/2`$.
In the vicinity of the second-order phase transition, the Eqs. (6 \- 7) is given by

$$\overline{m}^2=\frac{A_{1,(sq,sc)}(K,n,h)1}{A_{3,(sq,sc)}(K,n,h)}.$$
(8)
Since the magnetization ($`\overline{m}`$) goes to zero continuously, to obtain the second order phase transition line, in the $`Th`$ plane, we need solve the equations

$$A_{1,(sq,sc)}(K,n,h)=1,A_{3,(sq,sc)}(K,n,h)<0.$$
(9)
On the other hand, the r.h.s. of Eq. (8) must be positive; otherwise, the transition is interpreted to be of first-order. Unfortunately the first order transition line cannot to be analyzed in the frame of the differential operator technique, therefore we have confined our calculations only to the second order transitions, including the TCP. We obtain the TCP by solving the equations

$$A_{1,(sq,sc)}(K,n,h)=1,A_{3,(sq,sc)}(K,n,h)=0.$$
(10)
The numerical solution of Eq. (9) provide the second-order phase transition line which is shown in Figure 1. When the random field is zero, we obtain the well-known critical temperature of Heisenberg classical model of spin-$`1/2`$, $`K_{c,sq}^1=3.005`$ and $`K_{c,sc}^1=5.031`$. In Figure 1, the solid line indicate the second-order phase transitions; the black diamonds denote the position of the TCP at which phase transition changes from second to first order. The numerical solution of Equations (10) provide the values for the tricritical points $`(h_{TCP},K_{TCP}^1)=(1.800,0.958)`$ and $`(h_{TCP},K_{TCP}^1)=(2.785,2.389)`$ for the square and simple cubic lattice, respectivelly. We might conclude saying that EFT approximation is able to predict the presence of tricritical behavior for the Heisenberg classical model of spin-$`1/2`$ on a square and simple cubic lattice in a random magnetic field which takes on the random values $`\pm h`$ with equal probabilities. As can be seen, our results show the existence of the TCP for square lattice ($`z=4`$) in accordance with mean-field and Betthe-Peierls approximations. This occur because to obtain the first-order phase transition we must have at least the $`\overline{m}^5`$ term in the expansion of Eqs.(6 \- 7). Therefore, in the presente approach we cannot observe behavior tricritical for lattice with coordination numbers $`z<4`$.