text:
# FZJ-IKP(TH)-2001-02 Chiral Lagrangians at finite density

\[

## Abstract

The effective $`SU(2)`$ chiral Lagrangian with external sources is given in the presence of non-vanishing nucleon densities by calculating the in-medium contributions of the chiral pion-nucleon Lagrangian. As a by product, a relativistic quantum field theory for Fermi many-particle systems at zero temperature is directly derived from relativistic quantum field theory with functional methods.

PACS numbers: 12.39.Fe, 24.85.+p, 13.75.L, 13.75.Gx

\]

1. In the limit of massless up and down quarks the QCD Lagrangian is symmetric under the chiral group $`SU(2)_L\times SU(2)_R`$. One assumes that this symmetry is spontaneously broken to the diagonal subgroup $`SU(2)_{L+R}`$ giving rise to the appearance of 3 massless Goldstone bosons which finally acquire small masses due to the non-vanishing mass of the $`u`$ and $`d`$ quarks.

This symmetry breaking scenario constrains so much the interactions of the Goldstone bosons that the QCD Green functions can be calculated at low energies as an expansion in powers of momenta and quark masses. This is known as Chiral Perturbation Theory.

The extension of the theory to the case of low temperature at zero density was considered in ref.. In this article we study the case of small densities at zero temperature and derive the corresponding chiral Lagrangian by calculating the in-medium contributions due to the chiral pion-nucleon Lagrangian with functional methods. Although we focus our treatment to QCD, the relativistic many-body formalism here deduced for Fermi systems can be applied to processes governed by other dynamical theories, as the traditional non-relativistic zero temperature many-body  quantum theory which stresses the diagrammatic approach. Compared with standard quantum field theory at finite temperature $`T`$ in the grand canonical ensemble, one avoids the use of unknown chemical potentials which themselves have to be calculated in terms of the many-body forces. The former is accomplished by following quantum field theory at $`T=0`$ considering directly the change of the ground state from the vacuum to one with finite fermionic densities. In the same way one also avoids the non-trivial $`T0`$ limit due to the so called anomalous diagrams. The price to pay is to rely on the adiabatic hypothesis in order to determine the interacting ground state from that of the free case by turning on the interactions adiabatically.

2. Let us take first the case of symmetric and unpolarized nuclear matter, the extension of the formalism to the asymmetric and polarized case is straightforward and will be shown below. In the following we take the Heisenberg picture and, following closely scattering theory, we consider two ground states $`|\mathrm{\Omega }_{\text{out}}`$ and $`|\mathrm{\Omega }_{\text{in}}`$ which, under the action of any time dependent operator at asymptotic times $`t\pm \mathrm{}`$, respectively, behave as two symmetric Fermi seas of free protons and neutrons. The Fermi seas are filled up to the corresponding baryonic density, $`{\displaystyle \stackrel{N}{}}`$$`a^{}(𝕡_n)|0`$, where the label $`n`$ includes also the spin and isospin indices, $`N`$ is the number of momentum states inside the Fermi sea with Fermi momentum $`k_F=(3\pi ^2\rho /2)^{1/3}`$, $`\rho `$ is the total nuclear density and $`|0`$ is the vacuum. Our objective is to evaluate the generating functional $`𝒵[v,a,s,p]`$ in the presence of vector $`v_\mu `$, axial $`a_\mu `$, scalar $`s`$ and pseudoscalar $`p`$ external fields by working out the transition amplitude $`\mathrm{\Omega }_{\text{out}}|\mathrm{\Omega }_{\text{in}}_J`$, where the label $`J`$ just indicates the presence of the aforementioned external sources. In this way by taking functional derivatives of $`𝒵[v,a,s,p]`$ with respect to the external sources one evaluates the in-medium QCD connected Green functions (space-time averages at finite density of the quark currents coupled to the $`v`$, $`a`$, $`s`$ and $`p`$ sources). To do this we consider the effective chiral Lagrangians $`=_{\pi \pi }+_{\overline{\psi }\psi }+_{\overline{\psi }\psi \overline{\psi }\psi }+\mathrm{}`$ with increasing number of pairs of nucleon fields $`\psi (x)`$. We first restrict ourselves to the term with no nucleon fields $`_{\pi \pi }`$ and to that containing two of them $`_{\overline{\psi }\psi }=\overline{\psi }(x)D(x)\psi (x)`$, together with the previous external fields. We will discuss later a way to include perturbatively the contributions of Lagrangians with higher number of nucleons by considering them to arise from bilinear vertices through the exchange of an arbitrary heavy particle. Indeed, although we are talking about CHPT, the only thing that matters for the following derivations is that $`_{\overline{\psi }\psi }`$ is bilinear in the fermions. Consider now the transition amplitude for the ground states from $`t\mathrm{}`$ to $`t+\mathrm{}`$ in the presence of the previous external sources together with Grassmann sources $`\eta `$ and $`\eta ^{}`$, coupled to the nucleon fields:

$`\mathrm{\Omega }_{\text{out}}|\mathrm{\Omega }_{\text{in}}_{J,\eta ,\eta ^{}}={\displaystyle }[dU][d\psi ][d\psi ^{}]\mathrm{\Omega }_{\text{out}}|\psi (+\mathrm{})\times `$ (1)
$`e`$ $`{}_{}{}^{i{\scriptscriptstyle 𝑑x\left[_{\pi \pi }+\overline{\psi }D\psi +\eta ^{}\psi +\psi ^{}\eta \right]}}\psi (\mathrm{})|\mathrm{\Omega }_{\text{in}}`$ (2)
with the pion fields described by the $`2\times 2`$ unitary matrix $`U`$.

The ground state functional $`\psi (\pm \mathrm{})|\mathrm{\Omega }_{\begin{array}{c}\text{out}\\ \text{in}\end{array}}`$ can be expressed in terms of that of the vacuum by writing the $`a(𝕡_n)_{\begin{array}{c}\text{out}\\ \text{in}\end{array}}`$ operators as a function of $`\psi (x)`$ and its time derivative $`\dot{\psi }(x)`$:

$`a(𝕡_n)_{\begin{array}{c}\text{out}\\ \text{in}\end{array}}=\underset{t\pm \mathrm{}}{lim}{\displaystyle \frac{E(p)}{2m_N}}e^{iE(p)t}{\displaystyle }d𝕩e^{i𝕡𝕩}\overline{u}(𝕡_n)[\psi _n(x)`$ (5)
$`+{\displaystyle \frac{i}{E(p)}}\dot{\psi }_n(x)],`$ (6)
where as usual in scattering theory for $`t\pm \mathrm{}`$ the matrix elements are calculated as if there were no interactions. In the former expression, $`E(p)`$ is the energy of the nucleon with three-momentum $`𝕡`$, $`m_N`$ is the nucleon mass and $`u(𝕡_n)`$ is a Dirac spinor. Expressing $`\dot{\psi }(x)`$ in terms of $`\psi (x)`$ (one way is using the Dirac equation), taking into account that $`\psi (x)e^{i{\scriptscriptstyle 𝑑y\eta ^{}(y)\psi (y)}}=i\frac{\stackrel{}{\delta }}{\delta \eta ^{}(x)}e^{i{\scriptscriptstyle 𝑑y\eta ^{}(y)\psi (y)}}`$, the analogous expression for $`\psi ^{}(x)`$, and substituting all that in eq.(1), one obtains for $`\mathrm{\Omega }_{\text{out}}|\mathrm{\Omega }_{\text{in}}_{J,\eta ,\eta ^{}}`$:

$`\underset{\begin{array}{c}t+\mathrm{}\\ t^{}\mathrm{}\end{array}}{lim}{\displaystyle }[dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}({\displaystyle \underset{n}{\overset{N}{}}}{\displaystyle \frac{E(p_n)}{2m_N}}{\displaystyle }d𝕩_ne^{ip_nx_n}`$ (13)
$`\overline{u}(𝕡_n)[1{\displaystyle \frac{i}{E(p_n)}}(\gamma ^0{\displaystyle \underset{j=1}{\overset{3}{}}}\gamma ^j{\displaystyle \frac{}{x_n^j}}+i\gamma ^0m_N)]{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta ^{}(x_n)}})`$
$`e^{i{\scriptscriptstyle 𝑑x𝑑y\eta ^{}(x)D^1(x,y)\gamma ^0\eta (y)}}{\displaystyle \underset{m}{\overset{N}{}}}{\displaystyle \frac{E(q_m)}{2m_N}}{\displaystyle 𝑑𝕪_me^{iq_my_m}}`$
$`[{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta (y_m)}}{\displaystyle \frac{i}{E(q_m)}}({\displaystyle \underset{k=1}{\overset{3}{}}}{\displaystyle \frac{}{y_m^k}}{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta (y_m)}}\gamma ^k\gamma ^0+i{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta (y_m)}}`$
$`\gamma ^0m_N)]\gamma ^0u(𝕢_m)`$
where the integration over the nucleon fields and conjugate momenta is also done. Furthermore we define $`x_n=(t,𝕩_n)`$, $`y_m=(t^{},𝕪_m)`$, $`p_n=(E(p_n),𝕡_n)`$ and analogously for $`q_m`$. The action of the spatial derivatives can be readily taken into account by integrating by parts. In this way, they only act on the corresponding exponentials giving rise to three-momenta factors that can be further simplified by applying the Dirac equation on the Dirac spinors. In this way we have:

$`\overline{u}(𝕡_n)\left[1{\displaystyle \underset{j=1}{\overset{3}{}}}{\displaystyle \frac{p_n^j}{E(p_n)}}\gamma ^0\gamma _j+{\displaystyle \frac{m_N}{E(p_n)}}\gamma ^0\right]`$ $`=`$ $`{\displaystyle \frac{2m_N}{E(p)}}u^{}(𝕡_n)`$ (14)
$`\left[1{\displaystyle \underset{j=1}{\overset{3}{}}}{\displaystyle \frac{p_n^j}{E(p_n)}}\gamma ^0\gamma _j+{\displaystyle \frac{m_N}{E(p_n)}}\gamma ^0\right]\gamma ^0u(𝕡_n)`$ $`=`$ $`{\displaystyle \frac{2m_N}{E(p)}}u(𝕡_n)`$ (15)
Applying the previous results to eq.(13) it simplifies to:

$`\underset{\begin{array}{c}t+\mathrm{}\\ t^{}\mathrm{}\end{array}}{lim}{\displaystyle }[dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}({\displaystyle \underset{n}{\overset{N}{}}}{\displaystyle }d𝕩_ne^{ip_nx_n}u^{}(𝕡_n)`$ (18)
$`{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta ^{}(x_n)}})e^{i{\scriptscriptstyle 𝑑x𝑑y\eta ^{}(x)D^1(x,y)\gamma ^0\eta (y)}}{\displaystyle }_m^N{\displaystyle }d𝕪_me^{iq_my_m}`$ (19)
$`{\displaystyle \frac{\stackrel{}{\delta }}{\delta \eta (y_m)}}u(𝕢_m)`$ (20)
After acting with the left derivatives $`\stackrel{}{\delta }/\delta \eta ^{}(x_n)`$ on the exponential depending on the Grassmann sources, the former expression can be recast as:

$`\underset{\begin{array}{c}t+\mathrm{}\\ t^{}\mathrm{}\end{array}}{lim}{\displaystyle }[dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}({\displaystyle \underset{n}{\overset{N}{}}}{\displaystyle }d𝕩_ne^{ip_nx_n}{\displaystyle }dz_n^{}`$ (23)
$`u^{}(𝕡_n)D^1(x_n,z_n^{})\gamma ^0\eta (z_n^{}))e^{i{\scriptscriptstyle 𝑑x𝑑y\eta ^{}(x)D^1(x,y)\gamma ^0\eta (y)}}`$ (24)
$`{\displaystyle \underset{m}{\overset{N}{}}}{\displaystyle 𝑑𝕪_me^{iq_my_m}\frac{\stackrel{}{\delta }}{\delta \eta (y_m)}u(𝕢_m)}.`$ (25)
This result is equal to $`e^{i𝒵[v,a,s,p]}`$ when $`\eta 0`$ and $`\eta ^{}0`$.***Without taking the limit $`\eta ,\eta ^{}0`$, the generating functional contains baryonic sources and taking differential derivatives with respect to them one could directly evaluate in-medium baryonic Green functions, e.g. nucleon propagators. Nevertheless, since for our present purposes the baryons in the medium constitute just a background, we do not consider this case any further. Hence we can equalize to 1 the second exponential from the left and the remaining right derivatives and $`\eta (z_n)`$ sources have to be paired in order to finish with a non-vanishing result. Thus we have:

$`e^{i𝒵}=\underset{\begin{array}{c}t+\mathrm{}\\ t^{}\mathrm{}\end{array}}{lim}{\displaystyle [dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}\underset{\sigma }{}ϵ(\sigma )\underset{n}{\overset{N}{}}𝑑𝕩_n}`$ (28)
$`{\displaystyle 𝑑𝕩_{\sigma _n}e^{ip_nx_n}u^{}(𝕡_n)D^1(x_n,x_{\sigma _n})\gamma ^0e^{ip_{\sigma _n}x_{\sigma _n}}u(𝕡_{\sigma _n})}`$ (29)
with $`ϵ(\sigma )`$ the signature of the permutation $`\sigma `$ over all the indices –momenta, spin and isospin. In order to continue let us write the operator $`D(x)D_0(x)A(x)`$ with $`D_0(x)=i\gamma ^\mu _\mu m_N`$ the Dirac operator for the free motion of the nucleons. On the other hand, the operator $`A(x)`$ is completely general although in our case at hand, CHPT, it is subject to a chiral expansion of powers of soft three-momenta and quark masses. Furthermore, let us note that:

$`\underset{t+\mathrm{}}{lim}{\displaystyle 𝑑𝕩e^{ipx}u^{}(𝕡)D_0^1(x,x^{})}=ie^{ipx^{}}\overline{u}(𝕡)`$ (30)
$`\underset{t\mathrm{}}{lim}{\displaystyle 𝑑𝕩D_0^1(x^{},x)\gamma ^0e^{ipx}u(𝕡)}=ie^{ipx^{}}u(𝕡)`$ (31)
This result can be easily obtained by writing $`D_0^1(x,x^{})`$ in four-momentum space and then performing the integral over the temporal component of the momentum taking care of the imposed limits. For instance, let us take the first of the previous equations. Then we have:

$`\underset{t+\mathrm{}}{lim}{\displaystyle 𝑑𝕩e^{ipx}u^{}(𝕡)D_0^1(x,x^{})}=`$ (32)
$`\underset{t+\mathrm{}}{lim}{\displaystyle 𝑑𝕩e^{ipx}u^{}(𝕡)\frac{dR}{(2\pi )^4}\frac{e^{iR(xx^{})}(\overline{)𝑅}+m_N)}{R^2m_N^2+iϵ}},`$ (33)
with $`ϵ`$ a positive infinitesimal. Exchanging the order of the integrations, the spatial one gives rise to $`(2\pi )^3\delta (𝕡)`$ which fixes $``$. As a result one has:

$`\underset{t+\mathrm{}}{lim}{\displaystyle \frac{dR^0}{2\pi }u^{}(𝕡)\frac{e^{ip^0t}e^{iR^0(tx_0^{})}e^{i𝕡𝕩^{}}}{R_0^2𝕡^2m_N^2+iϵ}\left(R^0\gamma ^0𝕡\gamma +m_N\right)}`$ (34)
Since $`t+\mathrm{}`$ then $`tx_0^{}>0`$ and we close the integration contour over $`R^0`$ with a semi-circle of infinite radius on the lower half-plane picking up the pole at $`R^0=p^0iϵ=E(p)iϵ`$. Applying the Dirac equation to the result one arrives to eq.(30).

Then taking into account eq.(30) and the expansion $`D^1\gamma ^0=[D_0A]^1\gamma ^0=D_0^1\gamma ^0+D_0^1AD_0^1\gamma ^0+\mathrm{}`$, we can rewrite eq.(28) as:

$`e^{i𝒵}={\displaystyle [dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}\underset{\sigma }{}ϵ(\sigma )\underset{n}{\overset{N}{}}𝑑𝕩_n𝑑𝕩_{\sigma _n}}`$ (35)
$`e^{i𝕡_n𝕩_n}u^{}(𝕡_n)[\delta (𝕩_n𝕩_{\sigma _n})\delta _{n,\sigma _n}i{\displaystyle }dt{\displaystyle }dt^{}e^{iE(p_n)t}\gamma ^0`$ (36)
$`A(x_n,x_{\sigma _n})e^{iE(p_{\sigma _n})t^{}}i{\displaystyle 𝑑t𝑑t^{}𝑑z𝑑z^{}e^{iE(p_n)t}\gamma ^0}`$ (37)
$`A(x_n,z)D_0^1(z,z^{})A(z^{},x_{\sigma _n})e^{iE(p_{\sigma _n})t^{}}+\mathrm{}]`$ (38)
$`e^{i𝕡_{\sigma _n}𝕩_{\sigma _n}}u(𝕡_{\sigma _n})`$ (39)
where now both $`t`$ and $`t^{}`$ are integration variables being the time components of $`x_n`$ and $`x_{\sigma _n}`$, respectively. The dots just refer to those terms with an increasing number of insertions of the operator $`A`$ coming from the geometric expansion of $`D^1=\left[D_0A\right]^1`$ discussed above. Furthermore $`A(x,y)`$ is defined such that

$`{\displaystyle 𝑑x\overline{\psi }(x)A(x)\psi (x)}={\displaystyle 𝑑x𝑑y\overline{\psi }(x)A(x,y)\psi (y)}.`$ (40)
As a result of eq.(35) we can simply state that:

$$e^{i𝒵[v,a,s,p]}=[dU](\text{det}D)e^{i{\scriptscriptstyle 𝑑x_{\pi \pi }}}(\stackrel{~}{\text{det}})$$
(41)
where the tilde in $`\stackrel{~}{\text{det}}`$ indicates that the determinant has to be taken in the subspace of the Fermi sea states expanded by the basis functions $`e^{i𝕡_n𝕩}u(𝕡_n)`$ with $`|𝕡_n|<k_F`$. In this notation $``$ is given by:

$`I_3i{\displaystyle 𝑑t𝑑t^{}e^{iH_0t}\gamma ^0A\left[I_4D_0^1A\right]^1e^{iH_0t^{}}}`$ (42)
where $`I_3\delta (𝕩_n𝕩_m)\delta _{n,m}`$ and analogously $`I_4\delta (x_nx_m)\delta _{n,m}`$. On the other hand, $`e^{iH_0t^{}}e^{i𝕡_n𝕩^{}}u(𝕡_n)=e^{ip_nx^{}}u(𝕡_n)`$ and $`e^{i𝕡_n𝕩}u^{}(𝕡_n)e^{iH_0t}=e^{ip_nx}u^{}(𝕡_n)`$.

In order to obtain from eq.(41) the contributions of the surrounding medium to the generating functional it is convenient to exponentiate $`\stackrel{~}{\text{det}}`$ as $`\mathrm{exp}(\stackrel{~}{\text{Tr}}\mathrm{log})`$ (where the tilde has the same meaning as before). Then we have:

$`e^{i𝒵[v,a,s,p]}={\displaystyle }[dU](\text{det}D)\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }+{\displaystyle \underset{r=1}{\overset{2}{}}}{\displaystyle \underset{\alpha =1}{\overset{2}{}}}`$ (43)
$`{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}{\displaystyle 𝑑𝕩𝑑𝕪e^{i𝕡𝕩}u_{r\alpha }^{}(𝕡)\mathrm{log}[]}|_{(𝕩\alpha ,𝕪\alpha )}`$ (44)
$`e^{i𝕡𝕪}u_{r\alpha }(𝕡)\}`$ (45)
where we have indicated explicitly the spin $`r`$ and isospin $`\alpha `$ indices. It is very appropriate to stress at this point that eq.(43), although formal, is non-perturbative. From this result we can readily read out the new effective chiral Lagrangian density $`\stackrel{~}{}_{\pi \pi }`$ in the presence of a nuclear density just by equating the expression between curly brackets to $`i𝑑x\stackrel{~}{}_{\pi \pi }`$.

3. The perturbative theory is obtained by expanding $`\mathrm{log}`$ in eq.(43), so that $`e^{i𝒵}`$ can be written as:

$`{\displaystyle }[dU](\text{det}D)\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }i{\displaystyle \underset{r=1}{\overset{2}{}}}{\displaystyle \underset{\alpha =1}{\overset{2}{}}}{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}`$ (46)
$`{\displaystyle 𝑑x𝑑ye^{ipx}\overline{u}_{r\alpha }(𝕡)A[I_4D_0^1A]^1}|_{(x\alpha ,y\alpha )}e^{ipy}u_{r\alpha }(𝕡)`$ (47)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{r,s=1}{\overset{2}{}}}{\displaystyle \underset{\alpha ,\beta =1}{\overset{2}{}}}{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕢}{(2\pi )^32E(q)}}{\displaystyle 𝑑x𝑑x^{}𝑑y^{}𝑑y}`$ (48)
$`e^{ipx}\overline{u}_{r\alpha }(𝕡)A[I_4D_0^1A]^1|_{(x\alpha ,x^{}\beta )}e^{iqx^{}}u_{s\beta }(𝕢)e^{iqy^{}}\overline{u}_{s\beta }(𝕢)`$ (49)
$`A[I_4D_0^1A]^1|_{(y^{}\beta ,y\alpha )}e^{ipy}u_{r\alpha }(𝕡)+\mathrm{}\}.`$ (50)
Finally, taking into account the relation $`{\displaystyle \underset{r=1}{\overset{2}{}}}`$$`u_r(𝕡)\overline{u}_r(𝕡)=\overline{)𝑝}+m_N`$ we arrive at the following expression for the generating functional:

$`e^{i𝒵}={\displaystyle }[dU](\text{det}D)\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }i{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}`$ (51)
$`{\displaystyle 𝑑x𝑑ye^{ip(xy)}\text{Tr}\left(A[I_4D_0^1A]^1|_{(x,y)}(\overline{)𝑝}+m_N)\right)}+`$ (52)
$`{\displaystyle \frac{1}{2}}{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}{\displaystyle ^{k_F}}{\displaystyle \frac{d𝕢}{(2\pi )^32E(q)}}{\displaystyle 𝑑x𝑑x^{}𝑑y𝑑y^{}e^{ip(xy)}}`$ (53)
$`e^{iq(x^{}y^{})}\text{Tr}(A[I_4D_0^1A]^1|_{(x,x^{})}(\overline{)𝑞}+m_N)A[I_4`$ (54)
$`D_0^1A]^1|_{(y^{},y)}(\overline{)𝑝}+m_N))+\mathrm{}\}`$ (55)
where the trace refers both to the isospin and spinor indices. The previous formula implies a double expansion for obtaining the contributions of the nuclear medium to the chiral Lagrangian. One is the standard chiral expansion by expanding the vacuum operator $`A[I_4D_0^1A]^1=A+AD_0^1A+AD_0^1AD_0^1A+\mathrm{}`$ together with $`A`$ itself, $`A=A^{(1)}+A^{(2)}+A^{(3)}+\mathrm{}`$ in increasing powers of momenta and quark masses valid at low energies. Explicit expressions of $`A`$ up to $`𝒪(p^3)`$ can be obtained from ref.. The another is an expansion in the number of insertions of on-shell fermions belonging to the Fermi sea, schematically indicated in fig.1 by a thick solid line. Both expansions can be related by giving a chiral power couting to $`k_F`$ which can be naturally counted as $`𝒪(p)`$ since for nuclear saturation density $`k_F2M_\pi `$ with $`M_\pi `$ the pion mass. Moreover, the circles labeled by $`\mathrm{\Gamma }`$ correspond to the non-local operator $`iA[I_4D_0^1A]^1`$. Note that when inserting $`n`$ Fermi seas from the expansion of the logarithm one picks up a factor $`(1)^n/n`$ where the global minus sign appears due to the fermionic closed loop, $`n`$ is a combinatoric factor because any cyclic permutation in the trace of $`n`$ Fermi seas with their associated $`n`$ $`\mathrm{\Gamma }`$ operators gives the same result and finally the sign $`(1)^n`$ is a pure in-medium factor that one has to keep in mind and is already present in standard many-body theory.

A generally non-local vertex $`\mathrm{\Gamma }`$ comes from the iteration of the $`A`$ operator with intermmediate free baryon propagators $`D_0^1`$, obeying the usual Feynman rules, see fig.2. Notice that $`A`$ was defined from the Lagrangain $`\overline{\psi }D\psi `$ removing the free term $`\overline{\psi }D_0\psi `$ and changing the sign to the rest, this is why a minus sign appears in front of $`A`$ in fig.2. On the other hand $`iD_0^1(x,y)=i\frac{d^4p}{(2\pi )^4}\frac{\overline{)𝑝}+m_N}{p^2m_N^2+iϵ}`$ is the usual baryon propagator.

Hence a final diagram, when expanding $`\mathrm{\Gamma }`$ up to the required accuracy, will be a set of $`n1`$ Fermi-sea insertions, $`m0`$ free baryon-propagators and of $`m+n`$ vertices $`iA`$. First we include the $`1`$ global sign because of the fermionic closed loop, and the combinatoric factor $`1/n`$ together with the sign $`(1)^n`$. Then, following the diagram in the opposite sense to that of the fermionic arrows, write for each Fermi-sea an integral $`^{k_F}\frac{d𝕡(\overline{)𝑝}+m)}{(2\pi )^32E(p)}`$ with $`p^0=E(𝕡)`$, for each vaccum baryon propagator with free momentum $`p`$ write $`i\frac{dp}{(2\pi )^4}\frac{\overline{)𝑝}+m_N}{p^2m_N^2+iϵ}`$ and for a vertex in momentum space a term $`iA(2\pi )^4`$, keeping in mind the energy-momentum conservation at each vertex. Finally sum over the spin and isospin indices of the fermions. This defines explicitely the Feyman rules in momentum space in order to obtain $`i(2\pi )^4`$ times the desired connected graph accompied by the global delta of energy-momentum conservation. Analogous Feynman rules hold of course in configuration space, e.g., see eq.(51).

An equivalent way to state the previous rules, without including explicitely the integral symbols and factors $`2\pi `$, is to write the same sign-combinatoric factor $`(1)^n/n`$ and vertices $`iA`$ as before. Then for the free nucleons one has just the free propagator $`i\frac{\overline{)𝑝}+m}{p^2m^2+iϵ}`$ and for the Fermi-sea baryons the factor $`(2\pi )\delta (p^2m^2)\theta (p^0)(\overline{)𝑝}+m)\theta (k_F|𝕡|)`$. Finally sum over all the discrete indices attached to the fermions and integrate over all the free four-momenta with the measure $`d^4p/(2\pi )^4`$ after taking into account energy-momentum conservation at each vertex.

Fig.1 fixes the skeleton structure of standard in-medium vertices since still one has to consider the pion fields contained in $`A`$ over which one has to integrate in eqs.(43), (51) to finally obtain the generating functional. That is, from the vertices $`A`$ as well as from $`_{\pi \pi }`$, one can generate internal as well as external (coupled to the sources) pionic legs denoted by dashed lines in figs.3a and 3b. Here one has essentially the same Feynamn rules than in vaccum in order to proceed in a perturbative way, simply for each pionic line with four-momentum $`q`$, one writes the vacuum propagator $`i\frac{dq}{(2\pi )^4}\frac{1}{q^2M_\pi ^2+iϵ}`$, for the first version of the in-medium Feynman rules. For the second one has $`\frac{i}{q^2M_\pi ^2+iϵ}`$. The important remark to keep in mind is that a generalized vertex has analogous properties to those of a standard local quantum field theory one to the effects of determing the numerical factors accompanying the exchange of pion lines inside a given diagram. This can be seen just by applying standard perturbative techniques in path integrals to the action given between brackets in eq.(43) or more explicitely in eq.(51). Several examples are discussed in detail in ref..

4. Notice also that ultra-violet parts in the integration over running pionic momenta generate local multi-nucleon vertices and hence a proper treatment of pion loops can only be done by including simultaneously local nucleon interactions in order to reabsorb the divergences. This is denoted by the squares in figs. 3c and 3d. It is important to realize that while in fig.3a the pion is exchanged inside the same vertex (in the generalized sense of above) in fig.3b the pions are exchanged between two of them. This implies the presence of only one Dirac trace in fig.3c (the flow of the Dirac indices along the propagators on both sides of the square is indicated by the arrows of the open solid lines) and two Dirac traces in fig.3d. Thus, when writing the matrix elements corresponding to diagrams of the type of figs.3c and 3d (with local interactions) one has to keep track of the number of closed traces in the spinor indices since each of them will lead to a sign $`(1)`$ and its own combinatoric factor $`1/n_j`$ (with $`n_j`$ the number of Fermi-sea insertions in the closed Dirac loop) together with the sign $`(1)^{n_j}`$. One can also arrive to the same conclusions by considering that the dashed lines originate from the exchange of an arbitrary heavy meson (the only relevant point in our derivations is the bilinear character of the meson-$`\overline{N}N`$ interaction in the nucleon fields). In this way it is straightforward to realize the presence of a factor $`1/2`$ in front of fig.3d together with the factors $`_j(1)^{n_j+1}/n_j`$ which here is simply $`1`$. The relation between contact four nucleon interactions at low energies in effective field theories and its saturation by integrating-out heavy meson resonances has been established in ref..

There is still an important difference to be discussed when comparing figs.3c and 3d, which in fact is related to the presence of the factor (det$`D`$) in all the formulae from eq.(13) to eq.(51). The latter corresponds to contributions to the chiral Lagragian from closed fermion loops in the vacuum and in the spirit of the effective field theories these contributions, coming from states with masses close to or above the chiral scale $`\mathrm{\Lambda }_\chi M_\rho `$, are incorporated in the counterterms of the vacuum effectrive field theory. In this way we will set det$`D=1`$ in the following. Indeed on the right hand side of the square of fig.3c we can recognize a momentum loop flowing along the free baryon propagators without any insertion of a Fermi-sea. This simply means that since all the momenta that could go into this closed momentum loop are soft because they come from the circles, it just corresponds to vacuum renormalization from heavy particles and is reabsorbed in the low energy counterterms. This is schematically shown in fig.3e with a diamond corresponding to a higher order $`\pi N`$ counterterm. This is consistent as far as one is restricted to the low energy and momentum regime. However, when a Fermi-sea baryon line is present the momentum running though the baryon lines is of the form $`p_j=p+Q_j`$ with $`p=(E(p),𝕡)`$ and $`|𝕡|<k_F`$. Thus one is inserting a medium parameter $`k_F`$ and shifting upward up to around $`m_N`$ the energy level around which one is perturbating.

5. We now turn to the generalization of the formalism to the case of asymmetric nuclear matter, with different densities of neutrons, $`\rho _n`$, and protons, $`\rho _p`$, with Fermi momenta $`k_F^{(n)}=(3\pi ^2\rho _n)^{1/3}`$ and $`k_F^{(p)}=(3\pi ^2\rho _p)^{1/3}`$, respectively. Following the previous derivation of eq.(43) one can easily convince oneself that the only change is to remove the sum over isospin indices and to distinguish between $`\alpha =1`$ (proton) and $`\alpha =2`$ (neutron). In this way we will have:

$`e^{i𝒵[v,a,s,p]}={\displaystyle }[dU]\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }+{\displaystyle \underset{r=1}{\overset{2}{}}}`$ (56)
$`{\displaystyle ^{k_F^{(p)}}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)_1}}{\displaystyle 𝑑𝕩𝑑𝕪e^{i𝕡𝕩}u_{r\mathrm{\hspace{0.17em}1}}^{}(𝕡)\mathrm{log}[]}|_{(𝕩\mathrm{\hspace{0.17em}1},𝕪\mathrm{\hspace{0.17em}1})}`$ (57)
$`e^{i𝕡𝕪}u_{r\mathrm{\hspace{0.17em}1}}(𝕡)+{\displaystyle \underset{r=1}{\overset{2}{}}}{\displaystyle ^{k_F^{(n)}}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)_2}}{\displaystyle 𝑑𝕩𝑑𝕪e^{i𝕡𝕩}}`$ (58)
$`u_{r\mathrm{\hspace{0.17em}2}}^{}(𝕡)\mathrm{log}[]|_{(𝕩\mathrm{\hspace{0.17em}2},𝕪\mathrm{\hspace{0.17em}2})}e^{i𝕡𝕪}u_{r\mathrm{\hspace{0.17em}2}}(𝕡)\}`$ (59)
Notice that we have indicated separately the energies of protons $`E(p)_1`$ and neutrons $`E(p)_2`$ with three-momentum $`𝕡`$, since the previous equation is valid for the non-isospin limit as well. Nevertheless, in order to simplify the formulae, we will consider in the following the case with equal nucleon masses. We also introduce the $`2\times 2`$ matrix:

$`n(p)`$ $`=`$ $`\left(\begin{array}{cc}\theta (k_F^{(p)}p)& 0\\ 0& \theta (k_F^{(n)}p)\end{array}\right)`$ (62)
$``$ $`\left(\begin{array}{cc}n(p)_1& 0\\ 0& n(p)_2\end{array}\right)`$ (65)
$`=`$ $`\widehat{n}(p)I_2+\overline{n}(p)\tau _3`$ (66)
with $`I_2`$ the $`2\times 2`$ unity matrix, $`\tau _3`$ the usual Pauli matrix $`diag(1,1)`$, $`\widehat{n}(p)=(n(p)_1+n(p)_2)/2`$ and $`\overline{n}(p)=(n(p)_1n(p)_2)/2`$. Then eq.(56) can be rewritten as:

$`e^{i𝒵[v,a,s,p]}={\displaystyle }[dU]\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }+{\displaystyle }d𝕩{\displaystyle }d𝕪{\displaystyle }{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}`$ (73)
$`e^{i𝕡(𝕩𝕪)}\text{Tr}\left(n(p)\mathrm{log}[]|_{(𝕩,𝕪)}(\overline{)𝑝}+m)\right)\}=`$
$`={\displaystyle }[dU]\mathrm{exp}\{i{\displaystyle }dx_{\pi \pi }i{\displaystyle }{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}}{\displaystyle }dxdye^{ip(xy)}`$
$`\text{Tr}\left(A[I_4D_0^1A]^1|_{(x,y)}(\overline{)𝑝}+m_N)n(p)\right)+`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{d𝕡}{(2\pi )^32E(p)}\frac{d𝕢}{(2\pi )^32E(q)}𝑑x𝑑x^{}𝑑y𝑑y^{}e^{ip(xy)}}`$
$`e^{iq(x^{}y^{})}\text{Tr}(A[I_4D_0^1A]^1|_{(x,x^{})}(\overline{)𝑞}+m_N)n(q)`$
$`A[I_4D_0^1A]^1|_{(y^{},y)}(\overline{)𝑝}+m_N)n(p))+\mathrm{}\}.`$
Comparing this equation with eq.(51), the only difference in the Feynman rules is the inclusion of an isospin matrix $`n(p)`$ associated to every Fermi-sea insertion.

The case of a polarized nuclear matter can be treated in the same way just by doing in eq.(56) the replacement:

$`{\displaystyle \underset{r=1}{\overset{2}{}}}{\displaystyle ^{k_F}}𝑑𝕡u_{r\alpha }^{}(𝕡)\mathrm{}u_{r\alpha }(𝕡){\displaystyle ^{k_F}}𝑑𝕡u_{r_1\alpha }^{}(𝕡)\mathrm{}u_{r_1\alpha }(𝕡)`$ (74)
$`+{\displaystyle ^{k_F}}𝑑𝕡u_{r_2\alpha }^{}(𝕡)\mathrm{}u_{r_2\alpha }(𝕡)`$ (75)
since there are two spin states per momentum state.

6. The present formalism is applied in ref. to evaluating several quantities relevant for low energy QCD in the nuclear medium. There we will also address in detail the issue of the chiral counting from eq.(73) in the nuclear medium and the limitations of a plain perturbative treatment of CHPT at finite density. Just to mention that instead of a relativistic treatment of the baryons we could also have considered in the same way the non-relativistic case. This makes more straightforward the evaluation of baryon loops  although one has also to take into account recent developments in the field of effective field theories with propagation of relativistic heavy particles . In any case the main issue still to be addressed in the medium, as discussed in refs. , is to deal with the problem of the large S-wave scattering lengths in the nucleon-nucleon scattering which introduces a new extra scale of only $``$10 MeV alredy in the vaccum case, where consistent power counting schemes have been developed. Some interesting findings in this direction, requiring further consideration, can be found in ref. for a theory without pions.

To conclude, we have derived the $`SU(2)`$ chiral Lagrangian with extenal sources in the presence of non-zero nuclear density by explicitely working out in quantum field theory the in-medium contributions from the $`\pi N`$ chiral Lagrangian, eq.(51)(symmetric nuclear matter), eq.(73) (asymmetric nuclear matter) and eq.(74) (asymmetric and unpolarized nuclear matter). Then the perturbative theory is developed and the corresponding Feynman rules are given. Contrarily to standard many-body techniques, the rules and diagrams here derived for the general relativistc case are analogous to the usual ones from vacuum quantum field theory without modification of the baryon propagators establishing a neat separation between in-medium and vacuum contributions, as indicated in figs.1 and 2. Applications of this many-particle relativistic quantum field theory formalism to actual calculations can be found in ref..

I would like to thank Andreas Wirzba and Ulf-G. Meißner for useful discussions and for a critical reading of the manuscript. This work was supported in part by funds from DGICYT under contract PB96-0753 and from the EU TMR network Eurodaphne, contract no. ERBFMRX-CT98-0169.


format: text

text:
# A Measurement of H₀ from the Sunyaev-Zeldovich Effect

## 1 Introduction

The Sunyaev-Zeldovich effect (SZE) is a spectral distortion in the Cosmic Microwave Background (CMB) due to inverse Compton scattering of the CMB photons off hot electrons; at radio frequencies this distortion is manifested as a fractional decrement in intensity of order $`10^4`$. For over two decades it has been known (Silk & White, 1978; Cavaliere et al., 1979) that the combination of X-ray and SZE observations of rich galaxy clusters, under the assumption of spherical symmetry, yields a direct measurement of the cosmic distance scale. Only in the past few years have reliable and accurate applications of this method become possible (e.g. Birkinshaw et al., 1991; Herbig et al., 1995; Carlstrom et al., 1996; Grainge et al., 1996; Holzapfel et al., 1997). Due to the assumption of spherical symmetry, selection biases have been a great concern. To address this Myers et al. (1997) defined an X-ray flux-limited sample of 11 $`z<0.1`$ clusters and began a campaign to measure distances to these clusters. This yields an orientation-unbiased sample, and departures from spherical symmetry in individual clusters will average out in the determination of $`\mathrm{H}_0`$ from the sample as a whole. Myers et al. obtain a Hubble Constant of $`54\pm 14\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ from observations of four clusters. The accuracy of this result is limited by a 7% radio calibration uncertainty and estimated 15-30% X-ray model uncertainties for each cluster.

In this Letter we present an improved measurement of $`\mathrm{H}_0`$ based on the Herbig et al. and Myers et al. results on Coma, A478, A2142, and A2256, plus observations of three new clusters (A399, A401, A1651) from their complete sample. The results we present incorporate more accurate X-ray models (Mason & Myers, 2000, hereafter,MM2000) and a more accurate radio calibration (Mason et al., 1999). We also calculate the effect of intrinsic CMB anisotropies on our result and find them to be a limiting factor for a sample of this size. Throughout this paper we use $`\mathrm{H}_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and consider two cosmologies: SCDM ($`\mathrm{\Omega }_M=1,\mathrm{\Omega }_\mathrm{\Lambda }=0`$), and $`\mathrm{\Lambda }\mathrm{CDM}`$ ($`\mathrm{\Omega }_M=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$). Unless otherwise stated, we assume the SCDM cosmology.

## 2 The Sunyaev-Zeldovich Effect, Cluster Sample, and X-Ray Models

### 2.1 The Sunyaev-Zeldovich Effect

The observed fractional change in antenna temperature induced by a cloud of electrons with a beam-averaged Compton $`y`$parameter $`y_{obs}`$ is (Sunyaev & Zeldovich, 1980)

$$\frac{\mathrm{\Delta }T_{obs}}{T_{cmb}}=\frac{x^2e^x}{(e^x1)^2}\left[x\mathrm{coth}(x/2)4\right]\chi _{rel}^1\times y_{obs}.$$
(1)
Here $`x=h\nu /kT_{cmb}`$, and $`\chi _{rel}^1`$ is a relativistic correction factor which we compute using the analytic expression of Sazonov & Sunyaev (1998). We assume zero peculiar velocity, which introduces a $`5\%`$ uncertainty but will average out over the sample. Along a given line of sight, the Compton $`y`$parameter is

$$y=\frac{kT_e}{m_ec^2}\tau $$
(2)
where $`\tau `$ is the inverse Compton optical depth. At $`32.0\mathrm{GHz}(x=0.563)`$, Eq. 1 takes the form

$$\frac{\mathrm{\Delta }T_{obs}}{T_{cmb}}=1.897\times y_{obs}\chi _{rel}^1.$$
(3)
We follow the procedure of Myers et al. in correcting the observed Compton $`y`$parameters for relativistic effects rather than the models.

Given a model for the cluster density and temperature profiles from X-ray observations, we predict the SZE decrement in terms of the Compton $`y`$parameter averaged over the telescope beam and switching patterns:

$$y_{pred}=\frac{1}{\mathrm{\Omega }_{beam}}𝑑\mathrm{\Omega }y_{model}(\widehat{\mathrm{\Omega }})R_{N,sw}(\widehat{\mathrm{\Omega }}).$$
(4)
Here $`\mathrm{\Omega }_{beam}`$ is the solid angle of the telescope beam, $`y_{model}`$ is the predicted $`y`$ along a given line of sight, and $`R_{N,sw}(\widehat{\mathrm{\Omega }})`$ is the normalized beam response, including the effects of switching. It is easily shown that $`y_{pred}`$, as determined from X-ray observations and Eq. 4, is proportional to $`h^{1/2}`$, so that if $`y_{obs}=q`$ and $`y_{pred}=ph^{1/2}`$, the Hubble constant is given by

$$h=\left(\frac{p}{q}\right)^2.$$
(5)
Note that the $`\tau _{sw}`$ values presented in MM2000 are simply the model inverse Compton optical depths convolved with the telescope beam switching pattern as per Eq. 4. For an isothermal cluster, $`y_{pred}=\tau _{sw}kT_e/m_ec^2`$.

### 2.2 Cluster Sample

In MM2000 we define a larger, X-ray flux-limited cluster sample which expands upon that presented by Myers et al. . This is a 90% volume-complete sample selected from the XBAC catalog (Ebeling et al., 1996) with $`z<0.1`$ and $`F_X>1.0\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. At $`z=0.1`$, the volume-completeness criterion corresponds to $`L_X>1.13\times 10^{44}h^2\mathrm{erg}\mathrm{s}^1`$ (0.1 - 2.4 keV). The resulting set of 31 clusters contains the Myers et al. sample. While the 7 clusters for which we present measurements here are all members of the smaller Myers et al. sample, this work is part of an ongoing project to survey distances to the objects in our expanded sample. MM2000 discusses in detail the completeness of this XBAC-derived catalog.

### 2.3 X-Ray Models

MM2000 also presents X-ray models for the 22 clusters in our sample which had public ROSAT PSPC data as of May 1999. The primary focus of this analysis was to quantify the uncertainty in the beam-convolved inverse Compton optical depth, $`\tau _{sw}`$, a quantity which we find is robustly constrained by the ROSAT data. We adopt the Markevitch et al. (1998) measurements of the cluster gas temperature $`T_e`$, which— except for the Coma cluster— we assume to be isothermal. For the Coma cluster, following Hughes et al. (1988), we adopt a hybrid model which is isothermal inside a radius of $`500h^1\mathrm{kpc}`$, outside of which the temperature follows the density profile with an adiabatic index of $`\gamma =1.5`$. For A478 and A2142, we use the fits which have been corrected for the spectral bias induced by the cooling flow emission. Table 1 summarizes our cluster models and the resulting beam-averaged Compton $`y`$values, $`y_{pred}`$. In § 4 we discuss the impact on $`\mathrm{H}_0`$ of our choices for the cluster models.

## 3 SZE Observations

In this section we report on observations of 7 clusters from our sample. Two of these, A399 and A401, are in close proximity. The effects of this are included in our model predictions (see § 4).

### 3.1 New Observations

For the observations reported here, the OVRO 5.5-meter telescope (5-m) was outfitted with a HEMT receiver having a center frequency of 32 GHz and a bandwidth of 6.5 GHz. The HEMT input is Dicke switched every millisecond between two ambient-temperature corrugated feeds which give rise to two 7’.35 FWHM Gaussian beams \[each having a main-beam volume of $`\mathrm{\Omega }_{beam}=(5.21\pm 0.03)\times 10^6\mathrm{Sr}`$\] separated by 22’.16 in azimuth. To remove systematic effects due primarily to atmosphere and ground, we employ a triple-differencing technique. Two levels of differencing are provided by the Dicke switching between the target (ON) position and a reference (REF) position, and and by nodding the telescope in azimuth at a rate of $`0.1\mathrm{Hz}`$. The third level is provided by observing blank leading (LEAD) and trailing (TRAIL) fields over the same range of azimuth and elevation as the MAIN field, so that ground-based signals are cancelled in the difference. The cluster signal is given by the MAIN field signal minus the average of the LEAD and TRAIL signals, or MLT. The LEAD-TRAIL field difference (LTD) provides a diagnosis of possible residual signals due, e.g., to unsubtracted ground emission or intrinsic CMB fluctuations.

Flux density calibration is accomplished with reference to Cas A using an epoch 1998 flux of $`194\pm 4\mathrm{Jy}`$ (Mason et al., 1999), with the secular variation of Cas A modeled as per Baars et al. (1977). Data are edited by automatic filtering algorithms that use the scatter of the data and their standard deviations as rejection criteria. More discussion of our observing strategy and data filters is given in Readhead et al. (1989) and Myers et al. (1997); more details on the instrument are in Leitch et al. (2000).

#### 3.1.1 Discrete Source Removal

The brightest contaminating discrete sources in A399 and A1651, selected from 87GB (Gregory & Condon, 1991), were monitored at 18.5 GHz with the OVRO 40-meter telescope in November and December of 1996, concurrent with cluster SZE observations. In the Fall of 1997 we used the OVRO 40-m telescope to survey all NVSS (Condon et al., 1998) sources in these fields and those for A401 which, assuming a flat spectrum, would give a peak signal at 32 GHz $`>16\mu \mathrm{K}`$. We also observed all sources within $`12^{}`$ of any of the field centers and having 1.4 GHz flux densities $`>50\mathrm{mJy}`$. Sources which showed indications of variability were extrapolated to 32 GHz assuming a flat spectrum; otherwise the two-point spectral index between 1.4 and 32 GHz was used. The overall corrections in $`\mu \mathrm{K}`$ are given for each cluster in § 3.1.2.

#### 3.1.2 A399, A401, and A1651

Observations of A399 were taken from October 1996 through March 1997. For these observations the LEAD and TRAIL fields were separated from the MAIN by 26 minutes of Right Ascenscion. After statistical filters and weighting of the data, we acquired 40 hours of total integration time including LEAD and TRAIL fields, and time spent on the REF beams. To account for source contamination $`24\pm 5\mu \mathrm{K}`$ was subtracted from the observed signal. We then determine a decrement of $`\mathrm{\Delta }T_{obs}=164\pm 21\mu \mathrm{K}`$ (MLT). The LEAD-TRAIL difference is $`\mathrm{\Delta }T_{\mathrm{LTD}}=15\pm 21\mu \mathrm{K}`$. The observed $`\mathrm{\Delta }T`$ corresponds to $`y_{obs}=(3.28\pm 0.42)\times 10^5`$, including the relativistic correction $`\chi _{rel}=1.027`$.

A1651 was also observed during this period, with LEAD and TRAIL field separations of $`20^{\mathrm{min}}\mathrm{\hspace{0.17em}30}^{\mathrm{sec}}`$ from the MAIN field; the effective total integration time was 25 hours. We subtract $`5\pm 1\mu \mathrm{K}`$ from the observed decrement to correct for discrete sources. These are mostly in the REF beam of the LEAD field. With the source correction we determine a decrement $`\mathrm{\Delta }T_{obs}=247\pm 30\mu \mathrm{K}`$ (MLT). The LTD is $`\mathrm{\Delta }T_{\mathrm{LTD}}=92\pm 36\mu \mathrm{K}`$. The Compton-$`y`$ parameter is $`y_{obs}=(4.88\pm 0.59)\times 10^5`$ including the relativistic correction $`\chi _{rel}=1.023`$.

A401 was observed from October 1997 to March 1998, giving an effective total integration time of 43 hours. The LEAD and TRAIL field separations were $`16^{\mathrm{min}}\mathrm{\hspace{0.17em}36}^{\mathrm{sec}}`$. Discrete sources contribute $`16\pm 11\mu \mathrm{K}`$. The source-corrected decrement is $`\mathrm{\Delta }T_{obs}=338\pm 20\mu \mathrm{K}`$ (MLT); the LTD is $`\mathrm{\Delta }T_{\mathrm{LTD}}=133\pm 24\mu \mathrm{K}`$. This gives (with $`\chi _{rel}=1.031`$) $`y_{obs}=(6.33\pm 0.43)\times 10^5`$.

The statistically significant LTD’s for A401 and A1651 are discussed in § 4.

### 3.2  Coma Cluster

Herbig et al. (1995) used the OVRO 5-m telescope to determine the SZE decrement towards the Coma cluster (Abell 1656) during the observing seasons of 1992 and 1993, and found a Rayleigh-Jeans temperature decrement of $`\mathrm{\Delta }T_{obs}=302\pm 48\mu \mathrm{K}`$. These data were calibrated relative to DR21 assuming $`S_{DR21}=18.24\pm 0.55\mathrm{Jy}`$, and referenced to the telescope main beam assuming $`\mathrm{\Omega }_{Beam}=(5.16\pm 0.15)\times 10^6\mathrm{Sr}`$. From the Mason et al. flux density scale and $`\chi _{rel}=1.035`$, we determine a beam-averaged Compton $`y`$parameter of $`y_{sw}=(6.38\pm 1.01)\times 10^5`$.

### 3.3  A478, A2142, and A2256

A478, A2142, and A2256 were observed by Myers et al. (1997) from July 1993 to March 1994 resulting in decrements of $`375\pm 28,437\pm 25,`$ and $`243\pm 29\mu \mathrm{K}`$ respectively; these temperature are referred to the main beam, and have been corrected for point-source contamination. These data were calibrated using a brightness temperature of $`T_J=144\pm 8`$ K for Jupiter. After including the correction for the main beam area \[Myers et al. use $`(5.12\pm 0.14)\times 10^6\mathrm{Sr}`$\], we find an overall correction of $`f=1.037`$ for the Myers et al. data. The stated decrements yield switched, beam-averaged Compton-$`y`$ parameters of $`(7.25\pm 0.54)\times 10^5,(8.44\pm 0.48)\times 10^5,`$ and $`(4.70\pm 0.56)\times 10^5`$, respectively. We apply relativistic corrections ($`\chi _{rel}=1.033`$ , 1.037 , and 1.026) together with our calibration correction, and find Compton $`y`$values of $`(7.77\pm 0.58)\times 10^5`$ , $`(9.10\pm 0.52)\times 10^5`$ , and $`(5.00\pm 0.60)\times 10^5`$.

## 4 Interpretation

For SZE observations at our frequencies and angular scales intrinsic CMB anisotropies are a significant source of uncertainty. We use the RING5M measurements (Leitch et al., 2000) to determine the impact of this on our $`\mathrm{H}_0`$ results. In measurements taken with the OVRO 5-m telescope at 32 GHz, Leitch et al. determined an RMS temperature difference on $`22^{}`$ scales of $`\delta T_{22^{}}=79_{10}^{+11}\mu \mathrm{K}`$. In order to determine the degree to which parallactic angle averaging during a track on a cluster reduces this, we generated $`10^3`$ realizations of $`4\mathrm{deg}^2`$ patches of sky from a $`\mathrm{\Lambda }`$CDM power spectrum. Each realization was convolved with the 5-m main beam and the switching pattern characteristic of a typical SZE observation. We find (including the LEAD/TRAIL differencing) a residual CMB signal of $`81.6\%\times \delta T_{22^{}}`$, or $`64\mu \mathrm{K}`$ ($`\sigma _y=1.24\times 10^5`$) assuming the RING5M power level. The RMS of LTD predicted by this analysis is $`90\mu \mathrm{K}`$, close to the observed RMS of $`79\mu \mathrm{K}`$ in our 7 clusters. On this basis the non-zero LTD’s seen in A1651 and A401 are not unexpected.

While this noise level is not small compared to the signal in fainter clusters such as A399, we leave these clusters in our sample to avoid an orientation bias and perform a Maximum Likelihood analysis to determine $`\mathrm{H}_0`$. Assume that we have $`N`$ observed Compton-$`y`$ paramaters $`q_i`$, and $`N`$ predicted $`y`$-parameters (as per Eq. 5) $`p_i`$ determined from X-ray observations. The SZE observations are corrupted with random Gaussian noise $`\sigma _{q,i}`$ (including thermal noise and point-source subtraction uncertainties) and the predictions are corrupted with random Gaussian noise $`\sigma _{p,i}`$; the RMS residual CMB signal is $`\sigma _{cmb}=64\mu \mathrm{K}`$. In addition there are $`N`$ “true” predictions $`\widehat{p_i}`$ which are the predictions which we would have in the absence of errors in the X-ray models. Then we can determine $`h`$ by minimizing

$$\chi ^2=\underset{i=1}{\overset{N}{}}\frac{(q_ih^{1/2}\widehat{p_i})^2}{\sigma _{q,i}^2+\sigma _{cmb}^2}+\frac{(p_i\widehat{p_i})^2}{\sigma _{p,i}^2},$$
(6)
where $`\chi ^22lnL`$. The unconstrained $`\widehat{p_i}`$ must also be minimized in this process, but these can be projected out analytically. Setting the derivatives of $`\chi ^2`$ with respect to the $`\widehat{p_i}`$ equal to zero, we find

$$\widehat{p_i}=\frac{q_i+f_ih^{1/2}p_i}{h^{1/2}+f_ih^{1/2}}$$
(7)
with $`f_i=(\sigma _{q,i}^2+\sigma _{cmb}^2)/\sigma _{p,i}^2`$. This leaves a 1-dimensional parameter space. The value of $`h`$ which minimizes Eq. 6 is the most probable value; two-sided 68% confidence intervals are determined about this value by integrating the likelihood function $`L`$. We impose the prior $`0.01<h<5`$; our upper bound on $`\mathrm{H}_0`$ for A399 is slightly affected by our choice of prior, but other results are not. Monte Carlo tests with the estimator of Eq. 6 show that within the range of power levels allowed by RING5M, our result is not affected by the uncertainty in $`\sigma _{cmb}`$.

Table 2 shows the measured and predicted values of the Compton $`y`$parameter for our 7 clusters along with the Maximum-Likelihood values for and 68% confidence limits on $`\mathrm{H}_0`$. A399 and A401 are in close spatial proximity with the result that the observed decrement on either is reduced somewhat due to signal from the other in the reference beam; this is accounted for in the stated values of $`y_{pred}`$. Averaged over a track, we find a decrease in the predicted $`y`$ for A399 of $`(0.15\pm 0.05)\times 10^5h^{1/2}`$ and for A401 of $`(0.13\pm 0.05)\times 10^5h^{1/2}`$. For an SCDM cosmology, the sample average is $`64_{11}^{+14}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. The scatter of the $`\mathrm{H}_0`$ values gives an error in the mean of only $`8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ , which is less than that given by our Maximum Likelihood analysis. The minimum of $`\chi ^2`$ corresponds to $`\chi _\nu ^2=0.32`$ for $`\nu =6`$; there is a $`5\%`$ chance of obtaining a $`\chi _\nu ^2`$ value this low by chance for six degrees of freedom. We therefore think it likely that the scatter in these 7 clusters is fortuitously small. Since our observing and modelling errors have been well quantified, the Maximum Likelihood method gives a more reliable estimate of the uncertainty than the data scatter for a small sample such as ours. Our result is not sensitive to the model we choose for the cluster temperature profile: adopting hybrid models for all 5 clusters reduces the sample average by only 4%. Although more extreme models for the temperature profiles would have a larger effect on our result, such models are not motivated by current X-ray analyses (Irwin, Bregman & Evrard, 1999; White, 2000). For a $`\mathrm{\Lambda }\mathrm{CDM}`$ cosmology, the sample average is $`66_{11}^{+14}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. The calibration uncertainties are $`3\%`$ (radio) and $`8\%`$ (X-ray), and we estimate a $`10\%`$ uncertainty in the SZE predictions due to the possibility of substructure and non-isothermality in the ICM. Altogether we have a systematic error budget of $`\pm 14_{sys}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, or $`\pm 15_{sys}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ for $`\mathrm{\Lambda }\mathrm{CDM}`$.

The increase in $`\mathrm{H}_0`$ relative to the Myers et al. result is due primarily to differences in the updated X-ray models. As a check on our ROSAT-derived density models, we have compared the baryonic masses that we obtain inside $`500h^{}\mathrm{kpc}`$ to those obtained by Mohr et al. (1999) for 20 of the 22 clusters in common to our analyses: we find that the mean mass ratio is equal to unity at better than $`1\%`$ accuracy. See MM2000 for more discussion of this analysis. Also the $`T_e`$ values we have used tend to be higher than the Myers et al. values due to the correction for the spectral bias of the cooling flows. Using the Myers et al. temperatures increases the fractional scatter in the $`\mathrm{H}_0`$ values by over a factor of two and decreases the sample mean to $`55\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. The fact that the uncertainty in our result is comparable to that of Myers et al. — in spite of improved X-ray models, better calibration, and 3 more clusters— is due to our inclusion of intrinsic anisotropy in the analysis.

## 5 Conclusions

We have presented a determination of $`\mathrm{H}_0`$ resulting from measuring the SZE in 5 clusters from an unbiased, low-redshift sample. For an SCDM cosmology we find $`\mathrm{H}_0=64_{11}^{+14}\pm 14_{sys}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, while for a $`\mathrm{\Lambda }\mathrm{CDM}`$ cosmology we obtain $`\mathrm{H}_0=66_{11}^{+14}\pm 15_{sys}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. This result is in good agreement with other recent measurements. Madore et al. (1999), using Cepheid variables as distance indicators, find $`72\pm 5\pm 7_{sys}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Gravitational lens time delays give results consistent with these (e.g. Fassnacht et al., 1999).

Our value of $`\mathrm{H}_0`$ implies an age for the universe ranging from $`(10.2\pm 3.2)\times 10^9`$ years for an SCDM universe to $`(14.2\pm 4.5)\times 10^9`$ years for $`\mathrm{\Lambda }\mathrm{CDM}`$. These are both consistent with recent age determinations from main-sequence fitting, which give ages of $`(>12\pm 1)\times 10^9`$ years (Chaboyer et al., 1998). In the X-ray cluster analysis of MM2000 we find a mean ICM mass fraction $`f_{ICM}M_{bary}/M_{tot}=(7.02\pm 0.28)h^{3/2}\times 10^2`$ (within $`R_{500}1h^1\mathrm{Mpc}`$) in a sample of 22 nearbly clusters. We combine this result with the Big Bang Nucleosynthesis constraint $`\mathrm{\Omega }_Bh^2=0.019\pm 0.001`$ (Burles & Tytler, 1998) and our value for $`\mathrm{H}_0`$ to find a total matter density $`\mathrm{\Omega }_M=0.32\pm 0.05`$, arguing against SCDM cosmologies.

We are grateful to Russ Keeney for his many years of hard work on the 5-m and 40-m telescopes at OVRO, and to Erik Leitch for writing the software that was used to analyze the OVRO data. We thank Jonathan Sievers for assistance with the intrinsic CMB simulations and the 1997/1998 OVRO observations. We received suport from NSF grants AST 91-19847 and AST 94-19279. For part of the duration of this work BSM was supported by the Zacheus Daniels fund at the University of Pennsylvania. STM was supported by an Alfred R. Sloan Fellowship at the University of Pennsylvania. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.


text:
# Electron localization in the insulating state: Application to crystalline semiconductors

## I Introduction

A nonmetal is distinguished from a metal by its vanishing conductivity at low temperature and low frequency: we use here the term “insulator” to include any nonmetal, like the semiconducting materials which are the case studies actually addressed in this work.

Within classical physics, the qualitative difference between an insulator and a metal is attributed to the nature of the electronic charge: either “bound” (Lorentz model for insulators) or “free” (Drude model for metals). In other words, electrons are localized in insulators and delocalized in metals. In a milestone paper published in 1964, W. Kohn characterized the insulating state of matter in a way which is reminiscent of the classical picture: he gave evidence that the main feature determining the insulating behavior of matter is electron localization in the ground–state wavefunction. Although this work mainly addressed correlated many–electron systems, its message is very relevant even for materials where an independent–electron description is quite adequate, as the semiconductor crystals studied here. Recently a novel measure of electron localization—different from Kohn’s one—was proposed by Resta and Sorella, hereafter cited as RS. Their approach is deeply rooted into the modern theory of polarization.

Metals and insulators reveal their qualitative difference when static dielectric polarization is addressed. Suppose we expose a finite macroscopic sample to an electric field, say inserting it in a charged capacitor. In metals polarization is trivial: universal, material–independent, due to surface phenomena only (screening by free carriers). Therefore polarization in metals is not a bulk phenomenon. The opposite is true for insulators: macroscopic polarization is a nontrivial, material–dependent, bulk phenomenon. We can therefore phenomenologically characterize an insulator, in very general terms, as a material whose ground wavefunction sustains a nonzero bulk macroscopic polarization whenever the electronic Hamiltonian is non centrosymmetric. If the Hamiltonian is instead centrosymmetric, the polarization vanishes but remains a well defined bulk property, at variance with the metallic case. The phenomenological link between macroscopic polarization and insulating behavior was first pointed out and exploited—taking advantage of the modern theory of polarization —by RS in 1999. This approach is based on Berry phases and related concepts. Even the RS paper, like Kohn’s 1964 one, mostly concerns correlated systems. Furthermore, in order to keep the presentation simple and concise, most results are explicitly shown in one dimension, while the $`d`$–dimensional formulation is only sketched in the final paragraphs of RS. In the present paper we provide more details on how the RS theory of localization works in three dimensions, specializing to a system of noninteracting electrons, like the band insulators chosen as case studies here.

Some other important papers must be mentioned at this point. In 1997 Marzari and Vanderbilt, hereafter cited as MV, while not addressing metals at all (and hence their difference from insulators), establish nonetheless some results which are relevant to the present viewpoint. In a very recent comprehensive paper Souza, Wilkens, and Martin—hereafter cited as SWM—generalize and extend in various ways the main finding of RS: we adopt here some of their notations. Finally, after this work was completed, we became aware of Ref. , whose conclusions bear some implications for our results shown in Sect. VI.

The paper is organized as follows. In Sect. II we define the basic ingredients providing both polarization and localization, namely the expectation values of the many–body phase operators $`z_N^{(\alpha )}`$ for the three Cartesian coordinates, Eq. (4). In Sect. III, following RS, we show how the modulus of $`z_N^{(\alpha )}`$ defines a very meaningful quantity, the localization tensor, for which we adopt the SWM notations: such tensor is finite in insulators and diverges in metals. In Sect. IV we discuss the main properties of the localization tensor, and in Sect. V we present first–principle calculations for several elemental and binary semiconductors: the main trends are analyzed. In Sect. VI we calculate orbitals which are optimally localized in a given direction, and whose average quadratic spread coincides with the localization tensor. We also heuristically check the exponential localization of these orbitals, which we call “hermaphrodite orbitals”. In Sect. VII we draw our main conclusions. In the Appendix we consider a molecule or a cluster an we discuss our localization tensor therein, showing its relationship to some results of Boys localization theory, well known in quantum chemistry.

## II Many–body phase operators

We are addressing here, as it is done by MV, a crystalline system of independent–electrons, having in mind a Kohn–Sham scheme. The properties of interest, namely, macroscopic polarization and electron localization, are not properties of the individual KS orbitals: instead, they are global properties of the occupied KS manifold. As shown in Refs. , it proves formally convenient to deal with a many–body wavefunction $`\mathrm{\Psi }`$, obtained as a Slater determinant of occupied orbitals. This determinant is uniquely determined by the manifold of the occupied orbitals and is invariant by unitary transformation of these orbitals among themselves: for instance, in insulating crystals, an important transformation of this class converts the occupied Bloch orbitals into Wannier functions. Quantities which can be expressed solely in terms of $`\mathrm{\Psi }`$ are invariant in form under such transformations.

Throughout this work—with the exception of the Appendix—we adopt periodic Born–von–Kàrmàn boundary conditions (BvK) on a large cell, multiple of the crystalline elementary cell. The quantities of interest are intensive and have a well defined thermodynamic limit, while the wavefunction itself becomes an ill–defined mathematical object in that limit.

For the sake of simplicity, we assume a simple cubic cell of side $`a`$ and a large BvK cell of side $`L=Ma`$. More general structures can be dealt with using scaling, similarly to what shown e.g. in Ref.  or in SWM. The thermodynamic limit corresponds to $`M\mathrm{}`$, while practical calculations are performed at finite, and possibly large, $`M`$ values. The spinorbitals $`\psi `$ (spin–up) and $`\overline{\psi }`$ (spin–down) may be chosen of the Bloch form. In the finite system there are $`M^3`$ allowed Bloch vectors $`𝐪_s`$, arranged on a regular mesh in the unit reciprocal cell, where $`s(s_1,s_2,s_3)`$ and

$$𝐪_s=\frac{2\pi }{Ma}(s_1,s_2,s_3),s_\alpha =0,1,\mathrm{},M1.$$
(1)
We adopt a plane–wave–like normalization for the Bloch orbitals:

$$\psi _{n𝐪_s}|\psi _{n^{}𝐪_s^{}}=\frac{1}{L^3}_{\text{BvK cell}}𝑑𝐫\psi _{n𝐪_s}^{}(𝐫)\psi _{n^{}𝐪_s^{}}(𝐫)=\delta _{nn^{}}\delta _{ss^{}},$$
(2)
If the system is insulating with $`n_b`$ doubly occupied bands, there are $`N=2n_bM^3`$ independent spinorbitals, out of which we write a single–determinant many–body wavefunction for $`N`$ electrons:

$$\mathrm{\Psi }=𝖠\underset{n,s}{}\frac{1}{L^3}\psi _{n𝐪_s}\overline{\psi }_{n𝐪_s},$$
(3)
where the product runs over all occupied bands and all mesh points, $`𝖠`$ is the antisymmetrizer operator, and the factor ensures that the $`N`$–body wavefunction is normalized to one on the hypercube of side $`L`$. If, instead, the system is metallic, then the many–body wavefunction $`\mathrm{\Psi }`$ can still be written in the form of Eq. (3), but where not all the Bloch vectors of a given band are included in the product.

According to Refs. , the key quantities to deal with both macroscopic polarization and electron localization are expectation values of “many–body phase operators”. For a three–dimensional system there are three such operators, one for each Cartesian direction. We indicate as $`z_N^{(\alpha )}`$, where $`\alpha `$ is a Cartesian index, their ground–state expectation values:

$$z_N^{(x)}=\mathrm{\Psi }|\mathrm{e}^{i\frac{2\pi }{L}_{i=1}^Nx_i}|\mathrm{\Psi },$$
(4)
and analogously for $`y`$ and $`z`$ directions. This remarkably compact expression is very general and applies as it stands even to correlated and/or disordered systems: here we specialize to a crystalline system of independent electrons, whose wavefunction $`\mathrm{\Psi }`$ assumes the form of Eq. (3), where the product indices have to be differently specified in the insulating and metallic cases.

We may conveniently recast $`z_N^{(x)}`$ as an overlap

$$z_N^{(x)}=\mathrm{\Psi }|\stackrel{~}{\mathrm{\Psi }},$$
(5)
where $`\stackrel{~}{\mathrm{\Psi }}`$ is the Slater determinant of a different set of Bloch spinorbitals:

$$\stackrel{~}{\psi }_{n𝐪_s}(𝐫)=\mathrm{e}^{i\frac{2\pi }{Ma}x}\psi _{n𝐪_s}(𝐫),$$
(6)
and analogously for the bar (spin–down) ones. According to a well known theorem, the overlap between two single–determinant wavefunctions is equal to the determinant of the $`N\times N`$ overlap matrix built out of the occupied spinorbitals. Since the overlaps between different–spin spinorbitals vanish, and those between equal–spin ones are identical in pairs, we can write

$$z_N^{(x)}=(\text{det}𝒮)^2,$$
(7)
where $`𝒮`$ is the overlap matrix between spatial orbitals, having size $`N/2=n_bM^3`$. Its elements are:

$`𝒮_{n𝐪_s,n^{}𝐪_s^{}}={\displaystyle \frac{1}{L^3}}{\displaystyle _{\text{BvK cell}}}𝑑𝐫\psi _{n𝐪_s}^{}(𝐫)\stackrel{~}{\psi }_{n^{}𝐪_s^{}}(𝐫)`$ (8)
$`={\displaystyle \frac{1}{L^3}}{\displaystyle _{\text{BvK cell}}}𝑑𝐫u_{n𝐪_s}^{}(𝐫)u_{n^{}𝐪_s^{}}(𝐫)\mathrm{e}^{i(\frac{2\pi }{Ma}x+𝐪_s^{}𝐫𝐪_s𝐫)},`$ (9)
where the $`u`$’s are the periodic functions in the Bloch orbitals.

The matrix $`𝒮`$ is very sparse: in fact, given the geometry of the $`𝐪_s`$’s on the regular reciprocal mesh (see Eq. (1), the overlap integrals in Eq. (9) are nonvanishing only if $`s_1=s_1^{}+1`$, $`s_2=s_2^{}`$, and $`s_3=s_3^{}`$. We express the nonvanishing elements in terms of a small overlap matrix $`S`$, of size $`n_b\times n_b`$:

$$S_{nn^{}}(𝐪,𝐪^{})=u_{n𝐪}|u_{n^{}𝐪^{}}=\frac{1}{a^3}_{\text{cell}}𝑑𝐫u_{n𝐪}^{}(𝐫)u_{n^{}𝐪^{}}(𝐫).$$
(10)
Owing to the sparseness of $`𝒮`$, its determinant factors into products of determinants of small matrices $`S`$.

In the insulating case we use the wavefunction of Eq. (3), where all the Bloch vectors of a given band are occupied: the factorization is then

$$z_N^{(x)1/2}=\text{det}𝒮=\underset{s}{}\text{det}S(𝐪_{s_1+1,s_2,s_3},𝐪_{s_1,s_2,s_3}).$$
(11)
We get a more compact notation upon defining

$$\mathrm{\Delta }𝐪^{(x)}=\frac{2\pi }{L}(1,0,0)=\frac{2\pi }{Ma}(1,0,0),$$
(12)
which is the vector connecting nearest neighbor q points in the $`x`$ direction. We have then:

$$z_N^{(x)1/2}=\underset{s}{}\text{det}S(𝐪_s+\mathrm{\Delta }𝐪^{(x)},𝐪_s);$$
(13)
In the metallic case, instead, the $`z_N^{(\alpha )}`$’s are identically zero. This is easily understood by looking at the simple case with only one band. Suppressing the band index the small overlap matrix becomes a c–number $`S(𝐪,𝐪^{})`$, and Eq. (11) becomes a product of c–numbers, with no determinant to evaluate. In an insulator this product runs over the whole $`𝐪_s`$ mesh, and all factors are nonzero; in a metal the analogous product runs only on the $`𝐪_s`$’s within the Fermi surface. Looking at the definition of $`S(𝐪,𝐪^{})`$, Eqs. (9) and (10), it is clear that there exists at least one occupied $`𝐪_s`$, adjacent to the Fermi surface, such that $`S(𝐪_s,𝐪_s^{})`$ vanishes for all occupied $`s^{}`$. This is enough to imply that $`z_N^{(\alpha )}`$ vanishes as well.

The complex numbers $`z_N^{(\alpha )}`$ are ground–state expectation values, and do not access any spectral information. Yet they qualitatively discriminate between insulators and metals: they are in fact nonvanishing in the former materials, and vanishing in the latter ones. This shows, according to RS, that there is a qualitative difference in the organization of the electrons in their ground state. It is remarkable that, in the present case, such difference shows up already at finite $`N`$, before the thermodynamic limit is taken.

## III Electron localization

In centrosymmetric materials the expectation values $`z_N^{(\alpha )}`$ are real (provided the origin is chosen at a centrosymmetric site), while in noncentrosymmetric materials they are in general complex: their phases define then the Cartesian components of the macroscopic polarization in suitable units. In the metallic case the $`z_N^{(\alpha )}`$’s vanish and the polarization is ill defined, in agreement with the phenomenological viewpoint illustrated in Sect. I. We address electron localization using the moduli of these same $`z_N^{(\alpha )}`$’s. Following RS, electron localization is measured by a squared localization length in one dimension, and by a “localization tensor” in three dimensions. This tensor is an intensive quantity, has the dimensions of a squared length, and measures the localization of the many–electron system as a whole: in the present case, it is a global property of the occupied KS manifold. The localization tensor is finite for insulators and diverges for metals.

In the very recent SWM paper it is shown, among other things, that the RS localization tensor is related to the mean–square quantum fluctuation of the polarization: it is a second cumulant moment, which can be very elegantly extracted from a moment generating function. We adopt throughout notations inspired by SWM, and we indicate the localization tensor as $`r_\alpha r_\beta _\mathrm{c}`$, where the subscript stays for “cumulant”. For a material having cubic or tetrahedral symmetry, like the semiconductors considered in the present case studies, the localization tensor is isotropic: its only independent element is $`x^2_\mathrm{c}`$. Its expression is provided by RS, whose Eq. (18) we recast here as

$$x^2_\mathrm{c}=\frac{1}{N}\left(\frac{L}{2\pi }\right)^2\mathrm{ln}|z_N^{(x)}|^2,$$
(14)
and the thermodynamic limit is understood. For a metal $`z_N^{(x)}`$ vanishes and the localization tensor is formally infinite, even at finite $`N`$. For an insulator, whose wavefunction has the form of Eq. (3), we get from Eq. (13)

$$|z_N^{(x)}|=\underset{s}{}\text{det}S^{}(𝐪_s,𝐪_s+\mathrm{\Delta }𝐪^{(x)})S(𝐪_s,𝐪_s+\mathrm{\Delta }𝐪^{(x)});$$
(15)
$$x^2_\mathrm{c}=\left(\frac{a}{2\pi }\right)^2\frac{1}{2n_bM}\mathrm{ln}|z_N^{(x)}|^2.$$
(16)
Eqs. (15) and (16) are the typical expressions implemented in our test–case calculations discussed below. The thermodynamic limit is obtained as usual for $`M\mathrm{}`$ and takes, not surprisingly, the form of an integral performed over the reciprocal unit cell, or equivalently over the first Brillouin zone. The integral is:

$`x^2_\mathrm{c}={\displaystyle \frac{a^3}{n_b(2\pi )^3}}{\displaystyle }d𝐪({\displaystyle \underset{n}{}}{\displaystyle \frac{}{q_x}}u_{n𝐪}|{\displaystyle \frac{}{q_x}}u_{n𝐪}`$ (17)
$`{\displaystyle \underset{n,n^{}}{}}u_{n𝐪}|{\displaystyle \frac{}{q_x}}u_{n^{}𝐪}{\displaystyle \frac{}{q_x}}u_{n^{}𝐪}|u_{n𝐪}).`$ (18)
The proof is relatively straightforward, starting from Eq. (18) and discretizing integrals and derivatives on the mesh defined in Eq. (1).

Expressions such as Eq. (18) and similar ones had appeared in the literature before, in relationship to Wannier functions. By means of an expression of this kind, MV define a ground–state quantity $`\mathrm{\Omega }_\mathrm{I}`$ which sets a lower bound for the second (spherical) moments of the Wannier functions. More precisely, for an insulator with $`n_b`$ occupied bands (hence $`n_b`$ Wannier functions per cell) such second moment is no smaller in average than $`\mathrm{\Omega }_\mathrm{I}/n_b`$. It is worth mentioning at this point that the logic of the MV paper goes backwards with respect to the present approach: first they provide a continuum theory, and then they discretize for computational purposes. Their discretization is different from Eq. (16), which emerges naturally from the present formulation starting from the remarkably compact Eq. (14). Both discretizations obviously converge to the same $`M\mathrm{}`$ limit: their convergence properties are different, though.

Specializing MV to a cubic material, RS have found the simple relationship $`\mathrm{\Omega }_\mathrm{I}=3n_bx^2_\mathrm{c}`$: notice that $`x^2_\mathrm{c}`$ is intensive, while $`\mathrm{\Omega }_\mathrm{I}`$ is not such. Building upon MV’s work, we are now ready to generalize the localization tensor to materials of arbitrary symmetry:

$`r_\alpha r_\beta _\mathrm{c}={\displaystyle \frac{V_c}{n_b(2\pi )^3}}{\displaystyle }d𝐪({\displaystyle \underset{n}{}}{\displaystyle \frac{}{q_\alpha }}u_{n𝐪}|{\displaystyle \frac{}{q_\beta }}u_{n𝐪}`$ (19)
$`{\displaystyle \underset{n,n^{}}{}}u_{n𝐪}|{\displaystyle \frac{}{q_\alpha }}u_{n^{}𝐪}{\displaystyle \frac{}{q_\beta }}u_{n^{}𝐪}|u_{n𝐪}),`$ (20)
where $`V_c`$ is the cell volume. Notice that the imaginary part of the integrand in Eq. (20), being antisymmetric in $`𝐪`$, cancels in the integral, such that the localization tensor is real. Even the offdiagonal elements, as defined in Eq. (20), have a finite–$`N`$ counterpart in terms of many–body phase operators.

For an insulating crystal of arbitrary symmetry, $`\mathrm{\Omega }_\mathrm{I}`$ as defined by MV equals $`n_b`$ times the trace of our localization tensor $`r_\alpha r_\beta _\mathrm{c}`$. In a metal, expressions like Eqs. (18) and (20) do not make much sense, consistently with the fact that our finite–$`N`$ expression, Eq. (14), is formally infinite at any $`N`$ value.

## IV Properties of the localization tensor

We have already emphasized that the localization tensor is a property of the occupied KS manifold as a whole. The main quantity which defines such manifold is the (spin–integrated) single–particle density matrix $`\rho `$, which coincides with twice the projector $`P`$ over the occupied KS orbitals: this projector is invariant by unitary transformations of the orbitals. Using Bloch eigenfunctions the projector reads, for an insulator with $`n_b`$ occupied bands:

$$P(𝐫,𝐫^{})=\frac{1}{2}\rho (𝐫,𝐫^{})=\frac{1}{(2\pi )^3}\underset{n=1}{\overset{n_b}{}}𝑑𝐪\psi _{n𝐪}(𝐫)\psi _{n𝐪}^{}(𝐫^{}).$$
(21)
The localization tensor (in the thermodynamic limit) has been written as a Brillouin–zone integral in Eq. (20). This integral can be identically transformed into a particularly simple expression whose only ingredient is $`P`$:

$$r_\alpha r_\beta _\mathrm{c}=\frac{1}{2n_b}_{\text{cell}}𝑑𝐫_{\text{all space}}𝑑𝐫^{}(𝐫𝐫^{})_\alpha (𝐫𝐫^{})_\beta |P(𝐫,𝐫^{})|^2,$$
(22)
which is the second moment of the (squared) density matrix in the coordinate $`𝐫𝐫^{}`$. The proof of the equivalence between Eq. (22) and Eq. (20) can be worked out using the same algebra appearing in Ref. : for a different argument proving the same result, see the Appendix.

We have arrived at Eq. (22) considering an insulating crystal so far. In this case we know, under general arguments, that $`P(𝐫,𝐫^{})`$ is asymptotically exponential in the argument $`|𝐫𝐫^{}|`$: this confirms that the integral over all space in Eq. (22) converges and the localization tensor is therefore finite. At this point, it is worthwhile to apply the general form of Eq. (22) to the metallic case. For the simplest metal of all, the free electron gas, the density matrix is known exactly:

$$P(𝐫,𝐫^{})=\frac{1}{2}\rho (𝐫,𝐫^{})=\frac{3n_0}{2}\frac{j_1(k_\mathrm{F}|𝐫𝐫^{}|)}{k_\mathrm{F}|𝐫𝐫^{}|}.$$
(23)
Replacement of Eq. (23) into Eq. (22) results in a diverging integral, thus confirming that our localization tensor is formally infinite in this paradigmatic metal. Other, more realistic, metals feature this same divergence.

The fact that the density matrix $`\rho (𝐫,𝐫^{})`$ is short–range in the variable $`𝐫𝐫^{}`$ has been named “nearsightedness” by W. Kohn. The second moment expression in Eq. (22) shows that our localization tensor is indeed a meaningful quantitative measure of such nearsightedness. We are going to analyze below the major trends over an important class of materials: tetrahedral semiconductors. We mention at this point that a conceptually different measure of the nearsightedness of a given electronic ground state focuses instead on the exponent governing the exponential decay of $`\rho (𝐫,𝐫^{})`$ in insulators: some case studies have been recently investigated.

We have already observed that some of our findings are closely related to the previous work by MV. These authors’ main interest were the “optimally localized Wannier functions”, i.e. those localized orbitals which minimize the average spherical moment. They prove, among other things, that such moment is strictly larger than the trace of our localization tensor. Building on their results, it is straightforward to attribute a similar meaning to the tensor itself: for any transformation of the occupied orbitals into a set of unitarily equivalent ones, the second moment in a given direction can be no smaller than the localization tensor, projected on that direction.

Since we are going to apply our results to cubic materials only, we focus on those orbitals which minimize in average the quadratic spread (second moment) in the $`x`$ coordinate. The present formalism makes the definition of these orbitals particularly simple: they are in fact the eigenfunctions of the position operator $`x`$, projected over the occupied manifold. Calling $`\mathrm{\Xi }=PxP`$ this operator, its expression in the Schrödinger representation is:

$$\mathrm{\Xi }(𝐫,𝐫^{})=_{\text{all space}}𝑑𝐫^{\prime \prime }P(𝐫,𝐫^{\prime \prime })x^{\prime \prime }P(𝐫^{\prime \prime },𝐫^{}).$$
(24)
Notice that $`x`$ is incompatible with BvK boundary conditions and its matrix elements over Bloch states are ill defined; nonetheless, $`\mathrm{\Xi }`$ is—in insulators—a well defined operator, which maps any vector of the occupied manifold into another vector of the same manifold. This fact owes to the exponential localization of $`P`$ in Eq. (24).

The relationship between $`\mathrm{\Xi }`$ and the orbitals optimally localized in the $`x`$ direction is easily proved borrowing some results from MV; for a different argument leading to the same proof, see the Appendix. We also notice an important difference with respect to the three–dimensional localization explicitly considered by MV. While the trace of the localization tensor provides a lower bound for three–dimensional localization, its element $`x^2_\mathrm{c}`$ provides instead a genuine minimum for one–dimensional localization (in a cubic material). This qualitative difference owes to the fact that, while one can manifestly diagonalize $`PxP`$, one cannot simultaneously diagonalize $`PxP`$, $`PyP`$, and $`PzP`$.

We end this Section about the properties of the localization tensor with a most important issue: is $`r_\alpha r_\beta _\mathrm{c}`$ a measurable quantity? The answer, due to SWM, is “yes”. They prove the identity:

$$r_\alpha r_\beta _\mathrm{c}=\frac{\mathrm{}V_c}{2\pi e^2n_b}_0^{\mathrm{}}\frac{d\omega }{\omega }\mathrm{Re}\sigma _{\alpha \beta }(\omega ),$$
(25)
where $`\sigma _{\alpha \beta }`$ is the conductivity tensor. Notice that the left hand side, as emphasized throughout the present work, is a property of the electronic ground state, while the right hand side is a measurable property related to electronic excitations: therefore Eq. (25) must be regarded as a sum rule. The frequency integral in Eq. (25) diverges in metals and is finite in insulators, as obviously expected. Since in the latter materials there is a gap for electronic excitations, Eq. (25) immediately leads to the inequality:

$$r_\alpha r_\beta _\mathrm{c}<\frac{\mathrm{}V_c}{2\pi e^2n_b\epsilon _g}_0^{\mathrm{}}𝑑\omega \mathrm{Re}\sigma _{\alpha \beta }(\omega ),$$
(26)
where $`\epsilon _g`$ is the direct gap. Using then the oscillator–strength sum rule, Eq. (26) for a cubic material is cast as:

$$x^2_\mathrm{c}<\frac{\mathrm{}^2}{2m_e\epsilon _g}.$$
(27)
Below, we investigate the trends in both members of this inequality for our test–case materials.

## V Calculated localization tensors

We have studied several tetrahedrally coordinated crystalline materials, from the group IV, III–V, and II–VI, having the diamond and zincblende structure. The first-principle calculations have been performed within density-functional theory in the local-density approximation, using pseudopotentials  and plane waves. We implement a trivial extension of the formulas presented above, using a rectangular unit cell instead of a simple cubic one: we thus describe the diamond and zincblende structures by means of a tetragonal cell with a lattice constant $`a`$ in the basal plane and $`c=\sqrt{2}a`$. There are four atoms per unit cell, whose projections on the $`c`$ axis are equispaced; for the sake of consistency with the formal results, we take $`x`$ along the $`c`$ axis and $`yz`$ in the basal plane. We then use a BvK cell of sides $`M_xc`$, $`M_ya`$, and $`M_za`$, corresponding to a mesh in reciprocal space with $`M_x,M_y,M_z`$ points: this allows an easier control of convergence.

We start evaluating at the mesh points the Hermitian matrices

$$A_s=S^{}(𝐪_s,𝐪_s+\mathrm{\Delta }𝐪^{(x)})S(𝐪_s,𝐪_s+\mathrm{\Delta }𝐪^{(x)});$$
(28)
then Eqs. (15) and (16) are written as

$$x^2_\mathrm{c}=\frac{1}{M_yM_z}\underset{s_2=1}{\overset{M_y}{}}\underset{s_3=1}{\overset{M_z}{}}\left(\frac{c^2M_x}{4\pi ^2n_b}\underset{s_1=1}{\overset{M_x}{}}\text{ ln det }A_s\right).$$
(29)
In general convergence is fast in $`M_y`$, $`M_z`$, and slower in $`M_x`$. The expression in parenthesis in Eq. (29) is precisely the one–dimensional expression discussed in detail by RS, and the three dimensional one simply obtains from it as an average in the $`(q_y,q_z)`$ plane.

First we show in Fig. 1 the convergence of our expressions over a genuinely cubic grid, which coincides with the one used by MV in their evaluation of the quantity $`\mathrm{\Omega }_\mathrm{I}=3n_bx^2_\mathrm{c}`$. They use a different discretization of the same $`𝐤`$ space integral: both calculations converge to the same localization tensor, although our discretization, based on Eqs. (15) and (16), converges faster. All the following results have been obtained using noncubic grids, as in Eq. (29), in order to achieve faster convergence.

We have systematically calculated well converged localization tensors for several elemental and binary semiconductors. In Fig. 2 we plot the localization tensors versus the right–hand member of the inequality in Eq. (27), where for the gap $`\epsilon _g`$ we have used both (i) the theoretical and (ii) the experimental values. In case (i) the inequality owes to an exact sum rule and must be satisfied: we are therefore checking the internal consistency of the calculations. Also, it may be noticed that the inequality is very strongly verified. In case (ii) there is no a–priori guarantee that the inequality is verified, particularly given the fact that the experimental gap is systematically larger than the KS one. Nonetheless the localization tensor is obtained here as a pure ground state property, and it is well known that density-functional theory in the local-density approximation provides a good representation of the ground state, though not of the excitations. It is therefore interesting to verify that even for case (ii) the inequality in Eq. (27) is strongly verified for all the materials considered.

The localization tensor ranges roughly between 1 and 3 bohr<sup>2</sup> for all the materials considered, diamond being the most localized and germanium the most delocalized. The trend is qualitatively expected, in agreement with SWM’s statement that “the larger the gap, the more localized the electrons are”. However this is a trend more than a strict rule, and indeed a few materials show irregularities. Better trends are obtained when comparing families of materials: either isoelectronic series or isovalent series. In order to enhance such regularities, we have heuristically tried a few different laws. In Fig. 3 we plot our localization tensor versus $`1/\epsilon _g`$, using the minimum gaps instead of the direct ones: here monotonical trends are very perspicuous.

## VI Maximally localized hermaphrodite orbitals

We actually perform localizing transformations on the Bloch orbitals. At variance with the most standard approach, we focus on orbitals which are localized in one direction only, say $`x`$, while they are completely delocalized in the $`yz`$ directions. By analogy with the standard theory of Wannier functions, one obtains such orbitals by integration of the Bloch ones over one component only of the Bloch vector. The resulting orbitals are Wannier–like in one direction and Bloch–like in the other two: they can be therefore called “hermaphrodite orbitals”. Because of the same reasons as for ordinary Wannier functions, such hermaphrodite orbitals are nonunique: we focus here on those hermaphrodite orbitals which are optimally localized in the $`x`$ direction. It has been shown above that, in the thermodynamic limit, these orbitals are eigenfunctions of the operator $`\mathrm{\Xi }`$, Eq. (24), and their centers are the corresponding eigenvalues. It is expedient to consider the modified operator

$$\stackrel{~}{\mathrm{\Xi }}(𝐫,𝐫^{})=_{\text{all space}}𝑑𝐫^{\prime \prime }P(𝐫,𝐫^{\prime \prime })\mathrm{e}^{i\frac{2\pi }{M_xc}x^{\prime \prime }}P(𝐫^{\prime \prime },𝐫^{}),$$
(30)
which to leading order in $`1/M_x`$ has the same eigenfunctions as $`\mathrm{\Xi }`$, and simply related eigenvalues.

When considering a finite sample with BvK boundary conditions—or equivalently a discrete grid in the reciprocal unit cell—the operator $`\mathrm{\Xi }`$ as in Eq. (24) is useless because the operator $`x`$ therein becomes ill defined. Instead the operator $`\stackrel{~}{\mathrm{\Xi }}`$ is well defined, provided the value of $`M_x`$ is consistent with the choice of the grid. The integral in Eq. (30) is now performed over the BvK cell and not over all space; the projector projects over the finite occupied manifold, having dimension $`n_bM_xM_yM_z`$. Choosing the Bloch functions as the basis in the occupied manifold, the matrix elements of $`\stackrel{~}{\mathrm{\Xi }}`$ are nothing else than the matrix $`𝒮`$ defined in Eq. (9). Therefore in the discrete case the hermaphrodite orbitals which achieve optimal localization in the $`x`$ direction are simply obtained by diagonalizing the matrix $`𝒮`$. Since, as already observed, the matrix is already diagonal in $`s_2`$ and $`s_3`$, the problem is reduced to $`M_yM_z`$ independent diagonalizations of submatrices of size $`n_bM_x`$. We characterize our orbitals as $`w_{j,s_2,s_3}`$, where $`(s_2,s_3)`$ is a two–dimensional Bloch–index and $`j`$ is a one–dimensional Wannier–like index, with $`1jn_bM_x`$.

We are going to verify that these orbitals indeed minimize the average quadratic spread in one–dimension. If we define

$$z_{j,s_2,s_3}=_{\text{BvK cell}}𝑑𝐫|w_{j,s_2,s_3}(𝐫)|^2\mathrm{e}^{i\frac{2\pi }{M_xc}x},$$
(31)
then according to RS the quadratic spread of one given hermaphrodite orbital is:

$$w_{j,s_2,s_3}|x^2|w_{j,s_2,s_3}_\mathrm{c}=\frac{c^2M_x^2}{4\pi ^2}\mathrm{ln}|z_{j,s_2,s_3}|^2.$$
(32)
Taking now the average over all orbitals, and calling this quantity $`\lambda _{xx}^2`$, we get

$$\lambda _{xx}^2=\frac{1}{M_yM_z}\underset{s_2=1}{\overset{M_y}{}}\underset{s_3=1}{\overset{M_z}{}}\left(\frac{c^2M_x}{4\pi ^2n_b}\underset{j=1}{\overset{n_bM_x}{}}\mathrm{ln}|z_{j,s_2,s_3}|^2\right).$$
(33)
We then notice that $`w_{j,s_2,s_3}`$ are the eigenvectors of $`\stackrel{~}{\mathrm{\Xi }}`$, hence the expectation values $`z_{j,s_2,s_3}`$ are the corresponding eigenvalues. Since the product of the eigenvalues equals the determinant, standard manipulations prove that the average spread $`\lambda _{xx}^2`$ equals indeed the lower bound $`x^2_\mathrm{c}`$ as given in Eq. (29).

There is a subtlety about the diagonalization of the submatrices of $`𝒮`$, which are the projection over a certain finite dimensional manifold of the operator $`\mathrm{e}^{i\frac{2\pi }{M_xc}x}`$. Although the operator is unitary, its projection is not a unitary matrix, hence the eigenvectors are not exactly orthogonal, as instead honest localized orbitals must be: this is not a serious problem. In fact the larger is $`M_x`$, the closer to unitarity the matrix becomes: we know that the modulus of its determinant differs by one for a term of the order $`1/M_x`$, hence the modulus of each eigenvalue differs by one for a term of the order $`1/M_x^2`$. We recover exact orthonormality in the thermodynamic limit; in our calculations already at $`M_x20`$ deviations from orthogonality are hardly noticeable.

We have calculated the orbitals $`w_{j,s_2,s_3}`$ for several crystalline semiconductors: to the purpose of display, we call the $`yz`$ average of $`|w_{j,s_2,s_3}|^2`$ as $`n_{\mathrm{loc}}(x)`$, where the indices remain implicit. At fixed $`(s_2,s_3)`$ we have, given our double cell, $`8M_x`$ orbitals centered on a BvK period of length $`M_xc`$. There are however at most four different shapes, and one obtains all the functions upon translations in the $`x`$ direction (by multiples of $`c/2`$) of the four basic ones: this is not surprising, since the genuine unit cell is one half of our computational one. We find that the different shapes are actually always four, with the only exception of an elemental semiconductor at $`s_2=s_3=0`$. In this very special case the different shapes are only two, the orbitals are centered at the bond center, and their densities are centrosymmetric: the corresponding functions $`n_{\mathrm{loc}}(x)`$ are shown in Fig. 4 for the case of Si. The most general case is exemplified by Fig. 5: it shows the four different $`n_{\mathrm{loc}}(x)`$ for the case of GaAs, again at $`s_2=s_3=0`$. None of these four $`w`$ orbitals is therefore centered at a symmetry site, and none is centrosymmetrical, although they are obviously symmetrically related to each other. About the actual value of the quadratic spread (in the $`x`$ direction) of each of the $`w`$’s, we have found as a general feature that the least localized ones are those for $`s_2=s_3=0`$, i.e. at the $`\mathrm{\Gamma }`$ point in the two–dimensional reciprocal space.

We now address the long standing issue of exponential localization. Exact general results only exist for the genuine one–dimensional case, where W. Kohn has proved long ago that the Wannier functions which minimize the quadratic spread (i.e. are optimally localized) have an asymptotic exponential behavior. After the present work was completed, we became aware of Ref. , where the asymptotic exponential is shown to have a power–law prefactor. In three dimensions the problem is unsolved, with the exception of some very special cases. It has been conjectured by MV that their optimally localized Wannier functions enjoy three–dimensional exponential localization: an analytical proof looks very hard. Our hermaphrodite orbitals $`w`$ are optimally localized in the $`x`$ direction, and therefore in a sense they have one–dimensional character: nonetheless, they are have a genuine dependence on all three coordinates. Therefore an analytic proof along the lines of Refs.  and  is not easily extended to our case. Instead, it is simple to use the same arguments as given in Ref.  in order to prove that our $`w`$ orbitals decay in the $`x`$ direction faster that any inverse power of $`x`$.

Our very elongated BvK cells allow us to study the asymptotic behaviour heuristically on our calculated $`w`$’s. The nearest periodic replica of a given $`w`$ orbital is centered at a distance of $`M_xc`$ from it: therefore the interesting “asymptotic region” is accessible up to a distance somewhat smaller than $`M_xc/2`$, as it clearly appears from Fig. 6. The quantity of choice in order to “blow up” the exponential behavior is obviously the logarithm of $`n_{\mathrm{loc}}(x)`$, which we plot in Fig. 6 (thin solid lines) for the case of Ge and for two different $`(s_2,s_3)`$. It is seen that there is a wide region where the plots have an overall linear behavior, with superimposed oscillations having the crystal periodicity along the $`x`$ direction ($`c/2`$ in the present case). The slopes at different $`(s_2,s_3)`$ are very different, though; the $`n_{\mathrm{loc}}(x)`$ with the slowest decay corresponds to $`s_2=s_3=0`$ and therefore to the least localized, as we previously observed. Next, we filter the disturbing periodic oscillations using our favorite tool of the macroscopic average. We tried both ways: filtering $`n_{\mathrm{loc}}(x)`$ first and then taking the logarithm, or filtering $`\mathrm{ln}n_{\mathrm{loc}}(x)`$: the latter turns out to work best. The macroscopic filtering is also shown in Fig. 6 (thick solid lines): it is easily realized, expecially looking at the magnified plot in the lower panel, that there is a sizeable region, spanning several cells, where the plotted function looks accurately linear with $`x`$, hence $`n_{\mathrm{loc}}(x)\mathrm{exp}(\pm bx)`$. We therefore demonstrate “experimentally” the exponential localization of our $`w`$ orbitals. After we became aware of Ref. , we checked that the power–law prefactor suggested therein does not improve the quality of our fits. It is hard to assess whether this is due to a basic difference between our case and a genuine one–dimensional one, or to the limited resolution achievable in our selfconsistent three-dimensional finite–size calculation. Finally, in Fig. 7 we display some correlations between the localization length and the exponential decay length $`1/b`$ averaged over the two-dimensional mesh ($`s_2,s_3`$).

## VII Conclusions

In the present work we provide the three–dimensional formulation of the RS theory of electron localization, specializing it to the case of independent KS electrons; we discuss some of its relationships to the MV and SWM papers, and also (in the Appendix) how it relates to Boys theory of localization in molecules. We then implement the theory to several materials in the class of tetrahedrally coordinated semiconductors. Among the results, we find that in general the calculated localization length is a monotonical function of the gap, although a few materials show irregularities. The trend is more regular within a given family (isoelectronic or isovalent). Finally, we heuristically show that the orbitals which are optimally localized in a given direction (“hermaphrodite orbitals”) show exponential localization.

## Acknowledgments

R.R. acknowledges partial support by ONR grant N00014-96-1-0689.

## Relationship to Boys localization in molecules

We abandon here the BvK boundary conditions used throughout this work, and we consider an $`N`$–electron system which is bounded in space. Both the orbitals and the many–body wavefunction $`\mathrm{\Psi }`$ are therefore exponentially vanishing at large distances. Supposing that $`N`$ is even and the state is a singlet, for independent particles the wavefunction is the Slater determinant:

$$\mathrm{\Psi }=\frac{1}{N!}|\phi _1\overline{\phi }_1\phi _2\overline{\phi }_2\mathrm{}\phi _{N/2}\overline{\phi }_{N/2}|.$$
(34)
The orbitals enjoy no specific symmetry. Any unitary transformation of the orbitals produces the same many–body ground state (modulo an overall phase): a specific choice of the orbitals will be referred to as “choice of the gauge” in the following. Obviously all ground state properties are gauge–invariant. The density matrix is twice the projector over the occupied orbitals:

$$\rho (𝐫,𝐫^{})=2P(𝐫,𝐫^{})=2\underset{i=1}{\overset{N/2}{}}\phi _i(𝐫)\phi _i^{}(𝐫^{}).$$
(35)
We are interested in exploiting the gauge freedom in order to express the ground state in terms of localized orbitals. The standard Boys localization in molecules consists in minimizing spherical second moments, in perfect analogy with MV, which can be regarded as the solid–state analogue of Boys localization. Here instead we are mostly interested in localizing in one given direction, say $`x`$.

For any given choice of the single-particle orbitals $`\phi _i`$, the average quadratic spread in the $`x`$ direction is by definition:

$$\lambda _{xx}^2=\frac{2}{N}\underset{i=1}{\overset{N/2}{}}(\phi _i|x^2|\phi _i\phi _i|x|\phi _i^2).$$
(36)
We recast this identically as:

$`\lambda _{xx}^2={\displaystyle \frac{2}{N}}{\displaystyle \underset{i}{}}\phi _i|x(\mathrm{\hspace{0.17em}1}{\displaystyle \underset{j}{}}|\phi _j\phi _j|)x|\phi _i`$ (37)
$`+{\displaystyle \frac{2}{N}}{\displaystyle \underset{ij}{}}|\phi _i|x|\phi _j|^2.`$ (38)
The first term in Eq. (38) is gauge invariant, since we can identically write:

$$\lambda _{xx}^2=\frac{2}{N}\text{Tr }xPx(1P)+\frac{2}{N}\underset{ij}{}|\phi _i|x|\phi _j|^2,$$
(39)
where “Tr” indicates the trace on the electronic coordinate. The gauge–invariant term in Eq. (39) can be regarded as the $`xx`$ element of a more general tensor, which turns out to be the molecular analogue of our localization tensor. We use the same notation for molecules and for crystals:

$$r_\alpha r_\beta _\mathrm{c}=\frac{2}{N}\text{Tr }r_\alpha Pr_\beta (1P).$$
(40)
If we look for the orbitals which minimize the average spread in the $`x`$ direction, the solution, after Eq. (39), is provided by those orbitals which diagonalize the position operator $`x`$, projected over the occupied manifold. Obviously, a set of orthonormal orbitals which diagonalize $`PxP`$ can always be found, since $`PxP`$ is a Hermitian operator. The quadratic spread of these orbitals equals $`x^2_\mathrm{c}`$, the gauge–invariant part in Eq. (39). If we are interested instead in minimizing the spherical second moment, in general we cannot diagonalize simultaneously $`PxP`$, $`PyP`$, and $`PzP`$. Therefore the spherical spread will be in general strictly larger than the Cartesian trace of the localization tensor. This is a key feature in the work of Boys, and MV as well.

We have defined the localization tensor in Eq. (40). With an obvious generalization of the previous arguments, this tensor provides in general the maximum localizability in any given direction. An equivalent expression for the localization tensor is:

$$r_\alpha r_\beta _\mathrm{c}=\frac{1}{N}𝑑𝐫𝑑𝐫^{}(𝐫𝐫^{})_\alpha (𝐫𝐫^{})_\beta |P(𝐫,𝐫^{})|^2,$$
(41)
which has the meaning of the second moment of the (squared) density matrix in the coordinate $`𝐫𝐫^{}`$.

At this point, we may think of a crystalline solid as of a very large “molecule”, or a cluster, and take the thermodynamic limit. Since bulk properties must be independent of the choice of boundary conditions (either BvK or “free”), the density matrix and the localization tensor must be the same as the one previously found in this work. And indeed, a glance to Eq. (22) shows that it coincides with the thermodynamic limit of Eq. (41) in the insulating case. As for the metallic case, our previous findings bear an important message concerning Boys localization. For a cluster of finite size, no matter how large, one can doubtless build localized Boys orbitals. But our results prove that, in the large $`N`$ limit, the quadratic spread of these Boys orbitals diverges whenever the cluster is metallic.