text:
# Electronic structure of Millerite NiS

## I Introduction

There have been numerous electron spectroscopic investigations of the metastable hexagonal phase of nickel monosulphide (NiS), primarily to understand its unusual electronic phase transition . However, no electron spectroscopic studies have been undertaken on the millerite phase of NiS, which is the stable phase below 600 K. The resistivity measurements on the system show that it has a high metallic conductivity ($`2\times 10^4`$ Ohms-1-cm-1 at room temperature). It has been speculated  that the millerite phase is diamagnetic with the observed small magnetic susceptibility attributed to the presence of paramagnetic hexagonal modification of NiS. Millerite NiS crystallizes in the lower symmetry trigonal space group ($`R3m`$) with lattice parameters $`a=9.589`$ Å and $`c=3.165`$ Å . The crystal structure of millerite NiS is shown in Fig. 1(a). Here, the Ni atom has five nearest neighbor S atoms, occupying the corners of a square pyramid. In this pyramidal geometry (see Fig. 1(b)) the Ni atoms are displaced slightly out of the basal plane, towards the apical sulphur atom. Though the local coordination of Ni in millerite is very similar to that in another Ni-S system, BaNiS2, millerite has a three-dimensional structure in contrast to BaNiS2 which is highly two-dimensional in nature. The Ni-S bond in millerite NiS (shortest $`d_{NiS}=2.25`$ Å) is relatively shorter than that in hexagonal NiS ($`d_{NiS}=2.39`$ Å), NiS2 ($`d_{NiS}=2.40`$ Å) or BaNiS2 (shortest $`d_{NiS}=2.32`$ Å) and leads to stronger hybridization effects in this system.

There have been very few studies to understand the electronic structure of this compound. Band structure calculation performed for the millerite phase of NiS  showed the ground state to be metallic, in agreement with experimental results. However, no comparison between experimental and theoretical results exists and the details of the electronic structure in terms of the electronic structure parameters are still unknown.

In the present study, we investigate the electronic structure of millerite phase of NiS using x-ray photoemission (XP) and Ultra-Violet photoemission (UP) measurements in conjunction with ab inito band structure as well as parameterized many-body calculations. The present results provide a consistent and quantitative description of the electronic structure of the system. By comparing the hexagonal and millerite phases of NiS, we study the evolution of the spectral function with changing correlation effects in the metallic regime.

## II EXPERIMENTAL

For the preparation of the millerite samples, polycrystalline samples of hexagonal NiS was prepared first by solid state reaction . Hexagonal NiS sample was then sealed in quartz tube in high vacuum and was maintained at 573 K for about 2-3 weeks and was cooled slowly to room temperature over a period of 8 hours, to obtain the millerite phase. The $`x`$-ray powder diffraction measurements confirmed the phase purity of the sample. Spectroscopic measurements were carried out in a combined VSW spectrometer with a base pressure of 2$`\times 10^{10}`$ mbar equipped with a monochromatized Al K$`\alpha `$ $`x`$-ray source and a Helium discharge lamp. XP spectroscopic measurements were performed on the samples with an overall instrumental resolution of better than 0.8 eV, while UP measurements for He I and He II are performed with an instrumental resolution better than 90 meV and 120 meV, respectively. The sample surface was cleaned in situ periodically during the experiment by scraping with an alumina file and the surface cleanliness was monitored by recording the carbon 1$`s`$ and oxygen 1$`s`$ XP signal regions. The binding energy was calibrated to the instrumental Fermi-level which was determined by recording the Fermi-edge region from a clean silver sample.

## III CALCULATIONS

Scalar relativistic linearized muffin-tin orbital (LMTO) band structure calculations have been performed within the atomic sphere approximation (ASA) for millerite phase of NiS with the real crystal structure . Here, the rhombohedral unit cell consisting of 3 formula units was employed. Sphere radii used for Ni and S were 2.49 and 2.55 a.u., respectively. 17 empty spheres with sphere radii in the range of 0.9 to 1.65 a.u. were also used. Convergence was obtained with $`s`$, $`p`$ and $`d`$ orbitals at Ni and S atomic spheres and $`s`$ and $`p`$ for the empty spheres, with 220 $`k`$ points in the irreducible part of the Brillouin zone.

Core level and valence band (VB) spectra were calculated for millerite NiS, using a NiS5 cluster as found in the solid (see Fig. 1(b)), within a parameterized many-body model including orbital dependent electron-electron (multiplet) interactions. The calculational method has been described in detail elsewhere . The calculations were performed including all the transition metal 3$`d`$ and sulphur 3$`p`$ orbitals. In the VB calculation, Ni 3$`d`$ and S 3$`p`$ contributions to the valence band spectra were calculated within the same model and with the same parameter values. As the dimension of the Hamiltonian matrix is large, Lanczos method was used to evaluate the spectral function and the calculated one-electron removal spectra were appropriately broadened to simulate the experimental spectra. In the Ni 2$`p`$ core level calculation, Doniach-$`\stackrel{˘}{S}`$unji$`\stackrel{´}{c}`$ line shape function  was used for broadening the discrete energy spectrum of the cluster model, in order to represent the asymmetric line shape of core levels from these highly metallic compounds; similar asymmetric line shapes are also found in the other core levels in this system. In the case of VB calculations, energy dependent Lorentzian function was used for the lifetime broadening. Other broadening effects such as the resolution broadening and solid state effects were taken into account by convoluting the spectra with a Gaussian function. The broadening parameters were found to be consistent with values used for similar systems . Since the atomic cross-sections for the Ni 3$`d`$ and S 3$`p`$ states are vastly different, it is necessary to calculate a weighted average of these two contributions to the valence band. The atomic cross-section ratio  between S $`3p`$ and Ni $`3d`$ states ($``$ 0.17) is not appropriate in this context, since solid-state effects alter this ratio significantly . It was found that S 3$`p`$/Ni 3$`d`$ cross-section ratio of approximately 5.5 times that obtained from the atomic calculations gives the best result for the valence band calculations.

## IV RESULTS & DISCUSSIONS

The total as well as partial Ni $`d`$ and S $`p`$ densities of states (DOS) obtained from the LMTO band structure calculation are shown in Fig. 2. The thick solid line represents the total DOS, while thin solid line and dashed line show the Ni 3$`d`$ and S 3$`p`$ partial DOS, respectively. Our results are in good agreement with that of a previously published calculation . The overall features of the DOS for millerite phase is similar to that obtained for the paramagnetic hexagonal phase of NiS . However, in this case, the Fermi-level lies in the rising part of the DOS, instead of close to a minimum in DOS as in the case of hexagonal NiS where an instability in the Fermi-surface can open up a gap in certain directions of the Fermi-surface . This distinction may partially be responsible for the absence of any phase transition in the millerite phase. The DOS between -3.5 eV and 1.2 eV is dominated by the Ni 3$`d`$ contributions. Near the Fermi energy region, there is a substantial contribution from the S 3$`p`$ states, due to the strong covalency (Ni $`d`$-S $`p`$ interactions) in the system, forming the antibonding states. The DOS in the energy range -5.5 eV to -3.5 eV has dominant Ni 3$`d`$ and S 3$`p`$ contributions and represents the bonding states of the system. In the energy range, -8 eV to -5.5 eV, S 3$`p`$ contribution is dominant with a smaller contribution from the Ni 3$`d`$ states; these non-bonding states of the S 3$`p`$ are stabilized in energy compared to Ni $`d`$-S $`p`$ bonding states due to strong S-S interactions, similar to the case of hexagonal NiS and other sulphides . The peak at $`1.9`$ eV resulting from the Ni 3$`d`$ ($`t_{2g}`$-like states) is shifted to a higher energy compared to the $`t_{2g}`$ states of the hexagonal NiS. This is related to the formation of the direct Ni-Ni bonds in the millerite phase leading to the stabilization of this phase compared to the hexagonal phase .

The band structure results are compared with the experimental data in Fig. 3. The partial densities of states of Ni 3$`d`$ and S 3$`p`$ are broadened with the experimental XP spectroscopic resolution of 0.8 eV and are shown along with valence band spectra taken with 21.2 eV (He I), 40.8 eV (He II) and 1486.6 eV (Al K$`\alpha `$) photon energies. Such a comparison provides an experimental determination of the orbital characters of the various features in the experimental spectra, since the photoemission cross sections for Ni 3$`d`$ and S 3$`p`$ vary significantly with the photon energy. The decrease in the intensity of features B and C in going from 21.2 eV photon energy to 40.8 eV is explained by the Cooper minimum in the S 3$`p`$ photoionization cross section at around 50 eV, indicating the dominance of the S 3$`p`$ contributions in these two spectral features, with feature A dominated by the Ni 3$`d`$ contribution. From band structure results, the feature B has contributions from S 3$`p`$ as well as from Ni 3$`d`$, while the feature C is dominated by S 3$`p`$ states. Although the band structure calculations reproduce the features B and C rather well, the energy position of the feature A is not correctly reproduced, with the calculated peak appearing at about 0.4 eV higher binding energy compared to the experiment. This discrepancy between the experimental results and band structure calculations is attributed to the electron correlation effects within the Ni 3$`d`$ levels; thus, it appears that correlation effects continue to have important influence even for this highly metallic system. This has prompted us to go beyond the band structure theories and study the electronic structure of this system using a cluster-model, where the electron correlation effects are explicitly taken into account.

The Ni 2$`p`$ core level spectrum for the millerite phase of NiS is shown in the inset of Fig. 4 by solid circles. The spectrum consists of spin-orbit split 2$`p_{3/2}`$ (852.8 eV binding energy) and 2$`p_{1/2}`$ (870.2 eV binding energy) signals with strong satellite features at about 859.5 eV and 876 eV binding energies, corresponding to 2$`p_{3/2}`$ and 2$`p_{1/2}`$ signals, respectively. The intense satellite features in the Ni 2$`p`$ core level spectrum point to the presence of the electron correlations in the system. The 2$`p_{3/2}`$ and 2$`p_{1/2}`$ peaks show strong asymmetries, similar to hexagonal form of NiS . In order to determine the inelastic scattering background, we have performed electron energy loss spectroscopy (EELS) on these samples, with the same primary energy as that of the Ni $`2p`$ core level peak. Using a procedure that have been previously employed , the inelastic background function obtained for millerite NiS is shown in the inset of Fig. 4 as a solid line. The background function is a highly structured function with a broad plasmon loss feature around 876 eV, giving rise to an apparently stronger satellite intensity for the Ni 2$`p_{1/2}`$ peak compared to the satellite intensity accompanying the 2$`p_{3/2}`$ signal.

We have performed core and valence band calculations within a single model for a NiS5 cluster. The cluster fragment used for the calculation is shown in Fig. 1(b). Within this model, we consider only one Ni atom; effects due to the strong Ni-Ni interactions present in this compound could not be considered due to the prohibitively large basis set involved in such calculations. For Ni2+, the electron-electron interaction parameters, $`F_{dd}^2`$, $`F_{dd}^4`$, $`F_{pd}^2`$, $`G_{pd}^1`$ and $`G_{pd}^3`$ used were same as that for hexagonal form of NiS . In the main panel of Fig. 4, we show the calculated core level spectrum (solid line) including the experimentally determined inelastic scattering background, superimposed on the experimental spectrum (open circles) for the parameter set S-Ni $`(pd\sigma )`$ = -1.8 eV, $`\mathrm{\Delta }`$ = 1.0 eV and $`U_{dd}`$ = 4.0 eV. The calculated spectrum without any broadening effect is presented as a stick diagram at the bottom of the same figure. The calculated spectrum matches reasonably well with the experimental spectrum. However, there are some differences between the calculated spectrum and the experimental one; the rising edge of the Ni 2$`p_{3/2}`$ at $`852`$ eV and the satellite energy region of the 2$`p_{1/2}`$ peak are not accurately described by the calculation. Such discrepancies between the experimental and calculated spectrum may have its origin in the neglect of the strong metal-metal bonds in the system or due to the slight differences in the background function generated from the EELS spectrum and the actual background in the photoemission spectrum.

The $`(pd\sigma )`$ values obtained for millerite NiS (-1.8 eV) is considerably larger than that obtained for other Ni-S systems like, hexagonal NiS  ($`(pd\sigma )=1.4`$ eV), NiS2  ($`(pd\sigma )=1.5`$ eV) and BaNiS2  ($`(pd\sigma )=1.5`$ eV). This significant increase in the $`(pd\sigma )`$ for the millerite case is related to the shorter Ni-S bonds in the system compared to the other nickel sulphides, as mentioned earlier, leading to a high degree of covalency. The value of the charge transfer energy ($`\mathrm{\Delta }=1`$ eV) obtained for the millerite NiS is smaller than the hexagonal form of NiS  (2.5 eV) and NiS2  (2.0 eV). However, a similar value of $`\mathrm{\Delta }`$ has been observed in another Ni-S system, BaNiS2  ($`\mathrm{\Delta }=1.0`$ eV); significantly, BaNiS2 has the same local geometry around Ni atoms as that in the millerite NiS with NiS5 cluster in a pyramidal arrangement. It has been observed in the case of divalent nickel oxides  that as the dimensionality of the Ni-O connectivity decreases, the charge transfer energy decreases, via  a change in the Madelung potential. Thus, the decrease in the $`\mathrm{\Delta }`$ for millerite phase compared to hexagonal NiS is also possibly due to a change in the Madelung potential arising from the decreased local coordination of Ni atoms. The lower value of $`\mathrm{\Delta }`$ for the millerite phase compared to the hexagonal phase leads to a higher degree hybridization mixing of Ni $`d`$ and S $`p`$ states in the system, leading to an enhanced covalency. It turns out that the on-site Coulomb interaction strength, $`U`$, is very similar between millerite, hexagonal NiS  and NiS2 . This shows that $`U`$ is not very sensitive to the local co-ordination or to the metallic/insulating property within the series of nickel sulphides. According to the ZSA  phase diagram, these parameters (small $`\mathrm{\Delta }`$ and large hopping strength) place the millerite phase deep inside the $`pd`$-metallic regime.

Our calculations suggest that the ground state wavefunction of the millerite phase consists of 52%, 42.6% and 5.4% of $`d^8`$, $`d^9\underset{¯}{L}`$ and $`d^{10}\underset{¯}{L}^2`$ configurations with a high-spin $`S=1`$ state. The average $`d`$-occupancy ($`n_d`$) of the system is found to be 8.53, which is substantially larger than that obtained for hexagonal NiS ($`n_d=8.43`$). This establishes a strongly covalent ground state of the millerite phase, compared to other sulphides studied. In order to understand the origin of various features in the experimental core level spectrum, we have analyzed the characters of some typical final states, marked 1-10 in the main panel of Fig. 4, in terms of contributions from various configurations ($`d^8`$, $`d^9\underset{¯}{L}^1`$ and $`d^{10}\underset{¯}{L}^2`$) (see Table I). These features can be grouped into four regions; the main peak region, 852-855 eV (labelled 1-3 in Fig. 4), the weak satellite region between 856 eV and 858 eV (labelled 4 and 5), intense satellite region between 858 eV and 861 eV (labelled 6-8) and the high energy satellites beyond 863 eV (labelled 9 and 10). The first, third and fourth group of features have essentially similar characteristics as those in the hexagonal phase, with dominant $`d^9\underset{¯}{L}^1`$, significant contributions from all configurations, and dominant $`d^8`$ configuration, respectively . However, the second group of features, which is absent in the hexagonal phase of NiS, have almost 90% of the contributions coming from the well-screened $`d^9\underset{¯}{L}^1`$ configuration. Almost pure charge transfer nature of these states is possibly due to the shorter Ni-S distances and a lower value of $`\mathrm{\Delta }`$.

Within the same model, we have calculated the XP valence band spectrum of the millerite phase. The calculated spectrum (solid line) along with Ni 3$`d`$ (dashed line) and S 3$`p`$ (dot-dashed line) contributions for the millerite phase are shown superimposed on the experimental XP valence band spectrum (open circles) in Fig. 5. The inelastic scattering background is taken to have an integral energy dependence and is shown as a dotted line. The Ni 3$`d`$ contributions to the VB spectrum without any broadening effects are shown as a stick diagram. The parameter set used for the VB calculation is very similar to that used for the core level calculation; however, it was found that a $`(pd\sigma )`$ of -1.6 eV instead of -1.8 eV provides a better agreement with the experimental result. The agreement between the experimental and calculated spectra is reasonably good over the entire energy range. The features A and B are reproduced rather well, however, the feature C as well as the region close to the Fermi energy (lower binding energy side of the feature A) could not be described very accurately. This is however, not very surprising. The feature C is dominated by states arising from S $`p`$-S $`p`$ interactions, as shown by the band structure results (see Fig. 3). Since we do not take this interaction in to account, and also due to the intrinsic limitation of a cluster model to account for such band structure effects, the feature C is completely absent in the results of the cluster model. It is reasonable to expect the features close to the $`E_F`$ to be influenced by substantial Ni-Ni nearest-neighbor interactions present in this compound, explaining the discrepancy between the experiment and the results based on the cluster model neglecting such metal-metal interactions.

The analysis of the ground state wave function shows that the ground state has 55.2%, 40.4% and 4.4% of $`d^8`$, $`d^9\underset{¯}{L}^1`$ and $`d^{10}\underset{¯}{L}^2`$ configurations, respectively yielding an average occupancy of 8.49, with a high-spin configuration. The analysis of the final state configurations for a selected set (marked 1-10 in Fig. 5) of final states have been carried out and the results are given in table II. The various features can be grouped into different categories, namely the main peak region at 1-3 eV binding energy (labelled 1-3, corresponding to feature A), region B in the range of 4-7 eV (labelled 4-7) and weaker satellites beyond 7 eV (labelled 8-10 and corresponding to feature D). The first group of features has the dominant contribution coming from the well screened $`d^8\underset{¯}{L}^1`$ states with a non-negligible contribution from the poorly screened $`d^7`$ states and overscreened $`d^9\underset{¯}{L}^2`$ states. These contributions are found to be similar to those observed in the case of hexagonal NiS . The second group of features are predominantly contributed by $`d^8\underset{¯}{L}^1`$ and $`d^9\underset{¯}{L}^2`$ states, while in the case of hexagonal phase it was dominated by $`d^7`$ states. The third group of features have a high degree of $`d^7`$ character along with contributions from $`d^9\underset{¯}{L}^2`$ character, suggesting correlation effects within the Ni $`d`$ states. These high energy, weak intensity satellites manifest in the valence band spectrum as a weak tailing of the spectrum beyond about 7 eV (feature D in Fig. 5). It is to be noted here that the characters of these weak intensity features are nearly same as that of the features appearing between 6 and 8 eV binding energy in the VB spectrum of hexagonal NiS which have been attributed to the spectral signature of the lower Hubbard band, persisting in the $`pd`$-metallic regime. Hence, the lower Hubbard band features are shifted to higher binding energies as well as becoming weaker in intensity in going from hexagonal NiS to millerite NiS, primarily due to an enhanced hopping strength in the millerite phase driving the system deeper in to the metallic regime.

Although the presence of the Hubbard bands in the metallic regime has been predicted by the dynamical mean-field theoretical (DMFT) calculations, the evolution of these spectral features well inside the metallic regime has not been studied in detail as a function of various electronic interaction strengths. A recent DMFT calculation for the spectral functions of a charge-transfer insulator near the metal-insulator transition boundary  suggests that within the metallic regime, the spectral signatures of the Hubbard band moves towards the Fermi-level, collapsing with the coherent states, forming the metallic ground state of the system with decreasing $`\mathrm{\Delta }`$. On the other hand, our experimental results establish that well inside the $`pd`$-metallic regime, the Hubbard band is further stabilized, moving towards higher binding energy region with increasing metallicity driven by enhanced hopping interaction strength between Ni $`d`$ and S $`p`$ states. Further experimental and theoretical efforts are needed in this direction to have a better understanding of the evolution of the electronic structure of these systems.

In conclusion, the millerite phase of NiS has been studied by means of electron spectroscopic techniques, band structure and model Hamiltonian calculations. The band structure calculations were found to be more successful in describing the experimental valence band spectrum in comparison to the case of hexagonal NiS, suggesting a reduced effect of electron correlation in the millerite phase. This is consistent with the highly conducting ground state of the millerite phase, in contrast to the antiferromagnetic and poorly conducting ground state of hexagonal NiS. However, calculations including the electron correlation effects are found to be necessary for the description of certain features in the experimental spectra, indicating importance of such interaction effects even for such a highly metallic system. Thus, it appears that both band structure effects and correlation effects need to be treated on an equal footing for a complete description of such systems. The various electronic parameter strengths obtained from these calculations indicate that the millerite phase of NiS is a highly covalent metal ($`pd`$-metal). From the comparative study between hexagonal and millerite phases of NiS, the evolution of the spectral functions in a $`pd`$-metal as a function of the covalency is discussed.

## V Acknowledgments

The authors thank Professor C. N. R. Rao for continued support and the Department of Science and Technology, and the Board of Research in Nuclear Sciences, Government of India, for financial support. SRK thanks Dr. P. Mahadevan for helpful discussions. DDS thanks Dr. M. Methfessel, Dr. A. T. Paxton, and Dr. M. van Schiljgaarde for making the LMTO-ASA band structure program available. The authors also thank Professor S. Ramasesha and the Supercomputer Education and Research Center, Indian Institute of Science, for providing the computational facility.


format: text

text:
# Return mapping of phases and the analysis of the gravitational clustering hierarchy

## 1 Introduction

The morphology of the large-scale structure in the Universe is that of a complex hierarchy of nodes, filaments and sheets interlocking large voids. The Fourier-space description of such a pattern is dominated by the properties of the phases rather than the amplitudes of the Fourier modes \[Chiang 2001\]. According to the prevailing theoretical ideas this pattern developed by a process of gravitational instability from an amorphous pattern of density fluctuations characterized by a Gaussian field with random phases. Since the non-random phases of the present structure have grown from random-phase initial perturbations then there is strong motivation for understanding how phase information develops within this paradigm and to construct a statistical description of galaxy clustering that could be used as a test of the basic idea.

Unfortunately, quantifying the properties of Fourier phases is difficult for a number of technical reasons, so their use in statistical studies has so far been limited. Ryden & Gramann (1991), Soda & Suto (1992) and Jain & Bertschinger (1996) focused on the evolution of individual phases away from their initial values but since the initial phases are unknown these studies can not be used as the basis of a statistical descriptor. The pattern of association between phases is subtle and hard to visualize which makes a statistical test hard to construct a priori.

As the first step in a different approach towards quantifying phase information, Coles & Chiang (2000) proposed a colour representation method to visualize phase coupling that at least reveals qualitatively how phase information arises during the evolution of $`N`$-body experiments but does not in itself constitute a statistical descriptor. In a related study, Chiang & Coles (2000) quantified phase information using a statistic derived from the Shannon entropy of the distribution of successive phase differences. This study displayed interesting relationships between phase entropy and gravitational clustering but still did not provide a general statistical description.

In this paper we use a generalization of the concept of a return map \[May 1976, Chiang & Coles 2000\] to transfer the phases of different Fourier modes on to a bounded square upon which simple statistical tests can be applied. In this way, we build upon the earlier studies \[Chiang & Coles 2000, Coles & Chiang 2000\] to construct a method that allows us to transform the phase information in a clustering pattern into a more useful form.

## 2 Phase Coupling in the Nonlinear Regime

The mathematical description of an inhomogeneous Universe revolves around the dimensionless density contrast, $`\delta (𝒙)`$, which is obtained from the spatially-varying matter density $`\rho (𝒙)`$ via

$$\delta (𝒙)=\frac{\rho (𝒙)\rho _0}{\rho _0},$$
(1)
where $`\rho _0`$ is the global mean density. When the density perturbation is small, the evolution of the density contrast can be obtained analytically through linear perturbation theory from 3 coupled partial differential equations. They are the linearized continuity equation,

$$\frac{\delta }{t}=\frac{1}{a}\mathbf{}_x𝒗,$$
(2)
the linearized Euler equation

$$\frac{𝒗}{t}+\frac{\dot{a}}{a}𝒗=\frac{1}{\rho a}\mathbf{}_xp\frac{1}{a}\mathbf{}_x\varphi ,$$
(3)
and the linearized Poisson equation

$$\mathbf{}_{x}^{}{}_{}{}^{2}\varphi =4\pi Ga^2\rho _0\delta .$$
(4)
In these equations, $`a`$ is the expansion factor, $`p`$ is the pressure, $`\mathbf{}_x`$ denotes a derivative with respect to the comoving coordinates $`𝒙`$, $`𝒗=a\dot{𝒙}`$ is the peculiar velocity and $`\varphi (𝒙,t)`$ is the peculiar gravitational potential. From Eq. (2)-(4), and if one ignores pressure forces, it is easy to obtain an equation for the evolution of $`\delta `$:

$$\ddot{\delta }+2(\frac{\dot{a}}{a})\delta 4\pi G\rho _0\delta =0.$$
(5)
For a spatially flat universe dominated by pressureless matter, $`\rho _0(t)=1/6\pi Gt^2`$ and Eq. (5) admits two linearly independent power law solutions $`\delta (𝒙,t)=b_\pm (t)\delta _0(𝒙)`$, where $`\delta _0(𝒙)`$ is the initial condition, $`b_+(t)a(t)t^{2/3}`$ is the growing mode and $`b_{}(t)t^1`$ is the decaying mode.

It is useful to expand the density contrast in Fourier series, in which $`\delta `$ is treated as a superposition of plane waves:

$$\delta (𝒙)=\stackrel{~}{\delta }(𝒌)\mathrm{exp}(i𝒌𝒙).$$
(6)
The Fourier transform $`\stackrel{~}{\delta }(𝒌)`$ is complex and therefore possesses both amplitude $`|\stackrel{~}{\delta }(𝒌)|`$ and phase $`\varphi _𝒌`$ where

$$\stackrel{~}{\delta }(𝒌)=|\stackrel{~}{\delta }(𝒌)|\mathrm{exp}(i\varphi _𝒌).$$
(7)
In the standard picture of ‘gravitational instability’ model for the origin of cosmic structure, particularly those involving inflation, the initial perturbations are Gaussian \[Bardeen et al. 1986\]. The most relevant property of Gaussian random fields is that they possess Fourier modes whose real and imaginary parts are independently distributed. In other words, they have phase angles $`\varphi _𝒌`$ that are independently distributed and uniformly random on the interval $`[0,2\pi ]`$. Terms in the perturbative evolution equations for the Fourier modes that represent coupling between different waves are of second (or higher) order in $`\delta `$ and these are neglected in linear perturbation theory. When fluctuations are small, i.e., during the linear regime, the Fourier modes evolve independently (Eq.( 5) and (6)) and the Gaussian character is retained in the linear regime, where the phases remain independent and uniformly random. In the later stages of evolution, however, modes begin to couple together. In this non-linear regime that Fourier phases become non-random. For a thorough review of the theory and implications of non-linear evolution from the point of view of perturbation theory, see Bernardeau et al. (2002).

Standard methods of analysis proceed via the power-spectrum, $`P(k)`$, essentially proportional to $`|\stackrel{~}{\delta }(𝒌)|^2`$. The probabilistic properties of Gaussian random fields are completely specified by knowledge of $`P(k)`$. Higher-order quantities based on $`\stackrel{~}{\delta }(𝒌)`$ can also be defined, such as the bispectrum \[Peebles 1980, Matarrese, Verde & Heavens 1997, Scoccimarro et al. 1998, Scoccimarro, Couchman & Frieman 1999, Verde et al. 2000, Verde et al. 2001, Verde et al. 2002\], which vanishes for Gaussian fields, or quantities related to correlations of $`|\stackrel{~}{\delta }(𝒌)|^2`$ \[Stirling & Peacock 1996\]. Phase coupling results in a non-Gaussian field in which the bispectrum and higher-order polyspectra may be non-zero \[Watts & Coles 2002\]. Phase information is at the heart of non-linear galaxy clustering.

## 3 Directional Phase Mapping

There are two principal difficulties involved in constructing a statistic from Fourier phases. One is that because phases reflect the morphology, their values change according to the position of structural features \[Chiang 2001\]. For example, the phases of a peak in the form of Dirac-$`\delta `$ function $`\delta _D(xx_0)`$, $`\varphi _k=kx_0`$ suffer change in slope along the $`k`$-axis when there is shift of the peak $`\delta _D(xx_0x^{^{}})`$, the phases being $`\varphi _k=k(x_0+x^{^{}})`$. If a pattern is statistically homogeneous, any descriptor of it should be translation-invariant and this is manifestly not true of the phases themselves. The other problem is that phases are of circular measures and therefore defined modulo $`2\pi `$. Traditional measures of association, such as covariances of the form $`\varphi _i\varphi _j`$, are based on the assumption that the measure associated with the variable is linear and are therefore not appropriate to cases where values separated by $`2\pi `$ are in fact equal in value.

To address these problems, Chiang & Coles (2000) used the phase difference (or phase gradient), $`D_k`$, defined in one dimension by

$$D_k\varphi _{k+1}\varphi _k,$$
(8)
i.e. for neighbouring phases. In two or three dimensions differences can be taken in orthogonal directions. The quantity $`D_k`$ has the twin advantages that for random phases it is also random but upon translation $`x^{}`$ it changes by a constant $`x^{}`$ for all $`k`$. The statistical properties of the set of differences $`D_k`$ contain information about the correlations of neighbouring Fourier modes. Strong correlation of the neighbouring modes at large $`k`$ corresponds to the highest peak in the clustering pattern. Naselsky, Novikov and Silk (2002) used this characteristic to extract point sources in the CMB map. Moreover, Chiang et al. (2001) used the phase analysis for extracting the in-flight beam shape properties of CMB experiments. To construct a more general description of phase coupling we need to extend this method to modes that are not necessarily neighbours. We do this by constructing a directional phase map, based on the return maps used in non-linear dynamics \[May 1976\].

The basic idea is simple and based on a study of one-dimensional examples contained in Chiang & Coles (2000) which provides a useful illustration of the more general approach. With a set of phases $`\varphi _k`$ from the Fourier transform of a one-dimensional process, one can plot a map of $`\varphi _k`$ against $`\varphi _{k+1}`$ for each pair $`(\varphi _k,\varphi _{k+1})`$. If the phases are random this will be a scatter plot with points distributed randomly within the bounded square of side $`[0,2\pi ]`$. If there is association between neighbouring phases the plot will contain correlations; the quantity $`D_k`$ is sensitive to linear association. If the spatial pattern consists of a single high-amplitude peak the points display linear association and are mapped into a diagonal lines on the diagram.

In what follows, for illustration, we shall use two-dimensional examples based on numerical simulations with periodic boundaries, so we take $`\mathrm{\Delta }𝒌(m,n)`$ where $`m`$ and $`n`$ are integers. The simulations are done on a $`512^2`$ grid. In Fig. 1 we show 4 sets of such $`N`$–body simulations for initial power spectral index $`n=1`$, 0, 1 and 2 (see Chiang & Coles 2000 for details of the simulations). The evolutionary stages are characterized by an increasing scale of non-linearity defined by $`\left(\delta \rho /\rho \right)^2_{k_{\mathrm{NL}}}=b_+^2(t)_0^{k_{\mathrm{NL}}}P(k)d^2k=1`$, where $`b_+(t)`$ is the growing mode of the linear density contrast and $`P(k)`$ is the linear extrapolation of the initial power spectrum. This definition of $`k_{\mathrm{NL}}`$ identifies the corresponding scale $`2\pi k_{\mathrm{NL}}^1`$ as the boundary between linearity and non-linearity. The stages in Fig. 1 are chosen such that the scales of non-linearity $`k_{\mathrm{NL}}`$ between any two successive stages vary by a factor of 2. The levels of non-linearity of the stages are thus $`k_{\mathrm{NL}}=256`$, 128, 64, 32, 16, 8 and 4.111To avoid confusion with the panels in the captions, we re-name the stages as 1-7, which are originally named as stage $`a`$-$`g`$ in Chiang & Coles (2000).

We can use these simulations to illustrate how we extend the return mapping between neighbouring phases ($`\mathrm{\Delta }𝒌=1`$) to pairs of phases with any $`\mathrm{\Delta }𝒌`$ scales in $`k`$-space. We map all pairs of phases $`\varphi (i,j)`$ and $`\varphi (i+m,j+n)`$ onto the $`x`$ and $`y`$ values of the return map. The axes therefore range over $`[0,2\pi ]`$ for both $`\varphi (𝒌)`$ ($`x`$) and $`\varphi (𝒌+\mathrm{\Delta }𝒌)`$ ($`y`$) axes. For example, for $`(m,n)=(4,6)`$ we have points on the return map $`(\varphi (i,j),\varphi (i+4,j+6))`$ for all $`i[255,256]`$, $`j[1,256]`$, i.e., all points $`(\varphi (1,1),\varphi (5,7))`$, $`(\varphi (2,1),\varphi (6,7))`$, $`(\varphi (1,2),\varphi (5,8))\mathrm{}`$, from a 2D Fourier transform of a realisation. This represents the directional phase coupling for coupling scale $`\mathrm{\Delta }𝒌(\mathrm{\Delta }k_x,\mathrm{\Delta }k_y)(m,n)`$ in $`𝒌`$-space. The neighbouring phase differences in the $`k_x`$-direction and $`k_y`$-direction used by Chiang & Coles (2000) simply corresponds to $`(m,n)=(1,0)`$ and $`(0,1)`$ respectively.

In Fig. 2 (a) we show one example of phase mapping from the realization of stage 3, spectral index $`n=1`$ simulation. The particular $`\mathrm{\Delta }𝒌`$ in this example is $`(m,n)=(0,33)`$. This panel demonstrates how weak coupling between phase pairs with fixed scale $`\mathrm{\Delta }𝒌`$ can manifest itself in the return map as non-uniform density in the map plane.

This directional phase mapping approach circumvents the problem of the circular character of $`\varphi _𝒌`$ but does not attempt to condense all the related information into a single quantity. It, on the other hand, exploits all the information between all Fourier modes. For example, the phase coupling of a 1D distribution can be expressed in a (2D) return map. It therefore allows us to build simple statistics to test the significance of general non-randomness. The phase difference between any pair at a fixed scale becomes a single point on the return map for that scale. The circular characteristic of $`\varphi _𝒌`$ is transferred to a bounded square, topologically equivalent to a torus owing to the periodicity of $`x`$ and $`y`$ axes. The bands seen in Fig. 2 (a) therefore correspond to twisted linear features on this torus.

The key advantage of directional phase mapping is that, for a Gaussian random field, any directional phase mapping for any scale $`\mathrm{\Delta }𝒌(m,n)`$ should produce a random Poisson distribution. Weak phase coupling will produce correlations at large vectors $`(m,n)`$ while strong non-linearity will produce highly non-uniform phase maps at all scales $`(m,n)`$.

## 4 A $`\chi ^2`$ Test on Phase Maps

Once we have transferred the phase information onto a phase map like that shown in Fig. 2 (a), many different statistical tests can be used to analyse its properties. Here we outline a simple yet powerful method.

First we smooth the return map. Smoothing enhances our visualization of the pattern of phase coupling. In Fig. 2 (b), (c), and (d) we show the contour maps of the smoothed return map. In these simple illustrative experiments we divide the square of the return map into $`128^2`$ pixels, and we bin the 32 768 points so that, for a perfectly even distribution in the map plane the occupation of each pixel is $`p(i,j)=2`$. Then we smooth this $`128^2`$ mesh by

$$p(𝒙,R)=(\sqrt{2\pi }R)^2d^2𝒙^{^{}}p(𝒙^{^{}})\mathrm{exp}\left(\frac{|𝒙𝒙^{^{}}|^2}{2R^2}\right).$$
(9)
The contour levels are drawn upwards starting from the mean value. Panel (b) is the contour of a smoothed return map from a realization of random phases. Panel (c) is that from stage 2, $`n=1`$ of Fig. 1 with $`(m,n)=(2,44)`$. For comparison, the coupling scale $`\mathrm{\Delta }𝒌`$ for both (b) and (c) is set the same. We can see that even for the mildly non-linear regime represented by stage 2 (for which $`\delta 1.8`$), phase mapping does reveal the existence of coupling by starting to condense on the diagonal strip. Panel (d) is the smoothed version of (a) with $`(m,n)=(0,33)`$, in which the condensed strip from (a) is much clearer.

We define a mean $`\chi ^2`$ statistic as

$$\overline{\chi ^2}=\frac{1}{M}\underset{i,j}{}\frac{\left[p(i,j)\overline{p}\right]^2}{\overline{p}}$$
(10)
where $`M`$ is the number of pixels we assign on the return map and $`\overline{p}`$ is the mean value for each pixel.

In Fig. 2, the smoothing scale on the $`128^2`$ mesh is $`R=2`$ and the contour levels are drawn starting from the mean value. We use this pixel size and smoothing scale in all our calculations of $`\overline{\chi ^2}`$ in the following, but one can vary the scale as part of a statistical test.

## 5 Simulation Results

We now illustrate the results of this analysis using the 2D $`N`$-body simulations described earlier. First we Fourier transform the realizations of the $`N^2`$ mesh ($`N=512`$ in our simulations). Because of the reality of the original distribution, and the consequent Hermitian conjugate relations in the Fourier image, only half of the Fourier transform contains independent information. We end up with $`N^2/2`$ Fourier modes available, so we take $`k_{N/2}+1`$ to $`k_{N/2}`$ in $`k_x`$ axis and from 1 to $`k_{N/2}`$ in $`k_y`$ axis.

We carry out a phase mapping for each $`(m,n)`$. The directional phase mapping is performed for the vectors $`(m,n)`$, where $`m[49,50]`$ and $`n[1,50]`$ for the phase base $`\varphi (i,j)`$, where $`i[127,128]`$, and $`j[1,128]`$, a quarter of the available phases. The limited range of mapping vectors is chosen to ensure that the map can be constructed without running out of sensible wavenumbers. Thus, in such a case, for each return map of $`(m,n)`$ there are $`2\times 128^2`$ points. We set up $`M=128^2`$ pixels in this $`[0,2\pi ]`$ square and then smooth the return map to decrease the fluctuation. The $`\overline{\chi ^2}(m,n)`$ for each fixed scale $`(m,n)`$ is calculated by Eq. (10).

Figure 3 shows ‘supermaps’ of all the phase information contained in the $`\overline{\chi ^2}(m,n)`$ using a grey scale in the $`(m,n)`$ plane from various realisations of simulations. We are interested specifically on the mild non-linear regime, so in Fig. 3 the panels (a), (b), (c) and (d) correspond to stages 2, 3, 4 and 5 of $`n=1`$ in Fig. 1, respectively. The $`m`$-axis ranges from $`49`$ to $`50`$ and $`n`$-axis ranges from 1 to 50. Each pixel therefore displays the level of phase coupling on a certain scale $`(\mathrm{\Delta }k_x,\mathrm{\Delta }k_y)(m,n)`$ in terms of the $`\overline{\chi ^2}`$ statistics. The bright points in Figure 3 are direct indications of phase coupling for the corresponding $`\mathrm{\Delta }𝒌`$ scales.

The reason we present the specific scales $`(m,n)=(2,44)`$ for stage 2 and $`(m,n)=(0,33)`$ for stage 3, $`n=1`$ in Fig. 2 is that the calculation of the $`\overline{\chi ^2}`$ statistics shows them to be ‘hotspots’ of phase correlation. The scales of phase coupling for those two panels correspond to the brightest points, i.e. the maximum $`\overline{\chi ^2}`$ on the corresponding $`(m,n)`$ planes in Fig. 3.

In Fig. 4, 5 and 6 we plot the maximum $`\overline{\chi ^2}`$ against $`|\mathrm{\Delta }𝒌|`$ for all 4 sets of simulations. These 1D plots show that the maximal phase coupling from each ring $`|\mathrm{\Delta }𝒌|\sqrt{m^2+n^2}`$ of the $`(m,n)`$ plane. As we control the simulations by assigning the same set of initial random phases, These 1D plots will display how the scale $`|\mathrm{\Delta }𝒌|`$ of phase coupling is related to morphology.

The $`N`$-body simulations we have carried out is of self-similar nature, that is, a distribution function $`f(𝒙,t)`$ has the same statistical measure as the re-scaled one

$$f\lambda ^\alpha f(𝒙/\lambda ^\beta ,\lambda t),\text{ as }t\lambda t.$$
(11)
With the reciprocal-scaling property of Fourier transform,

$$f(\alpha 𝒙)=\frac{1}{|\alpha |}F\left(\frac{𝒌}{\alpha }\right),$$
(12)
for non-linear scales $`𝒙_{\mathrm{NL}}`$ increasing, i.e. $`\alpha >1`$, the corresponding scales in $`k`$ space are decreasing. Although phases do not possess a linear relationship owing to their circular nature, it can be understood qualitatively that, if there exists coupling between pairs of phases with fixed $`\mathrm{\Delta }𝒌`$, this scale has to decrease as gravitational clustering proceeds.

It is therefore clear that for the case such as $`n=1`$, where large-scale filaments are the prominent feature, phase coupling starts from low $`\mathrm{\Delta }𝒌`$, which also has cascade effect on to higher $`\mathrm{\Delta }𝒌`$, as phases strongly couple on any $`\mathrm{\Delta }𝒌`$ might also do so at multiples of $`\mathrm{\Delta }𝒌`$. For high $`n`$, on the other hand, where small clumps form first, phase coupling starts from large $`\mathrm{\Delta }𝒌`$. This also explains why the Shannon entropy from neighbouring phase difference can produce considerable results at early stages for $`n=1`$, but not for $`n=2`$ (see the fig.5 of Chiang & Coles 2000). On the other hand, the coupling between the amplitudes are enhanced by the factor $`1/\alpha `$ as clustering proceeds, so mode-mode coupling between amplitudes at early stages is not as obvious as between phases. For a discussion of how this relates to the development of phase correlations in mildly non-linear evolution, see Watts & Coles (2002).

The $`\overline{\chi ^2}`$ statistics shown on the $`(m,n)`$ plane and those 1D plots confirm the visualization of phase coupling presented by Coles & Chiang (2000). Phase coupling firstly appears on large $`\mathrm{\Delta }𝒌`$ when the the scale of non-linearity is small in real space, then it shifts on the $`(m,n)`$ plane to small $`\mathrm{\Delta }𝒌`$ in Fig. 3 (c) and (d) and finally dominates at neighbouring modes as seen in Fig. 6.

## 6 Conclusion

We have generalized a method based on phase mapping on the return map. This simple, easy-to-implement method can detect phase coupling at any scales $`\mathrm{\Delta }𝒌`$ in $`k`$ space. We apply this method to two-dimensional simulations of gravitational clustering and the result has shown that even when the evolution is in the mild non-linear regime, phase coupling on certain scale is revealed through the $`\overline{\chi ^2}`$ statistics on the $`(m,n)`$ plane.

In contrast to other methods, such as the Shannon entropy of the distribution of neighbouring phase differences \[Chiang & Coles 2000\], this method does not require large number of phases. Moreover, this approach can detect the scale of phase coupling through the phase mapping as shown in Fig. 2.

With the systematic $`N`$-body simulations shown in Fig. 1, we have also demonstrated in Fig. 4-6 that the scale of phase coupling differs according to the clustering morphology: modes between small $`\mathrm{\Delta }𝒌`$ for large-scale filaments, large for small clumps.

This method reveals a signature of non-linear gravitational instability, but also offers the opportunity to provide a general test of Gaussianity that could be applied to cosmic microwave background temperature maps. In future work we shall evaluate the effectiveness for such method.

## Acknowledgments

This paper was supported in part by Danmarks Grundforskningsfond through its support for the establishment of the Theoretical Astrophysics Center and by grants RFBR 17625. PC acknowledges support from PPARC.


text:
# References

HOPF ALGEBRA, THERMODYNAMICS AND ENTANGLEMENT
IN QUANTUM FIELD THEORY

A. IORIO1, G. LAMBIASE1 AND G. VITIELLO1,2

1Dipartimento di Fisica “E.R.Caianiello”, Università di Salerno, 84100 Salerno, Italy

and INFN, Gruppo Collegato di Salerno

2INFM, Sezione di Salerno

iorio@sa.infn.it

lambiase@sa.infn.it

vitiello@sa.infn.it

Abstract The quantum deformation of the Hopf algebra describes the skeleton of quantum field theory, namely its characterizing feature consisting in the existence of infinitely many unitarily inequivalent representations of the canonical commutation relations. From this we derive the thermal properties of quantum field theory, with the entropy playing the role of the generator of the non-unitary time evolution. The entanglement of the quantum vacuum appears to be robust against interaction with the environment: on cannot “unknot the knots” in the infinite volume limit.

1. Introduction

It is now long known that quantum deformations  of the enveloping algebras of Lie algebras have a Hopf algebra structure . It has been shown  that dissipative systems, as well as finite temperature non-equilibrium systems, are properly described in the frame of the q-deformed Hopf algebra. The analysis of such systems has revealed that the proper algebraic structure of quantum field theory (QFT) is the deformed Hopf algebra. The q-parameter acts as a label for the infinitely many unitarily inequivalent representations (uir) of the canonical commutation relations (ccr). Our task in this paper is to report about such results on the architecture of QFT and, by resorting to that, to show its intrinsic thermodynamic nature. We also want to comment on the robustness of the entanglement of the quantum vacuum against interaction with the environment .

Let us illustrate in a simple way that the Hopf algebra enters in the QFT formalism since the beginning. We observe that the introduction of the operator algebra necessary to set up QFT is usually limited to the introduction of the boson Weyl-Heisenberg (WH) algebra. The additivity of some observables such as the energy, the momentum and the angular momentum is such an obvious assumption that one does not even bother to spell it out. It is implicitly given as granted. However, if one is asked to express it explicitly and formally, then one realizes that the boson WH algebra is only a part of the full algebraic structure and it becomes natural to introduce the “coproduct” map, namely, for the “addition” of, e.g., the angular momentum $`J^\alpha `$, $`\alpha =1,2,3`$, of two particles, $`\mathrm{\Delta }J^\alpha =J^\alpha \mathrm{𝟏}+\mathrm{𝟏}J^\alpha J_1^\alpha +J_2^\alpha `$. It is, however, not only in this simple way that one needs the Hopf structure where the coproduct mapping plays a crucial role. The full algebraic structure which is needed has to take also into account  one of the very special features of QFT, namely the existence of infinitely many representations of the ccr  (this is a characterizing feature which makes QFT deeply different from quantum mechanics (QM). In QM all the representations of the ccr are unitary equivalent due to the von Neumann theorem). Then one is led to consider the q-deformation of the Hopf algebra. Let us see in the following how this goes.

2. Deformed Hopf algebra and the quantum field theory structure

In the following we shall focus on the case of bosons for simplicity. However, our conclusions can be extended also to fermions .

The coproduct is a homomorphism which duplicates the algebra, $`\mathrm{\Delta }:𝒜𝒜𝒜`$. The operator doubling implied by the coalgebra is a key ingredient of Hopf algebras. As already said, Lie-Hopf algebras are commonly used in the familiar addition of energy, momentum and angular momentum. Thus, the operational meaning of the coproduct is that it provides the prescription for operating on two modes. Associated to that, there is the doubling of the degrees of freedom of the system which is rich of physical meanings (e.g. near a black hole such a doubling perfectly describes the modes on the two sides of the horizon ).

The bosonic Hopf algebra for a single mode (the case of modes labelled by the momentum $`\kappa `$ is straightforward), also called $`h(1)`$, is generated by the set of operators $`\{a,a^{},H,N\}`$ with commutation relations:

$$[a,a^{}]=2H,[N,a]=a,[N,a^{}]=a^{},[H,]=0,$$
(1)
where $`a`$ ($`a^{}`$) is the generic annihilation (creation) operator and $`H`$ is a central operator, constant in each representation. The Casimir operator is given by $`𝒞=2NHa^{}a`$. In $`h(1)`$ the coproduct is defined by

$$\mathrm{\Delta }𝒪=𝒪\mathrm{𝟏}+\mathrm{𝟏}𝒪𝒪_1+𝒪_2,$$
(2)
where $`𝒪`$ stands for $`a`$, $`a^{}`$, $`H`$ and $`N`$. The $`q`$-deformation of $`h(1)`$ is the Hopf algebra $`h_q(1)`$:

$$[a_q,a_q^{}]=[2H]_q,[N,a_q]=a_q,[N,a_q^{}]=a_q^{},[H,]=0,$$
(3)
where $`N_qN`$ and $`H_qH`$ and $`[x]_q={\displaystyle \frac{q^xq^x}{qq^1}}`$. The Casimir operator is given by $`𝒞_q=N[2H]_qa_q^{}a_q`$. The coproduct for $`a_q`$ and $`a_q^{}`$ now changes to

$$\mathrm{\Delta }a_q=a_qq^H+q^Ha_q,\mathrm{\Delta }a_q^{}=a_q^{}q^H+q^Ha_q^{},$$
(4)
while, of course, stays the same for $`H`$ and $`N`$.

In the fundamental representation, obtained by setting $`H=1/2`$, $`𝒞=0`$, $`h(1)`$ and $`h_q(1)`$ coincide, as it happens for the spin-$`\frac{1}{2}`$ representation of $`su(2)`$ and $`su_q(2)`$. The differences appear in the coproduct and in the higher spin representations. We shall denote by $`_1`$ this representation space (single mode Fock space).

As customary, one requires $`a`$ and $`a^{}`$ ($`a_q`$ and $`a_{q}^{}{}_{}{}^{}`$) to be adjoint operators. This implies that $`q`$ can only be real or of modulus one. In the two mode Fock space $`_2_1_1`$, for $`|q|=1`$, the hermitian conjugation of the coproduct must be supplemented by the inversion of the two spaces for consistency with the coproduct homomorphism. Summarizing, on $`_2=_1_1`$ it can be written:

$$\mathrm{\Delta }a_q=a_1q^{1/2}+q^{1/2}a_2,\mathrm{\Delta }a_q^{}=a_1^{}q^{1/2}+q^{1/2}a_2^{},$$
(5)
$$\mathrm{\Delta }H=1,\mathrm{\Delta }N=N_1+N_2.$$
(6)
Note that $`[a_i,a_j]=[a_i,a_j^{}]=0,ij`$. It is now possible to show that the full set of infinitely many uir of the ccr in QFT are classified by use of the deformed Hopf algebra. To do that it is sufficient to show that the Bogolubov transformations are directly obtained by use of the deformed copodruct operation. As well known, indeed, the Bogolubov transformations relate different (i.e. unitary inequivalent) representations. We consider therefore the following operators (cf. (4) with $`H=1/2`$ and $`q(\theta )e^{2\theta }`$):

$$\alpha _{q(\theta )}\frac{\mathrm{\Delta }a_q}{\sqrt{[2]_q}}=\frac{1}{\sqrt{[2]_q}}(e^\theta a_1+e^\theta a_2),\beta _{q(\theta )}\frac{1}{\sqrt{[2]_q}}\frac{\delta }{\delta \theta }\mathrm{\Delta }a_q=\frac{1}{\sqrt{[2]_q}}(e^\theta a_1e^\theta a_2),$$
(7)
and h.c.. A set of commuting operators with canonical commutation relations is given by

$$\alpha (\theta )\frac{\sqrt{[2]_q}}{2\sqrt{2}}[\alpha _{q(\theta )}+\alpha _{q(\theta )}\beta _{q(\theta )}^{}+\beta _{q(\theta )}^{}],$$
(8)
$$\beta (\theta )\frac{\sqrt{[2]_q}}{2\sqrt{2}}[\beta _{q(\theta )}+\beta _{q(\theta )}\alpha _{q(\theta )}^{}+\alpha _{q(\theta )}^{}],$$
(9)
and h.c. One then introduces

$`A(\theta ){\displaystyle \frac{1}{\sqrt{2}}}\left(\alpha (\theta )+\beta (\theta )\right)`$ $`=`$ $`A\mathrm{cosh}\theta B^{}\mathrm{sinh}\theta ,`$ (10)
$`B(\theta ){\displaystyle \frac{1}{\sqrt{2}}}\left(\alpha (\theta )\beta (\theta )\right)`$ $`=`$ $`B\mathrm{cosh}\theta A^{}\mathrm{sinh}\theta ,`$ (11)
with $`[A(\theta ),A^{}(\theta )]=1,[B(\theta ),B^{}(\theta )]=1`$. All other commutators are equal to zero and $`A(\theta )`$ and $`B(\theta )`$ commute among themselves. For homogeneity of notation, we have used $`Aa_1`$ and $`Ba_2`$.

Eqs. (10) and (11) are nothing but the Bogolubov transformations for the $`(A,B)`$ pair. In other words, the Bogolubov-transformed operators $`A(\theta )`$ and $`B(\theta )`$ are linear combinations of the coproduct operators defined in terms of the deformation parameter $`q(\theta )`$ and of their $`\theta `$-derivatives. Notice that

$$i\frac{\delta }{\delta \theta }A(\theta )=[𝒢,A(\theta )],i\frac{\delta }{\delta \theta }B(\theta )=[𝒢,B(\theta )],$$
(12)
and h.c., where $`𝒢i(A^{}B^{}AB)`$ denotes the generator of (10) and (11). For a fixed value $`\overline{\theta }`$, we have

$$\mathrm{exp}(i\overline{\theta }p_\theta )A(\theta )=\mathrm{exp}(i\overline{\theta }𝒢)A(\theta )\mathrm{exp}(i\overline{\theta }𝒢)=A(\theta +\overline{\theta }),$$
(13)
and similar equations for $`B(\theta )`$. In eq.(13) the definition $`p_\theta i{\displaystyle \frac{\delta }{\delta \theta }}`$ has been used. It can be regarded as the momentum operator “conjugate” to the “degree of freedom” $`\theta `$, which thus acquires formal definiteness in the sense of the canonical formalism.

If $`|0`$ denotes the vacuum annihilated by $`A`$ and $`B`$, $`A|0=0=B|0`$, then, at finite volume $`V`$,

$$|0(\theta )=\mathrm{exp}(i\theta 𝒢)|0.$$
(14)
The vacuum $`|0(\theta )`$ turns out to be an $`SU(1,1)`$ generalized coherent state: coherence and the vacuum structure in QFT are thus intrinsically related to the deformed Hopf algebra.

In the continuum limit in the space of momenta, i.e. in the infinite volume limit, the number of degrees of freedom becomes uncountable infinite, hence we obtain $`0(\theta )|0(\theta ^{})0`$ as $`V\mathrm{},\theta ,\theta ^{},\theta \theta ^{}`$, thus the Hilbert spaces $`_\theta `$ and $`_\theta ^{}`$ become unitarily inequivalent. In this limit, the deformation parameter $`\theta ={\displaystyle \frac{1}{2}}\mathrm{ln}q`$ labels the set $`\{_\theta ,\theta \}`$ of the infinitely many uir of the ccr . In conclusion, the deformed Hopf algebra provides the typical structure one deals with in QFT.

Furthermore, in the case in which the deformation parameter is time-dependent, it turns out to be related to the so-called heat-term in dissipative systems. This can be seen by noticing that the Heisenberg equation for $`A(t,\theta (t))`$ is

$$i\dot{A}(t,\theta (t))=i\frac{\delta }{\delta t}A(t,\theta (t))i\frac{\delta \theta }{\delta t}\frac{\delta }{\delta \theta }A(t,\theta (t))=$$

$$[H,A(t,\theta (t))]+\frac{\delta \theta }{\delta t}[𝒢,A(t,\theta (t))]=[H+Q,A(t,\theta (t))],$$
(15)
where $`Q{\displaystyle \frac{\delta \theta }{\delta t}}𝒢`$ is the announced heat-term, and $`H`$ is the Hamiltonian responsible for the time variation in the explicit time dependence of $`A(t,\theta (t))`$. $`H+Q`$ is therefore to be identified with the free energy . Thus, the conclusion is that variations in time of the deformation parameter actually involve dissipation.

When the proper field description is taken into account, $`A`$ and $`B`$ depend on the momentum $`\kappa `$ and, as customary in QFT, one deals with the algebras $`{\displaystyle \underset{\kappa }{}}h_\kappa (1)`$.

3. Entropy, entanglement and environment

It is remarkable that the “conjugate momentum” $`p_\theta `$ generates transitions among the uir (in the infinite volume limit): $`\mathrm{exp}(i\overline{\theta }p_\theta )|0(\theta )>=|0(\theta +\overline{\theta })>`$. Use of eq. (14) shows that

$$\frac{}{\theta _k}|0(\theta )>=\left(\frac{1}{2}\frac{𝒮_𝒜}{\theta _k}\right)|0(\theta )>.$$
(16)
where, in full generality, we are using $`\theta =\theta _k`$ and

$$S_A\underset{\kappa }{}\left\{A_\kappa ^{}A_\kappa \mathrm{ln}\mathrm{sinh}^2\left(\theta _k\right)A_\kappa A_\kappa ^{}\mathrm{ln}\mathrm{cosh}^2\left(\theta _k\right)\right\}.$$
(17)
Thus $`i\left(\frac{1}{2}\mathrm{}\frac{𝒮_𝒜}{\theta }\right)`$ is the generator of translations in $`\theta `$ . The operator $`S_A`$ ($`S_B`$ for the $`B`$-modes) is recognized to be the entropy operator . In dissipative or unstable systems (or in thermal theories at non-equilibrium) the deformation parameter depends on time, $`\theta =\theta (t)`$, and the non-unitary time evolution is controlled by the entropy variations . We thus have found that entropy controls non-unitary time evolution: dissipation implies indeed the choice of a privileged direction in time evolution (the arrow of time) with the consequent breaking of time-reversal invariance. It can be also shown that $`S_AS_B`$ is a constant in time.

By introducing the free energy functional for the $`A`$-modes

$$_A<0(\theta )|\left(H_A\frac{1}{\beta }S_A\right)|0(\theta )>,$$
(18)
where $`H_A=_\kappa \mathrm{}\mathrm{\Omega }_\kappa A_\kappa ^{}A_\kappa `$, one shows that $`d_A=dE_A\frac{1}{\beta }d𝒮_A=0`$. This is the first principle of thermodynamics for a system coupled with environment at constant temperature and in absence of mechanical work. The change in time of particles condensed in the vacuum, $`d𝒩_A`$, turns into heat dissipation, defined as usual $`dQ=\frac{1}{\beta }dS`$, which recovers the conclusion of the previous section (cf. eq. (15)). These observations thus lead us to recognize the thermal features of $`\{|0(\theta (t))>\}`$.

It is also interesting to observe that the thermodynamical arrow of time, whose direction is defined by the increasing entropy direction, points in the same direction of the cosmological arrow of time, namely the inflating time direction for the expanding Universe. This can be shown by considering indeed the quantization of inflationary models . The concordance between the two arrows of time is not at all granted and is a subject of an ongoing debate (see, e.g., ).

Let us now turn our attention to the entanglement. It is convenient to consider the case in which we have two modes, say $`A`$ and $`\overline{A}`$ (as, for example, in the case of a complex scalar field). The corresponding doubled modes are $`B`$ and $`\overline{B}`$. For typographical simplicity, we denote them as $`AA^{(+)}`$, $`\overline{A}\overline{A}^{(+)}`$, $`BA^{()}`$ and $`\overline{B}\overline{A}^{()}`$. We also use $`\sigma =\pm `$. Then the vacuum can be expressed as a $`SU(1,1)\times SU(1,1)`$ generalized coherent state  of Cooper-like pairs

$$|0(\theta )=G(\theta )|0=\frac{1}{Z}\mathrm{exp}\left[\underset{\sigma }{}\underset{k}{}\mathrm{tanh}\theta _kA_k^{(\sigma )}\overline{A}_k^{(\sigma )}\right]|0,$$
(19)
where $`Z=_k\mathrm{cosh}^2\theta _k`$. The state $`|0(\theta )`$ in (19), can be rewritten as

$$|0(\theta )=\frac{1}{Z}\left[|0+\underset{k}{}\mathrm{tanh}\theta _k\left(|1_k^{(+)},\overline{0}|0,\overline{1}_k^{()}+|0,\overline{1}_k^{(+)}|1_k^{()},\overline{0}\right)+\mathrm{}\right],$$
(20)
where, we denote by $`|n_k^{(\sigma )},\overline{m}_k^{(\sigma )}`$ a state of $`n`$ particles and $`m`$ “antiparticles” in whichever sector $`(\sigma )`$. Note that for the generic $`n^{\mathrm{th}}`$ term, the state $`|n_k^{(\sigma )},\overline{0}|1_{k_1}^{(\sigma )},\mathrm{},1_{k_n}^{(\sigma )},\overline{0}`$, and similarly for antiparticles.

By introducing a well known notation, $``$ for a particle, and $``$ for an antiparticle, the two-particle state in (20) can be written as

$$|^{(+)}|^{()}+|^{(+)}|^{()},$$
(21)
which is an entangled state of particle and antiparticle living in the two sectors $`(\pm )`$. The generic $`n^{\mathrm{th}}`$ term in (20) shares exactly the same property as the two-particle state, but this time the $``$ describes a set of $`n`$ particles, and $``$ a set of $`n`$ antiparticles. The mechanism of the entanglement, induced by the q-deformation, takes place at all orders in the expansion, always by grouping particles and antiparticles into two sets. Thus the whole vacuum $`|0(\theta )`$ is an infinite superposition of entangled states111A similar structure also arises in the temperature-dependent vacuum of Thermo-Field Dynamics  (see also ).

$$|0(\theta )=\underset{n=0}{\overset{+\mathrm{}}{}}\sqrt{W_n}|\mathrm{Entangled}_n,$$
(22)
$$W_n=\underset{k}{}\frac{\mathrm{sinh}^{2n_k}\theta _k}{\mathrm{cosh}^{2(n_k+1)}\theta _k},$$
(23)
with $`0<W_n<1\mathrm{and}_{n=0}^+\mathrm{}W_n=1`$. The coefficients $`\sqrt{W_n}`$ of the expansion in Eq. (22) appear also in the discussion about the entropy of the black hole .

Of course, the probability of having entanglement of two sets of $`n`$ particles and $`n`$ antiparticles is $`W_n`$. At finite volume, being $`W_n`$ a decreasing monotonic function of $`n`$, the entanglement is suppressed for large $`n`$. It appears then that only a finite number of entangled terms in the expansion (22) is relevant. Nonetheless this is only true at finite volume (the QM limit), while the interesting case occurs in the infinite volume limit, which one has to perform in a QFT setting.

The entanglement is generated by $`G(\theta )`$, where the field modes in one sector $`(\sigma )`$ are coupled to the modes in the other sector $`(\sigma )`$ via the deformation parameter $`\theta _k`$. Since the deformation parameter may describe the background gravitational field (environment) , the temperature effects , the dissipative effects  and other effects of the environment in which the system is embedded, it appears that the origin of the entanglement is the environment, in contrast with the usual QM view, which attributes to the environment the loss of the entanglement. In the present treatment such an origin for the entanglement makes it quite robust.

One further reason for the robustness is that this entanglement is realized in the limit to the infinite volume once and for all since then there is no unitary evolution to disentangle the vacuum: at infinite volume one cannot ”unknot the knots”. Such a non-unitarity is only realized when all the terms in the series (22) are summed up, which indeed happens in the $`V\mathrm{}`$ limit.

As a final comment let us observe that the doubling implied by the deformed Hopf algebra has been useful in several applications, ranging from unstable particles , to coherence in quantum Brownian motion , squeezed states in quantum optics , topologically massive theories in the infrared region in 2+1 dimensions , the Chern-Simons-like dynamics of Bloch electrons in solids , the quantization of matter in curved background . These features are also common to two-dimensional gravity models  and to the study of quantization arising from the loss of information . Moreover, it has been applied  to the study of the memory capacity problem in the quantum model for the brain.


text:
# Exchange interaction parameters and adiabatic spin-wave spectra of ferromagnets: A “renormalized magnetic force theorem”

(24 July, 2002; revised 6 November 2002)

## Abstract

The “magnetic force theorem” is frequently used to compute exchange interaction parameters and adiabatic spin-wave spectra of ferromagnets. The interest of this approach is that it allows to obtain these results from a non-self-consistent calculation of the (single-electron) band energy only, which greatly reduces the computational effort. However, as discussed by various authors, this approach entail some systematic error. Here, a “renormalized magnetic force theorem” allowing to remove this systematic error without increasing significantly the computational effort is presented. For systems with one atom per unit cell, it amounts to a simple renormalization of the spin-wave spectrum. This renormalization greatly improves the agreement between calculated and experimental Curie temperatures of Fe and Ni.
published in: Phys. Rev. Lett. 90, 087205 (2003)

The ab initio study of magnetic interactions in magnetic metals and of interlayer exchange coupling in multilayers has been the subject of a large number of publications during the past few years. Many these studies rely explicitly or implicitly on the use of the “magnetic force theorem” (MFT), which allows to perform the calculations non-self-consistently, and by taking into account only single-particle energies. This is of great practical interest for it reduces the computational effort by several orders of magnitude.

The force theorem was first introduced for the case of non-magnetic systems Mackintosh1980 . Extensions to the case of magnetic systems have been published by Oswald et al. Oswald1985 , who focused on the case of magnetic impurities embedded in a non-magnetic host, and by Liechtenstein et al. Liechtenstein1984 , who addressed the case of exchange interactions and spin-wave spectra of ferromagnetic systems. However, as mentioned by Liechtenstein et al. in their 1984 paper, their MFT yields exact results only in the limit of infinite magnon wavelength. At finite magnon wavelength, their prescription entails some systematic error, as emphasized recently by various authors Bylander2000 ; Grotheer2001 .

In the present Letter I present a “renormalized MFT”, which corrects the systematic error entailed by the MFT of Liechtenstein et al. Liechtenstein1984  (hereafter called the “bare MFT”). In the case of periodic systems with a single site per unit cell, this results in an extremely simple renormalization of the spin-wave spectrum and Curie temperature. The Curie temperatures calculated with this procedure are found to be in excellent agreement with the experimental ones, both for Fe and Ni.

The problem we want to address here is that of calculating the energy required to create some static transverse deviation of the magnetic moments in a ferromagnet. The physical motivation for this is that, within the adiabatic approximation (i.e., provided that the characteristic time scale for the dynamics of such fluctuations is long as compared with the time scale for the dynamics of band electrons), it provides a good approximation to the energy of the low-lying magnetic excitations and therefore allows to address from an ab initio point of view the thermodynamics of itinerant ferromagnets. In practice this results into a mapping of the complicated itinerant ferromagnetic system onto an effective Heisenberg system having the same low-lying excitations, but much less degrees of freedom. Let $`n(𝐫)`$ and $`𝐦(𝐫)m(𝐫)𝐮(𝐫)`$ be, respectively, the charge- and spin-density. For transition-metal systems, it is usually justified to neglect (as is almost always done) intra-atomic non-collinearity and fluctuations (which are generally expected to cost a high energy). We therefore have $`𝐮(𝐫)=𝐮_𝐑`$ if $`𝐫`$ belongs to the atomic cell $`\mathrm{\Omega }_𝐑`$ around atom $`𝐑`$, and we wish to impose a prescribed direction to the unit vector $`𝐮_𝐑`$ of each site $`𝐑`$ generalization\_continuum .

The proper way of doing this relies on the “constrained density functional theory” of Dederichs et al. Dederichs1984 , which, in the present case, amounts to introduce some local external field $`𝐁_𝐑^{}`$, perpendicular to the local magnetic moment axis $`𝐮_𝐑`$, playing the role of Lagrange parameters. As we wish to constrain two components of the magnetic moment for each site, we have two Lagrange parameters per site, which are the magnitude and azimuthal angle of the local constraining field.

In (non-relativistic) spin-density functional theory Hohenberg1964 , the ground state energy $`E_0`$ is obtained by minimizing the Hohenberg-Kohn functional (HKF) $`_{\text{HK}}[\rho ]=𝒯_0[\rho ]+_{\text{ext}}[\rho ]+_H[\rho ]+_{\text{xc}}[\rho ],`$ with respect to the spinor-density $`\rho (n\sigma _0+𝐦𝝈)`$, where $`\sigma _0`$ is the unit spinor and $`𝝈`$ the vector-spinor whose components are the Pauli matrices. The various terms are respectively, the kinetic energy of a non-interacting system having the same spinor-density, the potential energy, the Hartree part of the Coulomb energy, and the exchange-correlation energy. The constrained ground state energy is found by introducing the new constrained HKF (CHKF) Dederichs1984  $`_{\text{HK}}[\rho ,𝐁^{}]_{\text{HK}}[\rho ]+_{\text{cons}}[\rho ,𝐁^{}],`$ where the constraint term is given by $`_{\text{cons}}[\rho ,𝐁^{}]_𝐑𝐁_𝐑^{}_{\mathrm{\Omega }_𝐑}d𝐫𝐦(𝐫),`$ and minimizing $`_{\text{HK}}[\rho ,𝐁^{}]`$ with respect to $`\rho (𝐫)`$, and to the set $`𝐁^{}\{𝐁_𝐑^{}\}`$. As we want to calculate the energy change associated with some infinitesimal transverse fluctuation $`\delta 𝐮_𝐑`$ away from the ferromagnetic configuration, we may attempt to use the variational properties of the CHKF in order to minimize the computational effort. The difficulty lies in the calculation of the kinetic energy term, which requires the knowledge of the effective one-electron spinor-potential $`w_{\text{eff}}[\rho ]V_{\text{eff}}\sigma _0𝐁_{\text{eff}}𝝈`$ of energy eigenvalues $`\epsilon _i(w_{\text{eff}}[\rho ])`$ (labeled in order of increasing energy) and eigenvectors $`\psi _i(𝐫)`$ which yields $`\rho `$ as output spinor-density, i.e, $`\rho (𝐫)=_{i=1}^{N_{\text{el}}}\left\{\left|\psi _i(𝐫)\right|^2\sigma _0+\left[\psi _i^{}(𝐫)𝝈\psi _i(𝐫)\right]𝝈\right\}.`$ The kinetic energy is then given by $`𝒯_0[\rho ]=_{i=1}^{N_{\text{el}}}\epsilon _i(w_{\text{eff}}[\rho ])d𝐫\left(nV_{\text{eff}}𝐦𝐁_{\text{eff}}\right).`$ This implicit dependence of $`w_{\text{eff}}[\rho ]`$ upon the spinor-density $`\rho `$ is the origin of our difficulties. In fact, starting from a trial input potential and magnetic field and solving the one-electron problem, we are only able to obtain the value of the HKF for the output spinor-density, which we did not know a priori, but not for some spinor-density chosen a priori. This problem can be circumvented by using an auxiliary energy functional, first introduced by Harris Harris1985 . The Harris functional (HF) $`_{\text{Harris}}[\rho ,𝐁^{}]`$, in the present context, has the same form as the CHKF $`_{\text{HK}}[\rho ,𝐁^{}]`$, with $`𝒯_0[\rho ]`$ replaced by

$$𝒯^{}[\rho ,𝐁^{}]\underset{i=1}{\overset{N_{\text{el}}}{}}\epsilon _i(w_{\text{eff}}^{}[\rho ,𝐁^{}])d𝐫\left(nV_{\text{eff}}^{}𝐦𝐁_{\text{eff}}^{}\right),$$
(1)
where the new effective spinor-potential $`w_{\text{eff}}^{}[\rho ,𝐁^{}]V_{\text{eff}}^{}\sigma _0𝐁_{\text{eff}}^{}𝝈`$ is defined explicitly in terms of $`\rho `$ and $`𝐁^{}`$ by

$$w_{\text{eff}}^{}[\rho ,𝐁^{}](𝐫)\frac{\delta \left(_{\text{ext}}[\rho ]+_H[\rho ]+_{\text{xc}}[\rho ]+_{\text{cons}}[\rho ,𝐁^{}]\right)}{\delta \rho (𝐫)},$$
(2)
where the functional derivative is taken at the a priori prescribed spinor density. This essential difference between $`w_{\text{eff}}^{}[\rho ,𝐁^{}]`$ and $`w_{\text{eff}}[\rho ]`$ allows us to calculate explicitly the value of the HF for a spinor-density $`\rho `$ and constrains $`𝐁^{}`$ chosen a priori (provided we know some suitable approximation of $`_{\text{xc}}[\rho ]`$), which constitutes a great advantage with respect to the CHKF. In view of the Kohn-Sham theorem Hohenberg1964 , the HF and CHKF obviously take the same value $`E_0`$ for the density $`\rho ^{}`$ and constrains $`𝐁^{}`$ corresponding to the constrained ground state (solution of the constrained Kohn-Sham equation). In addition, the HF can be shown to be stationary (but not necessarily minimal, in contrast to the CHKF) with respect to $`\rho `$ and $`𝐁^{}`$ in the vicinity of the constrained ground state $`(\rho ^{},𝐁^{})`$. The properties of the HF have been studied by various authors who found that, in fact, it often yields a better approximation of the ground state energy, in the vicinity of $`(\rho ^{},𝐁^{})`$, than the corresponding HKF Read1989 .

Starting from the ferromagnetic state (with $`𝐮_𝐑^0=𝐮^0`$ for all sites $`𝐑`$), for which we assume the self-consistent density $`\rho _0n_0\sigma _0+m_0𝐮_0𝝈`$ to be known, we perform some infinitesimal rotations: $`𝐮_𝐑=𝐮^0+\delta 𝐮_𝐑`$. For the ferromagnetic state, the constraining field vanishes everywhere. For the rotated state, we approximate the energy by using the HF $`_{\text{Harris}}[\rho ,𝐁^{}]`$, evaluated for a trial input density equal to $`\rho _{\text{in}}n_0\sigma _0+m_0(𝐮_\mathrm{𝟎}+\delta 𝐮)𝝈`$, and for some trial input constrain $`𝐁_{\text{in}}^{}`$ (to be specified later).

If we use the local density approximation (LDA) for the exchange correlation term (as is almost always done), the only term in our trial evaluation of the HF which varies with $`\delta 𝐮\{\delta 𝐮_𝐑\}`$ is the band energy (sum of eigenvalues $`\epsilon _i`$), so that the energy associated with the fluctuation $`\delta 𝐮`$ is

$`\mathrm{\Delta }E(\delta 𝐮)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_e}{}}}\left\{\epsilon _i(w_{\text{eff}}^{}[\rho _{\text{in}},𝐁_{\text{in}}^{}])\epsilon _i(w_{\text{eff}}^{}[\rho _0,0])\right\}`$ (3)
$`+𝒪_2(\delta n,\delta m,\delta 𝐁^{}),`$
where $`𝒪_p(x,y,\mathrm{})`$ is a quantity of order $`p`$ (and higher) in $`x,y,\mathrm{}`$, and $`\delta nn^{}n_0`$, $`\delta mm^{}m_0`$, $`\delta 𝐁^{}𝐁^{}𝐁_{\text{in}}^{}`$, respectively. We are looking for an expansion of $`\mathrm{\Delta }E(\delta 𝐮)`$ to second order in $`\delta 𝐮`$ of the form

$$\mathrm{\Delta }E(\delta 𝐮)=\underset{𝐑,𝐑^{}}{}A_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑\delta 𝐮_𝐑^{}+𝒪_4(\delta 𝐮),$$
(4)
with $`A_{\mathrm{𝐑𝐑}^{}}J_{\mathrm{𝐑𝐑}^{}}+\delta _{\mathrm{𝐑𝐑}^{}}\left(_{𝐑^{\prime \prime }}J_{\mathrm{𝐑𝐑}^{\prime \prime }}\right)`$, which defines the coupling parameters $`J_{\mathrm{𝐑𝐑}^{}}`$. The sum rule

$$\underset{𝐑}{}A_{\mathrm{𝐑𝐑}^{}}=\underset{𝐑^{}}{}A_{\mathrm{𝐑𝐑}^{}}=0$$
(5)
expresses the fact that the total energy remains invariant (within the non-relativistic theory) upon a uniform rotation of the magnetization. The definition (4) for $`A_{\mathrm{𝐑𝐑}^{}}`$ implies that it is related to the (static) transverse susceptibility $`\chi `$ by: $`2A_{\mathrm{𝐑𝐑}^{}}=M_𝐑\left(\chi ^1\right)_{\mathrm{𝐑𝐑}^{}}M_𝐑^{}`$.

The bare MFT of Liechtenstein et al. Liechtenstein1984  amounts to make the choice $`𝐁_{\text{in}}^{}0`$ for the trial input constraining fields. They then obtain

$$\mathrm{\Delta }E(\delta 𝐮)=\underset{𝐑,𝐑^{}}{}\stackrel{~}{A}_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑\delta 𝐮_𝐑^{}+𝒪_2(\delta n,\delta m,𝐁^{}),$$
(6)
with $`\stackrel{~}{A}_{\mathrm{𝐑𝐑}^{}}\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}}+\delta _{\mathrm{𝐑𝐑}^{}}\left(_{𝐑^{\prime \prime }}\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{\prime \prime }}\right)`$, and

$`\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}}`$ $``$ $`{\displaystyle \frac{1}{\pi }}\text{Im}{\displaystyle _{\mathrm{}}^{\epsilon _F}}d\epsilon {\displaystyle _{\mathrm{\Omega }_𝐑}}d𝐫{\displaystyle _{\mathrm{\Omega }_𝐑^{}}}d𝐫^{}B_{\text{xc}}(𝐫)G^{}(𝐫,𝐫^{})`$ (7)
$`\times B_{\text{xc}}(𝐫^{})G^{}(𝐫^{},𝐫).`$
Since $`\delta n`$ and $`\delta m`$ are even with respect to $`\delta 𝐮`$, they are generally of second order in $`\delta 𝐮`$; however, the constrains $`𝐁^{}`$ are odd with respect to $`\delta 𝐮`$, and therefore generally of first order in $`\delta 𝐮`$, so that $`𝒪_2(\delta n,\delta m,𝐁^{})=𝒪_2(\delta 𝐮)`$, although $`𝒪_2(\delta n,\delta m)=𝒪_4(\delta 𝐮)`$. As a consequence, the parameters $`\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}}`$ (hereafter called the bare exchange parameters) are not equal to the true exchange parameters $`J_{\mathrm{𝐑𝐑}^{}}`$ and entail some systematic error. Note, however, that this problem does not occur for the case considered by Oswald et al. Oswald1985  since the constrains vanish in the case they consider.

Clearly, to calculate correctly the true exchange parameters $`J_{\mathrm{𝐑𝐑}^{}}`$, we need to take the exact constrains $`𝐁^{}`$ (which are still unknown yet) as trial input in our estimate of the HF. If we do so, we obtain after some algebra,

$$\mathrm{\Delta }E(\delta 𝐮)=\underset{𝐑,𝐑^{}}{}\left[\stackrel{~}{A}_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑\delta 𝐮_𝐑^{}+(M_𝐑\delta _{\mathrm{𝐑𝐑}^{}}\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}})𝐁_𝐑^{}\delta 𝐮_𝐑^{}\frac{1}{2}\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}𝐁_𝐑^{}𝐁_𝐑^{}^{}\right]+𝒪_4(\delta 𝐮),$$
(8)
with $`M_𝐑_{\mathrm{\Omega }_𝐑}d𝐫m(𝐫)`$ is the magnetic moment of atom $`𝐑`$, and where the bare transverse susceptibility $`\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}`$ and the exchange-correlation response function $`\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}`$ are given by

$`\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}`$ $``$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _{\mathrm{}}^{\epsilon _F}}d\epsilon {\displaystyle _{\mathrm{\Omega }_𝐑}}d𝐫{\displaystyle _{\mathrm{\Omega }_𝐑^{}}}d𝐫^{}\text{Im}\left[G^{}(𝐫,𝐫^{})G^{}(𝐫^{},𝐫)\right],`$ (9)
$`\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}`$ $``$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^{\epsilon _F}}d\epsilon {\displaystyle _{\mathrm{\Omega }_𝐑}}d𝐫{\displaystyle _{\mathrm{\Omega }_𝐑^{}}}d𝐫^{}\text{Im}\left[G^{}(𝐫,𝐫^{})B_{\text{xc}}(𝐫^{})G^{}(𝐫^{},𝐫)+G^{}(𝐫,𝐫^{})B_{\text{xc}}(𝐫^{})G^{}(𝐫^{},𝐫)\right].`$ (10)
The constrains are obtained by expressing the transverse moments $`M_𝐑\delta 𝐮_𝐑`$ as resulting from the transverse exchange-correlation and constraining fields: $`M_𝐑\delta 𝐮_𝐑=_𝐑^{}\left(\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑^{}+\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}𝐁_𝐑^{}^{}\right).`$ In order to keep the expressions compact we introduce matrix notations: $`𝖠`$, $`\stackrel{~}{𝖠}`$, $`\stackrel{~}{𝖪}`$, $`\stackrel{~}{𝖷}`$, $`\stackrel{~}{𝖬}`$ are the matrices whose $`(\mathrm{𝐑𝐑}^{})`$-elements are, respectively, $`A_{\mathrm{𝐑𝐑}^{}}`$, $`\stackrel{~}{A}_{\mathrm{𝐑𝐑}^{}}`$, $`\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}`$, $`\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}`$, $`M_𝐑\delta _{\mathrm{𝐑𝐑}^{}}`$. Inserting the resulting expression of the constrains, $`𝐁_𝐑^{}=_𝐑^{}\left[𝖷^1(𝖬\stackrel{~}{𝖪})\right]_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑^{},`$ into the above expression of $`\mathrm{\Delta }E(\delta 𝐮)`$, we finally obtain the exact explicit expression of the renormalized exchange parameters:

$$𝖠=\stackrel{~}{𝖠}+\frac{1}{2}(𝖬\stackrel{~}{𝖪}^T)\stackrel{~}{𝖷}^1(𝖬\stackrel{~}{𝖪}).$$
(11)
One can easily prove the sum rule $`M_𝐑=_𝐑^{}\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}`$, which implies that the constrains $`𝐁_𝐑^{}`$ vanish for a coherent rotation of all magnetic moment, and that the sum rule (5) is indeed satisfied. The above result, Eq. (11), the “renormalized magnetic force theorem”, constitutes the main result of this Letter. Its importance is that it corrects the systematic error introduced by the “bare magnetic force theorem”, without increasing significantly the computational effort. It therefore provides a general method to calculate exactly the exchange parameters $`J_{\mathrm{𝐑𝐑}^{}}`$ which is several orders of magnitude faster than a fully self-consistent calculation.

The linearized equation of motion of the transverse fluctuations is

$$M_𝐑\frac{\mathrm{d}\delta 𝐮_𝐑}{\mathrm{d}t}=2\frac{\mathrm{\Delta }E}{\delta 𝐮_𝐑}\times 𝐮^0=4\underset{𝐑^{}}{}A_{\mathrm{𝐑𝐑}^{}}\delta 𝐮_𝐑^{}\times 𝐮^0,$$
(12)
and the spin-wave energies are given, as usual, by the eigenvalues of the symmetric matrix $`\mathrm{}\mathsf{\Omega }4𝖬^{1/2}\mathrm{𝖠𝖬}^{1/2}.`$

Some better physical insight into the nature of the above renormalization of the exchange parameters can be obtained if one performs a simple, but yet quite reasonable, approximation. Let us define

$$\mathrm{\Delta }_𝐑\frac{2}{M_𝐑}_{\mathrm{\Omega }_𝐑}d𝐫B_{\text{xc}}(𝐫)m(𝐫)=\frac{4}{M_𝐑}\underset{𝐑^{}}{}\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}},$$
(13)
which can be seen as some average of the exchange splitting on site $`𝐑`$. The second equality, in the above equation expresses a sum rule related to the invariance with respect to a global spin rotation. If the magnetization, within an atomic cell, is sufficiently rigid, i.e., if intra-atomic fluctuations of the spin-density have a large energy cost (as is usually the case in transition metals), one has approximately:

$$\stackrel{~}{K}_{\mathrm{𝐑𝐑}^{}}\frac{4\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}}}{\mathrm{\Delta }_𝐑},\text{ and }\stackrel{~}{\chi }_{\mathrm{𝐑𝐑}^{}}\frac{8\stackrel{~}{J}_{\mathrm{𝐑𝐑}^{}}}{\mathrm{\Delta }_𝐑\mathrm{\Delta }_𝐑^{}}.$$
(14)
Note that the above relations become exact if do not perform the discretization approximation $`𝐮(𝐫)𝐮_𝐑`$ (see footnote generalization\_continuum ). One then obtains the following simple and transparent renormalization for the exchange interaction parameters,

$$𝖠=\stackrel{~}{𝖠}\left(14𝖬^1\mathsf{\Delta }^1\stackrel{~}{𝖠}\right)^1,$$
(15)
and for the spin-wave matrix

$$\mathrm{}\mathsf{\Omega }=\mathrm{}\stackrel{~}{\mathsf{\Omega }}\left(1\mathsf{\Delta }^1\mathrm{}\stackrel{~}{\mathsf{\Omega }}\right)^1,$$
(16)
where $`(\mathsf{\Delta })_{\mathrm{𝐑𝐑}^{}}\mathrm{\Delta }_𝐑\delta _{\mathrm{𝐑𝐑}^{}}`$ and $`\mathrm{}\stackrel{~}{\mathsf{\Omega }}4𝖬^{1/2}\stackrel{~}{𝖠}𝖬^{1/2}`$ is the bare (unrenormalized) spin-wave matrix. Noting that the bare parameters $`\stackrel{~}{A}_{\mathrm{𝐑𝐑}^{}}`$ can be expressed in terms of the Stoner parameters $`\left(𝖨^{\text{xc}}\right)_{\mathrm{𝐑𝐑}^{}}\delta _{\mathrm{𝐑𝐑}^{}}\mathrm{\Delta }_𝐑/(2M_𝐑)`$ as $`2\stackrel{~}{𝖠}=\mathrm{𝖬𝖨}^{\text{xc}}\left(1\stackrel{~}{𝖷}𝖨^{\text{xc}}\right)𝖬`$, one easily shows that the renormalization (15) can be re-expressed as

$$2𝖠\mathrm{𝖬𝖷}^\mathrm{𝟣}𝖬=𝖬\left(\stackrel{~}{𝖷}^1𝖨^{\text{xc}}\right)𝖬.$$
(17)
The above expression has the familiar form of the random-phase approximation (RPA) result for the transverse susceptibility; it is important to realize, however, that in the present context (as shown by the above derivation) this result is exact within the LDA (except for the discretization approximation $`𝐮(𝐫)𝐮_𝐑`$, which as already indicated, can be removed easily). The present approach is therefore formally equivalent to approaches based on calculations of the transverse susceptibility Grotheer2001 , however without the need for self-consistent total energy calculations.

For the particular case of periodic lattices, the above equations are most conveniently solved in Fourier space. For systems with a single atom per unit cell, and using the approximation (14), the renormalization of the exchange parameters leads to the simple rescaling of the spin-wave spectrum

$$\mathrm{}\omega (𝐪)=\frac{\mathrm{}\stackrel{~}{\omega }(𝐪)}{1\mathrm{}\stackrel{~}{\omega }(𝐪)/\mathrm{\Delta }}.$$
(18)
The above result clearly shows that the bare MFT yields correct results only in the limit of long wavelength (in particular, it yields the correct spin-wave stiffness $`D`$), or if the spin-wave energies remain much smaller than the exchange splitting $`\mathrm{\Delta }`$. The above result is in agreement with the estimate of the error entailed by the bare MFT proposed by Grotheer Grotheer2001 : $`\mathrm{}\stackrel{~}{\omega }\mathrm{}\omega \mathrm{}\stackrel{~}{\omega }/(1\mathrm{}\stackrel{~}{\omega }/\mathrm{\Delta })^2`$.

The Curie temperature can be calculated by means of the random-phase-approximation (RPA) Green’s function method Pajda2001 . For periodic systems with a single atom per unit cell, one has

$$\frac{1}{k_BT_C^{\text{RPA}}}=\frac{6}{M}\frac{1}{N}\underset{𝐪}{}\frac{1}{\mathrm{}\omega (𝐪)},$$
(19)
so that, by using approximation (14), one therefore obtains an extremely simple renormalization of the (RPA) Curie temperature:

$$k_BT_C^{\text{RPA}}=k_B\stackrel{~}{T}_C^{\text{RPA}}\left(16\frac{k_B\stackrel{~}{T}_C^{\text{RPA}}}{M\mathrm{\Delta }}\right)^1,$$
(20)
where $`\stackrel{~}{T}_C^{\text{RPA}}`$ is the Curie temperature obtained from the bare exchange parameters. As seen from Table 1, the renormalization of the exchange parameters considerably improves the agreement between theoretical and experimental Curie temperatures of Fe and Ni.

The method discussed here provides a convenient and accurate approach to study the exchange interactions, spin-wave spectra and Curie temperature of complex systems such as disordered alloys, ultrathin films, nanostructures, dilute magnetic semiconductors (e.g., Ga<sub>1-x</sub>Mn<sub>x</sub>As), etc.

###### Acknowledgements.

I wish to express my thanks to Helmut Eschrig for some very stimulating discussions which motivated me to address this problem. I am grateful to Josef Kudrnovský, Leonid Sandratskii, and Valeri Stepanyuk for their helpful remarks on this work.