text:
# Casimir densities for two concentric spherical shells in the global monopole spacetime

## 1 Introduction

The Casimir effect is regarded as one of the most striking manifestation of vacuum fluctuations in quantum field theory. The presence of reflecting boundaries alters the zero-point modes of a quantized field, and results in the shifts in the vacuum expectation values of quantities quadratic in the field, such as the energy density and stresses. In particular, vacuum forces arise acting on constraining boundaries. The particular features of these forces depend on the nature of the quantum field, the type of spacetime manifold and its dimensionality, the boundary geometries and the specific boundary conditions imposed on the field. Since the original work by Casimir in 1948  many theoretical and experimental works have been done on this problem (see, e.g.,  and references therein). Many different approaches have been used: mode summation method with combination of the zeta function regularization technique, Green function formalism, multiple scattering expansions, heat-kernel series, etc. The Casimir effect can be viewed as a polarization of vacuum by boundary conditions. Another type of vacuum polarization arises in the case of external gravitational fields . In a previous paper  we have studied an example of a situation when both types of sources for the polarization are present. Namely, we have investigated the vacuum expectation values of the square of a scalar field and energy-momentum tensor induced by a spherical shell in the spacetime of a point-like global monopole.

Topological defects have attracted a great deal of attention because of their relevance to a number of different areas ranging from condensed matter to structure formation (for a review see ). In the context of hot big bang cosmology, the unified theories of the fundamental interactions predict that the universe passes through a sequence of phase transitions. These phase transitions might have given rise to several kinds of topological defects depending on the nature of the symmetry that is broken . If a global SO(3) symmetry of a triplet scalar field is broken, the point like defects called global monopoles are believed to be formed. The simplified global monopole was introduced by Sokolov and Starobinsky . The gravitational effects of the global monopole were studied in Ref. , where a solution is presented which describes a global monopole at large radial distances. The quantum vacuum effects of the matter fields on the global monopole background have been considered in . The effects produced by the non-zero temperature are investigated as well . The zeta function and the Casimir energy for spherical boundaries in this background are considered in .

In this paper, we continue the investigation started in  and consider the Casimir densities in the region between two concentric spherical shells on background of the $`D+1`$-dimensional spacetime of a point-like global monopole. The positive frequency Wightman function, the vacuum expectation values of the field square and energy-momentum tensor are investigated for a massive scalar field with general curvature coupling parameter $`\xi `$. In addition to describing the physical structure of the quantum field at a given point, the energy-momentum tensor acts as the source of gravity in the Einstein equations. It therefore plays an important role in modeling a self-consistent dynamics involving the gravitational field . We study the general case of Robin boundary conditions with different coefficients for the inner and outer spheres. The Robin boundary conditions are an extension of the ones imposed on perfectly conducting boundaries and may, in some geometries, be useful for depicting the finite penetration of the field into the boundary with the ”skin-depth” parameter related to the Robin coefficient . It is interesting to note that the quantum scalar field satisfying Robin condition on the boundary of a cavity violates the Bekenstein’s entropy-to-energy bound near certain points in the space of the parameter defining the boundary condition . This type of conditions also appear in considerations of the vacuum effects for a confined charged scalar field in external fields  and in quantum gravity . Mixed boundary conditions naturally arise for scalar and fermion bulk fields in the Randall-Sundrum model . In this model the bulk geometry is a slice of anti-de Sitter space and the corresponding Robin coefficient is related to the curvature scale of this space. For scalars with general curvature coupling the essential point is the relation between the mode sum energy, evaluated as a renormalized sum of the zero-point energies for each normal mode of frequency, and the volume integral of the renormalized energy density. In  it has been shown that in the discussion of this question for the Robin parallel plates it is necessary to include in the energy a surface term concentrated on the boundary (see  for similar issues in more complicated geometries and the discussion of the paper  in ).

We have organized the paper as follows. In the next section we derive a formula for the Wightman function in the region between two spheres. The reason for our choice of the Wightman function is that this function also determines the response of the particle detectors in a given state of motion. Following Refs. , to evaluate the bilinear field products we use the mode sum method in combination with the summation formulae from  (see also ). These formulae allow (i) to extract from vacuum expectation values the parts due to a single sphere on the global monopole background, and (ii) to present the ”interference” parts in terms of exponentially convergent integrals involving the modified Bessel functions. The vacuum expectation value of the field square is obtained computing the Wightman function in the coincidence limit, and is investigated in section 3. The various limiting cases are considered including the limit of the strong gravitational field corresponding to small values of the parameter associated with the solid angle deficit in the global monopole spacetime. Section 4 is devoted to the vacuum expectation values of the energy density and stresses. The formulae are derived for the interaction forces between the spheres. The behavior of these quantities is investigated in various limiting cases. Section 5 concludes the main results of the paper.

## 2 Wightman function

In this paper we consider a real scalar field $`\phi `$ with curvature coupling parameter $`\xi `$ in the $`D+1`$-dimensional spacetime of a point-like global monopole. In the hyperspherical polar coordinates $`(r,\vartheta ,\varphi )(r,\theta _1,\theta _2,\mathrm{}\theta _n,\varphi )`$, $`n=D2`$, the corresponding geometry is described by the line element

$$ds^2=dt^2dr^2\sigma ^2r^2d\mathrm{\Omega }_D^2,$$
(1)
where $`d\mathrm{\Omega }_D^2`$ is the line element on the surface of the unit sphere in $`D`$-dimensional Euclidean space, the parameter $`\sigma `$ is smaller than unity and is related to the symmetry breaking energy scale in the theory. In the spacetime given by line element (1) the solid angle deficit is $`(1\sigma ^2)S_D`$, with $`S_D=2\pi ^{D/2}/\mathrm{\Gamma }(D/2)`$ being the total area of the surface of the unit sphere in $`D`$-dimensional Euclidean space. The field equation has the form

$$\left(_i^i+m^2+\xi R\right)\phi =0,$$
(2)
where $`m`$ is the mass for the field quanta, $`_i`$ is the covariant derivative operator associated with the metric given by line element (1), and

$$R=n(n+1)\frac{\sigma ^21}{\sigma ^2r^2}$$
(3)
is the corresponding Ricci scalar (we adopt the convention of Birrell and Davies  for the curvature tensor). Note that for $`\sigma 1`$ the geometry is singular at the origin (point-like monopole), $`r=0`$. In (2) the values of the curvature coupling parameter $`\xi =0`$, and $`\xi =\xi _D`$ with $`\xi _D(D1)/4D`$ correspond to the minimal and conformal couplings, respectively.

In this paper we are interested in the vacuum expectation values (VEVs) of the field bilinear products on background of the geometry described by (1), assuming that the field satisfies the Robin boundary conditions

$$\left(\stackrel{~}{A}_r+\stackrel{~}{B}_r\frac{}{r}\right)\phi (x)=0,r=a,b,$$
(4)
on two spheres with radii $`a`$ and $`b`$, $`a<b`$, concentric with the monopole. Here the coefficients $`\stackrel{~}{A}_r`$ and $`\stackrel{~}{B}_r`$ are constants, in general, different for the inner and outer spheres. The imposition of this boundary condition on the quantum field $`\phi (x)`$ leads to the modification of the spectrum for the zero-point fluctuations and results in the shift in VEVs for physical quantities. In particular, vacuum forces arise acting on constraining boundary. This is the familiar Casimir effect. The VEVs for the physical quantities bilinear in the field can be evaluated if the corresponding Wightman function is known. For this reason we first concentrate on the positive frequency Wightman function $`G^+(x,x^{})=0|\phi (x)\phi (x^{})|0`$, where $`|0`$ is the amplitude for the corresponding vacuum state. The Wightman function also determines the response of the particle detectors in a given state of motion . By expanding the field operator over eigenfunctions and using the commutation rules one can see that

$$0|\phi (x)\phi (x^{})|0=\underset{\alpha }{}\phi _\alpha (x)\phi _\alpha ^{}(x^{}),$$
(5)
where $`\{\phi _\alpha (x),\phi _\alpha ^{}(x^{})\}`$ is a complete orthonormal set of positive and negative frequency solutions to the field equation with quantum numbers $`\alpha `$, satisfying boundary conditions (4). Note that for $`D=1`$ we have the standard Casimir-like geometry on Minkowski spacetime.

In the hyperspherical coordinates, for the region between two spheres the complete set of solutions to (2) with scalar curvature from (3), has the form

$$\phi _\alpha (x)=\beta _\alpha r^{n/2}g_{\nu _l}(\lambda a,\lambda r)Y(m_k;\vartheta ,\varphi )e^{i\omega t},l=0,1,2,\mathrm{},$$
(6)
where $`m_k=(m_0l,m_1,\mathrm{},m_n)`$, and $`m_1,m_2,\mathrm{},m_n`$ are integers such that

$$0m_{n1}m_{n2}\mathrm{}m_1l,m_{n1}m_nm_{n1},$$
(7)
$`Y(m_k;\vartheta ,\varphi )`$ is the surface harmonic of degree $`l`$ (see ). In (6)

$$g_{\nu _l}(\lambda a,\lambda r)J_{\nu _l}(\lambda r)\overline{Y}_{\nu _l}^{(a)}(\lambda a)\overline{J}_{\nu _l}^{(a)}(\lambda a)Y_{\nu _l}(\lambda r),\lambda =\sqrt{\omega ^2m^2},$$
(8)
$`J_\nu (z)`$ and $`Y_\nu (z)`$ are the Bessel and Neumann functions, and the functions with overbars are defined in accordance with

$$\overline{F}^{(\alpha )}(z)A_\alpha F(z)+B_\alpha zF^{}(z),A_\alpha =\stackrel{~}{A}_\alpha \stackrel{~}{B}_\alpha n/2\alpha ,B_\alpha =\stackrel{~}{B}_\alpha /\alpha ,\alpha =a,b.$$
(9)
Here and below we use the following notation

$$\nu _l=\frac{1}{\sigma }\left[\left(l+\frac{n}{2}\right)^2+(1\sigma ^2)n(n+1)\left(\xi \xi _{D1}\right)\right]^{1/2},$$
(10)
assuming that $`\nu _l^2`$ is non-negative. This corresponds to the restriction on the values of the curvature coupling parameter for $`n>0`$, given by the condition

$$\xi \frac{n\sigma ^2}{4(n+1)(1\sigma ^2)}.$$
(11)
This condition is satisfied by the most important special cases of the minimal and conformal couplings. The coefficients $`\beta _\alpha `$ in (6) can be found from the normalization condition

$$\left|\phi _\alpha (x)\right|^2\sqrt{g}𝑑V=\frac{1}{2\omega },$$
(12)
where the integration goes over the region between the spheres, $`arb`$. Substituting eigenfunctions (6), and using the relation

$$\left|Y(m_k;\vartheta ,\varphi )\right|^2𝑑\mathrm{\Omega }=N(m_k)$$
(13)
(the explicit form for $`N(m_k)`$ is given in  and will not be necessary for the following considerations in this paper) for the spherical harmonics and the value for the standard integral involving the square of a cylinder function , one finds

$$\beta _\alpha ^2=\frac{\pi ^2\lambda T_{\nu _l}^{ab}(b/a,\lambda a)}{4N(m_k)\omega a\sigma ^{D1}},$$
(14)
where we use the notation

$$T_\nu ^{ab}(\eta ,z)=z\left\{\frac{\overline{J}_\nu ^{(a)2}(z)}{\overline{J}_\nu ^{(b)2}(\eta z)}\left[A_b^2+B_b^2(\eta ^2z^2\nu ^2)\right]A_a^2B_a^2(z^2\nu ^2)\right\}^1,\eta =\frac{b}{a}.$$
(15)
The functions chosen in the form (8) satisfy the boundary condition on the sphere $`r=a`$. From the boundary condition on $`r=b`$ one obtains that the corresponding eigenmodes are solutions to the equation

$$C_{\nu _l}^{ab}(b/a,\lambda a)\overline{J}_\nu ^{(a)}(\lambda a)\overline{Y}_\nu ^{(b)}(\lambda b)\overline{J}_\nu ^{(b)}(\lambda b)\overline{Y}_\nu ^{(a)}(\lambda a)=0.$$
(16)
Below the roots to this equation will be denoted by $`\gamma _{\nu ,k}=\lambda a`$ , $`k=1,2,\mathrm{}`$. The corresponding eigenfrequencies $`\omega =\omega _{\nu _l,k}`$ are related to these zeros as $`\omega _{\nu _l,k}=\sqrt{\lambda _{\nu _l,k}^2/a^2+m^2}`$.

Substituting the eigenfunctions into the mode sum (5) and using the addition formula

$$\underset{m_k}{}\frac{1}{N(m_k)}Y(m_k;\vartheta ,\varphi )Y^{}(m_k;\vartheta ^{},\varphi ^{})=\frac{2l+n}{nS_D}C_l^{n/2}(\mathrm{cos}\theta ),$$
(17)
for the expectation value of the field product one finds

$$0|\phi (x)\phi (x^{})|0=\frac{\pi ^2(rr^{})^{n/2}}{4naS_D\sigma ^{D1}}\underset{l=0}{\overset{\mathrm{}}{}}(2l+n)C_l^{n/2}(\mathrm{cos}\theta )\underset{k=1}{\overset{\mathrm{}}{}}h(\gamma _{\nu _l,k})T_{\nu _l}^{ab}(b/a,\gamma _{\nu _l,k}),$$
(18)
with the function

$$h(z)=\frac{z}{\sqrt{z^2+m^2a^2}}g_{\nu _l}(z,zr/a)g_{\nu _l}(z,zr^{}/a)e^{i\sqrt{z^2/a^2+m^2}(t^{}t)},$$
(19)
and $`C_l^p(x)`$ is the Gegenbauer or ultraspherical polynomial of degree $`l`$ and order $`p`$. To sum over $`k`$ we will use the generalized Abel-Plana summation formula

$`{\displaystyle \frac{\pi ^2}{2}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}h(\gamma _{\nu ,k})T_\nu ^{ab}(\eta ,\gamma _{\nu ,k})`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{h(x)dx}{\overline{J}_\nu ^{(a)2}(x)+\overline{Y}_\nu ^{(a)2}(x)}}`$
$`{\displaystyle \frac{\pi }{4}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\overline{K}_\nu ^{(b)}(\eta x)}{\overline{K}_\nu ^{(a)}(x)}}{\displaystyle \frac{\left[h(xe^{\pi i/2})+h(xe^{\pi i/2})\right]dx}{\overline{K}_\nu ^{(a)}(x)\overline{I}_\nu ^{(b)}(\eta x)\overline{K}_\nu ^{(b)}(\eta x)\overline{I}_\nu ^{(a)}(x)}},`$
where $`I_\nu (x)`$ and $`K_\nu (x)`$ are the Bessel modified functions. The corresponding conditions for the function (19) are satisfied if $`r+r^{}+|tt^{}|<2b`$. Note that this is the case in the coincidence limit for the region under consideration. Note that we have assumed values $`A_\alpha `$ and $`B_\alpha `$ for which all zeros for (16) are real. In case of existence of purely imaginary zeros we have to include additional residue terms on the left of this formula (see ). Applying to the sum over $`k`$ in (18) formula (2), for the corresponding Wightman function one obtains

$$0|\phi (x)\phi (x^{})|0=\frac{\sigma ^{1D}}{2naS_D}\underset{l=0}{\overset{\mathrm{}}{}}\frac{2l+n}{(rr^{})^{n/2}}C_l^{n/2}(\mathrm{cos}\theta )\{_0^{\mathrm{}}\frac{h(z)dz}{\overline{J}_{\nu _l}^{(a)2}(z)+\overline{Y}_{\nu _l}^{(a)2}(z)}$$
(21)
$$\frac{2}{\pi }_{ma}^{\mathrm{}}\frac{zdz}{\sqrt{z^2a^2m^2}}\frac{\overline{K}_{\nu _l}^{(b)}(\eta z)}{\overline{K}_{\nu _l}^{(a)}(z)}\frac{G_{\nu _l}^{(a)}(z,zr/a)G_{\nu _l}^{(a)}(z,zr^{}/a)}{\overline{K}_{\nu _l}^{(a)}(z)\overline{I}_{\nu _l}^{(b)}(\eta z)\overline{K}_{\nu _l}^{(b)}(\eta z)\overline{I}_{\nu _l}^{(a)}(z)}\mathrm{cosh}\left[\sqrt{z^2/a^2m^2}(t^{}t)\right]\},$$

where we have introduced notations

$$G_\nu ^{(\alpha )}(z,y)=I_\nu (y)\overline{K}_\nu ^{(\alpha )}(z)\overline{I}_\nu ^{(\alpha )}(z)K_\nu (y),\alpha =a,b.$$
(22)
In the limit $`b\mathrm{}`$ the second integral on the right of (21) tends to zero (for large $`b/a`$ the subintegrand is proportional to $`e^{2bz/a}`$), whereas the first one does not depend on $`b`$. It follows from here that the term with the first integral in the figure braces corresponds to the Wightman function for the region outside a single sphere with radius $`a`$ on background of the global monopole geometry. As a result the Wightman function in the region between two spheres is presented in the form

$`0|\phi (x)\phi (x^{})|0`$ $`=`$ $`0^{(a)}|\phi (x)\phi (x^{})|0^{(a)}{\displaystyle \frac{\sigma ^{1D}}{\pi nS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2l+n}{(rr^{})^{n/2}}}C_l^{n/2}(\mathrm{cos}\theta )`$
$`\times `$ $`{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z\mathrm{\Omega }_{a\nu _l}(az,bz)}{\sqrt{z^2m^2}}}G_{\nu _l}^{(a)}(az,rz)G_{\nu _l}^{(a)}(az,r^{}z)\mathrm{cosh}\left[\sqrt{z^2m^2}(t^{}t)\right],`$
where

$$\mathrm{\Omega }_{a\nu }(az,bz)=\frac{\overline{K}_\nu ^{(b)}(bz)/\overline{K}_\nu ^{(a)}(az)}{\overline{K}_\nu ^{(a)}(az)\overline{I}_\nu ^{(b)}(bz)\overline{K}_\nu ^{(b)}(bz)\overline{I}_\nu ^{(a)}(az)},$$
(24)
and $`|0^{(a)}`$ is the amplitude for the vacuum state in the case of a single sphere with radius $`a`$. The Wightman function $`0^{(a)}|\phi (x)\phi (x^{})|0^{(a)}`$ is investigated in a previous paper , and below we will mainly concentrate on the terms induced by the presence of the second sphere. Note that in the coincidence limit, $`x^{}=x`$, the second summand on the right hand side of (2) will give a finite result for $`ar<b`$, and is divergent on the boundary $`r=b`$. In Ref.  the Wightman function for a single sphere in the global monopole spacetime is presented in the form

$$0^{(a)}|\phi (x)\phi (x^{})|0^{(a)}=0_\mathrm{m}|\phi (x)\phi (x^{})|0_\mathrm{m}+\phi (x)\phi (x^{})_b^{(a)},$$
(25)
where $`|0_\mathrm{m}`$ is the amplitude for the vacuum state in the case of the boundary-free global monopole geometry. The expressions for the boundary induced part $`\phi (x)\phi (x^{})_b^{(a)}`$ for both regions inside and outside a single shell are given in .

Using the identity

$`\mathrm{\Omega }_{\alpha \nu }(az,bz)G_\nu ^{(\alpha )}(\alpha z,rz)G_\nu ^{(\alpha )}(\alpha z,r^{}z)|_{\alpha =a}^{\alpha =b}`$ $`=`$ $`{\displaystyle \frac{\overline{I}_\nu ^{(a)}(az)}{\overline{K}_\nu ^{(a)}(az)}}K_\nu (rz)K_\nu (r^{}z)`$ (26)
$`{\displaystyle \frac{\overline{K}_\nu ^{(b)}(bz)}{\overline{I}_\nu ^{(b)}(bz)}}I_\nu (rz)I_\nu (r^{}z),`$
it can be seen that for the case of two spheres the Wightman function in the intermediate region can also be presented in the form

$`0|\phi (x)\phi (x^{})|0`$ $`=`$ $`0^{(b)}|\phi (x)\phi (x^{})|0^{(b)}{\displaystyle \frac{\sigma ^{1D}}{\pi nS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2l+n}{(rr^{})^{n/2}}}C_l^{n/2}(\mathrm{cos}\theta )`$
$`\times `$ $`{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z\mathrm{\Omega }_{b\nu _l}(az,bz)}{\sqrt{z^2m^2}}}G_{\nu _l}^{(b)}(bz,rz)G_{\nu _l}^{(b)}(bz,r^{}z)\mathrm{cosh}\left[\sqrt{z^2m^2}(t^{}t)\right],`$
with $`0^{(b)}|\phi (x)\phi (x^{})|0^{(b)}`$ being the Wightman function for the vacuum inside a single sphere with radius $`b`$, and

$$\mathrm{\Omega }_{b\nu }(az,bz)=\frac{\overline{I}_\nu ^{(a)}(az)/\overline{I}_\nu ^{(b)}(bz)}{\overline{K}_\nu ^{(a)}(az)\overline{I}_\nu ^{(b)}(bz)\overline{K}_\nu ^{(b)}(bz)\overline{I}_\nu ^{(a)}(az)}.$$
(28)
Note that formula (2) can be also derived by the procedure described above for (2), if in expression (6) for the eigenfunctions we replace the function $`g_{\nu _l}(\lambda a,\lambda r)`$ (given by (8)) by the function $`J_{\nu _l}(\lambda r)\overline{Y}_{\nu _l}^{(b)}(\lambda b)\overline{J}_{\nu _l}^{(b)}(\lambda b)Y_{\nu _l}(\lambda r)`$. In the coincidence limit, the second summand on the right of formula (2) is finite for $`a<rb`$ and diverges on the boundary $`r=a`$. It follows from here that if we write the regularized Wightman function in the form

$$0|\phi (x)\phi (x^{})|0=0_\mathrm{m}|\phi (x)\phi (x^{})|0_\mathrm{m}+\phi (x)\phi (x^{})_b^{(a)}+\phi (x)\phi (x^{})_b^{(b)}+\mathrm{\Delta }W(x,x^{}),$$
(29)
then in the coincidence limit the ”interference” term $`\mathrm{\Delta }W(x,x^{})`$ is finite for all values $`arb`$. Using the formula for the Wightman function for a single sphere (see ) and equation (2), it can be seen that this term may be presented as

$$\mathrm{\Delta }W(x,x^{})=\frac{\sigma ^{1D}}{\pi nS_D}\underset{l=0}{\overset{\mathrm{}}{}}\frac{2l+n}{(rr^{})^{n/2}}C_l^{n/2}(\mathrm{cos}\theta )_m^{\mathrm{}}\frac{zdz}{\sqrt{z^2m^2}}W^{(ab)}(r,r^{})\mathrm{cosh}\left[\sqrt{z^2m^2}(t^{}t)\right],$$
(30)
where

$$W^{(ab)}(r,r^{})=\frac{\overline{I}_{\nu _l}^{(a)}(az)}{\overline{I}_{\nu _l}^{(b)}(bz)}\frac{\overline{K}_{\nu _l}^{(b)}(bz)}{\overline{K}_{\nu _l}^{(a)}(az)}\left[\frac{G_{\nu _l}^{(a)}(az,zr)G_{\nu _l}^{(b)}(bz,zr^{})}{\overline{K}_{\nu _l}^{(a)}(az)\overline{I}_{\nu _l}^{(b)}(bz)\overline{K}_{\nu _l}^{(b)}(bz)\overline{I}_{\nu _l}^{(a)}(az)}I_{\nu _l}(zr^{})K_{\nu _l}(zr)\right].$$
(31)
Formula (29) presents the Wightman function in the region between two spheres, $`arb`$. In the regions $`ra`$ and $`rb`$, the Wightman functions are given by the formulae corresponding to a single sphere with radius $`a`$ and $`b`$ respectively, and are investigated in .

## 3 Vacuum expectation value for the field square

The VEV of the field square is obtained computing the Wightman function in the coincidence limit $`x^{}x`$. In the region between the spheres, from (2) and (2) we obtain two equivalent forms

$$0|\phi ^2(x)|0=0^{(\alpha )}|\phi ^2(x)|0^{(\alpha )}\frac{\sigma ^{1D}}{\pi nS_Dr^n}\underset{l=0}{\overset{\mathrm{}}{}}D_l_m^{\mathrm{}}𝑑z\frac{z\mathrm{\Omega }_{\alpha \nu _l}(az,bz)}{\sqrt{z^2m^2}}G_{\nu _l}^{(\alpha )2}(\alpha z,rz),$$
(32)
with $`\alpha =a,b`$, and

$$D_l=(2l+D2)\frac{\mathrm{\Gamma }(l+D2)}{\mathrm{\Gamma }(D1)l!}$$
(33)
being the degeneracy of each angular mode with given $`l`$, and $`\mathrm{\Gamma }(x)`$ is the gamma function. Using formula (29), the VEV for the field square can be also presented in the form

$$0|\phi ^2(x)|0=0_\mathrm{m}|\phi ^2(x)|0_\mathrm{m}+\phi ^2(x)_b^{(a)}+\phi ^2(x)_b^{(b)}+\phi ^2(x)^{(ab)},$$
(34)
where $`\phi ^2(x)_b^{(\alpha )}`$ is the VEV induced by a single sphere with radius $`\alpha `$, the ”interference” term is given by the formula

$`\phi ^2(x)^{(ab)}`$ $`=`$ $`{\displaystyle \frac{\sigma ^{1D}}{\pi nS_Dbr^n}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{zdz}{\sqrt{z^2m^2}}}{\displaystyle \frac{\overline{I}_{\nu _l}^{(a)}(az)\overline{K}_{\nu _l}^{(b)}(bz)}{\overline{K}_{\nu _l}^{(a)}(az)\overline{I}_{\nu _l}^{(b)}(bz)\overline{K}_{\nu _l}^{(b)}(bz)\overline{I}_{\nu _l}^{(a)}(az)}}`$ (35)
$`\times \left[{\displaystyle \frac{I_{\nu _l}(rz)}{\overline{I}_{\nu _l}^{(b)}(bz)}}G_{\nu _l}^{(b)}(bz,rz){\displaystyle \frac{K_{\nu _l}(rz)}{\overline{K}_{\nu _l}^{(a)}(az)}}G_{\nu _l}^{(a)}(az,rz)\right],`$
and is finite for all values $`arb`$. In the case of the Dirichlet scalar, by using the relation $`I_\nu (x)K_\nu (y)>I_\nu (y)K_\nu (x)`$ for $`y>x`$, it can be seen that $`\phi ^2(x)^{(ab)}>0`$. Note that for this boundary condition the both boundary induced terms $`\phi ^2(x)_b^{(\alpha )}`$, $`\alpha =a,b`$, inside and outside of a single spherical shell are negative .

Now let us consider the limiting cases of formula (35). In the limit $`a0`$ the subintegrand behaves as $`a^{2\nu _l}`$, and the dominant contribution to $`\phi ^2(x)^{(ab)}`$ comes from the summand with $`l=0`$. Using the standard formulae for the Bessel modified functions , noting that $`B_a=\stackrel{~}{B}_a/a`$, and assuming $`\nu _0>0`$, to the leading order one has

$`\phi ^2(x)^{(ab)}`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}}{\pi S_Dr^n\nu _0\mathrm{\Gamma }^2(\nu _0)}}{\displaystyle \frac{n+2\eta _a\nu _0}{n+2\nu _0}}\left({\displaystyle \frac{a}{2}}\right)^{2\nu _0}{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{z^{2\nu _0+1}dz}{\sqrt{z^2m^2}}}{\displaystyle \frac{\overline{K}_{\nu _0}^{(b)}(bz)}{\overline{I}_{\nu _0}^{(b)}(bz)}}`$ (36)
$`\times I_{\nu _0}(rz)[{\displaystyle \frac{\overline{K}_{\nu _0}^{(b)}(bz)}{\overline{I}_{\nu _0}^{(b)}(bz)}}I_{\nu _0}(rz)2K_{\nu _0}(rz)],a0,`$
where we use the notation

$$\eta _\alpha =\{\begin{array}{cc}\hfill 1,& \mathrm{for}\stackrel{~}{B}_\alpha =0\hfill \\ \hfill 1,& \mathrm{for}\stackrel{~}{B}_\alpha 0\hfill \end{array}.$$
(37)
For $`\nu _0=0`$ and $`a0`$, the dominant contribution behaves as $`\phi ^2(x)^{(ab)}1/\mathrm{ln}a`$. Note that for a minimally coupled scalar ($`\xi =0`$) one has $`\nu _0=D/21`$.

In the case of a massless scalar the asymptotic behavior of ”interference” part (35) for large values $`b`$ and fixed $`a`$ and $`r`$ can be obtained by changing the integration variable to $`y=zb`$ and expanding the subintegrand in terms of $`a/b`$ and $`r/b`$. The dominant contribution for the summand with a given $`l`$ has an order $`(a/b)^{2\nu _l}`$ (assuming that $`A_aB_a\nu _l`$) and the main contribution comes from the $`l=0`$ term. The leading term for the corresponding asymptotic expansion can be presented in the form

$`\phi ^2^{(ab)}`$ $``$ $`{\displaystyle \frac{2\sigma ^{1D}r^n}{\pi S_Da\mathrm{\Gamma }^2(\nu _0+1)}}\left({\displaystyle \frac{a}{2b}}\right)^{2\nu _0+1}{\displaystyle \frac{A_a+B_a\nu _0}{A_aB_a\nu _0}}\left[{\displaystyle \frac{A_a+B_a\nu _0}{A_aB_a\nu _0}}\left({\displaystyle \frac{a}{r}}\right)^{2\nu _0}2\right]`$ (38)
$`\times {\displaystyle _0^{\mathrm{}}}dzz^{2\nu _0}{\displaystyle \frac{nK_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}zK_{\nu _0}^{}(z)}{nI_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}zI_{\nu _0}^{}(z)}}.`$
Hence, the ”interference” part of the VEV of the square of the field operator vanishes in both limits $`a0`$ and $`b\mathrm{}`$.

Now we turn to the limit $`a,b\mathrm{}`$ for fixed $`ba`$ and $`\sigma `$. In this limit expression (35) diverges and the main contribution comes from large values of $`l\sigma [2(1a/b)]^1`$. Introducing in (35) a new integration variable $`y=bz/\nu _l`$, we can replace the Bessel modified functions by their uniform asymptotic expansions for large values of the order . To the leading order this gives

$$\phi ^2(x)^{(ab)}\frac{(b\sigma )^{1D}}{\pi S_D\mathrm{\Gamma }(D1)}\underset{l=0}{\overset{\mathrm{}}{}}l^{D2}_{m_l}^{\mathrm{}}𝑑y\frac{(y^2m_l^2)^{1/2}}{c_b(y)e^{2y}/c_a(y)1}F(y,r),$$
(39)
where we have introduced the notations

$$m_l=(ba)\sqrt{m^2+(l/b\sigma )^2},c_\alpha (y)=\frac{\stackrel{~}{A}_\alpha +y\stackrel{~}{B}_\alpha /(ba)}{\stackrel{~}{A}_\alpha y\stackrel{~}{B}_\alpha /(ba)},\alpha =a,b,$$
(40)
and

$$F(y,r)=\frac{1}{c_b(y)}\mathrm{exp}\left(2y\frac{rb}{ba}\right)+c_a(y)\mathrm{exp}\left(2y\frac{ar}{ba}\right)2.$$
(41)
In the limit under consideration in (39) we can make the replacement

$$\underset{l}{}l^{D2}f\left(\frac{l(ba)}{b\sigma }\right)\left(\frac{b\sigma }{ba}\right)^{D1}_0^{\mathrm{}}𝑑tt^{D2}f(t).$$
(42)
Further, introducing instead of $`y`$ a new integration variable $`u=\sqrt{y^2t^2m_l^2}`$ and converting to polar coordinates on the plane $`(u,t)`$, the angular part of the resulting integral is easily evaluated. Using the standard relations for the gamma function, one receives

$$\phi ^2(x)^{(ab)}\frac{(ba)^{1D}}{(4\pi )^{D/2}\mathrm{\Gamma }(D/2)}_{m_0}^{\mathrm{}}𝑑y\frac{(y^2m_0^2)^{D/21}}{c_b(y)e^{2y}/c_a(y)1}F(y,r),a,b\mathrm{},ba=\mathrm{const}.$$
(43)
This leading term of the corresponding asymptotic expansion does not depend on the parameter $`\sigma `$ and coincides with the corresponding quantity for two parallel plates with the separation $`ba`$, on background of the Minkowski spacetime.

And finally, consider the ”interference” term $`\phi ^2^{(ab)}`$ in the limit of strong gravitational field, $`\sigma 1`$, for fixed values $`a,b,r`$. For $`\xi >0`$, from (10) one has $`\nu _l1`$, and after introducing in (35) a new integration variable $`y=z/\nu _l`$, we can replace the modified Bessel function by their uniform asymptotic expansions for large values of the order. The integral over $`y`$ can be estimated by using the Laplace method. The main contribution to the sum over $`l`$ comes from the summand with $`l=0`$, and to the leading order we receive

$$\phi ^2(x)^{(ab)}\frac{\sigma ^{1D}\eta _a\eta _bexp[2\stackrel{~}{\nu }\mathrm{ln}(b/a)]}{(2\pi \stackrel{~}{\nu })^{1/2}r^nS_D\sqrt{b^2a^2}},\stackrel{~}{\nu }=\frac{1}{\sigma }\sqrt{n(n+1)\xi },\sigma 1.$$
(44)
Note that, using the corresponding expressions for single sphere parts , we can see that in the limit under consideration

$$\frac{\phi ^2(x)^{(ab)}}{\phi ^2(x)_b^{(\alpha )}}2\sqrt{\frac{|\alpha ^2r^2|}{b^2a^2}}\frac{\eta _a\eta _b}{\eta _\alpha }exp[2\stackrel{~}{\nu }|\mathrm{ln}(\alpha /r)|],\alpha =a,b,$$
(45)
and the ”interference” term is suppressed compared to the single sphere contribution. For $`\xi =0`$ and $`\sigma 1`$ for the terms with $`l0`$ one has $`\nu _l1`$. The corresponding contribution can be estimated by the way similar to that in the previous case. This contribution is exponentially small. For the summand with $`l=0`$ to the leading order over $`\sigma `$ we have $`\nu _l=n/2`$ in (35), and $`\phi ^2\sigma ^{1D}`$. As we see, the behavior of the ”interference” part in the VEV of the field square in the limit of the strong gravitational field is essentially different for minimally and and non-minimally coupled scalars. This behavior is clearly seen from figure 1 where we have plotted the dependence of the ”interference” term $`\phi ^2(x)^{(ab)}`$ on the ratio $`r/b`$ in the cases of conformally ($`\xi =\xi _D`$, left panel) and minimally ($`\xi =0`$, right panel) coupled massless Dirichlet scalars in $`D=3`$ dimensions for $`a/b=0.5`$. The separate curves correspond to the values $`\sigma =1`$ (a), $`\sigma =0.4`$ (b), $`\sigma =0.2`$ (c).

## 4 Vacuum expectation values of the energy-momentum tensor and the interaction forces between the spheres

In this section we will consider the VEVs for the energy momentum tensor operator in the region between two spheres on background of the global monopole spacetime. Using the standard classical expression, we can write the vacuum energy-momentum tensor as the coincidence limit

$$0|T_{ik}(x)|0=\underset{x^{}x}{lim}_i_k^{}0|\phi (x)\phi (x^{})|0+\left[\left(\xi \frac{1}{4}\right)g_{ik}_l^l\xi _i_k\xi R_{ik}\right]0|\phi ^2(x)|0,$$
(46)
where for the point-like global monopole spacetime the nonzero components of the Ricci tensor are given by expressions

$$R_2^2=R_3^3=\mathrm{}=R_D^D=n\frac{\sigma ^21}{\sigma ^2r^2},$$
(47)
with the indices $`2,3,\mathrm{},D`$ corresponding to the coordinates $`\theta _1,\theta _2,\mathrm{},\varphi `$ respectively. Substituting the Wightman function (2) into (46), we obtain that the vacuum energy-momentum tensor has the diagonal form (as expected by the symmetry of the model)

$$0|T_i^k|0=\mathrm{diag}(\epsilon ,p,p_{},\mathrm{},p_{}),$$
(48)
where the vacuum energy density $`\epsilon `$ and the effective pressures in radial, $`p`$, and azimuthal, $`p_{}`$, directions are functions of the radial coordinate only. Using the Wightman function from (2) and the VEV for the field square from (32), the components of the vacuum energy-momentum tensor can be presented in the form

$$q(a,b,r)=q(a,r)+q_a(a,b,r),a<r<b,q=\epsilon ,p,p_{},$$
(49)
where $`q(a,r)`$ are the corresponding functions for the vacuum outside a single sphere with radius $`a`$. In (49) the additional components are in the form

$`q_a(a,b,r)`$ $`=`$ $`{\displaystyle \frac{\sigma ^{1D}}{2\pi r^nS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z^3\mathrm{\Omega }_{a\nu _l}(az,bz)}{\sqrt{z^2m^2}}}F_{\nu _l}^{(q)}[G_{\nu _l}^{(a)}(az,y),G_{\nu _l}^{(a)}(az,y)]_{y=zr},`$ (50)
$`q`$ $`=`$ $`\epsilon ,p,p_{},`$
where for arbitrary functions $`f(y)`$ and $`g(y)`$ the functions $`F_\nu ^{(q)}[f(y),g(y)]`$ are defined by the relations

$`F_{\nu _l}^{(\epsilon )}[f(y),g(y)]`$ $`=`$ $`(14\xi )\left[f^{}g^{}{\displaystyle \frac{n}{2y}}(fg)^{}+\left({\displaystyle \frac{\nu _l^2}{y^2}}{\displaystyle \frac{1+4\xi 2(mr/y)^2}{14\xi }}\right)fg\right]`$ (51)
$`F_{\nu _l}^{(p)}[f(y),g(y)]`$ $`=`$ $`f^{}g^{}+{\displaystyle \frac{\stackrel{~}{\xi }}{2y}}(fg)^{}\left(1+{\displaystyle \frac{\nu _l^2+\stackrel{~}{\xi }n/2}{y^2}}\right)fg,\stackrel{~}{\xi }=4(n+1)\xi n`$ (52)
$`F_{\nu _l}^{(p_{})}[f(y),g(y)]`$ $`=`$ $`(4\xi 1)f^{}g^{}{\displaystyle \frac{\stackrel{~}{\xi }}{2y}}(fg)^{}+\left[4\xi 1+{\displaystyle \frac{\nu _l^2(1+\stackrel{~}{\xi })+\stackrel{~}{\xi }n/2}{(n+1)y^2}}\right]fg,`$ (53)
and the prime denotes the differentiation with respect to $`y`$. Quantities (50) are finite for $`ar<b`$ and diverge at the surface $`r=b`$. Similar to the Wightman function, the components of the vacuum energy-momentum tensor for a single sphere case can be presented in the form

$$q(a,r)=q_\mathrm{m}(r)+q_b(a,r),$$
(54)
where $`q_\mathrm{m}(r)`$ are the corresponding quantities for the boundary-free monopole geometry and the expressions for the sphere induced parts $`q_b(a,r)`$ are given in .

On the basis of formula (2), the vacuum energy-momentum tensor components may be written in another equivalent form:

$$q(a,b,r)=q(b,r)+q_b(a,b,r),a<r<b,q=\epsilon ,p,p_{},$$
(55)
with $`q(b,r)`$ being the corresponding components for the vacuum inside a single sphere with radius $`b`$. Here the additional components are given by the formula

$`q_b(a,b,r)`$ $`=`$ $`{\displaystyle \frac{\sigma ^{1D}}{2\pi r^nS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z^3\mathrm{\Omega }_{b\nu _l}(az,bz)}{\sqrt{z^2m^2}}}F_{\nu _l}^{(q)}[G_{\nu _l}^{(b)}(bz,y),G_{\nu _l}^{(b)}(bz,y)]_{y=zr},`$ (56)
$`q`$ $`=`$ $`\epsilon ,p,p_{}.`$ (57)
This expressions are finite for all $`a<rb`$ and diverge at the inner sphere surface $`r=a`$.

It follows from the above that if we present the vacuum energy-momentum tensor components in the form

$$q(a,b,r)=q_\mathrm{m}(r)+q_b(a,r)+q_b(b,r)+\mathrm{\Delta }q(a,b,r),a<r<b,$$
(58)
then the quantities

$$\mathrm{\Delta }q(a,b,r)=q_a(a,b,r)q(b,r)=q_b(a,b,r)q(a,r)$$
(59)
are finite for all $`arb`$. Near the surface $`r=a`$ it is suitable to use the first equality in (59), as for $`ra`$ both summands are finite. For the same reason the second equality is suitable for calculations near the outer surface $`r=b`$. Using formula (30) for the corresponding part of the Wightman function, it can be seen that the following formula takes place for the ”interference” parts

$`\mathrm{\Delta }q(a,b,r)`$ $`=`$ $`{\displaystyle \frac{\sigma ^{1D}}{2\pi r^nS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{z^3dz}{\sqrt{z^2m^2}}}{\displaystyle \frac{\overline{I}_{\nu _l}^{(a)}(az)\overline{K}_{\nu _l}^{(b)}(bz)}{\overline{I}_{\nu _l}^{(a)}(az)\overline{K}_{\nu _l}^{(b)}(bz)\overline{I}_{\nu _l}^{(b)}(bz)\overline{K}_{\nu _l}^{(a)}(az)}}`$ (60)
$`\times \left\{{\displaystyle \frac{F_{\nu _l}^{(q)}[I_\nu (y),G_{\nu _l}^{(b)}(bz,y)]}{\overline{I}_{\nu _l}^{(b)}(bz)}}{\displaystyle \frac{F_{\nu _l}^{(q)}[K_\nu (y),G_{\nu _l}^{(a)}(az,y)]}{\overline{K}_{\nu _l}^{(a)}(az)}}\right\}_{y=zr}.`$
It can be checked that this quantities satisfy the covariant continuity equation

$$r\frac{d\mathrm{\Delta }p}{dr}+(D1)(\mathrm{\Delta }p\mathrm{\Delta }p_{})=0.$$
(61)
Note that the ambiguities of the renormalization procedure for the VEV of the energy-momentum tensor in the form of an arbitrary mass scale (see ) are contained in the boundary-free parts $`q_\mathrm{m}(r)`$ of the corresponding components. The boundary induced parts $`q_b(\alpha ,r)`$, $`\alpha =a,b`$ and $`\mathrm{\Delta }q(a,b,r)`$ are unambiguously defined for $`a<r<b`$. In particular, for the massless conformally coupled scalar they contain no conformal anomalies and are traceless.

Now we turn to the investigation of limiting cases of formula (60). In the limit $`a0`$, the subintegrand behaves as $`a^{2\nu _l}`$ and the dominant contribution comes from the $`l=0`$ term:

$`\mathrm{\Delta }q(a,b,r)`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}(a/2b)^{2\nu _0}}{2\pi S_Dr^nb^3\nu _0\mathrm{\Gamma }^2(\nu _0)}}{\displaystyle \frac{n+2\eta _a\nu _0}{n+2\nu _0}}{\displaystyle _{mb}^{\mathrm{}}}{\displaystyle \frac{z^{2\nu _0+3}dz}{\sqrt{z^2m^2b^2}}}{\displaystyle \frac{\overline{K}_{\nu _0}^{(b)}(z)}{\overline{I}_{\nu _0}^{(b)}(z)}}`$ (62)
$`\times \{{\displaystyle \frac{\overline{K}_{\nu _0}^{(b)}(z)}{\overline{I}_{\nu _0}^{(b)}(z)}}F_{\nu _0}^q[I_{\nu _0}(y),I_{\nu _0}(y)]2F_{\nu _0}^q[I_{\nu _0}(y),K_{\nu _0}(y)]\}_{y=zr/b},a0,`$
assuming that $`\nu _0>0`$. For $`\nu _0=0`$ and $`a0`$ the ”interference” parts (60) behave as $`1/\mathrm{ln}a`$. Consider the limit $`b\mathrm{}`$ for fixed $`a`$ and $`r`$ in the case of a massless scalar field. By changing the integration variable to $`y=bz`$ and using the formula for the modified Bessel functions in the case of small values of the argument, we see that to the leading order the subintegrand in (60) behaves as $`(a/b)^{2\nu _l+1}`$ (assuming that $`\nu _l>0`$). The leading contribution is due to the $`l=0`$ term and one has:

$`\mathrm{\Delta }q(a,b,r)`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}(a/2b)^{2\nu _0+1}}{\pi S_Dar^D\mathrm{\Gamma }^2(\nu _0+1)}}{\displaystyle \frac{A_a+B_a\nu _0}{A_aB_a\nu _0}}\left[(2\nu _0+n)f_{1\nu _0}^{(q)}{\displaystyle \frac{A_a+B_a\nu _0}{A_aB_a\nu _0}}\left({\displaystyle \frac{a}{r}}\right)^{2\nu _0}+2f_{2\nu _0}^{(q)}\right]`$ (63)
$`\times {\displaystyle _0^{\mathrm{}}}dzz^{2\nu _0}{\displaystyle \frac{nK_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}zK_{\nu _0}^{}(z)}{nI_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}zI_{\nu _0}^{}(z)}},b\mathrm{},`$
with notations

$`f_{1\nu }^{(\epsilon )}=\nu (14\xi ),f_{1\nu }^{(p)}={\displaystyle \frac{\stackrel{~}{\xi }}{2}},f_{1\nu }^{(p_{})}={\displaystyle \frac{\nu +1/2}{n+1}}\stackrel{~}{\xi },`$ (64)
$`f_{2\nu }^{(\epsilon )}=0,f_{2\nu }^{(p)}=n1,f_{2\nu }^{(p_{})}=4n(n+1)\xi /\sigma ^2.`$ (65)
As we see in both limits $`a0`$ and $`b\mathrm{}`$ the ”interference” parts for the energy-momentum tensor components vanish.

In the limit $`a,b\mathrm{}`$ for fixed values $`ba`$ and $`\sigma `$, by the calculations similar to those for the ”interference” part of the field square, one receives

$$\mathrm{\Delta }q(a,b,r)\frac{(ba)^{D1}}{(4\pi )^{D/2}\mathrm{\Gamma }(D/2+1)}_{m_0}^{\mathrm{}}𝑑y\frac{(y^2m_0^2)^{D/2}}{c_b(y)e^{2y}/c_a(y)1}\mathrm{\Delta }F^q(y,r),$$
(66)
where

$$\mathrm{\Delta }F^\epsilon (y,r)=\mathrm{\Delta }F^p_{}(y,r)=1\frac{4D(\xi \xi _D)y^2m_0^2}{2(y^2m_0^2)}[F(y,r)+2],\mathrm{\Delta }F^p=1,$$
(67)
with the function $`F(y,r)`$ defined by (41). These expressions are exactly the same as the corresponding expressions for the geometry of two parallel plates on the Minkowski background investigated in  (note that in , the notations $`\stackrel{~}{B}_a/\stackrel{~}{A}_a=\beta _1`$ and $`\stackrel{~}{B}_b/\stackrel{~}{A}_b=\beta _2`$ are used) for a massless scalar and in Ref.  for the massive case.

And finally, let us consider the limit of the strong gravitational field, corresponding to $`\sigma 1`$. In this limit for $`\xi >0`$ one has $`\nu _l1`$. After the change of variable to $`y=bz/\nu _l`$, we replace in (60) the Bessel modified functions by their uniform asymptotic expansions for large values of the order. The leading contribution comes from the $`l=0`$ summand and we have the following limit for the ”interference” parts of the vacuum energy-momentum components:

$$\mathrm{\Delta }p\frac{\mathrm{\Delta }p_{}}{D1}\frac{\sigma ^{1D}\eta _a\eta _b\stackrel{~}{\nu }^{3/2}}{r^DS_D\sqrt{2\pi (b^2a^2)}}e^{2\stackrel{~}{\nu }\mathrm{ln}(b/a)},\mathrm{\Delta }\epsilon /\mathrm{\Delta }p\sigma ,$$
(68)
where $`\stackrel{~}{\nu }`$ is defined in (44). In this limit the ”interference” parts are exponentially suppressed with respect to single spheres contributions.

Now we turn to the interaction forces between the spheres. The vacuum force acting per unit surface of the sphere $`r=\alpha `$, $`\alpha =a,b`$, is determined by the $`{}_{1}{}^{}{}_{}{}^{1}`$–component of the vacuum energy-momentum tensor at this point. By virtue of relations (49) and (55), the corresponding effective pressures can be presented as a sum of the pressure for a single sphere with $`r=\alpha `$ when the second sphere is absent and the pressure induced by the presence of the second sphere, $`p_\alpha (a,b,r=\alpha )`$. The first term is divergent due to the well-known surface divergences and needs additional regularization. The second term is finite and can be termed as an interaction force between the spheres. This additional radial vacuum pressure on the sphere with $`r=\alpha `$, $`\alpha =a,b`$ due to the existence of the second sphere can be found from (50) and (56), respectively. Using the relations

$$G_\nu ^{(r)}(rz,rz)=B_r,rz\frac{}{y}G_\nu ^{(r)}(rz,y)_{y=rz}=A_r,r=a,b,$$
(69)
they can be presented in the form

$`p_\alpha (a,b,r=\alpha )`$ $`=`$ $`{\displaystyle \frac{\sigma ^{1D}}{2\pi \alpha ^DS_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{zdz}{\sqrt{z^2m^2}}}\mathrm{\Omega }_{\alpha \nu _l}(az,bz)`$
$`\times `$ $`\left\{\stackrel{~}{A}_\alpha ^24(D1)\xi \stackrel{~}{A}_\alpha B_\alpha \left[z^2\alpha ^2+\nu _l^2n^2/4\right]B_\alpha ^2\right\},\alpha =a,b.`$
Unlike the self-action forces, these quantities are finite for $`a<b`$ and need no further regularization. For the Dirichlet scalar one has $`\mathrm{\Omega }_{\alpha \nu _l}(az,bz)>0`$ and, hence, $`p_\alpha (a,b,r=\alpha )<0`$. This means that in this case the vacuum interaction forces between the spheres are attractive. Using the Wronskian relation for the Bessel modified functions, it can be seen that

$$\left[A_\alpha B_\alpha (z^2\alpha ^2+\nu ^2)\right]\mathrm{\Omega }_{\alpha \nu }(az,bz)=n_\alpha \alpha \frac{}{\alpha }\mathrm{ln}\left|1\frac{\overline{I}_\nu ^{(a)}(az)}{\overline{I}_\nu ^{(b)}(bz)}\frac{\overline{K}_\nu ^{(b)}(bz)}{\overline{K}_\nu ^{(a)}(az)}\right|,$$
(71)
where $`n_a=1,n_b=1`$. This allows us to write the expressions (4) for the interaction forces per unit surface in another equivalent form:

$`p_\alpha (a,b,r=\alpha )`$ $`=`$ $`{\displaystyle \frac{n_\alpha \sigma ^{1D}}{2\pi \alpha ^{D1}S_D}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}D_l{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{zdz}{\sqrt{z^2m^2}}}\left[1{\displaystyle \frac{\stackrel{~}{\xi }\stackrel{~}{A}_\alpha B_\alpha }{A_\alpha ^2B_\alpha ^2(z^2\alpha ^2+\nu _l^2)}}\right]`$
$`\times {\displaystyle \frac{}{\alpha }}\mathrm{ln}|1{\displaystyle \frac{\overline{I}_{\nu _l}^{(a)}(az)}{\overline{I}_{\nu _l}^{(b)}(bz)}}{\displaystyle \frac{\overline{K}_{\nu _l}^{(b)}(bz)}{\overline{K}_{\nu _l}^{(a)}(az)}}|,\alpha =a,b,`$
where $`\stackrel{~}{\xi }`$ is defined in (52). For Dirichlet and Neumann scalars the second term in the square brackets is zero.

Let us consider the interaction forces between the spheres in limiting cases. For small values of the radius of the inner sphere, $`a0`$, the leading contribution comes from the $`l=0`$ summand. The corresponding terms are given by formulae

$`p_a(a,b,r=a)`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}a^{2\nu _0D}}{2^{2\nu _01}\mathrm{\Gamma }^2(\nu _0)\pi S_D}}{\displaystyle \frac{n+2\eta _a\nu _0}{n+2\nu _0}}{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z^{2\nu _0+1}}{\sqrt{z^2m^2}}}{\displaystyle \frac{\overline{K}_{\nu _0}^{(b)}(bz)}{\overline{I}_{\nu _0}^{(b)}(bz)}},`$ (73)
$`p_b(a,b,r=b)`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}a^{2\nu _0}}{2^{2\nu _0}\nu _0\mathrm{\Gamma }^2(\nu _0)\pi S_Db^D}}{\displaystyle \frac{n+2\eta _a\nu _0}{n+2\nu _0}}{\displaystyle _m^{\mathrm{}}}𝑑z{\displaystyle \frac{z^{2\nu _0+1}}{\sqrt{z^2m^2}\overline{I}_{\nu _0}^{(b)2}(bz)}}`$ (74)
$`\times \left\{\stackrel{~}{A}_b^24(D1)\xi \stackrel{~}{A}_bB_b\left[z^2b^2+(1/\sigma ^21)n(n+1)\xi \right]B_b^2\right\},`$
where $`\eta _a`$ is defined in (37).

For large values of the radius of the outer sphere, $`b\mathrm{}`$, in the case of a massless scalar to the leading order we find

$`p_a(a,b,r=a)`$ $``$ $`{\displaystyle \frac{2\sigma ^{1D}(A_aB_a\nu _0)^2}{\pi a^DbS_D\mathrm{\Gamma }^2(\nu _0)}}\left({\displaystyle \frac{a}{2b}}\right)^{2\nu _0}\left[\stackrel{~}{A}_a^24(D1)\xi \stackrel{~}{A}_aB_a(\nu _0^2n^2/4)B_a^2\right]`$ (75)
$`\times {\displaystyle _0^{\mathrm{}}}dzz^{2\nu _0}{\displaystyle \frac{nK_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}K_{\nu _0}^{}(z)}{nI_{\nu _0}(z)2\delta _{0\stackrel{~}{A}_b}I_{\nu _0}^{}(z)}},`$
$`p_b(a,b,r=b)`$ $``$ $`{\displaystyle \frac{\sigma ^{1D}(a/2b)^{2\nu _0}}{\pi \nu _0\mathrm{\Gamma }^2(\nu _0)b^{D+1}S_D}}{\displaystyle \frac{A_a+B_a\nu _0}{A_aB_a\nu _0}}{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle \frac{z^{2\nu _0}}{I_{\nu _0}^2(z)}}.`$ (76)
In (76) we have assumed that $`\stackrel{~}{A}_b0`$. In the case $`\stackrel{~}{A}_b=0`$ (Neumann boundary condition on the outer sphere) the subintegrand is equal to $`(z^2+\nu _0^2n^2/4)z^{2\nu _0}[nI_{\nu _0}(z)/2zI_{\nu _0}^{}(z)]^2`$.

Now let us consider the limit $`\sigma 1`$. For $`\xi >0`$ one has $`\nu _l1`$ for all $`l`$, and to the leading order we have

$$p_\alpha (a,b,r=\alpha )\frac{\eta _a\eta _b\sigma ^{1D}\stackrel{~}{\nu }^{3/2}}{\sqrt{2\pi }\alpha ^DS_D\sqrt{b^2a^2}}e^{2\stackrel{~}{\nu }\mathrm{ln}(b/a)},\alpha =a,b.$$
(77)
In the case of a minimally coupled scalar, $`\xi =0`$, the main contribution is due to the summand $`l=0`$ with $`\nu _l=n/2`$ and is of an order $`\sigma ^{1D}`$. In this case the contributions coming from the $`l>0`$ summands are exponentially suppressed.

As in the cases of the VEVs for the field square and energy-momentum tensor, the dependence of the interaction forces on the parameter $`\sigma `$ is essentially different for minimally and non-minimally coupled scalars. This feature is illustrated in figures 2 and 3 where we have plotted the dependence of the interaction forces $`p_a(a,b,r=a)`$ (left panels) and $`p_b(a,b,r=b)`$ (right panels) between the spheres on $`a/b`$ for conformally and minimally coupled $`D=3`$ Dirichlet scalars on background of the global monopole background with $`\sigma =1`$ (a), $`\sigma =0.4`$ (b), $`\sigma =0.2`$ (c).

## 5 Conclusion

In this paper, we present the quantum vacuum effects produced by two concentric spherical shells in the $`D+1`$-dimensional point-like global monopole spacetime defined by the line element (1). The case of a massive scalar field with general curvature coupling parameter and satisfying the Robin boundary conditions on the spheres is considered. To derive formulae for the vacuum expectation values of the square of the field operator and the energy-momentum tensor, we first construct the positive frequency Wightmann function. This function is also important in considerations of the response of a particle detector at a given state of motion through the vacuum under consideration . The application of the generalized Abel-Plana formula to the mode sum over zeros of the combinations of the cylinder functions allows us to extract the part due to a single sphere on background of the global monopole geometry, formula (29). In this formula the ”interference” part $`\mathrm{\Delta }W(x,x^{})`$, given by formula (30), is finite in the coincidence limit for all values $`arb`$. The VEV of the square of the scalar field, $`0|\phi ^2(x)|0`$, is given by the evaluation of the Wightman function in the coincidence limit. The expectation values for the energy-momentum tensor are obtained by applying on the corresponding Wightman function a certain second-order differential operator and taking the coincidence limit. In both cases the expectation values can be presented as a sum of boundary-free global monopole, single sphere induced and ”interference” terms. The VEVs related to the spacetime of a point-like global monopole without boundaries are considered in , and the effects produced by a single sphere are investigated in a previous paper . Note that for the points not lying on the spheres, the boundary induced terms are unambiguously defined and the ambiguities in the renormalization procedure in the form of an arbitrary mass scale are contained in the boundary-free parts only. In particular, for the massless conformally coupled scalar the boundary induced vacuum energy-momentum tensor contains no conformal anomalies and is traceless. In this paper we concentrate on the ”interference” parts in the local Casimir densities and interaction forces between the spheres. The application of the generalized Abel-Plana formula allows us to derive closed expressions for these quantities, given by formula (35) in the case of the field square, and by formula (60) for the energy-momentum tensor. As bulk divergences are contained in the part corresponding to the global monopole geometry without boundaries, and the surface divergences are contained in the single sphere parts, the ”interference” parts are finite for all values $`arb`$. In particular, the integrals in the corresponding formulae are exponentially convergent and they are useful for numerical evaluations. We have considered various limiting cases of the formulae for the ”interference” parts. In the limits $`a0`$ or $`b\mathrm{}`$ for a fixed value of $`r`$, these parts vanish as $`a^{2\nu _0}`$ and $`b^{(2\nu _0+1)}`$ respectively, where $`\nu _0`$ is defined by relation (10) with $`l=0`$. In the limit $`a,b\mathrm{}`$ for a fixed value $`ba`$, the leading terms in the boundary produced parts do not depend on the parameter $`\sigma `$ associated with the solid angle deficit, and are exactly the same as the corresponding quantities for the geometry of two parallel Robin plates on the Minkowski background. We have also investigated the limit of strong gravitational fields corresponding to small values of the parameter $`\sigma `$, $`\sigma 1`$. In this limit the behaviour of the Casimir densities is drastically different for minimally ($`\xi =0`$) and non-minimally ($`\xi 0`$) coupled scalars. In the minimal coupling case the leading terms of the corresponding asymptotic expansions for both field square and the energy-momentum tensor VEVs behave as $`\sigma ^{1D}`$. For a non-minimally coupled scalar, the ”interference” parts behave as $`\phi ^2^{(ab)}\sigma ^{3/2D}\mathrm{exp}(\gamma /\sigma )`$ and $`\mathrm{\Delta }p\mathrm{\Delta }p_{}\mathrm{\Delta }\epsilon /\sigma \sigma ^{D1/2}\mathrm{exp}(\gamma /\sigma )`$, with $`\gamma =2\sqrt{n(n+1)\xi }\mathrm{ln}(b/a)`$. Note that in this case the ”interference” parts are exponentially suppressed with respect to single sphere contributions. The vacuum forces acting on spheres contain two terms. The first ones are the forces acting on a single sphere then the second boundary is absent. Due to the well–known surface divergences in the VEV’s of the energy-momentum tensor these forces are infinite and need an additional regularization. The another terms in the vacuum forces are finite and are induced by the presence of the second boundary. They correspond to the interaction forces per unit surface between the spheres and are determined by formula (4). For the Dirichlet scalar these forces are always attractive. In the limit of the strong gravitational field, $`\sigma 1`$, for the minimally coupled scalar field the interactions forces behave as $`\sigma ^{1D}`$, whereas for a non-minimally coupled scalar they are exponentially small, $`p_\alpha (a,b,r=\alpha )\sigma ^{D1/2}\mathrm{exp}(\gamma /\sigma )`$.

## Acknowledgement

We acknowledge support from the Research Project of the Kurdistan University. The work of AAS was supported in part by the Armenian Ministry of Education and Science (Grant No. 0887).


format: text

text:
# An ab-initio many-body method for electronic structure calculations of solids. II. Unscreened Hartree-Fock treatment for the 3d systems Fe, Co, Ni and Cu

## I Introduction

In the preceding paperSchnell2003-1 (1) we suggested a new method for ab-initio electronic-structure calculations of solids. The main steps of this procedure are:

1. Perform a conventional, self-consistent, band-structure calculation for an effective one-particle Hamiltonian, namely, the Hartree Hamiltonian, to obtain a suitable basis set of Bloch functions.
2. By taking into account only a finite number $`J`$ of bands one chooses a truncated one-particle Hilbert space. The Marzari-VanderbiltMV97 (2) algorithm is then used to construct a maximally localized set of Wannier functions, which span the same truncated one-particle Hilbert space.
3. All one-particle (tight-binding) and two-particle (Coulomb) matrix elements of the Hamiltonian within this Wannier function basis are calculated.
4. The resulting electronic many-body Hamiltonian in second quantization with parameters determined from first principles is studied within standard many-body approximations for lattice electron systems.

In this paper we apply this scheme to the 3d transition metals Fe, Co, Ni, and Cu. We use the “linear muffin-tin orbital” (LMTO) method within the “atomic-sphere approximation” (ASA)Skriver (3) to perform the band-structure calculation for the Hartree Hamiltonian in first quantization. The direct Coulomb matrix elements of the maximally localized Wannier basis are rather large, about 20 eV in magnitude. We then use the simplest possible many body approximation, the Hartree-Fock approximation (HFA), to study the second-quantized multiple-band Hamiltonian. Our results are compared with those obtained from a standard LSDA calculation.DrGross90 (4, 5) Although the 3d-bands and the 4s-band overlap in the L(S)DA approximation, our unscreened HFA calculations give 3d-bands that lie considerably lower (between 10 and 20 eV) than the 4s-band. The HFA correctly predicts ferromagnetism for the ferromagnetic metals Fe, Co, and Ni and no magnetism for Cu, but with a much larger exchange splitting between majority and minority 3d bands than obtained within LSDA and with a slightly larger magnetic moment per site than obtained experimentally or within LSDA. On the other hand, the total energy is lower in HFA than in LSDA. The LSDA results for metals are probably more reliable than our new HFA results, which lack important screening and correlation effects. In order for our method to go beyond LSDA we would need to use better many body methods than the (unscreened) HFA, which should be possible within our scheme.

Our purpose in presenting the HFA results is to demonstrate that the new ab-initio many-body method that we have proposed is feasible and can be applied to practical calculations of materials. In addition, because HFA is the assumed standard starting point before adding complex many-body correlations, and because it is the simplest many-body approximation (or the best mean-field one-particle approximation), it is useful to know what the HFA predicts for the 3d transition metals, which have been so heavily studied by other techniques. Comparisons of HFA results with experiments, other (higher order) theories, or established standard methods such as the L(S)DA should demonstrate the effects of correlation on electronic properties in d-electron systems.

To the best of our knowledge we do not know of any published HFA results (band structure, density of states, magnetism, magnetic moment, total energy, etc.) for the 3d ferromagnets Fe, Co and Ni, unless it was implicitly applied to these materials for schemes like the local ansatzStollhoff (6), where HFA results serve as an input to higher order calculations . The is not surprising since the HFA has, from very early on, been viewed as a poor approximation for metals. For example, when applied to the homogeneous electron gas (as the simplest model of an infinite metallic system), the HFA has well-known Fermi edge singularitiesMahan (7, 8). These lead, in particular, to a vanishing density of states (DOS) at the Fermi energy, which is, of course, unphysical. This unphysical feature usually prevails in actual HFA-calculations for real metalsMonkhorst79 (9), though sometimes this singularity is hard to see in actual HFA-resultsDovesiPisani82 (10). In our calculations the non-locality is handled through the calculation of expectation values (matrix elements of the density matrix), which makes HF calculations as easy as Hartree calculations. Furthermore, because of our localized Wannier basis, we only keep on-site and a few inter-site Coulomb and exchange matrix elements. Hence our calculations have an effective short-ranged Coulomb interaction. Although longer-range Coulomb matrix elements are small in our calculations, which is why we truncate them, it is possible that if all of them were kept to infinite distances that they could add up to give Fermi edge singularities (which are due to the long-ranged nature of the bare Coulomb interaction) and other standard anomalies. Correlation or screening would quickly kill these effects.

The paper is organized as follows. In Section II we describe the LMTO-Hartree calculations and the Hartree-results for the band structure and density of states, respectively. Section III describes some results obtained for the maximally localized Wannier functions obtained within the Marzari-Vanderbilt algorithmMV97 (2), in particular their localization properties. Results for the matrix elements, in particular the direct Coulomb and exchange matrix elements are given in Section IV; we also compare these results with calculations of the Slater integrals. The application of the (unscreened) HFA to the multiband many-body Hamiltonian in second quantization is the subject of Section V. For an interpretation of the results we compare the numerical HFA results obtained for the crystal with previous atomic HFA results and with numerical and analytical results for a simplified local (atomic or zero band width) model in Section VI. A comparison with the more standard LSDA-results follows in Section VII, before the paper closes with a short discussion.

## II LMTO Hartree calculation

For the four materials of interest (Fe, Co, Ni, Cu) we performed a selfconsistent Hartree band-structure calculation. Besides the nuclear charge we used the (experimentally) known results for the lattice structure (bcc for Fe, fcc otherwise; Co should actually be hexagonal) and for the lattice constant as input. For the band-structure calculation we used the LMTO methodSkriver (3, 11) within the atomic sphere approximation (ASA). The radius of the (overlapping) muffin-tin spheres (the Wigner-Seitz radius $`S`$) is determined by the condition that the sphere volume equals the volume of the unit cell. This yields the following values for the Wigner-Seitz radius: $`S=2.662a_0`$ for Fe, $`S=2.621a_0`$ for Co, $`S=2.602a_0`$ for Ni and $`S=2.669a_0`$ for Cu (Ref. Skriver, 3). Within the muffin-tin spheres the potential and wave functions are expanded in spherical harmonics with a cutoff $`l_{\mathrm{max}}=2`$, i.e., s, p, and d-orbitals are included.

Within the LMTO-ASA the eigenfunctions, i.e. the Bloch wave functions (inside a muffin-tin sphere), are given in terms of the solution to the radial Schrödinger equation $`\varphi _{\nu l}(r)`$ to some fixed energies $`E_{\nu l}`$ and its energy derivative $`\dot{\varphi }_{\nu l}(r)`$:

$$\mathrm{\Psi }_{n𝐤}(𝐫)=\underset{L}{}\left(\varphi _{\nu l}(r)A_L^{n𝐤}+\dot{\varphi }_{\nu l}(r)B_L^{n𝐤}\right)Y_L(\widehat{𝐫}),$$
(1)
where the $`Y_LY_l^m`$ denote the (complex) spherical harmonics and $`n`$ is the band index.

Figs. 1 \- 4 show the nonmagnetic (spin-degenerate) Hartree results for the density of states (DOS). One observes that the five narrow 3d-bands are similar in energy and therefore hybridize with the 4sp-bands. For Fe, Co, and Ni the Fermi level $`E_F`$ falls in the midst of the partially filled 3d-bands, and hence there is a large DOS at the Fermi level. For Cu $`E_F`$ is above the filled 3d bands, where the states have more 4sp-like character and the DOS is small (free electron like).

The Bloch functions obtained from this Hartree calculation are used as the basis of a truncated Hilbert space. In the band calculation, there is already a natural truncation due to the spherical harmonics expansion cutoff $`l_{max}=2`$, which limits the calculation to 9 bands for one atom per unit cell. The connection between the Bloch functions in the different muffin-tin spheres is given by the standard Bloch relation

$$\psi _{n𝐤}(𝐫+𝐑)=e^{i\mathrm{𝐤𝐑}}\psi _{n𝐤}(𝐫).$$
(2)
Hence the Bloch function in a single muffin-tin sphere determines the function for the whole crystal. However, this simple relation holds for Bloch functions only. Any other function (to be represented as linear combination of Bloch functions) can, nevertheless, be decomposed into its contributions within the different muffin-tins according to

$$\mathrm{\Phi }_\alpha (𝐫)=\underset{i}{}\mathrm{\Phi }_\alpha (𝐑_i;𝐫𝐑_i)$$
(3)
where $`\mathrm{\Phi }_\alpha (𝐑;𝐫)=0`$ for $`|𝐫|>S`$.

## III Maximally localized Wannier functions

We now turn to the construction and description of the Wannier functions and will present results for their main properties for the four materials under investigation. Because

$$w_{𝐑n}(𝐫)=w_{\mathrm{𝟎}n}(𝐫𝐑),$$
(4)
it is sufficient to investigate the Wannier functions only for one arbitrarily chosen lattice vector $`\mathrm{𝟎}`$: $`w_n(𝐫)w_{\mathrm{𝟎}n}(𝐫)`$. This Wannier function has contributions not only in the muffin-tin sphere around $`\mathrm{𝟎}`$ but also in other muffin-tin spheres, and we can use the decomposition (3) into the contributions from the different muffin-tin spheres

$$w_n(𝐫)=\underset{𝐑}{}w_n(𝐑;𝐫𝐑).$$
(5)
After application of the Marzari-Vanderbilt algorithm, which is briefly summarized in the preceding paperSchnell2003-1 (1) and in more detail in Ref. MV97, 2, the new set of bands that are used to calculate the maximally localized Wannier functions can no longer be classified by pure angular momentum quantum numbers. The Wannier functions are rather admixtures having different $`L`$-contributions (3d, 4s, 4p etc.). But, since the original Bloch functions from which the Wannier functions are constructed were given in terms of a spherical harmonics expansion, the new Wannier functions (and their contribution in each individual muffin-tin sphere) can also be decomposed into these spherical harmonics contributions

$$w_n(𝐑;𝐫)=\underset{L}{}\left\{\varphi _{\nu l}(r)A_L^{𝐑n}+\dot{\varphi }_{\nu l}(r)B_L^{𝐑n}\right\}Y_L(\widehat{𝐫}).$$
(6)
One can then calculate the weight of the contributions to the Wannier function (centered at $`\mathrm{𝟎}`$) within the different muffin-tin spheres

$$w_n|w_n_𝐑_𝐑d^3𝐫|w_n(𝐫)|^2=_\mathrm{𝟎}d^3𝐫|w_n(𝐑;𝐫)|^2,$$
(7)
and one can also decompose this into the different $`l`$-contributions according to:

$$w_n|w_n_𝐑=\underset{l}{}\underset{C_l^{𝐑n}}{\underset{}{\underset{m=l}{\overset{l}{}}\left\{|A_{lm}^{𝐑n}|^2+\dot{\varphi }_{\nu l}^2|B_{lm}^{𝐑n}|^2\right\}}}$$
(8)
For the 3d-system iron these quantities are tabulated in Table 1. The first line is the weight $`w_n|w_n_\mathrm{𝟎}`$ in the center muffin-tin. Between 88 and 98% of the total weight of the Wannier functions is to be found already within the center muffin-tin; this shows how well localized our Wannier functions are with the lowest five functions having values of more than 95%. Rows 2–4 in this table indicate the different $`l`$-contribution or $`l`$-character of the Wannier functions. One sees that the optimally localized Wannier functions are not pure within their $`l`$-character, but the lowest five Wannier functions (0-4) still have mainly $`l=2`$ (3d) character. Higher band-index states (which are slightly less well localized according to row 1) are admixtures that have mainly $`l=1`$ (4p) character (about 50 %), but also a considerable amount of $`l=0`$ (4s) and $`l=2`$ (3d) character. Corresponding results for the other 3d-systems Co, Ni, and Cu are similar and, therefore, not repeated here.

## IV One particle and Coulomb matrix elements

From the optimally localized Wannier functions we calculate the one-particle matrix elements

$$t_{12}=d^3𝐫w_1^{}(𝐫)\left(\frac{\mathrm{}^2}{2m}^2+V(𝐫)\right)w_2(𝐫)$$
(9)
and the Coulomb matrix elements of the Hamiltonian

$$W_{12,34}=d^3𝐫d^3𝐫^{}w_1^{}(𝐫)w_2^{}(𝐫^{})\frac{e^2}{|𝐫𝐫^{}|}w_3(𝐫^{})w_4(𝐫).$$
(10)
Here we use the abbreviated notation 1 to mean $`𝐑_1n_1`$ and 2 to mean for $`𝐑_2n_2`$, etc. We used two different numerical algorithms to calculate these Coulomb matrix elements, namely the FFT-algorithm briefly described in the previous paper and a spherical expansion algorithm described in some detail in Ref. SCA02, 12. The latter method makes use of the fact that (in each muffin-tin sphere) the Wannier functions are explicitly given as linear combinations of products of spherical harmonics and a radial wave function. The expansion

$$\frac{1}{|𝐫𝐫^{}|}=\underset{k=0}{\overset{\mathrm{}}{}}\frac{4\pi }{2k+1}\frac{r_<^k}{r_>^{k+1}}\underset{m=k}{\overset{k}{}}Y_K^{}(\widehat{𝐫}^{})Y_K(\widehat{𝐫})$$
(11)
($`K=\{k,m\}`$) makes it possible to express the on-site Coulomb integrals as one-dimensional integrals over products of the radial functions and the (tabulated) Gaunt coefficients. The results obtained by this algorithm and by the independent FFT-algorithm agree within the numerical errors (of at most 1%).

Results for the on-site direct and exchange Coulomb matrix elements between the optimally localized Wannier functions are given in Table 2 for iron (Fe). The direct Coulomb integrals $`U_{nm}=W_{nm,mn}`$ between the Wannier states with the lowest five band indices ($`n,m\{0,\mathrm{},4\}`$), which according to the discussion and table in the preceding section have mainly 3d-character, are rather large, up to 23 eV for Fe. Within the 3d-like bands the interband direct Coulomb matrix elements are of the same magnitude as the intraband matrix elements. The matrix elements between 3d-states and 4sp-states are considerably smaller, of the magnitude of 13 - 14 eV. For electrons in 4sp-states ($`n,m\{5,\mathrm{},9\}`$) the direct intraband Coulomb matrix elements are again of the order of 13 - 14 eV, but the interband matrix elements are slightly smaller, about 9 eV. The exchange matrix elements $`J_{nm}=W_{nm,nm}`$ are always much smaller, usually less than 1 eV (for $`nm`$). Again the corresponding results for the other 3d-systems investigated (Co, Ni and Cu) are very similar.

For the 5 states with predominant 3d-character we have calculated the averages of the on-site direct and exchange Coulomb matrix elements

$`U`$ $``$ $`{\displaystyle \frac{1}{25}}{\displaystyle \underset{mm^{}}{}}W_{mm^{}m^{}m}`$ (12)
$`J`$ $``$ $`{\displaystyle \frac{1}{20}}{\displaystyle \underset{mm^{}}{}}W_{mm^{}mm^{}},`$ (13)
as well as the averages of the absolute values of the nearest neighbor (NN) and next nearest neighbor (NNN) hopping matrix elements

$$t_{NN(N)}\frac{1}{25}\underset{n,m}{}|t_{𝐑nm}|.$$
(14)
The results obtained thereby for the 4 transition metals under consideration are shown in Table 3. The U-values vary between 21 eV for Fe and 25 eV for Cu, the J-values are smaller than 1 eV and the hopping matrix elements are of the magnitude 0.5 – 0.7 eV for nearest-neighbor (NN) and 0.1 – 0.2 eV for next-nearest-neighbor (NNN), and further on decrease with increasing distance.

We have also evaluated the Slater integralsS29 (13):

$$F^ke^2𝑑rr^2𝑑r^{}r^2|R_{l=2}(r)|^2\frac{r_<^k}{r_>^{k+1}}|R_{l=2}(r^{})|^2,$$
(15)
where $`R_{l=2}(r)`$ is a radial (atomic) d-wave function (obtained by solving the Schrödinger equation for a radial symmetric potential, for instance). Note that only the three integrals $`F^0`$, $`F^2`$ and $`F^4`$ are required to determine all the Coulomb $`d`$-matrix elements. Using the radial d-wave function obtained from the Hartree calculation we obtain the following values for the Slater integrals of the four 3d-systems: $`F^0=`$ 21.62 eV for Fe, 23.18 eV for Co, 24.69 eV for Ni, and 26.27 eV for Cu. This means, the Slater integrals $`F^0`$ are rather good estimates of our (averaged) Coulomb matrix elements. These values are also in agreement with older results obtained in calculations for 3d-atomsWatson59 (14). In Table 4 we show our $`F^k`$-values for the four 3d-crystals and compare them with corresponding atomic calculations from Ref. Watson59, 14. Obviously, there is fairly good agreement between these atomic and our results.

## V Unscreened Hartree-Fock approximation

After we have determined the matrix elements within our restricted basis set of 9 maximally localized Wannier functions (per site and spin), we have a Hamiltonian in second quantization of the form

$$H=\underset{12\sigma }{}t_{12}c_{1\sigma }^{}c_{2\sigma }+\frac{1}{2}\underset{1234\sigma \sigma ^{}}{}W_{12,34}c_{1\sigma }^{}c_{2\sigma ^{}}^{}c_{3\sigma ^{}}c_{4\sigma }$$
(16)
for which all the matrix elements are known from first principles. The simplest approximation one can now apply is the HFA, which replaces the many-body Hamiltonian by the effective one-particle Hamiltonian

$`H_{HF}`$ $`=`$ $`{\displaystyle \underset{12\sigma }{}}\left(t_{12}+\mathrm{\Sigma }_{12,\sigma }^{HF}\right)c_{1\sigma }^{}c_{2\sigma }`$ (17)
$`\text{with }\mathrm{\Sigma }_{12,\sigma }^{HF}`$ $`=`$ $`\mathrm{\Sigma }_{12}^{Hart}+\mathrm{\Sigma }_{12,\sigma }^{Fock}`$
$`=`$ $`{\displaystyle \underset{34\sigma ^{}}{}}\left[W_{13,42}\delta _{\sigma \sigma ^{}}W_{31,42}\right]c_{3\sigma ^{}}^{}c_{4\sigma ^{}}.`$
Here the expectation values $`c_{1\sigma }^{}c_{2\sigma }`$ have to be determined selfconsistently for the HF Hamiltonian (17). Note that the Fock (exchange) term is spin ($`\sigma `$) dependent and may, therefore, give rise to magnetic solutions.

The Hartree-Fock results for the four materials of interest are shown in Figs. 512. We show the effective HF band structure and its density of states (DOS). In our HFA calculations there are no singularities (or a vanishing DOS) at the Fermi level since we start from a localized description and consider the Coulomb matrix elements only locally (on-site and for a few neighbor shells). Therefore, we implicitely truncate the Coulomb interaction in real space and in practice work with an effective short-ranged interaction. Within HFA the main part of the 3d-bands lies between 18 and 22 eV below the Fermi level and is separated from the 4sp-bands. We find magnetism in HFA for Fe, Co, and Ni in agreement with experiment. The five majority spin d-bands are about 20 eV below the Fermi energy and are completely filled. But the partially filled minority d-bands have two (for Fe), three (for Co), and four (for Ni) filled bands between -18 and -15 eV, and the rest are around and above the Fermi level. This results in magnetic moments per atom (in units of the Bohr magneton $`\mu _B`$/atom) of 2.9 for iron, 1.9 for cobalt, and 0.76 for nickel. For copper no magnetism and exchange splitting of the 3d-bands is obtained, but the (spin degenerate) 3d-bands are at about 22 eV below the Fermi level and separated from the 4sp-bands. If we compare these results with the results of the simple Hartree approximation shown in Figs. 14 we see that the exchange term has two effects: it produces an exchange splitting and the possibility of magnetic solutions, and it draws the 3d-bands energetically down by an amount of about 20 eV. Compared with experiment the HFA overestimates magnetism and leads to overly large values for the magnetic moment per site (about 10–30% too large; experiment gives magnetic moments of about 2.2 for iron, 1.7 for Co, and 0.6 for Ni). This is consistent with Heisenberg or Ising model studies where the mean-field approximation HFA also has the tendency to overestimate magnetism and magnetic solutions.

However, the reason why the 3d-bands lie so far below the Fermi level and the 4sp-band in HFA has nothing to do with the existence and overestimation of magnetism. This can be seen already from the non-magnetic system Cu, for which the (fully occupied) 3d-bands also lie at about 22 eV below the Fermi level (see Figs. 11 and 12). To demonstrate this also for a system with a partially filled 3d-band we have done a non-magnetic Hartree-Fock-calculation for Co (by forcing equal occupation for both spin directions). The results for the band structure and the DOS are shown in Figs. 13,14. We observe again that the main part of the 3d-bands are well below the 4s-bands and Fermi level; note the hybridization gap caused by the unoccupied 3d-bands above the Fermi level.

## VI  Comparison with atomic Hartree-Fock results

We have seen in the previous section that one effect of the HFA calculation, when compared with the Hartree calculation, is the shift of the 3d-bands down (about 20 eV below the Fermi level and about 8–10 eV below the bottom of the 4sp-band). This shift of the d-bands is about the same energy as the Coulomb matrix elements $`U`$, and roughly agrees earlier atomic Hartree-Fock calculationsWatson59 (14, 15), where the 3d-states are also about 10 eV below the 4s-states.

Because the inter-site hopping matrix elements in Table 3 are much smaller than the U-values one may consider an expansion in $`t/U`$, with the zeroth order approximation to completely neglect hopping. Doing this, we have performed a quasi-atomic HFA calculation for Co, by including only the on-site one-particle and two-particle (Coulomb) matrix elements. The results are summarized in Fig. 15; the degeneracy of the different levels is also indicated. In the paramagnetic case,we find that the 3d-bands are below the 4s-bands (at the Fermi level) by about 6 to 7 eV, which is in rough agreement with the earlier atomic HFA resultsWatson59 (14, 15). The splitting between the occupied and unoccupied 3d-states is about 23 eV, which is the on-site $`U`$ for Co. Magnetic HFA solutions are also found in the atomic limit for Co, as shown in the right panel of Fig. 15. The majority-spin 3d states ($`T_{2g}`$ and $`E_g`$) are now completely filled and energetically lie lower than the corresponding non-magnetic HFA-states. But only the (3-fold degenerate) $`T_{2g}`$-states of the minority-spin electrons are filled whereas the $`E_g`$-states of the minority electrons are empty (and now even 26 eV above the occupied d-states). The additional energetical shifts between the occupied 3d-states in the paramagnetic and ferromagnetic atomic HFA solution are due to the exchange matrix elements $`J`$.

This behavior can qualitatively be understood within the framework of the following simple, analytically solvable model. Similar to the numerical HFA-results presented and discused above, we neglect all intersite one-particle (hopping) and interaction matrix elements. Furthermore, we assume that we have diagonalized the one-particle Hamiltonian, taking into account only the atomic 3d-levels and assuming that the on-site one-particle diagonal matrix elements $`\epsilon `$, the Coulomb matrix elements $`U`$, and the exchange matrix elements $`J`$ are equal, i.e., that the 3d-levels are degenerate in the atomic limit with no crystal-field effects. Then the atomic part of the many-body Hamiltonian can be written as

$`H={\displaystyle \underset{i\sigma }{}}\epsilon c_{i\sigma }^{}c_{i\sigma }`$ $`+`$ $`{\displaystyle \frac{U}{2}}{\displaystyle \underset{(i\sigma )(j\sigma ^{})}{}}c_{i\sigma }^{}c_{i\sigma }c_{j\sigma ^{}}^{}c_{j\sigma ^{}}`$ (19)
$`+`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{ij,\sigma \sigma ^{}}{}}c_{i\sigma }^{}c_{j\sigma ^{}}^{}c_{i\sigma ^{}}c_{j\sigma }`$
where $`i,j\{0,\mathrm{},4\}`$ denote the 5 (degenerate) 3d-states.

The standard Hartree-Fock decoupling leads to

$`H={\displaystyle \underset{i\sigma }{}}(\epsilon `$ $`+`$ $`U\left[{\displaystyle \underset{j\sigma ^{}}{}}c_{j\sigma ^{}}^{}c_{j\sigma ^{}}c_{i\sigma }^{}c_{i\sigma }\right]`$ (20)
$``$ $`J{\displaystyle \underset{ji}{}}c_{j\sigma }^{}c_{j\sigma })c_{i\sigma }^{}c_{i\sigma }.`$
Here we have assumed that the Hartree-Fock Hamiltonian has the same symmetry as the uncorrelated Hamiltonian, and hence off-diagonal expectation values $`c_{j\sigma ^{}}^{}c_{i\sigma }`$ for $`(i\sigma )(j\sigma ^{})`$ vanish. From this equation it is clear that the HF Hamiltonian can be written in terms of an effective one-particle energy

$$H=\underset{i\sigma }{}\epsilon _{i\sigma }^{\text{HFA}}c_{i\sigma }^{}c_{i\sigma }$$
(21)
where

$$\epsilon _{i\sigma }^{\text{HFA}}=\epsilon +U\left[\underset{j\sigma ^{}}{}c_{j\sigma ^{}}^{}c_{j\sigma ^{}}c_{i\sigma }^{}c_{i\sigma }\right]J\underset{ji}{}c_{j\sigma }^{}c_{j\sigma }.$$
(22)
In the simple Hartree approximation (HA) the exchange decouplings are neglected, which means that all the decoupling terms with the negative sign would not occur. Therefore, the corresponding Hartree one-particle energies are given by

$$\epsilon _{i\sigma }^{\text{HA}}=\epsilon +U\underset{j\sigma ^{}}{}c_{j\sigma ^{}}^{}c_{j\sigma ^{}}.$$
(23)
Comparing this result with the Hartree-Fock one-particle energies, we find that the HF occupied levels are shifted downwards by an amount of

$$Uc_{i\sigma }^{}c_{i\sigma }+J\underset{ji}{}c_{j\sigma }^{}c_{j\sigma }$$
(24)
relative to the Hartree levels. Momentarily setting $`J=0`$, we see that for $`N`$ occupied levels the Hartree approximation gives the one-particle energies

$$\epsilon _{i\sigma }^{\text{HA}}=\epsilon +NU$$
(25)
whereas the HFA yields

$$\epsilon _{i\sigma }^{\text{HFA}}=\epsilon +(N1)U.$$
(26)
The occupied Hartree-Fock one-particle energies are lower than the corresponding Hartree one-particle energies by $`U`$, which is a consequence of the artificial and unphysical self-interaction still present in the Hartree approximation that is exactly cancelled in Hartree-Fock. This also explains why the Hartree-Fock bands are shifted downwards from the Hartree bands by an energy of the amount $`U`$. One also sees from this simple atomic-limit Hartree-Fock model that the energy difference between the highest occupied and the lowest unoccupied effective Hartree-Fock one-particle energies is again essentially $`U`$, which is once more in agreement with our numerical results for the crystal and for the atom (cf. Fig.15). Note that we have ignored $`U_{sd}`$ interactions, which cause an additional shift of d-bands below the s-bands by about an additional 10 eV in the full HFA calculations.

Taking into account the exchange interaction $`J`$ again and denoting by $`N_\sigma `$ the number of occupied states with spin $`\sigma `$ (i.e. $`N=N_{}+N_{}`$) one obtains in HFA

$$\epsilon _\sigma ^{\text{HFA}}=\epsilon +(N1)U(N_\sigma 1)J.$$
(27)
Then the total energy in HFA is given by

$$E_{\text{tot}}=N\epsilon +\frac{N(N1)}{2}U\underset{\sigma }{}\frac{N_\sigma (N_\sigma 1)}{2}J.$$
(28)
For the total energy we have added the necessary correction term to the sum of the occupied energy levels (much like the double counting term that shows up in band-structure calculations). Now for partially filled 3d-shells the occupation of the different spin directions may be different. Denoting $`M=N_{}N_{}`$ we obtain for the total energy

$$E_{\text{tot}}=N\epsilon +\frac{N(N1)}{2}U\frac{N^2+M^2}{4}J+\frac{N}{2}J.$$
(29)
The magnetic ($`M0`$) total energy is lower than the nonmagnetic (consistent with Hund’s rules).

Take once more Co with 8 3d-electrons. The paramagnetic (nonmagnetic) state has the occupations $`N_{}=N_{}=4`$ ($`N`$=8 and $`M`$=0). For this configuration (corresponding to the left panel in Fig. 15) one obtains

$$E_{\text{tot}}^{(P)}=8\epsilon +28U12J.$$
(30)
The Hund’s rule magnetic solution has 3d-states of one spin-direction completely filled, i.e., $`N_{}=5`$ and $`N_{}=3`$ ($`N`$=8 and $`M`$=2). This gives

$$E_{\text{tot}}^{(M)}=8\epsilon +28U13J.$$
(31)
Therefore, the magnetic configuation (with a magnetic moment of 2 for the atom) is energetically more favorable by $`J`$. Note also the exchange splitting in the occupied energy eigenvalues

$$\epsilon _{}\epsilon _{}=2J$$
(32)
and that our model would predict the unoccupied minority spin $`E_{2g}`$ state to be $`U+J`$ higher in energy than the corresponding occupied majority spin state.

The simple model in this section differs from the results shown in Fig. 15 in that we have replaced the full matrix of $`U`$ and $`J`$ by scalar values for d-states only (ignoring s-d interactions, for example). However, it captures all of the important physics without attempting to be completely quantitative.

## VII  Comparison with LSDA results

For comparison with the HFA results described in Section V we have also performed a standard local spin-density approximation (LSDA), ab-initio, band-structure calculation with the LMTO-ASA method. The results for the DOS are shown in Figs. 1619. For the magnetic systems Fe, Co, and Ni one obtains also an exchange splitting and the prediction of magnetic solutions with magnetic moments of 2.18 for iron, 1.58 for Co, and 0.58 for Ni, which are in better agreement with experiment than the HFA results. This is consistent with the expectation that L(S)DA is tuned to accurately predict ground-state properties.

The energy spectra of the bands (DOS) are quite different from the HFA. For example, the 3d-bands now fall into the same energy region as the 4sp-bands, i.e., the LSDA-results are not so different from the Hartree-results. This means that the exchange-correlation energy leads only to a small shift of the 3d-bands downwards by at most a few eV and a smaller exchange splitting (also of the magnitude of 1 eV) in the case of magnetic solutions. On the other hand, in the LSDA calculations the self-interaction terms are not completely canceled, i.e., an (unrealistic) self-interaction is included, which may lead to 3d bands that lie energetically too high, as discussed for the atomic limit in the previous section.

To see the effect of correlations within LSDA, we have also performed an exchange-only calculation for Co. The result is shown in Fig. 20. Obviously the (majority) d-bands lie lower than the ones in the full LSDA-calculation shown in Fig. 17. But the shift is of the magnitude of 1 eV only, i.e., very minimal when compared with the large drop in the full HFA. On the other hand, this exchange-only LSDA-result also does not contain self-interaction corrections, which are responsible for the large shift downwards of the d-bands. Nevertheless, the LSDA result indicates that a possible effect of correlations is to shift the 3d-bands up relative to exchange-only calculations, and hence one would expect a similar effect if correlations could be added to the full HFA-calculations.

One can also calculate the total energy in the Hartree, HFA, and LSDA approximations. The results obtained for the four materials of interest are shown in Fig. 21. We see that the total energy is always significantly lower in HFA than in the Hartree approximation, which is expected because the HFA minimizes the total energy. The HFA total energy is also lower than the L(S)DA, and the LSDA result is lower than the simple Hartree result. Because of the unknown approximations that go into constructing L(S)DA, it is hard to guess ahead of time that this would be the case. However, it is well known that the L(S)DA approximation produces a bad total exchange-correlation energy; the reason why such good agreement with experiment is found is that relative exchange-correlation energies are nonetheless reasonably accurately calculated.

## VIII  Discussion and conclusion

In this paper we have presented the results of (unscreened) HFA calculations for the 3d transition metals Fe, Co, Ni and Cu. We obtain magnetic solutions for Fe, Co, and Ni with (slightly) too large magnetic moments when compared to experimental or LSDA results. The occupied HFA 3d-bands lie about 20 eV below the Fermi level (and the Hartree result), which is also the magnitude of the splitting between occupied and unoccupied 3d-bands and of the magnitude of the on-site Coulomb matrix element (the “Hubbard” U). This downwards shift of the HFA 3d-bands compared to the Hartree- and L(S)DA-3d-bands can be understood as due to the self-interaction correction of HFA.

One may argue that these results are not surprising and an artifact of using the unscreened HFA. Our ab-initio calculation of the direct Coulomb matrix elements yields large values of the magnitude of 20 eV. Therefore the Coulomb energies are large for these materials and HFA can be considered to be an approximation for the selfenergy which is correct only in linear order in the Coulomb interaction. But for these large values of the U-terms HFA is certainly not sufficient but one has to apply better many-body approximations. One should apply systematic extensions of HFA, which within the standard perturbational approach can be represented by (a resummation of an infinite series of) Feynman diagrams, or one can try to apply the recently so successful non-perturbational many-body schemes like “dynamical mean field theory” (DMFT)GKKR96 (16) or variational (Gutzwiller) approachesGebhardWeber (17). The simplest standard diagram series are the bubble diagrams leading essentially to the ”random phase approximation” (RPA). This means just a renormalization of the interaction line, i.e. the pure “naked” Coulomb interaction has to be replaced by a “dressed” interaction. Or in other words, the exchange (Fock) contribution has not to be calculated with the bare Coulomb matrix elements but with screened Coulomb matrix elements. Probably the non-perturbational schemes like DMFT are also only applicable for screened Coulomb matrix elements.

To summarize, the new many-body ab-initio approach proposed in Ref. Schnell2003-1, 1 is applicable to real materials. Maximally localized Wannier functions and their one- and two-particle matrix elements can be calculated from first principles so as to obtain a many-body multi-band Hamiltonian in second quantization. Standard methods of many-body theory can then be applied to this Hamiltonian. The simplest approximation (HFA) is unreliable, because the large Coulomb matrix elements are unscreened. Therefore, in the future better many-body approximations should be used. Within the standard Feynman diagram approach the most straightforward next step consists in a summation of bubble diagrams leading to a renormalized (screened) Coulomb interaction. This requires calculating the exchange contribution not with the bare but with a screened Coulomb interaction. To take into account the effects of screening one has to calculate the charge susceptibility and the (static) dielectric constant, which can be done within a generalized Lindhard theory, for instance.

###### Acknowledgements.

This work has been supported by a grant from the Deutsche Forschungsgemeinschaft No. Cz/31-12-1. It was also partially supported by the Department of Energy under contract W-7405-ENG-36. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.