Files
2021-04-05 09:47:49 -05:00

317 lines
18 KiB
TeX

\documentclass{article}
%other packages
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{physics}
\usepackage[
style=phys, articletitle=false, biblabel=brackets, chaptertitle=false, pageranges=false, url=true
]{biblatex}
\usepackage{graphicx}
\usepackage{todonotes}
\usepackage{siunitx}
\usepackage[compat=1.0.0]{tikz-feynman}
\usepackage{cleveref}
\title{Notes on the free energy calculation for Owen-Scalapino}
\date{2020-12-20}
\addbibresource{./bibliography.bib}
\graphicspath{{./figures/}}
\newcommand{\pf}{p_{\mathrm{F}}}
\newcommand{\vf}{v_{\mathrm{F}}}
\newcommand{\epsF}{\epsilon_{\mathrm{F}}}
\newcommand{\debye}{\omega_{\mathrm{D}}}
\newcommand{\corr}{\mu^{\ast}}
\newcommand{\dof}{\rho}
\DeclareMathOperator{\sgn}{sgn}
\begin{document}
\maketitle
We are interested in the free energy calculation because we want to be able to correctly limit the range of superconducting behaviour for Owen-Scalapino.
As seen in~\cite{OwenScalapino}, the superconducting state is restricted in its temperature range, with the traditional $T_c$ lowered and with an additional lower critical temperature below which the material returns to the normal state.
There are some subtleties with how we need to handle this calculation in the OS case, which we will discuss.
\section{Superconducting Free energy} \label{sec:scfreeen}
We begin with the standard BCS free energy expression:
\begin{equation}
F_s = 2 \sum_k \epsilon_k \left(f_k - 2 f_k v_k^2 + v_k^2\right) - \sum_{k, l} V_{kl} u_k v_k u_l v_l (1 - 2 f_k) (1 - 2 f_l) - TS \label{eq:bcs:fs}
\end{equation}
For the distribution function $f_k$, we have
\begin{equation}
f_k = \frac{1}{1 + \exp\frac{E_k - \corr}{T}},
\end{equation}
where $E_k^2 = \epsilon_k^2 + \Delta^2$.
Our coherence factors are unchanged by the additional $\corr$ constraint in our case, apart from the modified value to $\Delta_k$:
\begin{align}
\abs{u_k}^2 &= \frac12 \left(1 + \frac{\epsilon_k}{E_k} \right) \\
\abs{v_k}^2 &= \frac12 \left(1 - \frac{\epsilon_k}{E_k} \right).
\end{align}
These can be obtained most simply from the expression
\begin{equation}
\tan 2 \theta_k = - \frac{\Delta_k}{\epsilon_k}
\end{equation}
from the December 4 notes, as verification that the coherence factor does not change.
We'll also use
\begin{equation}
u_k^\ast v_k = \frac12 \frac{\Delta_k}{E_k}
\end{equation}
In \cref{eq:bcs:fs}, incidentally, there are some suppressed magnitude signs, as only one of $u_k$ and $v_k$ can always be made real.
This will be irrelevant here though, I believe, as at best the interaction term in the energy will have the phase cancel out.
We can begin reducing this for code by writing $F_s$ as $A - B - C$ with the three terms in \cref{eq:bcs:fs},
and simplifying each in turn:
\begin{align}
A &= 2 \sum_k \epsilon_k \left(f_k - 2 f_k v_k^2 + v_k^2 \right) \\
&= 2 \sum_k \epsilon_k \left(f_k - f_k\left(1 - \frac{\epsilon_k}{E_k} \right) + \frac12 \left(1 - \frac{\epsilon_k}{E_k} \right) \right) \\
&= 2 \sum_k \epsilon_k \left(f_k\frac{\epsilon_k}{E_k} + \frac12 - \frac12 \frac{\epsilon_k}{E_k} \right) \\
&= \sum_k \epsilon_k \left(1 + 2 f_k\frac{\epsilon_k}{E_k} - \frac{\epsilon_k}{E_k} \right) \\
\end{align}
Next,
\begin{align}
B &= \sum_{k, l} V_{kl} u_k v_k u_l v_l (1 - 2 f_k) (1 - 2 f_l) \\
&= \sum_{k, l} V_{kl} \frac12 \frac{\Delta_k}{E_k} \frac12 \frac{\Delta_l}{E_l} (1 - 2 f_k) (1 - 2 f_l) \\
&= \frac14 \sum_{k, l} V_{kl} \frac{\Delta_k}{E_k} \frac{\Delta_l}{E_l} (1 - 2 f_k) (1 - 2 f_l) \\
\end{align}
Finally, the entropy term:
\begin{align}
C &= TS \\
&= -2 T \sum_k f_k \log f_k + (1 - f_k) \log (1 - f_k)
\end{align}
Next, we make our assumptions about the momentum sum.
If we assume that our density of states is constant over the relevant interval, we can replace the integral sums with $\sum_k \rightarrow 2 \int_{0}^{\debye} \dof \dd{\epsilon_k}$.
The $2$ factor here comes from assuming that all of our sums are symmetric in $\epsilon_k$, which is also why our integral begins from $0$.
I don't know if this assumption is correct, but I believe it is (and it handles the intuition that $v_k$ represents our electronlike quasiparticles, which are relevant for $\epsilon_k > \epsF$)
This assumes that we don't consider quasiparticles with energy above $\debye$ as participating in this interaction, and that $\debye \ll \epsF$, which are both reasonable assumptions for this problem.
This means that our term $A$ is then
\begin{align}
A &= \int_{-\debye}^{\debye} \dof \dd{\epsilon_k} \epsilon_k \left(1 + 2 f_k\frac{\epsilon_k}{E_k} - \frac{\epsilon_k}{E_k} \right) \\
&= \dof \int_{-\debye}^{\debye} \dd{\epsilon_k} \epsilon_k \left(1 + 2 f_k\frac{\epsilon_k}{E_k} - \frac{\epsilon_k}{E_k} \right) \\
&= \dof \int_{-\debye}^{\debye} \dd{\epsilon_k} \epsilon_k \left(1 + 2 \left( \frac{1}{1 + \exp\frac{E_k - \corr}{T}} \right)\frac{\epsilon_k}{E_k} - \frac{\epsilon_k}{E_k} \right)
\end{align}
We drop the first term, as it is odd.
The next two terms are even, as they only depend on $\epsilon_k^2$.
\begin{align}
A &= 2 \dof \int_{0}^{\debye} \dd{\epsilon_k} 2 \epsilon_k \left( \frac{1}{1 + \exp\frac{E_k - \corr}{T}} \right)\frac{\epsilon_k}{E_k} - \frac{\epsilon_k^2}{E_k} \\
&= 2 \dof \int_{0}^{\debye} \dd{\epsilon_k} \frac{\epsilon_k^2}{E_k} \left( \frac{2}{1 + \exp\frac{E_k - \corr}{T}} - 1\right) \\;
&= 2 \dof \int_{0}^{\debye} \dd{\epsilon_k} \frac{\epsilon_k^2}{E_k} \frac{2}{1 + \exp\frac{E_k - \corr}{T}} - 2 \dof \int_0^{\debye} \dd{\epsilon_k} \frac{\epsilon_k^2}{\sqrt{k^2 + \Delta^2}} \\
&= 2 \dof \int_{0}^{\debye} \dd{\epsilon_k} \frac{\epsilon_k^2}{E_k} \frac{2}{1 + \exp\frac{E_k - \corr}{T}} - \dof \left( \debye \sqrt{\Delta^2 + \debye^2} + \Delta^2 \log \frac{\Delta}{\debye + \sqrt{\Delta^2 + \debye^2}} \right)
\end{align}
Similarly, we can assume that $V_{kl} \rightarrow - V$ and $\Delta_k \rightarrow \Delta$, as part of our standard weak coupling assumption set, giving us
\begin{align}
B &= V \int_{0}^{\debye} \dof \dd{\epsilon_k} \int_{0}^{\debye} \dof \dd{\epsilon_l} \frac{\Delta_k}{E_k} \frac{\Delta_l}{E_l} (1 - 2 f_k) (1 - 2 f_l) \\
&= V \left(\dof \int_{0}^{\debye} \dd{\epsilon_k} \frac{\Delta_k}{E_k} (1 - \frac{2}{1 + \exp\frac{E_k - \corr}{T}}) \right)^2 \\
&= V \left(\dof \Delta \int_{0}^{\debye} \frac{\dd{\epsilon_k}}{E_k} (1 - \frac{2}{1 + \exp\frac{E_k - \corr}{T}}) \right)^2.
\end{align}
Lastly,
\begin{align}
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \log (1 - f_k) \\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \log (1 - \frac{1}{1 + \exp\frac{E_k - \corr}{T}}) \\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \log (\frac{\exp\frac{E_k - \corr}{T}}{1 + \exp\frac{E_k - \corr}{T}}) \\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \left(\log(\exp\frac{E_k - \corr}{T}) - \log({1 + \exp\frac{E_k - \corr}{T}})\right) \\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \left(\frac{E_k - \corr}{T} - \log({1 + \exp\frac{E_k - \corr}{T}})\right)
\end{align}
We can simplify this slightly further:
\begin{align}
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} f_k \log f_k + (1 - f_k) \left(\frac{E_k - \corr}{T} + \log f_k\right) \\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} \log f_k + (1 - f_k) \frac{E_k - \corr}{T}\\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} -\log({1 + \exp\frac{E_k - \corr}{T}}) + (1 - f_k) \frac{E_k - \corr}{T}\\
C &= - 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} -\log({\exp\frac{-E_k + \corr}{T}} + 1) - \frac{E_k - \corr}{T} + (1 - f_k) \frac{E_k - \corr}{T}\\
C &= 4 T \int_{0}^{\debye} \dof \dd{\epsilon_k} \log({\exp\frac{-E_k + \corr}{T}} + 1) + f_k \frac{E_k - \corr}{T}\\
\end{align}
\section{Owen-Scalapino Gap Calculation}
We are interested in being able to calculate the gap for different quasiparticle densities $n$, as shown in Owen-Scalapino's figure 1b in~\cite{OwenScalapino}.
To do this, we use the modified gap equation:
\begin{equation}
\left[ N(0) V \right]^{-1} = \int_{- \omega_D}^{\omega_D} \frac{\dd{\epsilon_k}}{\sqrt{\Delta^2 + \epsilon_k^2}} \tanh{\frac{\sqrt{\Delta^2 + \epsilon_k^2} - \corr}{2 T}}. \label{eq:gap}
\end{equation}
OS's equation 7 is the definition of a parameter $n$ measuring the excess quasiparticle density in units of $4 N(0) \Delta_0$:
\begin{equation}
n = \frac{1}{\Delta_0} \int_0^\infty \dd{\epsilon} \left( \frac{1}{1 + \exp(\frac{\sqrt{\Delta^2 + \epsilon_k^2} - \corr}{T})} - \frac{1}{1 + \exp(\frac{\sqrt{\Delta^2 + \epsilon_k^2}}{T})} \right). \label{eq:n}
\end{equation}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{osGapVsTReproduction}
\caption{Reproduction of OS graph 1b, with parameters described in the text.} \label{fig:osRepro1b}
\end{figure}
First, to find $\Delta_0$, we solve the standard gap equation with $T=0$ and $\corr = 0$.
For \cref{fig:osRepro1b}, this is done for $N(0) V = 0.2$ and $\debye = 100$.
The values chosen have no special significance, but keep our energy scale at order 1, with $\Delta_0 = 1.34765$ and $T_c = 0.764083$.
Then, we numerically solve the pair of equations \cref{eq:gap} and \cref{eq:n} for the free parameters $\Delta$ and $\corr$ for our chosen value of $T$, keeping $N(0)$, $V$ and $\debye$ fixed.
Playing with the parameters slightly, I had success over the important part of the solution space starting at $(\Delta, \corr) = (\frac{\debye}{\sinh(\left(N(0) V\right)^{-1})}, n^2 \Delta_0)$.
A more sophisticated analysis could most likely be done to fix the starting space more accurately and avoid the messy behaviour in the region where the free energy is lower for the normal state.
However, as we discussed, for sufficiently small $n$ (n \lessapprox 0.1) this is unlikely to matter, as for small $T \rightarrow 0$ there is no return to the normal state, and for large $T$ the critical temperature is close to $T_c$ for $n = 0$ (OS find that the critical temperature only drops to around $0.9 T_{c, 0}$ for $n = 0.02$, judging by their plot).
The result of solving this pair of equations is the current correction to the chemical potential $\corr$ and the gap $\Delta$, which can be used in the Nam expression.
The numerical equations are not overly slow;
for the starting points given in the energetically favourable region the descent to the solution seems easy for Mathematica to handle.
\section{Nam usage}
Because the earlier result matches OS's curve, we are more confident that we can use this as a calculation for $\Delta$ in the Nam code.
However, the Nam code takes in the following parameters:
\begin{enumerate}
\item $\omega$ in \si{\per\s}
\item The Plasma frequency $\omega_p$, in \si{\per\s}
\item The impurity collision frequency $\tau$, in \si{\per\s}
\item $\vf$, in \si{\m\per\s}
\item A critical temperature $T_c$, in \si{\per\s}
\item The dipole moment $d$, in $\si{\coulomb \m}$
\end{enumerate}
The point here is the critical temperature $T_c$, which is essentially only used for calculating the gap $\Delta$ via $\Delta = 3.06 \sqrt{T_c \left( T_c - T\right)}$.
Because this is no longer the expression of interest, replacing this with the OS gap package should be enough to fulfill that role.
The OS gap package takes in its own parameters:
\begin{enumerate}
\item $T$, $\corr$ $\debye$ and the interaction $V$ in \si{\per\s}
\item $N(0)$ the density of states in \si{\s}.
\end{enumerate}
We can keep the same expressions as earlier with the same reduced dimensionless units, keeping the variable $\mu$:
\begin{align}
\xi &= \frac{\omega}{\Delta} \\
\xi' &= \frac{\omega'}{\Delta} \\
\nu &= \frac{1}{\tau \Delta} \\
\kappa &= \frac{q \vf}{\Delta} \\
t &= \frac{T}{\Delta} \\
\sigma_0 &= \frac{n e^2}{m \Delta}
\mu \leftarrow \frac{\mu}{\Delta}
\end{align}
And then our expression, with the changes becomes:
\begin{equation}
\sigma(\kappa, \xi) = -i \frac{3 \sigma_0}{4} \frac{1}{\xi}\left[\int_{1 - \xi}^{1}\dd{\xi} \tanh(\frac{\xi + \xi' - \mu}{2 t}) I_1 + \int_{1}^{\infty} \dd{\xi'} \left( \tanh(\frac{\xi + \xi' - \mu}{2t}) I_1 - \tanh(\frac{\xi' - \mu}{2t})I_2 \right) \right] \label{eq:nam}
\end{equation}
with
\begin{align}
I_1 &= F(\kappa, \Re[\sqrt{(\xi + \xi')^2 - 1} - \sqrt{\xi'^2 - 1}]) (g + 1) \nonumber\\
&\quad + F(\kappa, \Re[-\sqrt{(\xi + \xi')^2 - 1} - \sqrt{\xi'^2 - 1}]) (g - 1) \\
I_2 &= F(\kappa, \Re[\sqrt{(\xi + \xi')^2 - 1} - \sqrt{\xi'^2 - 1}]) (g + 1) \nonumber\\
&\quad + F(\kappa, \Re[\sqrt{(\xi + \xi')^2 - 1} + \sqrt{\xi'^2 - 1}]) (g - 1) \\
F(\kappa, E) &= \frac{1}{\kappa} \left[2 S(E) + (1 - S(E)^2)\ln(\frac{S(E) + 1}{S(E) - 1})\right] \\
S(\kappa, E) &= \frac{1}{\kappa} \left(E - i \left(\Im[\sqrt{(\xi + \xi')^2 - 1} + \sqrt{\xi'^2 - 1}] + 2 \nu \right) \right) \\
g &= \frac{\xi' \left( \xi + \xi'\right) + 1}{\sqrt{\xi'^2 - 1}\sqrt{(\xi + \xi')^2 - 1}}
\end{align}
\section{Mathematica implementation}
We want to dimensionally reduce these when we implement the code, by writing $\corr$, $\Delta$, $\omega_D$ and $T$ in units of $\Delta_0$:
\begin{align}
o_D &= \frac{\omega_D}{\Delta_0} \\
m &= \frac{\corr}{\Delta_0} \\
d &= \frac{\Delta}{\Delta_0} \\
t &= \frac{T}{\Delta_0}
\end{align}
giving us
\begin{equation}
\left[ N(0) V \right]^{-1} = \int_{- o_D}^{o_D} \frac{\dd{\varepsilon}}{\sqrt{d^2 + \varepsilon^2}} \tanh{\frac{\sqrt{d^2 + \varepsilon^2} - m}{2 t}}
\end{equation}
\begin{equation}
n = \int_0^\infty \dd{\epsilon} \left( \frac{1}{1 + \exp(\frac{\sqrt{d + \varepsilon^2} - m}{t})} - \frac{1}{1 + \exp(\frac{\sqrt{d^2 + \varepsilon^2}}{t})} \right)
\end{equation}
We can look at the Nam conductivity now:
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{realpartofconductivity}
\caption{$\frac{\Re[\sigma_{SC}]}{\sigma_N}$} \label{fig:real}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{imagpartofconductivity}
\caption{$\frac{\Im[\sigma_{SC}]}{\sigma_N}$} \label{fig:imag}
\end{figure}
As there is as yet no filtering based on the free energy, this is messy, so the comparison to \cref{fig:osRepro1b} is necessary to see which combinations of $T$ and $n$ are valid.
However, this is still quite numerically unstable, even with the dimensional reduction procedure described above for $n$ and $\corr$.
The information that higher $n$ leads to lower $\sigma$ is expected, and does suggest that the implementation is correct, although still very noisy.
\section{Noise calculation}
For our noise calculation, we assume a material parameterised by a Debye frequency $\debye$ and the interaction parameter $N(0) V$.
Because these aren't necessarily experimentally accessible, I tried to keep their values such that they lead to a physically reasonable $T_c$.
In order to keep $N(0) V$ small and $\debye$ bigger than the other parameters, I chose $N(0)V = 0.25$ and $\debye = \SI{1e13}{\per\s}$.
This leads to $T_{c0} = T_c(\mu = 0) = \SI{1.44e11}{\per\s}$, similar to the values used for the equilibrium Nam case.
Because that particular choice of $N(0) V$ and $\debye$ lead to some additional numerical noise for $T$ close to $T_c$, I changed the values slightly to $N(0) V = 0.2$ and $\debye = \SI{1e14}{\per\s}$, which leads to $T_{c0} = \SI{7.64e11}{\per\s}$, and ensures the graphs below look good for $T > 0.8 T_c$.
\begin{enumerate}
\item Use the Owen Scalapino coupled integral equations \cref{eq:gap,eq:n}, find $\mu$ and $\Delta$ for fixed $n$.
\item Find the expected gap from the approximation in OS, $T_c(n) \approx (1 - 4n) T_{c0}$.
If $T > T_c(n)$, then the calculation is skipped (a more complete handling would use either the Lindhard form or use a Nam expression that's been extended to $\Delta = 0$).
This is necessary because the coupled integral equations are very hard to solve
\item Using the modified Nam equations, calculate the dielectric function and create the approximated interpolation form, similar to the equilibrium case.
\item Calculate the noise as usual.
\end{enumerate}
\subsection{Figures}
Some examples of the noise are potrayed below, from \crefrange{fig:smallomeganoise}{fig:largerTnoise}.
There are some interesting features in the graphs:
\begin{itemize}
\item In the curves for constant omega, \crefrange{fig:smallomeganoise}{fig:bigomeganoise}, there seems to be a consistent dip, where initially as $T$ increases, the noise decreases.
That seems like it could be an artifact of the method used to calculate $\Delta$, but it is very odd (because nothing besides $\Delta$) should affect any other step.
\item In the constant omega graphs, for $n = .02$ and $n = .04$ the noise is greater than for $n = 0$ even in the region where it goes to the `normal state'.
This is confusing.
\item Most values of $n$ should not be superconducting for most of the range of $T$, as expected.
Near the critical temperature for a particular $n$, the noise increases as expected, showing a little gap.
See \cref{fig:mediumTnoise}.
\end{itemize}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constOmega1/1}
\caption{$T_1$ vs temperature, for very small omega} \label{fig:smallomeganoise}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constOmega1/10}
\caption{$T_1$ vs temperature, for medium omega} \label{fig:mediumomeganoise}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constOmega1/13}
\caption{$T_1$ vs temperature, for big omega} \label{fig:bigomeganoise}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constT1/43}
\caption{$T_1$ vs frequency, for a temperature where all the chosen $n$ values are allowed, showing the gap between $n = 0$ and $n > 0$, which seems potentially spurious.} \label{fig:smallTnoise}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constT1/35}
\caption{$T_1$ vs frequency, for a temperature where the $n = .04$ state is almost, but not quite surpressed.} \label{fig:mediumTnoise}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{constT1/15}
\caption{$T_1$ vs frequency, for temperature so large that only the $n = 0$ state is energetically allowed to be superconducting} \label{fig:largerTnoise}
\end{figure}
\printbibliography
\end{document}