format: text

text:
# The Antonov problem for rotating systems

## 1 Introduction

The name “Antonov problem” is usually referred to the study of the mean-field equilibrium state of self-gravitating $`N`$-body systems of classical particles in a finite volume. One considers the system with Hamiltonian

$$H_NH_N(\{𝒓_i\},\{𝒑_i\})=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2+\mathrm{\Phi }(\{𝒓_i\})$$
(1)
$$\mathrm{\Phi }(\{𝒓_i\})=G\underset{1i<jN}{}\frac{1}{|𝒓_i𝒓_j|}$$
(2)
enclosed in a finite three-dimensional spherical box of radius $`R`$, which is necessary in order to prevent particles from evaporating. In (1,2): $`𝒓_iV^3`$ and $`𝒑_i^3`$ denote, respectively, the position and momentum of the $`i`$-th particle, $`G`$ is the gravitational constant, and $`V`$ stands for the (volume of the) box containing the particles; we set the particles’ masses to one. The main problems are: (a) find the equilibrium configurations, i.e. the particle density profiles $`\rho (𝒓)`$ ($`_V\rho (𝒓)𝑑𝒓=N`$) that maximize the microcanonical entropy ($`\alpha =`$ constant)

$$S_NS_N(E)=\mathrm{log}\left[\frac{\alpha }{N!}\delta (H_NE)D𝒓D𝒑\right]$$
(3)
($`D𝒙=_{i=1}^Nd𝒙_i`$); and (b) analyze the corresponding thermodynamics. Clearly, the existence of such profiles and their actual structure will depend on the value of the total energy $`E`$. Intuitively, when $`E`$ is large and positive (that is when the kinetic term dominates) a gaseous-type of structure will be preferred, with particles uniformly spread over $`V`$, whereas when $`E`$ is sufficiently negative particles are expected to attract each other in a dense core embedded in a vapor. The situation presents however several subtleties. *Assuming that particles are point masses obeying Boltzmann statistics in space*, i.e. that they can come arbitrarily close to each other, it is found within the mean-field approximation  that

1. there is no global entropy maximum, therefore all maxima are local and correspond to metastable situations;
2. if $`\frac{RE}{GN^2}ϵ>ϵ_c0.335`$ (i.e. at sufficiently high energies, for we take $`R`$ to be fixed) there exist local entropy extrema where the particle density satisfies the condition

$$\rho (𝒓)=A\mathrm{exp}\left[\beta \frac{\rho (𝒓^{\mathbf{}})}{|𝒓𝒓^{\mathbf{}}|}𝑑𝒓^{\mathbf{}}\right]\beta =\frac{3/2}{\left[E\mathrm{\Phi }(\{𝒓_i\})\right]}$$
(4)
where $`A`$ is a constant. Notice that $`\beta `$ is itself a function of $`\rho `$. The physically admissible solutions of (4) are spherically symmetric (in fact, they can be mapped onto the solutions of the Emden equation $`y^{\prime \prime }+2y^{}/x=e^y`$, with $`y(0)=y^{}(0)=0`$) and are called microcanonical “isothermal spheres”;
3. in general, local extrema in the space of $`\rho `$’s can be maxima (metastable states) or saddle points (unstable). Microcanonical isothermal spheres are local maxima if the “density contrast” $`\rho (0)/\rho (R)709`$, otherwise they are at most saddle points of the entropy surface;
4. if $`ϵ`$ lies in the range between $`0.2`$ and $`ϵ_c`$ isothermal spheres have negative specific heat111As first pointed out by Thirring , in the microcanonical ensemble self-gravitating systems can have negative specific heat, at odds with the canonical setup, in which it is related to energy fluctuations and positive definite.: in such states the system heats up and shrinks upon decreasing energy;
5. finally, if $`ϵ<ϵ_c`$ (i.e. at sufficiently low energies) there are not even local entropy extrema.

The standard interpretation  is that, while at high energies an equilibrium configuration is achievable in the form of an isothermal sphere, below $`ϵ_c`$ the system collapses to a dense core and overheats, following the negative-specific-heat trend to its extreme consequences. Such a “transition”, occurring at the Antonov *point* $`ϵ=ϵ_c`$, is known as “gravothermal collapse” or “catastrophe” (see also  for a more formal treatment of this transition).

In this work we extend the above theory to rotating systems, by including the effects of the angular momentum $`𝑳=_{i=1}^N𝒓_i\times 𝒑_i`$ as a second conserved quantity besides $`E`$, and discuss the thermodynamics and the corresponding equilibrium states (which will now depend on both $`E`$ and $`𝑳`$). Our work is motivated mainly by the crucial role rotation is expected to play in astrophysical systems on different scales such as globular clusters and the cores of elliptic galaxies , both of which can be modeled by (1,2). (Hence our “particles” can be interpreted equally well as atoms and as stars.) A gravity-driven collapse can be seen as the prime cause for the equilibrium shaping of self-gravitating systems, both rotating and not. But the nature and the result of this collapse are expected to change significantly in the presence of rotation, and in particular for systems rotating sufficiently fast. Rotating systems were tackled in the past by different models (see e.g.  for a discussion of the stability of rotating gaseous cylinders, and the more recent works  and references therein), but a complete analysis of the Antonov problem with angular momentum is still lacking. Some progress was achieved recently for another class of systems where particles are assumed to obey a statistics that prohibits overlapping, the so-called Lynden-Bell statistics .

Indeed, we will show that the above picture (a)-(e) holds with minor modifications for slowly rotating systems, that is when $`|𝑳|=L`$ is sufficiently small. In particular, defining the reduced angular momentum $`\lambda =L/\sqrt{RGN^3}`$, we identify an Antonov *line* $`ϵ_c(\lambda )`$, with $`ϵ_c(0)=ϵ_c0.335`$ which plays in rotating systems the same role as the Antonov point plays in non-rotating ones. For $`\lambda `$ larger than a critical value $`\lambda _c0.45`$, however, the collapse occurs not to a single dense cluster but to a double cluster (a “binary star”), signaling a breaking of the (axial-)rotational symmetry of (1). It is remarkable that, strictly speaking, for $`\lambda >\lambda _c`$ there is no Antonov limit: double cluster solutions exist at all energies, albeit being locally unstable (they are saddle points of entropy in the space of $`\rho `$’s). The lower the energy, the more dense and point-like the two clusters. We calculate the global phase diagram of the model in the $`(ϵ,\lambda )`$ plane and analyze the thermodynamics by studying the phase transition that take place and the caloric curves. The equilibrium density profiles obtained in the different regions of the phase diagram show that, depending on the values of $`ϵ`$ and $`\lambda `$, the formation of “rings”, “single stars”, “double stars” and “disk”-like structures is possible.

From a broad statistical mechanics viewpoint, it is well known that self-gravitating systems, and more generally systems with long-range interaction potentials222By which we mean potentials decaying with the interparticle distance $`r`$ more slowly than $`r^{D\delta }`$ with $`\delta >0`$ in $`D`$ dimensions when $`r\mathrm{}`$., represent a subtle and delicate issue . For a recent critical discussion of the general statistical mechanics of non-trivial systems see . The main problem is that the usual thermodynamic limit $`(N\mathrm{},V\mathrm{},N/V\text{ finite})`$ for systems with potentials such as (2) does not exist . Their thermodynamics is evidently not extensive. Now extensive systems necessarily have positive specific heat, but this need not be true for non-extensive ones . Hence the canonical and microcanonical description of non-extensive systems do not necessarily coincide (they do only in “regular” phases with positive specific heat, and in the dilute limit $`(N\mathrm{},V\mathrm{},N/V^{1/3}`$ finite$`)`$ discussed at length in ). In the most general and physically interesting cases, self-gravitating systems should be studied in the microcanonical ensemble, a circumstance that obviously introduces several technical difficulties which are amplified by the introduction of rotation.

The choice of Boltzmann statistics, which is made by most authors , to describe the behaviour of particles contributes to complicate the situation further. The basic assumption is that all particles can eventually occupy the same point/cell in space, i.e. overlapping between particles is admitted. With Boltzmann statistics, *the system (1,2) has no ground state*, hence the entropy is not bounded. This is the ultimate origin of Antonov’s catastrophe, even (as we shall see) in rotating systems. Clearly, if one wants to describe a system where Newtonian gravity is the dominant macroscopic interaction then configurations with high particle density should be avoided, for at such short distances quantum effects and nuclear interactions become more important than gravity. One way around this problem is to introduce hard cores for particles. Another substantially equivalent way is to use Lynden-Bell statistics  instead of Boltzmann’s.

We shall proceed as follows. In Sec. 2, we will expound our microcanonical mean-field theory, deriving the analogous of Eq. (4) for rotating systems. In Sec. 3 we will show results from the numerical solution of the new equation, and discuss in particular the phase diagram in the $`(ϵ,\lambda )`$ plane and the caloric curves. A comparison between the Lynden-Bell scenario discussed in  and the Boltzmann theory is presented in Sec. 4. Finally, we formulate our conclusions with some final remarks .

## 2 Mean-field theory

The theoretical analysis of Antonov’s problem in a rotating framework can be carried out along the lines of the mean-field approximation described in , where the equilibrium properties of a self-gravitating and rotating system were studied assuming Lynden-Bell statistics for particles333A calculation similar to the one presented here (with Boltzmann statistics) was performed in , where however the resulting stationarity condition is studied in the context of the Thirring model , a simplified model of a rotating self-gravitating gas.. The microcanonical entropy in our case is given by

$$S_NS_N(E,𝑳)=\mathrm{log}\left[\frac{\alpha }{N!}\delta (H_NE)\delta (𝑳\underset{i=1}{\overset{N}{}}𝒓_i\times 𝒑_i)D𝒓D𝒑\right]$$
(5)
where $`\alpha `$ is a constant that makes the argument of the logarithm dimensionless. The integral over momenta, namely

$$F_N(\{𝒓_i\},K,𝑳)=\delta (K\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2)\delta (𝑳\underset{i=1}{\overset{N}{}}𝒓_i\times 𝒑_i)D𝒑$$
(6)
can be calculated via its Laplace transform (integrals resulting from insertion of the integral representation of the $`\delta `$-function are at most of Gaussian type, hence trivial )

$`\stackrel{~}{F}_N(\{𝒓_i\},s,𝑳)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}F_N(\{𝒓_i\},K,𝑳)e^{sK}𝑑K(\mathrm{}s>0)`$ (7)
$`=`$ $`(\text{det}𝕀)^{1/2}\left({\displaystyle \frac{2\pi }{s}}\right)^{\frac{3N3}{2}}e^{\frac{1}{2}s𝑳^T𝕀^1𝑳}`$ (8)
where $`𝕀𝕀(\{𝒓_i\})`$ denotes the inertia tensor, with elements ($`a,b=1,2,3`$)

$$I_{ab}(\{𝒓_i\})=\underset{i=1}{\overset{N}{}}(r_i^2\delta _{ab}r_{i,a}r_{i,b})$$
(9)
and $`𝑳^T𝕀^1𝑳=_{a,b}L_aI_{ab}^1L_b`$. Inversion of (8) gives

$$F_N(\{𝒓_i\},K,𝑳)=\{\begin{array}{cc}\frac{(2\pi )^{\frac{3N3}{2}}}{\mathrm{\Gamma }(\frac{3N3}{2})\sqrt{\text{det}𝕀}}(K\frac{1}{2}𝑳^T𝕀^1𝑳)^{\frac{3N5}{2}}\hfill & K>\frac{1}{2}𝑳^T𝕀^1𝑳\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(10)
hence one is left with

$$S_N=\mathrm{log}\left[\frac{\alpha A}{N!}\frac{[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }(\{𝒓_i\})]^{\frac{3N5}{2}}}{\sqrt{\text{det}𝕀}}D𝒓\right]$$
(11)
where $`A=(2\pi )^{\frac{3N3}{2}}/\mathrm{\Gamma }((3N3)/2)`$. If one is interested (as we are) in the large $`N`$ limit, one can keep only the terms of dominating order in $`N`$ in (11). This yields

$$S_N=\mathrm{log}\left[\frac{\alpha A}{N!}[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }(\{𝒓_i\})]^{\frac{3N}{2}}D𝒓\right]$$
(12)
and it remains to integrate over $`V^N`$.

To this aim, we partition the spherical box $`V`$ into $`K`$ identical cells labeled by the positions $`𝒓_k`$ of their centers ($`k=1,\mathrm{},K`$), and we let $`n(𝒓_k)`$ denote the number of particles inside the cell at position $`𝒓_k`$, with $`_{k=1}^Kn(𝒓_k)=N`$. In order to turn the integral into a sum over the cells, we adopt the following mean-field approximation for the potential $`\mathrm{\Phi }`$ and the components of the inertia tensor as functionals of the density profile $`\rho \rho (𝒓)=(K/V)n(𝒓)`$:

$$\mathrm{\Phi }(\{𝒓_i\})\mathrm{\Phi }[\rho ]=\frac{G}{2}\frac{\rho (𝒓)\rho (𝒓^{\mathbf{}})}{|𝒓𝒓^{\mathbf{}}|}𝑑𝒓𝑑𝒓^{\mathbf{}}$$
(13)
$$I_{ab}(\{𝒓_i\})I_{ab}[\rho ]=\rho (𝒓)\left(r^2\delta _{ab}r_ar_b\right)𝑑𝒓$$
(14)
This choice automatically excludes correlations from the theory. The microcanonical mean-field entropy can now be written as

$$S_N^{\text{mf}}=\mathrm{log}\left[\frac{\alpha A}{N!}[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[\rho ]]^{\frac{3N}{2}}P[\rho ]𝑑\rho (𝒓)\right]$$
(15)
where $`P[\rho ]`$ is the probability to observe a density profile $`\rho `$. Now $`P[\rho ]`$ is just proportional to the number of ways in which our $`N`$ particles can be distributed into the $`K`$ cells, and is thus determined by the statistical model by which one describes the behaviour of the particles. Assuming Boltzmann statistics, there is no spatial constraint for particles so that all of them can eventually occupy the same cell in $`V`$. This leads to

$$P[\rho ]\frac{N!}{n(𝒓_1)!\mathrm{}n(𝒓_K)}N!\mathrm{exp}\left[\frac{NK}{V}c(𝒓)(\mathrm{log}c(𝒓)1)𝑑𝒓\right]$$
(16)
where in obtaining the last expression we used Stirling’s approximation and defined the relative cell occupancy $`c(𝒓)=n(𝒓)/N=V\rho (𝒓)/(KN)`$. Introducing the dimensionless variable $`𝒙=𝒓/R`$ and neglecting irrelevant constants, one finally arrives at

$$S_N^{\text{mf}}=\mathrm{log}e^{N\mathrm{\Sigma }[c]}𝑑c(𝒙)$$
(17)
$$\mathrm{\Sigma }[c]=\frac{3}{2}\mathrm{log}\left[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[c]\right]\frac{1}{\mathrm{\Theta }}c(𝒙)(\mathrm{log}c(𝒙)1)𝑑𝒙$$
(18)
where $`\mathrm{\Theta }=c(𝒙)𝑑𝒙=\frac{V}{KR^3}`$.

For $`N`$ large though finite, $`S_N^{\text{mf}}`$ can be estimated by the steepest descent method, which for $`S=S_N^{\text{mf}}/N`$ yields

$$S=\underset{c(𝒙)}{\mathrm{max}}\left[\frac{3}{2}\mathrm{log}\left[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[c]\right]\frac{1}{\mathrm{\Theta }}c(𝒙)(\mathrm{log}c(𝒙)1)𝑑𝒙\right]$$
(19)
Notice that the stationarity properties of $`S_N`$ do not depend on $`\mathrm{\Theta }`$. This is most easily seen by setting $`f(𝒙)=c(𝒙)/\mathrm{\Theta }`$, which leads to

$$S=\underset{f(𝒙)}{\mathrm{max}}\left[\frac{3}{2}\mathrm{log}\left[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[f]\right]f(𝒙)(\mathrm{log}f(𝒙)1)𝑑𝒙+\mathrm{log}\mathrm{\Theta }\right]$$
(20)
We remark that $`𝕀𝕀[f]`$. In terms of $`f`$ and in units of $`GN^2/R`$ and $`NR^2`$, respectively, $`\mathrm{\Phi }`$ and $`I_{ab}`$ (see (13),(9)) become simply

$$\mathrm{\Phi }[f]=\frac{1}{2}\frac{f(𝒙)f(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙𝑑𝒙^{\mathbf{}}$$
(21)
$$I_{ab}[f]=f(𝒙)(x^2\delta _{ab}x_ax_b)𝑑𝒙$$
(22)
The stationarity condition $`\delta S/\delta f=0`$ implies the following equation for the extrema of the entropy:

$$f(𝒙)=\mathrm{exp}\left[\beta \frac{f(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙^{\mathbf{}}+\frac{\beta }{2}(𝝎\times 𝒙)^2+\mu \right]$$
(23)
where $`𝝎=𝕀^1𝑳`$ is the angular velocity and $`\mu `$ is a Lagrange multiplier that ensures the constraint $`f(𝒙)𝑑𝒙=1`$, and

$$\beta =\frac{3/2}{[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[f]]}$$
(24)
is the inverse kinetic energy. Notice that $`\beta \beta [f]`$.

Eq. (23) is the analogous of (4) for a rotating system. Unfortunately, it cannot be solved analytically. However, some simple qualitative considerations are possible. In absence of rotation, the solution of (23) is known to depend only on $`x=|𝒙|`$. The minimum expected effect of rotation is to create distortions that cannot be described by $`x`$ alone. Hence it is necessary to introduce angular variables, and it is convenient to pass to spherical coordinates $`(x,\theta ,\varphi )`$. If the system is slowly rotating and if, as we shall assume in the following, the angular momentum is parallel to the $`3`$-axis, one may neglect the $`\varphi `$-dependence and write formally the solution of (23) as a sum of Legendre polynomials weighted by functions of $`x`$ alone444When in the following we speak of axial rotational symmetry or for short of axial symmetry we refer always to the $`3`$-axis.:

$$f(𝒙)=\underset{l=0}{\overset{\mathrm{}}{}}f_l(x)P_l(\mathrm{cos}\theta )$$
(25)
This procedure is well-known , and has been recently applied to the Antonov problem for slowly rotating systems by Chavanis . The latter, however, do not display any serious qualitative difference from non rotating ones. The most interesting phenomena occur at high angular momenta, and in order to describe them it is necessary to include the dependence on $`\varphi `$. We thus write the solution of (23) as

$$f(𝒙)=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}f_{l,m}(x)Y_{lm}(\theta ,\varphi )$$
(26)
where $`Y_{lm}`$ are real spherical harmonics and $`f_{l,m}`$ are their weights (again functions of $`x`$ alone), which will have to be determined. Clearly, this expansion differs from the one in Legendre polynomials essentially for the terms with $`m0`$. We shall see that these terms are the ones that provide the relevant new physics at high angular momenta.

Using this expansion it is easy to show the following statements. The normalization $`f(𝒙)𝑑𝒙=1`$ is equivalent to

$$2\sqrt{\pi }x^2f_{0,0}(x)𝑑x=1$$
(27)
so that the $`(0,0)`$ harmonics is ensures this condition. The coordinates of the center of mass are given by

$$x_1^{CM}x_1f(𝒙)𝑑𝒙=2\sqrt{\frac{\pi }{3}}x^3f_{1,1}(x)𝑑x$$

$$x_2^{CM}x_2f(𝒙)𝑑𝒙=2\sqrt{\frac{\pi }{3}}x^3f_{1,1}(x)𝑑x$$
(28)
$$x_3^{CM}x_3f(𝒙)𝑑𝒙=2\sqrt{\frac{\pi }{3}}x^3f_{1,0}(x)𝑑x$$

so that harmonics with $`l=1`$ are responsible for fixing the center of mass. It is convenient to set the latter in the center of the spherical box ($`𝒙=\mathrm{𝟎}`$), which amounts to postulating $`f_{1,m}=0`$ for $`m\{1,0,1\}`$. The diagonal components of the inertia matrix read555We warn that indices like in $`I_{ab}`$ refer to spatial components, while indices like $`f_{l,m}`$ refer to harmonics.

$$I_{11}=\sqrt{\frac{16\pi }{9}}x^2f_{0,0}(x)𝑑x+\sqrt{\frac{4\pi }{45}}x^4f_{2,0}(x)𝑑x\sqrt{\frac{4\pi }{15}}x^4f_{2,2}(x)𝑑x$$
(29)
$$I_{22}=\sqrt{\frac{16\pi }{9}}x^2f_{0,0}(x)𝑑x+\sqrt{\frac{4\pi }{45}}x^4f_{2,0}(x)𝑑x+\sqrt{\frac{4\pi }{15}}x^4f_{2,2}(x)𝑑x$$
(30)
$$I_{33}=\sqrt{\frac{16\pi }{9}}x^2f_{0,0}(x)𝑑x\sqrt{\frac{16\pi }{45}}x^4f_{2,0}(x)𝑑x$$
(31)
while for the off-diagonal elements one gets

$$I_{12}=\sqrt{\frac{4\pi }{15}}x^4f_{2,2}(x)𝑑x$$
(32)
$$I_{13}=\sqrt{\frac{4\pi }{15}}x^4f_{2,1}(x)𝑑x$$
(33)
$$I_{23}=\sqrt{\frac{4\pi }{15}}x^4f_{2,1}(x)𝑑x$$
(34)
In particular, one has

$$I_{11}I_{22}=\sqrt{\frac{16\pi }{15}}x^4f_{2,2}(x)𝑑x$$
(35)
so that this parameter, which will play a foremost role in the following analysis, only depends on the $`(2,2)`$ harmonics.

Using spherical harmonics as a basis set, for the Newtonian potential one has :

$$\frac{1}{|𝒙𝒙^{\mathbf{}}|}=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}\frac{4\pi }{2l+1}\frac{\mathrm{min}\{x,x^{}\}^l}{\mathrm{max}\{x,x^{}\}^{l+1}}Y_{lm}(\theta ,\varphi )Y_{lm}(\theta ^{},\varphi ^{})$$
(36)
Now it is easy to see that the so far unknown functions $`f_{l,m}`$ satisfy the following integral equation (which can be found by multiplying both sides of (23) by $`Y_{lm}`$ and integrating over angular variables):

$$f_{l,m}(x)=Y_{lm}(\theta ,\varphi )e^{\beta _{l^{},m^{}}u_{l^{},m^{}}(x)Y_{l^{}m^{}}(\theta ,\varphi )+\frac{1}{2}\beta \omega ^2x^2\mathrm{sin}^2\theta }𝑑\theta 𝑑\varphi $$
(37)
with

$$u_{l,m}(x)=\frac{4\pi }{2l+1}\frac{\mathrm{min}\{x,x^{}\}^l}{\mathrm{max}\{x,x^{}\}^{l+1}}f_{l,m}(x^{})(x^{})^2𝑑x^{}$$
(38)
Despite their unfriendly look, Eqs (37,38), which should be solved together in a self-consistent way, can be dealt with numerically. Technical details are similar to those described in  in relation to the Lynden-Bell statistics case, and we will skip them here. Let it suffice to say that for the computation the series (36) was truncated at the order $`l_{\text{max}}=12`$ and that all harmonics (even and odd) up to $`l_{\text{max}}`$ were included, except those with $`l=1`$. Notice that this does not imply a truncation of the expansion (26), i.e. the $`f_{l,m}`$ exist also for orders higher than $`l_{\text{max}}`$. The effects of truncating the potential are well-known, and are discussed e.g. in . In short, when truncated at higher and higher order the series (36) differs more and more from the real potential in a narrower and narrower region around $`𝒙=𝒙^{\mathbf{}}`$, in a way that is reminiscent of the Gibbs phenomenon. As we said, the angular momentum was taken to lie along the $`3`$-direction. In this setup, solutions of (23) were calculated as functions of $`ϵ=\frac{ER}{GN^2}`$ and $`\lambda =\frac{L}{\sqrt{RGN^3}}`$; results of this analysis are given and discussed in the next section.

## 3 Results

### 3.1 Relevant quantities

The crucial quantity in the analysis that will follow is given by

$$\mathrm{Hes}(ϵ,\lambda )=\mathrm{det}\left(\begin{array}{cc}_ϵ^2S& _\lambda _ϵS\\ _ϵ_\lambda S& _\lambda ^2S\end{array}\right)=\kappa _1\kappa _2$$
(39)
where $`\kappa _{1,2}`$ are the eigenvalues of the matrix, ordered assuming $`\kappa _1>\kappa _2`$. The relevance of $`\mathrm{Hes}(ϵ,\lambda )`$ lies in the fact that it relates the thermostatistical equilibrium properties of systems in the microcanonical ensemble to the topology of the entropy surface in the space of conserved macroscopic quantities. In general, it can be shown  that pure, stable thermodynamic phases with positive heat capacity are characterized by having $`\kappa _1<0`$ (hence $`\mathrm{Hes}(ϵ,\lambda )>0`$). Here, $`S`$ as a function of both $`ϵ`$ and $`\lambda `$ is concave, and statistical ensembles (microcanonical and canonical) are equivalent. When at least $`\kappa _1>0`$, one has a phase coexistence region with negative specific heat, and statistical ensembles are inequivalent. $`S(ϵ,\lambda )`$ has a convex part that is ignored by the canonical ensemble, which considers the concave hull of $`S`$. The convex region corresponds to negative specific heat. For self-gravitating systems, (39) provides a powerful characterization of some of the relevant quantities. In particular, we shall concentrate on the following aspects of the phase diagram.

*Antonov limit.* For each fixed $`\lambda `$, the Antonov limit $`ϵ_c(\lambda )`$ is the minimum $`ϵ`$ for which a solution of (23) exists. In general, $`ϵ_c(\lambda )`$ can be *defined* by the condition

$$\underset{ϵϵ_c(\lambda )}{lim}\kappa _1=+\mathrm{}$$
(40)
(Strictly speaking, $`\kappa _1`$ can diverge only in the mean-field approximation in the $`N\mathrm{}`$ limit at fixed $`ϵ`$ and $`\lambda `$.) This is a straightforward generalization of the non-rotating case, where at the Antonov point $`ϵ_cϵ_c(0)`$ one has

$$\underset{ϵϵ_c}{lim}\frac{\beta }{ϵ}=+\mathrm{}$$
(41)
When two (or more) conserved quantities have to be taken into consideration, the above condition becomes a condition on the maximum eigenvalue of the matrix of second derivatives of $`S`$ with respect to the conserved quantities. Upon varying $`\lambda `$, the Antonov limit $`ϵ_c(\lambda )`$ describes a line in the $`(ϵ,\lambda )`$ plane. We call this line the Antonov line.

*Phase coexistence.* By virtue of the properties of (39), starting in the high-energy phase and decreasing $`ϵ`$ at fixed $`\lambda `$ one finds the beginning of a phase coexistence region when

$$\mathrm{Hes}(ϵ,\lambda )=0$$
(42)
is fulfilled. We shall denote by $`ϵ_0(\lambda )`$ the corresponding value of $`ϵ`$. In particular, at $`ϵ_0(\lambda )`$ the specific heat at constant $`\beta \omega =\frac{S}{\lambda }`$ ($`\omega =|𝝎|`$), i.e.

$$c_{\beta \omega }=\frac{\beta ^2}{\mathrm{Hes}(ϵ,\lambda )}\frac{^2S}{\lambda ^2}$$
(43)
has a pole and changes sign (from positive to negative) upon decreasing $`ϵ`$. Again, varying $`\lambda `$ the point $`ϵ_0(\lambda )`$ describes a line in the $`(ϵ,\lambda )`$ plane (with $`ϵ_0(0)0.2`$), which separates a pure thermodynamic phase from a phase coexistence region.

*Bifurcation to axially asymmetric solutions.* For sufficiently high $`\lambda `$, solutions of (23) without axial rotational symmetry (“binary stars”) bifurcate (continuously) from axially-symmetric solutions . This bifurcation can be detected with the aim of the order parameter (see (35))

$$D=|I_{11}I_{22}|=\left|\sqrt{\frac{16\pi }{15}}x^4f_{2,2}(x)𝑑x\right|$$
(44)
introduced in . It is easy to understand that solutions with $`D=0`$ are axially-symmetric, while solutions with $`D0`$ break axial rotational symmetry. If $`\lambda =0`$ the system is necessarily isotropic ($`I_{11}=I_{22}=I_{33}`$) and spherical symmetry cannot be broken. This reflects the fact that the solution of (23) depends only on $`x`$, so that the series expansion (26) is truncated after the first term, which involves only the $`Y_{00}`$ harmonics. When $`\lambda 0`$, anisotropies may occur ($`I_{33}I_{11},I_{22}`$) and one can have either solutions that are rotationally-symmetric around the $`3`$-axis ($`I_{11}=I_{22}`$ or $`D=0`$), which eventually depend on $`x`$ and $`\theta `$ to account for obvious rotation-induced distortions, or solutions that are asymmetric around the $`3`$-axis ($`I_{11}I_{22}`$ or $`D0`$). The latter must depend on $`x`$, $`\theta `$ and $`\varphi `$, and correspond to double clusters. However, it turns out that such solutions appear only when $`\lambda `$ is larger than a threshold value and $`ϵ`$ is sufficiently low (see also the discussion at the end of Sec. 4, and in particular Fig. 5).

*Stability of the double cluster solutions.* Finally, we shall be interested in the local stability properties of double-cluster solutions (high-$`ϵ`$, axially symmetric solutions are locally stable because in this limit the kinetic term of (1) dominates and the system behaves more and more as an ideal gas). This amounts to looking for the *marginal stability line* defined by the condition

$$\left[\frac{\delta ^2S}{\delta f(𝒙)\delta f(𝒙^{\mathbf{}})}\right]_{G[f^{}]=0}=0$$
(45)
where

$$G[f]=f(𝒙)\mathrm{exp}\left[\beta \frac{f(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙^{\mathbf{}}+\frac{\beta }{2}(𝝎\times 𝒙)^2+\mu \right]$$
(46)
and $`f^{}`$ is an axially asymmetric solution of (23). The marginal stability line marks the border between the region of the phase diagram where axially asymmetric solution are stable from that where they are at most metastable (saddle points in the entropy surface). In principle, at the most general level the local stability analysis requires solving the eigenvalue problem for the kernel $`K(𝒙,𝒙^{\mathbf{}})=\frac{\delta ^2S}{\delta f(𝒙)\delta f(𝒙^{\mathbf{}})}`$. We have not been able to carry out this task. However, it is possible to get numerical information about (45). Some technical details about this problem are given in , and we will not discuss them here.

### 3.2 Phase diagram

The phase diagram resulting from numerical solution of (23) is displayed in Fig. 1. In Fig. 2 we show instead a small sample of equilibrium density profiles. The situation for sufficiently low angular $`\lambda `$ resembles that of non-rotating systems. At high $`ϵ`$, a pure “gas” phase is found, where the equilibrium density profiles are axially-symmetric and particles are uniformly spread over $`V`$. As $`ϵ`$ is decreased, the Newtonian contribution becomes more and more important and at the $`ϵ_0(\lambda )`$ line, the specific heat becomes negative and one enters a phase coexistence region, labeled PC-I, where a collapsed configuration (see Fig. 2a) competes with the gas state (see Fig. 2b). Upon further decreasing $`ϵ`$ no more solutions of (23) are found and the Antonov limit is reached. In this region, we indeed find that

$$\underset{ϵϵ_c(\lambda )}{lim}\mathrm{Hes}(ϵ,\lambda )=\mathrm{}$$
(47)
which, comparing with (40), implies that $`\kappa _2<0`$. Rotation has the only effects of shifting both the onset of negative heat capacity and the Antonov limit to higher $`ϵ`$’s (as also argued in ), and of deforming the density profile (as clear from Fig. 2b).

When $`\lambda `$ becomes larger than the critical value $`\lambda _c0.45`$, however, the picture changes drastically. Here axially-symmetric solutions become unstable against fluctuations that break such a rotational symmetry and axially-asymmetric solutions, i.e. “double stars”, occur (see Fig. 2c). The picture is clear for high $`\lambda `$. Upon decreasing $`ϵ`$, one encounters the $`ϵ_0(\lambda )`$ line (which in that region coincides with the bifurcation line) and enters a phase coexistence region (labeled PC-II) with negative specific heat $`c_{\beta \omega }`$, Eq. (43)), where double-cluster and gas configurations compete. However, at odds with the low-$`\lambda `$ case, a further decrease in $`ϵ`$ makes double-clusters unstable. Remarkably, the solution for lower $`ϵ`$ still exists, hence there is no Antonov limit for $`\lambda >\lambda _c`$ strictly speaking. This issue will be further discussed in the next section.

The bifurcation line and the $`ϵ_0(\lambda )`$ line do not coincide in a small range close to $`\lambda _c`$ (the AB line in Fig. 1 gives the bifurcation line in this range). Here the bifurcation line separates the two mixed phases, i.e. upon crossing it by increasing $`\lambda `$ at fixed $`ϵ`$ double clusters are obtained. Transitions from PC-I to PC-II can also occur upon decreasing $`ϵ`$ at fixed $`\lambda `$.

It is important to stress that these results were obtained by truncating the potential via (36) to order $`l_{\text{max}}=12`$. While calculations performed with smaller $`l_{\text{max}}`$ in the lower half of the phase diagram ($`\lambda <\lambda _c`$) gave essentially the same results, for high angular momenta we observed some dependence on $`l_{\text{max}}`$, albeit weak. In particular, it is possible that upon increasing $`l_{\text{max}}`$ the marginal stability line is further inclined towards the $`ϵ_0(\lambda )`$ line. As a consequence, the PC-II phase would become thinner and thinner, so that double stars would be stable (if at all) only in a very tiny portion of the parameter space. This quantitative issue, which is strictly related to the use of Boltzmann statistics, certainly deserves further investigation. We believe that no other main (qualitative) changes to the picture presented here should occur.

Finally, moving on to thermodynamics, in Fig. 3 we show the caloric curves $`\beta `$ vs $`ϵ`$ at fixed low and, respectively, high $`\lambda `$. The peak in (b) corresponds to the bifurcation line.

## 4 Comparison with Lynden-Bell statistics

It is interesting to compare the behaviour of the self-gravitating and rotating gas of Boltzmann particles described above with that of a same system of Lynden Bell particles . Using Lynden-Bell statistics  instead of Boltzmann’s amounts to postulating a maximum capacity of $`n_0N`$ particles for each of the $`K`$ cells in which $`V`$ is subdivided in order to compute the integral (12), so that Eq. (16) would become

$$P[\rho ]\frac{N!}{n(𝒓_1)!\mathrm{}n(𝒓_K)!}\underset{\mathrm{cells}k}{}\frac{n_0!}{(n_0n(𝒓_k))!}$$

From the physical viewpoint, Lynden-Bell systems are leptodermous. This is not possible for Boltzmann particles, as seen by the fact that the clusters in Fig. 2 are sharply peaked, with the density decreasing smoothly from the maximum (Lynden-Bell particles can form stars with a dense, finite core and a fast-decreasing density ). Correspondingly, the stationary point condition for $`c=\mathrm{\Theta }f`$ reads

$$c(𝒙)=\left[1+\mathrm{exp}\left[\frac{\beta }{\mathrm{\Theta }}\frac{c(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙^{\mathbf{}}\frac{1}{2}\beta (𝝎\times 𝒙)^2+\mu \right]\right]^1$$

This case has been studied thoroughly in . Notice that, at odds with what happens in Boltzmann statistics (which is the $`n_0N`$ followed by $`N\mathrm{}`$ limit of the Lynden-Bell case), here the result depends on $`\mathrm{\Theta }`$, which is now given by $`\mathrm{\Theta }=\frac{NV}{n_0KR^3}`$, and is related to the total mass of the system.

The microcanonical phase diagram for this case is discussed in , where it is shown that the system has three pure phases (”gas”, single cluster and double-cluster) separated by a large phase coexistence region with negative specific heat. For historical reasons, in this paper we shall discuss the canonical phase diagram, which can be obtained by mapping the $`(ϵ,\lambda )`$ variables onto their conjugate intensive variables $`(T1/\beta ,\omega )`$. Evidently, in the canonical ensemble phase coexistence regions cannot be detected , as one must always have positive specific heat. This means that the collapse temperature is determined by the value of $`\beta `$ at which the derivative of $`\beta `$ versus $`ϵ`$ changes sign, and the region PC-I, corresponding to lower energies, is not accessible. For the present comparison purposes, however, this limited information is sufficient. Canonical phase diagrams are shown in Fig. 4.

In the Boltzmann case, one sees that for high temperatures the equilibrium state is an axially-symmetric, gas-like one. At a critical temperature $`T_c(\omega )`$ the system collapses into a single cluster (for $`\omega <\omega _c1.03`$) or into a double cluster (for $`\omega >\omega _c`$) and lower temperatures are forbidden. (For $`\omega =0`$ the known value $`T_c(0)0.3967`$ is recovered.) Lynden-Bell statistics turns out to shift the onset of both the collapse to a single cluster and the formation of the double cluster to lower temperatures. However, at lower temperatures the system can exist in a single-cluster state (low $`\omega `$) or in a double cluster state (high $`\omega `$). Notice how in both cases the phase coexistence regions with negative specific heat that are found in the microcanonical phase diagrams are lost in the canonical setup.

Nevertheless, the behaviour of the Boltzmann system beyond the marginal stability line is quite intriguing. As we said, Eq. (23) admits double-cluster type of solutions for all $`ϵ`$ if $`\lambda >\lambda _c0.45`$, so there is no Antonov limit. However, such solutions are locally unstable: there exist perturbations of the density profile which lead to an increase of entropy. A profile of the entropy density $`S`$ versus $`ϵ`$ at fixed $`\lambda =0.46`$ is shown in Fig. 5. At high energy the system is found in a stable, axially symmetric state: a “gas”. Upon decreasing $`ϵ`$, the specific heat $`c_{\beta \omega }`$ (43) stays positive until point C is reached, where it turns negative (the curvature $`\kappa _1`$ of $`S(ϵ,\lambda )`$ or equivalently the quantity $`\mathrm{Hes}(ϵ,\lambda )`$ is zero). One enters the phase coexistence region where collapsed configurations are possible. A further (small) decrease in $`ϵ`$ will drive us to the bifurcation point D, where axial-rotational symmetry breaks down in favor of such structures as “binary stars”. (The axially symmetric solution however persists, though locally unstable.) Decreasing $`ϵ`$ again one arrives at E, where double clusters become locally unstable, that is a saddle point of $`S[f]`$. Double cluster solutions (possibly more than one) however still exist left of E. In such states the two “stars” are extremely dense and almost point-like.

## 5 Conclusion and outlook

In summary, we have studied the equilibrium properties of a self-gravitating and rotating gas assuming that particles obey Boltzmann statistics in space. The physical properties of this system can be understood easily as deriving from the balance of (a) kinetic and gravitational energy at low angular momentum, and (b) gravitational and rotational energy at high angular momenta. In case (a), the effects of rotation are basically negligible except for trivial distortions of the equilibrium shapes. If the total energy is sufficiently high so that the kinetic term dominates the gravitational one, the equilibrium state of the system is axially-symmetric and gas-like. At low energies, instead, gravitation induces the system to collapse into a dense core in a low-density gas. In case (b), the high energy regime is similar but at low energy rotation is strong enough to hinder the collapse into a single cluster and drive the system to organize into a double cluster structure (a “binary star”). This simple, intuitive picture is confirmed entirely by our mean-field theory, which shows that in case (b) at least one main bifurcation (from axially symmetric to axially asymmetric) in the solutions of the basic Eq. (23) occurs, and possibly many others can be found upon varying the conserved quantities $`ϵ`$ and $`\lambda `$.

Our work leaves several directions for further research. The main question concerns probably the dynamics of these systems, which is an equally old  and widely-discussed theme in astrophysics, see e.g. . Among the many open issues and possible extensions for the equilibrium case (such as the effect of different masses, see ), we want to put forward the following few. First is the behaviour of the system left of the marginal stability line where double cluster solutions become unstable. (We mentioned in Sec. 3 that the location of this line turns out to depend on the maximum order of harmonics included in the calculation of the gravitational potential, $`l_{\text{max}}`$. Also the fact that no Antonov limit exists for $`\lambda >\lambda _c`$ could be an effect of truncation.) At low angular momenta, namely below Antonov limit, no statistical equilibrium state exists and dynamical methods should be used to grasp the physics of the system (see e.g.  and references therein). But for $`\lambda >\lambda _c`$ it is still possible to obtain (unstable) equilibrium states as saddle point of the entropy surface. This means that certain perturbations applied to such “singular” double-cluster structures (highly concentrated) will lead to an increase of entropy. It is not clear to us whether local entropy maxima (metastable states) exist in this region. Definitely, there are no global maxima for $`S`$ is everywhere unbounded. This issue certainly deserves to be studied further. Another important question that we didn’t explore concerns the relevance of mass as a further conserved quantity. If $`N`$ is fixed, it should be treated on the same level as $`ϵ`$ and $`\lambda `$, and this could lead to some qualitative changes in the picture presented here. As a final consideration, we stress the mean-field character of our solution. In spite of several claims being made of mean-field theory being exact for self-gravitating systems, we believe that it would be important to investigate the effects of particle-particle correlations, also at the merely quantitative level.

Acknowledgments. Useful discussions with O. Fliegans are gratefully acknowledged.


text:
# Center Symmetry and Abelian Projection at Finite Temperature

## Abstract

At finite temperature, there is an apparent conflict between Abelian projection and critical universality. For example, should the deconfinement transition of an $`SU(2)`$ gauge theory projected to $`U(1)`$ lie in the $`Z(2)`$ universality class of the parent $`SU(2)`$ theory or in the $`U(1)`$ universality class? I prove that the projected theory lies in the universality class of the parent gauge theory. The mechanism is shown to be non-local terms in the projected effective action involving Polyakov loops. I connect this to the recent work by Dunne et al. on the deconfinement transition in the 2+1 dimensional Georgi-Glashow model.

There is an apparent conflict between Abelian projection and critical universality. In its simplest form, critical universality states that the universality class of a phase transition depends only on the dimensionality of the system and the symmetry group of the order parameter. Universality thus tells us that the nature of the deconfinement transition for a pure gauge theory with gauge group G depends only on the dimensionality of space and the center of the gauge group C(G).

For example, $`SU(2)`$gauge theories are in the $`Z(2)`$universality class of the Ising model, while $`SU(3)`$ gauge theories are in Potts model, or $`Z(3)`$, universality class. In Abelian projection, the gauge group is reduced. This reduction can changes the center of the gauge group, and thus naively may change the universality class. A common example would be the reduction of $`SU(2)`$ to $`U(1)`$. The center of $`SU(2)`$is $`Z(2)`$, but the center of $`U(1)`$is $`U(1)`$itself. Thus the deconfinement transition of the projected theory appears to be in a different universality class from the original underlying gauge theory. However, this naive expectation is wrong: lattice simulations show that the critical exponents for $`SU(2)`$ projected to $`U(1)`$ are identical to the critical exponents of the underlying $`SU(2)`$ gauge theory . The explanation for this behavior is given below.

In lattice gauge theory, Abelian projection has a natural algorithmic formulation. For analytical purposes, I will use the formalism developed in , which uses gauge fields $`u_\mu (x)G`$, site-based gauge-fixing fields $`g(x)G`$, andAbelian gauge fields $`h_\mu (x)HG`$. Three different actions play a role: the gauge action

$$S_g=\frac{\beta }{2N}\underset{plaq}{}Tr\left(u_p+u_p^+\right),$$
(1)
the gauge-fixing action, which for $`SU(2)`$ is given by

$$S_{gf}=\lambda \underset{links}{}Tr\left[u_l\sigma _3u_l^+\sigma _3\right],$$
(2)
and the projection action

$$S_{proj}[u,h]=\underset{links}{}\left[\frac{p}{2N}Tr\left(u_l^+h_l+u_l^+g_l\right)\right]$$
(3)
As $`p,\lambda \mathrm{}`$, the usual lattice projection algorithm is formally obtained. The variables $`g`$and $`h`$ are quenched: $`u`$’s are generated by $`S_g`$, then $`g`$’s are generated by $`S_{gf}`$. Finally, $`h`$’s are generated from the gauge-fixed $`u`$’s. Expectation values are given by

$`O`$ $`=`$ $`{\displaystyle \frac{1}{Z_g}}{\displaystyle \left[du\right]e^{S_g\left[u\right]}\frac{1}{Z_{gf}[u]}\left[dg\right]e^{S_{gf}\left[\stackrel{~}{u}\right]}}`$ (4)
$`{\displaystyle \frac{1}{Z_{proj}[\stackrel{~}{u}]}}{\displaystyle \left[dh\right]e^{S_{proj}[\stackrel{~}{u},h]}O}.`$
Relations between the Green functions of the projected theory and the Green functions of the underlying non-Abelian gauge theory can be found using character expansions. The projection weight can be expanded as

$`\mathrm{exp}\left[{\displaystyle \frac{p}{2N}}Tr\left(u^+h+u^+g\right)\right]`$ (5)
$`={\displaystyle \underset{\alpha }{}}d_\alpha c_\alpha \left(p\right)\chi _\alpha \left(h^+u\right)`$
where $`\chi _\alpha `$ is a group character, $`d_\alpha `$ is the dimensionality of the representation, and $`c_a\left(p\right)`$ is a positive, p-dependent coefficient. I use $`\stackrel{~}{\chi }^\beta (h)`$ to denote the corresponding character of the projection subgroup. Note that $`\chi ^\alpha (h)`$ is in general a reducible representation of H. For example, the j=1 representation of SU(2) is the sum of three irreducible representations of U(1), given by m=+1,0,-1. To leading order in the character expansion for projection, and to order $`\lambda ^0`$ for gauge-fixing, we have

$`\stackrel{~}{\chi }^\beta (h_1..h_n)`$ (6)
$`={\displaystyle \underset{\alpha }{}}\left({\displaystyle \frac{c_\alpha (p)}{c_0(p)}}\right)^n{\displaystyle _H}(dh)\stackrel{~}{\chi }^\beta (h)\chi ^\alpha (h^+)`$
$`\chi ^\alpha (g_1..g_n)`$
for any closed Wilson loop, where $`n`$is the length of the loop. For a Polyakov loop $`u_P(\stackrel{}{x})`$, we have

$$Trh_P(\stackrel{}{x})\left(\frac{c_F(p)}{c_0(p)}\right)^{n_t}Tru_P(\stackrel{}{x})$$
(7)
Just above a second-order critical point, we have

$$Tru_P(\stackrel{}{x})\left(TT_c\right)^\beta $$
(8)
so we must have the same critical index $`\beta `$ in the projected theory as in the original. Similar results hold for higher correlation functions and thus for all critical indices. We conclude that the critical behavior of the projected theory must be that of the original gauge theory.

We can understand this result on the basis of center symmetry. A center symmetry transformation on the original gauge theory acts on all the timelike links on a given time slice as

$$u_0(\stackrel{}{x},t)zu_0(\stackrel{}{x},t)$$
(9)
where $`z`$is an element of the center C(G) of G. $`S_g\left[u\right]`$is invariant under this transformation but the Polyakov loop $`u_P\left(\stackrel{}{x}\right)`$ is not. The deconfinement transition is the spontaneous breakdown of center symmetry at high temperature. REF

We can define an effective action $`S_{eff}\left[h\right]`$ for the projected theory as

$`e^{S_{eff}\left[h\right]}`$ $`=`$ $`{\displaystyle \frac{1}{Z_g}}{\displaystyle \left[du\right]e^{S_g\left[u\right]+S_{gf}\left[u\right]+S_{proj}[u,h]}}`$ (10)
$`{\displaystyle \frac{1}{Z_{gf}[u]}}{\displaystyle \frac{1}{Z_{proj}[u]}}`$
Now suppose that w is an element of C(H). Then consider the transformation $`h_0(\stackrel{}{x},t)wh_0(\stackrel{}{x},t)`$ on a single time slice. It is easy to see that the replacement $`S_{eff}\left[h\right]S_{eff}\left[wh\right]`$ unless $`S_g\left[wu\right]=S_g\left[u\right]`$. Thus $`S_{eff}`$ is only invariant under $`C(G)`$, not $`C(H)`$.

The effective action can be constructed order by order in a strong-coupling expansion around $`\beta ,\lambda ,p=0`$. To leading order for SU(N) we have

$$S_{eff}\left[h\right]=\frac{\beta p^4}{32N^9}\underset{plaq}{}Tr\left(h_{plaq}+h_{plaq}^+\right)$$
(11)
This effective action is invariant under gauge transformations in $`H`$ and invariant under the global symmetry $`hwh`$. Note that contributions from Wilson loops of area $`A`$ are suppressed by a factor of $`\beta ^A`$, but such terms are also invariant under $`hwh`$.

Polyakov loops break the invariance of $`S_{eff}\left[h\right]`$ under $`hwh`$. The lowest such contribution comes from integration over $`N_t`$ temporal plaquettes, aligned to form a belt of plaquettes around the lattice in the temporal direction. Integration over the $`u`$ field leads to an interaction of the form

$`\left({\displaystyle \frac{\beta }{2N^2}}\right)^{N_t}\left({\displaystyle \frac{p}{2N^2}}\right)^{2N_t}`$
$`[Tr\left(h_P(\stackrel{}{x})\right)Tr\left(h_P^+(\stackrel{}{x}+\widehat{i})\right)+h.c.]`$ (12)
where $`h_P(\stackrel{}{x})`$ is a Polyakov loop formed from the $`h`$ fields:

$$h_P(\stackrel{}{x})=\underset{t}{}h_0(\stackrel{}{x},t).$$
(13)
This Polyakov loop term is not invariant under the change of variables $`h_Pwh_P`$ unless $`w`$ is in $`C(G)`$.

Consider the case of $`SU(2)`$ projected to $`U(1)`$. If we write

$$h_P(\stackrel{}{x})=\left(\begin{array}{cc}e^{i\theta (\stackrel{}{x})}\hfill & 0\hfill \\ 0\hfill & e^{i\theta (\stackrel{}{x})}\hfill \end{array}\right),$$
(14)
the interaction between Polyakov loops has the form

$$\mathrm{cos}\theta (\stackrel{}{x})\mathrm{cos}\theta (\stackrel{}{x}+\widehat{i}).$$
(15)
The transformation $`h_Pwh_P`$ induces a change $`\theta \theta +\alpha `$ which is a symmetry of the action only when $`\alpha =0`$ or $`\alpha =\pi `$, i.e., when $`wZ(2)`$. A true $`U(1)`$gauge theory, in comparison, could only have effective interactions like

$$2\mathrm{cos}\left[\theta (\stackrel{}{x})\theta (\stackrel{}{x}+\widehat{i})\right]$$
(16)
which is invariant for arbitrary $`\alpha `$. We thus see how $`S_{eff}\left[h\right]`$ includes Polyakov loop terms, non-local in Euclidean time, which force the global center symmetry to be $`C(G)`$ rather $`C(H)`$.

The need to include specific, non-local terms associated with Polyakov loops to effective actions is not limited to the case of Abelian projection. Similar behavior has recently been demonstrated by Dunne et al.  in the finite-temperature Georgi-Glashow model in $`2+1`$ dimensions. Adjoint representation scalars break an $`SU(2)`$ gauge symmetry to $`U(1)`$. The low-energy behavior of the model is that of a $`U(1)`$ theory, and thus the question of $`Z(2)`$ versus $`U(1)`$ center symmetry arises here as well.

This model has a semiclassical regime where confinement can be demonstrated as a consequence of instanton effects. Agasian and Zarembo  have shown that at finite temperature, these instantons are seen by Polyakov loops as vortices, of the type that drive Kosterlitz-Thouless phase transitions in $`2`$-dimensional systems. A detailed renormalization group analysis appeared to show that the deconfining phase transition in the $`d=2+1`$ Georgi-Glashow model was in the $`U(1)`$ universality class.

However, Dunne et al. have demonstrated that Polyakov loop effects associated with heavy $`W^\pm `$ gauge bosons must be explicitly included in any discussion of critical behavior. The effect on the finite-temperature, $`d=2+1`$ Georgi-Glashow model is precisely parallel to adding a symmetry-breaking term to a $`U(1)`$ spin model, reducing the global symmetry from $`U(1)`$ to $`Z(2)`$. It is possible to show explicitly that the critical behavior of the model lies in the $`Z(2)`$ universality class of the two-dimensional Ising model after all.

Abelian projection does not enlarge center symmetry. As a consequence, projection does not violate universality: the projected theory is in the same universality class as the original, underlying gauge theory. Any effective action for a projected theory must include non-local terms at non-zero temperatures to give the correct critical behavior at deconfinement. Comparison of projection at finite temperature with the Georgi-Glashow model at finite temperature shows the similarity of the issue in the two models, and the similar resolution. This suggests that we may have identified all the terms of the projected effective action which can be relevant in the renormalization group sense. If true, this is an important step in understanding how Abelian projection reproduces the behavior of the underlying gauge theory.


text:
# The spin angular gradient approximation in the density functional theory

## Abstract

A spin angular gradient approximation for the exchange correlation magnetic field in the density functional formalism is proposed. The usage of such corrections leads to a consistent spin dynamical approach beyond the local approximation. The proposed technique does not contain any approximations for the form of potential and can be used in modern full potential band structure methods. The obtained results indicate that the direct ’potential’ exchange in 3d magnets is rather small compared to the indirect ’kinetic’ exchange, thus justifies the dynamical aspect of the local density approximation in 3d metals.

exchange coupling, electronic structure  preprint: HEP/123-qed

In density functional theory (DFT) two magnetic terms, responsible for any spin ordering, coexist on both intraatomic and interatomic length scales and produce a very successful description of static properties in 3d magnets. One of these terms is the static non-local spin current, which is associated with the kinetic energy (’kinetic’ field (KF)), whereas another one is the large exchange-correlation magnetic field (EXF). While enjoying great success in weakly correlated materials, the standard local spin density approach (LSDA) for the exchange correlation energy is believed not appropriate for strongly correlated materials. However, the standard assumption of the LSDA - a proportionality of EXF and local magnetization $`𝐦(𝐫)`$ is, in turn appropriate only for strongly correlated systems (narrow bands), when the magnetic effects related to the ’kinetic’ exchange can be omitted. The consequence of this local approximation is that the large EXF is spin dynamically inactive in the LSDAPRB1 . This feature of LSDA restrains the possible range of directions for EXF, which immediately affects the spin-spin correlation function and, in turn a whole range of dynamic and thermal properties. Hence, it is highly desirable that both sources of internal magnetic field are treated with no constraints on their directions. A general non-local approach should be used to treat adequately the dynamic aspect of the interrelation between ’kinetic’ (indirect) and ’potential’ (direct) magnetic sources. The influence of non locality of EXF on spin dynamics (SD) of magnets is essentially unknown and was not studied in real materials. In this paper we undertake an attempt to overcame this problem and propose a spin angular gradient approximation (SAGA) for EXF and apply it to the calculations of spin wave (SW) stiffness in 3d magnetic metals.

The evolution of the many-electron inhomogeneous system in the presence of the external field is defined in a unique way by the time-dependent one-particle density matrix $`\rho _{\alpha \beta }(𝐫,𝐫,t)=\psi _\beta ^+(𝐫,t)\psi _\alpha (𝐫,t),`$ where $`\psi _\alpha (𝐫,t)`$ is the annihilation operator for the electron at the point $`𝐫`$ with spin projection $`\alpha `$ at the instant time $`t`$. Equivalently, one can introduce the charge $`n(𝐫,t)=Tr\rho (𝐫,𝐫,t)`$ and magnetization $`𝐦`$ $`(𝐫,t)=Tr\rho (𝐫,𝐫,t)\sigma `$ densities, where $`\sigma `$ are Pauli matrices. Starting from the Schroedinger equation for the many-electron system one can formally obtain the exact set of equations for these quantities:

$$\stackrel{}{𝐦}(𝐫,t)=\gamma 𝐦(𝐫,t)\times 𝐁_{ext}(𝐫,t)+\frac{i}{2}_{\mathrm{𝐫𝐫}^{}}^2(\rho _{\alpha \beta }(𝐫,𝐫^{},t)\sigma _{\beta \alpha }c.c.)_{𝐫^{}=𝐫},$$
(1)
$$\stackrel{}{n}(𝐫,t)=\frac{i}{2}_{\mathrm{𝐫𝐫}^{}}^2(\rho _{\alpha \alpha }(𝐫,𝐫^{},t)c.c.)_{𝐫^{}=𝐫}$$
(2)
where $`𝐁_{ext}(𝐫,t)`$ is the external magnetic field. These equations are not closed since, generally speaking, the quantity $`_𝐫^{}\rho _{\alpha \beta }(𝐫,𝐫^{},t)_{𝐫^{}=𝐫}`$ cannot be described directly in terms of $`\rho _{\alpha \beta }(𝐫,𝐫,t)`$ (see the discussion of the kinetic energy term in DFTPZ ). Introducing the wave functions of Kohn-Sham quasiparticles $`\phi _{\nu \alpha }`$ and making an adiabatic approximation which assumes they evolve rapidly in comparison to the spin degrees of freedom, we can obtain a closed set of equations describing charge and spin dynamics:

$`\stackrel{}{𝐦}(𝐫,t)`$ $`=`$ $`\gamma 𝐦(𝐫,t)\times 𝐁_{tot}(𝐫,t)+`$
$`{\displaystyle \frac{i}{2}}_𝐫({\displaystyle \underset{\nu }{\overset{occ}{}}}\phi _{\nu \alpha }^{}(𝐫,t)_𝐫\phi _{\nu \beta }(𝐫,t)\sigma _{\beta \alpha }c.c.),`$
$$\stackrel{}{n}(𝐫,t)=\frac{i}{2}_𝐫(\underset{\nu }{\overset{occ}{}}\phi _{\nu \alpha }^{}(𝐫,t)_𝐫\phi _{\nu \alpha }(𝐫,t)c.c.),$$
(4)
$`\left[\left({\displaystyle \frac{1}{2}}^2+V_{ext}+V_H\right)\delta _{\alpha \beta }𝐁_{ext}\sigma _{\alpha \beta }+{\displaystyle \frac{\delta E_{xc}}{\delta \rho _{\beta \alpha }}}\right]\phi _{\nu \beta }(𝐫,t)`$
$`=`$ $`ϵ_\nu \left(t\right)\phi _{\nu \alpha }(𝐫,t),`$
where the sum is over the occupied Kohn-Sham states; $`V_{ext}`$ and $`V_H`$ are the external and Hartree potentials, respectively, $`E_{xc}`$ is the exchange-correlation energy, $`𝐁_{tot}=𝐁_{ext}+𝐁_{xc}`$ is the total magnetic field including EXF $`𝐁_{xc}=\delta E_{xc}/\delta 𝐦.`$ The Kohn-Sham wave functions and the corresponding energies $`ϵ_\nu \left(t\right)`$ depend on time due to the time-dependence of the densities and external field (the latter is supposed to be weak relative to the characteristic electron energies).

The second term in Eq.(The spin angular gradient approximation in the density functional theory) can be conveniently cast in the following torque form

$`{\displaystyle \frac{d𝐦(𝐫,t)}{dt}}`$ $`=`$ $`\gamma 𝐦(𝐫,t)\times \left[𝐁_{kin}(𝐫)+𝐁_{xc}(𝐫)+𝐁_{ext}(𝐫)\right]`$
$`=`$ $`\gamma 𝐦(𝐫,t)\times 𝐁_{tot}(𝐫)`$
where the total field acting on the electronic spin at point $`𝐫`$ is

$$𝐁_{tot}(𝐫)=\frac{\delta E}{\delta 𝐦(𝐫)}=\frac{\delta T}{\delta 𝐦(𝐫)}\frac{\delta E_{xc}}{\delta 𝐦(𝐫)}+𝐁_{ext}(𝐫),$$
(7)
with $`𝐁_{kin}(𝐫)=\rho \sigma /m\rho `$, $`\rho =\varphi _\nu ^{}\varphi _\nu `$, $`d/dt=/t+𝐯`$ and $`𝐯`$ is a velocity. From Eq.(The spin angular gradient approximation in the density functional theory) directly follows that $`𝐦d𝐦/dt=0.`$Below we use $`𝐁_{ext}=0.`$

The basic assumption of LSDA is that $`E_{xc}`$ is obtained for the homogeneous electron gas (HEG) model for a given magnetization density $`m`$, and assumed to apply to real external potentials. In this local approximation $`𝐁_{xc}(𝐫)`$ $`𝐦(𝐫)`$, hence assuming very strong EXF compared to KF and smallness of dynamics associated with EXF. Let us now eliminate this requirement and generalize a shape of EXF for the spiral magnetic configuration (spin-density wave (SDW)) with the amplitude $`m`$ and the wave vector $`𝐐.`$ Hence, the term $`𝐦(𝐫,t)\times 𝐁_{xc}(𝐫)`$ will contribute to Eq.(The spin angular gradient approximation in the density functional theory).

The SDW state is characterized by anomalous averages $`s_𝐩=c_{𝐩+𝐐/2}^+c_{𝐩𝐐/2}`$ where $`c_{𝐩\sigma }^+`$,$`c_{𝐩\sigma }`$ are the operators of creation and annihilation of electrons with momentum $`𝐩`$ and spin projection $`\sigma .`$ Similar to the Gorkov-Nambu formalism in the theory of superconductivity (see, for instance, Ref.KURMAEV ) we introduce the spinor operator $`\psi _𝐩=(c_{𝐩+𝐐/2}^+,c_{𝐩𝐐/2})`$. Then the Hamiltonian can be expressed as

$$H=\underset{𝐩}{}\psi _𝐩h_𝐩\psi _𝐩+\frac{1}{2}\underset{𝐪0}{}\underset{\mathrm{𝐩𝐩}^{}}{}v_c\left(𝐪\right)\left(\psi _{𝐩+𝐪}^+\psi _𝐩\right)\left(\psi _{𝐩^{}𝐪}^+\psi _𝐩^{}\right)$$
(8)
where $`v_c\left(𝐪\right)=4\pi e^2/𝐪^2V,`$ $`V`$ is a volume and $`h_𝐩=\theta _𝐩+\tau _𝐩\sigma _z\mathrm{\Delta }_𝐩\sigma _x.`$ Here

$`\theta _𝐩`$ $`={\displaystyle \frac{1}{2}}\left(\epsilon _{𝐩+𝐐/2}+\epsilon _{𝐩𝐐/2}\right)=𝐩^2/2+𝐐^2/8\mu ,`$ (9)
$`\tau _𝐩`$ $`={\displaystyle \frac{1}{2}}\left(\epsilon _{𝐩+𝐐/2}\epsilon _{𝐩𝐐/2}\right)=\mathrm{𝐩𝐐}/2`$
where $`\epsilon _𝐩=𝐩^2/2\mu `$ is the energy of the free electron and $`2\mathrm{\Delta }_𝐩`$ is the antiferromagnetic gap, related to the formation of SDW. In the Hartree-Fock approximation (HFA) the latter is determined through the relation

$$\mathrm{\Delta }_𝐩=\underset{𝐩^{}}{}v_c(𝐩𝐩^{})s_𝐩^{}.$$
(10)
Now we will replace $`v_c`$ by the effective Stoner exchange splitting $`I=\left(V_{exc}^{}V_{exc}^{}\right)/\left(n_{}n_{}\right),`$ where $`V_{exc}^\sigma =\left(n\epsilon _{exc}\right)/n_\sigma .`$ Then, Eq.10 can be replaced by $`\mathrm{\Delta }=I\left(n_{}n_{}\right)/2`$, where $`\mathrm{\Delta }`$ does not depend on $`𝐩.`$

The ’bare’ Green function in the Matsubara representation has the form

$$G(i\omega _m,𝐩)=\frac{1}{i\omega _mh_𝐩}=\frac{i\omega _m\theta _𝐩+\tau _𝐩\sigma _z\mathrm{\Delta }_𝐩\sigma _x}{\left(i\omega _m\xi _𝐩_{}\right)\left(i\omega _m\xi _𝐩_{}\right)}$$
(11)
where $`\xi _{𝐩,}=\theta _𝐩E_𝐩`$ is a quasiparticle spectrum for SDW with $`E_𝐩=\sqrt{\tau _𝐩^2+\mathrm{\Delta }^2}`$. From Eq.11 one can find the occupation number matrix

$$2N_𝐩=\left(1+\frac{\tau _𝐩\sigma _z\mathrm{\Delta }\sigma _x}{E_𝐩}\right)f_𝐩+\left(1\frac{\tau _𝐩\sigma _z\mathrm{\Delta }\sigma _x}{E_𝐩}\right)f_𝐩$$
(12)
where $`f_{𝐩\sigma }=f\left(\xi _{𝐩\sigma }\right)`$ is a Fermi function. Then for the Fock contribution to the exchange-correlation energy we will have

$`E_{Fock}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐩𝐩}^{}}{}}v_c\left(𝐩𝐩^{}\right)Tr\left[N(𝐩)N(𝐩^{})\right]=E_{Fock}^{(1)}+E_{Fock}^{(2)},`$ (13)
$`E_{Fock}^{(1)}`$ $`={\displaystyle \frac{1}{4}}{\displaystyle \underset{\mathrm{𝐩𝐩}^{}\sigma }{}}v_c\left(𝐩𝐩^{}\right)f_{𝐩\sigma }f_{𝐩^{}\sigma }\left(1+{\displaystyle \frac{\tau _𝐩\tau _𝐩^{}+\mathrm{\Delta }^2}{E_𝐩E_𝐩^{}}}\right),`$
$`E_{Fock}^{(2)}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐩𝐩}^{}}{}}v_c\left(𝐩𝐩^{}\right)f_𝐩f_𝐩^{}\left(1{\displaystyle \frac{\tau _𝐩\tau _𝐩^{}+\mathrm{\Delta }^2}{E_𝐩E_𝐩^{}}}\right).`$
Expression (13) is suitable for the calculation of the exchange energy in SDW state for any $`𝐐`$ with the corresponding equation for the chemical potential $`\mu `$

$$N=Tr\underset{𝐩}{}N_𝐩=\underset{𝐩}{}f_{𝐩\sigma }.$$
(14)
Let us consider an important case of small $`𝐐`$. The expansion of Eq.14 gives a correction of order $`𝐐^2`$ for $`\stackrel{~}{\mu }=\mu 𝐐^2/8.`$ One has

$$\delta \stackrel{~}{\mu }=\stackrel{~}{\mu }_𝐐\stackrel{~}{\mu }_{𝐐=0}=\frac{𝐐^2}{8F(n_{},n_{})}$$
(15)
where $`F=(p_F+p_F)I(n_{},n_{})/2\pi ^2`$ is a dimensionless Stoner enhancement factor, $`p_{F\sigma }=(6\pi ^2n_\sigma )^{1/3}.`$

Expanding Eq.13 up to $`𝐐^2`$ we have the following expression

$$\frac{E_{Fock}}{V}=\frac{e^2}{8\pi ^3}\left\{\left(p_F^4+p_F^4\right)Q^2\left[\left(\frac{1}{2F}\frac{2}{3}\right)\left(p_F^2+p_F^2\right)+\frac{\left(p_F+p_F\right)^2}{12F^2}\right]\right\}.$$
(16)
Now let us consider the correlation effects. In the Gell-Mann-Brueckner theoryGMB  the correlations can be taken into account by using the polarization operator

$$\mathrm{\Pi }(i\omega ,𝐪)=Tr\underset{𝐩}{}T\underset{\epsilon _n}{}G(𝐩+𝐪,i\epsilon _n+i\omega _n)G(𝐩,i\epsilon _n).$$
(17)
With this function the correlation part of $`\mathrm{\Omega }`$-potential can be expressed as:

$$\mathrm{\Omega }_{corr}=\underset{𝐪}{}\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{d\omega }{4\pi }\left\{\mathrm{ln}\left[\frac{1+v_c\left(𝐪\right)\mathrm{\Pi }(i\omega ,𝐪)}{1+v_c\left(𝐪\right)\mathrm{\Pi }_{𝐐=0}(i\omega ,𝐪)}\right]v_c\left(𝐪\right)\left[\mathrm{\Pi }(i\omega ,𝐪)\mathrm{\Pi }_{𝐐=0}(i\omega ,𝐪)\right]\right\},$$
(18)
where we wrote only $`𝐐`$dependent part of the correlation energy. Substituting Eq.11 into Eq.17 we have

$`\mathrm{\Pi }(i\omega ,𝐪)`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐩}{}}\left(1+{\displaystyle \frac{\tau _𝐩\tau _{𝐩+𝐪}+\mathrm{\Delta }^2}{E_𝐩E_{𝐩+𝐪}}}\right){\displaystyle \underset{\sigma }{}}{\displaystyle \frac{f_{𝐩\sigma }f_{𝐩+𝐪\sigma }}{i\omega +\xi _{𝐩+𝐪\sigma }\xi _{𝐩\sigma }}}+`$ (19)
$`2\left(1{\displaystyle \frac{\tau _𝐩\tau _{𝐩+𝐪}+\mathrm{\Delta }^2}{E_𝐩E_{𝐩+𝐪}}}\right){\displaystyle \frac{f_𝐩f_{𝐩+𝐪}}{i\omega +\xi _{𝐩+𝐪}\xi _𝐩}}`$
with

$$\mathrm{\Pi }_{𝐐=0}(i\omega ,𝐪)=\underset{\sigma }{}\frac{p_{F\sigma }}{2\pi ^2}\left\{1+\frac{1}{2q}\left[p_{F\sigma }^2\left(\frac{i\omega }{q}+\frac{q}{2}\right)^2\right]\mathrm{ln}\frac{i\omega +q^2/2+qp_{F\sigma }}{i\omega +q^2/2qp_{F\sigma }}\frac{1}{2q}\left[p_{F\sigma }^2\left(\frac{i\omega }{q}\frac{q}{2}\right)^2\right]\mathrm{ln}\frac{i\omega q^2/2+qp_{F\sigma }}{i\omega q^2/2qp_{F\sigma }}\right\},$$
(20)
where $`\mathrm{ln}`$ means the main branch of the logarithm.

Now the corresponding exchange-correlation addition to the spin wave spectrum at finite $`𝐐`$ can be written as

$$\delta \omega _𝐐=\frac{4}{M}\left[E_{SDW}(𝐐)E_{SDW}(0)\right].$$
(21)
where $`E_{SDW}(𝐐)`$ is the total energy of the spin spiral and $`M`$ is a magnetic moment of the cell.

For the electronic gas of high density the main contribution into integral (18) comes from the region of small $`𝐪`$. At the same time from Eq.(19) one can see that the interband contribution contains additional factor $`𝐪^2`$ compared to the intraband one. Hence, it seems reasonable to take into account terms of $`𝐪^2`$ order for the correlation effects in intraband transitions.

Below we consider non-relativistic case only, when exchange-correlation energy is invariant with respect to the rotation in spin space. Then, for the weakly nonuniform distribution of spin density, one can write

$$E_{exc}=𝑑𝐫\left\{n\epsilon _{exc}(n_{},n_{})+\lambda (n_{},n_{})D\right\},$$
(22)
where $`D=\left(_\alpha e_\beta \right)\left(_\alpha e_\beta \right)=\left(\theta \right)^2+\mathrm{sin}^2\theta \left(\phi \right)^2`$ is the rotational invariant of lowest order. Here $`𝐞=𝐦/\left|𝐦\right|(\mathrm{sin}\theta \mathrm{cos}\phi ,\mathrm{sin}\theta \mathrm{sin}\phi ,\mathrm{cos}\theta ).`$

For the variation of $`E_{exc}`$ in the spiral structure with $`𝐞(𝐫)=(\mathrm{cos}\mathrm{𝐐𝐫},\mathrm{sin}\mathrm{𝐐𝐫},0)`$ one can write $`E_{exc}^QE_{exc}^{Q=0}=V\lambda Q^2`$. Now, one needs to estimate the energy of SDW with an amplitude $`\left|𝐦\right|=n_{}n_{}`$ and $`𝐐0`$. This problem has been considered by HerringHER  in HFA. The last one, however, seems not to be useful in the calculation of the total energy of the SDW due to an instability at $`Q2k_F`$SDW . At the same time, when screening is taken into account this instability is destroyed (see, for instance, Ref.AMUSIA ). Hence, it is essential to add screening effects to the Fock contribution to $`E_{exc}`$. So, the following expression for the function of $`\lambda `$ in Eq.22 can be obtained

$$\lambda (n_{},n_{})=\frac{e^2}{16\pi ^2}\left(\frac{1}{F}\frac{4}{3}\right)\left(V_{exc}^{}p_FV_{exc}^{}p_F\right)+\frac{e^2}{96\pi ^3F^2}.$$
(23)
This expression properly takes into account the main effect - shift of chemical potential (15) and exactly corresponds to Eq.16 in the Fock approximation.

Let us discuss the consequences of the addition (23) to exchange-correlation energy. First of all, EXF $`𝐁_{exc}(𝐫)=\delta E_{exc}/\delta 𝐦(𝐫)`$ is now noncollinear with respect to the local magnetization, so, the variation of Eq.22 leads to appearance of a new term

$$B_{exc}^\gamma (𝐫)=\frac{2}{m}\left(\delta _{\beta \gamma }e_\beta e_\gamma \right)_\alpha \left(\lambda _\alpha e^\beta (𝐫)\right).$$
(24)
The noncollinearity of $`𝐁_{exc}`$ with $`𝐞`$ corresponds to the ’direct’ exchange processes, which are absent in LSDA. This field can be directly included in any traditional full-potential LSDA technique.

Another important consequence of this non-locality is the appearance of an additional contribution for SWS of ferromagnons

$$D=\frac{4}{M}\left[\underset{𝐐0}{lim}\frac{E_{SDW}(𝐐)E_{SDW}(0)}{𝐐^2}\right].$$
(25)
For this contribution we have

$$\delta D=\frac{4}{M}𝑑r\lambda (n_{},n_{}),$$
(26)
with integration over the whole elementary cell. This correction represents an addition to the value of $`D`$ which was obtained in the random phase approximationstiff .

The numerical calculations were first performed for HEG model. The calculation of the $`𝐐`$ dependence of the correlation energy according to Eq.(16) revealed that this contribution is relatively small and does not have any peculiarities as a function of $`𝐐`$ and has monotonic behavior. In our approximation it could indicate that this term contributes to both intraatomic and interatomic range scales in similar fashion.

Further, we applied the technique described above for the studies of SD in ferromagnetic (FM) Fe and Ni. These systems are well studied and might be used as prototype systems for 3d magnetism research. Below we employ a non-collinear version of the full-potential linear augmented plane wave method (WIEN code) using both LSDA and generalized gradient correction approaches. In this approach the magnetization density is treated as a continuous vector quantity. The parameters of calculations (number of $`𝐤`$-points, plane waves and etc) were chosen to provide a convergence of SWS about 3-5%. Unfortunately, we could not compute the energy from Eq.(16) for arbitrary $`𝐐`$ with real DFT wave functions and our study of spin waves in Fe and Ni was performed only for small $`𝐐`$, i.e. for SWS from Eq.25. Our standard LSDA calculations (without SAGA) produced 239 meV/A2 and 692 meV/A2 for $`D`$ in Fe and Ni respectively (see, also Ref.Mark ). These results are in agreement with previous theoretical calculations and are close to the experimental numbers. The non-self-consistent addition of non-local correction from Eq.26 does not change significantly these results, producing a positive correction of 13 meV/A2 and 45 meV/A2 for Fe and Ni respectively (so, the total $`D`$ become 253 meV/A2 and 735 meV/A2). The addition of Eq.22 for the self-consistency procedure and the total energy calculations modifies these results. The self-consistent account of SAGA increases the total $`D`$ up to 271 meV/A$`^2,`$ with the increasing EXF contribution of 21 meV/A2. The effect of the increase of the total $`D`$ was also found in Ni where EXF contributes 94 meV/A2 to the total self-consistent value of 782 meV/A$`^2.`$ This result already indicates that SD in 3d metals exists mostly due to ’kinetic’ (indirect) exchange with relatively small contribution from ’potential’ (direct) exchange. Whereas the absolute value of the latter term is rather large, its non-local (dynamic) part is much weaker compared to ’kinetic’ spin current dynamics. Overall the addition of SAGA leads to an increase of the total $`D`$ from both the direct SAGA addition and the indirect modification of ’kinetic’ contribution (effects of self-consistency). Whereas this result is obtained only for small $`𝐐,`$ we expect that at larger $`𝐐`$ this effect will be more pronounced. Also one can argue that the addition of diagrams beyond the GW approximation can be important due to their influence on Stoner exchange splitting. Both LSDA and GW calculations produce much larger splitting than the one observed on experiment, hence one can expect that the other diagrams will correct it and as a result will generate a smaller SWS.

In summary, we estimated the non local part of the exchange correlation magnetic field in real magnets. Our results for Fe and Ni indicated that in spite of the large absolute value of the exchange correlation magnetic field, the non-local part of this field, which contributes to microscopical spin dynamics is rather small compared to the ’kinetic’ spin current term in 3d ferromagnetic metals. We expect that this correction can be important in the interstitials of weak magnets, systems with strong interaction between the local moment and conduction electrons, and magnets with non-collinear ordering.

This work was carried out at the Ames Laboratory, which is operated for the U.S.Department of Energy by Iowa State University under Contract No. W-7405-82. This work was supported by the Director for Energy Research, Office of Basic Energy Sciences of the U.S.Department of Energy.