# Dissipation by surface states in superconducting RF cavities

Sean Deyo<EMAIL_ADDRESS>Michelle Kelley  Nathan Sitaraman  Danilo B.
Liarte  Tomas Arias  James P. Sethna  Laboratory of Atomic and Solid State
Physics, Cornell University, Ithaca, NY, USA  Thomas Oseroff  Matthias Liepe
Cornell Laboratory for Accelerator-Based Sciences and Education, Cornell
University, Ithaca, NY, USA

###### Abstract

Recent experiments on superconducting cavities have found that under large
radio-frequency (RF) electromagnetic fields the quality factor can improve
with increasing field amplitude, a so-called “anti-$Q$ slope.” Linear theories
of dissipation break down under these extreme conditions and are unable to
explain this behavior. We numerically solve the Bogoliubov-de Gennes equations
at the surface of a superconductor in a parallel AC magnetic field, finding
that at large fields there are quasiparticle surface states with energies
below the bulk value of the superconducting gap. As the field oscillates, such
states emerge and disappear with every cycle. We consider the dissipation
resulting from inelastic quasiparticle-phonon scattering into these states and
investigate the ability of this mechanism to explain features of the
experimental observations, including the field dependence of the quality
factor. We find that this mechanism is likely not the dominant source of
dissipation and does not produce an anti-$Q$ slope by itself; however, we
demonstrate in a modified two-fluid model how these bound states can play a
role in producing an anti-$Q$ slope.

##  I Introduction

Superconducting radio-frequency (SRF) cavities are useful in a variety of
modern applications, such as free-electron lasers and particle colliders [1].
The primary advantage of SRF cavities over their normal-conducting
alternatives is a lower surface resistance, which has the obvious benefit of a
lower energy footprint [2]. Despite decades of steady improvements, allowing
for quality factors in the neighborhood of $10^{10}$ [1] and accelerating
fields up to $25-45$ MV/m [2], some unanswered questions remain. At high
fields there is often a “$Q$ slope,” meaning the quality factor decreases—or,
equivalently, surface resistance increases—as field strength increases. More
confounding, recently some cavities, in particular niobium cavities doped with
nitrogen, have exhibited an “anti-$Q$ slope,” meaning that the quality factor
actually increases with field [1]. Increasing the resonant frequency of the
cavity has been associated with stronger anti-$Q$ slopes [3]. Since inelastic
scattering rates [4] can be comparable to typical SRF frequencies [5], this
suggests inelastic scattering may play an important role in determining the
field dependence of surface resistance.

Definitive theoretical explanations for the anti-$Q$ slope remain elusive.
Conventional theories of AC dissipation in superconductors predict no $Q$
slope (constant quality factor). Mattis and Bardeen [6] employed linear
response methods along with BCS theory [7] to calculate the complex
conductivity of a superconductor subject to a magnetic field. They found a
surface resistance that roughly goes as $\omega^{2}\exp(-\Delta/k_{B}T)$ at
low temperatures, where $\Delta$ is the superconducting gap and $\omega$ is
the angular frequency of the field. Central to these calculations is the idea
that a quasiparticle in a state with energy $E_{1}$ can absorb a photon of
energy $\hbar\omega$, and transition to a state of energy
$E_{2}=E_{1}+\hbar\omega$. (Here we ignore other sources of dissipation such
as trapped magnetic vortices oscillating near the cavity surface [8, 9, 10]
and within grain boundaries [11, 12].) In order for the photon energy to be
well defined the absorption must happen coherently over many cycles of the AC
perturbation.

Extensions of the linear-response theory of [6] have been routinely invoked to
interpret experimental data showing an anti-$Q$ slope. Gurevich applied the
Keldysh formalism in the context of non-equilibrium Green’s functions to
rederive the conductivity and evaluate the surface resistance at low
frequencies and mean free paths, and high magnetic fields [13]. The anti-$Q$
slope there originates in the smearing of the density of states due to
oscillating superflow. At high enough fields, the decrease of the
superconducting gap and subsequent smearing of the density of states weakens
the zero-frequency singularity of the Mattis-Bardeen resistance, leading to a
field dependence. This calculation involves an approximation to determine the
nonequilibrium quasiparticle distribution function, a key quantity that is
very difficult to find starting from a fundamental theory. Goldie and
Withington obtained non-linear solutions of the kinetic equations for the
coupled quasiparticle and phonon systems [14]. Their solutions for the non-
thermal quasiparticle distribution function were later combined with the
Mattis-Bardeen theory by de Visser et al., who proposed a mechanism for
microwave suppression on superconducting aluminum resonators [15].

Though some of these extensions can be used in the regimes of strong fields
[13, 16], they maintain the notion of $E_{1}\to E_{2}=E_{1}+\hbar\omega$. At
large fields, however, we shall see that there are quasiparticle surface
states whose energies change dramatically during each cycle. The situation is
reminiscent of the smeared density of states in [13], except that we
explicitly include depth and time dependence of both the superconducting gap
and the quasiparticle states. In particular, we find that $E_{2}(t)-E_{1}(t)$
is not constant for such states, making the equality $E_{2}=E_{1}+\hbar\omega$
questionable. This discrepancy demands a new framework to handle such states
in strong AC fields. Rather than consider transitions occurring coherently
over many cycles, we take the opposite limit and consider scattering events
occurring within a single cycle: solving the same quantum dissipation problem,
but starting from the adiabatic limit rather than the weak-coupling limit.

A simpler phenomenological model relying on the two-fluid idea [17] yields
explicit analytic formulas for the surface resistance but lacks essential
physical ingredients present in the Mattis-Bardeen calculation, such as
coherent effects on the transition probabilities [18]. We briefly discuss how
to incorporate the quasiparticle states with changing energies into a modified
version of the two-fluid model.

In Section II we introduce the Bogoliubov-de Gennes equations and describe
quasiparticle surface states whose energies exhibit strong field dependence.
In Section III we consider how quasiparticle-phonon scattering changes the
occupation of these surface states during the RF cycle and can lead to
dissipation. In Section IV we compute the field and frequency dependence of
our dissipation mechanism and discuss its relevance to the anti-$Q$ slope.
Finally, in Section V we offer concluding remarks and possibilities for
further study.

##  II Bogoliubov-de Gennes formalism

To find quasiparticle states we solve the Bogoliubov-de Gennes (BdG) equations
[19]:

$\displaystyle(H_{e}+U)u+\Delta v$ $\displaystyle=Eu$ (1)
$\displaystyle-(H_{e}^{*}+U)v+\Delta^{*}u$ $\displaystyle=Ev$

with

$H_{e}(\mathbf{r})=(-i\hbar\nabla-e\mathbf{A}(\mathbf{r})/c)^{2}/2m+U_{0}(\mathbf{r})-E_{F},$

where $e$ is the fundamental charge, $c$ is the speed of light, and $m$ is the
mass of an electron. These equations must be solved self-consistently with the
potentials given by

$\displaystyle U$
$\displaystyle=-V\sum_{n}|u_{n}|^{2}f_{n}+|v_{n}|^{2}(1-f_{n})$ (2)
$\displaystyle\Delta$ $\displaystyle=V\sum_{n}u_{n}v_{n}^{*}(1-2f_{n}),$

where $f_{n}$ is the occupation of state $n$ and $V$ is a constant describing
the strength of the interaction. To accomplish this numerically, we make an
initial guess for the potentials, find all the solutions $(u_{n},v_{n})$, use
them to compute a refined guess for the potentials, and iterate until the
process settles on a set of $(u_{n},v_{n})$ consistent with $U$ and $\Delta$.
Note that Eqns. (1) possess a symmetry: If $(u,v)$ is a solution with energy
$E$, $(v^{*},-u^{*})$ is a solution with energy $-E$. If one takes advantage
of this symmetry one must also replace the state’s occupation $f$ with $1-f$
in order to preserve Eqns. (2). In doing so, one treats the quasiparticle
state as a ‘quasi-hole’ state with opposite energy and occupation.

We solve Eqns. (1) at the surface of a superconductor in a parallel magnetic
field: $\mathbf{A}=A_{0}\sin(\omega t)e^{-z/\lambda}\hat{y}$, where $\lambda$
is the London penetration depth of the superconductor. For simplicity we
impose the vector potential externally, rather than using the current density
and Maxwell’s equations to solve for a new vector potential after each
iteration of the self-consistency process. However, we can confirm that our
current density, after reaching self-consistency with $U$ and $\Delta$, is
also consistent with our vector potential. This accords with observations that
$\lambda$ does not change significantly for RF fields below $H_{c1}$ [20].

For the sake of simplicity and generality we neglect the atomic potential
$U_{0}$ and assume a spherical Fermi surface. To avoid simulating a semi-
infinite half-space, we limit the domain to a depth $L$. Provided
$L\gg\lambda$, this truncation does not affect the energies of states
localized near the surface. We also take the occupation fractions to be in
equilibrium, $f_{n}=\left(1+e^{\beta E_{n}}\right)^{-1}$, when self-
consistently finding the potentials, having verified numerically that this
approximation is of little effect.

Factoring out the dependence in the directions parallel to the surface,
$u(\mathbf{r})\to e^{i(k_{x}x+k_{y}y)}u(z)$ and likewise for $v$, we obtain a
coupled system of differential equations:

$\displaystyle\left\\{\frac{\hbar^{2}}{2m}\left[k_{x}^{2}+\left(k_{y}-\frac{eA_{0}}{\hbar
c}\sin(\omega
t)e^{-z/\lambda}\right)^{2}-\frac{d}{dz}^{2}\right]-E_{F}+U(z)\right\\}u(z)+\Delta(z)v(z)$
$\displaystyle=Eu(z),$ (3)
$\displaystyle-\left\\{\frac{\hbar^{2}}{2m}\left[k_{x}^{2}+\left(k_{y}+\frac{eA_{0}}{\hbar
c}\sin(\omega
t)e^{-z/\lambda}\right)^{2}-\frac{d}{dz}^{2}\right]-E_{F}+U(z)\right\\}v(z)+\Delta(z)u(z)$
$\displaystyle=Ev(z),$

with $k_{x}$ and $k_{y}$ as parameters. Since the time dependence of the field
$A_{0}(t)$ is much slower than quantum relaxation times (apart from inelastic
scattering), we obtain the time dependence of the eigenstates simply by
solving these equations at a series of times.

For boundary conditions we set $u(0)=u(L)=v(0)=v(L)=0$ to confine the
quasiparticles to the slab. One can see from Eqns. (2) that this forces the
potentials to vanish at the surface, which might seem to be a problem given
that, for instance, $\Delta$ should be a nonzero constant if the field is
zero. However, the relevant length scale for features of $\Delta$ is the
correlation length, while the sums in Eqns. (2) include terms that oscillate
on the much shorter length scale of $1/k_{F}$. For niobium $1/k_{F}=0.08$ nm
[21], about $500$ times smaller than the correlation length of nearly $40$ nm
[5]. Thus, any oscillations due to forcing $\Delta(0)=0$ are only present
within a tiny fraction of a penetration depth from the surface, reminiscent of
the Friedel oscillations in electron density near a surface [22]. Others who
have studied the Bogoliubov-de Gennes equations near a surface observe the
same phenomenon [23], or a similar phenomenon for the BdG equations with the
tight-binding model [24]. Authors studying superconductivity in thin films and
wires have also noted oscillations in the energy gap as a function of
thickness [25, 26]. It is important to distinguish between oscillations in the
energy gap as a function of thickness in a thin film and oscillations in the
pair potential as a function of position in a thick slab or half-space, but
both effects ultimately result from boundary conditions. Alternative boundary
conditions, as for instance in the jellium model [27], could soften the
oscillations in the pair potential, but we would not expect them to disappear
entirely.

##  III Dissipation mechanism

In the absence of a field, there is a uniform pair potential
$\Delta(z)=\Delta_{0}$ and all quasiparticle energies must be larger than
$\Delta_{0}$. Once the field is applied, states with energies less than
$\Delta_{0}$ appear. These states are localized near the surface. Examples are
plotted in Figure 1.

Figure 1: The BdG quasiparticle bound states for $k_{x}=k_{F}$, $k_{y}=0$, at
a field of $40$ MV/m. For each state the wave functions $u$ and $v$ are
plotted with a vertical shift according to their energy. The situation is
somewhat analogous to a particle in a potential $\Delta(z)$ with a hard wall
on the left.

There is also a continuum of states, not localized to the surface, whose
energies change negligibly during the cycle. Because elastic scattering is
quick—mean-free path estimates suggest an elastic scattering time in the
femtosecond range, several orders of magnitude quicker than typical SRF field
frequencies in the GHz range [28]—we can assume that a state stays in
equilibrium as long as it can scatter elastically with this reservoir of
continuum states. Bound states with energies below $\Delta_{0}$ cannot scatter
elastically with continuum states, so they can only be filled by the much
slower inelastic processes. If the inelastic scattering rate is comparable to
or slower than $\omega$ there will be disequilibrium.

The fact that the scattering rate can change during the cycle complicates the
dynamics. For any fermion state it is true that

$\frac{df}{dt}=\frac{1-f}{\tau_{h}}-\frac{f}{\tau_{p}}$ (4)

where $\tau_{h}$ and $\tau_{p}$ are the hole and particle lifetimes,
respectively. Setting Eqn. (4) to zero allows us to relate the equilibrium
occupation $f_{0}$ to the two lifetimes:

$(1-f_{0})\tau_{p}=f_{0}\tau_{h}\qquad\text{or}\qquad
f_{0}=\frac{\tau_{p}}{\tau_{p}+\tau_{h}}.$ (5)

We can use these relations to eliminate either of the two lifetimes:

$\frac{df}{dt}=\frac{f_{0}-f}{f_{0}\tau_{h}}=\frac{f_{0}-f}{(1-f_{0})\tau_{p}}.$
(6)

At least one of the lifetimes must be time-dependent if the energy is time-
dependent, since $f_{0}=(1+e^{\beta E})^{-1}$. The simplest approximation is
to take one of the lifetimes to be constant and let the other be determined by
Eqn. (5). For a bound state at energy $E<\Delta_{0}$, $\tau_{p}$ is determined
in large part by the ability to absorb a phonon of energy $\Delta_{0}-E$. The
occupation numbers for such phonons depend exponentially on $\Delta_{0}-E$,
making $\tau_{p}$ strongly time-dependent. $\tau_{h}$, on the other hand, is
determined largely by spontaneous emission of phonons, making the occupation
numbers irrelevant. Thus, constant $\tau_{h}$ is the more physical choice for
a first approximation. Measurements in bulk niobium indicate an inelastic
scattering lifetime of $\tau_{p}=16$ ns at $2$ K for a quasiparticle at
$E=\Delta_{0}$ [4], which corresponds to $\tau_{h}=100$ µs.

To refine the approximation, let us consider what might make the hole lifetime
not constant. There is a simple phase space argument for the energy dependence
of $\tau_{h}$ in the low-temperature regime: Kaplan et al. provide well-known
formulas for the inelastic scattering rates of quasiparticles in terms of the
Eliashberg spectral function, $\alpha^{2}F(\Omega)$, which is approximately
quadratic in the phonon energy $\Omega$ at low frequencies [4]. Consider one
quasiparticle bound state with energy $E(t)<\Delta_{0}$: the typical
scattering processes at low temperatures are a quasiparticle in a continuum
state with energy $\Delta_{0}+k_{B}T$ emitting a phonon and entering the bound
state, and a quasiparticle in the bound state absorbing a phonon and entering
a continuum state with energy $\Delta_{0}+k_{B}T$. These rates are
proportional to

$\alpha^{2}F(\Omega)\propto\Omega^{2}\propto(\Delta_{0}+k_{B}T-E(t))^{2}.$

Furthermore, because bound state wave functions are concentrated within one
penetration depth of the surface, the scattering rate between a bound state
and a continuum state is roughly proportional to the local density of the
continuum state at the surface. We can compute the average density of
continuum quasiparticles within one penetration depth of the surface,
$\rho_{1}$, as a function of field and compare it to the average density at
zero field, $\rho_{0}$. Combining the phase space argument with the surface
density effect, we conclude that the scattering rate is approximately
proportional to $\alpha^{2}F(\Omega)\cdot\rho_{1}/\rho_{0}$. Thus, our refined
approximation for the hole lifetime is

$\tau_{h}(t)=\tau_{0}\frac{\rho_{0}}{\rho_{1}(t)}\left(\frac{k_{B}T_{c}}{\Delta_{0}+k_{B}T-E(t)}\right)^{2},$
(7)

where $T_{c}$ is the superconducting critical temperature and the constant
$\tau_{0}$ is set by demanding $\tau_{h}=100$ µs at $2$ K for a state at
$E=\Delta_{0}$.

Even with these refinements, our model of $\tau_{h}$ is still an
approximation. We have not explicitly considered possible effects of surface
ordering in enhancing $T_{c}$ at the surface [29]. This could be an important
point to consider for a more detailed model. Merely changing $T_{c}$ in Eqn.
(7) would have little effect, given that $T_{c}$ only appears as a
multiplicative constant and the coefficient $\tau_{0}$ fixes the overall
scaling of $\tau_{h}$; however, the physics of surface ordering (among other
things) could demand a depth-dependent expression for $\tau_{h}$.

With these scattering rates we can calculate $f(t)$. Every time the field goes
to zero the bound states return to the continuum, where they scatter
elastically with continuum states and return almost instantaneously to
$f=f_{0}$, dumping their excess occupation into the reservoir of continuum
states. The process repeats when the field starts to increase again and bound
states reemerge. Thus, for each bound state we can find the net work done by
the field in a single half cycle and multiply by $\omega/\pi$ to obtain the
power dissipated:

$P_{\text{diss}}=\frac{\omega}{\pi}\int_{0}^{T/2}f(t)\frac{dE}{dt}dt.$ (8)

In this form it is clear that the dissipation comes from the changing energies
of occupied states. Ultimately the excess energy is dumped via phonons, and to
illustrate this one can integrate by parts to obtain an equivalent formula:

$P_{\text{diss}}=\frac{\omega}{\pi}\int_{0}^{T/2}\frac{df}{dt}\left(E_{0}-E(t)\right)dt.$
(9)

The rate of phonon emission is $df/dt$, and the phonon energy is $E_{0}-E(t)$.
The formula holds for any $E_{0}$; that is, any energy the quasiparticle may
have come from before emitting the phonon. When using this formula one must be
careful to handle the delta function in $df/dt$ that occurs when the bound
state returns to the continuum at the end of the half cycle.

At large fields some bound states reach negative energies. When this happens
one can use the ‘quasi-hole’ symmetry described in Section II. If the
quasiparticle energy drops below $-\Delta_{0}$ or, equivalently, the quasi-
hole reaches an energy above $\Delta_{0}$, the state can briefly scatter
elastically with continuum states until it reenters the interval
$(-\Delta_{0},\Delta_{0})$.

Figure 2: Top: quasiparticle scattering from a bulk state into a bound state,
releasing a phonon in the process. As the bound state energy rises, the field
does work on the quasiparticle, resulting in net dissipation. Bottom: the
macroscopic picture tracking the average occupation of an ensemble of such
states for several values of $\omega\tau$. Dissipation is maximized when
${\omega\tau\approx 1}$.

Figure 2 illustrates the dynamics. For simplicity, the figure was made
assuming ${df/dt=-(f-f_{0})/\tau}$, which is the low-field limit of Eqn. (6)
with $\tau=\tau_{p}$. In this limit we can explicitly solve for $f(t)$:

$f(t)=e^{-t/\tau}\left[f_{0}(E(0))+\frac{1}{\tau}\int_{0}^{t}e^{s/\tau}f_{0}(E(s))ds\right].$

One can see that $f$ is larger in the second half of the plot because it lags
behind $f_{0}$. The second half is when $dE/dt$ is positive, so this imbalance
is the reason the field does positive net work. The size of the imbalance
depends on $\tau$: If inelastic scattering is much quicker than the field
(${\omega\tau\ll 1}$), the state stays near equilibrium the whole time ($f$ is
close to $f_{0}$) and little dissipation occurs. If inelastic scattering is
slow (${\omega\tau\gg 1}$), the state hardly ever gets filled ($f$ barely
changes) and again little dissipation occurs. The maximum imbalance occurs for
${\omega\tau\approx 1}$. At typical SRF frequencies of a few GHz, niobium’s
${\tau_{p}=16\text{ ns}}$ puts us in the ${\omega\tau\gg 1}$ regime.

We stress that this mechanism does not apply to unbound states, whose wave
functions oscillate throughout the slab and are not concentrated near the
surface. We discuss the dissipation due to these states in Section IV, using a
modified two-fluid formalism [18].

After computing the dissipation over an infinitesimal area $ds$, we find the
surface resistance $R_{\text{s}}$ from

$P_{\text{diss}}=\frac{1}{2}R_{\text{s}}\int|H|^{2}\,ds$ (10)

and the quality factor from

$Q=\frac{\omega
U}{P_{\text{diss}}}=\frac{\omega\mu_{0}\int|H|^{2}\,dv}{R_{\text{s}}\int|H|^{2}\,ds}=\frac{G}{R_{\text{s}}}$
(11)

where $G=\frac{\omega\mu_{0}\int|H|^{2}\,dv}{\int|H|^{2}\,ds}$ is the geometry
factor, a parameter independent of the frequency and size of the SRF cavity.
In order to, say, double the frequency of a cylindrical cavity, one must halve
the radius of the cavity in order to keep it in the same wave guide mode, so
$G$ remains unchanged and is in fact purely geometrical.

We use parameters realistic for an elliptical SRF cavity made from niobium: a
Fermi energy of $E_{F}=5.32$ eV [21], a bulk gap of $\Delta_{0}=1.5$ meV, a
penetration depth of $\lambda=40$ nm, a temperature of $2$ K, a geometry
factor of $G=270$ $\Omega$, a frequency of $1.3$ GHz, and accelerating
gradients up to $30$ MV/m [5]. Such gradients correspond to magnetic fields of
more than half of niobium’s lower critical field $H_{c1}$ at $2$ K [30].

Additionally, we can use the Bogoliubov-de Gennes states we calculate as
inputs into the two-fluid model, with superconducting electrons treated as an
inductive channel and normal conducting electrons as a resistive channel. With
these two channels in parallel, and in the limit that the inductive channel
has a much lower impedance, the AC dissipation is proportional to the density
of the normal fluid:

$R_{\text{s}}\approx\delta^{-1}\sigma_{1}/\sigma_{2}^{2}$ (12)

with

$\displaystyle\sigma_{1}$ $\displaystyle=n_{n}e^{2}\tau_{n}/m$ (13)
$\displaystyle\sigma_{2}$ $\displaystyle=n_{s}e^{2}/m\omega$

where $\delta\approx c\,(4\pi\omega\sigma_{2})^{-1/2}$ is the skin depth,
$n_{n}$ and $n_{s}$ are the densities of the normal and superconducting
electrons, and $\tau_{n}$ is the relaxation time of the normal electrons [18].
Note that $\tau_{n}$ is not the same as the inelastic scattering time; it is a
phenomenological parameter one can calculate from the electrical resistivity,
which tends to be dominated by elastic scattering with impurities at low
temperatures [21, 31].

To estimate the local fraction of normal electrons in the superconducting
phase, we compute the proportion of electrons relative to the normal phase
available for conduction

$\frac{n_{n}}{n}(z)=\frac{1}{N(0)}\sum_{n}\left(-\frac{\partial
f_{n}}{\partial E_{n}}\right)\left(u_{n}^{2}(z)+v_{n}^{2}(z)\right),$ (14)

where $N(0)$ denotes the electronic density of states at the Fermi level. This
quantity, and other normal-phase properties, are readily computed within the
Bogoliubov-de Gennes framework by forcing $\Delta=0$ in Eqns. (1). For a
homogeneous superconductor, Eqn. (14) reduces to the standard BCS form [18].

If the bound states have a different relaxation time than the continuum
states, we treat them as two parallel resistive channels with conductivity

$\sigma_{1}=\left(n_{b}\tau_{b}+n_{c}\tau_{c}\right)e^{2}/m,$ (15)

where $n_{b}$ and $n_{c}$ are the densities of bound and continuum states, and
$\tau_{b}$ and $\tau_{c}$ are their respective relaxation times. These
densities are calculated using Eqn. (14) but restricting the summation over
states with ${E_{n}<\Delta_{0}}$ or ${E_{n}\geq\Delta_{0}}$. We call this
modification the “three-fluid” model.

A decrease in $\sigma_{1}$ causes more current to flow through the inductive
channel, ultimately reducing $R_{\text{s}}$. Near the surface, as the field
increases the density of bound states increases and the density of continuum
states decreases. Accordingly, any source of scattering that makes
$\tau_{b}\ll\tau_{c}$ can cause a decrease in $\sigma_{1}$ and thus an
anti-$Q$ slope in our three-fluid model.

##  IV Results

Figure 3: Quality factor vs. field amplitude for a niobium cavity at frequency
$1.3$ GHz. The dashed lines show the quality factor computed from our
inelastic scattering disequilibrium mechanism, using either a constant
$\tau_{h}$ or the time-dependent $\tau_{h}(t)$ in Eqn. (7). The dash-dot lines
show the quality factor from this mechanism plus the three-fluid model as
described in the text with $\tau_{c}=1$ ps. Markers show experimental data
from Grassellino et al. [32] for a selection of electro-polished (EP) or
nitrogen-doped cavities.

Figure 3 plots the quality factor vs. field amplitude for a constant
$\tau_{h}=100$ $\mu$s and for the time-dependent $\tau_{h}(t)$ in Eqn. (7),
along with several experimental measurements from Grassellino et al. [32] for
comparison. Constant $\tau_{h}$ yields an anti-$Q$ slope, while the refined
$\tau_{h}$ has $Q$ declining as the field strength increases. Both
approximations yield quality factors well above the experimental values,
indicating that the mechanism of inelastic scattering into bound states is
likely not the dominant source of dissipation, but may become important for
cavities with quality factors approaching $10^{11}$. Exactly how important
depends in large part on the magnitude of $\tau_{h}$.

As noted earlier, the parameters for a typical niobium cavity put us in the
slow-scattering regime, wherein dissipation would increase if the scattering
rates were to increase. A tenfold increase in scattering rates would make the
quality factor ten times lower. Given the rough nature of our approximations
and the strong influence of $\tau_{h}$ on the dissipation, more precise
calculations of both the magnitude and field dependence of scattering rates
are likely warranted.

Another important consideration is the frequency dependence of the quality
factor. In the $\omega\tau\gg 1$ regime, dissipation is limited by the rarity
of scattering events that fill a given bound state. Doubling the frequency
leads to half as many scattering events and thus half as much dissipation per
cycle, but with twice as many cycles per second the energy dissipated per
second remains constant. On the other hand, in the $\omega\tau\ll 1$ regime
the dissipation is limited by the fact that the bound states stay close to
equilibrium. Increasing the frequency leads to more disequilibrium and more
cycles per second, so the result is that dissipation is proportional to
frequency squared. Numerical computation verifies that dissipation goes with
$\omega^{2}$ when $\omega\tau\ll 1$ and flattens out once $\omega\tau\gg 1$.
Inelastic scattering rates in bulk niobium suggest $\omega\tau\gg 1$, so
unless the filling of surface states is much faster than what we have inferred
from bulk estimates, dissipation from this mechanism is independent of
frequency within the typical SRF range.

We also compute the dissipation resulting from the three-fluid model, using
our Bogoliubov-de Gennes states to calculate the densities of the two normal
fluids. The relative sizes of $\tau_{c}$ and $\tau_{b}$ depend on the nature
of the scattering. Accurate descriptions of surface scattering will depend on
factors such as the cavity’s impurity profiles, surface roughness [33] or
coatings such as the Nb3Sn coating often applied to the surface of SRF
cavities [12]. When we model the elastic scattering of quasiparticles with
point-like impurities using a three-dimensional delta function as the
perturbing potential, we estimate using Fermi’s golden rule that $\tau_{b}$ is
shorter than $\tau_{c}$ for states near the surface by roughly a factor of
two. When we model an extended defect such as a surface nanohydride with a
one-dimensional delta function $\delta(z-z_{0})$, we find values for
$\tau_{b}$ as much as six times smaller than $\tau_{c}$.

Given the uncertainty about the ratio $\tau_{b}/\tau_{c}$, in Figure 3 we
assume $\tau_{c}=1$ ps [18] and plot the resulting quality factor for two
limiting cases, $\tau_{b}=0$ (yellow dash-dot curve) and $\tau_{b}=\tau_{c}$
(red dash-dot-dot curve), to show the range of plausible $Q$ slopes. The
results for $\tau_{b}=0$ roughly align with the anti-$Q$ slope in the
experimental data for the first $15$ MV/m, whereas the $\tau_{b}=\tau_{c}$
curve is similar to that of the electropolished cavity with no anti-$Q$ slope.

As the field strength increases, the density of the bound state fluid near the
surface grows while the density of the continuum state fluid shrinks. The
increase in bound states exceeds the decrease in continuum states, so if
$\tau_{b}=\tau_{c}$ the net effect is to increase $\sigma_{1}$, drawing more
current through the resistive channel and dissipating more energy. However, if
$\tau_{b}$ is negligible compared to $\tau_{c}$, the net effect is to decrease
$\sigma_{1}$, causing less dissipation and yielding an anti-Q slope. We
estimate an anti-$Q$ slope to appear approximately when
$\tau_{b}<\tau_{c}/{5}$. The magnitude of the quality factor could change
significantly depending on the value of $\tau_{c}$, but the main point is the
difference in slope: These results suggest that one could manipulate the $Q$
slope by modulating any scattering sources that affect the bound states more
than the continuum states.

It is worth noting that in the two- and three-fluid models, as in the
conventional linear response theory, the dissipation has an $\omega^{2}$
dependence. If a mechanism with such dependence produces an anti-$Q$ slope,
the addition of a mechanism with flat frequency dependence and no anti-$Q$
slope—such as our inelastic scattering disequilibrium mechanism—would explain
why the anti-$Q$ slope becomes more pronounced at higher frequencies. While
this argument is speculative, it nonetheless demonstrates the potential for
bound states to contribute to the anti-$Q$ slope.

##  V Conclusion

For superconductors in large AC fields there are quasiparticle states for
which a linear response approach to dissipation is inadequate. These states
are bound to the surface with energies below the value of the bulk
superconducting gap, and their energies change throughout an AC cycle. By
focusing on these fundamentally nonperturbative bound states we have made a
stark departure from the conventional theories derived in a weak field
setting. For the range of parameters we expect in a typical niobium SRF
cavity, the dissipation from inelastic quasiparticle-phonon scattering into
these states is independent of frequency and is likely not the dominant source
of dissipation. However, if inelastic scattering is much stronger at the
surface of the superconductor than it is in bulk, the dissipation from this
mechanism could become important. More refined calculations of scattering
rates are a sensible next step to build on these results.

Moreover, we argue that when the relaxation time of the bound states differs
from that of the continuum states, the two-fluid model should be modified to
consider the continuum and bound quasiparticles as two separate normal-
conducting fluids in addition to the superconducting fluid. Depending on
whether the relaxation time of the bound state fluid is comparable to or much
smaller than that of the continuum state fluid, the resulting quality factor
could either increase or decrease with field strength. Our estimations suggest
that some kinds of scattering can indeed affect the bound states more than the
continuum states. If one can identify and control the concentration of such
scattering sources in an SRF cavity, in principle one could use them to
activate the anti-$Q$ slope.

##  VI Acknowledgments

This work was supported by the U.S. National Science Foundation under Award
PHY-1549132, the Center for Bright Beams.

## References

  * Dhakal [2020] P. Dhakal, Physics Open 5, 100034 (2020).
  * Padamsee [2017] H. Padamsee, Superconductor Science and Technology 30, 053003 (2017).
  * Martinello  _et al._ [2018] M. Martinello, M. Checchin, A. Romanenko, A. Grassellino, S. Aderhold, S. K. Chandrasekaran, O. Melnychuk, S. Posen, and D. A. Sergatskov, Physical Review Letters 121, 224801 (2018).
  * Kaplan  _et al._ [1976] S. B. Kaplan, C. C. Chi, D. N. Langenberg, J. J. Chang, S. Jafarey, and D. J. Scalapino, Physical Review B 14, 4854 (1976).
  * Gonnella [2016] D. A. Gonnella,  _The Fundamental Science of Nitrogen-Doping of Niobium Superconducting Cavities_ , Ph.D. thesis, Cornell University (2016).
  * Mattis and Bardeen [1958] D. C. Mattis and J. Bardeen, Phys. Rev. 111, 412 (1958).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * Gurevich and Ciovati [2013] A. Gurevich and G. Ciovati, Phys. Rev. B 87, 054502 (2013).
  * Liarte  _et al._ [2018] D. B. Liarte, D. Hall, P. N. Koufalis, A. Miyazaki, A. Senanian, M. Liepe, and J. P. Sethna, Phys. Rev. Applied 10, 054057 (2018).
  * Checchin  _et al._ [2018] M. Checchin, M. Martinello, A. Grassellino, S. Aderhold, S. K. Chandrasekaran, O. S. Melnychuk, S. Posen, A. Romanenko, and D. A. Sergatskov, Applied Physics Letters 112, 072601 (2018), https://doi.org/10.1063/1.5016525 .
  * Sheikhzada and Gurevich [2017] A. Sheikhzada and A. Gurevich, Phys. Rev. B 95, 214507 (2017).
  * Carlson  _et al._ [2021] J. Carlson, A. Pack, M. K. Transtrum, J. Lee, D. N. Seidman, D. B. Liarte, N. S. Sitaraman, A. Senanian, M. M. Kelley, J. P. Sethna, T. Arias, and S. Posen, Phys. Rev. B 103, 024516 (2021).
  * Gurevich [2014] A. Gurevich, Physical Review Letters 113, 087001 (2014).
  * Goldie and Withington [2013] D. J. Goldie and S. Withington, Superconductor Science and Technology 26, 015004 (2013).
  * de Visser  _et al._ [2014] P. J. de Visser, D. J. Goldie, P. Diener, S. Withington, J. J. A. Baselmans, and T. M. Klapwijk, Phys. Rev. Lett. 112, 047004 (2014).
  * Kubo and Gurevich [2019] T. Kubo and A. Gurevich, Physical Review B 100, 064522 (2019).
  * Halbritter [1974] J. Halbritter, Zeitschrift für Physik 266, 209 (1974).
  * Tinkham [2013] M. Tinkham,  _Introduction to Superconductivity_ (Dover Publications, 2013).
  * de Gennes [2018] P.-G. de Gennes,  _Superconductivity of Metals and Alloys_ (CRC Press, 2018).
  * Carlson and McLean [1968] D. E. Carlson and W. L. McLean, Physical Review 176, 538 (1968).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Harcourt, 1976).
  * Crommie  _et al._ [1993] M. Crommie, C. P. Lutz, and D. Eigler, Nature 363, 524 (1993).
  * Troy and Dorsey [1995] R. J. Troy and A. T. Dorsey, Physical Review B 51, 11729 (1995).
  * Croitoru  _et al._ [2020] M. Croitoru, A. Shanenko, Y. Chen, A. Vagov, and J. A. Aguiar, Physical Review B 102, 054513 (2020).
  * Blatt and Thompson [1963] J. M. Blatt and C. J. Thompson, Physical Review Letters 10, 332 (1963).
  * Shanenko and Croitoru [2006] A. Shanenko and M. Croitoru, Physical Review B 73, 012510 (2006).
  * Lang and Kohn [1970] N. Lang and W. Kohn, Physical Review B 1, 4555 (1970).
  * Maniscalco  _et al._ [2017] J. T. Maniscalco, D. Gonnella, and M. Liepe, Journal of Applied Physics 121, 043910 (2017).
  * Ginzburg [1964] V. Ginzburg, Phys. Letters 13, 101 (1964).
  * French [1968] R. A. French, Cryogenics 8, 301 (1968).
  * Webb [1969] G. W. Webb, Physical Review 181, 1127 (1969).
  * Grassellino  _et al._ [2013] A. Grassellino, A. Romanenko, D. Sergatskov, O. Melnychuk, Y. Trenikhina, A. Crawford, A. Rowe, M. Wong, T. Khabiboulline, and F. Barkov, Superconductor Science and Technology 26, 102001 (2013).
  * Jacob  _et al._ [1990] U. Jacob, V. Vancea, and H. Hoffman, Physical Review B 41, 11852 (1990).


# Automated robotic intraoperative ultrasound for brain surgery

Michael Dyck*1,2, Alistair Weld*3, Julian Klodmann1, Alexander Kirst1, Giulio
Anichini4,
Luke Dixon4, Sophie Camp4, Stamatia Giannarou3, Alin Albu-Schäffer1,2
1Institute of Robotics and Mechatronics, German Aerospace Center, 2Department
of Informatics, Technical University of Munich, 3Hamlyn Centre for Robotic
Surgery, Imperial College London, 4Department of Neurosurgery, Charing Cross
Hospital, Imperial College London, UK*These authors contributed equally to the
work -<EMAIL_ADDRESS><EMAIL_ADDRESS>This project was
supported by the International Graduate School of Science and Engineering
(IGSSE); TUM Graduate School, and the UK Research and Innovation (UKRI) Centre
for Doctoral Training in AI for Healthcare (EP/S023283/1), the Royal Society
(URF$\setminus$R$\setminus$2 01014]), the NIHR Imperial Biomedical Research
Centre. The authors want to thank ImFusion GmbH, and especially Dr. Marco
Esposito, for providing us an academic license to their software and constant
support throughout our research.

##  I Introduction

During brain tumour resection, localising cancerous tissue and delineating
healthy and pathological borders is challenging, even for experienced
neurosurgeons and neuroradiologists [1]. Intraoperative imaging is commonly
employed for determining and updating surgical plans in the operating room.
Ultrasound (US) has presented itself a suitable tool for this task, owing to
its ease of integration into the operating room and surgical procedure.
However, widespread establishment of this tool has been limited because of the
difficulty of anatomy localisation and data interpretation.

Publications on various applications of robotic US scanning have been produced
in recent years. Robotic ultrasound scanning of the spine is a well-posed
tracking task. [2] uses a semantic segmentation network to segment the spine,
and in conjunction with a depth map, generates a robotic trajectory. [3]
presents an automated method for lung diagnosis, using a pose estimation
network to identify regions on the chest for path defining. [4] uses image
quality optimization and interaction force adjustment to optimize probe pose
and surface coupling.

In this work, we present a robotic framework designed and tested on a soft-
tissue-mimicking brain phantom, simulating intraoperative US (iUS) scanning
during brain tumour surgery. Our work builds on what was introduced in [5, 6].

Figure 1: Experimental setup showing the robotic arm with it’s attached probe
and camera holder. Scanning a custom-made brain phantom.

##  II Methods

The workflow of our robotic iUS tissue scanning platform is depicted in Fig.
2.

Figure 2: Workflow for our iUS tissue scanning platform. Automatic detection
of the object’s location and RGB-D reconstruction of its surface with a stereo
camera results in a triangular mesh of the phantom. The impedance controller
incorporates the geometry of this mesh into the real-time control loop,
guiding the US probe to establish contact and execute scanning trajectories.
Graphic based on [6].

###  II-A Soft-Tissue-Mimicking Brain Phantoms

A custom-made soft-tissue mimicking phantoms was made using GELITA® GELATINE
Type Ballistic 3 with a brain-like silicone mold. To mimic soft tissue and
extend the durability of the phantom, the phantom contained
water:glycerine:gelatine at 45:45:10 wt.%. Inspired by [7], different objects
were placed inside the phantom, including: olives, grapes, blueberries,
screws. This ensures the presence of random features within the US recordings
of the phantom.

###  II-B Localisation and Surface Reconstruction

For localisation, four ArUco markers are first manually placed around the
area. Then the camera is positioned above the area of interest, such that all
markers are visible within the image. This is enough to initialise the
automatic localisation and reconstruction routine. Using the OpenCV [8] ArUco
library, the four fiducial markers can be identified. Using these points, a
plane - defined by a centre point and normal vector - containing the phantom
is calculated. Using the detected plane, the robot automatically aligns the
camera such that it’s z-axis points towards the plane’s centre at an angle of
45°and a distance of 30 cm. Using the RGB-D reconstruction algorithm - from
the ImFusion Suite software solution by ImFusion GmbH - and executing a simple
rotational trajectory - rotating the camera around the normal of the
identified plane - a triangle mesh of the tissue phantom is extracted.

###  II-C Impedance-Controlled US Scanning

We use the impedance controller introduced in [9]. Surface-specific
coordinates are defined on the surface mesh. These coordinates are the
distance $d\in\mathbb{R}$ between transducer and surface, and the three
orientation coordinates $\boldsymbol{\epsilon}\in\mathbb{R}^{3\times 1}$ that
align the probe axis and surface normal. Two additional coordinates
$(s_{1},s_{2})\in\mathcal{X}\subset\mathbb{R}^{2}$ control the position of the
transducer on the surface. We stack all coordinates in a vector
$\boldsymbol{\rho}=(s_{1},s_{2},d,\epsilon_{1},\epsilon_{2},\epsilon_{3})^{T}\in\mathbb{R}^{6\times
1}$ (cmp. Fig. 2) and implement the unified impedance control framework by
[10]

$\boldsymbol{\tau}=\boldsymbol{J}_{\rho}^{T}(\boldsymbol{q})[\boldsymbol{K}_{\rho}(\boldsymbol{\rho}_{d}-\boldsymbol{\rho}(\boldsymbol{q}))+\boldsymbol{D}_{\rho}(\boldsymbol{q})(\boldsymbol{\dot{\rho}}_{d}-\boldsymbol{\dot{\rho}})].$
(1)

Here, $\boldsymbol{K}_{\rho},\boldsymbol{D}_{\rho}$ represent the symmetric,
positive-definite stiffness and damping matrices, respectively.
$\boldsymbol{q}\in\mathbb{R}^{7\times 1}$ are the generalized configuration
coordinates of the robot,
($\boldsymbol{\rho}_{d},\boldsymbol{\dot{\rho}}_{d}$) are the desired position
and velocity in the coordinates, $\boldsymbol{J}_{\rho}\in\mathbb{R}^{6\times
6}$ represents the Jacobian matrix mapping joint velocities
$\boldsymbol{\dot{q}}$ to velocities $\boldsymbol{\dot{\rho}}$. Controlling
the distance results in passive interaction forces between US transducer and
tissue, that depend on controller and object stiffnesses.

##  III Results

The experimental setup can be seen in Fig. 1. A seven degrees of freedom (DoF)
DLR Miro surgical robotic arm is used. Attached is an Intel® RealSense™ Depth
Camera D435i, and a GE 12L-RS linear US probe. The US transducer is connected
to a GE LOGIQ e US machine. An ATI Mini45 F/T-sensor measures six DoF
interaction forces during scanning. Results of the marker detection and
surface 3D reconstruction can be seen in Fig. 2. Fig. 3 shows the interaction
force between the US probe and phantom. The US probe was aligned normal to the
surface ($\boldsymbol{\varepsilon}_{d}=\boldsymbol{0}_{3\times 1}$). The
setpoint for the distance coordinate ($d_{d}$) was continuously decreased
(increasing penetration depth) and held constant at a manually chosen value
based on qualitative US image quality. The example US images were recorded
during the experiment and show different stages of contact establishment -
ranging from very little contact on the left to good probe-tissue coupling on
the right. By controlling the distance, rather than the force, safe coupling
of probe and surface can be achieved - invariant to tissue properties and
geometry.

Figure 3: Interaction force between US transducer and brain phantom during
contact establishment. The US images below were recorded during the
interaction and show qualitatively different stages of US contact.

##  IV Conclusion

This work presents a framework for automated robotic iUS brain scanning. Where
testing on a complex custom-made brain phantom shows that site localisation
and surface reconstruction can be achieved and enables probe servoing with low
interaction force. Future work will concern visual servoing for probe pose and
US quality optimisation, while minimising required interaction forces.

## References

  * [1] L. Dixon, A. Lim, M. Grech-Sollars, D. Nandi, and S. J. Camp, “Intraoperative ultrasound in brain tumor surgery: A review and implementation guide,” _Neurosurgical Review_ , vol. 45, pp. 2503 – 2515, 2022.
  * [2] M. Victorova, M. K.-S. Lee, D. Navarro-Alarcon, and Y. P. Zheng, “Follow the curve: Robotic ultrasound navigation with learning based localization of spinous processes for scoliosis assessment,” _IEEE Access_ , vol. PP, pp. 1–1, 2021.
  * [3] X. Ma, Z. Zhang, and H. K. Zhang, “Autonomous scanning target localization for robotic lung ultrasound imaging,” _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 9467–9474, 2021.
  * [4] Z. Jiang, M. Grimm, M. Zhou, Y. Hu, J. Esteban, and N. Navab, “Automatic force-based probe positioning for precise robotic ultrasound acquisition,” _IEEE Transactions on Industrial Electronics_ , vol. 68, no. 11, pp. 11 200–11 211, 2021.
  * [5] A. Weld, M. Dyck, J. Klodmann, G. Anichini, L. Dixon, S. J. Camp, A. Albu-Schäffer, and S. Giannarou, “Collaborative robotic ultrasound tissue scanning for surgical resection guidance in neurosurgery,” _Proceedings of The 14th Hamlyn Symposium on Medical Robotics_ , vol. abs/2301.08174, 2022.
  * [6] M. Dyck, A. Weld, J. Klodmann, A. Kirst, L. Dixon, G. Anichini, S. J. Camp, S. Giannarou, and A. Albu-Schäffer, “Collaborative robotic ultrasound tissue scanning for surgical resection guidance in neurosurgery,” _IEEE Transactions on Medical Robotics and Bionics_ , 2022, submitted for publication.
  * [7] H. Morehouse, H. P. Thaker, and C. Persaud, “Addition of metamucil to gelatin for a realistic breast biopsy phantom,” _Journal of Ultrasound in Medicine_ , vol. 26, 2007.
  * [8] Itseez, “Open source computer vision library,” https://github.com/itseez/opencv, 2015.
  * [9] M. Dyck, A. Sachtler, J. Klodmann, and A. Albu-Schäffer, “Impedance control on arbitrary surfaces for ultrasound scanning using discrete differential geometry,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 7738–7746, 2022.
  * [10] A. Albu-Schäffer, C. Ott, and G. Hirzinger, “A unified passivity-based control framework for position, torque and impedance control of flexible joint robots,” _The international journal of robotics research_ , vol. 26, no. 1, pp. 23–39, 2007.