text:
# Density Functional Theory for a Confined Fermi System with Short-Range Interaction

## I Introduction

Calculations of bulk observables for medium to heavy nuclei typically rely on nonrelativistic (Skyrme) or covariant (QHD) mean-field models. How can one systematically go beyond such models? One possibility is to view nuclear mean-field approaches as approximate implementations of Kohn-Sham density functional theory (DFT) KOHN65 ; PARR89 ; DREIZLER90 ; KOHN99 ; ARGAMAN00 , which is widely used in condensed matter and quantum chemistry applications. To date, the refinement of DFT methods has focused almost exclusively on Coulomb systems and DFT has had little explicit impact on nuclear structure phenomenology PETKOV91 . Effective field theory (EFT), however, offers a systematic approach for describing low-energy nuclear physics that can provide a framework for nuclear DFT. EFT approaches have been making steady progress on two- and three-body nuclear systems and certain halo nuclei, and ab initio shell model methods should be able to extend the calculations to a wide range of light nuclei EFT98 ; EFT99 ; Birareview ; BEANE99 ; Furnstahl:2001hs . Our ultimate goal is to systematically describe *heavier* nuclei by applying EFT methods to Kohn-Sham DFT. In this paper we take the first steps toward this goal.

Density functional theory provides a calculational framework that greatly extends the range of many-body calculations in finite systems. In Kohn-Sham DFT, the coordinate-space ground-state density for $`A`$ fermions is found by simply summing the squared wavefunctions for a set of single-particle orbitals. These wavefunctions are solutions to a Schrödinger equation featuring the Kohn-Sham single-particle potential, which is local and energy independent. At zero temperature, for nondegenerate ground states (e.g., closed shells), the sum for the density is over the orbitals with the lowest $`A`$ eigenvalues, equally weighted (i.e., “occupation numbers” are either one or zero). The Kohn-Sham potential is itself a functional of the orbitals, so we have a self-consistent problem that is solvable by iteration. When the iterations have converged, the orbitals can be plugged into an energy functional to find the ground-state energy. While this procedure is very similar to the self-consistent Hartree approximation, it incorporates *all* correlations *if* the correct Kohn-Sham functional is used.

The advantages of the Kohn-Sham DFT procedure are clear. One is that the external potential, if present, appears in a very simple way in a separate term in the functional, with the rest of the functional being universal. Another is that solving the Schrödinger equation with local potentials for eigenvalues is relatively simple and fast. The computational advantages would be lost, of course, if the construction and evaluation of the Kohn-Sham potential is itself too difficult or expensive. The experience for Coulomb systems is that local density approximations (LDA) for the potential part of the energy functional and their improvement with gradient expansions work very well in many (although not all) systems PERDEW96 ; PERDEW99 ; ARGAMAN00 . That is, the most important nonlocality to treat explicitly is the kinetic energy, which is why the Kohn-Sham framework is generally superior to the Thomas-Fermi approximation, in which the kinetic energy is given by an LDA (see Sect. IV.4 below). The gradient expansion approximations greatly simplify the evaluation of the Kohn-Sham potential. We anticipate that the convergence of derivative expansions for nuclear systems will be rapid in general, but this will need to be confirmed.

Our strategy is to apply an effective action formalism COLEMAN88 ; PESKIN95 ; WEINBERG96  to calculate the Kohn-Sham potential and energy functional order-by-order in an EFT expansion, and to use EFT power counting to organize and justify a derivative expansion of the functional. In the present work, we use a dilute, confined Fermi system with short-range interactions as a laboratory to explore how EFT can be used to carry out systematic DFT calculations. In this case, the explicit expansion parameter is the local Fermi momentum times the scattering length (and other effective range parameters). We assume a gradient expansion parameter that justifies a local density approximation, but the verification of this assumption is postponed to future work. Ultimately we are interested in calculating self-bound systems (e.g., nuclei), with spin- and isospin-dependent interactions and long-range forces (e.g., pion exchange). These are all significant but well-defined extensions of the model described here. In the meantime, the model provides a prototype for more complex systems and also has a physical realization in recent and forthcoming experiments on fermionic atoms in optical traps GRANADE02 .

The Kohn-Sham approach to DFT was proposed in Ref. KOHN65 . Since then, the literature of DFT applications has grown exponentially, primarily in the areas of quantum chemistry and electronic structure ARGAMAN00 . A general introduction to density functional theory as conventionally applied is provided in the books by Dreizler and Gross DREIZLER90  and Parr and Young PARR89 , while Ref. KOCH00  is a practitioners guide to DFT for quantum chemists. The connection of DFT to nonrelativistic mean-field approaches to nuclei (e.g., Skyrme models) was pointed out in Ref. BRACK85  (and no doubt elsewhere) and was explored for covariant nuclear mean-field models in Refs. SCHMID95 ; SCHMID95a . However, it has not led, to our knowledge, to new or systematically improved mean-field-type functionals for nuclei.

The use of functional Legendre transformations for DFT with the effective action formalism was first detailed by Fukuda and collaborators FUKUDA94 ; FUKUDA95 , who also discuss the inversion and auxiliary field methods of constructing the effective action. The connection to Kohn-Sham DFT was shown by Valiev and Fernando VALIEV96 ; VALIEV97 ; VALIEV97b ; RASAMNY98  and later by other authors in Refs. FAUSSURIER00 ; CHITRA00 ; CHITRA01 . Recent work by Polonyi and Sailer applies renormalization group methods and a cluster expansion to an effective-action formulation of generalized DFT for Coulomb systems Polonyi:2001uc . To our knowledge, however, there is no prior work on merging the Kohn-Sham density functional approach and effective field theory.

The plan of the paper is as follows. In Sect. II, we review effective field theory for a dilute system of fermions. In Sect. III, the effective action approach for determining a Kohn-Sham energy functional through a Legendre transformation is reviewed. A systematic approximation procedure for constructing the energy functional, the inversion method, is presented. In Sect. IV, the formalism is applied to a dilute Fermi gas in a harmonic trap and results are presented through third order (NNLO) in the dilute EFT expansion using an LDA. Section V summarizes our results and future plans.

## II EFT for Infinite Dilute Fermi Systems

### II.1 Background

Effective field theory (EFT) provides a powerful framework to study low-energy phenomena in a model-independent way LEPAGE89 ; Birareview ; BEANE99 . The EFT approach is grounded in some very general physical principles LEPAGE89 . If a system is probed or interacts at low energies, resolution is also low, and fine details of what happens at short distances or in high-energy intermediate states are not resolved. Therefore, the short-distance structure can be replaced by something simpler without distorting the low-energy observables. This is analogous to a multipole expansion, in which a complicated, charge or current distribution is replaced for long-wavelength probes by a series of point multipoles. EFT uses local Lagrangian field theory as a framework for carrying out this program in a complete and systematic way.111Note that conventional nuclear phenomenology also relies on these principles in using potentials cut off at short distances. However, cutoff *independence* is a goal of EFT that is not usually achieved in phenomenological approaches. The uncertainty principle implies that high-energy intermediate states are highly virtual and only last for a short time, so their effects are not distinguishable from those of local operators LEPAGE89 . This physics can then be systematically absorbed into the coefficients of these operators using renormalization.

The effective degrees of freedom (dof’s) in an EFT depend on a separation or resolution momentum scale $`\mathrm{\Lambda }`$, which sets the radius of convergence of an EFT. Long-range dof’s with respect to $`\mathrm{\Lambda }`$ must be treated explicitly while short-range physics is encoded in the coefficients of the local operators. (See Ref. BEANE99  for details of how the pion can be considered either a short- or long-distance degree of freedom in the two-nucleon problem, depending on the resolution scale.) The hierarchy of scales in the system is exploited to provide expansion parameters. For example, if the typical momenta $`k`$ are small compared to the inverse range of the interaction $`1/R`$, we can take $`\mathrm{\Lambda }1/R`$ and low-energy observables can be described by a controlled expansion in $`kR`$. All short-distance effects are systematically absorbed into low-energy constants through renormalization. As a result, the EFT approach allows for accurate calculations of low-energy processes and properties with well-defined error estimates (based on the order of truncation) EFT98 ; EFT99 ; BEANE99 ; Birareview .

The application of EFT methods to many-body problems promises a consistent organization of many-body corrections, with reliable error estimates, and insight into the analytic structure of observables (see, for example, the identification using renormalization group methods of logarithmic contributions to the energy of dilute systems in Refs. BRAATEN97  and HAMMER00 ). The EFT provides a model-independent description of finite-density observables in terms of parameters that can be fixed from scattering in the vacuum or from a subset of finite density properties. One can also exploit the freedom in an EFT of using different regulators and renormalization schemes to find simplifications and clarifications HAMMER00 .

While EFT has shown early promise in applications to basic many-body problems (e.g., Refs. BRAATEN97  and HAMMER00 ), there are formidable challenges in carrying out many-body calculations, particularly for finite, non-uniform systems. For sufficiently small numbers of fermions, the many-body Schrödinger equation for finite systems can be solved directly, for example by Green’s function Monte Carlo methods. However, the computational cost of these methods grows as a power of the number of particles (or faster), which prevents their application to very many systems of interest in condensed matter and quantum chemistry (Coulomb systems) and nuclear physics (medium to heavy nuclei).

The purpose of the present work is to address these challenges systematically. We seek to merge the organizational advantages and insight provided by EFT with the calculational power and relative ease of DFT for finite systems. We will build the basic EFT formalism of this merger for a dilute gas of identical fermions with short-range interactions. We review below the results from the calculation of the energy density in the case of a uniform system HAMMER00 . These results will be the starting point of the DFT calculation of the energy density of a dilute Fermi gas confined by an external harmonic potential. They will also serve as the basis for the LDA to the energy functional.

### II.2 Lagrangian and Energy Density for Uniform System

We consider a general local Lagrangian for a nonrelativistic fermion field that is invariant under Galilean, parity, and time-reversal transformations:

$``$ $`=`$ $`\psi ^{}\left[i_t+\mu +{\displaystyle \frac{\stackrel{}{}^{\mathrm{\hspace{0.17em}2}}}{2M}}\right]\psi {\displaystyle \frac{C_0}{2}}(\psi ^{}\psi )^2+{\displaystyle \frac{C_2}{16}}\left[(\psi \psi )^{}(\psi \stackrel{}{}^2\psi )+\text{ h.c.}\right]`$ (1)
$`+{\displaystyle \frac{C_2^{}}{8}}(\psi \stackrel{}{}\psi )^{}(\psi \stackrel{}{}\psi )+\mathrm{},`$
where $`\stackrel{}{}=\stackrel{}{}\stackrel{}{}`$ is the Galilean invariant derivative and h.c. denotes the Hermitian conjugate. The terms proportional to $`C_2`$ and $`C_2^{}`$ contribute to $`s`$-wave and $`p`$-wave scattering, respectively, while the dots represent terms with more derivatives and/or more fields. The Lagrangian Eq. (1) represents a particular canonical form, which can be reached via field redefinitions. For example, higher-order terms with time derivatives are omitted, as they can be eliminated in favor of terms with spatial derivatives HAMMER00 .

To reproduce the results in Ref. HAMMER00 , we can write a conventional generating functional with the Lagrangian of Eq. (1) and Grassmann sources coupled to $`\psi ^{}`$ and $`\psi `$, respectively NEGELE88 . The non-quadratic part of the Lagrangian is removed in favor of functional derivatives with respect to the Grassmann sources and the remaining quadratic part is evaluated in terms of a non-interacting Green’s function times the sources. Perturbative expansions for Green’s functions (and subsequently S–matrix elements) follow by taking successive functional derivatives, and the ground state energy density follows by applying the linked cluster theorem (see Ref. NEGELE88  for details). In calculating the energy, finite density boundary conditions at $`T=0`$ can be incorporated into the non-interacting Green’s function using the chemical potential $`\mu `$ or by including them by hand with a non-interacting chemical potential. The latter approach is simplest at $`T=0`$ and was adopted in Ref. HAMMER00 .

The coefficients $`C_0`$, $`C_2`$, and $`C_2^{}`$ can be obtained from matching the EFT to a more fundamental theory or to (at least) three independent pieces of experimental data. We follow the regularization and renormalization prescription described in Ref. HAMMER00 , namely dimensional regularization with minimal subtraction, which is particularly convenient for the dilute, natural system. By matching to the effective-range expansion for low-energy fermion-fermion scattering, we can express the $`C_{2i}`$ in terms of the effective-range parameters:

$$C_0=\frac{4\pi a_s}{M},C_2=C_0\frac{a_sr_s}{2},\text{and}C_2^{}=\frac{4\pi a_p^3}{M},$$
(2)
where $`a_s`$ ($`a_p`$) are the $`s`$-wave ($`p`$-wave) scattering length and $`r_s`$ is the $`s`$-wave effective range, respectively.

The energy density for the uniform dilute fermion system with natural scattering length HAMMER00  is calculated as a perturbative expansion in $`k_\mathrm{F}/\mathrm{\Lambda }`$, where $`k_\mathrm{F}`$ is the Fermi momentum and $`\mathrm{\Lambda }`$ is the resolution scale (e.g., $`\mathrm{\Lambda }1/R`$ for hard spheres). The non-interacting energy density at zero temperature for $`A`$ particles with spin-degeneracy $`g`$ in volume $`V`$ can be written as

$$_0=\rho \frac{3}{5}\frac{k_\mathrm{F}^2}{2M},$$
(3)
where the density $`\rho `$ is

$$\rho =\frac{A}{V}=g\frac{d^3k}{(2\pi )^3}\theta (k_\mathrm{F}k)=\frac{gk_F^3}{6\pi ^2}.$$
(4)
The order-by-order corrections to Eq. (3) due to interactions can be represented by the Hugenholtz diagrams given in Fig. 1. The Feynman rules for calculating these graphs along with details of the renormalization needed to render divergent graphs finite are given in Ref. HAMMER00 . Here we will only need the final results and we simply quote them up to next-to-next-to-leading order (NNLO).

The LO diagram of order $`k_\mathrm{F}^6`$ \[Fig. 1(a)\], which represents the Hartree-Fock result, contributes

$$_1=\rho (g1)\frac{k_\mathrm{F}^2}{2M}\frac{2}{3\pi }k_\mathrm{F}a_s,$$
(5)
where Eqs. (2) and (4) have been applied. At order $`k_\mathrm{F}^7`$ (NLO) there are two diagrams. The three-loop diagram \[Fig. 1(c)\] is an example of an “anomalous diagram,” which vanishes identically for a uniform system in the zero-temperature formalism but is nonzero when calculated in the zero-temperature limit. In the latter case, its contribution is precisely canceled by the shift between the noninteracting and interacting chemical potentials, as dictated by the Kohn-Luttinger-Ward theorem KOHN60 ; FETTER71 ; NEGELE88 . We will find an analogous cancellation in the diagrammatic expansion of the Kohn-Sham functional. The other diagram at this order \[the “beach ball” diagram of Fig. 1(b)\], makes the contribution

$$_2=\rho (g1)\frac{k_\mathrm{F}^2}{2M}(k_\mathrm{F}a_s)^2\frac{4}{35\pi ^2}\left(112\mathrm{ln}2\right).$$
(6)
Finally, we have the graphs of order $`k_\mathrm{F}^8`$ (NNLO). The first three are anomalous diagrams. In addition to graphs containing the $`C_0`$ vertex \[Figs. 1(g) and (h)\], there also graphs that contain $`C_2`$ and $`C_2^{}`$. Altogether these graphs give the correction

$`_3`$ $`=`$ $`\rho {\displaystyle \frac{k_\mathrm{F}^2}{2M}}[(g1){\displaystyle \frac{1}{10\pi }}(k_\mathrm{F}a_s)^2k_\mathrm{F}r_s+(g+1){\displaystyle \frac{1}{5\pi }}(k_\mathrm{F}a_p)^3`$ (7)
$`+(g1)\{(0.07550\pm 0.00003)+(g3)(0.05741\pm 0.00002)\}(k_\mathrm{F}a_s)^3],`$
where the integrals for Figs. 1(g) and (h) have been evaluated numerically. The expansion can be continued systematically, but this is as far as we will need here.

## III Effective Actions and The Inversion Method

The energy density given in the last section illustrates that the ground state energy of a *uniform*, interacting many-body system may be written as a function of the constant density. The DFT theorems of Hohenberg and Kohn (HK) HK64  formally prove the existence of a generalization to finite, *nonuniform* systems: the ground state energy can be obtained from a functional of the local density alone. Specifically, for a system with an external potential $`v(𝐱)`$ coupled to the density $`\rho (𝐱)`$, there exists an energy functional $`E[\rho ]`$ that can be decomposed as

$$E[\rho (𝐱)]=F_{HK}[\rho (𝐱)]+d^3𝐱v(𝐱)\rho (𝐱),$$
(8)
where the functional $`F_{HK}[\rho ]`$, which is known as the HK free energy, is universal (i.e., independent of the potential $`v`$). A variational principle ensures that the functional $`E[\rho ]`$ is a minimum equal to the ground state energy when evaluated at the exact ground state density. The practical problem of density functional theory is to find accurate and tractable approximations to $`F_{HK}[\rho ]`$ DREIZLER90 ; PARR89 ; ARGAMAN00 .

In applications to Coulomb systems, $`F_{HK}[\rho ]`$ is typically decomposed into a noninteracting kinetic energy $`F_{ni}[\rho ]`$, the Hartree term $`E_H[\rho ]`$, and everything else, which is defined as the exchange-correlation energy $`E_{xc}[\rho ]`$ ARGAMAN00 :

$$F_{HK}[\rho ]=F_{ni}[\rho ]+E_H[\rho ]+E_{xc}[\rho ].$$
(9)
In the Thomas-Fermi approximation, the noninteracting functional is evaluated within an LDA. This approximation is accurate at high densities when the characteristic scale over which the fermion density changes is small compared to the Fermi wavelength. The approximation is inadequate for most ordinary Coulomb systems and improves only when some account is taken of the nonlocality in $`F_{ni}`$. Kohn and Sham KOHN65  introduced a method involving auxiliary orbits by which $`F_{ni}`$ could be treated exactly, which has become the basis for most practical DFT applications ARGAMAN00 ; KOHN99 .

The utility of DFT then rests on finding an explicit expression, exact or approximate, for $`E_{xc}`$. For Coulomb systems, the conventional procedure is to approximate $`E_{xc}`$ in the LDA based on a fit to Monte Carlo results for the energy of a uniform electron gas as a function of the density, and then to include semi-phenomenological gradient corrections ARGAMAN00 . For applications to other systems, such as nuclei or trapped atoms, we adopt an approximation scheme based on effective action methods. Here we follow the work of Fukuda et al. FUKUDA94 ; FUKUDA95  and its further development by Valiev and Fernando VALIEV97 . We extend this work by merging it with an effective field theory expansion. For consistency with the EFT treatment of dilute systems in Ref. HAMMER00 , we work at zero temperature in Minkowski space.

We begin with the system described by the Lagrangian in Eq. (1), to which we add a term for an external potential $`v(𝐱)`$ coupled to the density operator, $`v(𝐱)\psi ^{}\psi `$.222In general, the potential will be coupled to more operators than simply $`\psi ^{}\psi `$. Field redefinitions, which leave observables unchanged, could be used to eliminate couplings to these operators, which would be higher order in the EFT expansion FURNSTAHL01 . However, such field transformations would also induce energy-dependent terms in the Lagrangian. Since we have already used field transformations to achieve a canonical, energy-independent Lagrangian, we will assume that omitted higher-order couplings, which should be suppressed numerically according to power counting, can be treated perturbatively. Here we take the external potential to be an isotropic harmonic confining potential

$$v(𝐱)=\frac{1}{2}m\omega ^2|𝐱|^2,$$
(10)
as might be appropriate for some atomic traps BRUUN98 , but the discussion holds for a general external potential. We also introduce a c-number source, $`J(x)`$, coupled to the composite density operator and write down the generating functional using the path integral formulation,

$$Z[J]=e^{iW[J]}=D\psi D\psi ^{}e^{i{\scriptscriptstyle d^4x[+J(x)\psi ^{}(x)\psi (x)]}}.$$
(11)
For simplicity, normalization factors are considered to be implicit in the functional integration measure. (See Refs. FUKUDA94 ; FUKUDA95  for a more careful treatment of the path integrals.) Using the definition in Eq. (11), we see that the density (in the presence of $`J`$) is

$$\rho (x)\psi ^{}(x)\psi (x)_J=\frac{\delta W[J]}{\delta J(x)}.$$
(12)
While it is possible to absorb $`v(𝐱)`$ into the definition of $`J(x)`$, we find it more convenient to recover our original system in the limit $`J0`$.

The effective action is defined through the functional Legendre transformation

$$\mathrm{\Gamma }[\rho ]=W[J]d^4xJ(x)\rho (x).$$
(13)
This transformation ensures that $`\mathrm{\Gamma }`$ has no explicit dependence on $`J`$. The existence of such a functional follows from the concavity of $`W[J]`$ which guarantees that Eq. (12) can be inverted to give $`J=J[\rho ]`$. The functional dependence between $`J`$ and $`\rho `$ can be used to write Eq. (13) entirely in terms of $`\rho `$. The proof that $`W[J]`$ is strictly concave is given in Ref. VALIEV97 .

Since we are interested here in describing finite many-body ground states at $`T=0`$, it is most convenient to work with functions of the particle number $`A`$ rather than the chemical potential $`\mu `$. This can be achieved via a conventional Legendre transformation on $`W`$ or $`\mathrm{\Gamma }`$, but is more simply carried out implicitly, by choosing appropriate finite-density boundary conditions that enforce a given $`A`$ by hand (the actual procedure is detailed below). In the following, we will assume this has been done. *Thus, $`\mathrm{\Gamma }`$ and $`W`$ are functions of $`A`$ and variations over $`\rho (x)`$ conserve $`A`$.* In addition, we restrict the discussion to time independent sources and densities. In this case the effective action acquires a factor that corresponds to the time interval over which the source is acting, which we indicate schematically as

$$\mathrm{\Gamma }[\rho ]=E[\rho ]\times _{\mathrm{}}^{\mathrm{}}𝑑t,$$
(14)
as in Ref. FUKUDA94 . To avoid overly cluttered notation, we will divide out this ubiquitous time factor everywhere it appears and write

$$\stackrel{~}{\mathrm{\Gamma }}[\rho ]\mathrm{\Gamma }[\rho ]\times \left[_{\mathrm{}}^{\mathrm{}}𝑑t\right]^1=E[\rho ].$$
(15)
and similarly with $`W[J]`$ and the expansions below. (We will continue to use $`\stackrel{~}{\mathrm{\Gamma }}`$ rather than $`E`$ in this section.)

In conventional treatments (e.g., see Refs. WEINBERG96 ; PESKIN95 ), an effective action is derived from a Legendre transformation with respect to a source coupled to one of the fields in the Lagrangian, rather than to a composite operator as in the present case. However, the usual advantages of working with an effective action are also present here. The effective action has extrema at the possible quantum ground states of the system, and when evaluated at the minimum is proportional (at zero temperature) to the ground state energy FUKUDA94 ; FUKUDA95 ; HU96 . In particular, Eq. (14) defines an energy functional $`E[\rho ]`$ equal to the ground-state energy when evaluated with the ground-state density $`\rho `$.

The extremization condition is shown as follows. Combining Eq. (12) and Eq. (13) we find

$$\frac{\delta \stackrel{~}{W}[J]}{\delta J(𝐱)}=d^3𝐲\left(\frac{\delta \stackrel{~}{\mathrm{\Gamma }}[\rho ]}{\delta \rho (𝐲)}\right)\left(\frac{\delta \rho (𝐲)}{\delta J(𝐱)}\right)+\rho (𝐱)+d^3𝐲\left(\frac{\delta \rho (𝐲)}{\delta J(𝐱)}\right)J(𝐲)$$
(16)
or

$$d^3𝐲\left(\frac{\delta \stackrel{~}{\mathrm{\Gamma }}[\rho ]}{\rho (𝐲)}+J(𝐲)\right)\left(\frac{\delta \rho (𝐲)}{\delta J(𝐱)}\right)=0.$$
(17)
The invertibility of (12) implies

$$\frac{\delta \rho (𝐲)}{\delta J(𝐱)}0,$$
(18)
so we must have

$$\frac{\delta \stackrel{~}{\mathrm{\Gamma }}[\rho ]}{\delta \rho (𝐱)}=J(𝐱).$$
(19)
The above equation tell us that when $`J(𝐱)=0`$ the effective action is extremized, which is a statement of the second HK theorem HK64 . The strict concavity of $`\stackrel{~}{W}[J]`$ implies the strict concavity of $`\stackrel{~}{\mathrm{\Gamma }}[\rho ]`$ and so the extremum is a maximum. Since $`J(𝐱)=0`$ corresponds to the original system we see that the energy functional is minimized when evaluated at the exact expectation value of the density.

We observe that the separation of a $`v(𝐱)`$-dependent part of the DFT energy functional (or $`\stackrel{~}{\mathrm{\Gamma }}`$) from a universal part follows directly in the effective action formalism FUKUDA94 . From the definitions of $`Z`$ and $`W`$, it follows that

$$\stackrel{~}{W}_{v=0}[J]=\stackrel{~}{W}[J+v]$$
(20)
for any $`J(𝐱)`$. If we designate $`J_\rho (𝐱)`$ the inversion of $`\delta \stackrel{~}{W}/\delta J=\rho `$ and $`J_\rho ^0(𝐱)`$ the inversion of $`\delta \stackrel{~}{W}_{v=0}/\delta J=\rho `$ for the same density, then

$$\frac{\delta \stackrel{~}{W}[J_\rho ]}{\delta J(𝐱)}=\frac{\delta \stackrel{~}{W}_{v=0}[J_\rho ^0]}{\delta J(𝐱)}=\frac{\delta \stackrel{~}{W}[J_\rho ^0+v]}{\delta J(𝐱)},$$
(21)
which implies

$$J_\rho =J_\rho ^0+v.$$
(22)
Upon substituting Eq. (22) into Eq. (13), the effective action becomes

$`\stackrel{~}{\mathrm{\Gamma }}[\rho ]`$ $`=`$ $`\stackrel{~}{W}[J_\rho ]{\displaystyle d^3𝐱J_\rho (𝐱)\rho (𝐱)}`$ (23)
$`=`$ $`\stackrel{~}{W}_{v=0}[J_\rho ^0]{\displaystyle d^3𝐱\left(J_\rho ^0(𝐱)+v(𝐱)\right)\rho (𝐱)}`$
$`=`$ $`\left[\stackrel{~}{W}_{v=0}[J_\rho ^0]{\displaystyle d^3𝐱J_\rho ^0\rho (𝐱)}\right]{\displaystyle d^3𝐱v(𝐱)\rho (𝐱)}`$
$`=`$ $`\stackrel{~}{\mathrm{\Gamma }}_{v=0}[\rho ]{\displaystyle d^3𝐱v(𝐱)\rho (𝐱)},`$
which is the promised result \[Eq. (8)\] with an overall minus sign.

Now consider a system that can be characterized by an effective field theory power-counting parameter, which we label $`\lambda `$. This parameter may be dimensionful and can appear to all orders. We will not exhibit it explicitly here but indicate by subscripts the order in $`\lambda `$ of a given function or functional. In previous discussions of the inversion method FUKUDA94 ; FUKUDA95 , the parameter $`\lambda `$ was always a coupling constant (e.g., $`e^2`$ for the Coulomb interaction). In contrast, we associate $`\lambda `$ with an appropriate EFT expansion parameter. For example, $`\lambda `$ could be $`1/\mathrm{\Lambda }`$ in the dilute expansion (which we use here) or $`1/N`$ in a large $`N`$ expansion, where “$`N`$” is the spin degeneracy $`g`$ HAMMER00 . The effective action functional will depend on $`\lambda `$ but we treat $`\rho `$ and $`\lambda `$ as independent variables:

$$\stackrel{~}{\mathrm{\Gamma }}=\stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ].$$
(24)
However, the ground state expectation value $`\rho _g(𝐱)`$ will naturally depend on $`\lambda `$ as determined by

$$\frac{\delta \stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]}{\delta \rho (𝐱)}|_{\rho =\rho _g}=0.$$
(25)
That is, if $`\lambda `$ is changed, a different $`\rho _g`$ will be necessary to satisfy this equation.

The Legendre transformation defining $`\stackrel{~}{\mathrm{\Gamma }}`$ is

$$\stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]=\stackrel{~}{W}[J,\lambda ]J(\mathrm{𝟏})\rho (\mathrm{𝟏}),$$
(26)
where $`J`$ is a functional of $`\rho `$ and $`\lambda `$ as well, as dictated by

$$\frac{\delta \stackrel{~}{W}[J,\lambda ]}{\delta J(\mathrm{𝟏})}=\rho (\mathrm{𝟏}).$$
(27)
Here we have introduced the convenient shorthand notation

$$A(\mathrm{𝟏})B(\mathrm{𝟏})d^3𝐱_1A(𝐱_1)B(𝐱_1).$$
(28)
As $`\lambda `$ is changed, $`J`$ must be adjusted so that the same $`\rho `$ is obtained when taking this derivative. This is how the dependence of $`J`$ on $`\lambda `$ arises; clearly this dependence can become quite complicated.

The inversion method now proceeds by expanding each of the quantities that depend on $`\lambda `$ in Eq. (26) in a Taylor series in $`\lambda `$:

$`J[\rho ,\lambda ]`$ $`=`$ $`J_0[\rho ]+J_1[\rho ]+J_2[\rho ]+\mathrm{},`$ (29)
$`\stackrel{~}{W}[J,\lambda ]`$ $`=`$ $`\stackrel{~}{W}_0[J]+\stackrel{~}{W}_1[J]+\stackrel{~}{W}_2[J]+\mathrm{},`$ (30)
$`\stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]`$ $`=`$ $`\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]+\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]+\stackrel{~}{\mathrm{\Gamma }}_2[\rho ]+\mathrm{},`$ (31)
where, as advertised, the power of $`\lambda `$ associated with each function or functional is indicated by the subscript. We can substitute the expansion for $`J`$ into the expansion for $`W`$ and do a functional Taylor expansion of $`W[J]`$ about $`J_0`$; this makes the $`\lambda `$ dependence manifest. Equating equal powers of $`\lambda `$ gives ($`l=0,1,2,\mathrm{}`$)

$`\stackrel{~}{\mathrm{\Gamma }}_l[\rho ]`$ $`=`$ $`\stackrel{~}{W}_l[J_0]J_l(\mathrm{𝟏})\rho (\mathrm{𝟏})+{\displaystyle \underset{k=1}{\overset{l}{}}}{\displaystyle \frac{\delta \stackrel{~}{W}_{lk}[J_0]}{\delta J_0(\mathrm{𝟏})}}J_k(1)`$ (32)
$`+{\displaystyle \underset{m=2}{\overset{l}{}}}{\displaystyle \frac{1}{m!}}{\displaystyle \underset{k_1,\mathrm{},k_m1}{\overset{k_1+\mathrm{}+k_ml}{}}}{\displaystyle \frac{\delta ^m\stackrel{~}{W}_{l(k_1+\mathrm{}+k_m)}[J_0]}{\delta J_0(\mathrm{𝟏})\mathrm{}\delta J_0(𝐦)}}J_{k_1}(\mathrm{𝟏})\mathrm{}J_{k_m}(𝐦).`$
Since $`\rho `$ is independent of $`\lambda `$ OKUMURA96 , each $`J_k[\rho ]`$ follows from each $`\stackrel{~}{\mathrm{\Gamma }}_k`$:

$$J_k(\mathrm{𝟏})=\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_k[\rho ]}{\delta \rho (\mathrm{𝟏})}.$$
(33)
We reiterate that all of the $`J_l`$’s as defined here are functionals of $`\rho `$.

Starting with the zeroth order expression,

$$\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]=\stackrel{~}{W}_0[J_0]J_0(\mathrm{𝟏})\rho (\mathrm{𝟏}),$$
(34)
we take its functional derivative with respect to $`\rho `$:

$$\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_0[\rho ]}{\delta \rho (\mathrm{𝟏})}=J_0(\mathrm{𝟏})=\frac{\delta \stackrel{~}{W}_0[J_0]}{\delta J_0(\mathrm{𝟏}^{})}\frac{\delta J_0(\mathrm{𝟏}^{})}{\delta \rho (\mathrm{𝟏})}J_0(\mathrm{𝟏})\rho (\mathrm{𝟏}^{})\frac{\delta J_0(\mathrm{𝟏}^{})}{\delta \rho (\mathrm{𝟏})}.$$
(35)
Rearranging,

$$\left(\frac{\delta \stackrel{~}{W}_0[J_0]}{\delta J_0(\mathrm{𝟏}^{})}\rho (\mathrm{𝟏}^{})\right)\frac{\delta J_0(\mathrm{𝟏}^{})}{\delta \rho (\mathrm{𝟏})}=0,$$
(36)
which implies

$$\rho (\mathrm{𝟏})=\frac{\delta \stackrel{~}{W}_0[J_0]}{\delta J_0(\mathrm{𝟏})},$$
(37)
since the strict concavity of $`\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]`$ prohibits $`\delta J_0(\mathrm{𝟏}^{})/\delta \rho (\mathrm{𝟏})`$ from having zero eigenvalues.

Equation (37) says that $`J_0(𝐱)`$ is the source (or potential, see below) that generates the expectation value $`\rho `$ from the noninteracting system (that is, the system defined by $`\lambda =0`$, which includes the external potential and $`J_0(𝐱)`$ but no interactions). $`J_0(𝐱)`$ is not an arbitrary function, but the *particular* one that has this property. The existence of a $`J_0(𝐱)`$ with this property is the cornerstone of the Kohn-Sham formalism.

Equation (37) also implies that the second term in Eq. (32) cancels with the $`k=l`$ term of the first sum for all $`l>0`$, and thus $`\stackrel{~}{\mathrm{\Gamma }}_l`$ simplifies to

$`\stackrel{~}{\mathrm{\Gamma }}_l[\rho ]`$ $`=`$ $`\stackrel{~}{W}_l[J_0]\delta _{l,0}J_l(\mathrm{𝟏})\rho (\mathrm{𝟏})+{\displaystyle \underset{k=1}{\overset{l1}{}}}{\displaystyle \frac{\delta \stackrel{~}{W}_{lk}[J_0]}{\delta J_0(\mathrm{𝟏})}}J_k(\mathrm{𝟏})`$ (38)
$`+{\displaystyle \underset{m=2}{\overset{l}{}}}{\displaystyle \frac{1}{m!}}{\displaystyle \underset{k_1,\mathrm{},k_m1}{\overset{k_1+\mathrm{}+k_ml}{}}}{\displaystyle \frac{\delta ^m\stackrel{~}{W}_{l(k_1+\mathrm{}+k_m)}[J_0]}{\delta J_0(\mathrm{𝟏})\mathrm{}\delta J_0(𝐦)}}J_{k_1}(\mathrm{𝟏})\mathrm{}J_{k_m}(𝐦).`$
These equations allow us to build the $`\stackrel{~}{\mathrm{\Gamma }}_l`$’s recursively. Note that the $`\stackrel{~}{W}_k`$ functionals have the same diagrammatic expansion as in Fig. 1, but the fermion lines are evaluated with Kohn-Sham (KS) propagators (see below). For a given $`l`$, we only need $`\stackrel{~}{W}_k`$’s with $`k`$ less than or equal to $`l`$ and $`J_k`$’s with $`k`$ smaller than $`l`$ (which means lower-order $`\stackrel{~}{\mathrm{\Gamma }}_k`$’s). We will illustrate the procedure by constructing the first few orders.

Since the lowest-order term in $`\stackrel{~}{W}[J]`$ corresponds to the system without interactions between the fermions, we can write $`\stackrel{~}{W}_0[J]`$ explicitly by introducing normalized single-particle orbitals that satisfy the equation

$$\left(\frac{^2}{2M}+v(𝐱)J_0(𝐱)\right)\phi _i(𝐱)=\epsilon _i\phi _i(𝐱).$$
(39)
The index $`i`$ represents all quantum numbers except for the spin (we consider only spin-independent interactions here). $`\stackrel{~}{W}_0[J_0]`$ is then (minus) the sum of the single-particle eigenvalues up to the Fermi energy $`\epsilon _\mathrm{F}`$ (which is equal to the chemical potential)

$$\stackrel{~}{W}_0[J_0]=g\underset{\epsilon _i<\epsilon _\mathrm{F}}{}\epsilon _i.$$
(40)
\[Equation (40) can be derived by evaluating $`W_0[J_0]\mathrm{Tr}\mathrm{ln}(G_{\mathrm{ks}}^0)^1`$ using the Kohn-Sham Green’s function $`G_{\mathrm{ks}}^0`$ defined below in Eq. (48).\] In practice, $`\epsilon _\mathrm{F}`$ is determined by simply counting orbitals until the $`A`$ lowest are filled (accounting for the spin degeneracy $`g`$). Note that the $`\epsilon _i`$’s are functionals of $`J_0`$ through Eq. (39); using the normalization of the $`\phi _i`$’s, we find

$$\frac{\delta \epsilon _i}{\delta J_0(𝐲)}=\frac{\delta }{\delta J_0(𝐲)}d^3𝐱\phi _i^{}(𝐱)\left(\frac{^2}{2M}+v(𝐱)J_0(𝐱)\right)\phi _i(𝐱)=\phi _i^{}(𝐲)\phi _i(𝐲).$$
(41)
Equations (37) and (40) show that the density may be written as

$$\rho (𝐱)=g\underset{i}{\overset{\mathrm{occ}.}{}}\frac{\delta \epsilon _i}{\delta J_0(𝐲)}=g\underset{i}{\overset{\mathrm{occ}.}{}}\phi _i^{}(𝐱)\phi _i(𝐱),$$
(42)
where the sum is over occupied (“occ.”) states. Equation (42) corresponds to the famous result of Kohn and Sham, which gives the exact ground state density in terms of the orbitals of a non-interacting system KOHN65 .

With the above results, the lowest order effective action is

$$\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]=g\underset{i}{\overset{\mathrm{occ}.}{}}\epsilon _id^3𝐱J_0(𝐱)\rho (𝐱).$$
(43)
We can also use Eq. (39) to eliminate $`\epsilon _i`$ from Eq. (43) so that it reads

$$\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]=T_s[\rho ]d^3𝐱v(𝐱)\rho (𝐱),$$
(44)
where

$$T_s[\rho ]=g\underset{i}{\overset{\mathrm{occ}.}{}}d^3𝐱\phi _i^{}(𝐱)\left(\frac{^2}{2M}\right)\phi _i(𝐱)$$
(45)
is the total kinetic energy of the KS non-interacting system. If $`\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]`$ were given as an explicit functional of $`\rho `$, then $`J_0[\rho ]`$ could be determined by taking a functional derivative according to Eq. (35). However, taking the functional derivative of the expression in Eq. (43) merely reproduces the result of Eq. (37). Instead, we follow Ref. VALIEV97  and determine $`J_0`$ from the *interacting* effective action, which we now construct.

From Eq. (38) we can find $`\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]`$ since

$$\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]=\stackrel{~}{W}_1[J_0[\rho ]].$$
(46)
For the dilute Fermi system, this is easily calculated. We first introduce the Green’s function of the KS non-interacting system, $`G_{\mathrm{ks}}^0(x,x^{})`$, which satisfies

$$\left(i_t+\frac{^2}{2M}v(𝐱)+J_0(𝐱)\right)G_{\mathrm{ks}}^0(𝐱t,𝐱^{}t^{})=\delta ^3(𝐱𝐱^{})\delta (tt^{})$$
(47)
with finite density boundary conditions FETTER71 . The KS Green’s function has the usual spectral decomposition in terms of the orbitals of Eq. (39):

$$iG_{\mathrm{ks}}^0(𝐱t,𝐱^{}t^{})=\underset{i}{}\phi _i(𝐱)\phi _i^{}(𝐱^{})e^{i\epsilon _i(tt^{})}[\theta (tt^{})\theta (\epsilon _i\epsilon _\mathrm{F})\theta (t^{}t)\theta (\epsilon _\mathrm{F}\epsilon _i)].$$
(48)
The Feynman rules in position space follow conventionally from Eq. (1) NEGELE88 ; HAMMER00  and we have from Fig. 1(a) \[with fermion lines representing $`iG_{\mathrm{ks}}^0`$\]

$$\stackrel{~}{W}_1[J_0]=\frac{1}{2}g(g1)C_0d^3𝐱G_{\mathrm{ks}}^0(x,x^+)G_{\mathrm{ks}}^0(x,x^+),$$
(49)
where the right side is independent of $`x_0`$ by Eq. (48). The Green’s function with equal arguments can be directly expressed in terms of the density,

$$\rho (𝐱)=igG_{\mathrm{ks}}^0(x,x^+).$$
(50)
Using this result and Eq. (46), we have

$$\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]=\frac{1}{2}\frac{(g1)}{g}C_0d^3𝐱|\rho (𝐱)|^2.$$
(51)
Since the dependence on $`\rho (𝐱)`$ is explicit in Eq. (51), we can directly take the functional derivative with respect to $`\rho `$ to obtain

$$J_1(𝐱)=\frac{C_0(g1)}{g}\rho (𝐱).$$
(52)
Direct functional derivatives with respect to $`\rho `$ will not be possible at higher order. However, we can find functional derivatives with respect to $`J_0`$. An alternative path to $`J_1[\rho ]`$ from $`\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]`$ is

$$J_1(𝐱)=\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_1[\rho ]}{\delta \rho (𝐱)}=d^3𝐲\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_1[\rho ]}{\delta J_0(𝐲)}\frac{\delta J_0(𝐲)}{\delta \rho (𝐱)}=d^3𝐲D^1(𝐱,𝐲)\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐲)},$$
(53)
which defines the inverse “density-density” correlator

$$D^1(𝐱,𝐲)\frac{\delta J_0(𝐲)}{\delta \rho (𝐱)}=\left(\frac{\delta ^2\stackrel{~}{W}_0[J_0]}{\delta J_0(𝐱)\delta J_0(𝐲)}\right)^1.$$
(54)
We can find an expression for $`D^1`$ by taking the functional derivative with respect to $`J_0`$ in Eq. (49). Remembering that $`G_{\mathrm{ks}}^0`$ is a functional of $`J_0`$ it is straightforward to show that

$$\frac{\delta G_{\mathrm{ks}}^0(x_1,x_2)}{\delta J_0(𝐱)}=G_{\mathrm{ks}}^0(x_1,x)G_{\mathrm{ks}}^0(x,x_2)𝑑x_0.$$
(55)
Using Eq. (55), the derivative of $`\stackrel{~}{W}_1`$ becomes

$`{\displaystyle \frac{\delta \stackrel{~}{W}_1[J]}{\delta J_0(𝐱)}}`$ $`=`$ $`g(g1)C_0{\displaystyle d^4y𝑑x_0G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(x,y^+)G_{\mathrm{ks}}^0(y,y^+)}`$ (56)
$`=`$ $`i(g1)C_0{\displaystyle d^4y𝑑x_0G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(x,y^+)\rho (𝐲)}.`$
Comparing Eq. (56) to Eq. (53) we find

$$D^1(𝐱,𝐲)=\frac{i}{g}\left[𝑑y_0G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)\right]^1.$$
(57)
This correlator will appear in all higher-order contributions.

Having determined $`J_1[\rho ]`$, we can find $`\mathrm{\Gamma }_2[\rho ]`$ from

$`\stackrel{~}{\mathrm{\Gamma }}_2[\rho ]`$ $`=`$ $`\stackrel{~}{W}_2[J_0]+{\displaystyle d^3𝐱\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐱)}J_1(𝐱)}+{\displaystyle \frac{1}{2}}{\displaystyle d^3𝐱d^3𝐲\frac{\delta ^2\stackrel{~}{W}_0[J_0]}{\delta J_0(𝐱)\delta J_0(𝐲)}J_1(𝐱)J_1(𝐲)}`$ (58)
$`=`$ $`\stackrel{~}{W}_2[J_0]+{\displaystyle \frac{1}{2}}{\displaystyle d^3𝐱d^3𝐲\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐱)}D^1(𝐱,𝐲)\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐲)}}.`$
$`\stackrel{~}{W}_2[J_0]`$ is calculated from the graphs Figs. 1(b) and (c):

$`\stackrel{~}{W}_2[J_0]`$ $`=`$ $`ig(g1){\displaystyle \frac{C_0^2}{4}}{\displaystyle d^4xd^4yG_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(y,x)}`$ (59)
$`ig(g1)^2{\displaystyle \frac{C_0^2}{2}}{\displaystyle d^4xd^4yG_{\mathrm{ks}}^0(x,x^+)G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(y,y^+)}.`$
By using Eqs. (56) and (57), we can show explicitly that the second term in the expression for $`\stackrel{~}{\mathrm{\Gamma }}_2`$ exactly cancels the second term in $`\stackrel{~}{W}_2`$. First note that

$`{\displaystyle d^3𝐲D^1(𝐱,𝐲)\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐲)}}`$ $`=`$ $`iC_0(g1){\displaystyle d^3𝐲\left[𝑑y_0G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)\right]^1}`$ (60)
$`\times {\displaystyle d^4w𝑑y_0G_{\mathrm{ks}}^0(w,y)G_{\mathrm{ks}}^0(y,w^+)G_{\mathrm{ks}}^0(w,w^+)}`$
$`=`$ $`iC_0(g1)G_{\mathrm{ks}}(x,x^+).`$
We now see that

$`{\displaystyle \frac{1}{2}}{\displaystyle d^3𝐱d^3𝐲\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐱)}D^1(𝐱,𝐲)\frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(𝐲)}}`$ (61)
$`=ig(g1)^2{\displaystyle \frac{C_0^2}{2}}{\displaystyle d^4xd^4yG_{\mathrm{ks}}^0(x,x^+)G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(y,y^+)},`$
which is the negative of the second term in Eq. (59). This cancellation is shown in Fig. 2, where $`D^1`$ is represented by a double line. Thus, Eq. (58) reduces to

$$\stackrel{~}{\mathrm{\Gamma }}_2[\rho ]=ig(g1)\frac{C_0^2}{4}d^4xd^4yG_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(x,y)G_{\mathrm{ks}}^0(y,x)G_{\mathrm{ks}}^0(y,x).$$
(62)
The second term in $`W_2`$ corresponds to the “anomalous” graph (c) in Fig. 1. This graph vanishes identically in the uniform system at zero temperature. However, it does not vanish at finite temperature or in a finite system (at any temperature). Rather, the cancellation exhibited above is analogous to the Kohn-Luttinger-Ward theorem mentioned earlier. Similar cancellations, of the type illustrated diagrammatically in Fig. 3, completely eliminate contributions of $`D^1`$ to the effective action up to N3LO in the EFT expansion for short-range forces. This complete cancellation does *not* occur for long-range forces, or if the zero-range delta functions at the $`C_0`$ vertices are regulated by a cutoff rather than by dimensional regularization, as used here.

We will illustrate the construction of $`\stackrel{~}{\mathrm{\Gamma }}_3`$ graphically to indicate how the cancellations occur at NNLO. We start with the explicit functional expression:

$`\stackrel{~}{\mathrm{\Gamma }}_3`$ $`=`$ $`\stackrel{~}{W}_3[J_0]+{\displaystyle \frac{\delta \stackrel{~}{W}_2[J_0]}{\delta J_0(\mathrm{𝟏})}}J_1(\mathrm{𝟏})+{\displaystyle \frac{\delta \stackrel{~}{W}_1[J_0]}{\delta J_0(\mathrm{𝟏})}}J_2(\mathrm{𝟏})+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta ^2\stackrel{~}{W}_1[J_0]}{\delta J_0(\mathrm{𝟏})\delta J_0(\mathrm{𝟐})}}J_1(\mathrm{𝟏})J_1(\mathrm{𝟐})`$ (63)
$`+{\displaystyle \frac{\delta ^2\stackrel{~}{W}_0[J_0]}{\delta J_0(\mathrm{𝟏})\delta J_0(\mathrm{𝟐})}}J_1(\mathrm{𝟏})J_2(\mathrm{𝟐})+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\delta ^3\stackrel{~}{W}_1[J_0]}{\delta J_0(\mathrm{𝟏})\delta J_0(\mathrm{𝟐})\delta J_0(\mathrm{𝟑})}}J_1(\mathrm{𝟏})J_1(\mathrm{𝟐})J_1(\mathrm{𝟑}).`$
Figure 5 illustrates that the two terms with $`J_2`$ factors cancel with each other. Every functional differentiation with respect to $`J_0`$, which inserts a density operator, is represented by a small dot (the large dot is a $`C_0`$ vertex). $`\stackrel{~}{W}_1`$ corresponds to Fig. 1(a), so the functional derivative in Fig. 5(a) yields the diagram shown; an explicit diagram for $`J_2`$ is not needed. Similarly, the second term is represented in Fig. 5(b). As shown, substituting $`J_1`$ (see Fig. 7) yields minus the first contribution, so the sum is zero and $`J_2`$ does not appear in $`\stackrel{~}{\mathrm{\Gamma }}_3`$.

Thus, $`\stackrel{~}{\mathrm{\Gamma }}_3`$ reduces to

$`\stackrel{~}{\mathrm{\Gamma }}_3`$ $`=`$ $`\stackrel{~}{W}_3[J_0]+{\displaystyle \frac{\delta \stackrel{~}{W}_2[J_0]}{\delta J_0(\mathrm{𝟏})}}J_1(\mathrm{𝟏})+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta ^2\stackrel{~}{W}_1[J_0]}{\delta J_0(\mathrm{𝟏})\delta J_0(\mathrm{𝟐})}}J_1(\mathrm{𝟏})J_1(\mathrm{𝟐})`$ (64)
$`+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\delta ^3\stackrel{~}{W}_1[J_0]}{\delta J_0(\mathrm{𝟏})\delta J_0(\mathrm{𝟐})\delta J_0(\mathrm{𝟑})}}J_1(\mathrm{𝟏})J_1(\mathrm{𝟐})J_1(\mathrm{𝟑}).`$
$`\stackrel{~}{W}_3`$ is given by the sum of Hugenholtz graphs in Fig. 1(d) through (j). But the sum of the other terms in Eq. (64) exactly cancel the diagrams in Fig. 1(d), (e), and (f), as shown in Fig. 6. The factors in front of the Feynman diagrams indicate additional multiplicative factors beyond those prescribed by the Feynman rules. These factors conspire to precisely subtract the anomalous diagrams in $`\stackrel{~}{W}_3`$, leaving only Fig. 1(g) through (j).

All higher orders in $`\stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]`$ are determined in a similar manner. Direct calculation becomes cumbersome, but one can formulate Feynman rules for $`\stackrel{~}{\mathrm{\Gamma }}`$, which dictate how to make appropriate insertions of $`D^1`$; these are given in Refs. FUKUDA95 ; YOKOJIMA95 ; OKUMURA96 ; VALIEV97 . It is important to note that the Kohn-Sham potential $`J_0`$ completely determines each order in the expansion of $`\stackrel{~}{\mathrm{\Gamma }}`$.

The cancellations exhibited here at low orders are expected from a comparison of the DFT/EFT expansion to perturbation theory. In perturbation theory for the confined system, the energy is calculated by evaluating the diagrams in Fig. 1 using the noninteracting propagator in the presence of the external potential for the fermion lines. (This propagator would take the form of Eq. (48) with harmonic oscillator wave functions for the orbitals.) Using the self-consistent $`G_{\mathrm{ks}}`$ instead sums an infinite class of higher-order diagrams at each order. The cancellation of anomalous diagrams from $`\stackrel{~}{\mathrm{\Gamma }}`$ corresponds to the removal of contributions already included through self-consistency (e.g., tadpoles). Although the nonperturbative contributions for the dilute, natural system are not required by power counting, the self-consistent calculation of the energy and density together is actually easier in practice than the purely perturbative calculation.

The appearance of the inverse density-density correlator $`D^1`$ can be understood by comparison to the effective actions for local fields and non-local composite fields. In the former case, the Legendre transformation removes one-particle intermediate states (leaving only one-particle-irreducible diagrams), while in the latter case, the Legendre transformation removes two-particle intermediate states (leaving two-particle-irreducible diagrams). Thus we infer that the role of $`D^1`$ is to remove intermediate states created by $`\psi ^{}\psi `$. The difference here is that we cannot write a closed-form expression for the effective action, which is possible in the other cases. However, we have seen that for short-range interactions, the extra diagrams at low orders cancel against anomalous diagrams, which is a great simplification.

To find an expression for $`J_0`$, we apply the variational principle satisfied by $`\stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]`$ to its expansion:

$`{\displaystyle \frac{\delta \stackrel{~}{\mathrm{\Gamma }}[\rho ,\lambda ]}{\delta \rho (\mathrm{𝟏})}}|_{\rho =\rho _{\mathrm{gs}}}`$ $`=`$ $`0`$ (65)
$`=`$ $`{\displaystyle \frac{\delta (\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]+\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ])}{\delta \rho (\mathrm{𝟏})}}|_{\rho =\rho _{\mathrm{gs}}}=J_0(\mathrm{𝟏})|_{\rho _{\mathrm{gs}}}+{\displaystyle \frac{\delta \stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]}{\delta \rho (\mathrm{𝟏})}}|_{\rho =\rho _{\mathrm{gs}}},`$
where the interaction effective action is

$$\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]\underset{i=1}{}\stackrel{~}{\mathrm{\Gamma }}_i[\rho ],$$
(66)
or

$$J_0(\mathrm{𝟏})|_{\rho =\rho _{\mathrm{gs}}}=\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]}{\delta \rho (\mathrm{𝟏})}|_{\rho =\rho _{\mathrm{gs}}}=D^1(\mathrm{𝟏},\mathrm{𝟏}^{})\frac{\delta \stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]}{\delta J_0(\mathrm{𝟏}^{})}|_{\rho =\rho _{\mathrm{gs}}}.$$
(67)
We stress that these relations hold only when we are solving for the Kohn-Sham potential corresponding to the ground-state density $`\rho _{\mathrm{gs}}(𝐱)`$.

The second equality in Eq. (67), which is shown diagrammatically through NLO in Fig. 7, is the key to the general Kohn-Sham self-consistent procedure:

1. Choose an approximation for $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]`$ by truncating the expansion in $`\lambda `$ at some order.
2. Make a reasonable initial guess for the Kohn-Sham potential $`J_0(𝐱)`$.
3. Calculate $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}[\rho ]`$ starting from $`J_0`$.
4. Use Eq. (67) to determine a new Kohn-Sham potential $`J_0(𝐱)`$.
5. Repeat the last two steps until self-consistency is reached (i.e., until some measure of the change in $`J_0`$ is less than a given tolerance).

In the next section, we will apply an LDA to $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{int}}`$, which means that it will given explicitly as a functional of $`\rho `$. In this case, the second equality in Eq. (67) is superfluous. One can also avoid explicitly evaluating Eq. (67) by adjusting $`J_0(𝐱)`$ using a steepest-descent approach VALIEV97 .

To this point we have neglected to mention that expressions such as $`\stackrel{~}{\mathrm{\Gamma }}_2`$ in Eq. (62) are divergent. For a uniform system, $`J_0`$ is constant. In this case, $`\stackrel{~}{W}_2`$ and subsequently $`\stackrel{~}{\mathrm{\Gamma }}_2`$ are renormalized by dimensional regularization and minimal subtraction as in Ref. HAMMER00 . In a finite system, with $`J_0`$ a function of $`𝐱`$, the ultraviolet linear divergence in Eq. (62) is renormalized by the same counterterm COLLINS86 , but it is computationally awkward to renormalize in the finite system. By using a derivative expansion, we can perform all renormalizations in the uniform system. This will be carried out explicitly in future work. With the LDA truncation applied below, we can simply use the renormalized expressions for the energy density from Ref. HAMMER00  to calculate the renormalized Kohn-Sham potential and energy functional.

## IV Results for Dilute Fermi System in a Trap

In this section, we present representative numerical results for the dilute Fermi system defined in Sect. II when confined in a harmonic oscillator trap. Our principal goal here is to illustrate the nature of the convergence of the EFT in a finite system both qualitatively and quantitatively. We consider effective range parameters corresponding to both attractive and repulsive underlying interactions (but we do not allow for pairing). To emphasize the difference between Kohn-Sham (KS) and Thomas-Fermi (TF) approaches and since we are ultimately interested in nuclear systems, we consider relatively small numbers of trapped fermions. Current experiments with trapped atoms use $`10^5`$$`10^6`$ atoms GRANADE02 ; as noted below, for such large systems the differences between KS and TF get washed out.

### IV.1 The Local Density Approximation (LDA)

To carry out the Kohn-Sham self-consistent procedure, we need to evaluate the expressions in the expansion of $`\stackrel{~}{\mathrm{\Gamma }}[\rho ]`$, so that we can use Eq. (67) to find $`J_0(𝐱)`$. In applications to Coulomb systems, the Kohn-Sham energy functional is conventionally written as

$$E[\rho ]=T_s[\rho ]+d^3𝐱v(𝐱)\rho (𝐱)+E_\mathrm{H}[\rho ]+E_{xc}[\rho ],$$
(68)
where $`T_s`$ is given in Eq. (45), $`E_\mathrm{H}`$ is the Hartree energy and $`E_{xc}[\rho ]`$ is known as the exchange-correlation energy. The Hartree energy is singled out because it can be explicitly written in terms of the density. When we have contact interactions only, the Fock term has the same dependence on the density as the Hartree term, and so we can also include it explicitly and replace $`E_\mathrm{H}`$ by $`E_{\mathrm{HF}}`$, redefining $`E_{xc}`$ appropriately.

This decomposition corresponds to writing the effective action as VALIEV97

$$\stackrel{~}{\mathrm{\Gamma }}[\rho ]=\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]+\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]+\underset{i=2}{\overset{\mathrm{}}{}}\stackrel{~}{\mathrm{\Gamma }}_i[\rho ],$$
(69)
up to an overall minus sign. The expression for $`\stackrel{~}{\mathrm{\Gamma }}_0`$ given in Eq. (44) shows that it has the same form as the first two terms in the energy functional in Eq. (68). In general, the explicit functional dependence of $`T_s`$ on the density is unknown since it enters implicitly through the orbitals. In practice it is easiest to calculate the first two terms by writing them as in Eq. (43):

$$\stackrel{~}{\mathrm{\Gamma }}_0[\rho ]=g\underset{i}{}\epsilon _id^3𝐱J_0(𝐱)\rho (𝐱).$$
(70)
The other terms in $`\stackrel{~}{\mathrm{\Gamma }}`$ may be written using the Green’s function for the KS non-interacting system as shown in Sec. III. The KS Green’s functions are explicit functionals of $`J_0`$ not $`\rho `$, however, and therefore almost all of $`\stackrel{~}{\mathrm{\Gamma }}`$ is not given as an explicit functional of the density; the general exception is the Hartree term and here the Fock term since we have only contact interactions. Furthermore, the actual expressions are quite difficult to evaluate in a finite system.

As a first approximation, we use the LDA, which may be considered the lowest-order term in a derivative expansion of the energy functional. The idea is to expand around the uniform system where (as was shown in Section II.2) the energy functional can be written as an explicit function of $`\rho `$. The LDA prescription is

$$E_{xc}^{\mathrm{LDA}}[\rho (𝐱)]d^3𝐱_{2+}(\rho _0)|_{\rho _0\rho (𝐱)},$$
(71)
where $`_{2+}`$ is the EFT energy density of the uniform system, including terms at second order and higher.

By combining results from Secs. II.2 and III, we may write the exchange-correlation energy, to third order (NNLO) in the EFT expansion, as

$`E_{xc}^{\mathrm{LDA}}[\rho (𝐱)]`$ $`=`$ $`{\displaystyle d^3𝐱\left\{_2(\rho (𝐱))+_3(\rho (𝐱))+\mathrm{}\right\}}`$ (72)
$`=`$ $`b_1{\displaystyle \frac{a_s^2}{2M}}{\displaystyle d^3𝐱[\rho (𝐱)]^{7/3}}`$
$`+\left(b_2a_s^2r_s+b_3a_p^3+b_4a_s^3\right){\displaystyle \frac{1}{2M}}{\displaystyle d^3𝐱[\rho (𝐱)]^{8/3}}+\mathrm{},`$
where the dimensionless $`b_i`$ are

$`b_1`$ $`=`$ $`{\displaystyle \frac{4}{35\pi ^2}}(g1)\left({\displaystyle \frac{6\pi ^2}{g}}\right)^{4/3}(112\mathrm{ln}2),`$
$`b_2`$ $`=`$ $`{\displaystyle \frac{1}{10\pi }}(g1)\left({\displaystyle \frac{6\pi ^2}{g}}\right)^{5/3},`$
$`b_3`$ $`=`$ $`{\displaystyle \frac{1}{5\pi }}(g+1)\left({\displaystyle \frac{6\pi ^2}{g}}\right)^{5/3},`$
$`b_4`$ $`=`$ $`\left({\displaystyle \frac{6\pi ^2}{g}}\right)^{5/3}\left(0.0755(g1)+0.0574(g1)(g3)\right).`$ (73)
In order to solve for the orbitals in Eq. (39) and to calculate the energy we need the expression for $`J_0(𝐱)`$. In the LDA, this is simple since

$$J_0(𝐱)=\frac{\delta }{\delta \rho (𝐱)}\left(\stackrel{~}{\mathrm{\Gamma }}_1[\rho ]+\underset{i=2}{\overset{\mathrm{}}{}}\stackrel{~}{\mathrm{\Gamma }}_i[\rho ]\right)=\frac{\delta }{\delta \rho (𝐱)}\left(E_{\mathrm{HF}}[\rho ]+E_{xc}[\rho ]\right).$$
(74)
To NNLO we find:

$$J_0(𝐱)=\frac{(g1)}{g}\frac{4\pi a_s}{M}\rho (𝐱)\frac{7}{3}b_1\frac{a_s^2}{2M}[\rho (𝐱)]^{4/3}\frac{8}{3}\left(b_2a_s^2r_s+b_3a_p^3+b_4a_s^3\right)\frac{1}{2M}[\rho (𝐱)]^{5/3}.$$
(75)
A convenient expression for the total binding energy (through NNLO) follows by substituting for $`J_0(𝐱)`$ and combining terms:

$`E[\rho (𝐱)]`$ $`=`$ $`g{\displaystyle \underset{i}{\overset{\mathrm{occ}.}{}}}\epsilon _i{\displaystyle \frac{1}{2}}{\displaystyle \frac{(g1)}{g}}{\displaystyle \frac{4\pi a_s}{M}}{\displaystyle d^3𝐱[\rho (𝐱)]^2}{\displaystyle \frac{4}{3}}b_1{\displaystyle \frac{a_s^2}{2M}}{\displaystyle d^3𝐱[\rho (𝐱)]^{7/3}}`$ (76)
$`{\displaystyle \frac{5}{3}}\left(b_2a_s^2r_s+b_3a_p^3+b_4a_s^3\right){\displaystyle \frac{1}{2M}}{\displaystyle d^3𝐱[\rho (𝐱)]^{8/3}}.`$
In the numerical calculations given below, the noninteracting case ($`C_0=0`$) uses $`J_0(𝐱)0`$ and the first term in Eq. (76), LO uses the first term in Eq. (75) and the first two terms in Eq. (76), and so on for NLO and NNLO.

### IV.2 Kohn-Sham Self-Consistent Procedure

Here we describe the numerical procedure used to find the Kohn-Sham orbitals. We assume closed shells, so the density and potentials are functions only of the radial coordinate $`r|𝐱|`$. (This restriction is straightforward to relax.) The Kohn-Sham iteration procedure is as follows:

1. Guess an initial density profile $`\rho (r)`$. If the system is particularly nonlinear this initial choice may be critical. More generally, there may be a metastable state (e.g., if the system ultimately collapses). For nuclear systems, the experience is that a crude caricature of the true density (such as a Woods-Saxon shape) is adequate. In the present case, the non-interacting harmonic oscillator density is sufficient.
2. Evaluate the local single-particle potential

$$v_s[\rho (r)]v_s(r)v(r)J_0(r)$$
(77)
using Eqs. (75) at the chosen level of approximation (e.g., NLO). Beyond the LDA $`v_s`$ could be a functional not just of the density but of the Kohn-Sham orbitals individually.
3. Solve the Schrödinger equation for the lowest $`A`$ states (including degeneracies), to find a set of orbitals and Kohn-Sham eigenvalues $`\{\phi _\alpha ,\epsilon _\alpha \}`$:

$$\left[\frac{^2}{2M}+v_s(r)\right]\phi _\alpha (𝐱)=\epsilon _\alpha \phi _\alpha (𝐱).$$
(78)
4. Compute a new density from the orbitals:

$$\rho (r)=\underset{\alpha =1}{\overset{A}{}}|\phi _\alpha (𝐱)|^2.$$
(79)
All other ground-state observables are functionals of $`\{\phi _\alpha ,\epsilon _\alpha \}`$.
5. Repeat steps 2.–4. until changes are acceptably small (“self-consistency”). In practice, the changes in the density are “damped” by using a weighted average of the densities from the $`(n1)`$th and $`n`$th iterations:

$$\rho (r)=\beta \rho _{n1}(r)+(1\beta )\rho _n(r),$$
(80)
with $`0<\beta 1`$.

This procedure has been implemented for dilute fermions in a trap using computer codes written both in C and in Mathematica. Two methods for carrying out step 3. were tested. The Kohn-Sham single-particle equations are solved in one approach by direct integration of the differential equations via the Numerov method BLATT67  and in the other approach by diagonalization of the single-particle Hamiltonian in a truncated basis of unperturbed harmonic oscillator wavefunctions. The same results are obtained to high accuracy. For closed shells, either method is efficient and easy to code.

### IV.3 Fermions in a Harmonic Trap

The interaction through NNLO is specified in terms of the three effective range parameters $`a_s`$, $`r_s`$, and $`a_p`$. For the numerical calculations presented here, we consider two cases in the dilute limit, $`a_s\{r_s,a_p\}0`$ with both signs for $`a_s`$, and hard sphere repulsion with radius $`R`$, in which case $`a_s=a_p=R`$ and $`r_s=2R/3`$.

Lengths are measured in units of the oscillator parameter $`b\sqrt{\mathrm{}/M\omega }`$, masses in terms of the fermion mass $`M`$, and $`\mathrm{}=1`$. In these units, $`\mathrm{}\omega `$ for the oscillator is unity and the Fermi energy of a non-interacting gas with filled shells up to $`N_F`$ is $`E_F=(N_F+3/2)`$. The total number of fermions $`A`$ is related to $`N_F`$ by

$$A=\frac{g}{6}(N_F+1)(N_F+2)(N_F+3).$$
(81)
Since we have only considered spin-independent interactions, our results are independent of whether the spin degeneracy $`g`$ actually originates from spin, isospin, or some flavor index.

With interactions included, single-particle states are labeled by a radial quantum number $`n`$, an orbital angular momentum $`l`$ with $`z`$-component $`m_l`$, and the spin projection. The radial functions depend only on $`n`$ and $`l`$, so the degeneracy of each level is $`g\times (2l+1)`$. Excluding spin, the solutions are of the form

$$\phi _{nlm_l}(𝐱)=\frac{u_{nl}(r)}{r}Y_{lm_l}(\mathrm{\Omega }),$$
(82)
where the radial function $`u_{nl}(r)`$ satisfies

$$\left[\frac{1}{2}\frac{d^2}{dr^2}+v_s(r)+\frac{l(l+1)}{2r^2}\right]u_{nl}(r)=\epsilon _{nl}u_{nl}(r).$$
(83)
The $`u_{nl}`$’s are normalized according to

$$_0^{\mathrm{}}|u_{nl}(r)|^2𝑑r=1.$$
(84)
Thus the density is given by:

$$\rho (r)=g\underset{i}{\overset{\mathrm{occ}.}{}}\frac{|u_i(r)|^2}{4\pi r^2}=g\underset{nl}{\overset{\mathrm{occ}.}{}}\frac{(2l+1)}{4\pi r^2}|u_{nl}(r)|^2$$
(85)
The interactions are sufficiently weak that the occupied states are in one-to-one correspondence with those occupied in the non-interacting harmonic oscillator potential.

To check our numerical calculations, we first compared to density distributions at zero temperature from Ref. BRUUN98 , which used a contact interaction with strength corresponding to $`a_s=0.16`$ in our units. In that work, non-interacting and mean-field (Hartree-Fock) results for the normal state (for possible comparison to superfluid solutions) were presented for systems with 240 and 330 atoms. The corresponding densities in our EFT expansion are the curves labeled “$`C_0=0`$ exact” and “Kohn-Sham LO” in Figs. 9 and 13. We also show densities for the same systems but with $`a_s=+0.16`$ in Figs. 9 and 13. As one might expect, the attractive interaction pulls the density in while the repulsive interaction pushes it out relative to the non-interacting density. Values for the energy per particle, the average Fermi momentum, and the rms radius for each calculation are given in Table 1. Averages are defined as

$$f(𝐱)\frac{1}{A}d^3𝐱f(𝐱)\rho (𝐱).$$
(86)
In each of these figures, the other effective range parameters, $`r_s`$ and $`a_p`$, are set to zero. We have also shown results for a representative system with spin degeneracy $`g=4`$ with weaker attractive ($`a_s=0.10`$) and repulsive ($`a_s=+0.10`$) interactions in Figs. 11 and 11. To illustrate the impact of $`r_s`$ and $`a_p`$ on the NNLO results, we have also included two calculations where the underlying interaction is hard-sphere repulsion with $`R=0.16`$ ($`g=2`$) and $`R=0.10`$ ($`g=4`$) in Figs. 15 and 15.

From Figs. 9 and 9, we see what seems to be good convergence of the density by NNLO. Note that the LO results are not well converged in either case. From Table 1, we find that the average Fermi momentum

$$k_\mathrm{F}\frac{1}{A}d^3𝐱\left[\frac{6\pi ^2\rho (𝐱)}{g}\right]^{1/3}\rho (𝐱)$$
(87)
is such that the expansion parameter in an LDA sense, $`k_\mathrm{F}a_s`$, is equal to or greater than one-half, which implies poor convergence. In fact, the convergence is misleading, as seen by comparison to the hard-sphere repulsive case in Fig. 15. (See the discussion of energy convergence in Sec. IV.5.)

The systems with $`g=4`$ have $`k_\mathrm{F}a_s1/4\text{}1/3`$, and the convergence of even the hard-sphere case is good. We note that the instability for large $`g`$ discussed in Ref. Furnstahl:2002gt  is manifested as a collapse of the density iteration by iteration when the stability criterion is violated by the interior density.

### IV.4 Comparison to Thomas-Fermi

In this section, we compare the LDA Kohn-Sham (KS) results to those from a Thomas-Fermi (TF) approximation. By Thomas-Fermi, we mean that the entire energy, including the kinetic energy, is calculated in a local density approximation. In particular, the Thomas-Fermi kinetic energy functional is ARGAMAN00

$$T_{TF}[\rho ]=\frac{3}{10M}\left(\frac{6\pi ^2}{g}\right)^{2/3}d^3r[\rho (𝐫)]^{5/3},$$
(88)
rather than being computed as the sum of the kinetic energy of Kohn-Sham orbitals. We define the Thomas-Fermi potential energy functional to have the same form as the Kohn-Sham potential energy functional in the LDA. Historically, the Thomas-Fermi approach, as applied to atoms and molecules, was the first attempt to use the (electron) density as the basic variable rather than solving for the wavefunction. In its original form, only the Hartree term and the external nuclear-electron attractive potential were included in the potential energy. We will generalize to define LO, NLO, and NNLO Thomas-Fermi approximations.

The idea behind Thomas-Fermi is that in each volume element $`dV`$ we have chemical equilibrium, with chemical potential $`\mu `$. Each volume element is labeled by $`r`$ (we assume spherical symmetry for convenience). Local eigenvalues $`E_k(r)`$ are computed at each $`r`$ as if in a uniform system:

$$E_k(r)=\frac{k^2}{2M}+v_s(r),$$
(89)
and levels are filled until

$$E_{k_\mathrm{F}}(r)=\mu ,$$
(90)
which defines the local Fermi momentum $`k_\mathrm{F}(r)`$. This is turn defines the Thomas-Fermi density $`\rho _{TF}(r,\mu )`$ for a given chemical potential $`\mu `$ as:

$$\rho _{TF}(r,\mu )=\{\begin{array}{cc}\frac{g}{6\pi ^2}\left(2M[\mu v_s(r)]\right)^{3/2}& \text{if }\mu >v_s(r)\hfill \\ 0& \text{if }\mu v_s(r)\hfill \end{array},$$
(91)
where the potential $`v_s(r)`$ includes the external trap potential and the LDA Kohn-Sham potential \[see Eq. (77)\]. The conventional TF procedure combines this equation with an equation for the potential (e.g., a Poisson equation in the Coulomb case) by substituting for the potential. Here we solve it in a two-step process closely analogous to the Kohn-Sham solution procedure.

For a given choice of $`\mu `$, the number of fermions $`A_{TF}`$ is given by:

$$A_{TF}(\mu )=4\pi _0^{\mathrm{}}r^2\rho _{TF}(r,\mu )𝑑r,$$
(92)
where we’ve assumed a spherically symmetric distribution. We find the correct value of $`\mu `$ for an input value of $`A`$ using a root finding program, which finds the zero of

$$f(\mu )AA_{TF}(\mu )$$
(93)
in the interval $`\mu _{\mathrm{min}}<\mu <\mu _{\mathrm{max}}`$. The procedure to find the self-consistent $`\rho _{TF}`$ is iterative:

1. Guess an initial $`\rho _{TF}(r)`$ (for example, the unperturbed density);
2. given $`\rho _{TF}(r)`$, compute $`v_s(r)`$ for the noninteracting, LO, NLO, or NNLO case;
3. using this $`v_s(r)`$ in Eq. (91), find $`\mu `$ so that $`f(\mu )=0`$;
4. with the new value of $`\mu `$, calculate a new $`\rho _{TF}(r,\mu )`$ from Eq. (91);
5. return to step 2., continuing until $`\rho _{TF}`$ and $`\mu `$ do not change within a prescribed tolerance.

This procedure is rather simple and scales very well with the number of particles. There are deficiencies to the Thomas-Fermi approach, however, for many problems of interest. These deficiencies are most evident for the Coulomb problem, where it is proved that molecules do not bind (the separate atoms always have lower energy) TELLER62 . Gradients expansions can improve the approximation but have not had a quantitative impact compared to Kohn-Sham approaches.

The kinetic energy contribution is a leading source of nonlocality in the energy functional. One of the chief virtues of the Kohn-Sham approach is that it treats this component much more effectively than a derivative expansion would. The comparison with Thomas-Fermi calculations highlights this difference. The most visible consequences are the absence of shell structure in ground state densities in the Thomas-Fermi approach and the incorrect treatment of the low-density asymptotic region (dominated by the last Kohn-Sham orbitals). In contrast, the last Kohn-Sham orbital has the correct energy (the ionization energies in the exact and Kohn-Sham systems are equal), so the exponential tail of the distribution is correct. These deficiencies of Thomas-Fermi are relevant for the calculation of nuclear charge and matter densities.

In Figs. 17 and 17, Thomas-Fermi and Kohn-Sham densities are compared for a small system of 20 atoms (the other parameters are given in the figures). Both non-interacting and NNLO curves are shown in each figure. The shell structure is evident in the non-interacting density and is only slightly damped by the interactions in the Kohn-Sham approach. In contrast, the Thomas-Fermi curves are featureless. The same comparisons but with an order of magnitude more atoms ($`A=240`$) are shown in Figs. 19 and 19. In these case the differences are much smaller. For a small number of particles, the Kohn-Sham procedure is a comparable computation to Thomas-Fermi. With thousands of atoms, the Thomas-Fermi approximation should be accurate and efficient.

In Table 2, the energies per particle and other properties are compared for Kohn-Sham and Thomas-Fermi solutions. With $`A=20`$ atoms, the energy is reproduced at about the 2% level and the average Fermi momentum at the 1% level. With $`A=240`$ atoms, the energy is reproduced better than 1% and the average Fermi momentum to a third of a percent. Thus the Thomas-Fermi density will be quite adequate in making error estimates, which we turn to next.

### IV.5 Power Counting and Convergence

A major motivation for the use of effective field theory for many-body problems is the promise of error estimates. We face two challenges in making such estimates for a Kohn-Sham DFT. There are usually new low-energy constants (LEC’s) at each successive order in the EFT expansion. We first need to estimate the size of the LEC’s in the omitted orders. Second, we need to estimate their numerical impact on the functional and subsequently on observables in finite systems. In the present case, where we have a perturbative low-density expansion, we can do both directly.

Naive dimensional analysis, or NDA, is frequently used to estimate unknown LEC’s in an EFT Lagrangian. The idea is to identify the relevant dimensional momentum and mass scales in the problem and to rescale each term in the Lagrangian as an appropriate combination of scales times a dimensionless coefficient. An estimate of the truncation error follows by assuming the dimensionless coefficient is order unity (e.g., from 1/3 to 3).

In Ref. FURNSTAHL99b , NDA appropriate to low-energy chiral effective field theories of QCD was applied to covariant energy functionals for nuclei. While these functionals were viewed as approximate Kohn-Sham energy functionals, their form was directly derived from a chiral Lagrangian in the Hartree approximation. That is, there was a one-to-one correspondence between terms in the Lagrangian and terms in the corresponding energy functional. As a result, estimates of coefficients in the Lagrangian were immediately translated into error estimates in the functional, which led to estimates in specific finite nuclei through local density approximations FURNSTAHL99b .

In the present case, the Hartree (actually Hartree-Fock) terms again provide an immediate connection between NDA estimates in the Lagrangian and estimates of the energy per particle in terms of the fermion density. Using LDA estimates of the average density, we get energy estimates for specific systems of trapped atoms. EFT power counting then provides estimates for the higher-order terms using the normalization established by the Hartree terms. For a purely short-ranged interaction, the NDA is simple: the estimate of the Hartree-Fock energy contribution from a given term in the Lagrangian is found by replacing $`\psi ^{}\psi `$ by the average density (and including an appropriate spin factor). Any reasonable density (e.g., the Thomas-Fermi result) can be used to find the average density, since we already allow for a much larger uncertainty in the coefficients. Local density estimates for derivatives of densities are also sufficiently accurate.

In Fig. 20, the contributions to the energy per particle from LO, NLO, and NNLO terms in the energy per particle are shown as round symbols for three of the systems from Table 1. So, for example, the NLO contribution is from $`|E_{\mathrm{NLO}}E_{\mathrm{LO}}|/A`$. The square symbol denote estimates based on naive dimensional analysis, with error bars indicating a 1/2 to 2 uncertainty in the coefficients. In particular, the LO contribution is an estimate of the Hartree-Fock term and the NLO and NNLO estimates were found by multiplying the LO result by $`k_\mathrm{F}a_s`$ and $`(k_\mathrm{F}a_s)^2`$, respectively. With one exception, the NDA estimates are good predictors of the actual contribution. Therefore, we can use these estimates to predict the uncertainty in the energy per particle from higher orders. The exception is the NNLO estimate for the $`g=2`$ system, which greatly overestimates the actual contribution at that order. The reason is evident from the $`b_4`$ coefficient in Eq. (73), which determines the size of this contribution. Each of the two terms in $`b_4`$ are the size expected from NDA, but for $`g=2`$ they happen to largely cancel. There is always this possibility of unnaturally small coefficients (and subsequent contributions) because of accidental cancellations. If we compare instead the NNLO energy from the hard-sphere-repulsion case, which has additional contributions at NNLO, we find that it is close to the NDA estimate.

## V Summary and Conclusions

In this paper, we construct a Kohn-Sham density functional for a confined, dilute Fermi gas as an order-by-order effective field theory (EFT) expansion. The starting point is a generating functional with a source $`J(x)`$ coupled to the composite density operator $`\psi ^{}\psi `$. A functional Legendre transformation with respect to the source leads to an effective action of the density. In general we might expect complications with such an effective action because of the renormalization of composite operators YOKOJIMA95 ; OKUMURA96 , but these are avoided in the present case because the fermion number is a conserved charge COLLINS86 .

A conventional effective action is a functional of the expectation value of elementary fields in the Lagrangian. It defines a classical field theory that contains all of the quantum effects of the interacting quantum field theory associated with the original Lagrangian. That is, one reproduces results of the full field theory from tree level calculations based on the propagators and vertices of the effective action COLEMAN88 ; PESKIN95 ; WEINBERG96 . In the present case, the effective action of the density can be used to calculate the ground state energy, including all correlations, with what looks like a Hartree calculation (which is the many-body equivalent of tree level). This is density functional theory.

The calculation is carried out by adapting the inversion method proposed in Refs. FUKUDA94 ; FUKUDA95  and the Kohn-Sham procedure of Ref. VALIEV97  to an EFT treatment of the dilute Fermi gas. Instead of organizing the perturbative inversion in terms of a coupling constant (e.g., the electron charge squared), we use the EFT expansion parameter for the natural dilute system, which is $`1/\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is the resolution scale. In the uniform system, the ultimate dimensionless expansion parameters are products of the Fermi momentum $`k_\mathrm{F}`$ and parameters of the effective range expansion (e.g., the scattering length $`a_s`$, which is of order $`1/\mathrm{\Lambda }`$ for a natural system). In a finite system, the density-weighted average of a local Fermi momentum times the effective range parameters controls one expansion, with another expansion involving gradients of the density.

In the present work, we used the local density approximation (LDA), which means that only the first expansion was tested. The observed convergence of the density and energy in sample systems confirms this expansion for a finite system. An error plot of contributions to the energy per particle versus the order of the calculation shows that we can reliably estimate the truncation error in a finite system. The next step is to develop a derivative expansion and test it for convergence. Work is in progress to adapt to the present problem the methods that have been used to derive derivative expansions for one-loop effective actions.

The EFT expansion for the uniform system using dimensional regularization and minimal subtraction (DR/MS) offered many simplifications over conventional treatments HAMMER00 . The derivative expansion approach would preserve all of these advantages. In particular, renormalization of the derivative expansion coefficient functions would be carried out entirely in the uniform system using DR/MS, so that a given diagram contributes to only a single prescribed order in the expansion.

An alternative to the inversion method that might be preferred for some systems is to introduce an auxiliary field (or fields) and to carry out the Legendre transformation conventionally with respect to that field FUKUDA94 ; FUKUDA95 . This could be applied to the large-N EFT expansion discussed in Ref. Furnstahl:2002gt . The construction of Kohn-Sham DFT in the auxiliary field formulation is discussed in Refs. FAUSSURIER00 ; VALIEV96 .

Many of the phenomena of greatest interest in experiments on trapped fermion atoms involve tuning the system so that the $`s`$-wave scattering length is large. For example, one would like to study superfluid transitions at low temperature GRANADE02 ; HOLLAND01 ; TIMMERMANS01 . systematic solution to the large scattering length problem at finite density is not yet available, even for a uniform system. Furthermore, with *any* attractive interaction we expect a pairing instability. Extending the DFT/EFT expansion procedure to include pairing will be an important challenge. The inversion method has been applied to conventional BCS superconductivity in Ref. INAGAKI92  using a source coupled to the pair creation and destruction operator. Work to adapt this procedure to the EFT is in progress.

These extensions will be directly relevant for nuclear applications as well. In addition, to adapt the density functional procedure to chiral effective field theories with explicit pions, we will need to extend the discussion to include long-range forces (“long range” means compared to $`1/\mathrm{\Lambda }`$). Other extensions include DFT with alternative source terms (e.g., spin-density DFT) and time-dependent DFT in the EFT formalism. Studies of each of these extensions are in progress.

###### Acknowledgements.

We thank P. Bedaque, H.-W. Hammer, A. Schwenk and B.D. Serot for useful comments. This work was supported in part by the National Science Foundation under Grant No. PHY–0098645.


format: text

text:
# Chapter 1 Investigating Galaxy Clusters through 𝛾-ray Emission

## Chapter 1 Investigating Galaxy Clusters through $`\gamma `$-ray Emission

### Abstract

We address the role of $`\gamma `$-ray astronomy in the investigation of nonthermal processes in the large scale structure of the universe. Based on EGRET upper limits on nearby galaxy clusters (GCs) we constrain the acceleration efficiency of CR electrons at intergalactic shocks to $`1\%`$ than the shock ram pressure. That implies a contribution to the cosmic $`\gamma `$-ray background from intergalactic shocks of order 25 % of the measured level. We model spatial and spectral properties of nonthermal $`\gamma `$-ray emission due to shock accelerated cosmic-rays (CRs) in GCs and emphasize the importance of imaging capability of upcoming $`\gamma `$-ray facilities for a correct interpretation of the observational results. GC observations at this photon energy will help us understand the origin of the radio emitting particles, the possible level of CR pressure and the strength of magnetic fields in intracluster environment and possibly will allow us detect the accretion shocks.

### 1. Introduction

Cosmic shocks emerge during structure formation in the universe, due to gravitationally driven supersonic gas infall onto collapsing objects. Astrophysical shocks are collisionless and capable of accelerating cosmic-rays (CRs) via first order Fermi mechanism . CR electrons suffer severe energy losses dominated, in typical intra-cluster environment, by inverse Compton (IC) for energies above $`150`$ MeV and Coulomb losses below that. Low energy, sub-relativistic CR protons also have short lifetimes due to Coulomb losses. However, for relativistic protons the dominant energy loss mechanism up to the energy threshold for photo-pair production is p-p inelastic collisions with a timescale longer than a Hubble time. Thus, once accelerated CR protons accumulate within large scale structures (LSSs) confined there by $`\mu `$G strong turbulent magnetic fields .

It is of interest to investigate CR acceleration at LSS shocks for several reasons. Firstly, GCs exhibit non-thermal radiation. This mainly consists of diffuse radio emission extending on Mpc scales. It is thought to be synchrotron radiation but there is no consensus as to the origin of the emitting electrons. Also, excess of radiation with respect to thermal emission has been reported at hard X-ray  and, more controversially, extreme UV photon energies . Secondly, if shocks acceleration operates efficiently, the proton component could bear a significant fraction of the total gas pressure, affecting the dynamics of cosmic structures . Finally, CR electrons accelerated at inter-galactic shocks could contribute to the cosmic $`\gamma `$-ray background (CGB) . This contribution addresses the role of $`\gamma `$-ray observations for investigating CR acceleration at LSS shocks. We describe the method in §LABEL:fm:num.sec, the results in LABEL:fm:cgb.sec-LABEL:fm:coma.sec and conclude in §LABEL:fm:disc.sec

### 2. Numerical Simulation

We a perform numerical simulation that models simultaneously the formation of the LSS and the evolution of three CR species: namely shock accelerated protons and (primary) electrons and secondary e± generated in p-p inelastic collisions of the simulated CR protons and the thermal gas. For the LSS we adopt a canonical, flat $`\mathrm{\Lambda }`$CDM model with total mass density $`\mathrm{\Omega }_m=0.3`$, vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, Hubble constant $`H_0=67`$ km s-1 Mpc-1, baryonic mass density, $`\mathrm{\Omega }_b=0.04`$, and spectral index and cluster normalization for initial density perturbations, $`n_s=1`$ and $`\sigma _8=0.9`$ respectively (see ref.  for full details).

The CR dynamics, includes: injection at shocks, acceleration, energy losses and spatial transport. It is computed numerically through the code COSMOCR . According to the adopted thermal leakage injection scheme, a fraction about $`10^4`$ of the protons crossing the shock is injected as CRs. For shocks with Mach number, $`M=410`$, responsible for most of the heating of the intergalactic gas , this roughly corresponds to converting 30% of the shock ram pressure into CR proton pressure. Further, we assume that the ratio of accelerated primary CR electrons to protons at relativistic energies is given by a parameter $`R_{e/p}10^2`$ . In accord to the test particle limit of diffusive shock acceleration theory , the accelerated particles are distributed in momentum as a power-law with slope $`q=4/(1M^2)`$. All relevant loss mechanisms are accounted for . We follow CR protons in the momentum range between $`0.1`$ GeV/c and $`10^6`$ GeV/c, and CR electrons and e± between 15 MeV/c and 20 TeV/c.

### 3. Contribution to the Cosmic Gamma-ray Background

Fig. LABEL:fm:gamma.fig (left) shows the contribution to the CGB from CRs accelerated at LSS shocks . It is computed as

$$\epsilon ^2J(\epsilon )=\epsilon \frac{c}{4\pi H_0}_0^{z_{max}}\frac{e^{\tau _{\gamma \gamma }}}{[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}}\frac{j[\epsilon (1+z),z]}{(1+z)^4}𝑑z$$
(1)
where $`\epsilon `$ is the photon energy, $`j(\epsilon ,z)`$ is the computational-box-averaged spectral emissivity in units ‘photons cm-3 s-1’ at red-shift, $`z`$, $`\tau _{\gamma \gamma }`$ is an attenuation factor due to photo-pair creation, $`\gamma \gamma e^\pm `$, $`c`$ is the speed of light and $`z_{max}`$ an upper limit of integration. EGRET observational data (solid dots)  are also shown for comparison. We consider the following emission processes: IC emission of CR electrons scattering off cosmic microwave background photons (dot line), decay of neutral pions ($`\pi ^0\gamma \gamma `$) produced in p-p inelastic collisions (dash line) and IC emission from secondary e± (dot-dash line). In fig. LABEL:fm:gamma.fig the total flux (solid line) corresponds roughly to a constant value at the level of 0.2 keV cm-2 s-1 sr-1 throughout the spectrum. It is dominated by IC emission from primary electrons. Fractions of order 30% and 10 % are produced by $`\pi ^0`$-decay and IC emission from secondary e± respectively.

All three components produce the same flat spectrum, consistent with the observations. The result is a reflection of the fact that the CRs distributions were generated in strong shocks . However, the computed flux is only $``$ 15 % of the observed CGB. It is difficult to imagine a higher contribution from $`\pi ^0`$-decay and IC emission from e± . In fact, if more CR protons were produced at shocks, CR-induced shock modifications would actually reduce the population of $`\gamma `$-ray emitting protons (and e± ). On the other hand, the fraction, $`\eta `$, of shock ram pressure converted into CR electrons, can be constrained by comparing the simulated clusters’ $`\gamma `$-ray photon luminosity above 100 MeV to the upper limits set by the EGRET  for nearby GCs. This is done in fig. LABEL:fm:gamma.fig (right panel). The simulation data are best-fit by the curve (solid line):

$$\mathrm{L}_\gamma (>100\mathrm{MeV})=8.7\times 10^{43}\left(\frac{\eta }{4\times 10^3}\right)\left(\frac{\mathrm{T}_\mathrm{x}}{\mathrm{keV}}\right)^{2.6}\mathrm{ph}\mathrm{s}^1.$$
(2)
According to our calculations , the EGRET upper limits require that $`\eta 0.8\%`$. This implies an upper limit on the computed $`\gamma `$-ray flux of about 0.35 keV cm-2 s-1 sr-1 or a fraction of order $`25`$ % of the CGB.

### 4. Synthetic emission maps and spectrum for a Coma-like prototype

In this section we explore the spectral and spatial properties of $`\gamma `$-ray radiation between 10 keV and 10 TeV due to shock accelerated CRs in GCs. As non-thermal bremsstrahlung turns out unimportant , we consider the following emission processes: $`\pi ^0`$-decay and IC emission from both primary CR electrons and secondary e± . We first compute the emission spectrum for the largest virialized object in the computational box which has an X-ray core temperature $`T_x4`$ keV. The various emission components are then renormalized to the case of a Coma-like cluster. In particular, IC emission from electrons is rescaled according to equation LABEL:fm:gamfit.eq. The total number of e± is rescaled assuming that these particles are responsible the for synchrotron emission of Coma radio halo ***This assumption is not necessary but it is assumed for simplicity. For the purpose we took a radio flux at 1.4 GHz $`S_{1.4GHz}=640`$ mJy  and a two values for the volume average magnetic field $`B`$, namely 1 and 3 $`\mu `$G .

In the left quadrant of fig. LABEL:fm:figc4.fig we present synthetic maps of the integrated photon flux above 100 keV (top panels) and 100 MeV (bottom panels). They show that the emission from $`\pi ^0`$-decay (bottom left) and e± (top left) is confined to the cluster core where it creates a diffuse halo which rapidly fades with distance from the center. In fact, e± and $`\pi ^0`$ are produced at the highest rate in the densest regions where both the parent CR ions and target nuclei are most numerous. On the other hand, because of severe energy losses, $`\gamma `$-ray emitting primary electrons are only found in the vicinity of strong shocks where they are accelerated. The emission (right panels) is extended and with a rich morphology reflecting the complex “web” of accretion shocks surrounding GCs . The corresponding synthetic spectra extracted from a core (top; with a $`0.5^o`$ radius) and an outskirts region (bottom; a ring with inner and outer radii of $`0.5^o`$ and $`1.5^o`$ respectively) are shown in the central quadrant of fig. LABEL:fm:figc4.fig. They shows that the emitted radiation in the outskirts region is strongly dominated by IC emission from primary primary e-. Conversely in the core region $`\pi ^0`$-decay (solid thin line) dominates at high photon energy ($`>`$ 100 MeV) (top panel). Below $``$ 100 MeV the relative contribution from primary and secondary e± depends on the strength of the magnetic field. Finally, in the right quadrant the predicted radiation spectra in various spatial regions are compared to sensitivity limitsThese are for point sources and will be worse for extended sources considered here. of planned $`\gamma `$-ray observatories.

### 5. Conclusions

Clusters observations at $`\gamma `$-ray energy will provide important information about GCs. First, the extended emission is produced at the location of accretion shocks which, unlike weak merging shocks, have yet to be observed directly. Thus, detection and imaging of IC emission from primary electrons would provide an opportunity for their direct observation. The flux about 100 keV together with radio measurements, will allow an estimate of the intracluster magnetic field. It is important, however, to separate out the contribution from the CR electrons accelerated at the outlining shocks. In fact, because the magnetic field strength is expected to drop at the cluster outskirts, these electrons likely generate only a weak radio emission. Nonetheless, they produce a strong IC flux which can easily dominate the total soft $`\gamma `$-ray emission (cf. fig. LABEL:fm:figc4.fig). Measuring the $`\gamma `$-ray flux at and above 100 MeV will allow us to confirm or rule out secondary e± models for radio halos in GC  and to estimate the non-thermal CR pressure there . In this perspective several authors estimated for nearby clusters the $`\gamma `$-ray flux from $`\pi ^0`$-decay. However, radiation flux from IC emission can be comparable . Therefore, for a correct interpretation of the measurements the contributions from these two components will need to be separated as outlined here.

Measuring the spectrum of IC $`\gamma `$-ray from primary CR electrons will also be instrumental for addressing their contribution to the CGB , the details of their acceleration mechanism and the physical conditions of the shocks.

### 6. Acknowledgments

This work is partially supported by the European Community under the contract HPRN-CT2000-00126 RG29185 (The Physics of the Intergalactic Medium). The computational work was carried out at the the Rechenzentrum in Garching.

1. Allen, G. E., Petre, R., & Gotthelf, E. V. 2001, ApJ, 558, 739

2. Blandford, R. D. & Eichler, D. 1987, Phys. Rep. , 154, 1

3. Blasi, P. & Colafrancesco, S. 1999, Astropart. Phys. , 12, 169

4. Bowyer, S., Berghoefer, T. W., & Korpela, E. 1999, ApJ, 526, 592

5. Deiss, B. M., Reich, W., Lesch, H., & Wielebinski, R. 1997, A&A, 321, 55

6. Dolag, K. & Enßlin, T. 2000, A&A, 362, 151

7. Fusco-Femiano, R., et al. 1999, ApJ, 513, L21

8. Keshet, U., et al. 2002, ApJ, in press (astro-ph/0202318)

9. Kim, K.-T., et al. 1990, ApJ, 355, 29

10. Lieu, R., et al. 1996, Science , 274, 1335

11. Loeb, A. & Waxman, E. 2000, Nature (London) , 405, 156

12. Müller, D. & et al. 1995, in Int. Cosmic Ray Conference, Vol. 3, Rome, 13

13. Miniati, F. 2001, Comp. Phys. Comm. , 141, 17

14. Miniati, F. 2002, MNRAS, 337, 199

15. Miniati, F. 2002, MNRAS, submitted

16. Miniati, F., Jones, T. W., Kang, H., & Ryu, D. 2001, ApJ, 562, 233

17. Miniati, F., et al. 2000, ApJ, 542, 608

18. Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559,

19. Reimer, O. 1999, in ICRC, Vol. 4, ed. D. Kieda, et al., Salt Lake City, 89

20. Sreekumar, et al. 1996, ApJ, 464, 628

21. Sreekumar, et al. 1998, ApJ, 494, 523

22. Völk, H. J., Aharonian, F. A., & Breitschwerdt, D. 1996, Sp.Sci.Rev., 75, 279