nam_paper/paper.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{cleveref}

\title{EWJN from a BCS Superconductor}

\addbibresource{./bibliography.bib}

\graphicspath{{./figures/}}

\newcommand{\vf}{v_{\mathrm{F}}}
\newcommand{\qf}{q_{\mathrm{F}}}

\begin{document}

\maketitle

\section{Introduction \label{sec:intro}}
\begin{itemize}
	\item Motivate with~\cite{Tenberg2019} and~\cite{Kolkowitz2015}
	\item Existing work with Lindhard expression~\cite{QubitRelax}
	\item Variety of existing superconducting dielectric functions like in~\cite{AGD, llv9, Zimmermann1991, Mattis, Tinkham}.
	Advantage of Nam expression\cite{Nam1967} is that it's got the worked out dependence on momentum and impurities, both of which we need.
\end{itemize}
\section{Methods \label{sec:methods}}

Physically, we are interested in the relaxation time of a qubit near the surface of the metal.
Sufficiently close to the metal, such that the separation between qubit and metal is much smaller than the shortest dimension of the metallic body, we can approximate the metal as a half space.
This defines a natural coordinate system, and we can allow the metal to take up all points with $z$-coordinate less than zero, extending to infinity in the $x$ and $y$ directions.

For a charge qubit with level separation $\omega$ and dipole moment $\vec{d}$, the relaxation rate $\frac{1}{T_1}$ depends on the qubit's distance from the surface $z$, as well as its orientation $i$.
The vacuum wavelength $\lambda = \frac{c}{\omega}$ is a natural unit for this distance $z$, so we wil measure $z$ in units of $\lambda$.
The electromagnetic field fluctuations that contribute to qubit relaxation have been described in~\cite{QubitRelax} and~\cite{Henkel1999}.
Both use Fermi's golden rule and the fluctuation-dissipation theorem to relate the relaxation rate to the spectral density of the field fluctuations, and obtain the following expression:
\begin{equation}
	\frac{1}{T_1} = \frac{d^2}{\epsilon_0} \frac{\omega^3}{c^3} \chi_{i}^{(E)}(z, \omega) \coth\frac{\omega}{2 T}.
\end{equation}
%\todo{All the Nam stuff is in Gaussian units, so should pick one unit system and stick with it. Doesn't %affect results so far, as $\chi$ is unitless and only depends on quantities that are the same in SI / %Gaussian. Still bad though.}
%Here and elsewhere we take $\hbar = k_{\mathrm{B}} = 1$.

Similarly, for spin qubits with dipole moment $\vec{\mu}$, both~\cite{QubitRelax} and~\cite{Henkel1999} have a similar expression with a different spectral density expression:
\begin{equation}
	\frac{1}{T_1} = \frac{\mu^2}{\epsilon_0} \frac{\omega^3}{c^3} \chi_{i}^{(B)}(z, \omega) \coth\frac{\omega}{2 T}.
\end{equation}

The spectral densities are in general functions of the geometry and electrodynamic response of the metal object.
A standard limiting case for the response is that of local, linear electrodynamics, where the relationship between the electric displacement $\vec{D}$ and electric field $\vec{E}$ is a proportionality:
\begin{equation}
	\vec{D}(\vec{r}) = \epsilon \vec{E}(\vec{r})
\end{equation}
For a metal with conductivity $\sigma$, this dielectric function will typically have a large imaginary value as determined by the Drude expression in Gaussian units:
\begin{equation}
	\epsilon = 1 + i\frac{4 \pi \sigma}{\omega}
\end{equation}

For a half space, these spectral densities can be written in terms of the surface's reflection coefficients, as derived by\cite{Henkel1999} for local electrodynamics.
Considering for now a qubit pointing in the direction perpendicular to the half-space, we can write
\begin{align}
	\chi_\perp^E(z, \omega) &= \Re \int_0^{+\infty} \dd{u} \frac{u^3}{v} e^{2 i z v} r_p(u). \label{eq:chi}
\end{align}
The integration variable $u$ effectively represents a momentum in units of $\flatfrac{1}{\lambda}$, with $v = \sqrt{1 - u^2}$.
If $v \geq 1$, we take the positive square root $v = i \sqrt{u^2 - 1}$.
The magnetic spectral density is the same, except with an additional factor of $\flatfrac{1}{c^2}$ and with $r_s$ instead of $r_p$.
In the local limit, for dielectric constant $\epsilon$, the reflection coefficients will be the Fresnel $r_p$ and $r_s$:
\begin{align}
	r_p &= \frac{\epsilon v - \sqrt{\epsilon - u^2}}{\epsilon v + \sqrt{\epsilon - u^2}} \\
	r_s &= \frac{v - \sqrt{\epsilon - u^2}}{ v + \sqrt{\epsilon - u^2}}
\end{align}

However, as discussed in~\cite{QubitRelax} and~\cite{Henkel2006}, this expression no longer remains accurate for arbitrarily small distances, with \eqref{eq:chi} diverging as $\frac{1}{z^3}$ as $z \rightarrow 0$.
The divergence here stems from the unphysicality of local electrodynamics at very small scales;
for $z$ smaller than the electromagnetic coherence length of the metal the full response function defined by
\begin{equation}
	\vec{D}(\vec{r, t}) = \int \dd{r'} \dd{t'} \epsilon(r, r', t, t') \vec{E}(\vec{r', t'})
\end{equation}
becomes necessary.
The experiment described by Kolkowitz et al.\cite{Kolkowitz2015} gives an example of a physical situation where the local description is insufficient, with a qubit probe at a distance from a silver film less than the film's mean free path.
Given that the coherence lengths of typical BCS superconductors can be larger than the distances used in that experiment, we should expect that the purely local description will not work for the superconducting case.

We can keep~\eqref{eq:chi} accurate for sufficiently small $z$ if we make appropriate changes to our reflection coefficients\cite{QubitRelax,Henkel2006}, so we use the reflection coefficients in Ford and Weber~\cite{Ford1984}:
\begin{align}
	r_p(u) &= \frac{\pi v - \zeta_p(u)}{\pi v + \zeta_p(u)} \\
	r_s(u) &= \frac{\zeta_s(u) - \frac{\pi}{v}}{\zeta_s(u) + \frac{\pi}{v}} \\
	\zeta_p(u) &= 2i \int_0^\infty \dd{y} \frac{1}{\kappa^2} \left( \frac{y^2}{\epsilon(\frac{\omega}{c}\kappa, \omega) - \kappa^2} + \frac{u^2}{\epsilon_(\frac{\omega}{c}\kappa, \omega)} \right) \label{eq:zp} \\
	\zeta_s(u) &= 2i \int_0^\infty \dd{y} \frac{1}{\epsilon(\frac{\omega}{c}\kappa, \omega) - \kappa^2} \label{eq:zs} \\
	\kappa^2 &= u^2 + y^2
\end{align}
The quantities $\zeta_p$ and $\zeta_s$ are proportional to the surface impedances of the metal, derived for this case by Ford and Weber\cite{Ford1984}, and with more thorough discussion in~\cite{LandauLifshitzElectrodynamics}.
The treatment in~\cite{QubitRelax} compares the difference between these expressions and the simpler Fresnel reflection coefficients.

The dielectric function $\epsilon(q, \omega)$, then, contains the information needed to describe the electromagnetic properties of the surface near the qubit.
For the normal state, the dielectric function derived by Lindhard~\cite{Lindhard} used in~\cite{QubitRelax} describes the non-local electromagnetic response of a metal.
Using the form described by Solyom\cite{SolyomV3},
\begin{equation}
	\epsilon_{\mathrm{Lindhard}}(\vec{q}, \omega) = 1 + \frac{q_{TF}^2}{q^2}\frac{\displaystyle 1 + \frac{\omega + \flatfrac{i}{\tau}}{2 \vf q} \ln(\frac{\omega - \vf q + \flatfrac{i}{\tau}}{\omega + \vf q + \flatfrac{i}{\tau}})}{\displaystyle 1 + \frac{\flatfrac{i}{\tau}}{2 \vf q} \ln(\frac{\omega - \vf q + \flatfrac{i}{\tau}}{\omega + \vf q + \flatfrac{i}{\tau}})}. \label{eq:lindhardsolyom}
\end{equation}
Here, $q_{TF}$ is the Thomas-Fermi wavevector $q_{TF}^2 = 3 \flatfrac{\omega_p^2}{\vf^2}$, $\omega_p$ is the plasma frequency $\sqrt{\flatfrac{4 \pi n e^2}{m}}$, $\tau$ is the collision time and $\vf$ is the Fermi velocity.
The branch cut for the logarithm is chosen here such that their imaginary parts lie between $\pm i \pi$.
This can be shown to reduce to the Drude dielectric function in the $q \rightarrow 0$ limit.

The Lindhard dielectric function reflects the important property of having an imaginary part that vanishes for $q$ such that $\abs{\varepsilon_q - \omega} >  q\vf$.
This is a very generic feature of these types of response functions, and occurs because there are no points on the assumed spherical Fermi surface further than $2 q_{\mathrm{F}}$ apart.
Thus, there are no available quasiparticle-hole excitations available for energy dissipation (cf discussion in \cite{AGD}, \cite{FetterWalecka} or \cite{SolyomV3}).
This is a general argument, and it should be expected that a superconducting dielectric function should also have zero imaginary part above some momentum on the order of the Fermi momentum.
Because the Lindhard dielectric function's imaginary part vanishes above a cutoff $q_c\left(\omega\right)$ and has real part that goes as $\Re \epsilon_{\mathrm{Lindhard}} - 1 \sim \frac{1}{q^2}$, all of the integrals in $\eqref{eq:chi}, \eqref{eq:zp} and \eqref{eq:zs}$ converge.

We use the expressions from Nam in~\cite{Nam1967} to represent the superconducting response function.
This extends the previous models by Mattis and Bardeen~\cite{Mattis} and Abrikosov, Dzyaloshinskii and Gorkov\cite{AGD} to give expressions that allow for broader ranges of impurity values.\todo{Including the full expressions from Nam here is a bit space-prohibitive, but it may be important to show exactly what our assumptions encode to in his notation? ex: by assuming no magnetic impurities, our renormalisation factor becomes simpler.}

Here, we look at Nam's expressions in the weak coupling limit, for no magnetic impurities and an isotropic material.
\begin{equation}
	\epsilon_{\mathrm{Nam}}(\vec{q}, \omega) = 1 + \frac{3 \pi}{\omega^2} \frac{n e^2}{m} \left[\int_{\Delta - \omega}^{\Delta}\dd{\omega'} \tanh(\frac{\omega + \omega'}{2 T}) I_1 + \int_{\Delta}^{\infty} \dd{\omega'} \left( \tanh(\frac{\omega + \omega'}{2 T}) I_1  - \tanh(\frac{\omega'}{2 T})I_2 \right) \right], \label{eq:eps}
\end{equation}
with
\begin{align}
	I_1 &= F(q, \Re[\sqrt{(\omega + \omega')^2 - \Delta^2} - \sqrt{\omega'^2 - \Delta^2}]) (g + 1) \nonumber\\
	&\quad + F(q, \Re[-\sqrt{(\omega + \omega')^2 - \Delta^2} - \sqrt{\omega'^2 - \Delta^2}]) (g - 1) \\
	I_2 &= F(q, \Re[\sqrt{(\omega + \omega')^2 - \Delta^2} - \sqrt{\omega'^2 - \Delta^2}]) (g + 1) \nonumber\\
	&\quad + F(q, \Re[\sqrt{(\omega +  \omega')^2 - \Delta^2} + \sqrt{\omega'^2 - \Delta^2}]) (g - 1) \\
	F(q, E) &= \frac{1}{q \vf} \left[2 S(E) + (1 - S(E)^2)\ln(\frac{S(E) + 1}{S(E) - 1})\right] \label{eq:NamF} \\
	S(q, E) &= \frac{1}{q \vf} \left( E - i \left(\Im[\sqrt{(\omega + \omega')^2 - \Delta^2} + \sqrt{\omega'^2 - \Delta^2}] + \frac{2}{\tau} \right) \right) \\
	g &= \frac{\omega' \left(\omega + \omega'\right) + \Delta^2}{\sqrt{\omega'^2 - \Delta^2}\sqrt{(\omega + \omega')^2 - \Delta^2}}.
\end{align}
The assumption of isotropy suppresses the $q$ dependence for $\Delta$, which then is just a function of temperature, and can be described using the well-known BCS expression $\Delta \approx 3.06 \sqrt{T_c(T_c - T)}$ (see for example \cite{Tinkham}).

\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{Cond1Re}
	\caption{$\Re[\epsilon(q)]$ for $\omega = 1$, $\tau = 0.5$, $\omega_p = 10$, $\vf = 1$, $T = .9999 T_c$, $T_c = 3$} \label{fig:cond1Re}
\end{figure}

\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{Cond1Im}
	\caption{$\Im[\epsilon(q)]$ for $\omega = 1$, $\tau = 0.5$, $\omega_p = 10$, $\vf = 1$, $T = .9999 T_c$, $T_c = 3$} \label{fig:cond1Im}
\end{figure}

The Lindhard and Nam dielectric constants are compared in \cref{fig:cond1Re} and \cref{fig:cond1Im}, plotting the real and imaginary part for small representative values.
In this regime, $\omega_p > T_c > \omega$, as is typical for the frequency regime of interest, while $\tau$ is chosen to be smaller than $\omega$.
For a typical metal in this description, the Fermi wavevector $\qf$ is around the same order as $\sqrt{3}\frac{\omega_p}{\vf}$ (see discussion on this point in Solyom\cite{SolyomV3}).
We can see in \cref{fig:cond1Im} that the Lindhard dielectric function goes to zero for $q < \qf \approx 10 \sqrt{3}$.

\section{Numerical Techniques \label{sec:technical}}

The noise integral \eqref{eq:chi} can be calculated numerically, with proper care taken to handle the integrand's behaviour across the entire range.
For small momenta, where $\vf q \ll \omega$, both \eqref{eq:lindhardsolyom} and \eqref{eq:eps} can be series expanded to give explicit expressions.
The Lindhard dielectric function, up to $\mathcal{O}(q^2)$, becomes
\begin{gather}
	\epsilon_{\mathrm{Lindhard}}(\vec{q}, \omega) = 1 - \frac{\omega_p^2}{\omega^2} \left(\frac{\omega}{(\omega + \frac{i}{\tau})} + (\vf q)^2  \frac{9 \omega + 5 \frac{i}{\tau}}{15 (\omega + \frac{i}{\tau})^3} \right). \label{eq:lindhardsmallkseries}
\end{gather}
As expected for a description of the normal state, for $q \rightarrow 0$ this reduces to the Drude expression.

All of the momentum dependence in the Nam expression is contained within the $F(q, E)$ function in \eqref{eq:NamF}.
Expanding this out to second order in the momentum gives
\begin{align}
	F = \frac43 \frac{1}{\eta} + (\vf q)^2\frac{4}{15} \frac{1}{\eta^3},
\end{align}
where
\begin{equation}
	\eta = E - i \left(\Im[\sqrt{(\omega + \omega')^2 - \Delta^2} + \sqrt{\omega'^2 - \Delta^2}] + 2 \frac{1}{\tau} \right).
\end{equation}
This, and other limiting forms, were stated by Nam as well~\cite{Nam1967}.
Inserting this in $\eqref{eq:NamF}$ suffices to obtain the small $q$ values in a more numerically stable way.

By comparison to the $q \rightarrow 0$ case, the large momentum dependence is more involved to correctly handle.
For sufficiently large $q$, the logarithm in $\eqref{eq:NamF}$ goes to $i \pi$, and $F \rightarrow \frac{i \pi}{q \vf}$.
This behaviour leads to an unphysical divergence in \eqref{eq:chi} as $z \rightarrow 0$.
Tinkham~\cite{Tinkham} and Abrikosov, Dzyaloshinskii and Gorkov~\cite{AGD} both describe expressions with this same $\frac{1}{q}$ dependence.
The root cause of inaccuracies above $q \gg 2 q_{\mathrm{F}}$ is that the Green's function method used to derive the superconducting response function makes the assumption that $q$ is sufficiently close to the Fermi surface.

However, it is insufficient to only account for the $\frac{1}{q}$ dependence, as a stronger condition is necessary to ensure convergence for arbitrarily small $z$.
In the large $q$ regime, \eqref{eq:chi} can be reduced to
\begin{align}
	\chi_\perp^E(z, \omega) &= \Re \int_0^{+\infty} \dd{u} \frac{u^3}{i u} e^{-2 z u} r_p(u),
\end{align}
which means that for $z = 0$
\begin{align}
	\chi_\perp^E(z, \omega) &= \int_0^{+\infty} \dd{u} u^2 \Im r_p(u).
\end{align}
For a dielectric function that asymptotically scales as $\frac{A + i B}{q}$, for real $A$ and $B$, it can be shown that
\begin{equation}
	\Im r_p(u) \sim \frac{B}{u},
\end{equation}
which clearly leads to a divergent $\chi_\perp^E$.

This is the same divergence discussed in Langsjoen et al.~\cite{QubitRelax}, where for the normal state this problem does not arise if the Lindhard dielectric function is used.
Ultimately, this is because for $q > 2 q_{\mathrm{F}}$ the imaginary part of the Lindhard function goes to zero.

This can be handled in two coarse approximations.
The key is that the integral in \eqref{eq:chi} picks out values around $u = \frac{c}{\omega} \sim \frac{1}{z}$ over most of its range, because of the $u^2 e^{-2 z u}$ factor for $u \gg 1$.
If we cut off the imaginary part of the Nam response function above some momentum $q_{\mathrm{cutoff}}$, then if $\frac{1}{z} \ll \frac{\omega}{c} q_{\mathrm{cutoff}}$, the imprecision of such a low order approximation will not greatly affect the final noise integral.
For a metal with $\vf \sim \SI{e6}{\m\per\s}$, this corresponds to $z \sim \SI{1}{\nm}$.
At this length scale, inter-atomic spacing would itself be enough to invalidate these results already.

As a further correction, the Nam expression's $\frac{1}{q}$ dependence causes an overestimation in $\chi_\perp^E$ if $\frac{1}{z}$ is above the momentum at which superconductivity should play a role, on the order of the Debye momentum.
Above this point, the Lindhard result suffices as a description of the metal's response.
This can be incorporated into a numerical model by interpolating the results;
one crude way of accomplishing this is to simply define a new effective reflection coefficient
\begin{equation}
	r_{p, \mathrm{effective}}(u) = \min_{\Im}\left[r_{p, \mathrm{Lindhard}}\left(u\right), r_{p, \mathrm{Nam}}\left(u\right)\right],
\end{equation}
where the right hand side should be interpreted to choose the argument with the smaller imaginary part.
This effectively forces the reflection coefficient to choose the Lindhard result above some crossover $q_{\mathrm{crossover}}$.
Once again, unless $\frac{1}{z} \sim q_{\mathrm{crossover}}$, the inaccuracies this causes should be exponentially damped in $\chi_\perp^E$.

\section{Experiments \label{sec:experiments}}
\begin{itemize}
	\item Show that expression is of right order of magnitude to describe~\cite{Tenberg2019}.
	\item Look at temperature at which~\cite{Tenberg2019} should show SC effects (not far off!)
	\item Compare to~\cite{Kolkowitz2015}?
	Much further from relevant temperature range.
	\item Metal choices?
	Obviously large $T_c$ is most important to get noise reduction benefits.
	What else is most relevant for usable metals?
\end{itemize}
\section{Conclusions \label{sec:conclusions}}

\printbibliography

\end{document}