View source for PHY6937

'''Welcome to Phy 6937 Superconductivity and superfluidity'''

PHY6937 is a one semester advanced graduate level course. Its aim is to introduce concepts and theoretical techniques for the description of superconductors and superfluids. This course is a natural continuation of the "many-body" course PHY5670 and will build on the logical framework introduced therein, i.e. broken symmetry and adiabatic continuity. The course will cover a range of topics, such as the connection between the phenomenological Ginzburg-Landau and the microscpic BCS theory, Migdal-Eliashberg treatment of phonon mediated superconductivity, unconventional superconductivity, superfluidity in He-4 and He-3, and Kosterlitz-Thouless theory of two dimensional superfluids. 


The key component of the course is the collaborative student contribution to the course Wiki-textbook. Each team of students is responsible for BOTH writing the assigned chapter AND editing chapters of others.

'''Team assignments:''' [[Phy6937_Spring11_teams|Spring 2011 student teams]]


----
'''Outline of the course:'''
== Pairing Hamiltonian and BCS instability ==

To see the origins of superconductivity, it is helpful to look at a toy system, which we already know will give us superconducting behavior. This is useful because the toy system is only a simple change to a non-interacting electron gas. By adding in some small ''attractive'' interaction, we will arrive at a superconducting system! This interaction need only occur between two electrons occupying the same position in space (and necessarily having opposite spin!). Additionally, we still find the interesting behaviour regardless of the size of the interaction; the only requirement is that it be non-zero!

We can write the Hamiltonian of the system as:

<math>H=\sum_\vec{r}[\psi_\sigma^\dagger (\vec{r})(\epsilon_\vec{p}-\mu)\psi_\sigma^\dagger (\vec{r}) +g\psi_\uparrow^\dagger (\vec{r})\psi_\downarrow^\dagger (\vec{r})\psi_\downarrow (\vec{r})\psi_\uparrow (\vec{r})]</math>

where <math>\ g<0</math> and <math>\ |g|<<\epsilon_{F}</math>.

For this system, the partition function is:

<math>Z=\int D[\psi_\sigma ^{*} (\tau, \vec{r}), \psi_\sigma  (\tau, \vec{r})]e^{-S_{BCS}}</math>

where <math>S_{BCS}=\int_0^\beta d\tau \sum_\vec{r}[\psi_\sigma^\dagger (\tau, \vec{r})(\partial _\tau+\epsilon_\vec{p}-\mu)\psi_\sigma^\dagger (\vec{r}) +g\psi_\uparrow^\dagger (\vec{r})\psi_\downarrow^\dagger (\vec{r})\psi_\downarrow (\vec{r})\psi_\uparrow (\vec{r})]</math>

It doesn't matter to multiply partition function by a constant:

<math>Z\rightarrow Z=\int D[\psi_\sigma ^{*} (\tau, \vec{r}), \psi_\sigma  (\tau, \vec{r})] D[\Delta^{*}(\tau, \vec{r}),\Delta (\tau, \vec{r})] e^{-S_{BCS}-S_{\Delta}}</math>

where <math>S_\Delta=-\int_0^\beta d\tau\sum_{\vec{r}}\frac{1}{g}\Delta^*(\tau,\vec{r})\Delta(\tau,\vec{r})</math>

<math>\psi^\dagger</math> and <math>\ \psi</math> are grassmann numbers.
<math>\ \Delta^*</math> and <math>\ \Delta</math> are constant.
<math>\psi_\uparrow\psi_\downarrow</math> and <math>\psi_\downarrow\psi_\uparrow</math> behave like constant.

Let's make a shift of the constant:

<math>\Delta \rightarrow \Delta+g\psi_\uparrow\psi_\downarrow</math>

<math>\Delta^*\rightarrow \Delta^*+g\psi^\dagger_\downarrow\psi^\dagger_\uparrow</math>

Then, <math>S_\Delta=-\int_0^\beta d\tau \sum_{\vec{r}}{\{\frac{1}{g}\Delta^*\Delta + \Delta^*\psi_\uparrow \psi_\downarrow + \Delta\psi^\dagger_\downarrow \psi^\dagger_\uparrow+g\psi^\dagger_\downarrow \psi^\dagger_\uparrow \psi_\uparrow \psi_\downarrow}\}</math>

<math>\begin{align}S=&S_{BCS}+S_{\Delta}\\
=&\int_0^\beta d\tau \sum_{\vec{r}}\{  \psi_\sigma^\dagger(\tau, \vec{r})(\partial _\tau+\epsilon_\vec{p}-\mu)\psi_\sigma^\dagger (\tau, \vec{r})  \ \ \ \ \ \ \ \ \ \ \rightarrow  S_0 \\
&+\Delta^*(\tau, \vec{r})\psi_\uparrow (\tau, \vec{r})\psi_\downarrow (\tau, \vec{r}) \Delta (\tau, \vec{r})\psi^\dagger_\downarrow (\tau, \vec{r})\psi^\dagger_\uparrow (\tau, \vec{r})  \rightarrow S_{int}\\
&-\frac{1}{g}\Delta^* (\tau, \vec{r})\Delta (\tau, \vec{r}) \}   \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \   \rightarrow S_{\Delta}
\end{align}</math>

then, <math>Z=\int D[\psi_{\sigma}^{*}(\tau,\mathbf{r}),\psi_{\sigma}(\tau,\mathbf{r})]D[\Delta^{*}(\tau,\mathbf{r}),\Delta(\tau,\mathbf{r})]e^{-(S_{0}+S_{int.}+S_{\Delta})}</math>.

<math>\left\langle e^{-S_{int.}}\right\rangle _{0}\cong exp[\frac{1}{2}\left\langle S_{int.}^{2}\right\rangle _{0}+\frac{1}{4!}(\left\langle S_{int.}^{4}\right\rangle _{0}-3\left\langle S_{int.}^{2}\right\rangle _{0}^{2})]</math> by cumulant expansion, which guarantees that until the 2nd order, it is accurate.

Use Matsubara's Method

<math>\psi_{\sigma}(\tau,\mathbf{r})=\frac{1}{\beta}\underset{\omega_{n}}{\sum}\underset{\mathbf{k}}{\sum}e^{i\mathbf{k}\cdot\mathbf{r}}e^{-i\omega_{n}\tau}\psi_{\sigma}(i\omega_{n},\mathbf{k}), \omega_{n}=(2n+1)\frac{\pi}{\beta};</math>

<math>\Delta(\tau,\mathbf{r})=\frac{1}{\beta}\underset{\Omega_{n}}{\sum}\underset{\mathbf{k}}{\sum}e^{i\mathbf{k}\cdot\mathbf{r}}e^{-i\Omega_{n}\tau}\Delta_{\mathbf{k}}(i\Omega_{n}), \omega_{n}=2n\frac{\pi}{\beta}.</math>

Then <math>S_{0}=\frac{L^{D}}{\beta}\underset{\omega_{n}}{\sum}\underset{\mathbf{k}}{\sum}[-i\omega_{n}+\varepsilon_{\mathbf{k}}-\mu]\psi_{\sigma}^{\dagger}(i\omega_{n},\mathbf{k})\psi_{\sigma}(i\omega_{n},\mathbf{k}).</math>

<math>S_{int.}=\frac{L^{D}}{\beta^{2}}\underset{\omega_{n},\Omega_{n}}{\sum}\underset{\mathbf{k},\mathbf{q}}{\sum}[\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\psi_{\uparrow}(i\Omega_{n}-i\omega_{n},\mathbf{\mathbf{q}-k})\psi_{\downarrow}(i\omega_{n},\mathbf{k})+\Delta_{\mathbf{q}}(i\Omega_{n})\psi_{\downarrow}^{\dagger}(i\omega_{n},\mathbf{k})\psi_{\uparrow}^{\dagger}(i\Omega_{n}-i\omega_{n},\mathbf{\mathbf{q}-k})].</math>

The Fourier transform of 1 body Green's function is (<math>i=1,2</math> mean <math>{\mathbf{r}_{i},\tau_{i}}</math>) <math>G(1-2)=\left\langle \psi(1)\psi^{*}(2)\right\rangle =\frac{1}{\beta}\underset{\omega_{n}}{\sum}\frac{1}{L^{D}}\underset{\mathbf{k}}{\sum}e^{-i\omega_{n}(\tau_{1}-\tau_{2})}e^{i\mathbf{k}\cdot(\mathbf{r}_{1}-\mathbf{r}_{2})}\frac{1}{-i\omega_{n}+\varepsilon_{\mathbf{k}}-\mu}</math>,

so <math>G_{\sigma}^{0}(i\omega_{n},\mathbf{k})=\left\langle \psi_{\sigma}(i\omega_{n},\mathbf{k})\psi_{\sigma}^{\dagger}(i\omega_{n},\mathbf{k})\right\rangle _{0}=\frac{\beta}{L^{D}}\frac{1}{-i\omega_{n}+\varepsilon_{\mathbf{k}}-\mu}</math>.

Then <math>\left\langle S_{int.}^{2}\right\rangle _{0}=\frac{2L^{2D}}{\beta^{4}}\underset{\omega_{n},\Omega_{n}}{\sum}\underset{\mathbf{k},\mathbf{q}}{\sum}[G_{\uparrow}^{0}(i\omega_{n},\mathbf{k})G_{\downarrow}^{0}(i\Omega_{n}-i\omega_{n},\mathbf{q}-\mathbf{k})]\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})=L^{D}\frac{2}{\beta}\underset{\Omega_{n},\mathbf{q}}{\sum}\chi_{p}(\mathbf{q},i\Omega_{n})\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})</math>,

where <math>\chi_{p}(\mathbf{q},i\Omega_{n})=\frac{L^{D}}{\beta^{3}}\underset{\omega_{n},\mathbf{k}}{\sum}G_{\uparrow}^{0}(i\omega_{n},\mathbf{k})G_{\downarrow}^{0}(i\Omega_{n}-i\omega_{n},\mathbf{q}-\mathbf{k})</math> is called pairing susceptibility.

Let's calculate it:

<math>\chi_{p}(\mathbf{q},i\Omega_{n})=\frac{L^{D}}{\beta^{3}}\underset{\omega_{n},\mathbf{k}}{\sum}G_{\uparrow}^{0}(i\omega_{n},\mathbf{k})G_{\downarrow}^{0}(i\Omega_{n}-i\omega_{n},\mathbf{q}-\mathbf{k})=\frac{1}{L^{D}}\frac{1}{\beta}\underset{\omega_{n},\mathbf{k}}{\sum}\frac{-1}{i\omega_{n}-\varepsilon_{\mathbf{k}}+\mu}\times\frac{1}{i\omega_{n}-i\Omega_{n}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\mu}</math>,

[[Image:Fig.1.jpg]]

<math>\Rightarrow=\frac{1}{L^{D}}\frac{1}{\beta}\underset{\mathbf{k}}{\sum}\oint_{c}\frac{dz}{2\pi i}\frac{-1}{z-\varepsilon_{\mathbf{k}}+\mu}\times\frac{1}{z-i\Omega_{n}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\mu}\frac{1}{e^{\beta z}+1}</math>.

Since <math>\frac{-1}{z-\varepsilon_{\mathbf{k}}+\mu}\times\frac{1}{z-i\Omega_{n}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\mu}=\frac{1}{\varepsilon_{\mathbf{q}-\mathbf{k}}+\varepsilon_{\mathbf{k}}-2\mu-i\Omega_{n}}[\frac{1}{z-\varepsilon_{\mathbf{q}}+\mu}-\frac{1}{z-i\Omega_{n}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\mu}]</math>,

and change the integral path to 

[[Image:Fig.2.jpg]]

<math>\Rightarrow=-\frac{1}{L^{D}}\frac{1}{\beta}\underset{\mathbf{k}}{\sum}\frac{1}{\varepsilon_{\mathbf{q}-\mathbf{k}}+\varepsilon_{\mathbf{k}}-2\mu-i\Omega_{n}}[\frac{1}{e^{\beta(\varepsilon_{\mathbf{q}}-\mu)}+1}-\frac{1}{e^{\beta(-\varepsilon_{\mathbf{q}-\mathbf{k}}+\mu)}+1}]=\int\frac{d^{D}k}{(2\pi)^{D}}\frac{1}{\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{q}-\mathbf{k}}-2\mu-i\Omega_{n}}[1-f(\varepsilon_{\mathbf{k}})-f(\varepsilon_{\mathbf{q}-\mathbf{k}})].</math>

In the static (<math>\ \Omega_{n}=0</math>) and uniform (<math>\mathbf{q}=0</math>) limit,<math> 1-2f(\varepsilon_{\mathbf{k}})=Tanh[\frac{\beta}{2}(\varepsilon_{\mathbf{k}}-\mu)]</math>.

Then <math>\chi_{p}(0,0)=\int\frac{d^{D}k}{(2\pi)^{D}}\frac{Tanh[\frac{\beta}{2}(\varepsilon_{\mathbf{k}}-\mu)]}{2(\varepsilon_{\mathbf{k}}-\mu)}</math>.

In low energy, integrate the energy in the shell near Fermi energy:

[[Image:Fig.3.jpg]]

<math>\Rightarrow\chi_{p}(0,0)\cong N(0)\int_{\hbar\omega_{D}}^{-\hbar\omega_{D}}d\xi\frac{Tanh[\xi\beta/2]}{2\xi}\cong N(0)\int_{0}^{-\hbar\omega_{D}}d\xi\frac{Tanh[\xi\beta/2]}{\xi}=N(0)ln[\frac{2\hbar\omega_{D}e^{\gamma}}{\pi k_{B}T}].</math>

Then <math>\frac{1}{2}\left\langle S_{int.}^{2}\right\rangle _{0}=L^{D}\frac{1}{\beta}\chi_{p}(0,0)\underset{\Omega_{n},\mathbf{q}}{\sum}\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})</math>.

If we ignore the higher order in the cumulant expansion,

<math>S_{eff}=-\underset{\mathbf{r}}{\sum}\int_{0}^{\beta}d\tau\frac{1}{g}\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})-\frac{1}{2}\left\langle S_{int.}^{2}\right\rangle _{0}=\underset{\mathbf{r}}{\sum}\int_{0}^{\beta}d\tau[\frac{1}{\left|g\right|}-N(0)ln(\frac{2\hbar\omega_{D}e^{\gamma}}{\pi k_{B}T})]\Delta^{*}(\tau,\mathbf{r})\Delta(\tau,\mathbf{r})</math>.

Because the partition function <math>Z=\int D\Delta^{*}D\Delta e^{-S_{eff}(\Delta)}</math>, if we only consider the <math>\Delta</math> related factors.

The superconductivity phase transition temperature is the temperature makes <math>\frac{1}{\left|g\right|}-N(0)ln(\frac{2\hbar\omega_{D}e^{\gamma}}{\pi k_{B}T})=0</math>, which is <math>T_{c}=\frac{\hbar\omega_{D}}{k_{B}}\frac{2}{\pi}e^{\gamma}e^{-\frac{1}{N(0)\left|g\right|}}=1.134\frac{\hbar\omega_{D}}{k_{B}}e^{-\frac{1}{N(0)\left|g\right|}}</math>.

Beyond the critical temperature, the <math>\Delta</math> related factors in the partition function is just <math>1</math>, the same as no cooper pair, which is normal state; below the critical temperature, the <math>\Delta</math> related factors in the partition function will diverge, which means superconductivity phase transition.

Finite <math>\vec{q}</math> (small) <math>\ (\Omega_n=0)</math>

<math>\chi_p (q,0)-\chi_p (0,0)=\frac{1}{L^D} \sum_k \frac{1}{\beta} \sum_{i\omega_n}\frac{-1}{i\omega_n-\epsilon_k+\mu}(\frac{1}{i\omega_n+\epsilon_{q-k}-\mu}-\frac{1}{i\omega_n+\epsilon_{-k}-\mu}) </math>

for small <math>\vec{q}</math>,
<math>\epsilon_{q-k}=\epsilon_{-k}+q\frac{\partial \epsilon_\rho}{\partial \rho}|_{\rho=-k}=\epsilon_\vec{k}+\vec{q}\cdot \vec{v}_{-k}=\epsilon_k-\vec{q}\cdot\vec{v}_k</math>

and <math>\frac{1}{i\omega_n+\epsilon_k-\mu-q\cdot v_k}=\frac{1}{i\omega_n+\epsilon_k-\mu}+\frac{q\cdot v_k}{(i\omega_n+\epsilon_k-\mu)^2}+\frac{(q\cdot v_k)^2}{(i\omega_n+\epsilon_k-\mu)^3}</math>

Thus, 

<math>\begin{align}\chi_p(q,0)-\chi_p(0,0)
&=\frac{1}{L^D}\sum_k\frac{1}{\beta}\sum_{i\omega_n}\frac{-1}{i\omega_n-\epsilon_k+\mu}\left(\frac{q\cdot v_k}{(i\omega_n+\epsilon_k-\mu)^2}+\frac{(q\cdot v_k)^2}{(i\omega_n+\epsilon_k-\mu)^3}+...\right) \\
&=\frac{-1}{\beta}\sum_i\omega_n \int \frac{d^D k}{(2\pi)^D}\frac{(\vec{q}\cdot \vec{v}_k)^2}{(i\omega_n-\epsilon_k+\mu)(i\omega_n+\epsilon_k-\mu)^3} 
\end{align}</math>

Consider the states near the <math>\hbar \omega_D</math> shell near fermi surface, we have

<math>\chi_p(q,0)-\chi_p(0,0) =\frac{1}{\beta}\sum_{i\omega_n}\int\frac{d\Omega_{F.S.}}{\Omega_D}(q\cdot v_F)^2\int_{-\infty}^{+\infty} d\xi N(\xi+\mu)\frac{1}{(\xi-i\omega_n)(\xi+i\omega_n)^3}</math>

where, <math>\ \xi=\epsilon_k-\mu</math>

and 

<math>
\begin{align}
\int_{-\infty}^{+\infty} d\xi \frac{1}{(\xi-i\omega_n)(\xi+i\omega_n)^3}
&=\frac{2\pi i}{(2i\omega_n)^3}\theta(\omega_n)-\frac{2\pi i}{(2i\omega_n)^3}\theta(-\omega_n)\\
&=\frac{2\pi i}{(2i|\omega|)^3}
\end{align}
</math>

So, 

<math>
\begin{align}
\chi_p(q,0)-\chi_p(0,0)
&=\frac{1}{\beta}\sum_{i\omega_n}N(0)\int\frac{d\Omega}{\Omega_D}(q\cdot v_F)^2\frac{2\pi i}{(2i|\omega|)^3}\\
&=N(0)v_F^2|\vec{q}|^2\int\frac{d\Omega}{\Omega_D}(q\cdot v_F)^2 \frac{1}{\beta}\sum_{i\omega}\frac{2\pi i}{-i8|\frac{(2n+1)\pi}{\beta}|^3}\\
&=-\frac{1}{4}N(0)v_F^2q^2(<(\hat{q}\cdot \hat{v_F})>_{F.S.})\frac{\beta^2}{\pi^2}(\sum_{N=-\infty}^{+\infty}\frac{1}{|2n+1|^3})
\end{align}
</math>

<math>
\sum_{N=-\infty}^{+\infty}\frac{1}{|2n+1|^3}=\sum_{n=0}^\infty\frac{2}{(2n+1)^3}=\frac{2}{\pi}\frac{7\zeta(3)}{8}
</math>

where, <math>\ \zeta(3)</math> is Riemann zeta function.

<math>
<(\hat{q}\cdot \hat{v}_F)^2>_{F.S.}=\frac{1}{D}
</math>

For spherical F.S. in 3D,

<math>
\int\frac{d\Omega}{\Omega_D}(\hat{q}\cdot\hat{v}_F)^2=\frac{2\pi}{4\pi}\int_{-1}^{1}dcos\theta cos^2\theta = \frac{1}{3}
</math>

For circular F.S. in 2D,

<math>
\int\frac{d\Omega}{\Omega_D}(\hat{q}\cdot\hat{v}_F)^2=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta cos^2\theta = \frac{1}{2}
</math>

Then

<math>
\begin{align}
\chi_p(q,0)-\chi_p(0,0)
&=-\frac{1}{4}N(0)v_{F}^{2}q^{2}\frac{1}{D}\frac{\beta^{2}}{\pi^{2}}\frac{2}{\pi}\frac{7\zeta(3)}{8} \\
&=-N(0)\frac{7\zeta(3)}{16D\pi^{2}}q^{2}\frac{1}{\pi \hbar^{2}}\left(\frac{\hbar v_{F}}{k_{B}T}\right)^{2} \\ &\equiv-N(0)q^{2}\xi^{2}
\end{align}
</math>


So 

<math>
\begin{align}
\frac{1}{2}\left\langle S_{int.}^{2}\right\rangle _{0}&=L^{D}\frac{1}{\beta}\underset{\Omega_{n},\mathbf{q}}{\sum}\chi_{p}(q,0)\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n}) \\
&=N(0)ln[\frac{2\hbar\omega_{D}e^{\gamma}}{\pi k_{B}T}]L^{D}\frac{1}{\beta}\underset{\Omega_{n},\mathbf{q}}{\sum}\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})-L^{D}\frac{1}{\beta}\underset{\Omega_{n},\mathbf{q}}{\sum}N(0)q^{2}\xi^{2}\Delta_{\mathbf{q}}^{*}(i\Omega_{n})\Delta_{\mathbf{q}}(i\Omega_{n})
\end{align}
</math>.

<math>S_{eff}=\underset{\mathbf{r}}{\sum}\int_{0}^{\beta}d\tau\left[\left(\frac{1}{\left|g\right|}-N(0)ln(\frac{2\hbar\omega_{D}e^{\gamma}}{\pi k_{B}T})\right)\Delta^{*}(\tau,\mathbf{r})\Delta(\tau,\mathbf{r})-N(0)\xi^{2}(\nabla\cdot\Delta^{*}(\tau,\mathbf{r}))(\nabla\cdot\Delta(\tau,\mathbf{r}))\right]</math>.

Note that the last term in the expression tells us that <math> S_{eff} </math> would increase if gradient of <math> \Delta </math> is not zero.

Note that the above expression has a one-one correspondant to the Giznburg-Landau functional:

<math> F=\int d^{D}r\left[ \alpha (T-T_{c}) |\Psi(\vec{r})|^{2}+\frac{\hbar^{2}}{2m^{*}}|\nabla \Psi(\vec{r})|^{2} \right] </math>,

here <math>\ \Psi(\vec{r}) </math> corresponds to <math>\ \frac{\Delta(\tau,\vec{r})}{|g|N(0)a_{0}} </math> in <math>S_{eff}</math>. 


=== Little Parks experiment ===

[[Image:Layout48.gif]]

Refer to the fig, a thin shell of superconductor with radius R is shown and a small uniform magnetic field is passing through the hollow center of the cylinder. The experiment intends to show the variation of the critical temperature with change of the magnetic field passing through the hollow superconductor cylinder. 

Before showing it, we first have to rewrite the Giznburg-Landau functional to make it taken the presence of magnetic field into account. Hamiltonian for a free electron moving in a magnetic field can be written as:

<math> \frac{1}{2m}(p+\frac{eA}{c})^{2}\psi + V\psi = E\psi </math>

The physical observable magnetic field B would remain the same if we choose a different vector potential
<math> A\rightarrow A+ \nabla \chi </math> (ie perform gauge transformation). To maintain the same eigen-energy E which is observable, the wave function have to undergo a phase change:
<math> \psi \rightarrow e^{i\phi}\psi </math> where <math> \phi=\frac{e}{c\hbar}\chi </math>

Now in our Hamiltonian, the wave function is arranged as  
<math> \Delta^{*}(\tau,\vec{r})\psi_\uparrow (\tau,\vec{r})\psi_\downarrow (\tau,\vec{r}) + \Delta(\tau,\vec{r}) \psi_\downarrow^\dagger (\tau,\vec{r})\psi_\uparrow^\dagger (\tau,\vec{r}) </math>

since <math> \psi \rightarrow e^{i\phi}\psi </math>, so if we want the Hamiltonian to remind the same, <math>\ \Delta </math> has to transform as <math> \Delta \rightarrow e^{-2i\phi}\Delta  </math>

Since <math>\ \Delta </math> corresponds to <math>\ \Psi </math> in the Giznburg-Landau functional, so the Giznburg-Landau functional is modified as 

<math> F=\int d^{D}r\left[ \alpha (T-T_{c}) |\Psi(\vec{r})|^{2}+\frac{1}{2m^{*}}| ( \frac{\hbar \nabla}{i} - \frac{2e}{c}A(\vec{r}) ) \Psi(\vec{r})|^{2} \right] </math>

choose symmetric gauge:
<math> \vec{A}=\frac{1}{2}\vec{H}\times\vec{r}=\frac{1}{2}Hr\hat{\phi} </math>

In cylindrical coordinate: <math> \vec{\nabla}=\hat{r}\frac{\partial}{\partial r} + \frac{\hat{\phi}}{r}\frac{\partial}{\partial \phi} + \hat{z}\frac{\partial}{\partial z} </math>

define unit flux as <math> \Phi_{0}=\frac{hc}{2e} </math> 

define fluxoid as <math> \Phi(R) = \pi HR^{2}\ </math>, so we have 

<math> \begin{align} F&=\int d^{D}r\left[ \alpha (T-T_{c})|\Psi(\vec{r})|^{2} 
+\frac{\hbar^{2}}{2m^{*}}| (\frac{1}{R}\frac{\partial}{\partial \phi} - \frac{ie}{\hbar c} HR )\Psi(\vec{r}) |^{2}+
\frac{\hbar^{2}}{2m^{*}}|\frac{\partial}{\partial z} \Psi(\vec{r}) |^{2} \right] \\
&=\int d^{D}r\left[ \alpha (T-T_{c})|\Psi(\vec{r})|^{2} 
+\frac{\hbar^{2}}{2m^{*}R^{2}}| (\frac{\partial}{\partial \phi} - \frac{i\Phi}{\Phi_{0}} )\Psi(\vec{r}) |^{2} \right] \\
 \end{align} </math>

When <math> \Phi = N\Phi_{0}\ </math>, the critical temperature will remain the same and the phase of <math> \Psi\ </math> is changed as <math> \Psi \rightarrow e^{iN\phi} \Psi </math>. When  <math> \Phi \neq N\Phi_{0}\ </math>, the critical temperature is found to vary as
<math> T_{c}^{new}=T_{c}- \frac{\hbar^{2}}{2m^{*}R^{2}\alpha}\left (N-\frac{\Phi}{\Phi_{0}}\right )^{2}</math>. See the fig.

[[Image:Graph49.gif]]

== Microscopic derivation of the Giznburg-Landau functional ==

<math>Z=Z_{0}< e^{-S_{int}} ></math>

where
<math>Z_{0}=\int D\psi ^{*} D\psi D\Delta ^{*} D\Delta  e^{-(S_{\Delta} +S_{0})}</math>

we can expand this average for small<math>\Delta</math> near<math>T_{c}</math>, for this perpose we can assume  asecond order phase transition 
so that it increases continiously from zero to finite number after <math>T_{c}</math>

we need to calculate the average of <math>e^{-s_{int}}</math>which can be calculated by Tylor expansion:

<math>e^{-S_{int}}=<-S_{int}+\frac{1}{2}S_{int}^{2}-\frac{1}{3}S_{int}^{3}+\frac{1}{4!}S_{int}^{4}+...></math>

=<math>1-<S_{int}>+\frac{1}{2} < S_{int}^{2}> -\frac{1}{3!}< S_{int}^{2}> +\frac{1}{4!}< S_{int}^{4}> +...</math>

.......................the odd power terms are zero because <math><\psi _{\uparrow}(r,\tau )\psi _{\downarrow}(r,\tau ) > =0 </math>

<math>=e^{\frac{1}{2}< S_{int}^{2}>}e^{\frac{1}{4!}< S_{int}^{4}>-\lambda }</math>

if we expand these two terms in to the second order the following expression can be got:

<math>(1+\frac{1}{2} < S_{int}^{2}>+\frac{1}{2}(<\frac{1}{2} S_{int}^{2}>)^{2} +...)(1+\frac{1}{4!}< S_{int}^{4}>+...)</math>

<math>=1+\frac{1}{2} < S_{int}^{2}>+\frac{1}{8}(< S_{int}^{2}>)^{2} +...)+\frac{1}{4!}< S_{int}^{4}>-\lambda +...</math>

<math>\lambda</math> can be choosed in such a way .......

so,
<math>\lambda =\frac{1}{8}< S_{int}^{2}> ^{2} </math>

<math>=e^{\frac{1}{2} < S_{int}^{2}>+\frac{1}{4!}(< S_{int}^{4}>-3<S_{int}^{2}>^{2})+...}</math>

according to the expression we got before:


<math>S_{int}=\frac{L^{D}}{\beta ^{2}}\sum_{\omega _{n},\Omega _{n}}\sum _{k,q}[\Delta ^{*}_{q}(i\Omega _{n})\psi_{\downarrow}(i\Omega _{n}-i\omega _{n}),\vec{q}-\vec{k})\psi_{\uparrow}(i\omega _{n},k)+\Delta _{q}(i\Omega_{n})\psi_{\uparrow}^{\dagger }(i\omega _{n},k)\psi_{\downarrow}^{\dagger }(i\Omega _{n}-i\omega _{n}),\vec{q}-\vec{k})]</math>

let's write <math>S_{int}</math> in terems od <math>a</math> for simplification. where

<math>a=\int \Delta ^{*}(1)\psi_{\downarrow}(1) \psi_{\uparrow}(1)+\Delta (1)\psi_{\downarrow}^{*}(1) \psi_{\uparrow}^{*}(1)</math>

<math>a_{1}</math> is a couple grassman number, so we do not need to be worry about the sign when these terms comute with other terms.

<math>< S_{int}^{4}> =\int_{1234} < (a_{1}^{*}+a_{1})(a_{2}^{*}+a_{2})(a_{3}^{*}+a_{3})(a_{4}^{*}+a_{4})></math> 

<math>=(< a_{1}^{*}a_{2}^{*}a_{3}a_{4}>+ < a_{1}^{*}a_{2}a_{3}^{*}a_{4}>+< a_{1}^{*}a_{2}a_{3}a_{4}^{*}>+< a_{1}a_{2}^{*}a_{3}^{*}a_{4}> +< a_{1}a_{2}^{*}a_{3}a_{4}^{*}>)</math>

<math>=6< a_{1}^{*}a_{2}^{*}a_{3}a_{4}>=6\int _{1234}\Delta ^{*}(1)\Delta ^{*}(2)\Delta(3)\Delta(4)< \psi_{\downarrow}(1)\psi_{\uparrow}(1)\psi_{\downarrow}(2)\psi_{\uparrow}(2)\psi_{\downarrow}^{*}(3)\psi_{\uparrow}^{*}(3)\psi_{\downarrow}^{*}(4)\psi_{\uparrow}^{*}(4)> </math>

<math>3< S_{int}^{2}> ^{2}=3\int_{1,2}2< a_{1}^{*}a_{2}> \int_{3,4}2< a_{3}^{*}a_{4}>=12\int_{1,2,3,4}\Delta ^{*}(1)\Delta ^{*}(2)\Delta(3)\Delta(4)< \psi_{\downarrow}(1)\psi_{\uparrow}(1)\psi_{\uparrow}^{*}(3)\psi_{\downarrow}^{*}(3)>< \psi_{\downarrow}(2)\psi_{\uparrow}(2)\psi_{\uparrow}^{*}(4)\psi_{\downarrow}^{*}(4)></math> 



<math>-2G(2-3) G(2-4)G(1-4)G(1-3)=-12\int_{1,2,3,4}\Delta ^{*}(1)\Delta _{*}(2)\Delta (3)\Delta (4)G(2-3)G(2-4)G(1-4)G(1-3)</math>

Recall 

<math>G(2-3)=< \psi (r_{2},\tau _{2})\psi ^{*}(r_{3},\tau _{3})>=\frac{1}{\beta } \sum_{\omega _{n}}\frac{1}{L^{D}}\sum_{k}e^{-i\omega _{n}(\tau _{2}-\tau _{3})}e^{ik.(r_{2}-r_{3})}\frac{1}{-i\omega _{n}+\epsilon _{k}-\mu }</math>

see the solution which are <math>\tau</math>independent

<math>\frac{1}{\beta ^{4}}\sum_{\omega _{{n}_{1}}}...\sum_{\omega _{{n}_{4}}}\int_{0}^\beta {d\tau_{1}}\int_{0}^\beta {d\tau_{2}} \int_{0}^\beta {d\tau_{3}} \int_{0}^\beta {d\tau_{4}}e^{-i\omega _{{n}_{1}}(\tau _{1}-\tau _{3})} e^{-i\omega _{{n}_{2}}(\tau _{1}-\tau _{4})}e^{-i\omega _{{n}_{3}}(\tau _{2}-\tau _{3})}e^{-i\omega _{{n}_{4}}(\tau _{2}-\tau _{4})} G(i\omega _{{n}_{1}},r_{1}-r_{3})G(i\omega _{{n}_{2}},r_{1}-r_{4})G(i\omega _{{n}_{3}},r_{2}-r_{3})G(i\omega _{{n}_{4}},r_{2}-r_{4}) 
</math>

after getting integration over <math>\tau_{1}</math> we will get <math>\beta \delta (\omega _{n_{1}},-\omega _{n_{2}}) </math> and similarly by getting integration over <math>\tau_{2}</math> we have <math>\beta \delta (\omega _{n_{3}},-\omega _{n_{4}}) </math>

So, the final result can be written:
<math>\sum_{i\omega _{n}}G(i\omega _{n},r_{1}-r_{3}G(-i\omega _{n},r_{2}-r_{3})G(-i\omega _{n},r_{1}-r_{4})G(i\omega _{n},r_{2}-r_{4}) </math>

Now, we wish to perform gradiant expansion:

<math>\Delta ^{\ast }(1)=\Delta ^{\ast }(\frac{r_{1}+r_{4}}{2}+\frac{r_{1}-r_{4}}{2})=\Delta ^{\ast }(\frac{r_{1}+r_{4}}{2})+(\frac{r_{1}-r_{4}}{2})\nabla\Delta^{\ast }(\frac{r_{1}+r_{4}}{2})+... </math>

<math>\Delta ^{\ast }(2)=\Delta ^{\ast }(\frac{r_{2}+r_{3}}{2}+\frac{r_{2}-r_{3}}{2})=\Delta ^{\ast }(\frac{r_{2}+r_{3}}{2})+(\frac{r_{2}-r_{3}}{2})\nabla\Delta^{\ast }(\frac{r_{2}+r_{3}}{2})+...  </math>
<math>-12\int d^{D}R_{1,4}d^{D}R_{2,3}d^{D}\mu _{1,4}d^{D}\mu_{2,3}\Delta ^{\ast }(R_{1,4})\Delta ^{\ast }(R_{2,3})\Delta (R_{2,3}) \Delta (R_{1,4})\sum_{\omega_{n}}G(i\omega_{n},R_{1,4}-R_{2,3}+\frac{1}{2}(\mu _{1,4}+\mu _{2,3}))</math>

where:

<math>R=\frac{1}{2}(R_{1,4}+R_{2,3})</math>
          
<math>\mu =R_{1,4}-R_{2,3} </math>

<math>r_{1}-r_{3}=\frac{1}{2}(\mu _{1,4}+\mu_{2,3})+\mu</math>

<math>r_{2}-r_{4}=\frac{1}{2}(\mu _{1,4}+\mu_{2,3})-\mu</math>




Starting from the microscopic model, we found that <math>Z\backsimeq Z_{0}\int D\Delta*D\Delta e^{-S_{eff}} </math>,
where the <math>4^{th}</math> order in <math>\Delta</math>, and keeping only quadratic qradient terms, we have:

<math>S_{eff}=\frac{1}{k_{B}T}\sum_{r}\left[\underset{A}{\underbrace{\left(\frac{1}{|g|}-N(0)In\left[\frac{2\hbar\omega_{D}e^{\gamma_{E}}}{\pi k_{B}T}\right]\right)}}|\Delta(r)|^{2}+N(0)\xi^{2}(\nabla\Delta*(r)).(\nabla\Delta(r))+\frac{1}{2}\underset{B}{\underbrace{\frac{7\zeta(3)N(0)}{8\pi^{2}k_{B}^{2}T^{2}}}}|\Delta(r)|^{4}\right]</math>

We can use this expression to make quantitative experimental predictions. The path integral over <math>\Delta</math> is still imposible to carry out exactly, despite our approximations for <math>S_{eff}</math>, because <math>S_{eff}</math>contains quartic terms in <math>\Delta</math> and so we are not dealing with a Gaussian integral. The approximation strategy whic we will pursue is called saddle point approxiation, which in our contetxt means that we will expand teh integrand about a solution which minimizes <nowiki>S_{eff}</nowiki> with respect to <math>\Delta</math>. What we end up doing is replacing Z with
 <math>Z\sim Z_{0}e^{-S_{eff}[\Delta_{min}]}</math>,
 where <math>\Delta_{min}</math> is determined from<math> \left|\frac{\delta S_{eff}}{\delta\Delta^{*}}\right|_{\Delta=\Delta_{min}}=0</math> 
At this point, let's seek uniform solutions to their equations, in whcih case we can drop the gradient terms in <math>S_{eff}</math> :
<math>\frac{\delta S_{eff}}{\delta\Delta^{*}}=\frac{L^{D}}{k_{B}T}\left(A\Delta+B\Delta\Delta^{*}\Delta\right)</math>
where: 
<math>A=\frac{1}{|g|}-N(0)ln\left(\frac{\hbar\omega_{D}}{k_{B}T}\frac{2e\gamma}{\pi}\right)</math> and 
<math>B=\frac{7\zeta(3)}{8\pi}\frac{N(0)}{k_{B}^{2}T^{2}}</math>
Note that for <math>T>T_{C}\qquad A>0,B>0</math>

and <math>T<T\qquad A<0,B>0</math>

So  <math>A\Delta_{min}+B\Delta_{min}^{3}=0</math>

<math>\Delta_{min}=0,\: T>T_{C}</math>

and <math>\Delta_{min}=\sqrt{\frac{-A}{B}},\: T<T_{C} So, S_{eff}[\Delta_{min}]=\frac{L^{D}}{k_{B}T}\left(A\Delta_{min}^{2}+\frac{1}{2}B\Delta_{min}^{4}\right)</math>

<math>T>T_{C},\quad S_{eff}[\Delta_{min}]=0T<T_{C},\quad S_{eff}[\Delta_{min}]=\left(A\frac{-A}{B}+\frac{1}{2}B\frac{A^{2}}{B^{2}}\right)\frac{L^{D}}{k_{B}T}=\frac{-A^{2}}{2B}\frac{L^{D}}{k_{B}T}</math>

Since, we now have the approximate expression for the partition function we can calculate thermodynamic physical properties. the one we will focus on is the specific heat. Recall that, <math>Z=e^{-\beta F}=\sum_{n}e^{-\beta E_{n}}\Longrightarrow<E>=\frac{1}{Z}\sum_{n}E_{n}e^{-\beta E_{n}}=-\frac{\partial}{\partial\beta}lnZ =-\frac{\partial}{\partial\beta}lne^{-\beta F}=\frac{\partial}{\partial\beta}\left(\beta F\right)=F+\beta\frac{\partial F}{\partial\beta}</math>

== Effects of an applied magnetic field; Type I and Type II superconductivity ==

=== Derivation of the Ginzburg-Landau equations ===

Our starting point will be the Ginzburg-Landau (GL) free energy in the presence of an external magnetic field,

<math>F=\int d^d \vec{r} \left [\alpha(T-T_c)|\Psi(\vec{r})|^2+\tfrac{1}{2}b|\Psi(\vec{r})|^4+\frac{\hbar^2}{2m}\left |\left (\nabla-i\frac{2e}{\hbar c}\vec{A}(\vec{r})\right )\Psi(\vec{r})\right |^2+\frac{1}{8\pi}(\nabla\times\vec{A}(\vec{r}))^2-\frac{1}{c}\vec{J}_{\text{ext}}(\vec{r})\cdot\vec{A}(\vec{r})\right ],</math>

where <math>\vec{A}</math> is the total vector potential and <math>\vec{J}_{\text{ext}}</math> is an external current density, assumed to be controlled experimentally.  This current satisfies

<math>\nabla\times\vec{H}=\frac{4\pi}{c}\vec{J}_{\text{ext}},</math>

where <math>\vec{H}</math> is the external magnetic field.  The expression is the sum of the energy due to the superconducting order parameter, with the magnetic field introduced via the gauge invariance argument given above, the energy of the magnetic field alone, and the work done by the superconductor to maintain the external current at a constant value.

Let us first derive the "saddle point" equations satisfied by the magnetic field in the normal state.  In this case, we set <math>\Psi</math> to zero everywhere and set

<math>\left. \frac{\delta F}{\delta\vec{A}(\vec{r})}\right |_{\vec{A}=\vec{A}_\text{min}}=0.</math>

We will find this derivative by first finding the variation <math>\delta F</math> in the free energy for this case, which is

<math>\delta F=\int d^d \vec{r}' \left [\frac{1}{4\pi}(\nabla\times\vec{A}(\vec{r}'))\cdot(\nabla\times\delta\vec{A}(\vec{r}'))-\frac{1}{c}\vec{J}_{\text{ext}}(\vec{r}')\cdot\delta\vec{A}(\vec{r}')\right ],</math>

where <math>\delta\vec{A}</math> is a small variation in the vector potential; we assume that it vanishes on the "surface" of our system.  We now transform the first term using the identity,

<math>(\nabla\times\vec{A})\cdot(\nabla\times\vec{B})=\nabla\cdot[\vec{A}\times(\nabla\times\vec{B})]+\vec{A}\cdot[\nabla\times(\nabla\times\vec{B})],</math>

obtaining

<math>\delta F=\int d^d \vec{r}' \left [\frac{1}{4\pi}\nabla\cdot[\delta\vec{A}(\vec{r}')\times(\nabla\times\vec{A}(\vec{r}'))]+\frac{1}{4\pi}\delta\vec{A}(\vec{r}')\cdot[\nabla\times(\nabla\times\vec{A}(\vec{r}'))]-\frac{1}{c}\vec{J}_{\text{ext}}(\vec{r}')\cdot\delta\vec{A}(\vec{r}')\right ].</math>

The first term is a "surface" term; since we assumed that <math>\delta\vec{A}</math> vanishes everywhere on the "surface", we are left with just

<math>\delta F=\int d^d \vec{r} \left [\frac{1}{4\pi}[\nabla\times(\nabla\times\vec{A}(\vec{r}'))]-\frac{1}{c}\vec{J}_{\text{ext}}(\vec{r}')\right ]\cdot\delta\vec{A}(\vec{r}').</math>

We conclude that the variational derivative that we are interested in is

<math>\frac{\delta F}{\delta\vec{A}(\vec{r})}=\frac{1}{4\pi}[\nabla\times(\nabla\times\vec{A}(\vec{r}'))]-\frac{1}{c}\vec{J}_{\text{ext}}(\vec{r}').</math>

At the "saddle point", this derivative is zero, so we obtain the equation,

<math>\nabla\times(\nabla\times\vec{A})=\frac{4\pi}{c}\vec{J}_{\text{ext}}.</math>

We may introduce the total magnetic field <math>\vec{B}=\nabla\times\vec{A}</math>, thus obtaining

<math>\nabla\times\vec{B}=\frac{4\pi}{c}\vec{J}_{\text{ext}}.</math>

Comparing this to the definition of <math>\vec{J}_{\text{ext}}</math> given above, we conclude that <math>\vec{B}=\vec{H}</math> in the normal state.  In reality, this will only be approximately true due to para- or diamagnetic effects in the metal, but these effects will be small in comparison to those due to superconductivity, which we will now derive.

First, we will apply the "saddle point" condition for the superconducting order parameter, <math>\Psi</math>, which is

<math>\left. \frac{\delta F}{\delta\Psi^{*}(\vec{r})}\right |_{\Psi=\Psi_{\text{min}}}=0.</math>

Again, we start by finding the variation in the free energy in terms of a small variation <math>\delta\Psi^{*}</math> in the order parameter:

<math>\delta F=\int d^d \vec{r}' \left [\alpha(T-T_c)\Psi(\vec{r}')\,\delta\Psi^{*}(\vec{r}')+b|\Psi(\vec{r}')|^2\Psi(\vec{r}')\,\delta\Psi^{*}(\vec{r}')-\frac{e}{mc}\vec{A}(\vec{r}')\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')\,\delta\Psi^{*}(\vec{r}')-\frac{1}{2m}\frac{\hbar}{i}\nabla\delta\Psi^{*}(\vec{r}')\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')\right ]</math>

The last term is equal to

<math>-\frac{1}{2m}\left\{\frac{\hbar}{i}\nabla\cdot\left [\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')\,\delta\Psi^{*}(\vec{r}')\right ]-\left [\frac{\hbar}{i}\nabla\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}
')\right ]\delta\Psi^{*}(\vec{r}')\right\}.</math>

The second term in this expression is a "surface" term.  If we assume that <math>\delta\Psi^{*}</math> is zero on the "surface", then this term vanishes, leaving us with

<math>\delta F=\int d^d \vec{r}' \left\{\alpha(T-T_c)\Psi(\vec{r}')+b|\Psi(\vec{r}')|^2\Psi(\vec{r}')-\frac{e}{mc}\vec{A}(\vec{r}')\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')+\frac{1}{2m}\left [\frac{\hbar}{i}\nabla\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')\right ]\right\}\delta\Psi^{*}(\vec{r}').</math>

We can now immediately write down the variational derivative, which, upon being set to zero, gives us the first GL equation,

<math>\frac{1}{2m}\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}\right )^2\Psi+\alpha(T-T_c)\Psi+b|\Psi|^2\Psi=0.</math>

We also need to minimize the free energy with respect to the magnetic field.  We have already done this for the normal case, and there is only one more term that we need to consider in the superconducting case; we will therefore only treat this term.  We can quickly write down the variation in the superconducting part of the free energy <math>F_{SC}</math>, which is

<math>\delta F_{SC}=i\frac{e\hbar}{mc}\int d^d \vec{r}' \left [\Psi^{*}(\vec{r}')\left (\nabla-i\frac{2e}{\hbar c}\vec{A}(\vec{r}')\right )\Psi(\vec{r}')-\Psi(\vec{r}')\left (\nabla+i\frac{2e}{\hbar c}\vec{A}(\vec{r}')\right )\Psi^{*}(\vec{r}')\right ]\cdot\delta\vec{A}(\vec{r}').</math>

Combining this result with the previous result for the normal metal, we obtain the second GL equation,

<math>\frac{1}{4\pi}\nabla\times(\nabla\times\vec{A})-\frac{1}{c}\vec{J}_{\text{ext}}-\frac{e}{mc}\left (\Psi^{*}\frac{\hbar}{i}\nabla\Psi-\Psi\frac{\hbar}{i}\nabla\Psi^{*}\right )+\frac{4e^2}{mc^2}|\Psi|^2\vec{A}=0,</math>

or, introducing <math>\vec{B}</math> and <math>\vec{H}</math>,

<math>\frac{e}{m}\left (\Psi^{*}\frac{\hbar}{i}\nabla\Psi-\Psi\frac{\hbar}{i}\nabla\Psi^{*}\right )-\frac{4e^2}{mc}|\Psi|^2\vec{A}=\frac{c}{4\pi}\nabla\times(\vec{B}-\vec{H}).</math>

Given the definition of <math>\vec{H}</math> and the Maxwell equation (assuming static fields),

<math>\nabla\times\vec{B}=\frac{4\pi}{c}\vec{J},</math>

where <math>\vec{J}</math> is the total current density, we conclude that the left-hand side of this equation is the current density induced inside the superconductor.

Let us now suppose that we do not assume that <math>\delta\Psi^{*}</math> vanishes on the surface.  It may then be shown that the following boundary condition holds on the surface (see P. G. de Gennes, <i>Superconductivity in Metals and Alloys</i>):

<math>\hat{n}\cdot\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}\right )\Psi=\frac{i\hbar}{b_{dG}}\Psi.</math>

This relation holds for a superconductor-metal interface; for a superconductor-insulator interface, <math>b_{dG}\rightarrow\infty</math>.  We may show that this condition implies that the normal component of the current density on the surface vanishes.  If we multiply the above condition by <math>\Psi^{*}</math> on both sides, we obtain

<math>\hat{n}\cdot\Psi^{*}\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}\right )\Psi=\frac{i\hbar}{b_{dG}}\Psi^{*}\Psi.</math>

Taking the complex conjugate of both sides gives us

<math>\hat{n}\cdot\Psi\left (-\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}\right )\Psi^{*}=-\frac{i\hbar}{b_{dG}}\Psi^{*}\Psi.</math>

Adding these two equations together gives us

<math>\hat{n}\cdot\left [\left (\Psi^{*}\frac{\hbar}{i}\nabla\Psi-\Psi\frac{\hbar}{i}\nabla\Psi^{*}\right )-\frac{4e}{c}|\Psi|^2\vec{A}\right ]=0.</math>

The left-hand side is proportional to the normal component of the current density inside the superconductor.

=== The GL Equations in Dimensionless Form ===

We will find it convenient to introduce dimensionless variables when working with the GL equations.  We start by introducing a dimensionless order parameter, <math>\psi=\frac{\Psi}{\Psi_0}</math>, where

<math>\Psi_0^2=\frac{\alpha(T_c-T)}{b}.</math>

We may rewrite the first GL equation in terms of this parameter as

<math>\frac{1}{2mb\Psi_0^2}\left (\frac{\hbar}{i}\nabla-\frac{2e}{c}\vec{A}\right )^2\psi-(\left |\psi\right |^2-1)\psi=0,</math>

and the second as

<math>\frac{e\Psi_0^2}{mc}\left (\psi^{\ast}\frac{\hbar}{i}\nabla\psi-\psi\frac{\hbar}{i}\nabla\psi^{\ast}\right )-\frac{4e^2\Psi_0^2}{mc^2}\left |\psi\right |^2\vec{A}=\frac{1}{4\pi}\nabla\times[\nabla\times(\vec{A}-\vec{A}_0)],</math>

where we re-introduced <math>\vec{A}</math> into the right-hand side and also introduced <math>\vec{A}_0</math>, defined as

<math>\vec{H}=\nabla\times\vec{A}_0.</math>

Next, we introduce a dimensionless position vector,

<math>\tilde{r}=\frac{1}{\lambda}\vec{r},</math>

where <math>\lambda=\sqrt{\frac{mc^2}{16\pi e^2\Psi_0^2}}</math> is known as the penetration depth of the superconductor; we will see where this name comes from shortly.  In terms of this vector, the first GL equation becomes

<math>\left (\frac{1}{\Psi_0\lambda\sqrt{2mb}}\frac{\hbar}{i}\tilde{\nabla}-\frac{1}{\Psi_0\sqrt{2mb}}\frac{2e}{c}\vec{A}\right )^2\psi+(\left |\psi\right |^2-1)\psi=0</math>

and the second becomes

<math>\frac{4\pi e\lambda\Psi_0^2}{mc}\left (\psi^{\ast}\frac{\hbar}{i}\tilde{\nabla}\psi-\psi\frac{\hbar}{i}\tilde{\nabla}\psi^{\ast}\right )-\left |\psi\right |^2\vec{A}=\tilde{\nabla}\times[\tilde{\nabla}\times(\vec{A}-\vec{A}_0)].</math>

Finally, we introduce a dimensionless vector potential,

<math>\tilde{A}=\frac{1}{\Psi_0\sqrt{2mb}}\frac{2e}{c}\vec{A}</math>

and the dimensionless parameter,

<math>\kappa=\frac{\Psi_0\lambda\sqrt{2mb}}{\hbar}.</math>

In terms of these, the first GL equation becomes

<math>\left (-\frac{i}{\kappa}\tilde{\nabla}-\tilde{A}\right )^2\psi+(\left |\psi\right |^2-1)\psi=0</math>

and the second becomes

<math>\frac{1}{2\kappa}\left (\psi^{\ast}\frac{\tilde{\nabla}}{i}\psi-\psi\frac{\tilde{\nabla}}{i}\psi^{\ast}\right )-\left |\psi\right |^2\vec{A}=\tilde{\nabla}\times[\tilde{\nabla}\times(\tilde{A}-\tilde{A}_0)].</math>

We see that our theory has a dimensionless parameter in it, namely <math>\kappa</math>, which is known as the Ginzburg-Landau parameter.  We may write this parameter as

<math>\kappa=\frac{\lambda\sqrt{2m\alpha(T_c-T)}}{\hbar}=\frac{\lambda}{\xi_\text{GL}},</math>

where

<math>\xi_\text{GL}=\frac{\hbar}{2m\alpha(T_c-T)}</math>

is the GL coherence length.  This tells us that <math>\kappa</math> is the ratio of two length scales associated with the superconductor, namely the scale over which the order parameter "heals" (the coherence length <math>\xi_\text{GL}</math>) and that over which the magnetic field dies out (the penetration depth <math>\lambda</math>, as we will demonstrate shortly).  It also turns out that this parameter decides what type of superconductor we are dealing with.  If <math>\kappa < \tfrac{1}{\sqrt{2}}</math>, then we have a Type I superconductor, while, if <math>\kappa > \tfrac{1}{\sqrt{2}}</math>, then we have a Type II superconductor.

We may now find the value of this parameter in the microscopic model we considered earlier.  In that case, we found that

<math>\frac{\hbar^2}{2m}=N(\mu)\xi^2=N(\mu)\frac{7\zeta(3)}{16\pi^2 d}l_T^2,</math>

where <math>N(\mu)</math> is the density of states at the Fermi level, <math>\xi</math> is the coherence length, <math>d</math> is the number of dimensions that we are working in, and <math>l_T</math> is the thermal wavelength.  We will state the result for <math>d=3</math>.  Given that

<math>l_T=\frac{\hbar v_F}{k_B T}</math>

and that, in this case,

<math>N(\mu)=\frac{1}{2\pi^2}\frac{mk_F}{\hbar^2},</math>

we find that

<math>\kappa=\sqrt{\frac{18\pi^3}{7\zeta(3)}}\frac{k_B T_c}{\sqrt{mc^2}\sqrt{e^2 k_F}}\left (\frac{c}{v_F}\right )^2.</math>

Note that we set <math>T=T_c</math> in the expression for <math>l_T</math>; this is because the GL theory is only valid just below the transition temperature.  We may also express this in terms of the Fermi energy,

<math>E_F=\frac{p_F^2}{2m}=\tfrac{1}{2}p_F v_F.</math>

Doing so, we obtain

<math>\kappa=\sqrt{\frac{18\pi^3}{7\zeta(3)}}\frac{k_B T_c}{\sqrt{mc^2}\sqrt{\alpha}\sqrt{2E_F}}\left (\frac{c}{v_F}\right )^{3/2}=8\times 10^{-6}\cdot\frac{T_c\,[\text{K}]}{\sqrt{E_F\,[\text{eV}]}}\left (\frac{c}{v_F}\right )^{3/2}.</math>

In a typical metal, <math>v_F\approx 10^5 - 10^6\,\tfrac{\text{m}}{\text{s}},</math> so

<math>\kappa=(0.04 - 1.3)\cdot\frac{T_c\,[\text{K}]}{\sqrt{E_F\,[\text{eV}]}}.</math>

=== A Simple Example - The Strongly Type-I Superconductor ===

As a simple demonstration of the solution of the GL equations, let us consider a strongly Type I (<math>\kappa\ll 1</math>) superconductor with a planar boundary between it and an insulator.  Let us set up our coordinate system so that the boundary is at <math>x=0</math>.

[[Image:SC_I_Boundary.PNG]]

We apply a magnetic field along the <math>z</math> axis,

<math>\vec{H}=H\hat{z}.</math>

We expect by symmetry that the total magnetic field <math>\vec{B}=B(x)\hat{z}</math>.  We will choose our gauge such that

<math>\vec{A}=A_y(x)\hat{y}.</math>

We also take the order parameter to depend only on <math>x</math>.  The first GL equation becomes

<math>-\frac{1}{\kappa^2}\frac{d^2\psi}{d\tilde{x}^2}+A_y^2\psi-\psi+\psi^3=0.</math>

Since we are taking <math>\kappa</math> to be small, the derivative term dominates, and we may therefore approximate this equation as

<math>\frac{d^2\psi}{d\tilde{x}^2}=0,</math>

so that <math>\psi(\tilde{x})=c_1+c_2\tilde{x}</math>.  Our boundary condition states that

<math>\left. \frac{1}{i\kappa}\frac{d\psi}{d\tilde{x}}\right |_{\tilde{x}=0}=0,</math>

so that <math>c_2=0</math>.  Since <math>\Psi=\Psi_0</math> in the bulk, we conclude that <math>\psi(\tilde{x})=1</math> for <math>\tilde{x}<0</math>.  Similarly, <math>\Psi=0</math> deep into the insulating region, so that <math>\psi(\tilde{x})=0</math> for <math>\tilde{x}>0</math>.

Now we consider the second equation.  In this case, it becomes, for <math>\tilde{x}<0</math>,

<math>-\tilde{A}=\tilde{\nabla}\times(\tilde{B}-\tilde{H})=\tilde{\nabla}\times\tilde{B},</math>

or

<math>-\tilde{B}=\tilde{\nabla}\times(\tilde{\nabla}\times\tilde{B}).</math>

The right-hand side is just

<math>\tilde{\nabla}\times(\tilde{\nabla}\times\tilde{B})=-\frac{d^2\tilde{B}}{d\tilde{x}^2},</math>

so that the equation is now

<math>\frac{d^2\tilde{B}}{d\tilde{x}^2}=\tilde{B}.</math>

The solution to the equation in simply <math>\tilde{B}(\tilde{x})=\tilde{B}_{+}e^{\tilde{x}}+\tilde{B}_{-}e^{-\tilde{x}}</math>, or, in terms of dimensional quantities,

<math>B(x)=B_{+}e^{x/\lambda}+B_{-}e^{-x/\lambda}.</math>

Since our superconductor is in the region <math>x<0</math>, we must take <math>B_{-}=0</math>.  Furthermore, the field must equal the applied field at <math>x=0</math>, so

<math>B(x)=He^{x/\lambda}.</math>

For <math>x>0</math>, the second GL equation becomes

<math>\frac{d^2\tilde{B}}{d\tilde{x}^2}=0.</math>

The solution, in terms of dimensional quantities, is <math>B(x)=B_0+B_1 x</math>.  We must set <math>B_1=0</math> so that the field does not increase indefinitely as we move away from the superconductor.  Since <math>\vec{B}=\vec{H}</math> in the normal state, we conclude that <math>B(x)=H</math> for <math>x>0</math>.

We have now shown why we called <math>\lambda</math> the penetration depth; it sets the length scale over which the magnetic field tends to zero inside the superconductor.  We have also illustrated the expulsion of applied magnetic fields from the interior of a Type I superconductor; this is known as the Meissner effect.

=== Thermodynamics of Type-I Superconductors in Magnetic Fields ===

In a bulk superconductor, surface effects are unimportant; for now, we will assume that the order parameter <math>\Psi</math> is constant everywhere in the superconductor and that magnetic fields are completely expelled.  In this case, the free energy per unit volume of the superconductor is

<math>f_s=\alpha(T-T_c)\Psi_0^2+\tfrac{1}{2}b\Psi_0^4=-\frac{\alpha^2}{2b}(T-T_c)^2.</math>

==Going below <math>T_c</math> with the Saddle Point Approximation==

So, previously all of this work has shown us the behavious of a superconducting system near <math>T_c</math> only. If we want to go into lower temperatures, we will have to make a careful saddle-point approximation, following Bardeen, Cooper, and Schrieffer (BCS). Once again, we can start from our [[#Pairing Hamiltonian and BCS instability|microscopic 'toy' Hamiltonian]], and gain useful information. 

Recall that the partition function can be written, 

<math> \mathbb{Z} = \int{d\Delta^* d\Delta} \left[ \int{D\psi^* D\psi \ e^{-S_{BCS} - S_{\Delta}}} \right] </math>

Where

<math> S_{\Delta} = \frac{-1}{g} \int_{0}^{\beta}{d\tau} \int{d^3r}\Delta^*(\vec{r},\tau) \Delta(\vec{r},\tau) </math>

and

<math> S_{BCS} = S_0 + S_{int} \;\; \text{where} \;\; S_0 = \int_{0}^{\beta}{d\tau}\int{d^3r}\left[ \psi_{\sigma}^*(\vec{r},\tau) \left( \frac{\partial}{\partial \tau} + \epsilon_p - \mu \right) \psi_{\sigma}(\vec{r},\tau) \right]  </math>
<math> \text{and} \;\; S_{int} = \int_{0}^{\beta}{d\tau}\int{d^3r} \left[\Delta^*(\vec{r},\tau)\psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau) + \Delta(\vec{r},\tau)\psi_{\downarrow}(\vec{r},\tau)\psi_{\uparrow}(\vec{r}.\tau) \right]</math>

Previously, we used a cumulant expansion around <math>T = T_c</math> to find the Ginzburg-Landau equations, along with the assumption that <math>\Delta</math> was small. Now, if we throw away this assumption, pray that fluctuations are small, and that there is only one saddle point, we will be able to successfully describe the superconducting state deep below <math>T_c</math> by demanding:

<math>\frac{\partial S_{eff}[\Delta]}{\partial \Delta*} = 0 </math>

===Self-Consistency Equation===

The solution of this functional derivative equation will give the value of <math>\Delta = \Delta_{sp}</math> at the saddle point (and also the self-consistency equation for this mean-field theory.)

To that end: 
 
: <math> S_{eff}[\Delta] = S_{\Delta} - \ln \left( \int{D\psi D\psi^*} \;\; e^{-S_{BCS}}\right) </math>

: <math> \frac{\partial S_{\Delta}}{\partial \Delta^*(\vec{r},\tau)} = \frac{-\Delta(\vec{r},\tau)}{g}</math>

: <math> \frac{\partial \left[ -\ln\left(\int{D\psi D\psi^*} \;\; e^{-S_{BCS}}\right) \right] }{\partial \Delta^*(\vec{r},\tau)} = \frac{-1}{\int{D\psi D\psi^*} \;\; e^{-S_{BCS}}} \int{D\psi D\psi*} \frac{\partial e^{-S_{BCS}}}{\partial \Delta^*} = \frac{\int{D\psi D\psi^*} \;\; e^{-S_{BCS}} \frac{\partial S_{BCS}}{\partial \Delta^*}}{\int{D\psi D\psi^*} \;\; e^{-S_{BCS}}} </math>

and, since 

: <math>\frac{\partial S_{BCS}}{\partial \Delta^*} = \psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau)</math>

We arrive at 

: <math>\frac{\partial \left[ -\ln\left(\int{D\psi D\psi^*} \;\; e^{-S_{BCS}}\right) \right] }{\partial \Delta^*(\vec{r},\tau)} = \frac{\int{D\psi D\psi^*} \;\; e^{-S_{BCS}} \psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau)}{\int{D\psi D\psi^*} \;\; e^{-S_{BCS}}} = \langle \psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau)\rangle_{BCS} </math>

So that, from the saddle-point condition, we find the Self-Consistency Equation: 

 <math> \frac{\Delta_{sp}(\vec{r},\tau)}{g} = \langle\psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau)\rangle_{BCS} </math>

This contains the same information as the Ginzburg-Landau equations, but also more, as we can now go far below <math>T_c</math>. Evaluation of this cam be done formally, but is difficult and not terribly enlightening. Instead, we will search for a solution in which <math\psi</math> is independent of <math>\vec{r}</math> and <math>\tau</math>, similar to our solution near <math>T_c<math>. In this case, we also expect that <math>\Delta</math> will be independent of <math>\vec{r}</math> and <math>\tau</math>. 

First, it is useful to write down the action in momentum and frequency space, before evaluating the correlator.

Start with <math> \psi(\vec{r},\tau) = \frac{1}{\beta}\sum_{\omega_n}\sum_{\vec{k}}\frac{1}{\sqrt{V}} e^{i\vec{k}\cdot\vec{r}}\psi_{k}(i\omega_n)</math>

For the three terms in <math>S_{BCS}</math>, we have:

<math> S_0 = \int_{0}^{\beta}{d\tau}\int{d^3r}\left[ \psi_{\sigma}^*(\vec{r},\tau) \left( \frac{\partial}{\partial \tau} + \epsilon - \mu \right) \psi_{\sigma}(\vec{r},\tau) \right] </math>
: <math> = \frac{1}{V}\sum_{\vec{k},\vec{k'}} \frac{1}{\beta^2}\sum_{\omega_n,\nu_m}\int_{0}^{\beta}{d\tau}\int{d^3r}e^{i\nu_{m}\tau}e^{-i\vec{k'}\cdot\vec{r}} \left(-i\omega_n + \epsilon_p - \mu \right) e^{-i\omega_n\tau}e^{i\vec{k}\cdot\vec{r}}\psi_{\sigma,\vec{k'}}^*(i\nu_m)\psi_{\sigma, \vec{k}}(i\omega_n)</math>

The integral over all space give <math>V \delta_{\vec{k},\vec{k'}}</math>, and the integral over imaginary time gives a factor of <math>\beta \delta(\omega_n - \nu_m)</math>, so that we find

<math>S_0 = \frac{1}{\beta}\sum_{\omega_n}\sum_{\vec{k}} \psi_{\sigma, \vec{k}}^*(i\omega_n)\left( -i\omega_n + \epsilon_p - \mu \right)\psi_{\sigma,\vec{k}}(i\omega_n) </math>

Now, for the 2nd (pairing) term: 

<math> \Delta^*\int_{0}^{\beta}{d\tau} \int{d^3r} \psi_{\downarrow}(\vec{r},\tau) \psi_{\uparrow}(\vec{r},\tau) = \Delta^* \frac{1}{V}\sum_{\vec{k},\vec{k'}} \frac{1}{\beta^2}\sum_{\omega_n,\nu_m}\int_{0}^{\beta}{d\tau}\int{d^3r} e^{-i(\nu_m + \omega_n)\tau) e^{i(\vec{k} + \vec{k'})\cdot\vec{r}} \psi_\downarrow,\vec{k'}}(i\nu_m) \psi_{\uparrow,\vec{k}}(i\omega_n) </math>

This time, the integral over real space gives <math>V \delta(\vec{k}+\vec{k'})</math>, and the integral over imaginary time gives a factor of <math>\beta \delta(\omega_n + \nu_m)</math>, so the 2nd term becomes: 

: <math> = \frac{\Delta^*}{\beta}\sum_{\vec{k}}\sum_{\omega_n}\psi_{\downarrow,\vec{k}}(i\omega_n)\psi_{\uparrow,-\vec{k}}(-i\omega_n)</math>

The hermitian conjugate of the above (the 3rd term in <math>S_{BCS}</math>) gives

:<math> \;\; \frac{\Delta}{\beta}\sum_{\vec{k}}\sum_{\omega_n}\psi_{\uparrow,-\vec{k}}^*(-i\omega_n)\psi_{\downarrow,\vec{k}}^*(i\omega_n)</math>

So that, for the entire <math>S_{BCS}</math>, we have found:

<math>S_{BCS} = \frac{1}{\beta}\sum_{\omega_n}\sum_{\vec{k}}\left[ \left( i\omega_n + \epsilon_{-\vec{k}} - \mu \right)\psi_{\uparrow, -\vec{k}}^*(-i\omega_n)\psi_{\uparrow,-\vec{k}}(-i\omega_n) + \left( -i\omega_n + \epsilon_{\vec{k}} - \mu \right)\psi_{\downarrow, \vec{k}}^*(i\omega_n)\psi_{\downarrow,\vec{k}}(i\omega_n) \right] </math>
::<math> + \frac{1}{\beta}\sum_{\omega_n}\sum_{\vec{k}}\left[ \Delta^*\psi_{\downarrow,\vec{k}}(i\omega_n)\psi_{\uparrow,-\vec{k}}(-i\omega_n) + \Delta\psi_{\uparrow,-\vec{k}}^*(-i\omega_n)\psi_{\downarrow,\vec{k}}^*(i\omega_n) \right] </math>

Now, to evaluate the Gaussian integrals in the correlator, it is extremely beneficial to write <math>S_{BCS}</math> like a matrix. To this end, we construct the so-called Nambu Spinors: 

<math> \Psi_{\vec{k}}(i\omega_n) = \begin{pmatrix} \psi_{-\vec{k},\uparrow}(-i\omega_n) \\ \psi^*_{\vec{k},\downarrow}(i\omega_n) \end{pmatrix} </math> 

<math> \Psi^*_{\vec{k}}(i\omega_n) = \begin{pmatrix} \psi^*_{-\vec{k},\uparrow}(-i\omega_n) & \psi^*_{\vec{k},\downarrow}(i\omega_n) \end{pmatrix} </math>

So now, we can write the BCS action as:

<math>
S_{BCS} = \frac{1}{\beta}\sum_{\omega_N}\sum_{\vec{k}}\Psi^*_{\vec{k}}(i\omega_n) 
\begin{pmatrix} 
i\omega_n + \epsilon_{-\vec{k}} - \mu & \Delta \\
\Delta^* & i\omega_n - \epsilon_{\vec{k}} + \mu 
\end{pmatrix}
\Psi_{\vec{k}}(i\omega_n)
</math>

Now, we can examine the correlator we found on the right-hand-side of the self-consistency equation: 

<math> \langle \psi_{\downarrow}(\vec{r},\tau)\psi_{\uparrow}(\vec{r},\tau)\rangle = \frac{1}{\beta^2 V}\sum_{\omega_n,\nu_m}e^{-i(\omega_n + \nu_M)\tau} \sum_{\vec{k},\vec{k'}} e^{i(\vec{k} + \vec{k'})\cdot\vec{r}}\langle \psi_{\vec{k},\downarrow}(i\omega_n) \psi_{\vec{k'},\uparrow}(i\nu_m) \rangle </math>
:<math> = \frac{1}{\beta^2 V}\sum_{\omega_n,\nu_m}e^{-i(\omega_n + \nu_M)\tau} \sum_{\vec{k},\vec{k'}} e^{i(\vec{k} + \vec{k'})\cdot\vec{r}} \langle \Psi^*_{\vec{k}}(i\omega_n)_2 \Psi_{-\vec{k'}}(-i\nu_m)_1 \rangle </math>

So, we need to write down a generic matrix for <math>\langle \Psi^*_{\vec{k}}(i\omega_n)_{\mu} \Psi_{-\vec{k'}}(-i\nu_m)_{\lambda} \rangle</math>, and take element <math> (2,1) </math> for our result. The gaussian integrals require that <math>\vec{k} = -\vec{k'}</math> and <math> i\omega_n = -i\nu_m </math> for convergence. 
Using our [[Phy5670#Gaussian_Integrals| technology from last semester]], we can show that:

<math>\langle \Psi^*_{\vec{k}}(i\omega_n)_{\mu} \Psi_{-\vec{k'}}(-i\nu_m)_{\lambda} \rangle = \delta(\vec{k} + \vec{k'})\beta \delta(\omega_n + \nu_m) \left[ 
\begin{pmatrix}
i\omega_n + \epsilon_{\vec{k}} - \mu & \Delta \\
\Delta^* & i\omega_n - \epsilon_{\vec{k}} + \mu 
\end{pmatrix}^{-1} \right]_{\mu\lambda}
</math>

<math>  \Psi^*_{\vec{k}}(i\omega_n)_2 \Psi_{-\vec{k'}}(-i\nu_m)_1 \rangle = \delta(\vec{k} + \vec{k'})\beta \delta(\omega_n + \nu_m)\frac{\Delta}{-\omega_n^2 - (\epsilon_{\vec{k}} - \mu)^2 - \Delta^* \Delta } </math>

Now, we can easily write down the correlator:

<math>\langle \psi_{\downarrow}(\vec{r},\tau)\psi_{\uparrow}(\vec{r},\tau)\rangle = \frac{1}{\beta}\sum_{\omega_n}\frac{1}{V}\sum_{\vec{k}}\frac{\Delta}{-\omega_n^2 - (\epsilon_{\vec{k}} - \mu)^2 - \Delta^* \Delta } </math>

This reshapes our consistency equation to:

 <math>\frac{\Delta}{g} = \frac{1}{\beta}\sum_{\omega_n}\frac{1}{V}\sum_{\vec{k}}\frac{\Delta}{-\omega_n^2 - (\epsilon_{\vec{k}} - \mu)^2 - \Delta^* \Delta } </math>

Which has a trivial, <math> \Delta = 0 </math> solution, and far more interesting solution where

<math>\frac{1}{g} = \frac{1}{\beta}\sum_{\omega_n}\frac{1}{V}\sum_{\vec{k}}\frac{1}{-\omega_n^2 - (\epsilon_{\vec{k}} - \mu)^2 - |\Delta|^2 } </math>

Note that this is only possible if <math> g < 0 </math>!

== Beyond saddle-point approximation, collective modes and response in the broken symmetry state==
Recall we can write our partition function as

<math>Z=\int D\Delta^*D\Delta e^{-S_{eff}[\Delta]}</math>

or

<math>Z=\int D(\Re e\Delta)D(\Im m\Delta)e^{-S_{eff}[\Re e\Delta,\Im m\Delta]}</math>

where

<math>S_{eff}[\Delta]=\frac{1}{|g|}\int_0^\beta d\tau\int d^Dr\left[(\Re e\Delta(r,\tau))^2+(\Im m\Delta(r,\tau))^2\right]-\ln\left[\int D\psi^*D\psi e^{-S_0-S_{int}}\right]</math>

and

<math>\begin{align} S_0&=\int_0^\beta d\tau\int d^Dr \psi^*(r,\tau)(\frac{\partial}{\partial\tau}+\epsilon_p-\mu)\psi(r,\tau)\\
S_{int}&=\int_0^\beta d\tau\int d^Dr(\Delta^*(r,\tau)\psi_\downarrow(r,\tau)\psi_{\uparrow}(r,\tau)+\Delta(r,\tau)\psi^*_\uparrow(r,\tau)\psi^*_{\downarrow}(r,\tau))
\end{align}</math>

We can rewrite the interaction term in the action as

<math>S_{int}=\int_0^\beta d\tau\int d^Dr\left[\Re e(\Delta(r,\tau))(\psi_{\downarrow}\psi_{\uparrow}+\psi_{\uparrow}^*\psi_{\downarrow}^*)(r,\tau)+i*\Im m(\Delta(r,\tau))(\psi_{\uparrow}^*\psi_{\downarrow}^*-\psi_{\downarrow}\psi_{\uparrow})(r,\tau)\right]</math>

Consider now functional derivatives of <math>S_{eff}[\Delta]</math>:

<math>
\frac{\delta S_{eff}[\Delta]}{\delta\Re e\Delta(r,\tau)}=\frac{2}{|g|}\Re e\Delta(r,\tau)-\frac{1}{\int D\psi^*D\psi e^{-S_0-S_{int}}} \int D\psi^*D\psi e^{-S_0-S_{int}}\left(-\frac{\delta S_{int}}{\delta \Re e\Delta(r,\tau)}\right)
</math>

The functional derivative of <math>S_{int}</math> with <math>\Re e\Delta(r,\tau)</math> is

<math>\frac{\delta S_{int}}{\delta \Re e\Delta(r,\tau)}=\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi_{\uparrow}^*(r,\tau)\psi_{\downarrow}^*(r,\tau)</math>

Therefore, we have

<math>\frac{\delta S_{eff}[\Delta]}{\delta \Re e\Delta(r,\tau)}=\frac{2}{|g|}\Re e\Delta(r,\tau)+\langle\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi_{\uparrow}^*(r,\tau)\psi_{\downarrow}^*(r,\tau)\rangle\ \ (1)</math>

Similarly

<math>\frac{\delta S_{eff}[\Delta]}{\delta \Im m\Delta(r,\tau)}=\frac{2}{|g|}\Im m\Delta(r,\tau)+i\langle\psi_{\uparrow}^*(r,\tau)\psi_{\downarrow}^*(r,\tau)-\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)\rangle\ \ (2)</math>

If we were to set the LHS of the above two equations to zero, we would obtain our self-consistency conditions. 

The strategy is to take

<math>S_{eff}[\Delta(r,\tau)]=S_{eff}\left[\Delta_{sp}+(\Delta(r,\tau)-\Delta_{sp})\right]</math>

where <math>\Delta_{sp}</math> solves (1) and (2) with LHS set to zero and expand in  powers of <math>\Delta(r,\tau)-\Delta_{sp}</math>. So

<math>\begin{align}
&S_{eff}[\Delta_{sp}+(\Delta(r,\tau)-\Delta_{sp})]\\
&=S_{eff}[\Delta_{sp}]+\int_0^\beta d\tau\int d^3r\Re e(\Delta(r,\tau)-\Delta_{sp})\frac{\delta S_{eff}}{\delta\Re e\Delta(r,\tau)}|_{\Delta_{sp}}+\int_0^\beta d\tau\int d^3r\Im m(\Delta(r,\tau)-\Delta_{sp})\frac{\delta S_{eff}}{\delta\Im m\Delta(r,\tau)}|_{\Delta_{sp}}\\
&+\frac{1}{2}\int_0^\beta d\tau\int_0^\beta d\tau'\int d^3r\int d^3r'\Re e(\Delta(r,\tau)-\Delta_{sp})\Re e(\Delta(r',\tau')-\Delta_{sp})\frac{\delta^2 S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}|_{\Delta_{sp}}\\
&+\frac{1}{2}\int_0^\beta d\tau\int_0^\beta d\tau'\int d^3r\int d^3r'\Im m(\Delta(r,\tau)-\Delta_{sp})\Im m(\Delta(r',\tau')-\Delta_{sp})\frac{\delta^2 S_{eff}}{\delta\Im m\Delta(r,\tau)\delta\Im m\Delta(r',\tau')}|_{\Delta_{sp}}\\
&+\int_0^\beta d\tau\int_0^\beta d\tau'\int d^3r\int d^3r'\Re e(\Delta(r,\tau)-\Delta_{sp})\Im m(\Delta(r',\tau')-\Delta_{sp})\frac{\delta^2 S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Im m\Delta(r',\tau')}|_{\Delta_{sp}}+...
\end{align}</math>

Since by definition of <math>\Delta_{sp}</math> we have

<math>\frac{\delta S_{eff}}{\delta\Re e\Delta(r,\tau)}|_{\Delta_{sp}}=\frac{\delta S_{eff}}{\delta\Im m\Delta(r,\tau)}|_{\Delta_{sp}}=0</math>

So only the 0th and 2nd order terms contribute. The 0th order term gave us the condensation energy, and the 2nd order term will give us information about collective modes (in the broken symmetry phase).

<math>\begin{align}
&\frac{\delta^2 S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}\\
&=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')+\frac{\delta}{\delta\Re e\Delta(r,\tau)}\left[\frac{1}{\int D\psi^*D\psi e^{-S_0-S_{int}}}\int D\psi^*D\psi e^{-S_0-S_{int}}\left(\frac{\delta S_{int}}{\delta \Re e\Delta(r',\tau')}\right)\right]\\
&=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')-\frac{\int D\psi^*D\psi e^{-S_0-S_{int}}\left(\frac{\delta S_{int}}{\delta \Re e\Delta(r,\tau)}\right)\left(\frac{\delta S_{int}}{\delta \Re e\Delta(r',\tau')}\right)}{\int D\psi^*D\psi e^{-S_0-S_{int}}}\\
&+\left(\frac{\int D\psi^*D\psi e^{-S_0-S_{int}}\left(\frac{\delta S_{int}}{\delta \Re e\Delta(r,\tau)}\right)}{\int D\psi^*D\psi e^{-S_0-S_{int}}}\right)\left(\frac{\int D\psi^*D\psi e^{-S_0-S_{int}}\left(\frac{\delta S_{int}}{\delta \Re e\Delta(r',\tau')}\right)}{\int D\psi^*D\psi e^{-S_0-S_{int}}}\right)
\end{align}</math>

So

<math>\begin{align}
&\frac{\delta^2 S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')\\
&-\langle\left(\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi^*_{\uparrow}(r,\tau)\psi^*_{\downarrow}(r,\tau)\right)\left(\psi_{\downarrow}(r',\tau')\psi_{\uparrow}(r',\tau')+\psi^*_{\uparrow}(r',\tau')\psi^*_{\downarrow}(r',\tau')\right)\rangle_{S_0+S_{int}}\\
&+\langle\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi^*_{\uparrow}(r,\tau)\psi^*_{\downarrow}(r,\tau)\rangle\langle\psi_{\downarrow}(r',\tau')\psi_{\uparrow}(r',\tau')+\psi^*_{\uparrow}(r',\tau')\psi^*_{\downarrow}(r',\tau')\rangle_{S_0+S_{int}}
\end{align}</math>

i.e.

<math>\begin{align}
&\frac{\delta^2 S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')\\
&-\langle\left(\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi^*_{\uparrow}(r,\tau)\psi^*_{\downarrow}(r,\tau)\right)\left(\psi_{\downarrow}(r',\tau')\psi_{\uparrow}(r',\tau')+\psi^*_{\uparrow}(r',\tau')\psi^*_{\downarrow}(r',\tau')\right)\rangle_{S_0+S_{int}}^{\mbox{con}}
\end{align}</math>

Similarly,

<math>\begin{align}
&\frac{\delta^2 S_{eff}}{\delta\Im m\Delta(r,\tau)\delta\Im m\Delta(r',\tau')}=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')\\
&+\langle\left(\psi^*_{\uparrow}(r,\tau)\psi^*_{\downarrow}(r,\tau)-\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)\right)\left(\psi^*_{\uparrow}(r',\tau')\psi^*_{\downarrow}(r',\tau')-\psi_{\downarrow}(r',\tau')\psi_{\uparrow}(r',\tau')\right)\rangle_{S_0+S_{int}}^{\mbox{con}}
\end{align}</math>

And

<math>
\frac{\delta^2 S_{eff}}{\delta\Im m\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}=-i\langle\left(\psi^*_{\uparrow}(r,\tau)\psi^*_{\downarrow}(r,\tau)-\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)\right)\left(\psi_{\downarrow}(r',\tau')\psi_{\uparrow}(r',\tau')+\psi^*_{\uparrow}(r',\tau')\psi^*_{\downarrow}(r',\tau')\right)\rangle_{S_0+S_{int}}^{\mbox{con}}
</math>

If we evaluate these functional derivatives at the saddle point <math>\Delta_{sp}</math> we have

<math>\begin{align}
S_0+S_{int}&\rightarrow\int_0^\beta d\tau\int d^3r\sum_{\sigma=\uparrow,\downarrow}\psi^*_\sigma(r,\tau)\left(\frac{\partial}{\partial\tau}+\epsilon_p-\mu\right)\psi_{\sigma}(r,\tau)\\
&+\Delta_{sp}\int_0^\beta d\tau\int d^3r\left(\psi_{\downarrow}(r,\tau)\psi_{\uparrow}(r,\tau)+\psi_{\uparrow}^*(r,\tau)\psi_{\downarrow}^*(r,\tau)\right)
\end{align}</math>

where we take the saddle point solution <math>\Delta_{sp}</math> to be purely real. Arranging the Grassman fields into the Nambu spinor we have:

<math>(\psi^*_{\uparrow}\ \psi_\downarrow)\left(\begin{align}&\frac{\partial}{\partial\tau}+\epsilon_p-\mu&\Delta_{sp}\\ &\Delta_{sp}&\frac{\partial}{\partial\tau}-\epsilon_p+\mu\end{align}\right) \left(\begin{align}&\psi_\uparrow\\&\psi^*_\downarrow\end{align}\right)</math>

where the Nambu spinor is defined as

<math>\Psi=\left(\begin{align}&\psi_\uparrow\\&\psi^*_\downarrow\end{align}\right)</math>

So

<math>\begin{align}
&\psi^*_\uparrow\psi^*_\downarrow+\psi_\downarrow\psi_\uparrow=\Psi^*\left(\begin{align}&0&1\\&1&0\end{align}\right)\Psi =\Psi^*\sigma_x\Psi\\
&\psi^*_\uparrow\psi^*_\downarrow-\psi_\downarrow\psi_\uparrow=\Psi^*\left(\begin{align}&0&1\\-&1&0\end{align}\right)\Psi =\Psi^*i\sigma_y\Psi
\end{align}</math>

<math>\begin{align}
\frac{\delta^2S_{eff}}{\delta\Re e\Delta(r,\tau)\delta\Re e\Delta(r',\tau')}|_{\Delta_{sp}}&=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')-\langle\Psi^*(r,\tau)\sigma_x\Psi(r,\tau)\Psi^*(r',\tau')\sigma_x\Psi(r',\tau')\rangle^{\mbox{con}}\\
&=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')+\mbox{Tr}\left[\sigma_xG(r-r',\tau-\tau')\sigma_xG(r'-r,\tau'-\tau)\right]
\end{align}</math>

where the Green's functions are 2*2 matrices, 

<math>G^{-1}_k(i\omega_n)=\left(\begin{align}&-i\omega_n+\epsilon_k-\mu&\Delta_{sp}\\&\Delta_{sp}&-i\omega_n-\epsilon_k+\mu\end{align}\right)</math> 

Notice that this is a function of <math>\tau-\tau'</math> and <math>r-r'</math>. Let's call it <math>\Pi_{++}(r-r',\tau-\tau')</math>. This will give rise to

<math>\begin{align}
&\frac{1}{2}\int_0^\beta d\tau\int_0^\beta d\tau'\int d^3r\int d^3r'\Re e(\Delta(r,\tau)-\Delta_{sp})\Pi_{++}(r-r',\tau-\tau')\Re e(\Delta(r',\tau')-\Delta_{sp})\\
&=\frac{1}{2}\int_0^\beta d\tau\int_0^\beta d\tau'\int d^3r\int d^3r'\delta\Delta_+(r,\tau)\Pi_{++}(r-r',\tau-\tau')\delta\Delta_+(r',\tau')
\end{align}</math>

Fourier transforming we find the contribution to <math>S_{eff}</math>

<math>\frac{1}{2}\frac{1}{\beta}\sum_{\Omega_n}\sum_{q}\delta\Delta_+(-q,-i\Omega_n)\Pi_{++}(q,i\Omega_n)\delta\Delta_+(q,i\Omega_n)</math>

Similarly

<math>\begin{align}
\Pi_{--}(r-r',\tau-\tau')&=\frac{2}{|g|}\delta(\tau-\tau')\delta(r-r')+\mbox{Tr}\left[\sigma_yG(r-r',\tau-\tau')\sigma_yG(r'-r,\tau'-\tau)\right]\\
\Pi_{-+}(r-r',\tau-\tau')&=\langle\Psi^*\sigma_y\Psi(r,\tau)\Psi^*\sigma_x\Psi(r',\tau')\rangle^{\mbox{con}}=-\mbox{Tr}\left[\sigma_yG(r-r',\tau-\tau')\sigma_xG(r'-r,\tau'-\tau)\right]
\end{align}</math>

Then

<math>\begin{align}
&S_{eff}[\Delta]\approx S_{eff}[\Delta_{sp}]+\\
&\frac{1}{2}\frac{1}{\beta}\sum_{\Omega_n}\sum_q\left[\delta\Delta_+(-q,-i\Omega_n)\Pi_{++}(q,i\Omega_n)\delta\Delta_+(q,i\Omega_n)+\delta\Delta_-(-q,-i\Omega_n)\Pi_{--}(q,i\Omega_n)\delta\Delta_-(q,i\Omega_n)\right.\\
&\left.+2\delta\Delta_-(-q,-i\Omega_n)\Pi_{-+}(q,i\Omega_n)\delta\Delta_+(q,i\Omega_n)\right]+...
\end{align}</math>

To proceed with the evaluation of <math>\Pi</math>'s, note that

<math>\begin{align}
&G^{-1}_k(i\omega_n)=-i\omega_n\mathbb{I}+(\epsilon_k-\mu)\sigma_3+\Delta_{sp}\sigma_1\\
\Rightarrow&G_k(i\omega_n)=\frac{i\omega_n\mathbb{I}+(\epsilon_k-\mu)\sigma_3+\Delta_{sp}\sigma_1}{\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2}
\end{align}</math>

where <math>\omega_n=\frac{(2n+1)\pi}{\beta}</math>. After Fourier transform,

<math>\begin{align}
\Pi_{++}(q,i\Omega_n&=\frac{2}{|g|}+\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}\mbox{Tr}\left[\sigma_1G_{k+q}(i\omega_n+i\Omega_n)\sigma_1G_k(i\omega_n)\right]\\
\Pi_{--}(q,i\Omega_n&=\frac{2}{|g|}+\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}\mbox{Tr}\left[\sigma_2G_{k+q}(i\omega_n+i\Omega_n)\sigma_2G_k(i\omega_n)\right]\\
\Pi_{-+}(q,i\Omega_n&=-\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}\mbox{Tr}\left[\sigma_2G_{k+q}(i\omega_n+i\Omega_n)\sigma_1G_k(i\omega_n)\right]\\
\end{align}</math>

Consider now

<math>\begin{align}
&\mbox{Tr}\left[\sigma_{1,2}\frac{(i\omega_n+i\Omega_n)\mathbb{I}+(\epsilon_{k+q}-\mu)\sigma_3+\Delta_{sp}\sigma_1}{(\omega_n+\Omega_n)^2+(\epsilon_{k+q}-\mu)^2+\Delta_{sp}^2}\sigma_{1,2}\frac{i\omega_n\mathbb{I}+(\epsilon_k-\mu)\sigma_3+\Delta_{sp}\sigma_1}{\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2}\right]\\
=&2\frac{(i\omega_n+i\Omega_n)i\omega_n-(\epsilon_{k+q}-\mu)(\epsilon_k-\mu)\pm\Delta^2_{sp}}{\left[(\omega_n+\Omega_n)^2+(\epsilon_{k+q}-\mu)^2+\Delta_{sp}^2\right]\left[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2\right]}
\end{align}</math>

and

<math>\begin{align}
&\mbox{Tr}\left[\sigma_{2}\frac{(i\omega_n+i\Omega_n)\mathbb{I}+(\epsilon_{k+q}-\mu)\sigma_3+\Delta_{sp}\sigma_1}{(\omega_n+\Omega_n)^2+(\epsilon_{k+q}-\mu)^2+\Delta_{sp}^2}\sigma_{1}\frac{i\omega_n\mathbb{I}+(\epsilon_k-\mu)\sigma_3+\Delta_{sp}\sigma_1}{\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2}\right]\\
=&\frac{(i\omega_n+i\Omega_n)(\epsilon_k-\mu)\mbox{Tr}[\sigma_2\sigma_1\sigma_3]+i\omega_n(\epsilon_{k+q}-\mu)\mbox{Tr}[\sigma_2\sigma_3\sigma_1]}{\left[(\omega_n+\Omega_n)^2+(\epsilon_{k+q}-\mu)^2+\Delta_{sp}^2\right]\left[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2\right]}\\
=&\frac{(i\omega_n+i\Omega_n)(\epsilon_k-\mu)(-2i)+i\omega_n(\epsilon_{k+q}-\mu)2i}{\left[(\omega_n+\Omega_n)^2+(\epsilon_{k+q}-\mu)^2+\Delta_{sp}^2\right]\left[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2\right]}
\end{align}</math>

Now, note that at <math>q=0,\Omega_n=0</math>,

<math>\Pi_{-+}(0,0)=0</math>

because the numerator vanishes for all <math>k</math> and <math>\omega_n</math>. Also note that by self-consistency condition

<math>\begin{align}
\Pi_{--}(0,0)&=\frac{2}{|g|}-\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}2\frac{1}{[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2]}\\
&=\frac{2}{|g|}+2\int\frac{d^3k}{(2\pi)^3}\frac{1}{\beta}\sum_{\omega_n}\frac{1}{(i\omega_n-E_k)(i\omega_n+E_k)}\\
&=\frac{2}{|g|}+2\int\frac{d^3k}{(2\pi)^3}\left(\frac{n_F(E_k)}{2E_k}-\frac{n_F(-E_k)}{2E_k}\right)\\
&=2\left(\frac{1}{|g|}-\int\frac{d^3k}{(2\pi)^3}\frac{1-2n_F(E_k)}{2E_k}\right)\\
&=0
\end{align}</math>

However,

<math>\begin{align}
\Pi_{++}(0,0)&=\frac{2}{|g|}+\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}2\frac{-\omega_n^2-(\epsilon_k-\mu)^2+\Delta_{sp}^2}{[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2]^2}\\
&=4\Delta_{sp}^2\frac{1}{\beta}\sum_{\omega_n}\int\frac{d^3k}{(2\pi)^3}\frac{1}{[\omega_n^2+(\epsilon_k-\mu)^2+\Delta_{sp}^2]^2}\\
&=4\Delta_{sp}^2N_0\frac{1}{\beta}\sum_{\omega_n}\int_{-\infty}^\infty\frac{d\xi}{[\xi^2+\omega_n^2+\Delta_{sp}^2]^2}\\
&=4\Delta_{sp}^2N_0\frac{1}{\beta}\sum_{\omega_n}\frac{\pi}{2}\frac{1}{(\omega_n^2+\Delta_{sp}^2)^{\frac{3}{2}}}\\
&=2\pi N_0\Delta_{sp}^2\frac{1}{\pi^3T^2}\sum_{n=-\infty}^\infty\frac{1}{[(2n+1)^2+(\frac{\Delta_{sp}}{\pi T})^2]^{\frac{3}{2}}}
\end{align}</math>

This sum is slowly convergent. To evaluate it efficiently, we note that

<math>\frac{1}{a^{\frac{3}{2}}}=\frac{2}{\sqrt{\pi}}\int_0^\infty d\lambda\sqrt{\lambda}e^{-\lambda a}\ \ \mbox{for }a>0</math>

So

<math>\begin{align}
\sum_{n=-\infty}^\infty\frac{1}{[(2n+1)^2+(\frac{\Delta_{sp}}{\pi T})^2]^\frac{3}{2}}&=\frac{2}{\sqrt{\pi}} \int_0^\infty d\lambda\sqrt{\lambda}e^{-\lambda(\frac{\Delta_{sp}}{\pi T})^2}\left(\sum_{n=-\infty}^\infty e^{-\lambda(2n+1)^2}\right)\\
&=\frac{2}{\sqrt{\pi}} \int_0^\infty d\lambda\sqrt{\lambda}e^{-\lambda(\frac{\Delta_{sp}}{\pi T})^2}\Theta_2(0,e^{-4\lambda})\\
&=\frac{1}{4\sqrt{\pi}} \int_0^\infty d\xi\sqrt{\xi}e^{-\xi(\frac{\Delta_{sp}}{2\pi T})^2}\Theta_2(0,e^{-\xi})
\end{align}</math>

where <math>\Theta_2</math> is Jacobi elliptic theta function. Then

<math>\begin{align}
\Pi_{++}(0,0)&=2\pi N_0\Delta^2_{sp}\frac{1}{\pi^3T^2}\sum_{n=-\infty}^\infty\frac{1}{[(2n+1)^2+(\frac{\Delta_{sp}}{\pi T})^2]^\frac{3}{2}}\\
&=2N_0\Phi(\frac{\Delta_{sp}}{2\pi T})
\end{align}</math>

where

<math>\Phi(x)=x^2\frac{1}{\sqrt{\pi}}\int_0^\infty d\xi\sqrt{\xi}e^{-\xi x^2}\Theta_2(0,e^{-\xi})</math>

Note that

<math>\Pi_{++}(0,0)\rightarrow 2N_0\frac{7\zeta(3)}{4\pi^2T^2}\Delta^2_{sp}</math>

as <math>\Delta_{sp}\rightarrow 0</math>, i.e. as <math>T\rightarrow T_C</math>.

That is precisely the curvature of the new minimum in the Ginzburg-Landau free energy we found before. So, at 2nd order our effective action corresponds to the action for two real free bosons, <math>\delta\Delta_+</math> and <math>\delta\Delta_-</math>. At <math>q=0, i\Omega_n=0</math>, <math>\delta\Delta_+</math> mode is gapped (massive), but <math>\delta\Delta_-</math> mode is not gapped (massless). Physically, <math>\delta\Delta_+</math> corresponds to the fluctuations of the order parameter amplitude (because we chose <math>\Delta_{sp}</math> to be real). 

Recall our discussion from many-body course:

Amplitude fluctuations are not hydrodynamic modes since they do not correspond to either conserved or to broken symmetry variable. We have to extend our approach to higher order in <math>\delta\Delta_\pm</math> to describe its (rapid) decay.

<math>\delta\Delta_-</math> corresponds to fluctuations along the direction of the minimum of the double well potential, where there is no barrier. It corresponds to (part of) a "phase" mode.

To determine the kinematics (of our collective modes) we need to expand <math>\Pi_{\mu\nu}(q, i\Omega_n)</math> in powers of <math>\vec{q}</math> and <math>\Omega_n</math>. Our small expansion parameters are

<math>\Omega_n \ll \Delta_{sp}</math> and <math>q \ll \frac{\Delta_{sp}}{\hbar v_F}</math>

(Obviously <math>q \ll k_F</math>)

Start with <math>\Pi_{++}</math> :
What we need is

<math>\Pi_{++}(q, i\Omega_n) - \Pi_{++}(0, 0)</math>

<math>= \frac{2}{\beta}\sum_{\omega_n}{\int{\frac{d^3k}{(2\pi)^3}\left[\frac{1}{\omega_n^2 + E_k^2}\right]}}</math>