%\documentclass[12pt,aps,floats,superscriptaddress]{revtex4}
\documentclass[12pt,aps]{revtex4}
\usepackage{amsmath}

\begin{document}

\begin{center}
{\Large Solution for Homework 2}
\end{center}

\vspace{.5cm}
\noindent
Problem 1.

We consider the symmetric ($\epsilon_f=-{\frac{U}{2}}$) 
Anderson impurity model:
$$
H=\sum_{k\sigma}\epsilon_k c_{k\sigma}^{\dagger}c_{k\sigma}
+\sum_{k\sigma} (V_k c_{k\sigma}^{\dagger}f_{\sigma}+
V_k^*f_{\sigma}^{\dagger}c_{k\sigma})
+\epsilon_f(n_{\uparrow}+n_{\downarrow})+Un_{\uparrow}n_{\downarrow}
$$
Evaluating the impurity Greens function 
$G_{\sigma}(\tau)=
\langle f_{\sigma}(\tau) f_{sigma}^{\dagger}(0)\rangle$
we can integrate out the conduction electrons $c_{k\sigma}$ to
end up with:
$$
G(i\omega)=\frac{1}{i\omega-\epsilon_f-\Delta(i\omega)-\Sigma(i\omega)}
$$
where the the hybridization function is
$$
\Delta(i\omega_n)=
\sum_k{|V_k|^2\over i\omega_n-\epsilon_k}=-i \Gamma{\rm sign}(\omega_n)
$$
and $\Gamma=V^2\rho=1$. In the symmetric case the Hartrey term in
the self energy cancels $\epsilon_f$ and we only need to consider 
the second order term, so:
\begin{equation}
G(i\omega)=\frac{1}{G_0^{-1}(i\omega)-\Sigma_2(i\omega)}
\label{greensf}
\end{equation}
where $G_0^{-1}(i\omega)=i\omega+i\Gamma{\rm sign}(\omega_n)$.
In bare perturbation theory the self energy can be written 
as a convolution type integral in frequency representation:
\begin{equation}
\Sigma_2(i\omega)=U^2 T^2\sum_{i\omega',i\Omega}
G_0(i\omega+i\Omega)G_0(i\omega')G_0(i\omega'+i\Omega)
\label{sef}
\end{equation}
or as a product of Green's function in time representation:
\begin{equation}
\Sigma_2(\tau)=U^2 G_0^3(\tau)
\label{set}
\end{equation}
In selfconsistent perturbation theory the free Green's function
$G_0$ in the above two equations should be substituted by 
the full Green's function defined in Eq.(\ref{greensf}) and
equations solved selfconsistently.

When solving equations numerically we discretize variables
and store functions as arrays (lists). The analog of the
Fourier transformation of the continues function $a(\tau)$
is the discrete Fourier transformation of the discrete
function $a_r$ will be:
\begin{eqnarray}
b_s&=&\frac{c}{N}\sum_{r=-n}^{r=n}a_r {\rm e}^{2\pi i r s/N}
\nonumber\\
a_r&=&c^{-1}\sum_{s=-n}^{s=n}b_s {\rm e}^{-2\pi i s r/N}
\label{dft}
\end{eqnarray}
where the second line is the inverse transformation, 
$c$ is an arbitrary constant and $N=2n+1$. To make
correspondence to Matsubara frequencies we set:
$$
\omega_s=\pi T s,\,\,\,\,\, \tau_r=2\beta r/N,\,\,\,\,\, c=\beta
$$
In the limit $N\to\infty$ Eq.(\ref{dft}) becomes:
\begin{eqnarray}
a(i\omega)&=&\frac{1}{2}\int_{-\beta}^{\beta}d\tau a(\tau) 
{\rm e}^{i\omega_n\tau}
\nonumber\\
a(\tau)&=&T\sum_{i\omega_n} a(i\omega){\rm e}^{-i\omega_n\tau}
\label{matsubara}
\end{eqnarray}
Frequencies $s=\pm 1,\, \pm 3,\,...$ are relevant in Fermi statistics,
while $s=0,\, \pm 2,\, \pm 4,\,...$ are relevant in Bose statistics.

Two-point functions in perturbation theory have (slowly decaying)
asymptotics: $C(i\omega)~1/(i\omega)$. On the other hand, 
number of matsubara frequencies which we can sum up
using a computer is quite finite. To handle this problem we
use the knowledge of the function's asymptotics to separate
it into a slowly decaying part, which can be summed analytically,
and a fast decaying part, which can be easily summed numerically.
For example, to Fourier transform a Green's function we do:
$$
G(i\omega)=\{G(i\omega)-(i\omega)^{-1}\}+(i\omega)^{-1}
$$
The first term (curly brackets) has asymptotics $~\omega^{-2}$ and
can be Fourier transformed numerically. For the last term we have:
$$
T\sum_{i\omega_n} \frac{1}{i\omega}{\rm e}^{-i\omega_n\tau}=
-\frac{1}{2}\Theta(\tau)
$$
In case of the real frequencies, similarly:
$$
G(i\omega)=\{G(\omega)-\omega^{-1}\}+\omega^{-1}
$$
$$
\mathcal{P}\int \frac{d\omega}{2\pi}\omega^{-1}{\rm e}^{-i\omega\tau}=
-\frac{i}{2}\Theta(\tau)
$$

In general, a convolution type transformation, like Fourier transformation,
if implemented effectively, can be performed in $O(N \ln{N})$ time, much
faster than a naive $O(N^2)$ estimate. If one computes self energy
straitforwardly in frequency domain, like prescribed by Eq.(\ref{sef}),
one would need $O(N^2)$ time. Therefor we chose to resort to time domain 
Eq.(\ref{set}), price for that is Fourier transformations, which can be 
done in $O(N \ln{N})$. Here is the algorithm:
$$
G(i\omega)\xrightarrow{\rm InvFT}G(\tau)\xrightarrow{}\Sigma(\tau)
\xrightarrow{\rm FT}\Sigma(i\omega)\xrightarrow{}G(i\omega)
$$
In the bare perturbation theory we start with $G=G_0$ and iterate
this loop once, while in selfconsistent approach the loop is iterated
until convergence is reached.

To compute the spectral function $A(\omega)$ at $T=0$ 
we compute Green's function in real frequency $G(\omega)$ 
and take the imaginary part:
\begin{equation}
A(\omega)=-{\rm Im}G_R(\omega)=-{\rm Im}G(\omega+i\delta)
\label{spectralf}
\end{equation}
The computation of finite $T$ spectral function is more involved.
We first derive the expression of the self energy in terms of
the spectral functions. We use the following representation
of the Green's function:
$$
G(\omega)=-\frac{1}{\pi}\int d\omega'\frac{A(\omega')}{\omega'-\omega}
$$
%
%where the spectral function defined as
%\begin{equation}
%A(\omega)=-{\rm Im}G(\omega+i\delta)
%\label{spectralf}
%\end{equation}
to write:
$$
\Sigma_2(i\omega)=-U^2 T^2\sum_{i\omega',i\Omega}
\int \frac{dx dy dz}{\pi^3} 
\frac{1}{(x-i\omega')(y-i\omega'-i\Omega)(z-i\omega-i\Omega)}
$$
Performing standard Matsubara frequency summation, continuing
to real frequencies $i\omega\to\omega+i\delta$ and taking
the imaginary part of the expression we find:
\begin{equation}
{\rm Im}\Sigma(\omega+i\delta)=
-U^2\int\frac{dx dy}{\pi^2}[b(y)+f(y+\omega)][f(x)-f(x+y)]
A(x)A(x+y)A(y+\omega)
\label{sesf}
\end{equation}
where $f(x)$ and $b(x)$ are Fermi and Bose functions.
In principle, we could use the Kramers-Kronig relation 
to compute $\Sigma(\omega)$ from its imaginary part
and then use Eq.(\ref{spectralf}) to close the loop.
However direct use of Eq.(\ref{sesf}) or Kramers Kronig
requires $O(N^2)$ of computational time. We will
show how one can do all calculations in $O(N\ln{N})$.

Using relation $b(y)[f(x)-f(x+y)]=f(-x) f(x+y)$ we rewrite
Eq.(\ref{sesf}) as:
\begin{eqnarray}
&&{\rm Im}\Sigma(\omega+i\delta)=-U^2\int\frac{dx dy}{\pi^2}
\nonumber\\
&&[A_2(x)A_1(x+y)A(y+\omega)+A_1(x)A(x+y)A_1(y+\omega)-
A(x)A_1(x+y)A_1(y+\omega)]
\label{seconv}
\end{eqnarray}
where 
$$
A_1(x)=A(x)f(x),\,\,\,\,\, A_2(x)=A(x)f(-x)
$$
In time representation Eq.(\ref{seconv}) becomes
$$
{\rm Im}\Sigma(t)=-4 U^2 [A_2(t)A_1(-t)A(t)+A_1(t)A(-t)A_1(t)-
A(t)A_1(-t)A_1(t)]
$$
We should how to compute ${\rm Im}\Sigma$ by using Fourier
transform and regular multiplication. Instead of using Kramers
Kronig we can compute real part of self energy in the
same way as we did for real part. For that we use a slightly
different representation for one out of three Green's functions
in Eq.(\ref{sef}):
$$
G(\omega)=-\frac{i}{\pi}\int d\omega'\frac{B(\omega')}{\omega'-\omega}
$$
where $B(\omega)$ is defined as
$$
B(\omega)={\rm Re}G(\omega+i\delta)
$$
Then for the real part of $\Sigma$ we find:
\begin{eqnarray}
&&{\rm Re}\Sigma(\omega)=U^2\int\frac{dx dy}{\pi^2}
\nonumber\\
&&[A_2(x)A_1(x+y)B(y+\omega)-A_1(x)B(x+y)A_1(y+\omega)-
B(x)A_1(x+y)A_1(y+\omega)]
\label{seconvreal}
\end{eqnarray} 
or in time representation:
$$
{\rm Re}\Sigma(t)=4 U^2 [A_2(t)A_1(-t)B(t)+A_1(t)B(-t)A_1(t)-
B(t)A_1(-t)A_1(t)]
$$
This completes calculation of $\Sigma={\rm Re}\Sigma+i{\rm Im}\Sigma$.

In the attached Mathematica notebook file we implement algorithm
we outlined above. All computations are done in $O(N\ln{N})$ time.


\vspace{.5cm}

\end{document}
