text:
# Current-density functional for disordered systems

## I Introduction

The currently used method to obtain observables for disordered systems is to perform a quenched averaging, an averaging of the Green functions over random external field . The computational algorithm boils down either to the analytical continuation in the number of replicas , or to the introduction of fictious particles related to the real ones by super-transformations  or to the use of the Keldysh contour in computing loop-integrals . The impact of weak disorder on the conductivity has nicely been captured by resumming the dominant graphs  in a self-consistent manner . But the whole scheme is based on perturbation expansion and it remains unclear how to treat strongly disordered systems and understand the transition to localization. This feature appears to be a serious limiting factor in the description of strongly correlated systems, as well, where the interactions are supposed to be non-perturbative. The goal of this paper is to demonstrate that methods, well established in other domains of physics, in particular the density functional scheme and the functional extension of the renormalization group might be useful to tackle this problem. Namely, a new method is proposed for the systemical, non-perturbative computation of the queneched average of the one-particle irreducible (1PI) vertex functions.

The density functional  turned out to be a powerful formalism in handle the ground states. Its virtue is that the variational setting gives a large degrees of flexibility for improvements of the density functional. The extension to include the current in the variational techniques has already been suggested, as well . This method is improved on the following points. First, the density-current functional is generalized for time-dependent density and current configurations in order to encompass the dynamics and to arrive at a closed set of equations to determine this functional. Second, a constructive definition of the functional is used which is based on the Legendre transform of the generator functional for the connected Green functions for the density and current. In fact, it has been realized lately that such an effective action can be identified with the density functional . Third, the functional generalization  of the traditional renormalization group procedure  is used in the internal space  to obtain the effective action. The evolution equation describes the changes generated by the gradual increase of the amplitude of the fluctuations. The only approximation committed in this scheme is the truncation of the effective action, the projection of the functional renormalization group equation into a constrained functional space. One hopes to carry out this truncation in a manner which is independent of the amplitude of the fluctuations and the resulting scheme will be non-perturbative. We shall apply gradient expansion and retain arbitrarily high orders which seems to be a consistent scheme in determing the kinetic transport coefficients.

The functional renormalization group method has already been applied for disordered systems, namely to the pinning of an elastic system in a random potential . A random external potential introduced in an elastic system generates infinitely many relevant or marginal coupling constants and their treatment required the functional extension of the renormalization group procedure. The U.V. cut-off was gradually decreased and the resulting renormalization group equation was solved within the framework of the replica method. In the present work we do not touch the cut-off, instead we follow the evolution of the system when a coupling constant is gradually ’turned on’. Such a modification of the rules of the renormalization group method allows us to embark directly the issue of quenched averaging and offers an alternative to the replica method.

The initial condition is usually chosen to be a theory with such parameters which allow an approximative determination of the effective action. The evolution equation generates the ’renormalized trajectory’ along which the physics changes since the quantum/thermal fluctuations are turned on gradually. The integration of the evolution equation from the initial condition with small fluctuations to the strongly coupled theory under consideration represents a general purpose algorithm to solve strongly interactive systems.

It is shown below that the quenched average of the generator functional for the connected Green functions for the density satisfies a simple evolution equation when the strength of the disorder is gradually changed. This is not surprising for weak disorder since the perturbative averaging over the impurity potential introduces correlations which are represented by introducing ’interaction vertices’ for the density and the corresponding perturbation series involves Green functions with increasing number of density insertions. The evolution equation is a differential equation for the generator functional for the Green functions which expresses the change induced by the infinitesimal increase of the strength of the disorder in terms of functional derivatives with respect to the external source coupled to the density. In case of an interactive system the coupling constants of the annealed interactions are evolved, too. The equation is obtained without any reference to small parameters in the dynamics and is valid non-perturbatively. The evolution equation, being a functional differential equation represents a wonderful mathematical problem. In lacking theorems and other general support one has to truncate this equation by projecting it into a restricted functional space. The ordinary functions characterizing the functionals in this space, the running coupling strengths in the language of the renormalization group, will then be evolved by differential equations. It is rather natural that local functionals are easier to handle in a truncation scheme. Therefore the evolution equation is recasted in terms of a local functional, the effective action for the density.

The generator functionals are well suited to the linear response formalism since the introduction of the external sources can easily be achieved by functional derivation with respect to appropriately chosen external sources. Therefore we include the electric current besides the density in our formalism in order to keep track of issues related to charge transport and will use the generator functional for the connected Green functions for the current and density operators and their effective action. The formal evolution equation will be presented for the case of electrons with Coulomb interaction but the numerical solution will be discussed for the non-interactive case only.

We start in section II with the introduction of the model and its generator functionals. The average of the generator functional for the connected Green functions for the current and the density operators over the impurity field is defined in the next step. Its functional Legendre transform, the effective action or the current-density functional is the main object of this paper. Connection to the density functional theory is briefly commented, as well. It is furthermore mentioned that the generator functional is well suited to the linear response computations except that the causal propagators have to be replaced by the retarded or advanced ones in order to implement the physical boundary conditions in time. This gives a simple short-cut to re-derive the Kubo expressions for the transport coefficients.

The evolution equation method is discussed in section III. Though the complete interactive case is covered, the equation is studied in a more detailed manner for non-interactive particles only. The functional differential equation is approximated by truncating either the powers of the fields or the power of the space-time gradients acting on the fields. The former is useful to establish a contact with the resummation methods of the perturbation series, the latter is more appropriate for the computation of the transport coefficients.

Before solving the evolution equation we have to construct the initial condition. This is the subject of section IV. The initial condition is defined in the absence of disorder where the effective action is obtained by retaining infinitely many orders of the gradient expansion. First we compute the one-loop two-point functions for the current and density and construct the generator functional. After that we perform the functional Legendre transformation to find the effective action.

The brief survey of the numerical solution of the evolution equation is presented in section V. The extrapolation of the conductivity to the homogeneous case found out to be rather difficult for weak disorder and should be studied in a more detailed manner. But strong disorder suppresses the long range modes and the finite size dependence and allows to study the conductivity on lattices up to the size $`320^3`$. The preliminary results presented in this paper show a discontinuity developing in the static limit, in particular the d.c. conductivity drops to zero at a certain strength of the disorder. It is natural to identify this point with the onset of localization.

Section VI is for the conclusion. The technical details of the functional calculus, a formal relation between the evolution equation and the replica method and the perturbative evaluation of the two-point functions for the current and density are presented in the appendices together with the derivation of the Ward identities which are needed to assure the finiteness of the d.c. conductivity.

## II Current-density functional

The model will be introduced in this section together with the generator functional for its connected Green functions, its quenched average and functional Legendre transform, the effective action.

### A The model

We consider spinless electrons propagating on an external electromagnetic potentials $`V_x^{\mathrm{e}xt}`$ and $`A_x^{\mathrm{e}xt}`$ and eventually interacting with the Coulomb potential. The path integral representation of the vacuum-to-vacuum amplitude is

$$Z=𝒟[\psi ]𝒟[\psi ^{}]𝒟[u]e^{\frac{i}{\mathrm{}}S[\psi ^{},\psi ,u]},$$
(1)
with

$$S[\psi ^{},\psi ,u]=_x\left\{\psi _x^{}\left[i\mathrm{}_t\frac{1}{2m}\left(\frac{\mathrm{}}{i}+A_x^{\mathrm{e}xt}\right)^2+\mu +eV_x^{\mathrm{e}xt}+eu_x\right]\psi _x\frac{1}{2}(u_x)^2\right\},$$
(2)
c.f. Appendix A for the notations. We shall use a space-time lattice with spacings $`a_j=a`$, $`j=1,2,3`$ and $`a_0=a_\tau `$. There several reasons to introduce lattice regularization. The most important one is that the non-relativistic models are usually non-renormalizable and need cut-off for their definition. But even if we are satisfied by considering non-interacting electrons, as will be the case in most of this paper, which are U.V. finite the functions appearing in the description such as the Green functions develop I.R. singularities. This requires a special care in the numerical treatment of the one-loop integral of the evolution equation. The breakdown of this loop integral into a large but finite sum can be done in a singularity free manner only if the gauge invariance properties are not spoiled in the finite sum. The only known method to achieve this is to use lattice regularization. The regulated action will be

$$S[\psi ^{},\psi ,u]=a^6a_\tau ^2\underset{x^{},x}{}\psi _x^{}^{}G_{x^{},x}^1\psi _x+\frac{a^3a_\tau }{2}\underset{x}{}(u_x)^2,$$
(3)
in terms of the inverse propagator

$`G_{x^{},x}^1[V^{\mathrm{e}xt},A^{\mathrm{e}xt}]`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{a_\tau }}\left(\delta _{x^{},x}^K\delta _{x^{},x+\widehat{0}}^KU_0(t,𝐱)\right)`$ (5)
$`+{\displaystyle \frac{\mathrm{}^2}{2ma^2}}{\displaystyle \underset{j>0}{}}\left(\delta _{x^{}+\widehat{j},x}^KU_j^{}(x\widehat{j})+\delta _{x^{}\widehat{j},x}^KU_j(x)2\delta _{x^{},x}^K\right),`$
where $`\widehat{\mu }=a_\mu e_\mu `$ and the link variables

$$U_0(x)=e^{\frac{i}{\mathrm{}}a_\tau (\mu +eV_x^{\mathrm{e}xt}+eu_x)},U_j(x)=e^{\frac{i}{\mathrm{}}eaA_{j,x}^{\mathrm{e}xt}}$$
(6)
are used to assure gauge invariance .

The generator functional for the connected Green functions of local operators $`𝒪_{\alpha ,x}`$,$`\stackrel{~}{𝒪}_{\stackrel{~}{\alpha },x}`$ coupled to the external sources $`\sigma _{\alpha ,x}`$, $`\stackrel{~}{\sigma }_{\stackrel{~}{\alpha },x}`$, $`\alpha =0,1,2,3`$, $`\stackrel{~}{\alpha }=1,2,3`$, is

$$e^{\frac{i}{\mathrm{}}W_0[\sigma ]}=𝒟[u]𝒟[\psi ]𝒟[\psi ^{}]e^{\frac{i}{\mathrm{}}\psi ^{}(G^1[V^{\mathrm{e}xt},A^{\mathrm{e}xt}]+\sigma /)\psi +\frac{i}{2\mathrm{}}uu},$$
(7)
where $`\psi ^{}\sigma /\psi =\sigma _0𝒪_0+\sigma _j𝒪_j+\stackrel{~}{\sigma }_j\stackrel{~}{𝒪}_j`$. We shall choose

$`\psi ^{}𝒪_{0,x}\psi `$ $`=`$ $`{\displaystyle \frac{\delta S}{\delta eV_x^{\mathrm{e}xt}}}=i\psi _{x+\widehat{0}}^{}U_0(x)\psi _x,`$ (8)
$`\psi ^{}𝒪_{j,x}\psi `$ $`=`$ $`{\displaystyle \frac{\delta S}{\delta eA_{j,x}^{\mathrm{e}xt}}}={\displaystyle \frac{i\mathrm{}}{2ma}}{\displaystyle \underset{j>0}{}}\left[\psi _{x+\widehat{j}}^{}U_j(x)\psi _x\psi _x^{}U_j^{}(x)\psi _{x+\widehat{j}}\right],`$ (9)
$`\delta _{j,k}\delta _{x,y}\psi ^{}\stackrel{~}{𝒪}_{j,x}\psi `$ $`=`$ $`m{\displaystyle \frac{\delta \psi ^{}𝒪_{j,x}\psi }{\delta eA_{k,y}^{\mathrm{e}xt}}}={\displaystyle \frac{1}{2}}\delta _{j,k}\delta _{x,y}\left[\psi _{x+\widehat{j}}^{}U_j(x)\psi _x+\psi _x^{}U_j^{}(x)\psi _{x+\widehat{j}}\right].`$ (10)
The operators $`\stackrel{~}{𝒪}_{j,x}`$ are introduced in order to derive simplify the formal steps in deriving the Ward identities.

### B Averaging over the disorder

Let us suppose that the presence of the impurities can be described by a static potential $`v_𝐱`$, with Gaussian distribution . The quenched average of observables corresponding to a typical static impurity configuration can be constructed by means of the generator functional for averaged connected Green functions,

$$W[\sigma ]=\frac{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}W_0[\sigma _0+v,\sigma _1,\mathrm{}]}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}.$$
(11)
Two remarks are in order at this point. First, local impurity potential distribution is assumed for the sake of simplicity only, each step below can be reproduced with any other Gaussian distribution. Second, such a simple averaging is available for the Green functions of the density only, the operator which couples linearly to the impurity potential.

Since the impurity potential appears as an additional external source coupled to the density one expects that its modification of the dynamics can be obtained by means of the external potential $`\sigma _0`$. A simple explicit equation which realizes this possibility can immediately be obtained by rewriting (11) as

$$W[\sigma ]=W_0[\sigma _0+\frac{\delta }{\delta j},\sigma _1,\mathrm{}]e_{|j=0}^{\frac{g}{2}_𝐱j_𝐱^2}.$$
(12)
This equation can be used to generate the usual perturbation series .

Another variant of Eq. (11) which will be useful for the derivation of the evolution equation is established by using the ’equation of motion’ for the impurity potential. The equation of motion expresses the invariance of the path integral with respect to the infinitesimal shift of the integral variable,

$$0=𝒟[\varphi ]e^{\frac{i}{\mathrm{}}S[\varphi ]}\frac{\delta S[\varphi ]}{\delta \varphi _x}.$$
(13)
The finite shift $`v_𝐩v_𝐩i\sigma _{\omega =0,𝐩}/T`$ where where $`T`$ is the Euclidean time extent of the system carried out for the Fourier transform gives

$$W[\sigma ]=\frac{𝒟[v]e^{\frac{1}{2g}_𝐩(v_p\frac{1}{T}\sigma _{0,\omega =0,p})(v_𝐩\frac{1}{T}\sigma _{0,\omega =0,𝐩})}W_0^v[\sigma ]}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}.$$
(14)
where shifted functional is

$$W_0^v[\sigma _0,\sigma _1,\mathrm{}]=W_0[\sigma _{0,\omega =0}=Tv,\sigma _{0,\omega 0},\sigma _1,\mathrm{}].$$
(15)
This equation explains the use of the internal space blocking method for the computation of quenched averages. In this blocking procedure a parameter which controls the amplitude of the fluctuations, say the mass for relativistic scalar theories, is driving the evolution. In the present case $`g`$, the average amplitude square of the impurity field, plays a similar role. In fact, the impurity field fluctuations just smear out the generator functional within an interval of the size $`𝒪(\sqrt{g})`$ in the dependence on the external source. The evolution equation introduced below follows as the increase of $`g`$ gradually ’opens up’ a large interval for the source in which the disorder acts.

The functional $`W[\sigma ]`$ will be useful to express the change of the dynamics under an infinitesimal change of a parameter of the dynamics but leads to complications when such an evolution equation is truncated. The complications arise from the fact that $`W[\sigma ]`$ belongs to the class of non-local functionals which are very hard to characterize and handle. The way out form this problem is the introduction of the current-density functional, the effective action for the composite operators $`𝒪_\alpha `$ as the Legendre transform of $`W[\sigma ]`$,

$$\mathrm{\Gamma }[\rho ]=W[\sigma ]+\sigma _\alpha \rho _\alpha +\stackrel{~}{\sigma }_{\stackrel{~}{\alpha }}\rho _{\stackrel{~}{\alpha }},\rho _{\alpha ,x}=\frac{\delta W[\sigma ]}{\delta \sigma _{\alpha ,x}},\stackrel{~}{\rho }_{\stackrel{~}{\alpha },x}=\frac{\delta W[\sigma ]}{\delta \stackrel{~}{\sigma }_{\stackrel{~}{\alpha },x}}.$$
(16)
This functional is local and offers useful truncation schemes.

A relation which will later be used a number of occasions is

$$_x^{}\frac{\delta ^2W[\sigma ]}{\delta \sigma _{\alpha ,x}\delta \sigma _{\gamma ,x^{}}}\frac{\delta ^2\mathrm{\Gamma }[\rho ]}{\delta \rho _{\gamma ,x^{}}\delta \rho _{\beta ,x^{\prime \prime }}}=\delta _{x,x^{\prime \prime }}\delta _{\alpha ,\beta }.$$
(17)
### C Variational method

The Hohenberg-Kohn theorems  assure that (i) the density functional $`\mathrm{\Gamma }[\rho ]`$ which gives the ground state energy in the sector of the Fock space belonging to a given charge density $`\rho _x`$ is a well defined functional and (ii) the minimum of the density functional is reached at the density profile of the true ground state where its value reproduces the energy of the ground state. We generalize the density functional first to include the current $`\mathrm{\Gamma }[\rho ]\mathrm{\Gamma }[\rho _\mu ]`$, where $`\rho _0`$ and $`\rho _j`$ ($`j=1,2,3`$) is the expectation value of the density and the current, respectively. After that we allow the density and the current to become space and time dependent.

This generalization is realized by the effective action (16). In fact, the inverse transformation is carried out by the change of variable

$$\sigma _{\alpha ,x}=\frac{\delta \mathrm{\Gamma }[\rho ]}{\delta \rho _{\alpha ,x}},\stackrel{~}{\sigma }_{\stackrel{~}{\alpha },x}=\frac{\delta \mathrm{\Gamma }[\rho ]}{\delta \stackrel{~}{\rho }_{\stackrel{~}{\alpha },x}},$$
(18)
and this relation indicates that the ground state corresponds to the stationary points of the effective action. The transformation (16) can be used to argue that the effective action gives the quantum action ($`i\mathrm{}`$ times the logarithm of the transition amplitude) for real time or the free energy for imaginary time. Therefore the Euclidean effective action is bounded from below and the ground state is found at its minimum. Therefore the effective action plays the role of the current-density functional whose minimization leads us to the ground state.

### D Linear response

As mentioned in the Introduction, the generator functional is well suited to the linear response formalism. To demonstrate this point let us consider the system in the presence of an additional external electromagnetic field $`A_\mu `$ imposed for $`t<0`$ by making the replacement $`V_t^{\mathrm{e}xt}V_t^{}=V_t^{\mathrm{e}xt}+A_{0,t}`$ and $`A_t^{\mathrm{e}xt}A_t^{}=A_t^{\mathrm{e}xt}+A_t`$ in (5) where the eventual time dependence arising from the external fields is shown explicitly. For the computation of the conductivity $`j_k=\sigma _{k,\mathrm{}}E_{\mathrm{}}`$, we use the electric field $`𝐄=_t𝐀`$ and the interaction Hamiltonian $`H^{}=e_xA_{k,x}j_{k,x}`$. The ground state $`|0_t`$, following the external fields in an adiabatic manner might show time dependence, as well. The expectation value $`0_t|j_{k,t,𝐱}|0_t`$ of the induced gauge invariant current for $`t=0`$ is

$`0_0|j_{k,0,𝐱}|0_0`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}e}{2ma}}0_0|\psi _{s,0,𝐱+k}^{}U_k^{}(x)\psi _{s,0,𝐱}\psi _{s,0,𝐱}^{}U_k^{{}_{}{}^{}}(x)\psi _{s,0,𝐱+k}|0_0`$ (19)
$``$ $`e0_0|𝒪_{k,0,𝐱}|0_0{\displaystyle \frac{e^2}{2m}}A_{k,0,𝐱}0_0|\psi _{s,0,𝐱+k}^{}U_k(x)\psi _{s,0,𝐱}+\psi _{s,0,𝐱}^{}U_k^{}(x)\psi _{s,0,𝐱+k}|0_0`$ (20)
$`=`$ $`e{\displaystyle _y}\mathrm{\Theta }(y^0){\displaystyle \frac{\delta 0_0|𝒪_{k,0,𝐱}|0_0}{\delta A_{\mathrm{},y}}}A_{\mathrm{},y}{\displaystyle \frac{e^2}{m}}A_{k,0,𝐱}\stackrel{~}{\rho }^{}.`$ (21)
Since the bra $`0_0|`$ is the time reversed version of the ket $`|0_0`$ in the absence of spin the functional derivative in the last equation contains the contributions of the bra and the ket with a relative negative sign which corresponds to the commutator between the external source and the observable in the Kubo formula. The time inversion can formally be realized in the framework of the path integral formalism by using external field which is symmetrized with respect to the time inversion and by taking the complex conjugate of the integrand for $`t>0`$. This latter step brings the product of the Heavyside function and the functional derivative into the retarded current-current Green function,

$$0_0|j_{k,0,𝐱}|0_0=e^2_y\stackrel{~}{D}_{(0,𝐱),y}^{Rk,\mathrm{}}A_{\mathrm{},y}\frac{e^2}{m}A_{k,0,𝐱}\stackrel{~}{\rho }^{},$$
(22)
and we have

$$\sigma _{k,\mathrm{}}=\frac{ie^2}{\omega }\left(\frac{\stackrel{~}{\rho }^{}}{m}\delta ^{k,\mathrm{}}+\stackrel{~}{D}_{\omega ,𝐪}^{Rk,\mathrm{}}\right)=\frac{ie^2}{\omega }\left(\frac{\stackrel{~}{\rho }^{}}{m}\delta ^{k,\mathrm{}}+Im\stackrel{~}{D}_{\omega +i0^+,𝐪}^{k,\mathrm{}}\right)$$
(23)
in agreement with the Kubo formula. In the last equation we have the imaginary part of the causal propagator because the current operator is hermitian.

## III Evolution equation

We introduce the evolution equation and discuss its structure for non-interactive systems first. The generalization which includes interaction, the Coulomb potential in particular is presented briefly later.

In the traditional renormalization group method one orders the degrees of freedom $`\varphi _1,\varphi _2,\mathrm{},\varphi _N`$ according to their momenta, $`p_1<p_2<\mathrm{}<p_N`$ and eliminates them one-by one, starting with the perturbative, high momentum modes. The corresponding blocking for the action is

$$e^{iS_{N1}(\varphi _1,\varphi _2,\mathrm{},\varphi _{N1})}=𝑑\varphi _Ne^{iS_N(\varphi _1,\varphi _2,\mathrm{},\varphi _N)}.$$
(24)
Such a method may be called external space renormalization group since the momentum is a scale in the external space, i.e. in the space-time. As an alternative, one may order the fluctuations according to their amplitude, the scale parameter in the internal space of quantum field theory . The evolution from $`\lambda =0`$ to $`\lambda =\mathrm{}`$ of the partition function

$$Z_\lambda =_\lambda ^\lambda 𝑑\varphi _1\mathrm{}_\lambda ^\lambda 𝑑\varphi _Ne^{iS_N(\varphi _1,\varphi _2,\mathrm{},\varphi _N)}$$
(25)
takes into account the fluctuations in the order of increasing amplitude. A more practical, smooth cut-off in the internal space can be implemented by introducing a suppression term in the action $`S[\varphi ]S[\varphi ]+S_\lambda [\varphi ]`$, such as

$$S_\lambda [\varphi ]=\frac{\lambda }{2}\varphi \varphi .$$
(26)
In relativistic field theory this term generates mass and the realizes the Callan-Symanzik scheme .

The evolution equation can easiest be derived by starting with the rather trivial identity

$$\frac{d}{d\lambda }Z[j]=\frac{d}{d\lambda }𝒟[\varphi ]e^{\frac{i}{\mathrm{}}S[\varphi ]+\frac{i}{\mathrm{}}\varphi j}=\frac{i}{\mathrm{}}𝒟[\varphi ]e^{\frac{i}{\mathrm{}}S[\varphi ]+\frac{i}{\mathrm{}}\varphi j}\frac{dS[\varphi ]}{d\lambda }.$$
(27)
In the second step the expectation value standing in the right hand side is expressed by acting with the functional derivatives with respect to the source $`j`$ on the generator functional $`Z[j]`$ and the result is a functional differential equation for $`Z[j]`$. The harmless looking first step is actually a generalization of the Hellman-Feynman theorem  for the time evolution operator. Different choices for $`\lambda `$ yield important, well known equations. The choice of $`\lambda `$ as a field variable $`\varphi (x)`$ leads to the Schwinger-Dyson equation. When $`\lambda `$ is chosen to be a parameter of a transformation of the field variable which leaves a part or the full action invariant then the resulting equations are the Ward identities.

The disadvantage of the internal space renormalization group method is that the renormalized trajectory does not belong a to fixed physics anymore and the evolution of the control parameter $`\lambda `$ has little to do with the scale dependence in the theory. But the advantage is the flexibility in choosing $`\lambda `$. In fact, by identifying $`\lambda `$ with the Planck constant or a gauge coupling one can construct evolution schemes where the semiclassical structure remains invariant or the gauge invariance is preserved explicitely, respectively. The solution of the evolution from a perturbative initial condition with small amplitude fluctuations to the desired level of fluctuations provides us a systematical non-perturbative algorithm to compute irreducible Green functions.

### A Non-interacting particles

It was mentioned in Section II B that the effects of the disorder and the external source $`\sigma _0`$ are related. We follow up this remark by showing that the infinitesimal change of the strength of the disorder can be reproduced by acting on $`W[\sigma ]`$ with functional derivatives with respect to $`\sigma _0`$. In order to keep the measure of the strength of the disorder, $`g`$, explicit in the equations we make the replacement $`gg\lambda `$ and look for the evolution equation in $`\lambda `$,

$`\dot{\mathrm{\Gamma }}[\rho ]`$ $`=`$ $`\dot{W}[\sigma ]`$ (28)
$`=`$ $`{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g\lambda }_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}W_0^v[\sigma ]_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}{2g\lambda ^2𝒟[v]e^{\frac{1}{2g\lambda }_𝐱v_𝐱^2}}}`$ (30)
$`+{\displaystyle \frac{1}{2g\lambda ^2}}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g\lambda }_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}W_0^v[\sigma ]}{𝒟[v]e^{\frac{1}{2g\lambda }_𝐱v_𝐱^2}}}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g\lambda }_𝐱v_𝐱^2}_𝐱v_𝐱^2}{𝒟[v]e^{\frac{1}{2g\lambda }_𝐱v_𝐱^2}}}`$
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐩}{\displaystyle \frac{\delta ^2W[\sigma ]}{\delta \sigma _{0,\omega =0,p}\delta \sigma _{0,\omega =0,𝐩}}}`$ (31)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐩}\left[{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ]}{\delta \rho \delta \rho }}\right]_{(0,\omega =0,p),(0,\omega =0,𝐩)}^1,`$ (32)
where the dot stands for $`_\lambda `$. We shall impose the initial conditions at $`\lambda =0`$. It is shown in Appendix B that this formalism is consistent with the replica method so long the latter is applicable.

The functional differential equations are too difficult to handle and the step of central importance is to find a truncation scheme which renders this equation useful. We consider two truncations schemes following the Landau-Ginzburg strategy of a double expansion in the amplitude and the gradient of the fields. The truncation in the amplitude has more pedagogical value, it provides a relation between the resummation of the perturbation series and the evolution equation. The truncation in the powers of the momenta, the gradient expansion, seems to be a natural approximation scheme for the effective action since the averaging over the impurity field leaves the ground state translation invariant and the conductivity, as other transport parameters, are closely related to the Taylor expansion of the 1PI Green functions in the external momentum.

### B Coulomb interaction

We considered the system of non-interacting electrons so far. But the evolution equation method can be used for interactive systems, too. The non-trivial point is that both the annealed and the quenched averages can be computed in a parallel manner by generalizing the method of Ref.  for disordered systems.

The generator functional for the connected Green functions for real time involving the rescaled photon field $`uu/e`$ is

$$e^{\frac{i}{\mathrm{}}W_0[\sigma ,j]}=𝒟[u]𝒟[\psi ]𝒟[\psi ^{}]e^{\frac{i}{\mathrm{}}\psi _s^{}(G^1+\sigma /)\psi _s+\frac{i}{2\mathrm{}e^2}uu+\frac{i}{\mathrm{}}(j+\psi _s^{}\psi _s)u}$$
(33)
We make the rescaling $`e^2e^2\lambda `$ and find the evolution equation for the annealed average

$`\dot{W}_0[\sigma ,j]`$ $`=`$ $`{\displaystyle \frac{1}{2e^2\lambda ^2}}{\displaystyle _x}\mathrm{\Delta }_xu_xu_y_{|x=y}`$ (34)
$`=`$ $`{\displaystyle \frac{1}{2e^2\lambda ^2}}\left(i\mathrm{}\mathrm{T}r\mathrm{\Delta }{\displaystyle \frac{\delta ^2W_0[\sigma ,j]}{\delta j\delta j}}{\displaystyle \frac{\delta W_0[\sigma ,j]}{\delta j}}\mathrm{\Delta }{\displaystyle \frac{\delta W_0[\sigma ,j]}{\delta j}}\right),`$ (35)
in the electric charge. The functional we need is effective action corresponding to the quenched average of the connected Green functions,

$$\mathrm{\Gamma }[\rho ,w]=W[\sigma ,j]+\sigma _\alpha \rho _\alpha +jw,\rho _{\alpha ,x}=\frac{\delta W[\sigma ,j]}{\delta \sigma _{\alpha ,x}},w_x=\frac{\delta W[\sigma ,j]}{\delta j_x},$$
(36)
where

$$W[\sigma ,j]=\frac{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}W_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}.$$
(37)
In order to perform the quenched averaging by evolving the system we make the rescaling $`gg\lambda `$ in the distribution of the impurity potential and find

$`\dot{\mathrm{\Gamma }}[\rho ,w]`$ $`=`$ $`\dot{W}[\sigma ,j]`$ (38)
$`=`$ $`{\displaystyle \frac{1}{2g\lambda ^2}}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g}_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}W_0^v[\sigma ,j]_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}}`$ (41)
$`+{\displaystyle \frac{V}{2\lambda }}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g}_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}W_0^v[\sigma ,j]}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}}`$
$`+{\displaystyle \frac{1}{2e^2\lambda ^2}}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}\left(i\mathrm{}\mathrm{T}r\mathrm{\Delta }\frac{\delta ^2W_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{\delta j\delta j}\frac{\delta W_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{\delta j}\mathrm{\Delta }\frac{\delta W_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{\delta j}\right)}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}}.`$
The last line, the piece which generates the Coulomb interaction is represented graphically in Fig. 1.

The formal difference between the annealed and the quenched averages appears in the manner the disconnected parts are treated in the last two terms on the right hand side of the second equation. In order to linearize the right hand side in the functionals which facilitates the computation of the quenched average we introduce a new functional

$$K_0[\sigma ,j]=\frac{\delta W_0[\sigma ,j]}{\delta j}\mathrm{\Delta }\frac{\delta W_0[\sigma ,j]}{\delta j}$$
(42)
whose quenched average

$$K[\sigma ,j]=\frac{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}K_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}}$$
(43)
satisfies the evolution equation

$`\dot{K}[\sigma ,j]`$ $`=`$ $`{\displaystyle \frac{1}{2g\lambda ^2}}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g}_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}K_0^v[\sigma ,j]_𝐩|v_p\frac{1}{T}\sigma _{0,\omega =0,p}|^2}{𝒟[v]e^{\frac{1}{2g\lambda }_xv_x^2}}}`$ (45)
$`+{\displaystyle \frac{V}{2\lambda }}{\displaystyle \frac{𝒟[v]e^{\frac{1}{2g\lambda }_xv_x^2}K_0[\sigma _0+v,\sigma _1,\mathrm{},j]}{𝒟[v]e^{\frac{1}{2g\lambda }_xv_x^2}}}`$
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐳}{\displaystyle \frac{\delta ^2K[\sigma ,j]}{\delta \sigma _{0,0,𝐳}\delta \sigma _{0,0,𝐳}}}.`$ (46)
We write $`L[\rho ,w]=K[\sigma ,j]`$ and the corresponding evolution equation is

$`\dot{L}[\rho ,w]`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x}}{\displaystyle \frac{\delta ^2\rho _x}{\delta \sigma _{0,0,𝐳}\delta \sigma _{0,0,𝐳}}}{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x,y}}{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho _x\delta \rho _y}}{\displaystyle \frac{\delta \rho _x}{\delta \sigma _{0,0,𝐳}}}{\displaystyle \frac{\delta \rho _y}{\delta \sigma _{0,0,𝐳}}}`$ (47)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x}}{\displaystyle \frac{\delta ^3W[\sigma ,j]}{\delta \sigma _x\delta \sigma _{0,0,𝐳}\delta \sigma _{0,0,𝐳}}}{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x,y}}{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho _x\delta \rho _y}}{\displaystyle \frac{\delta ^2W[\sigma ,j]}{\delta \sigma _x\delta \sigma _{0,0,𝐳}}}{\displaystyle \frac{\delta ^2W[\sigma ,j]}{\delta \sigma _y\delta \sigma _{0,0,𝐳}}}`$ (48)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x}}{\displaystyle \frac{\delta }{\delta \sigma _x}}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{(0,0,𝐳),(0,0,𝐳)}^1{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}`$ (50)
$`{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x,y}}{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho _x\delta \rho _y}}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{x,(0,0,𝐳)}^1\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{y,(0,0,𝐳)}^1`$
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _{𝐳,x,y,u}}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{(0,0,𝐳),y}^1{\displaystyle \frac{\delta }{\delta \sigma _x}}{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho _y\delta \rho _u}}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{u,(0,0,𝐳)}^1{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}`$ (52)
$`{\displaystyle \frac{g}{2}}{\displaystyle _𝐳}\left[\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)^1{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho \delta \rho }}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)^1\right]_{(0,0,𝐳),(0,0,𝐳)}`$
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐳}\left[\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)^1\left({\displaystyle _x}{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}{\displaystyle \frac{\delta }{\delta \sigma _x}}{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}+{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho \delta \rho }}\right)\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)^1\right]_{(0,0,𝐳),(0,0,𝐳)}`$ (53)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐳}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)^1\left({\displaystyle _{x,y}}{\displaystyle \frac{\delta L[\rho ,w]}{\delta \rho _x}}\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{x,y}^1{\displaystyle \frac{\delta ^3\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho \delta \rho _y}}+{\displaystyle \frac{\delta ^2L[\rho ,w]}{\delta \rho \delta \rho }}\right)\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right)_{(0,0,𝐳),(0,0,𝐳)}^1`$ (54)
which is shown graphically in Fig. 2.

Therefore the (38) can be written as

$`\dot{\mathrm{\Gamma }}[\rho ,w]`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}{\displaystyle \frac{\delta ^2W[\sigma ]}{\delta \sigma _{0,0,𝐱}}}+{\displaystyle \frac{1}{2e^2\lambda ^2}}\left\{i\mathrm{}\mathrm{T}r\left[\mathrm{\Delta }\left({\displaystyle \frac{\delta ^2W[\sigma ,j]}{\delta j\delta j}}K[\sigma ,j]\right)\right]\right\}`$ (55)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐳}\left[{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta \rho \delta \rho }}\right]_{(0,0,𝐳),(0,0,𝐳)}^1+{\displaystyle \frac{1}{2e^2\lambda ^2}}\left\{i\mathrm{}\mathrm{T}r\mathrm{\Delta }\left[\left({\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ,w]}{\delta w\delta w}}\right)^1L[\rho ,w]\right]\right\}.`$ (56)
The system of equations (47) and (55) is closed and describes the evolution of the interactive system.

### C $`𝒪(\rho ^4)`$ truncation

In the rest of this Section we return to the non-interacting case and consider two different approximations for the evolution equation. First we use the ansatz

$$\mathrm{\Gamma }[\rho ]=\mathrm{\Gamma }+\mathrm{\Gamma }_{\widehat{a}}\rho _{\widehat{a}}+\frac{1}{2}\mathrm{\Gamma }_{\widehat{a},\widehat{b}}\rho _{\widehat{a}}\rho _{\widehat{b}}+\frac{1}{3!}\mathrm{\Gamma }_{\widehat{a},\widehat{b},\widehat{c}}\rho _{\widehat{a}}\rho _{\widehat{b}}\rho _{\widehat{c}}+\frac{1}{4!}\mathrm{\Gamma }_{\widehat{a},\widehat{b},\widehat{c},\widehat{d}}\rho _{\widehat{a}}\rho _{\widehat{b}}\rho _{\widehat{c}}\rho _{\widehat{d}}$$
(57)
where the combined ’super-index’ $`\widehat{a}=(\alpha ,x_a)`$ has been introduced and summation/integration is assumed for double indices. The inverse of

$$\frac{\delta ^2\mathrm{\Gamma }[\rho ]}{\delta \rho _{\widehat{a}}\delta \rho _{\widehat{b}}}=(F^1)_{\widehat{a},\widehat{b}}+\mathrm{\Gamma }_{\widehat{a},\widehat{b}}^{(3)}[\rho ]+\mathrm{\Gamma }_{\widehat{a},\widehat{b}}^{(4)}[\rho ]$$
(58)
where

$`(F^1)_{\widehat{a},\widehat{b}}`$ $`=`$ $`\mathrm{\Gamma }_{\widehat{a},\widehat{b}}`$ (59)
$`\mathrm{\Gamma }_{\widehat{a},\widehat{b}}^{(3)}[\rho ]`$ $`=`$ $`\mathrm{\Gamma }_{\widehat{a},\widehat{b},\widehat{c}}\rho _{\widehat{c}}`$ (60)
$`\mathrm{\Gamma }_{\widehat{a},\widehat{b}}^{(4)}[\rho ]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }_{\widehat{a},\widehat{b},\widehat{c},\widehat{d}}\rho _{\widehat{c}}\rho _{\widehat{d}}`$ (61)
on the right hand side can be obtained by the Neumann-series,

$$\dot{\mathrm{\Gamma }}[\rho ]=\frac{g}{2}_𝐱\left(FF(\mathrm{\Gamma }^{(3)}[\rho ]+\mathrm{\Gamma }^{(4)}[\rho ])F+F(\mathrm{\Gamma }^{(3)}[\rho ]+\mathrm{\Gamma }^{(4)}[\rho ])F(\mathrm{\Gamma }^{(3)}[\rho ]+\mathrm{\Gamma }^{(4)}[\rho ])F\mathrm{}\right)_{(0,0,𝐱),(0,0,𝐱)}.$$
(62)
The identification of the coefficients of the same $`\rho `$ powers gives the equations

$`\dot{\mathrm{\Gamma }}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}F_{(0,0,𝐱),(0,0,𝐱)}`$ (63)
$`\dot{\mathrm{\Gamma }}_{\widehat{a}}\rho _{\widehat{a}}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}\left(F\mathrm{\Gamma }^{(3)}[\rho ]F\right)_{(0,0,𝐱),(0,0,𝐱)}`$ (64)
$`\dot{\mathrm{\Gamma }}_{\widehat{a},\widehat{b}}\rho _{\widehat{a}}\rho _{\widehat{b}}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}\left[F\left(\mathrm{\Gamma }^{(4)}[\rho ]\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]\right)F\right]_{(0,0,𝐱),(0,0,𝐱)}`$ (65)
$`\dot{\mathrm{\Gamma }}_{\widehat{a},\widehat{b},\widehat{c}}\rho _{\widehat{a}}\rho _{\widehat{b}}\rho _{\widehat{c}}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}\left[F\left(2\mathrm{\Gamma }^{(4)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]\right)F\right]_{(0,0,𝐱),(0,0,𝐱)}`$ (66)
$`\dot{\mathrm{\Gamma }}_{\widehat{a},\widehat{b},\widehat{c},\widehat{d}}\rho _{\widehat{a}}\rho _{\widehat{b}}\rho _{\widehat{c}}\rho _{\widehat{d}}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐱}[F(\mathrm{\Gamma }^{(4)}[\rho ]F\mathrm{\Gamma }^{(4)}[\rho ]\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(4)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]`$ (68)
$`2\mathrm{\Gamma }^{(4)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]+\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ]F\mathrm{\Gamma }^{(3)}[\rho ])F]_{(0,0,𝐱),(0,0,𝐱)},`$
shown on Fig. 3. It is easy to understand the structure of these equations. They express the principle that an infinitesimal change of the strength of the disorder, $`gg+\delta g`$, can be realized in the order $`𝒪(\delta g)`$ by visiting each correlation line of the graphs and change its coefficient, represented by the cross on the correlation lines of the graphs. The equations express this process in a hierarchical manner, in terms of the one particle irreducible vertices.

The resummation of the perturbation series, implicit in the solution of the evolution equation, can easily be seen by iteration,

$$W_\lambda ^{(n+1)}[\sigma ]=W_0^{(n)}[\sigma ]+\frac{g}{2}_0^\lambda 𝑑\lambda ^{}_𝐩\frac{\delta ^2W_\lambda ^{}^{(n)}[\sigma ]}{\delta \sigma _{0,\omega =0,p}\delta \sigma _{0,\omega =0,p}},$$
(69)
where each step increases the order of the interactions/correlations present by one and the integration over $`\lambda ^{}`$ inserts the necessary symmetry factors. In the numerical integration of the evolution equation each step $`\lambda \lambda +\delta \lambda `$ increases the order of the interactions/correlations by one and all graphs are resummed as $`\delta \lambda 0`$. Naturally this resummation is performed in the given truncation of the irreducible vertex functions in the momentum space. In particular, the the self-energy insertions, the ladder and crossed diagrams of a particle-hole loop  are resummed among other contributions by the first term on the right hand side of the third equation in (63).

### D Gradient expansion

We consider now the ansatz

$$\mathrm{\Gamma }[\rho ]=_x\left[\rho _{\alpha ,x}V_{\alpha ,\beta }(\rho _x,i_x)\rho _{\beta ,x}U(\rho _x)\right]$$
(70)
for the non-interactive case where the space-time derivatives act on the function $`\rho _{\beta ,x}`$ only,

$$V_{\alpha ,\beta }(\rho ,p)=\underset{A}{}V_{\alpha ,A,\beta }(\rho )p^A.$$
(71)
The combined index $`A=(\nu _0,\nu _1,\nu _2,\nu _3)`$ defines the partial derivative $`^A=_0^{\nu _0}_1^{\nu _1}_2^{\nu _2}_3^{\nu _3}`$. Since static disorder appears through the $`\rho _0`$ dependence in the evolution equation we retain the $`𝒪((\rho _j)^2)`$ terms only which are needed for the conductivity. Therefore we shall use the expression

$`\mathrm{\Gamma }[\rho ]`$ $`=`$ $`{\displaystyle _z}[{\displaystyle \frac{1}{2}}\rho _{0,z}\mathrm{\Gamma }^{tt}(\rho _{0,z},\mathrm{sign}(Imi_0)i_0,\mathrm{\Delta })\rho _{0,z}`$ (74)
$`+{\displaystyle \frac{1}{2}}\rho _{0,z}\mathrm{\Gamma }^{ts}(\rho _{0,z},\mathrm{sign}(Imi_0)i_0,\mathrm{\Delta })_k\rho _{k,z}+{\displaystyle \frac{1}{2}}\rho _{k,z}\mathrm{\Gamma }^{ts}(\rho _{0,z},\mathrm{sign}(Imi_0)i_0,\mathrm{\Delta })_k\rho _{0,z}`$
$`+{\displaystyle \frac{1}{2}}\rho _{j,z}(\mathrm{\Gamma }^T(\rho _{0,z},\mathrm{sign}(Imi_0)i_0,\mathrm{\Delta })T_{j,k}()+\mathrm{\Gamma }^L(\rho _{0,z},\mathrm{sign}(Imi_0)i_0,\mathrm{\Delta })L_{j,k}())\rho _{k,z}U(\rho _z)]`$
with

$$U(\rho )=\frac{1}{2}\mathrm{\Gamma }^t\rho _j^2+V(\rho _0)$$
(75)
and

$$L()=\frac{}{\mathrm{\Delta }},T()=1L().$$
(76)
The coefficient functions satisfy the relations $`\mathrm{\Gamma }^{ts}()=\mathrm{\Gamma }^{ts}()`$, $`\mathrm{\Gamma }^{tt}()=\mathrm{\Gamma }^{tt}()`$, $`\mathrm{\Gamma }^T()=\mathrm{\Gamma }^T()`$, $`\mathrm{\Gamma }^L()=\mathrm{\Gamma }^L()`$. The effective action contains derivatives in arbitrary high order the truncation is that they act on a single $`\rho `$ only. The computation of the second functional derivative used in the evolution equation is sketched in Appendix D.

The evolution of the functions parameterizing the effective action can the easiest be obtained by writing $`\rho _{\alpha ,x}=\rho _\alpha +\delta \rho _{\alpha ,x}`$ with $`_x\delta \rho _{\alpha ,x}=0`$ and expanding $`\mathrm{\Gamma }[\rho +\delta \rho ]`$ in $`\delta \rho `$ up to terms $`𝒪(\delta \rho ^2)`$. We shall use the notation $`_k^j=\rho _0_{\rho _0}^{j+1}+k_{\rho _0}^j`$ with the property $`_{\rho _0}_k^j=_{k+1}^{j+1}`$, the $`\gamma `$-functions $`\gamma _p^{tt}=\mathrm{\Gamma }^t+_1^0\mathrm{\Gamma }^{tt}(p)`$, $`\gamma _p^{ts}=_2^0\mathrm{\Gamma }^{ts}(p)`$, $`\gamma _p^L=\mathrm{\Gamma }^s+\mathrm{\Gamma }^L(p)`$, $`\gamma _p^T=\mathrm{\Gamma }^s+\mathrm{\Gamma }^T(p)`$, and write the left hand side of the evolution equation as

$$\dot{\mathrm{\Gamma }}[\rho ]+_p\left\{\frac{1}{2}\delta \rho _{0,p}\dot{\gamma }^{tt}\delta \rho _{0,p}+i\delta \rho _{0,p}\dot{\gamma }^{ts}(p)p_j\delta \rho _{j,p}+\frac{1}{2}\delta \rho _{j,p}\left[\dot{\gamma }^T(p)T_p+\dot{\gamma }^L(p)L_p\right]\delta \rho _{k,p}\right\}$$
(77)
The right hand side gives

$`{\displaystyle \frac{g}{2}}{\displaystyle _p}\delta _{p_0,0}V(\mathrm{\Gamma }^{(0)1})_{p,p}^{0,0}+{\displaystyle \frac{g}{4}}{\displaystyle _p}\delta _{p_0,0}(\mathrm{\Gamma }^{(0)1})_{p,p}^{0,\beta }{\displaystyle _r}\delta \rho _{\gamma ,r}\delta \rho _{\rho ,r}\mathrm{\Gamma }_{p,p,r,r}^{\beta ,\kappa ,\gamma ,\rho }(\mathrm{\Gamma }^{(0)1})_{p,p}^{\kappa ,0}`$ (78)
$`{\displaystyle \frac{g}{2}}{\displaystyle _p}\delta _{p_0,0}(\mathrm{\Gamma }^{(0)1})_{p,p}^{0,\beta }{\displaystyle _r}\delta \rho _{\gamma ,r}\delta \rho _{\sigma ,r}\mathrm{\Gamma }_{p,p+r,r}^{\beta ,\rho ,\gamma }(\mathrm{\Gamma }^{(0)1})_{pr,p+r}^{\rho ,\kappa }\mathrm{\Gamma }_{pr,p,r}^{\kappa ,\mu ,\sigma }(\mathrm{\Gamma }^{(0)1})_{p,p}^{\mu ,0}+𝒪(\delta \rho ^3),`$ (79)
where $`F_r^{\rho ,\kappa }`$, $`\mathrm{\Gamma }_{t,t,r,r}^{\beta ,\kappa ,\gamma ,\rho }`$ and $`\mathrm{\Gamma }_{r,t,rt}^{\kappa ,\mu ,\sigma }`$ are recorded in Appendix D. The second term resums the impurity correlations attached to the particle-hole loop, the self-energy insertions together with the vertex corrections. The third term contributes at vanishing energy only and will be neglected.

The $`𝒪(\delta \rho ^0)`$ part of the evolution equation is

$$\dot{U}=\frac{g}{2}_𝐩\mathrm{\Gamma }^{1(0,0)})_𝐩=\frac{g}{2}_𝐩\frac{\gamma _p^L}{\gamma _p^{tt}\gamma _p^L+\widehat{𝐩}^2(\gamma _p^{ts})^2}$$
(80)
The $`𝒪(\delta \rho ^2)`$ part gives

$`\dot{\gamma }_r^{tt}`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle _𝐩}{\displaystyle \frac{[_{\rho _0}^2\gamma _p^{tt}+_{\rho _0}^2\gamma _r^{tt}](\gamma _p^L)^2}{[\gamma _p^{tt}\gamma _p^L+\widehat{𝐩}^2(\gamma _p^{ts})^2]^2}}`$ (81)
$`\dot{\gamma }_r^{ts}`$ $`=`$ $`{\displaystyle \frac{g}{4}}_{\rho _0}^2\gamma _r^{ts}{\displaystyle _𝐩}{\displaystyle \frac{(\gamma _p^L)^2}{[\gamma _p^{tt}\gamma _p^L+\widehat{𝐩}^2(\gamma _p^{ts})^2]^2}}`$ (82)
$`\dot{\gamma }_r^T`$ $`=`$ $`{\displaystyle \frac{g}{4}}_0^2\gamma _r^T{\displaystyle _𝐩}{\displaystyle \frac{(\gamma _p^L)^2}{[\gamma _p^{tt}\gamma _p^L+\widehat{𝐩}^2(\gamma _p^{ts})^2]^2}}`$ (83)
$`\dot{\gamma }_r^L`$ $`=`$ $`{\displaystyle \frac{g}{4}}_0^2\gamma _r^L{\displaystyle _𝐩}{\displaystyle \frac{(\gamma _p^L)^2}{[\gamma _p^{tt}\gamma _p^L+\widehat{𝐩}^2(\gamma _p^{ts})^2]^2}}.`$ (84)
One can verify as a consistency check that the $`𝒪(r^0)`$ parts of the equations for $`\gamma ^L`$ and $`\gamma ^T`$ agree with (80). Eqs. (81) have the form of a one-dimensional diffusion equation with $`(\rho _0,g)`$ as ’space-time’ coordinates with a ’space’-dependent diffusion constant given by the one-loop integrals.

The evolution equations need two boundary conditions in $`\rho _0`$. One can be obtained by noting that the Green functions $`\stackrel{~}{D}^{\mu ,\nu }`$ are symmetrical under charge conjugation $`\mu 6\mathrm{}^2/ma^2\mu `$. Therefore $`_{\rho _0}\gamma =0`$ will be chosen as a boundary condition at half-filling, $`\mu =3\mathrm{}^2/ma^2`$, $`\rho _0=1/2a^3`$ and the density will be restricted into the interval $`0\rho _01/2a^3`$. The other boundary condition is imposed at the lowest density point which is chosen at sufficiently small chemical potential in such a manner that only one particle is found in the system. By assuming that the Green functions increase linearly with the density in this low density regime the boundary condition $`_{\rho _0}^2\gamma =2\gamma /\rho _0^2`$ will be imposed.

## IV Initial conditions

After having obtained the evolution equations (80) we need their initial conditions, imposed in the perturbative regime. One might choose either $`\lambda =0`$ and use the free effective action or take a small but non-vanishing $`\lambda `$ and compute the initial conditions perturbatively. For a computation carried out at vanishing frequency and on infinitely large system one has to start with small, non-vanishing strength of disorder. The reason is that we integrate a differential equation in this scheme and have to assume that the $`\lambda `$-dependence is differentiable. But non-interacting particles move ballistically in the absence of disorder and an arbitrarily weak disorder already induces important changes in the transport coefficients. In fact, the system becomes localised in one or two dimensions at arbitrarily weak disorder. Sufficiently dilute three dimensional systems should be localised at weak impurity field and the construction of the initial conditions at weak disorder is rather non-trivial .

In a more formal language one may say that the difference between the initial conditions imposed at $`\lambda =0`$ and $`\lambda 0`$ lies in the way time reversal symmetry is broken. In order to arrive at finite conductivity one needs the breakdown of the time reversal and the space-translation invariance. Neither of these breakdowns takes place in the absence of interactions between the electrons and their environment. The interaction with static impurities breaks translation invariance. The interactions break the time reversal invariance dynamically in infinite volume. The advantage of the perturbative initial conditions imposed at small but non-vanishing $`\lambda `$ is that they contain the desired symmetry breaking pattern.

But we have to bear in mind that the numerical solution of the evolution equation can be obtained for a finite system only. Furthermore, in order to avoid the singular chemical potential dependence in finite systems it is necessary to perform the computation at finite frequency and to study the limit $`\omega 0`$ numerically. The $`\lambda `$-dependence of the Green functions of a finite system or at finite frequency is regular at $`\lambda =0`$ and no difficulty is expected during the integration of the evolution equation. Therefore the initial conditions will be constructed in the absence of impurities, at $`\lambda =0`$. The time reversal invariance is then broken externally, by the non-vanishing frequency and the transition from the ballistic to the diffusive regime will be induced smoothly during the evolution.

It is worthwhile noting that the issue of localization of non-interacting particles is far from being a one-particle problem and the density-current functional is a highly non-trivial object even for non-interacting electrons without disorder, at $`\lambda =0`$. This is because the Pauli blocking of the occupied states represents a strong topological many-particle interaction which obeys no small parameter. This circumstance is crucial in expecting a singularity in the ground state in function of the disorder strength.

### A Current and density Green functions

The first step is the construction of the generator functional for the connected Green functions of the density and the current. For this end we write $`\sigma _{\alpha ,x}=\sigma _\alpha +\delta \sigma _{\alpha ,x}`$ where $`_x\delta \sigma _{\alpha ,x}=0`$ and find the following generator functional

$$W[\sigma +\delta \sigma ]=i\mathrm{}\mathrm{T}r\mathrm{log}[G^1+\sigma /+\delta \sigma /]=\mathrm{}\mathrm{T}r\mathrm{log}D^1+\rho ^{}\delta \sigma \frac{1}{2}\delta \sigma \stackrel{~}{D}\delta \sigma +𝒪(\delta \sigma ^3)$$
(85)
where $`D_q^1=G_q^1+\sigma /`$, $`\rho _{\alpha ,z}^{}=i\mathrm{}\mathrm{T}r\left[D𝒪_{\alpha ,z}\right]=\rho _\alpha ^{}`$, and $`\stackrel{~}{D}_{(\alpha ,x),(\beta ,y)}=i\mathrm{}\mathrm{T}r\left[D𝒪_{\alpha ,x}D𝒪_{\beta ,x}\right]`$. Due to rotational invariance one has the form

$$\stackrel{~}{D}_{\omega ,𝐪}^{\mu ,\nu }=\left(\begin{array}{cc}\stackrel{~}{D}^{tt}(\omega ,𝐪^2)& i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\stackrel{~}{D}^{ts}(\omega ,𝐪^2)\\ i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\stackrel{~}{D}^{ts}(\omega ,𝐪^2)& \stackrel{~}{D}^{j,k}(\omega ,𝐪)\end{array}\right)$$
(86)
with $`\stackrel{~}{D}^{j,k}(\omega ,𝐪)=T^{j,k}(𝐪)\stackrel{~}{D}^T(\omega ,𝐪^2)+L^{j,k}(𝐪)\stackrel{~}{D}^L(\omega ,𝐪^2)`$, $`\widehat{𝐪}=\frac{2}{a}\mathrm{sin}\frac{a𝐪}{2}`$, $`L=\widehat{𝐪}\widehat{𝐪}/\widehat{𝐪}^2`$ and $`T=1L`$. But on a lattice with finite lattice spacing there is no rotational symmetry and the transverse piece becomes more involved. This complication will be neglected below. A more detailed expression for the propagator $`\stackrel{~}{D}^{\mu ,\nu }`$ is presented in Appendix E.

### B Variable transformation

The second step in constructing the density-current functional for the initial condition is the change of variables $`\sigma \rho `$. For this end we have to invert the relation

$$\rho _{\alpha ,x}=(\rho _\alpha +\delta \rho _{\alpha ,x})=i\mathrm{}\mathrm{T}r[\frac{1}{D^1+\sigma /+\delta \sigma /}𝒪_{\alpha ,x}]=i\mathrm{}\mathrm{T}r[D𝒪_{\alpha ,x}]+i\mathrm{}\mathrm{T}r[D\delta \sigma /D𝒪_{\alpha ,x}]+𝒪(\rho ^2).$$
(87)
We write $`\mu +\rho _0\mu `$ set $`\rho _j=\sigma _j=0`$ in the absence of Coulomb interaction and find

$$\delta \sigma =\stackrel{~}{D}^1\delta \rho $$
(88)
### C Effective action

Finally, the effective action can be written as

$$\mathrm{\Gamma }[\rho ]=i\mathrm{}\mathrm{T}r\mathrm{log}D^1+\sigma \rho \frac{1}{2}\delta \rho \stackrel{~}{D}^1\delta \rho +𝒪(\delta \rho ^3).$$
(89)
The identification of the part $`𝒪(\delta \rho ^0)`$ yields

$$U=\sigma _\alpha \rho _\alpha \frac{i\mathrm{}}{V}\mathrm{T}r\mathrm{log}D^1.$$
(90)
The $`𝒪(\delta \rho ^2)`$ order, the second functional derivative gives

$$\left(\begin{array}{cc}\gamma _q^{tt}& i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\gamma _q^{ts}\\ i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\gamma _q^{ts}& T\gamma _q^T+L\gamma _q^L\end{array}\right)=\left(\begin{array}{cc}\stackrel{~}{D}^{tt}& i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\stackrel{~}{D}^{ts}\\ i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\stackrel{~}{D}^{ts}& T\stackrel{~}{D}_T+L\stackrel{~}{D}_L\end{array}\right)^1=\left(\begin{array}{cc}\frac{\stackrel{~}{D}_L}{D}& i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\frac{\stackrel{~}{D}^{ts}}{D}\\ i\widehat{𝐪}e^{\frac{i}{2}a𝐪}\frac{\stackrel{~}{D}^{ts}}{D}& \frac{1}{\stackrel{~}{D}_T}T+\frac{\stackrel{~}{D}^{tt}}{D}L\end{array}\right),$$
(91)
where

$$D(\omega ,𝐪^2)=\stackrel{~}{D}^{tt}(\omega ,𝐪^2)\stackrel{~}{D}^L(\omega ,𝐪^2)+\widehat{𝐪}^2(\stackrel{~}{D}^{ts}(\omega ,𝐪^2))^2.$$
(92)
Therefore

$`\gamma ^{tt}(\omega ,𝐪^2)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{D}^L(\omega ,𝐪^2)}{D(\omega ,𝐪^2)}}`$ (93)
$`\gamma ^T(\omega ,𝐪^2)`$ $`=`$ $`{\displaystyle \frac{1}{D^T(\omega ,𝐪^2)}}`$ (94)
$`\gamma ^L(\omega ,𝐪^2)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{D}^{tt}(\omega ,𝐪^2)}{D(\omega ,𝐪^2)}}`$ (95)
$`\gamma ^{ts}(𝐪^2)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{D}^{ts}(\omega ,𝐪^2)}{D(\omega ,𝐪^2)}}.`$ (96)
The Euclidean Ward identities require

$$\stackrel{~}{D}^{ts}(\omega ,𝐪)=\frac{1}{\omega }\left(\stackrel{~}{D}^{j,k}(\omega ,𝐪)+\frac{1}{m}\delta ^{j,k}\stackrel{~}{\rho }^{}\right)$$
(97)
and

$$\stackrel{~}{D}^{tt}(\omega ,𝐪^2)=\frac{\widehat{𝐪}^2}{\omega }\stackrel{~}{D}^{ts}(\omega ,𝐪^2),$$
(98)
c.f. Appendix C. It is easy to see that these equations restrict $`\gamma _p^L`$ to be momentum independent.

## V Numerical results

There are two points where the numerical integration can be made faster. First, the one-loop integral in the evolution equation requires the knowledge of the $`\gamma `$ functions in the whole Brillouin zone. It was assumed that the $`\gamma `$ functions depend on $`\widehat{𝐩}^2`$, momentum square only. This assumption was found to be acceptable in the initial condition. Note in this respect that $`\stackrel{~}{D}^L`$ and $`\gamma ^L`$ remain well defined in a non-isotrop Brillouine zone. The other improvement on the speed of the algorithm is the runtime adjustment of the step size $`\mathrm{\Delta }g`$ during the numerical integration. This was carried out by keeping the largest relative increment in the differential equations under $`1\%`$. Lattices up to size $`320^3`$ were used but the results presented below correspond to the volume $`80^3`$ unless it is stated otherwise. The units $`a=m=\mathrm{}=1`$ were used to express dimensional quantities. The integration of the evolution equation started at $`g=0`$ with $`\mathrm{\Delta }g=10^6`$ and ended at $`g10^5`$, where $`\mathrm{\Delta }g10^4`$ was used. 60 division points were used to monitor the density dependence of the $`\gamma `$ functions in the interval $`2\times 10^6<\rho _0<0.5`$. The Euclidean frequency was restricted by the bounds $`10^5<\omega <335`$.

The Ward identities allow us to express the longitudinal conductivity as

$$\frac{\sigma (\omega ,𝐩^2)}{ie^2}=\stackrel{~}{D}^{ts}(\omega ,𝐩^2)=\frac{\omega }{𝐩^2\gamma ^{tt}(\omega ,𝐩^2)+\omega ^2\gamma ^L(\omega )},$$
(99)
indicating that the d.c. conductivity in the homogeneous $`𝐩=0`$ limit corresponds to an I.R. cut-off, a ’mass gap’ in the propagator $`\stackrel{~}{D}^{ts}(\omega ,𝐩^2)`$ and it can only be non-vanishing due to an $`𝒪(1/𝐩^2)`$ singularity in $`\gamma ^{tt}`$ functions ($`\gamma ^L`$ remains finite as $`\omega 0`$). A rough, qualitative picture of what happens as the disorder is increased is given by the analogy with the diffusion equation. The $`\gamma `$-functions of the initial condition have a singular peak at $`\rho _00`$. Therefore the diffusive evolution smears this peak out and tries to ’fill up’ the $`\gamma `$ functions at large densities where they take smaller values. Since the propagator $`\stackrel{~}{D}^{tt}`$ is $`𝒪(1/V)`$ for a system of finite volume $`V`$ the diffusion is rather slow at low densities. Therefore the main $`g`$-dependence will be an increase of the $`\gamma `$ functions for intermediate and large densities. As this increase diffuses ’up’ in the density it increases the 1PI functions and suppresses the conductivity.

Let us start with the $`g`$-dependence of $`\sigma ^1(\omega ,𝐩^2)`$, shown in Fig. 4 (a) at $`\omega =2\times 10^5`$ and $`𝐩=𝐩_1=(0,0,\pi /80)`$ for densities $`7.6\times 10^3<\rho _0<0.5`$. The inverse conductivity displays a plateau at weak disorder and starts to increase at $`g0.1`$ only. In order to see clearer what happens at this crossover point $`\mathrm{\Delta }\sigma ^1(g)=\sigma ^1(g)\sigma ^1(0)`$, the contribution due to the disorder, is shown in Fig. 4 (b). For some densities one finds a slight drop in $`\sigma ^1`$ at the beginning of the evolution, the figure shows the values with $`\mathrm{\Delta }\sigma ^1(g)>0`$ only. We believe that this weak non-monotonic $`g`$-dependence is a finite size effect. It is induced by the evolution equation when the density dependence in the initial conditions is distorted by a finite size effect. The latter is due to the constraint $`\mathrm{\Theta }(E_𝐪)`$ in the integrand of the zero temperature Green functions which induces an artificial step function-like chemical potential dependence in the initial conditions, (E17). When this finite size effect is ignored the lesson of Fig. 4 (b) is that the change of the nature of $`\sigma ^1`$ at $`g1`$ in Fig. 4 (a) is not a singularity. What happens is that the increase of $`\sigma ^1`$ due to disorder becomes comparable with the initial value and disorder starts to play a more important role. $`\mathrm{\Delta }\sigma ^1(g)`$ plotted for different values of the momentum $`𝐩`$ shows similar behavior. One finds the scaling $`\mathrm{\Delta }\sigma ^1(g)g^k`$ with $`k=1`$ for $`g<<1`$ and $`k=1/2`$ as $`g>>1`$ without any non-analicity in the $`g`$ dependence. The behavior $`\mathrm{\Delta }\sigma ^1(g)g`$ characterizes weak disorder. It is easy to see that the condition of the simplification that the $`g`$-dependence generated by the evolution equations (81) factorizes as an overall multiplicative factor is just the asymptotic scaling $`\mathrm{\Delta }\sigma ^1(g)\sqrt{g}`$. The possibility of any sudden transition at $`\mathrm{p}=0`$ developing in the thermodynamical and the $`\omega 0`$ limit between the weak and strong disorder region is still hidden in the momentum dependence, the issue we turn to now.

In order to explore the I.R. regime the extrapolation of the conductivity to $`𝐩=0`$ was performed by fitting the $`𝐩^2`$ dependence of $`\sigma ^1(\omega ,𝐩^2)`$ by a polynomial up to $`𝐩^{12}`$. The result is dominated by finite size effects at weak disorder. In fact, $`g=0`$, a perfectly ordered system should have vanishing $`\sigma ^1`$ at $`𝐩=0`$, i.e. long range Green function $`\stackrel{~}{D}^{ts}`$. In more physical terms, the density of state is a sum of Dirac-delta peaks for finite systems and one needs rather large volume in order to approach a continuous density of state which is needed for the Bragg reflections on the lattice and for finite conductivity. The dependence of the inverse conductivity on the strength of disorder, depicted in Fig. 5 (a) supports this expectation. Furthermore one can see that at low densities where the density of states is lower the conductivity is strongly suppressed and $`\sigma ^1`$ is large and positive for any $`g`$. When the negative parts of $`\sigma ^1`$ are discarded as in Fig. 5 (b) one sees the sudden drop of the conductivity as the strength of disorder is increased. $`\sigma ^1`$ is plotted against the density in Fig. 6 showing that the density dependence of the a.c. conductivity is strongly suppressed for strong disorder.

For the identification of the sudden increase of $`\sigma ^1`$ with the onset of the localized phase we have to perform the limit $`\omega 0`$. The frequency dependence of $`\sigma ^1`$ is shown in Fig. 7 (a). One recovers the usual $`\sigma 1/\omega `$ behavior at any momentum $`𝐩`$ at high frequencies except a short flattening around $`\omega 10`$ where the disorder decouples. In fact, the evolution is suppressed at high frequencies due to the smallness of propagators in the initial conditions, c.f. (E17), and the scaling $`\sigma 1/\omega `$ is recovered at higher frequencies with a proportionality constant given approximately by the initial conditions. As we turn towards low frequencies we find two qualitatively different behaviors for the conductivity extrapolated to $`𝐩=0`$. At weaker disorder $`\sigma ^1`$ decreases and becomes negative, more precisely is dominated by finite size effects and should not be taken seriously when computed at the present size lattices. On the contrary, for stronger disorder the computation is free of finite size effects and $`\sigma `$ decreases as $`\omega 0`$. Motivated by the Mott form

$$\sigma _f=C\left(\frac{\omega }{\omega _0}\right)^\kappa \left|\mathrm{log}\frac{\omega }{\omega _0}\right|^\eta $$
(100)
with $`(\kappa ,\eta )=(2,4)`$, the values of $`\kappa `$ $`\eta `$, $`\omega _0`$ and $`C`$ were fitted in Eq. (100) for strong disorder. Mott’s result alone compares rather poorly with the numerical results as shown in Fig. 7 (b). The simplest choice, $`(\kappa ,\eta )=(1,0)`$ represents a better approximation but the fit with the results $`(\kappa ,\eta )=(1.57,3.18)`$ shows the presence of logarithmic corrections and gives a further improvement.

There is a discontinuity at $`g=g_{cr}`$ separating finite and diverging values of $`\sigma ^1`$ at weak and strong disorder, respectively. The critical strength of the disorder, $`g_{cr}`$ can be estimated by finding where $`\sigma ^1`$ changes sign first as $`g`$ is decreased, c.f. Fig. 8. This is how the localization-delocalization transition appears in this scheme.

## VI Conclusions

A new, non-perturbative scheme is presented in this paper for the description of disordered system. The basic idea is to express the effects of an infinitesimal increase of the disorder by changing the density. No small parameter is required to obtain the evolution equation but truncation is needed to convert it into a useful numerical algorithm. We believe that the gradient expansion is a well suited and consistent truncation scheme for the computation of the transport coefficients.

There is a rather detailed space-time picture for weak localization drawn from the partial resummation of the perturbation expansion . The mechanism proposed in this paper is more formal, it is a diffusion process where the density and disorder strength are considered as space and time variables. The 1PI functions take large values at low density in the initial conditions. The conductivity decreases when the disorder is becoming stronger due to the diffusion of the low density peak of the 1PI functions towards higher densities. Localization appears as a special, correlated increase of the 1PI functions at different momentum values.

The numerical results obtained in solving the evolution equation indicates that weak localization is rather difficult to establish due to strong finite size effects. In fact, the finite conductivity arises from an $`𝒪(1/𝐩^2)`$ I.R. singularity of a 1PI function whose reproduction requires unusually large lattices. But strong disorder suppresses the long range correlations and renders the numerical solution more reliable.

It was found that the conductivity is vanishing as $`\omega 0`$ when the strength of disorder exceeds a certain threshold. But the fit of the numerical results yields fractional powers of $`\omega `$ and $`\mathrm{log}\omega `$, indicating that several terms contribute in the frequency interval considered. It is natural to identify this threshold with the localization transition.

The results presented here are preliminary, obviously more careful study of the finite size effects is needed to establish contact with perturbation expansion and to identify the weak localization regime. Furthermore, a more extended frequency interval must be studied in order to identify the low frequency scaling form of the conductivity in the localized regime. Finally, the most important step is the inclusion of the Coulomb interaction in the numerical solution.

## ACKNOWLEDGMENTS

I thank János Hajdu, Shimon Levit and Peter Wölfle for illuminating discussions and encouragement.

## A Notations

We summarize here briefly the notations used in the paper. The space-time integrals and the Fourier transform are

$$_{T/2}^{T/2}𝑑t=_t,d^dx=a^d\underset{x}{}=_x,_xf_xg_x=fg,\frac{d^dp}{(2\pi )^d}=\frac{1}{V}\underset{p}{}=_p,f_p=_xe^{ipx}f_x,f_x=_pe^{ipx}f_p,$$
(A1)
with $`x=(𝐱,t)`$, $`p=(𝐩,\omega )`$, and $`xp=\mathrm{𝐱𝐩}t\omega `$. The Dirac-deltas $`\delta _{x,y}=a^da_\tau \delta _{x,y}^K`$ and $`\delta _{p,q}=VT\delta _{p,q}^K`$ are expressed in terms of the Kronecker-deltas $`\delta _{x,y}^K`$ and $`\delta _{p,q}^K`$. Identities for the functional derivatives are

$`{\displaystyle \frac{\delta }{\delta f_x}}`$ $`=`$ $`{\displaystyle \frac{}{a^3a_\tau f_x}}={\displaystyle \underset{p}{}}e^{ipx}{\displaystyle \frac{}{f_p}}={\displaystyle \frac{1}{VT}}{\displaystyle \underset{p}{}}e^{ipx}{\displaystyle \frac{\delta }{\delta f_p}}={\displaystyle _p}e^{ipx}{\displaystyle \frac{\delta }{\delta f_p}},`$ (A2)
$`{\displaystyle \frac{\delta }{\delta f_p}}`$ $`=`$ $`VT{\displaystyle \frac{}{f_p}}={\displaystyle \underset{x}{}}e^{ipx}{\displaystyle \frac{}{f_x}}={\displaystyle _x}e^{ipx}{\displaystyle \frac{\delta }{\delta f_x}}.`$ (A3)
Greek and Latin indices will be used for the range $`0,1,2,3`$ and $`1,2,3`$, respectively.

## B Annealed v.s. quenched evolution

It is shown that the annealed and the quenched averages correspond to simple evolution equation. Let us take a function $`f(w)`$ and introduce the generator functional

$$f(W[\sigma ])=\frac{𝒟[v]e^{S_d[v]}f(W_0[\sigma iv])}{𝒟[v]e^{S_d[v]}}$$
(B1)
where

$$S_d[v]=\frac{1}{2}vKv.$$
(B2)
We make the replacement $`vv\sigma `$ $`(M)`$, $`vvi\sigma `$ $`(E)`$, i.e. use the equation of motion for $`v`$ to write

$$f(W[\sigma ])=\frac{𝒟[v]e^{S_d[vi\sigma ]}f(W_0[iv])}{𝒟[v]e^{S_d[v]}}.$$
(B3)
We therefore have the relations

$`\dot{W}[\sigma ]f^{}(W[\sigma ])`$ $`=`$ $`{\displaystyle \frac{𝒟[v](vi\sigma )\dot{K}(vi\sigma )e^{\frac{1}{2}(vi\sigma )K(vi\sigma )}f(W_0[iv])}{2𝒟[v]e^{\frac{1}{2}vKv}}}+{\displaystyle \frac{1}{2}}f(W[\sigma ])\mathrm{T}r\dot{K}K^1,`$ (B4)
$`{\displaystyle \frac{\delta W[\sigma ]}{\delta \sigma _x}}f^{}(W[\sigma ])`$ $`=`$ $`{\displaystyle \frac{𝒟[v]\frac{\delta S_d[vi\sigma ]}{\delta \sigma _x}e^{S_d[v\sigma ]}f(W_0[iv])}{𝒟[v]e^{S_d[v]}}},`$ (B5)
$`{\displaystyle \frac{\delta }{\delta \sigma _y}}\left({\displaystyle \frac{\delta W^M[\sigma ]}{\delta \sigma _x}}f^{}(W^E[\sigma ])\right)`$ $`=`$ $`{\displaystyle \frac{𝒟[v]\left((K(vi\sigma ))_x(K(vi\sigma ))_yK_{x,y}\right)e^{\frac{1}{2}(vi\sigma )K(vi\sigma )}f(W_0^E[iv])}{𝒟[v]e^{\frac{1}{2}vKv}}}`$ (B6)
which allow us to write the evolution equation as

$$\dot{W}[\sigma ]=\frac{1}{2}\mathrm{T}r\dot{K}^1\left(\frac{\delta ^2W[\sigma ]}{\delta \sigma \delta \sigma }+\frac{\delta W[\sigma ]}{\delta \sigma }\frac{f^{\prime \prime }(W[\sigma ])}{f^{}(W[\sigma ])}\frac{\delta W[\sigma ]}{\delta \sigma }\right).$$
(B7)
The simplest equation arises from the assumption

$$f^{\prime \prime }(W[\sigma ])=\frac{n}{\mathrm{}}f^{}(W[\sigma ])$$
(B8)
where $`n`$ is a constant. The solution,

$$f(x)=\{\begin{array}{cc}x\hfill & n=0\hfill \\ e^{nx/\mathrm{}}\hfill & n0\hfill \end{array}$$
(B9)
corresponds to the quenched average for $`n=0`$ and to $`n`$ annealed replicas when $`n0`$. The evolution equation reads as

$$\dot{W}_q[\sigma ]=\frac{1}{2}\mathrm{T}r\dot{K}^1\left(\frac{\delta ^2W_a[\sigma ]}{\delta \sigma \delta \sigma }+\frac{n}{\mathrm{}}\frac{\delta W_a[\sigma ]}{\delta \sigma }\frac{\delta W_a[\sigma ]}{\delta \sigma }\right)$$
(B10)
One finds

$$\dot{\mathrm{\Gamma }}[\sigma ]=\frac{1}{2}\mathrm{T}r\dot{K}^1\left[\left(\frac{\delta ^2\mathrm{\Gamma }_a^E[\rho ]}{\delta \rho \delta \rho }\right)^1+\frac{n}{\mathrm{}}\rho \rho \right]$$
(B11)
in terms of the effective action $`\mathrm{\Gamma }[\sigma ]`$. The condition for the applicability of the replica method is the assumption that the solution of this differential equation not only converges as $`n0`$ but the $`n`$-dependence is continuous, i.e. the limit agrees with the solution obtained with $`n=0`$. According to Eq. (B10) $`n`$ should approaches zero at least with the inverse volume in the thermodynamical limit. This requirement, well known from elementary considerations, may pose problem in phase transitions or when some symmetry is broken spontaneously.

## C Ward identity

The Kubo formula for the conductivity describes the way external electromagnetic potential generates the motion of the charges. One needs special care in using this formalism when the charges in question are treated in the second quantized formalism. The point is that one should not use approximations which violate the particle number. In fact, an uncontrolled number of particles created by an approximation renders the transport coefficients divergent when the IR limit is considered, whatever weak error is committed in the approximation. It is easy to control the total number of particles by global symmetries. But the complications start when the particles carry a conserved charge which can take positive as well negative values. In this case the conservation of the algebraic sum of the charge is not enough any more to protect the transport coefficient against IR divergences. In fact, imagine an approximation scheme which though preserves the total electric charge exactly but leads to the creation of an uncontrolled number of neutral pairs. The conductivity will obviously be IR divergent in such a scheme due to the polarization cloud generated by the erroneous approximation. Charged particles are usually handled by gauge theories and the approximations which leads to the creation of neutral pairs without violating the charge conservation are those which are invariant under global symmetry transformations but violate local gauge transformations.

The best way to make sure that a symmetry is preserved is to check if the equation of motion in the given regularization and approximation scheme remains symmetrical. In the derivation of the equation of motion, (13), then one uses infinitesimal shift of the dynamical field variables which corresponds to a symmetry transformation. The action remains invariant and any contribution to the equation of motion must come from non-symmetrical terms, i.e. non-symmetrical parts of the dynamics are the only sources contributing to such a special combination of the equation of motion, called Ward identity.

In our case we generate the shift of the field variables by an infinitesimal gauge transformation,

$$\psi _xe^{\frac{i}{\mathrm{}}ϵ_x}\psi _x.$$
(C1)
Since the electromagnetic fields appear as an external source only without dynamics they do not participate in this transformation. Nevertheless they appear in the discussion because the transformation (C1) is equivalent with

$$U_\mu (x)U_\mu (x)e^{\frac{i}{\mathrm{}}a_\mu _\mu ^+ϵ_x}$$
(C2)
where

$$_\mu ^\pm f_x=\pm \frac{f_{x\pm \widehat{\mu }}f_x}{a_\mu }$$
(C3)
as long as the dynamics is gauge invariant. The transformation (C3) induces the change

$`\delta S`$ $`=`$ $`a^3a_\tau {\displaystyle \underset{x}{}}_\mu ^+ϵ_x\psi ^{}𝒪_{\mu ,x}\psi +𝒪(ϵ^2),`$ (C4)
$`\delta \psi ^{}𝒪_{0,x}\psi `$ $`=`$ $`{\displaystyle \frac{ia_\tau }{\mathrm{}}}_0^+ϵ_x\psi ^{}𝒪_{0,x}\psi +𝒪(ϵ^2),`$ (C5)
$`\delta \psi ^{}𝒪_{j,x}\psi `$ $`=`$ $`{\displaystyle \frac{1}{m}}_j^+ϵ_x\psi ^{}\stackrel{~}{𝒪}_{j,x}\psi +𝒪(ϵ^2),`$ (C6)
$`\delta \psi ^{}\stackrel{~}{𝒪}_{j,x}\psi `$ $`=`$ $`{\displaystyle \frac{ma^2}{\mathrm{}^2}}_j^+ϵ_x\psi ^{}𝒪_{j,x}\psi +𝒪(ϵ^2).`$ (C7)
Using this transformation as a reparametrization of the functional integral including the operators $`𝒪_\mu `$ and $`\stackrel{~}{𝒪}_\mu `$ with a corresponding source we find

$$0=\underset{x}{}\left\{_0^+ϵ_x\left(1+\frac{ia_\tau }{\mathrm{}}\sigma _{0,x}\right)\frac{\delta W_0[\sigma ,\stackrel{~}{\sigma }]}{\psi _{0,x}}+\underset{j}{}_j^+ϵ_x\left[\left(1+\frac{a^2}{\mathrm{}^2}\stackrel{~}{\sigma }_{j,x}\right)\frac{\delta W_0[\sigma ,\stackrel{~}{\sigma }]}{\psi _{j,x}}\frac{1}{m}\sigma _{j,x}\frac{\delta W_0[\sigma ,\stackrel{~}{\sigma }]}{\delta \stackrel{~}{\sigma }_{j,x}}\right]\right\},$$
(C8)
After integration in part, omitting the boundary terms and averaging over the disorder we have

$`0`$ $`=`$ $`{\displaystyle }𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}\{_0^{}\left[(1+{\displaystyle \frac{ia_\tau }{\mathrm{}}}(\sigma _{0,x}iv_𝐱)){\displaystyle \frac{\delta W_0[\sigma _0iv,\sigma _j,\stackrel{~}{\sigma }]}{\psi _{0,x}}}\right]`$ (C10)
$`+{\displaystyle \underset{j}{}}_j^{}[(1+{\displaystyle \frac{a^2}{\mathrm{}^2}}\stackrel{~}{\sigma }_{j,x}){\displaystyle \frac{\delta W_0[\sigma _0iv,\sigma _j,\stackrel{~}{\sigma }]}{\psi _{j,x}}}{\displaystyle \frac{1}{m}}\sigma _{j,x}{\displaystyle \frac{\delta W_0[\sigma _0iv,\sigma _j,\stackrel{~}{\sigma }]}{\delta \stackrel{~}{\sigma }_{j,x}}}]\}.`$
By taking the functional derivative of the ’equation of motion’ for the impurity averaging,

$$0=𝒟[v]e^{\frac{1}{2g}_𝐱v_𝐱^2}\left(i_t\frac{\delta W_0[\sigma _0iv,\sigma _j,\stackrel{~}{\sigma }]}{\psi _{0,t,𝐱}}+\frac{1}{g}v_𝐱W_0[\sigma _0+v,\sigma _j,\stackrel{~}{\sigma }]\right)$$
(C11)
we can write find the Ward identity

$$0=_0^{}\left[\left(1+\frac{ia_\tau }{\mathrm{}}\sigma _{0,x}\right)\frac{\delta W[\sigma ,\stackrel{~}{\sigma }]}{\psi _{0,x}}\frac{ga_\tau }{\mathrm{}}_t\frac{\delta ^2W[\sigma ,\stackrel{~}{\sigma }]}{\psi _{0,t,𝐱}\psi _{0,x}}\right]+\underset{j}{}_j^{}\left[\left(1+\frac{a^2}{\mathrm{}^2}\stackrel{~}{\sigma }_{j,x}\right)\frac{\delta W[\sigma ,\stackrel{~}{\sigma }]}{\psi _{j,x}}\frac{1}{m}\sigma _{j,x}\frac{\delta W[\sigma ,\stackrel{~}{\sigma }]}{\delta \stackrel{~}{\sigma }_{j,x}}\right].$$
(C12)
By assuming the form

$$W[\sigma ]=W_0+\sigma \rho ^{}\frac{1}{2}\sigma \stackrel{~}{D}\sigma +\frac{1}{3!}W_{\stackrel{~}{\alpha },\stackrel{~}{\beta },\stackrel{~}{\gamma }}\sigma _{\stackrel{~}{\alpha }}\sigma _{\stackrel{~}{\beta }}\sigma _{\stackrel{~}{\gamma }}+𝒪(\sigma ^4)$$
(C13)
the first functional derivative of (C12) with respect to $`\sigma `$ taken at $`\sigma =\stackrel{~}{\sigma }=0`$ yields the conserved currents

$$J_{0,x}=\stackrel{~}{D}_{y,x}^{\kappa ,0}\frac{ia_\tau }{\mathrm{}}\delta ^{\kappa ,0}\delta _{y,x}\rho ^{}+\frac{ga_\tau }{\mathrm{}}_tW_{(0,t,𝐱),(0,x),(\kappa ,y)},J_{j,x}=\stackrel{~}{D}_{y,x}^{\kappa ,j}+\frac{1}{m}\delta ^{\kappa ,j}\delta _{y,x}\stackrel{~}{\rho }_j^{}$$
(C14)
for arbitrary $`\kappa `$ and $`y`$. For U.V. finite model or when the U.V. divergences are logarithmic only the terms $`𝒪(a_\tau )`$ can be neglected.

## D Functional derivatives

We give few details of the computation of the second functional derivative of the effective action (70). The effective action is supposed to be the sum of terms like

$$\gamma [\rho ]=_xV_{A,\mu }(\rho _x)^A\rho _{\mu ,x}.$$
(D1)
and the first four functional derivatives are the following

$`{\displaystyle \frac{\delta \gamma [\rho ]}{\delta \rho _{\alpha ,a}}}`$ $`=`$ $`{\displaystyle _z}\left[\delta _{a,z}_{\rho _\alpha }V_{A,\mu }^A\rho _{\mu ,z}+V_{A,\alpha }^A\delta _{a,z}\right]`$ (D2)
$`{\displaystyle \frac{\delta ^2\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}}}`$ $`=`$ $`{\displaystyle _z}\left[\delta _{a,z}\delta _{b,z}_{\rho _\beta }_{\rho _\alpha }V_{A,\mu }^A\rho _{\mu ,z}+\delta _{a,z}_{\rho _\alpha }V_{A,\beta }^A\delta _{b,z}+\delta _{b,z}_{\rho _\beta }V_{A,\alpha }^A\delta _{a,z}\right]`$ (D3)
$`{\displaystyle \frac{\delta ^3\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}\delta _{c,z}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{A,\mu }^A\rho _{\mu ,z}+\delta _{a,z}\delta _{b,z}_{\rho _\beta }_{\rho _\alpha }V_{A,\gamma }^A\delta _{c,z}`$ (D5)
$`+\delta _{a,z}\delta _{c,z}_{\rho _\gamma }_{\rho _\alpha }V_{A,\beta }^A\delta _{b,z}+\delta _{b,z}\delta _{c,z}_{\rho _\gamma }_{\rho _\beta }V_{A,\alpha }^A\delta _{a,z}]`$
$`{\displaystyle \frac{\delta ^4\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}\delta \rho _{ϵ,e}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}\delta _{c,z}\delta _{e,z}_{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{A,\mu }^A\rho _{\mu ,z}+\delta _{a,z}\delta _{b,z}\delta _{c,z}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{A,ϵ}^A\delta _{e,z}`$ (D8)
$`+\delta _{a,z}\delta _{b,z}\delta _{e,z}_{\rho _ϵ}_{\rho _\beta }_{\rho _\alpha }V_{A,\gamma }^A\delta _{c,z}+\delta _{a,z}\delta _{c,z}\delta _{e,z}_{\rho _ϵ}_{\rho _\gamma }_{\rho _\alpha }V_{A,\beta }^A\delta _{b,z}`$
$`+\delta _{b,z}\delta _{c,z}\delta _{e,z}_{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }V_{A,\alpha }^A\delta _{a,z}]`$
We write $`V_{A,\sigma }=\rho _\kappa V_{\kappa ,A,\sigma }`$ and find

$`_{\rho _\alpha }V_{A,\sigma }`$ $`=`$ $`V_{\alpha ,A,\sigma }+\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\sigma }`$ (D10)
$`_{\rho _\beta }_{\rho _\alpha }V_{A,\sigma }`$ $`=`$ $`_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _\alpha }V_{\beta ,A,\sigma }+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma }`$ (D11)
$`_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{A,\sigma }`$ $`=`$ $`_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,\sigma }+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,\sigma }+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma }`$ (D12)
$`_{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{A,\sigma }`$ $`=`$ $`_{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _ϵ}_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,\sigma }+_{\rho _ϵ}_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,\sigma }`$ (D14)
$`+_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{ϵ,A,\sigma }+\rho _\kappa _{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma }`$
and

$`{\displaystyle \frac{\delta \gamma [\rho ]}{\delta \rho _{\alpha ,a}}}`$ $`=`$ $`{\displaystyle _z}\left[\delta _{a,z}V_{\alpha ,A,\sigma }^A\rho _{\sigma ,z}+\delta _{a,z}\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\sigma }^A\rho _{\sigma ,z}+\rho _\kappa V_{\kappa ,A,\alpha }^A\delta _{a,z}\right]`$ (D15)
$`=`$ $`{\displaystyle _z}\left[\delta _{a,z}V_{\alpha ,A,\sigma }^A\rho _{\sigma ,z}+\delta _{a,z}\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\sigma }^A\rho _{\sigma ,z}+(1)^{|A|}\delta _{a,z}^A(\rho _\kappa V_{\kappa ,A,\alpha })\right]`$ (D16)
$`{\displaystyle \frac{\delta ^2\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}(_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _\alpha }V_{\beta ,A,\sigma }+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma })^A\rho _{\sigma ,z}`$ (D18)
$`+\delta _{a,z}(V_{\alpha ,A,\beta }+\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\beta })^A\delta _{b,z}+\delta _{b,z}(V_{\beta ,A,\alpha }+\rho _\kappa _{\rho _\beta }V_{\kappa ,A,\alpha })^A\delta _{a,z}]`$
$`{\displaystyle \frac{\delta ^3\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}\delta _{c,z}(_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,\sigma }+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,\sigma }+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma })^A\rho _{\sigma ,z}`$ (D20)
$`+(\delta _{a,z}\delta _{b,z}(_{\rho _\beta }V_{\alpha ,A,\gamma }+_{\rho _\alpha }V_{\beta ,A,\gamma }+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\gamma })^A\delta _{c,z}+perms.)]`$
$`{\displaystyle \frac{\delta ^4\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}\delta \rho _{ϵ,e}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}\delta _{c,z}\delta _{e,z}(_{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,\sigma }+_{\rho _ϵ}_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,\sigma }+_{\rho _ϵ}_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,\sigma }+_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{ϵ,A,\sigma }`$ (D23)
$`+\rho _\kappa _{\rho _ϵ}_{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\sigma })^A\rho _{\sigma ,z}+(\delta _{a,z}\delta _{b,z}\delta _{c,z}(_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,ϵ}+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,ϵ}`$
$`+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,ϵ}+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,ϵ})^A\delta _{e,z}+perms.],`$
where $`|A|=|\nu _0|+\mathrm{}`$. For homogeneous configurations $`\rho _{\alpha ,x}=\rho _\alpha `$ we find

$`{\displaystyle \frac{\delta \gamma [\rho ]}{\delta \rho _{\alpha ,a}}}`$ $`=`$ $`0`$ (D24)
$`{\displaystyle \frac{\delta ^2\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}}}`$ $`=`$ $`{\displaystyle _z}\left[\delta _{a,z}(V_{\alpha ,A,\beta }+\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\beta })^A\delta _{b,z}+\delta _{b,z}(V_{\beta ,A,\alpha }+\rho _\kappa _{\rho _\beta }V_{\kappa ,A,\alpha })^A\delta _{a,z}\right]`$ (D25)
$`{\displaystyle \frac{\delta ^3\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}(_{\rho _\beta }V_{\alpha ,A,\gamma }+_{\rho _\alpha }V_{\beta ,A,\gamma }+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\gamma })^A\delta _{c,z}+perms.]`$ (D26)
$`{\displaystyle \frac{\delta ^4\gamma [\rho ]}{\delta \rho _{\alpha ,a}\delta \rho _{\beta ,b}\delta \rho _{\gamma ,c}\delta \rho _{ϵ,e}}}`$ $`=`$ $`{\displaystyle _z}[\delta _{a,z}\delta _{b,z}\delta _{c,z}(_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,ϵ}+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,ϵ}+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,ϵ}+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,ϵ})^A`$ (D28)
$`+perms.]`$
or in Fourier space

$`{\displaystyle \frac{\delta ^2\gamma [\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}}}`$ $`=`$ $`{\displaystyle _{a,b,z}}e^{iap+ibq}\delta _{a,z}\delta _{b,z}\left[(V_{\alpha ,A,\beta }+\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\beta })(iq)^A+(V_{\beta ,A,\alpha }+\rho _\kappa _{\rho _\beta }V_{\kappa ,A,\alpha })(ip)^A\right]`$ (D29)
$`=`$ $`\delta _{p+q,0}\left[(V_{\alpha ,A,\beta }+\rho _\kappa _{\rho _\alpha }V_{\kappa ,A,\beta })(iq)^A+(V_{\beta ,A,\alpha }+\rho _\kappa _{\rho _\beta }V_{\kappa ,A,\alpha })(ip)^A\right]`$ (D30)
$`{\displaystyle \frac{\delta ^3\gamma [\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}\delta \rho _{\gamma ,r}}}`$ $`=`$ $`\delta _{p+q+r,0}[(_{\rho _\beta }V_{\alpha ,A,\gamma }+_{\rho _\alpha }V_{\beta ,A,\gamma }+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,\gamma })(ir)^A+perms.]`$ (D31)
$`{\displaystyle \frac{\delta ^4\gamma [\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}\delta \rho _{\gamma ,r}\delta \rho _{ϵ,s}}}`$ $`=`$ $`\delta _{p+q+r+s,0}[(_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,A,ϵ}+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,A,ϵ}+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,A,ϵ}+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,A,ϵ})(is)^A`$ (D33)
$`+perms.].`$
The first four functional derivatives of the full effective action on homogeneous background are given by

$`{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}}}`$ $`=`$ $`\delta _{p+q,0}\mathrm{\Gamma }_{p,q}^{(0)\alpha ,\beta }`$ (D34)
$`{\displaystyle \frac{\delta ^3\mathrm{\Gamma }[\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}\delta \rho _{\gamma ,r}}}`$ $`=`$ $`\mathrm{\Gamma }_{p,q,r}^{\alpha ,\beta ,\gamma }\delta _{p+q+r,0},`$ (D35)
$`{\displaystyle \frac{\delta ^4\mathrm{\Gamma }[\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}\delta \rho _{\gamma ,r}\delta \rho _{ϵ,s}}}`$ $`=`$ $`\mathrm{\Gamma }_{p,q,r,s}^{\alpha ,\beta ,\gamma ,ϵ}\delta _{p+q+r+s,0},`$ (D36)
where

$`\mathrm{\Gamma }_{p,q}^{(0)\alpha ,\beta }`$ $`=`$ $`V_{\alpha ,\beta }(q)+\rho _\kappa _{\rho _\alpha }V_{\kappa ,\beta }(q)+V_{\beta ,\alpha }(p)+\rho _\kappa _{\rho _\beta }V_{\kappa ,\alpha }(p)_{\rho _\alpha }_{\rho _\beta }U`$ (D37)
$`=`$ $`V_{\alpha ,\beta }(q)+V_{\beta ,\alpha }(q)+\rho _\kappa _{\rho _\alpha }V_{\kappa ,\beta }(q)+\rho _\kappa _{\rho _\beta }V_{\kappa ,\alpha }(q)_{\rho _\alpha }_{\rho _\beta }U`$ (D38)
$`\mathrm{\Gamma }_{p,q,r}^{\alpha ,\beta ,\gamma }`$ $`=`$ $`(_{\rho _\beta }V_{\alpha ,\gamma }(r)+_{\rho _\alpha }V_{\beta ,\gamma }(r)+\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,\gamma }(r)+perms.)_{\rho _\alpha }_{\rho _\beta }_{\rho _\gamma }U,`$ (D39)
$`=`$ $`\rho _\kappa _{\rho _\beta }_{\rho _\alpha }V_{\kappa ,\gamma }(r)+_{\rho _\beta }V_{\alpha ,\gamma }(r)+_{\rho _\alpha }V_{\beta ,\gamma }(r)+\rho _\kappa _{\rho _\alpha }_{\rho _\gamma }V_{\kappa ,\beta }(q)+_{\rho _\gamma }V_{\alpha ,\beta }(q)+_{\rho _\alpha }V_{\gamma ,\beta }(q)`$ (D41)
$`+\rho _\kappa _{\rho _\beta }_{\rho _\gamma }V_{\kappa ,\alpha }(p)+_{\rho _\gamma }V_{\beta ,\alpha }(p)+_{\rho _\beta }V_{\gamma ,\alpha }(p)_{\rho _\alpha }_{\rho _\beta }_{\rho _\gamma }U,`$
$`\mathrm{\Gamma }_{p,q,r,s}^{\alpha ,\beta ,\gamma ,ϵ}`$ $`=`$ $`_{\rho _\gamma }_{\rho _\beta }V_{\alpha ,ϵ}(s)+_{\rho _\gamma }_{\rho _\alpha }V_{\beta ,ϵ}(s)+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,ϵ}(s)+\rho _\kappa _{\rho _\gamma }_{\rho _\beta }_{\rho _\alpha }V_{\kappa ,ϵ}(s)+perms._{\rho _\alpha }_{\rho _\beta }_{\rho _\gamma }_{\rho _ϵ}U`$ (D42)
$`=`$ $`\rho _\kappa _{\rho _\beta }_{\rho _\alpha }_{\rho _\gamma }V_{\kappa ,ϵ}(s)+_{\rho _\alpha }_{\rho _\gamma }V_{\beta ,ϵ}(s)+_{\rho _\beta }_{\rho _\gamma }V_{\alpha ,ϵ}(s)+_{\rho _\beta }_{\rho _\alpha }V_{\gamma ,ϵ}(s)`$ (D46)
$`+(\rho _\kappa _{\rho _\beta }_{\rho _\alpha }_{\rho _ϵ}V_{\kappa ,\gamma }(r)+_{\rho _\alpha }_{\rho _ϵ}V_{\beta ,\gamma }(r)+_{\rho _\beta }_{\rho _ϵ}V_{\alpha ,\gamma }(r)+_{\rho _\beta }_{\rho _\alpha }V_{ϵ,\gamma }(r)`$
$`+(\rho _\kappa _{\rho _\beta }_{\rho _\gamma }_{\rho _ϵ}V_{\kappa ,\alpha }(p)+_{\rho _\gamma }_{\rho _ϵ}V_{\beta ,\alpha }(p)+_{\rho _\beta }_{\rho _ϵ}V_{\gamma ,\alpha }(p)+_{\rho _\beta }_{\rho _\gamma }V_{ϵ,\alpha }(p)`$
$`+(\rho _\kappa _{\rho _\alpha }_{\rho _\gamma }_{\rho _ϵ}V_{\kappa ,\beta }(q)+_{\rho _\gamma }_{\rho _ϵ}V_{\alpha ,\beta }(q)+_{\rho _\alpha }_{\rho _ϵ}V_{\gamma ,\beta }(q)+_{\rho _\alpha }_{\rho _\gamma }V_{ϵ,\beta }(q)_{\rho _\alpha }_{\rho _\beta }_{\rho _\gamma }_{\rho _ϵ}U.`$
and the dependence on $`\rho `$ is suppresed. Therefore the second functional derivative up to $`𝒪(\delta \rho ^2)`$ is

$`{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[\rho ]}{\delta \rho _{\alpha ,p}\delta \rho _{\beta ,q}}}`$ $`=`$ $`\delta _{p+q,0}\mathrm{\Gamma }_{p,q}^{(0)\alpha ,\beta }+\mathrm{\Gamma }_{p,q}^{(1)\alpha ,\beta }+{\displaystyle \frac{1}{2}}\mathrm{\Gamma }_{p,q}^{(2)\alpha ,\beta }+𝒪(\delta \rho ^3)`$ (D47)
$`\mathrm{\Gamma }_{p,q}^{(1)\alpha ,\beta }`$ $`=`$ $`{\displaystyle _r}\delta _{p+q+r,0}\delta \rho _{\gamma ,r}\mathrm{\Gamma }_{p,q,r}^{\alpha ,\beta ,\gamma },`$ (D48)
$`\mathrm{\Gamma }_{p,q}^{(2)\alpha ,\beta }`$ $`=`$ $`{\displaystyle _{r,s}}\delta _{p+q+r+s,0}\delta \rho _{\gamma ,r}\delta \rho _{\rho ,s}\mathrm{\Gamma }_{p,q,r,s}^{\alpha ,\beta ,\gamma ,\rho }.`$ (D49)
The appplication of these formulae for the effective action (74) yields the second derivative

$$\mathrm{\Gamma }_{p,q}^{(0)\alpha ,\beta }=\left(\begin{array}{cc}\mathrm{\Gamma }^t+\frac{1}{2}_1^0\mathrm{\Gamma }^{tt}(p)+\frac{1}{2}_1^0\mathrm{\Gamma }^{tt}(q)& i𝐩\frac{1}{2}\mathrm{\Gamma }^{st}(p)+i𝐪\frac{1}{2}_1^0\mathrm{\Gamma }^{ts}(q)\\ i𝐪\frac{1}{2}\mathrm{\Gamma }^{st}(q)+i𝐩\frac{1}{2}_1^0\mathrm{\Gamma }^{ts}(p)& \frac{1}{2}\mathrm{\Gamma }^T(q)T_q+\frac{1}{2}\mathrm{\Gamma }^L(q)L_q+\frac{1}{2}\mathrm{\Gamma }^T(p)T_p+\frac{1}{2}\mathrm{\Gamma }^L(p)L_p\end{array}\right)$$
(D50)
where $`\mathrm{\Gamma }^t=_0^2U`$, $`\mathrm{\Gamma }^s=_1^2U`$. The propagator is

$$\delta _{p+q,0}\mathrm{\Gamma }_{p,q}^{(0)1\alpha ,\beta }=\left(\begin{array}{cc}\frac{\mathrm{\Gamma }^L}{\mathrm{\Gamma }_D}& i𝐪\frac{1}{2}\frac{_2^0\mathrm{\Gamma }^{ts}}{\mathrm{\Gamma }_D}\\ i𝐪\frac{1}{2}\frac{_2^0\mathrm{\Gamma }^{ts}}{\mathrm{\Gamma }_D}& \frac{1}{\mathrm{\Gamma }^T}T+\frac{\mathrm{\Gamma }^t+_1^0\mathrm{\Gamma }^{tt}}{\mathrm{\Gamma }_D}L\end{array}\right)$$
(D51)
with

$$\mathrm{\Gamma }_D=(\mathrm{\Gamma }^t+_1^0\mathrm{\Gamma }^{tt})\mathrm{\Gamma }^L+\frac{𝐪^2}{4}(_2^0\mathrm{\Gamma }^{ts})^2.$$
(D52)
The third and fourth derivatives are of the form

$`\mathrm{\Gamma }_{p,q,r}^{0,0,0}`$ $`=`$ $`_2^1V_{0,0}(p)+_2^1V_{0,0}(q)+_2^1V_{0,0}(r)+_0\mathrm{\Gamma }^t`$ (D53)
$`=`$ $`{\displaystyle \frac{1}{2}}_2^1\mathrm{\Gamma }^{tt}(p)+{\displaystyle \frac{1}{2}}_2^1\mathrm{\Gamma }^{tt}(q)+{\displaystyle \frac{1}{2}}_2^1\mathrm{\Gamma }^{tt}(r)+_0\mathrm{\Gamma }^t`$ (D54)
$`\mathrm{\Gamma }_{p,q,r}^{0,0,\mathrm{}}`$ $`=`$ $`_2^1V_{0,\mathrm{}}(r)+_0V_{\mathrm{},0}(p)+_0V_{\mathrm{},0}(q)+_{\mathrm{}}\mathrm{\Gamma }^t`$ (D55)
$`=`$ $`{\displaystyle \frac{i}{2}}_2^1\mathrm{\Gamma }^{ts}(r)r_{\mathrm{}}+{\displaystyle \frac{i}{2}}_0\mathrm{\Gamma }^{st}(p)p_{\mathrm{}}+{\displaystyle \frac{i}{2}}_0\mathrm{\Gamma }^{st}(q)q_{\mathrm{}}+_{\mathrm{}}\mathrm{\Gamma }^t`$ (D56)
$`\mathrm{\Gamma }_{p,q,r}^{0,k,\mathrm{}}`$ $`=`$ $`_0V_{k,\mathrm{}}(r)+_0V_{\mathrm{},k}(q)+\delta _{k,\mathrm{}}_0\mathrm{\Gamma }^s`$ (D57)
$`=`$ $`{\displaystyle \frac{1}{2}}_0\mathrm{\Gamma }^T(r)T_r^{k,\mathrm{}}+{\displaystyle \frac{1}{2}}_0\mathrm{\Gamma }^T(q)T_q^{\mathrm{},k}+{\displaystyle \frac{1}{2}}_0\mathrm{\Gamma }^L(r)L_r^{k,\mathrm{}}+{\displaystyle \frac{1}{2}}_0\mathrm{\Gamma }^L(q)L_q^{\mathrm{},k}+\delta _{k,\mathrm{}}_0\mathrm{\Gamma }^s`$ (D58)
$`\mathrm{\Gamma }_{p,q,r,s}^{0,0,0,0}`$ $`=`$ $`_3^2V_{0,0}(p)+_3^2V_{0,0}(q)+_3^2V_{0,0}(r)+_3^2V_{0,0}(s)+_0^2\mathrm{\Gamma }^t`$ (D59)
$`=`$ $`{\displaystyle \frac{1}{2}}_3^2\mathrm{\Gamma }^{tt}(p)+{\displaystyle \frac{1}{2}}_3^2\mathrm{\Gamma }^{tt}(q)+{\displaystyle \frac{1}{2}}_3^2\mathrm{\Gamma }^{tt}(r)+{\displaystyle \frac{1}{2}}_3^2\mathrm{\Gamma }^{tt}(s)+_0^2\mathrm{\Gamma }^t`$ (D60)
$`\mathrm{\Gamma }_{p,q,r,s}^{0,0,0,k}`$ $`=`$ $`_3^2V_{0,k}(s)+_0^2V_{k,0}(p)+_0^2V_{k,0}(q)+_0^2V_{k,0}(r)+_0_k\mathrm{\Gamma }^t`$ (D61)
$`=`$ $`{\displaystyle \frac{i}{2}}_3^2\mathrm{\Gamma }^{ts}(s)s^k+{\displaystyle \frac{i}{2}}_0^2\mathrm{\Gamma }^{st}(p)p^k+{\displaystyle \frac{i}{2}}_0^2\mathrm{\Gamma }^{st}(q)q^k+{\displaystyle \frac{i}{2}}_0^2\mathrm{\Gamma }^{st}(r)r^k+_0_k\mathrm{\Gamma }^t`$ (D62)
$`\mathrm{\Gamma }_{p,q,r,s}^{0,0,k,\mathrm{}}`$ $`=`$ $`_0^2V_{k,\mathrm{}}(s)+_0^2V_{\mathrm{},k}(r)+\delta ^{k,\mathrm{}}_0^2\mathrm{\Gamma }^s`$ (D63)
$`=`$ $`{\displaystyle \frac{1}{2}}_0^2\mathrm{\Gamma }^T(s)T_s^{k,\mathrm{}}+{\displaystyle \frac{1}{2}}_0^2\mathrm{\Gamma }^L(s)L_s^{k,\mathrm{}}+{\displaystyle \frac{1}{2}}_0^2\mathrm{\Gamma }^T(r)T_r^{k,\mathrm{}}+{\displaystyle \frac{1}{2}}_0^2\mathrm{\Gamma }^L(r)L_r^{k,\mathrm{}}+\delta ^{k,\mathrm{}}_0^2\mathrm{\Gamma }^s.`$ (D64)
## E One-loop $`W[\sigma ]`$

This Appendix we briefly record the expressions needed for the computation of the Green functions occuring in the one-loop approximation of the generator funtional $`W[\sigma ]`$. This generator functional is non-trivial even for free electrons due to the exchange interaction.

### 1 Densities

We start with the expresions for the densities

$$\rho ^{}=\frac{i\mathrm{}}{2}G_{t,𝐱,t+\delta ,𝐱}+c.c.=_𝐪\mathrm{\Theta }(E_𝐪)$$
(E1)
and

$$\stackrel{~}{\rho }^{}=\frac{1}{3}\underset{j}{}\stackrel{~}{\rho }_j^{}=\frac{i\mathrm{}}{6}\underset{j>0}{}G_{t,𝐱,t+\delta ,𝐱+j}+c.c.=_𝐪\mathrm{\Theta }(E_𝐪)\mathrm{cos}q_1a$$
(E2)
in terms of the electron propagator $`G_p`$. The local potential is given by

$$\frac{\mathrm{}}{Vi}\mathrm{T}r\mathrm{log}D^1=i_{𝐪,E}\mathrm{log}[EE_𝐪+iϵ_𝐪].$$
(E3)
Charge conjugation acts as $`\mu \mu ^c=6\mathrm{}^2/ma^2\mu `$, $`\rho ^{}\rho ^c=_\rho \rho ^{}`$ and $`\stackrel{~}{\rho }^{}\stackrel{~}{\rho }^c=\stackrel{~}{\rho }^{}`$ where the band contribution is

$$_\rho =_𝐪1$$
(E4)
### 2 Two-point functions

We need

$`\stackrel{~}{D}_p^{0,0}`$ $`=`$ $`i\mathrm{}{\displaystyle _r}G_rG_{p+r},`$ (E5)
$`\stackrel{~}{D}_p^{0,j}`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}^2}{ma}}e^{i\frac{p_ja}{2}}{\displaystyle _r}G_rG_{p+r}\mathrm{sin}a\left(r_j+{\displaystyle \frac{p_j}{2}}\right)`$ (E6)
$`\stackrel{~}{D}_p^{j,k}`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}^3}{m^2a^2}}{\displaystyle _r}G_rG_{p+r}\mathrm{sin}a\left(r_j+{\displaystyle \frac{p_j}{2}}\right)\mathrm{sin}a\left(r_k+{\displaystyle \frac{p_k}{2}}\right),`$ (E7)
where the Fourier transform is defined as

$$\stackrel{~}{D}_{x,y}^{\alpha ,\beta }=_p\stackrel{~}{D}_p^{\alpha ,\beta }e^{ip(xy)}.$$
(E8)
It is straightforward to perform the energy integrals. The results are shown for the external momentum $`\widehat{𝐩}=(0,0,\widehat{p})`$ after having carried out the replacement $`\mu \mu \sigma _0`$,

$`\stackrel{~}{D}^{0,0}(\omega ,𝐩^2)`$ $`=`$ $`2𝒫{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2\mathrm{}^2\omega ^2}}+i\pi {\displaystyle _𝐪}\mathrm{\Theta }(E_𝐪)\mathrm{\Theta }(E_{𝐪+𝐩})[\delta (\mathrm{\Delta }E_{𝐪,𝐩}\mathrm{}\omega )+\delta (\mathrm{\Delta }E_{𝐪,𝐩}+\mathrm{}\omega )]`$ (E9)
$`\stackrel{~}{D}^{0,j}(\omega ,𝐩)`$ $`=`$ $`\delta ^{j,3}{\displaystyle \frac{2\omega \mathrm{}}{am}}𝒫{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}a\left(q_3+\frac{p}{2}\right)}{(\mathrm{\Delta }E_{𝐪,𝐩})^2\mathrm{}^2\omega ^2}}`$ (E11)
$`{\displaystyle \frac{i\mathrm{}\pi }{am}}{\displaystyle _𝐪}\mathrm{\Theta }(E_𝐪)\mathrm{\Theta }(E_{𝐪+𝐩})\mathrm{sin}a\left(q_j+{\displaystyle \frac{p_j}{2}}\right)[\delta (\mathrm{\Delta }E_{𝐪,𝐩}\mathrm{}\omega )\delta (\mathrm{\Delta }E_{𝐪,𝐩}+\mathrm{}\omega )]`$
$`=`$ $`\stackrel{~}{D}^{0,j}(\omega ,𝐩)=\stackrel{~}{D}^{0,j}(\omega ,𝐩)=\stackrel{~}{D}^{j,0}(\omega ,𝐩)`$ (E12)
$`\stackrel{~}{D}^T(\omega ,𝐩^2)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{a^2m^2}}𝒫{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}^2aq_1\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2\mathrm{}^2\omega ^2}}`$ (E14)
$`+{\displaystyle \frac{i\pi \mathrm{}^2}{a^2m^2}}{\displaystyle _𝐪}\mathrm{\Theta }(E_𝐪)\mathrm{\Theta }(E_{𝐪+𝐩})\mathrm{sin}^2aq_1[\delta (\mathrm{\Delta }E_{𝐪,𝐩}\mathrm{}\omega )+\delta (\mathrm{\Delta }E_{𝐪,𝐩}+\mathrm{}\omega )]`$
$`\stackrel{~}{D}^L(\omega ,𝐩^2)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{a^2m^2}}𝒫{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}^2a\left(q_3+\frac{p}{2}\right)\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2\mathrm{}^2\omega ^2}}`$ (E16)
$`+{\displaystyle \frac{i\pi \mathrm{}^2}{a^2m^2}}{\displaystyle _𝐪}\mathrm{\Theta }(E_𝐪)\mathrm{\Theta }(E_{𝐪+𝐩})\mathrm{sin}^2a\left(q_1+{\displaystyle \frac{p}{2}}\right)[\delta (\mathrm{\Delta }E_{𝐪,𝐩}\mathrm{}\omega )+\delta (\mathrm{\Delta }E_{𝐪,𝐩}+\mathrm{}\omega )],`$
where $`\mathrm{\Delta }E_{𝐪,𝐩}=E_𝐪E_{𝐪,𝐩}`$ and $`𝒫`$ denotes the principal value prescription. The singularity and the Dirac delta in the integrand render these expressions unpractical for numerical evaluation. We used Wick rotated integrals and considered the dependence in Euclidean energies in the numerical work. The corresponding expressions are

$`\stackrel{~}{D}^{0,0}(\omega ,𝐩^2)`$ $`=`$ $`2{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2+\mathrm{}^2\omega ^2}}`$ (E17)
$`\stackrel{~}{D}^{0,j}(\omega ,𝐩)`$ $`=`$ $`{\displaystyle \frac{2i\omega \mathrm{}^2}{am}}{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}a\left(q_3+\frac{p}{2}\right)1}{(\mathrm{\Delta }E_{𝐪,𝐩})^2+\mathrm{}^2\omega ^2}}`$ (E18)
$`\stackrel{~}{D}^T(\omega ,𝐩^2)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{a^2m^2}}{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}^2aq_1\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2+\mathrm{}^2\omega ^2}}`$ (E19)
$`\stackrel{~}{D}^L(\omega ,𝐩^2)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{}^2}{a^2m^2}}{\displaystyle _𝐪}{\displaystyle \frac{\mathrm{\Theta }(E_𝐪)\mathrm{sin}^2a\left(q_3+\frac{p}{2}\right)\mathrm{\Delta }E_{𝐪,𝐩}}{(\mathrm{\Delta }E_{𝐪,𝐩})^2+\mathrm{}^2\omega ^2}}.`$ (E20)
These expressions show clearly the need of keeping a non-vanishing energy $`\mathrm{}\omega `$ in any computation carried out in a finite system. In fact, the limit $`\omega 0`$ generates a Dirac delta in the integrands which yields singular results in a finite system where the density of states is the sum of Dirac delta peaks.