# Why dust pressure matters in debris discs

Elliot M. Lynch,1 Joshua B. Lovell,2 and Antranik A. Sefilian3
1 Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France
2 Centre for Astrophysics, Harvard & Smithsonian, 60 Garden St, Cambridge, MA
02138, United States
3Astrophysikalisches Institut und Universitätssternwarte, Friedrich-Schiller-
Universität Jena, Schillergäßchen 2–3, D-07745 Jena, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

There is a common assumption in the particulate disc community that the
pressure in particulate discs is essentially zero and that the disc
streamlines follow Keplerian orbits, in the absence of self-gravity or
external perturbations. It is also often assumed that the fluid description of
particulate discs is not valid in the presence of crossing orbits (e.g. from
nonzero free eccentricities). These stem from the misconception that fluid
pressure arises due to the (typically rare) collisions between particles and
that the velocity of particles in fluids are single-valued in space. In
reality, pressure is a statistical property of the particle distribution
function which arises precisely because there is a distribution of velocities
at a given position. In this letter we demonstrate, with simple examples, that
pressure in particulate discs is non-zero and is related to the inclination
and free eccentricity distributions of the constituent particles in the discs.
This means many common models of debris discs implicitly assume a nonzero, and
potentially quite significant, dust pressure. We shall also demonstrate that
the bulk motion of the dust is not the same as the particle motion and that
the presence of pressure gradients can lead to strong departures from
Keplerian motion.

###### keywords:

planetary systems - circumstellar matter - submillimetre: planetary systems -
celestial mechanics

††pubyear: 2022††pagerange: Why dust pressure matters in debris discs–Why dust
pressure matters in debris discs

##  1 Introduction

Debris discs are circumstellar discs of dust and other solid particles
typically left over from the dispersal of a gaseous protoplanetary disc, along
with the collisional grind-down of planetesimal belts. Dynamically particles
in the disc interact with the large bodies in the system through gravity and,
in the case of stars, through the stellar radiation field and winds. While
collisions between particles in the disc are important for setting the
planetesimal and grain size distribution, they are typically not frequent
enough to strongly affect the particle orbits (Wyatt et al., 1999). Thus, the
only important collective effect, on individual particle orbits, is that due
to self-gravity (see Wyatt, 2008; Matthews et al., 2014; Hughes et al., 2018).

Observations, however, are not sensitive to individual particle orbits but are
in fact sensitive to bulk (i.e. averaged) dust properties such as the dust
surface density. The bulk motion of the dust can depart strongly from the
motion of the individual dust particles as a result of the dust pressure
gradients, which, as we shall show, can be important even in completely
collisionless system. This can be particularly important in very narrow rings
such as Fomalhaut, as imaged by Kalas et al. (2005), Boley et al. (2012),
MacGregor et al. (2017) and Gáspár et al. (2023); HD53143 (MacGregor et al.,
2022; Stark et al., 2023); HR4796 (Kennedy et al., 2018; Schneider et al.,
2018; Olofsson et al., 2019) and HD202628 (Faramaz et al., 2019). The
narrowness of the ring means that even a small dust pressure can result in
order unity corrections to the bulk motion from the large pressure gradients.

In this letter we demonstrate the connection between dust pressure and the
free eccentricity and inclination distributions in dust discs. We shall show
how the presence of dust pressure can lead to the bulk dust motion strongly
departing from Keplerian motion, even when every particle is constrained to
follow a Keplerian ellipse. Finally, we shall discuss the importance of dust
pressure in dynamical and morphological models of continuous debris discs
(e.g. Davydenkova & Rafikov, 2018; Marino et al., 2019; Sefilian et al., 2021,
2023; Lovell et al., 2021, 2022; Lynch & Lovell, 2022; Lovell & Lynch, 2023).

The importance of pressure to the dynamics of completely collisionless systems
is a well studied phenomena in galactic dynamics and has been known about
since Jeans (1915). The results we present here are a consequence of Jeans
Equations and the dust pressure has a close relationship to the velocity
ellipsoid used to characterise the distribution functions of galaxies (Binney
& Tremaine, 2008). In addition to Jeans equations, there is existing work, in
the galactic dynamics literature, linking velocity dispersion to gravitational
softening (Miller, 1971; Tremaine, 2001; Touma, 2002). While the
correspondence is not exact (velocity dispersion effects also occur in the
absence of 2-body interactions), Tremaine (2001) has shown that velocity
dispersion, gas pressure and potential softening have similar effects on
spiral density waves propagating in particle discs.

In this letter we shall only consider dust pressure in the absence of gas.
Dust pressure is also important in the presence of gas, particularly for
larger grains or for time varying or turbulent flows. The presence of gas
(particularly turbulent gas) makes the problem significantly more complicated
and is unnecessary for understanding the main effects of dust pressure in gas-
poor debris discs.

In Section 2 we shall show how dust pressure arises in a circular disc with an
inclination distribution and how this leads to sub-Keplerian orbital motion.
In Section 3 we shall provide a brief sketch of how this effect works in a
coplanar disc with nonzero free eccentricity. Finally Section 4 discusses our
results in the context of existing models and Section 5 presents our
conclusions.

##  2 Inclination Distribution

The simplest setup that can demonstrate the importance of dust pressure is to
consider a disc where all the particles are on circular orbits with a constant
inclination, $\iota$, with respect to the midplane and a random longitude of
ascending node. Consider a point, $p$, on the midplane with (spherical) radius
$r$, azimuthal angle $\phi$ and polar angle $\theta=\pi/2$. The particles
which pass through this point will be split into two families, family 1 with
the longitude of ascending node aligned with $\phi$ and family 2 where it is
anti-aligned (see Figure 1). Note that both of these families orbit in a
prograde direction.

By simple geometric arguments, one can show that the velocity (in spherical
polars) of particles in family 1, as they pass through p, is $\Omega\equiv
v^{\phi}=\Omega_{K}\cos\iota$ and $v^{\theta}=\Omega_{K}\sin\iota$, where
$\Omega_{K}$ is the circular Keplerian angular velocity at radius $r$. The
velocity in family 2 is $v^{\phi}=\Omega_{K}\cos\iota$ and
$v^{\theta}=-\Omega_{K}\sin\iota$, with both families having $v^{r}=0$.
Assuming the total mass of dust particles in family 1 & 2 are the same then
the bulk (i.e. mass weighted average) vertical velocity $u^{\theta}=0$, and
the bulk azimuthal velocity $u^{\phi}=\Omega_{K}\cos\iota$. Thus, the bulk
motion of the dust in the midplane consists of circular sub-Keplerian motion,
which only matches the particle motion when $\iota=0$.

Figure 1: Two families of inclined orbits in a circular disc with constant
inclination, $\iota$, and anti-aligned longitude of ascending nodes. At a
given point on the midplane, if there is the same mass in family 1 and family
2, the averaged vertical velocity is zero and the average azimuthal velocity
(as defined in Section 2) is $u^{\phi}=\Omega_{K}\cos\iota$.

We now show how dust pressure gradients are responsible for the sub-Keplerian
character of the velocity. As particles in family 1 all have the same orbital
elements, the distribution function for these particles can be obtained as a
product of delta distributions. The distribution for family 2 can be
constructed similarly. Assuming the masses in family 1 and family 2 are the
same, and the mass crossing the midplane on the spherical shell with radius,
$r$, is $M_{r}(r)$ then the distribution function, $f$, for the particles on
the midplane $(\theta=\pi/2)$ is

$\displaystyle\begin{split}f=&\frac{1}{2}\Sigma^{\circ}(r)\delta(v^{r})\delta(v^{\phi}-\Omega_{K}\cos\iota)\\\
&\times\left[\delta(v^{\theta}-\Omega_{K}\sin\iota)+\delta(v^{\theta}+\Omega_{K}\sin\iota)\right],\end{split}$
(1)

where $\Sigma^{\circ}=\frac{M_{r}}{2\pi r}$ and $\delta(x)$ is the Dirac-delta
distribution. One can take velocity moments of Equation 1 to calculate the
bulk properties of the particulate disc (which is typically what observations
are sensitive to). The zeroth moment $\Sigma=\int
f\,dv^{r}\,dv^{\phi}\,dv^{\theta}=\Sigma^{\circ}$ is the surface density on
the midplane. The bulk velocity components, $u^{i}$, are given by
$u^{i}=\frac{1}{\Sigma}\int v^{i}f\,dv^{r}\,dv^{\phi}\,dv^{\theta}$, with
$i=\\{r,\phi,\theta\\}$, which, as expected, yield the same velocities
described in the previous paragraph.

We can also compute the second velocity moment of the distribution function,
the dust pressure tensor, $P^{ij}$, which arises irrespective of the
collisionality of the particles,

$P^{ij}=\int(v^{i}-u^{i})(v^{j}-u^{j})f\,dv^{r}\,dv^{\phi}\,dv^{\theta}.$ (2)

As all particles have $v^{r}=u^{r}=0$ and $v^{\phi}=u^{\phi}$, the components
of the pressure tensor which involve $r$ or $\phi$ are zero and the only
nonzero pressure tensor component is $P^{\theta\theta}$ as $u^{\theta}=0$.
However, the individual orbits have velocities of
$v^{\theta}=\pm\Omega_{K}\sin\iota$, which results in the following pressure
tensor component,

$\displaystyle\begin{split}P^{\theta\theta}&=\int(v^{\theta})^{2}f\,dv^{r}\,dv^{\phi}\,dv^{\theta},\\\
&=\Sigma^{\circ}\Omega_{K}^{2}\sin^{2}\iota.\end{split}$ (3)

The momentum equation for the bulk motion of the dust, which arises from the
first velocity moment of the (collisionless) Boltzmann equation, is

$\Sigma D\mathbf{u}=-\Sigma\nabla\Phi-(\nabla\cdot\mathbf{P}),$ (4)

where $\Phi$ is the gravitational potential and
$D=\partial_{t}+\mathbf{u}\cdot\nabla$ is the Lagrangian time derivative with
respect to the dust flow. In galactic dynamics, Equation (4) is known as Jeans
equation (Jeans, 1915; Binney & Tremaine, 2008). The last term,
$-(\nabla\cdot\mathbf{P})$, is the divergence of the pressure tensor. In a gas
disc, this corresponds to the usual pressure gradient term $-\nabla P$. In
dust discs (e.g. gas-poor debris discs) this term is usually neglected with
the justification that the only force on the particle is due to the point mass
potential of the central star. However, such a justification is incorrect as
$-(\nabla\cdot\mathbf{P})$ is not due to the forces on the particles but is
instead a statistical property that arises as a consequence of the particle
velocities being multivalued at any given position in the disc. For the
pressure tensor associated with our distribution function (Equation 1) this is

$(\nabla\cdot\mathbf{P})=-rP^{\theta\theta}\hat{\mathbf{e}}_{r},$ (5)

close to the midplane, where $\hat{\mathbf{e}}_{r}$ is the radial unit vector.
Looking for the steady solutions to the momentum equation (Equation 4), near
the midplane, one finds:

$\displaystyle-r[\Omega^{2}+(u^{\theta})^{2}]$
$\displaystyle=-r\Omega_{K}^{2}+r\frac{P^{\theta\theta}}{\Sigma^{\circ}},$ (6)
$\displaystyle r^{-1}u^{r}\Omega$ $\displaystyle=0,$ (7) $\displaystyle
r^{-1}u^{r}u^{\theta}$ $\displaystyle=0.$ (8)

By Equation (7), $u^{r}=0$ and by symmetry we require $u^{\theta}=0$.
Substituting the expression for $P^{\theta\theta}$ into Equation (6) we arrive
at $\Omega^{2}=\Omega_{K}^{2}\cos^{2}\iota$, which matches that derived
directly from the distribution function. Thus, we see that the pressure tensor
gradients are necessary to obtain the correct velocity for the bulk motion,
which is sub-Keplerian due to pressure support, similar to gas discs. If one
were to instead neglect the pressure gradients in the momentum equation, one
would instead obtain $\Omega=\Omega_{K}$ and incorrectly arrive at a Keplerian
bulk dust motion.

This result will carry over to a more realistic inclination distribution (e.g.
Brown, 2001; Matrà et al., 2019) with the added complication of evaluating the
integral in Equation 3. Finally, in discs where only one of the two families
exist, there are no intersections between particle orbits and the velocity
remains single valued. Thus, the pressure term vanishes, $\mathbf{P}=0$, and,
by Equation (7) one has $\Omega^{2}+(u^{\theta})^{2}=\Omega_{K}^{2}$
corresponding to an inclined Keplerian disc.

##  3 Free Eccentricity Distribution

Figure 2: In a disc with only free eccentricity, the bulk motion of the dust
consists of circular, non-Keplerian motion.

Another situation where dust pressure is important is in axisymmetric discs
with nonzero free eccentricity (see Figure 2). For simplicity, one can
consider the limiting case where the semimajor axis distribution is
infinitesimally narrow. Here, the radial velocity of the particle orbits
averages to zero. Thus, the bulk fluid motion consists of circular, non-
Keplerian, motion. That the bulk dust motion is non-Keplerian should be
obvious from the fact that the inner disc edge is located at the pericentre of
the particle orbits and, therefore, the particle velocities must be higher
than the circular Keplerian values (similarly, the orbits are slower than the
circular Keplerian value on the outermost orbit).

As with the inclination distribution, if one considers a point, $p$, with
cylindrical radius, $r$, and azimuthal angle, $\phi$, then there are two
orbits which pass through this point (as depicted Figure 2); one with a
velocity of $v^{r}=\sqrt{\frac{GM}{\lambda}}e\sin\psi(r,\phi)$ and
$v^{\phi}=\sqrt{\frac{GM}{\lambda^{3}}}[1+e\cos\psi(r,\phi)]^{2}$; and the
second with the same $v^{\phi}$ but
$v^{r}=-\sqrt{\frac{GM}{\lambda}}e\sin\psi(r,\phi)$ (where $\lambda$ and
$\psi(r,\phi)$ are the semi-latus rectum and true anomaly of the particle,
respectively). This results in the (mass-weighted) mean velocity of $u^{r}=0$
and $u^{\phi}=\Omega_{K}\sqrt{\frac{\lambda}{r}}$. Figure 3 shows how this
angular velocity compares with the Keplerian value for different
eccentricities. It is evident that the angular velocity is super-Keplerian in
the inner region of the disc and sub-Keplerian (potentially highly so) in the
outer region.

Figure 3: Angular velocity in a disc with free eccentricity from $e_{p}=0.1$
(solid black line) to $e_{p}=0.9$ (dashed black line) in intervals of $0.1$.
$a_{0}$ is the, arbitrary, semimajor axis of the orbits and $w=2a_{0}e_{p}$ is
the ring width. The discs are super-Keplerian at small radii and sub-Keplerian
at large radii. For the most eccentric discs the angular velocity is limited
to being less than twice the Keplerian value in order for the particles to
remain bound, however, the outer regions of the disc can become almost
entirely pressure supported with an arbitrarily low rotation rate. In galactic
dynamics this effect is a well known consequence of the virial theorem (Binney
& Tremaine, 2008).

We shall not repeat the same calculation we did with the circular disc with an
inclination distribution (i.e. Equation 2), as the integration to obtain the
pressure tensor is more involved. This is because, even if we consider the
simplest case of a set of a disc with a constant free eccentricity and an
infinitesimally narrow semimajor axis distribution (so that each individual
ring has a constant surface density, Statler, 1999), these rings cannot be
directly summed in phase space in the orbital coordinate system, as the
orbital coordinates are functions of both position and velocity. Regardless of
the difficulty in analytically evaluating the pressure tensor, the bulk
velocity within the disc will satisfy the momentum equation (Equation 4),
provided the pressure gradient term is included.

##  4 Validity of morphological and dynamical models

###  4.1 Implication for particle debris disc models

Debris disc observations are often interpreted using Monte-Carlo particle
based-methods, which sample the orbits in the disc and typically include
constant free and forced eccentricities (e.g. MacGregor et al., 2017; Kennedy,
2020; MacGregor et al., 2022). By sampling the particle distribution directly,
this method implicitly accounts for the dust pressure. This will, however,
mean the bulk motion around the ring is non-Keplerian, particularly when the
ring is narrow due to the large pressure gradients present (for instance, the
bulk motion in the disc models of Fomalhaut and HD53143 by MacGregor et al.
(2017) and MacGregor et al. (2022), respectively, have strong departures from
Keplerian motion). One challenge for Monte-Carlo models is their slow
convergence at high resolution (Rafikov, 2023; Lovell & Lynch, 2023). This,
however, will be further compounded if the velocity distribution needs to be
adequately sampled to account for the effects of dust pressure.

###  4.2 Implication for parameterised debris disc models

One question that the fluid description may help to answer is the question of
the validity of simplified models for debris disc morphology in the presence
of orbital intersections, such as those developed by Marino et al. (2019) and
Lynch & Lovell (2022) (and pioneered by Tremaine (1995); Statler et al.
(1999); Statler (1999); Peiris & Tremaine (2003) in the galactic dynamics
context). The model of Lynch & Lovell (2022) postulates that the density in
the debris disc satisfies the continuity equation for a fluid flow which
follows a set of non-overlapping, confocal, elliptical Keplerian orbits -
similar to the eccentric gas disc theories of Ogilvie (2001), Ogilvie & Barker
(2014) and Ogilvie & Lynch (2019). The validity of this assumption is thus
likely the same as that for a fluid disc, i.e. they should provide a good
description of the bulk dust properties on an orbital timescale when the
forces due to pressure gradients are weak enough, relative to the Keplerian
terms, to act on secular timescales. This is likely to be the case for
geometrically thin (scale height $H\ll R$ the radial lengthscale), but not for
radially wide debris discs (with radial lengthscales much longer than $H$)
with small to moderate free eccentricities. As such, the model presented in
Lynch & Lovell (2022) will break down for geometrically thick, or sufficiently
narrow debris discs (i.e. discs with $H\sim R$ or radial variations on
lengthscales comparable with $H$), in the same way that the secular theories
for eccentric gas discs break down under similar circumstances (Ogilvie, 2001;
Ogilvie & Barker, 2014; Ogilvie & Lynch, 2019).

The model of Pan et al. (2016) models eccentric discs using 1D line densities,
and predicts the disc surface brightness at long wavelengths to be $\propto
T_{\rm eff}/v_{\rm orb}$, where $T_{\rm eff}$ is equilibrium temperature of
the disc and $v_{\rm orb}$ is the orbital velocity. This is on a less secure
footing as the dust does not obey the zero-th moment equation. Lynch & Lovell
(2022) showed that the model of Pan et al. (2016) gives the correct behaviour
in the low (observational) resolution limit of the more complete theory. Thus,
the requirements for the Pan et al. (2016) model to be valid are very
restrictive in that the dust pressure must be low, whilst the beam size must
be much larger than the ring width, but not so large that it is comparable
with the disc radius. The issue with applying the model of Pan et al. (2016)
to a very narrow disc is that a relatively small free eccentricity, or
inclination distribution, will lead to non-Keplerian bulk dust motion around
the ring. We note that in situations where the model of Pan et al. (2016) is
valid, so too must be the models of Lynch & Lovell (2022). Finally, we note
that Lovell & Lynch (2023) recently presented a parameterised model which
includes free eccentricity and thus implicitly includes dust pressure, similar
to the particle based models.

###  4.3 Implication for dynamical debris disc models

For dynamical models of debris discs (e.g. Wyatt et al., 1999; Davydenkova &
Rafikov, 2018; Sende & Löhne, 2019; Sefilian et al., 2021, 2023), the presence
of dust pressure is potentially even more important. The morphological models,
discussed in Section 4.2, are primarily interested in describing dust
properties on orbital timescales, in order to obtain the dust surface density,
and thus require very large pressure gradients to be strongly modified by the
dust pressure. Dynamical models, however, are often concerned with describing
the evolution of the disc on secular timescales. This means pressure can be
much weaker and still significantly affect the dynamics on the timescales of
interest. Notably, for lower mass discs, it is likely that pressure, not self-
gravity, is the dominant collective effect and may need to be accounted for in
models like those in Davydenkova & Rafikov (2018), Sefilian et al. (2021) and
Sefilian et al. (2023). In fact, a simple consideration of the virial balance
implies that the correction to the (bulk) angular velocity for a typical disc
inclination distribution can greatly exceed that due to self-gravity.

To demonstrate the above, we consider a pessimistic scenario (a thin but
massive disc), where the disc mass is $M_{d}=3\times 10^{-3}M_{*}$ (this
includes the mass of the planetesimals which contain the large majority of the
total disc mass) and a mean particle inclination above/below the midplane of
$\iota\approx H/R=0.02$. We can then estimate the angular velocity correction
due to self gravity to be

$\displaystyle\begin{split}|\omega_{\rm
SG}/\Omega_{K}|&\approx\left(\frac{GM_{d}}{R^{3}}\right)^{1/2}\left(\frac{GM_{*}}{R^{3}}\right)^{-1/2}\\\
&=\sqrt{\frac{M_{\rm d}}{M_{\rm*}}}\lesssim 5\times 10^{-2}\left(\frac{M_{\rm
d}}{3\times 10^{-3}M_{*}}\right)^{1/2}.\end{split}$ (9)

Similarly, based on the results from Section 2, the angular velocity
correction due to dust pressure is approximately

$|\omega_{\rm p}/\Omega_{K}|\approx\sqrt{1-\cos^{2}\iota}\approx\iota\gtrsim
2\times 10^{-2}\left(\frac{H/R}{0.02}\right),$ (10)

where we have neglected any contribution from free eccentricity. Thus we see
that, even in this pessimistic scenario, the angular velocity correction from
dust pressure is comparable to that from self-gravity. For thicker or less
massive discs the correction due to pressure can potentially exceed that due
to self-gravity by several orders of magnitude (for similar reasons, the
effects of dust pressure can also exceed that due to some planetary
perturbations). There is, however, a crucial difference between dust pressure
and self-gravity in that self gravity is able to change the orbits of the
particles, whereas pressure in the collisionless/non-interacting limit is not.
It is primarily through its interaction with other dynamical processes that
the influence of pressure on the dynamics of the disc will be felt. The dust
pressure is likely the physical mechanism by which the disc self-gravity is
softened (where softened gravity is required) when the finite size of the
particles is unimportant (Tremaine (2001); see also Sefilian & Rafikov (2019)
for a detailed discussion of various softening prescriptions employed in
nearly Keplerian discs.

A dynamical model where dust pressure is implicitly present is Wyatt et al.
(1999). Wyatt et al. (1999) has shown that the (linear) perturbations from an
eccentric planet result in an eccentric debris disc where the complex
eccentricity of the particles is a sum of free and forced eccentricity. In the
fluid picture, the planet produces a force in the momentum equation which
affects the dust bulk velocity (resulting in a forced eccentricity), along
with heating the dust fluid. The fluid thermal motion takes the form of the
distribution of free eccentricities in the disc. There is no distinction
between dynamical heating and thermodynamic heating in the large particle
number limit. Softened N-ring models (Touma, 2002; Hahn, 2003; Sefilian et
al., 2023) also imply the presence of the dust pressure, but only when the
rings overlap (although the softening is motivated by velocity dispersion à la
Tremaine (2001) and Hahn (2003)). This results in a subtle difference with
Wyatt et al. (1999) in that free eccentricity in a single ring in Wyatt et al.
(1999) results in an axisymmetric disc with dust pressure (similar to that
considered in Section 3), while in Sefilian et al. (2023) one instead has a
zero pressure precessing eccentric ring. The dynamics of a (softened) N-ring
model when the rings share a common semimajor axis has been considered by
Touma (2002); Touma & Sridhar (2012); Touma & Tremaine (2014). Notably Touma &
Tremaine (2014) showed that, when self gravity is included, the thermal
axisymmetric eccentricity distribution (as considered in Section 3) becomes
unstable to developing an asymmetry, with azimuthal wavenumber $m=1$, if the
particle eccentricities are low enough.

Finally, one consequence of the nonzero dust pressure in debris discs is the
existence of dust sound/density waves. These sound waves are closely linked to
the possibility of having an arbitrary initial mass distribution around an
orbit. In relatively wide debris discs, free sound waves quickly become
extremely tightly wound as a result of phase mixing on $\sim 10$’s of orbits
and have little effect on observations or dynamics. However, if the disc’s
spread in semimajor axes is small (making phase mixing slow), or the waves are
being actively excited (e.g. by a planet), these waves may have important
dynamical and observational effects that are worth exploring further. It is
plausible that the spiral feature seen in the simulation of Sende & Löhne
(2019) can be interpreted as a form of spiral-density wave, excited by the
planet, but propagating in the disc as a consequence of dust pressure, rather
than the more usually considered self-gravity (with the added complication of
radiation pressure affecting the orbits of small particles and artificial
diffusion in orbital element space required for stability).

In conclusion, it should be possible to rederive several well-known results of
debris disc dynamics from the fluid perspective (e.g. the shaping of an
eccentric ring by a companion as described in Wyatt et al. (1999)). Of course,
unlike fluid discs, there is no guarantee that one can find closed form
expressions for the evolution of the dust pressure111There does exists an
exact closure for the Keplerian test particle problem when the collisionless
Boltzmann equation is written in orbital elements space which might prove
useful for this., so it remains to be seen whether there are aspects of debris
disc dynamics which are easier to treat with the fluid description. A closure,
when two body relaxation is included was obtained by Touma & Tremaine (2014).
We leave the task of finding a closed form fluid description of debris discs
for future work. We note that one can sidestep the closure problem by directly
working with the distribution function, as the fluid description is just a
tool for interpreting and/or solving the collisionless Boltzmann equation.

##  5 Conclusion

We have shown, with two simple examples, the importance of dust pressure in
dust discs, even when the individual particles are on exact Keplerian orbits.
Our findings have important implications for the study of gas-poor debris
discs. The pressure in debris discs is not related to additional forces on
particles, or collisions between particles, but is instead related to the non-
zero velocity dispersion present in debris discs. One can interpret the
inclination and free eccentricity distributions expected in typical debris
discs to be equivalent to (anisotropic) dust pressure. The presence of dust
pressure in debris discs means that the bulk dust motion does not match the
particle orbits and is not, generally, Keplerian. Finally, observations are
primarily sensitive to bulk dust properties not individual particle orbits, so
dust pressure is likely important in their interpretation.

## Acknowledgements

We thank the reviewer for their suggestions to improve the manuscript clarity.
This research was supported by the ERC through the CoG project PODCAST No
864965. This project has received funding from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement No 823823. J. B. Lovell would like to thank the Smithsonian
Institute for a Submilimeter Array Fellowship. A. A. Sefilian acknowledges
support by the Alexander von Humboldt Foundation through a Humboldt Research
Fellowship for postdoctoral researchers.

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Boley et al. (2012) Boley A. C., Payne M. J., Corder S., Dent W. R. F., Ford E. B., Shabram M., 2012, ApJ, 750, L21
  * Brown (2001) Brown M. E., 2001, AJ, 121, 2804
  * Davydenkova & Rafikov (2018) Davydenkova I., Rafikov R. R., 2018, ApJ, 864, 74
  * Faramaz et al. (2019) Faramaz V., et al., 2019, AJ, 158, 162
  * Gáspár et al. (2023) Gáspár A., et al., 2023, Nature Astronomy, 7, 790
  * Hahn (2003) Hahn J. M., 2003, ApJ, 595, 531
  * Hughes et al. (2018) Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541
  * Jeans (1915) Jeans J. H., 1915, MNRAS, 76, 70
  * Kalas et al. (2005) Kalas P., Graham J. R., Clampin M., 2005, Nature, 435, 1067
  * Kennedy (2020) Kennedy G. M., 2020, Royal Society Open Science, 7, 200063
  * Kennedy et al. (2018) Kennedy G. M., Marino S., Matrà L., Panić O., Wilner D., Wyatt M. C., Yelverton B., 2018, MNRAS, 475, 4924
  * Lovell & Lynch (2023) Lovell J. B., Lynch E. M., 2023, MNRAS,
  * Lovell et al. (2021) Lovell J. B., et al., 2021, MNRAS, 506, 1978
  * Lovell et al. (2022) Lovell J. B., et al., 2022, MNRAS, 517, 2546
  * Lynch & Lovell (2022) Lynch E. M., Lovell J. B., 2022, MNRAS, 510, 2538
  * MacGregor et al. (2017) MacGregor M. A., et al., 2017, ApJ, 842, 8
  * MacGregor et al. (2022) MacGregor M. A., et al., 2022, ApJ, 933, L1
  * Marino et al. (2019) Marino S., Yelverton B., Booth M., Faramaz V., Kennedy G. M., Matrà L., Wyatt M. C., 2019, MNRAS, 484, 1257
  * Matrà et al. (2019) Matrà L., Wyatt M. C., Wilner D. J., Dent W. R. F., Marino S., Kennedy G. M., Milli J., 2019, AJ, 157, 135
  * Matthews et al. (2014) Matthews B. C., Krivov A. V., Wyatt M. C., Bryden G., Eiroa C., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 521–544 (arXiv:1401.0743), doi:10.2458/azu_uapress_9780816531240-ch023
  * Miller (1971) Miller R. H., 1971, Ap&SS, 14, 73
  * Ogilvie (2001) Ogilvie G. I., 2001, MNRAS, 325, 231
  * Ogilvie & Barker (2014) Ogilvie G. I., Barker A. J., 2014, MNRAS, 445, 2621
  * Ogilvie & Lynch (2019) Ogilvie G. I., Lynch E. M., 2019, MNRAS, 483, 4453
  * Olofsson et al. (2019) Olofsson J., et al., 2019, A&A, 630, A142
  * Pan et al. (2016) Pan M., Nesvold E. R., Kuchner M. J., 2016, ApJ, 832, 81
  * Peiris & Tremaine (2003) Peiris H. V., Tremaine S., 2003, ApJ, 599, 237
  * Rafikov (2023) Rafikov R. R., 2023, MNRAS, 519, 5607
  * Schneider et al. (2018) Schneider G., et al., 2018, The Astronomical Journal, 155, 77
  * Sefilian & Rafikov (2019) Sefilian A. A., Rafikov R. R., 2019, MNRAS, 489, 4176
  * Sefilian et al. (2021) Sefilian A. A., Rafikov R. R., Wyatt M. C., 2021, ApJ, 910, 13
  * Sefilian et al. (2023) Sefilian A. A., Rafikov R. R., Wyatt M. C., 2023, arXiv e-prints, p. arXiv:2305.00951
  * Sende & Löhne (2019) Sende J. A., Löhne T., 2019, A&A, 631, A141
  * Stark et al. (2023) Stark C. C., Ren B., MacGregor M. A., Howard W. S., Hurt S. A., Weinberger A. J., Schneider G., Choquet E., 2023, ApJ, 945, 131
  * Statler (1999) Statler T. S., 1999, ApJ, 524, L87
  * Statler et al. (1999) Statler T. S., King I. R., Crane P., Jedrzejewski R. I., 1999, AJ, 117, 894
  * Touma (2002) Touma J. R., 2002, MNRAS, 333, 583
  * Touma & Sridhar (2012) Touma J. R., Sridhar S., 2012, MNRAS, 423, 2083
  * Touma & Tremaine (2014) Touma J., Tremaine S., 2014, arXiv e-prints, p. arXiv:1401.5534
  * Tremaine (1995) Tremaine S., 1995, AJ, 110, 628
  * Tremaine (2001) Tremaine S., 2001, AJ, 121, 1776
  * Wyatt (2008) Wyatt M. C., 2008, ARA&A, 46, 339
  * Wyatt et al. (1999) Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., Holmes E. K., Piña R. K., 1999, ApJ, 527, 918


Now at ]Max Planck Institute for Chemical Physics of Solids

# Local magnetic response of superconducting Sr2RuO4 thin films and rings

G. M. Ferguson  Department of Physics, Laboratory of Atomic and Solid State
Physics, Cornell University, Ithaca, New York 14853, USA  [  Hari P. Nair
Department of Materials Science and Engineering, Cornell University, Ithaca,
New York 14853, USA  Nathaniel J. Schreiber  Department of Materials Science
and Engineering, Cornell University, Ithaca, New York 14853, USA  Ludi Miao
Department of Physics, Laboratory of Atomic and Solid State Physics, Cornell
University, Ithaca, New York 14853, USA  Department of Physics, New Mexico
State University, Las Cruces, NM 88003  Kyle M. Shen  Department of Physics,
Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New
York 14853, USA  Kavli Institute at Cornell for Nanoscale Science, Ithaca, New
York 14853, USA  Darrell G. Schlom  Department of Materials Science and
Engineering, Cornell University, Ithaca, New York 14853, USA  Kavli Institute
at Cornell for Nanoscale Science, Ithaca, New York 14853, USA  Leibniz-
Institut für Kristallzüchtung, Max-Born-Str. 2, 12489 Berlin, Germany  Katja
C. Nowack<EMAIL_ADDRESS>Department of Physics, Laboratory of Atomic and
Solid State Physics, Cornell University, Ithaca, New York 14853, USA  Kavli
Institute at Cornell for Nanoscale Science, Ithaca, New York 14853, USA

###### Abstract

We conduct local magnetic measurements on superconducting thin-film samples of
Sr2RuO4 using scanning Superconducting Quantum Interference Device (SQUID)
susceptometry. From the diamagnetic response, we extract the magnetic
penetration depth, $\lambda$, which exhibits a quadratic temperature
dependence at low temperatures. Although a quadratic dependence in high-purity
bulk samples has been attributed to non-local electrodynamics, our analysis
suggests that in our thin-film samples the presence of scattering is the
origin of the quadratic dependence. While we observe micron-scale variations
in the diamagnetic response and superconducting transition temperature, the
form of the temperature dependence of $\lambda$ is independent of position.
Finally, we characterize flux trapping in superconducting rings
lithographically fabricated from the thin films, paving the way to systematic
device-based tests of the superconducting order parameter in Sr2RuO4.

††preprint: APS/123-QED

Since the discovery of superconductivity in Sr2RuO4 [1], a substantial
research effort has been dedicated to the determination of the superconducting
order parameter [2, 3, 4, 5, 6]. Although early experiments were interpreted
in favor of an odd-parity, time reversal symmetry breaking order parameter,
recent experiments have called this picture into question [7, 8, 9, 10, 11].
High-quality superconducting thin films of Sr2RuO4 have recently been grown by
molecular beam epitaxy [12], raising the prospect of device-based tests of the
order parameter symmetry [13]. Given the extreme sensitivity of the
superconducting state in Sr2RuO4 to disorder [14], a direct comparison between
the superconducting properties of thin films and bulk single crystals is
desirable.

The dependence on temperature, $T$, of the magnetic penetration depth,
$\lambda$, contains information about the gap structure of superconductors. A
nodal superconducting gap gives rise to a power-law dependence of $\lambda(T)$
at low temperature, while a fully-gapped superconductor exhibits an
exponential temperature dependence. In Sr2RuO4, the observation of
$\lambda(T)\sim T^{2}$ in high-quality bulk single crystals provided early
evidence for the presence of nodes in the superconducting gap function [15].

Here, we study the superconducting penetration depth in a
$30\text{\,}\mathrm{nm}\text{/}$ thick Sr2RuO4 thin film grown by molecular-
beam expitaxy. The substrate is
$(\mathrm{LaAlO_{3}})_{\mathrm{0.29}}-(\mathrm{SrAl}_{0.5}\mathrm{Ta}_{0.5}\mathrm{O}_{3})_{0.71}$
(LSAT) with the c-axis oriented perpendicular to the substrate surface. The
LSAT substrate preserves the tetragonal symmetry of Sr2RuO4 and induces a
small tensile strain of approximately $0.045\%$ when the sample is cooled to
low temperatures. Details of the growth of these films are described in Ref.
[12] and in Appendix C.Micrometer-scale rings were fabricated using
photolithography and argon ion milling from a second thin film with thickness
$30\text{\,}\mathrm{nm}\text{/}$ grown on LSAT (see Appendix D for details on
the fabrication).

Our experimental approach is illustrated schematically in Fig. 1a. We use a
scanning SQUID susceptometer [16, 17] to detect the diamagnetic response of
the superconducting film to a local magnetic field produced by a current
$I_{FC}$ sourced through an integrated field coil with radius $a$. In the
superconducting state, the Sr2RuO4 film generates supercurrents that screen
the magnetic field produced by the field coil, effectively reducing the mutual
inductance $M$ between the field coil and the SQUID pickup loop.

In Fig. 1b we show the change in the measured mutual inductance $\delta
M/M_{FC}$ at a fixed temperature, as the susceptometer approaches the sample
surface. Here, we set $\delta M=0$ when the SQUID is far from the sample, and
$M_{FC}$ is the bare mutual inductance between the SQUID pickup loop and field
coil, which we measured to be
$330\text{\,}\mathrm{\Phi_{0}}\text{/}\mathrm{A}$ far from the sample surface.
With these choices, $\delta M/M_{FC}=-1$ would correspond to a complete
screening of the magnetic field in the plane of the pickup loop. The signals
that we observe correspond to a fraction of this value as expected for a thin
film. In the Pearl limit $\lambda>>d$ with $d$ the film thickness and
neglecting the finite thickness and detailed geometry of the pickup loop and
field coil, Kogan and Kirtley et al. derived a model for the dependence of
$\delta M/M_{FC}$ on the height $z$ above the sample:

$\frac{\delta
M(z,T)}{M_{FC}}=-\frac{ad}{2\lambda^{2}(T)}\left[1-\frac{2z}{\sqrt{a^{2}+4z^{2}}}\right].$
(1)

Figure 1: (a) Schematic of the SQUID pickup loop and field coil pair close to
the surface of an Sr2RuO4 thin film grown on an LSAT substrate. A current
$I_{FC}$ applied to the field coil with an inner diameter of
$8\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$ induces a screening current
$I_{SC}$ in the superconducting Sr2RuO4. The sample response to $I_{FC}$ is
detected via the SQUID pickup loop, with an inner diameter of
$1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$. (b) Change in mutual
inductance $\delta M$ between the SQUID and field coil as the voltage $V_{z}$
applied to the z-positioner is increased, bringing the susceptometer towards
the sample surface. The solid line is a fit to Eq. 1

.

In Fig. 1b we include a fit to Eq. 1 of the height dependence of $\delta
M/M_{FC}$. We constrain the fit parameters to specific values:
$a=$4\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$$,
$d=$30\text{\,}\mathrm{nm}\text{/}$$ and
$z_{0}=$4\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$$, where $z_{0}$ is
the height of the susceptometer above the sample surface when the corner of
the SQUID chip first makes mechanical contact with the sample. The estimated
value of $z_{0}$ is based on our alignment and the SQUID geometry. The fitting
procedure yields $\lambda=$350\text{\,}\mathrm{nm}\text{/}$$ and
$\alpha=$120\text{\,}\mathrm{nm}\text{/}\mathrm{V}$$. Here, $\alpha$
characterizes the change in SQUID height per volt applied to the piezo
positioner. These values agree reasonably well with previous measurements of
$\alpha$ for our scanner and estimates of the low-temperature value of
$\lambda$ for films with a $T_{c}$ of $\sim$900\text{\,}\mathrm{mK}\text{/}$$
(See Appendix A for details). The absolute value of $\lambda$ is affected by
uncertainty in the fit parameters, which are strongly correlated with
$\lambda$, making a meaningful determination of the absolute value of
$\lambda$ challenging (see Appendix A: Scanning SQUID Susceptometry for
details).

Fortunately, an important feature of Eq. 1 is that its sole temperature
dependence originates from $\lambda(T)$, which appears in a simple prefactor.
This implies that by measuring $M$ as a function of temperature at a constant
height above the sample, we can accurately extract relative changes in
$\lambda(T)$, even in the presence of uncertainties in the geometric factors.
To determine the temperature dependence of $\lambda$, we measure $M$ at a
fixed position as a function of temperature, while maintaining light
mechanical contact with the sample. Assuming that $M$ measured at
$50\text{\,}\mathrm{mK}\text{/}$ reasonably approximates
$\lambda_{0}=\lambda(T=0)$ at a given position, we can obtain
$\lambda(T)/\lambda_{0}-1=\delta\lambda(T)\lambda_{0}$ from
$\frac{M(z_{o},T)}{M(z_{o},T=0)}=\frac{\lambda_{0}^{2}}{\lambda^{2}(T)}$.

Figure 2: (a) Spatial map at $20\text{\,}\mathrm{mK}\text{/}$ of the local
diamagnetic response. The two dots indicate the locations where detailed
temperature dependent data were acquired. The left-most dot (blue) had a
higher local superconducting transition temperature
$T_{c}=$0.90\text{\,}\mathrm{K}\text{/}$$ than the right dot (purple), with
$T_{c}=$0.85\text{\,}\mathrm{K}\text{/}$$. (b)
$\delta\lambda\left(T\right)=\lambda(T)-\lambda(0)$, normalized by
$\lambda(0)$ at the positions indicated in (a). (c)
$\delta\lambda(T)/\lambda(0)$ plotted against $\left(T/T_{c}\right)^{2}$.

Fig. 2a shows an image of the local magnetic response acquired by scanning
several microns above the sample surface while maintaining the sample
temperature at $20\text{\,}\mathrm{mK}\text{/}$. The image reveals micrometer-
scale variations in the local magnetic susceptiblity, and accordingly
$\lambda$. Next, we conduct detailed temperature-dependent measurements of the
penetration depth at two positions on the sample surface marked in Fig. 2a. In
Fig. 2b, we show $\delta\lambda(T)/\lambda_{0}$ versus temperature for these
two positions. We find a stronger diamagnetic response at base temperature and
a higher local critical temperature, $T_{c}$ at Position 1 than at Position 2.
For both positions, we observe no significant difference between data
collected upon warming and cooling indicating that we sweep temperature
sufficiently slow for the sample to thermalize.

In Fig. 2c, we plot $\delta\lambda(T)/\lambda_{0}$ as a function of
$T^{2}/T_{c}^{2}$ for both positions. The data fall on a straight line,
indicating that for $T<0.4T_{c}$, $\lambda(T)$ exhibits a $T^{2}$ temperature
dependence. A temperature dependence of $\lambda(T)\sim T^{2}$ has been
observed previously in high-quality single crystals Sr2RuO4 using a tunnel
diode oscillator technique [15], and more recently using scanning SQUID
microscopy on single crystals under uniaxial strain [18]. The power-law
temperature dependence observed in our study and in the single crystal work
provides evidence for nodes in the superconducting gap function of Sr2RuO4.

For nodal superconducting order parameters in the clean limit, a linear
temperature dependence, $\delta\lambda(T)\sim T$ is expected. In Sr2RuO4
single crystals, the $\delta\lambda(T)\approx T^{2}$ temperature dependence
was attributed to non-local electrodynamics of the nodal quasiparticles [19,
15, 18]. Impurity scattering can also change the temperature dependence of
$\lambda(T)$ superconductors with nodes in the gap function. For example, for
an order parameter with $d_{x^{2}-y^{2}}$ symmetry, strong impurity scattering
is expected to give $\delta\lambda(T)\sim T^{2}$ [20]. In our thin-film
samples, $T_{c}\approx$1\text{\,}\mathrm{K}\text{/}$$ is suppressed compared
to the highest $T_{c}\approx$1.5\text{\,}\mathrm{K}\text{/}$$ observed in
single crystals, and it is plausible that impurity scattering plays a larger
role in determining the superconducting properties of the thin films than in
bulk crystals.

To assess if scattering can explain the temperature dependence of the
penetration depth, we examine the superfluid density across the full
temperature range. The normalized superfluid density is directly related to
the penetration depth through
${\rho_{s}(T)}/{\rho_{0}}={\lambda_{0}^{2}}/{\lambda^{2}(T)}$, where
$\rho_{0}$ is its zero temperature value. Figure 3 shows the normalized
superfluid density for both positions measured in Figure 2. Although we found
local variations in both the magnitude of $\rho_{s}$ and $T_{c}$, we find that
the normalized superfluid density collapses onto a single curve for both
positions. Additional measurements collected on other areas of the sample,
show the same temperature dependence of the normalized superfluid density (See
Supplementary Information for Details).

In general, the temperature dependence of $\rho_{s}$ depends among other
factors on the gap structure of the superconductor, the underlying band
structure and scattering in the sample. Nevertheless, it is useful to compare
the temperature dependence of $\rho_{s}$ to expectations from different
simplified phenomenological models. In Fig. 3 we directly compare our
measurements to calculations of the superfluid density for different
superconducting gap functions assuming a single circular Fermi surface in the
weak-coupling BCS limit (see Appendix: Calculation of the Superfluid Density
for details). We find that neither the fully gapped model nor the d-wave
models without scattering closely resemble our data.

To include the effect of scattering on the superfluid density, we estimate the
strength of pair-breaking scattering in our samples using the theory of
Abrikosov and Gor’kov [21], which has been used previously to successfully
describe the dependence of $T_{c}$ on disorder in Sr2RuO4 in both single
crystals and epitaxial thin films [14]. Within this theory, $T_{c}$ satisfies

$\ln\left(\frac{T_{c0}}{T_{c}}\right)=\Psi\left(\frac{1}{2}+\frac{\Gamma_{N}}{2\pi
T_{c}}\right)-\Psi\left(\frac{1}{2}\right),$ (2)

where $\Psi$ is the digamma function, $\Gamma_{N}=\hbar/2\tau k_{B}$
characterizes the strength of pair-breaking scattering with $\tau$ the
corresponding scattering time, $k_{B}$ the Boltzmann constant, $\hbar$ the
reduced Planck’s constant, and $T_{c0}$ is the zero-disorder limit of $T_{c}$.
Here we take $T_{c0}=$1.5\text{\,}\mathrm{K}\text{/}$$, the highest value of
$T_{c}$ observed in bulk crystals [22]. In the thin film grown on LSAT studied
here, we observe a suppressed $T_{c}\approx$0.9\text{\,}\mathrm{K}\text{/}$$,
corresponding to $\Gamma_{N}\approx$0.7\text{\,}\mathrm{K}\text{/}$$. Here we
assume that the Abrikosov-Gor’kov theory is applicable to both thin-film and
bulk samples [23], and that both types of sample have the same maximum
$T_{c0}$. We note that Sr2RuO4 thin films host defects that are not present in
bulk samples. For example, step edges in the substrate due to a slight
misorientation between the substrate surface and the (001) LSAT axis can cause
out-of-phase boundaries that extend through a significant fraction of the film
thickness [12, 24].

Using the estimated value for $\Gamma_{N}$, we calculate the normalized
superfluid density assuming a d-wave gap function with vertical line nodes, a
single circular Fermi surface and including the effects of impurity scattering
by following the approach of [25] (details are provided in the Appendix). The
result is included as line labeled “disordered nodal” in Fig. 3, and provides
much better agreement with our experimental data. We emphasize that this curve
is not a fit to the experimental data, but a model for a disordered d-wave
superconductor with the scattering rate estimated by the experimentally
measured $T_{c}$.

The agreement between this model and the data is remarkable given the
simplicity of the Fermi surface and gap structure used in the calculations.
Further refinement of this model by introducing multiband effects, the
experimentally determined Fermi surface and a realistic pairing potential may
improve the agreement between the model and data further.

Figure 3: Normalized superfluid density $\lambda^{2}(0)/\lambda^{2}(T)$
measured at the two positions indicated in Fig. 2a. Theoretical curves
assuming a weak-coupling BCS superconductor with a single cylindrical Fermi
surface are included for comparison.

We next discuss characterization of micrometer-scale rings fabricated from a
Sr2RuO4 thin film. Lithographically defined devices are enabled by thin film
samples of Sr2RuO4 and promising for future device-based tests of the
superconducting order parameter symmetry. We patterned an array of rings on a
second Sr2RuO4 thin film sample grown on an LSAT substrate using standard
photolithography and argon ion milling. We present characterization of rings
with an inner diameter of $4\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$
and an outer diameter of $10\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$.

Figure 4: (a) Magnetic susceptibility image of three superconducting Sr2RuO4
rings. (b) SQUID signal as a function of field coil current acquired at
several temperatures below $T_{c}$ with the susceptometer positioned above the
right-most ring in (a). The response of the ring is hysteretic with sharp
jumps indicating of a discrete change in the ring fluxoid number superimposed
on a diamagnetic background. (c) Size of the fluxoid jumps plotted against the
diamagnetic susceptibility extracted from (b). Symbols indicate the direction
of the field coil ramp for each jump. Color scale is the same as (b). As the
temperature decreases, both the fluxoid jump size and the strength of the ring
diamagnetic response increase.

In Fig. 4a we show a magnetic susceptibility image of three representative
Sr2RuO4 rings. To characterize the magnetic response of the rings in detail,
we position the susceptometer above the ring and apply a current oscillating
sinusoidally at a few Hz to the field coil while monitoring the ring response
with the SQUID. We average the signal from up to a few hundred oscillations of
the field coil current in order to achieve a high signal-to-noise ratio. A
background signal, measured with the SQUID retracted from the sample at each
temperature, is subtracted from the data. Fig. 4b shows the magnetic response
of the right-most ring in Fig. 4a at several temperatures below $T_{c}$. As
the field coil current is swept, we observe a linear response which is
periodically interrupted by abrupt jumps in the SQUID signal. At each
temperature measured, the response of the ring to the field coil current is
hysteretic. The size of the hysteresis loop, the slope of the ring response,
and the size of the jumps in the signal all decrease as the temperature of the
sample is raised towards $T_{c}$.

To understand the magnetic response of the rings, we turn to the following
London equation which determines the current distribution in the ring,

$\vec{j}=-\frac{1}{\mu_{0}\lambda^{2}}\left[\frac{\Phi_{0}}{2\pi}\left(\nabla\theta\right)+\vec{A}\right],$

where $\vec{j}$ is the supercurrent density, $\mu_{0}$ is the vacuum
permeability, $\Phi_{0}=\frac{h}{2e}$ is the superconducting flux quantum,
$\theta$ is the order parameter phase and $\vec{A}$ is the magnetic vector
potential. Integrating over the film thickness and introducing cylindrical
coordinates $r$ and $\phi$ with respect to the center of the ring gives,

$g_{\phi}(r)=\frac{d}{\mu_{o}\lambda^{2}}\left[\frac{\Phi_{o}}{2\pi}\frac{N}{r}-A_{\phi}(r)\right],$
(3)

where we have enforced single-valuedness of the superconducting order
parameter with the integer $N$ the number of vortices trapped inside the ring,
and $g_{\phi}(r)$ the sheet current density flowing along the azimuth $\phi$
of the ring. $A_{\phi}$ is the azimuthal component of the magnetic vector
potential.

This model captures the main qualitative features of our data. We interpret
the sudden jumps in Fig. 4b as changes in the vortex number $N$ via the
transport of a superconducting vortex accross the wall of the ring. Changes in
$N$ generate discontinuous changes in the current flowing through the ring,
which produce jumps in the flux coupled into the SQUID pickup loop. Away from
these transitions, the flux threading the ring changes linearly with the field
coil current corresponding to a linear change in $A_{\phi}$ which increases
the diamagnetic screening current in the ring and consequently the flux
coupled into the SQUID pickup loop.

A quantitative prediction of the jump height requires a self-consistent
solution of $g_{\rho}(r)$ in the presence of the $A_{\phi}$ generated by the
field coil excitaion. For our samples and susceptometer, where the field coil
dimensions are comparable to the ring dimensions, these calculations must be
carried out numerically [26]. Eq. 3 shows that both the size of the jumps when
$N$ changes by one and the strength of the diamagnetic screening at a given
temperature $T$ are proportional to $1/\lambda^{2}$. In Fig. 4c we show the
size of the jumps at each $T$ plotted against the slope of the response away
from the jumps, which we attribute to diamagnetic screening. Consistent with
Eq. 3, we find that the jump height and the strength of diamagnetic screening
increase monotonically and proportional to each other as the temperature is
lowered and $1/\lambda^{2}$ increases. We find that the number of vortex
transitions goes down as the temperature is lowered, consistent with vortex
motion accross the ring wall. As the superfluid density in the ring increases,
the energetic barrier for a vortex tunneling accross the ring increases. The
field coil current required to induce a transition therefore increases and the
hysteresis in the ring signal becomes more pronounced at lower temperatures.
We investigate the temperature range between
$0.650\text{\,}\mathrm{K}\text{/}$ and $0.550\text{\,}\mathrm{K}\text{/}$,
near the superconducting transition of the rings
$T_{c}=$0.725\text{\,}\mathrm{K}\text{/}$$, where fluxoid transitions are
readily accessed with field coil currents of a few mA.

In summary, we find that the penetration depth in thin-film Sr2RuO4 shows a
$\lambda(T)\propto T^{2}$ dependence at low temperatures similar to reports in
single crystal bulk samples. This suggests, as is expected, that the gap
structure in our thin films is comparable to single crystal samples and
features nodes. Analysis of the full temperature dependence of the superfluid
density shows that it can be explained by a nodal superconducting gap combined
with scattering in our samples. While this does not exclude non-local effects
in the Meissner screening, scattering is in our case the simplest explanation.
In addition, we provide the first magnetic characterization of mesoscopic
superconducting rings lithographically patterned from a Sr2RuO4 thin-film.
Looking forward, our work shows that the delicate superconducting state of
Sr2RuO4 can be preserved in high-quality thin film devices, which is an
important first step towards device-based tests of the order parameter
symmetry. Furthermore, the effect of biaxial strain on the superconducting
state in Sr2RuO4 may be studied by choosng suitable substrates. Uniaxial
strain has proven to be an effective tuning parameter for superconductivity in
bulk Sr2RuO4 [7, 8]. Biaxially strained thin-film samples [27] may exhibit a
qualitatively different ground state than their uniaxially strained
counterparts [28]. The experimental approach demonstrated in this work will
provide direct access to new superconducting or magnetic states that emerge in
such samples.

###### Acknowledgements.

Scanning SQUID measurements were primarily supported through the Cornell
University Materials Research Science and Engineering Center DMR-1719875. Thin
film growth was supported by the National Science Foundation Platform for the
Accelerated Realization, Analysis, and Discovery of Interface Materials
(PARADIM) under Cooperative Agreement No. DMR-2039380. In addition, this
research was funded in part by the Gordon and Betty Moore Foundation’s EPiQS
Initiative through Grant Nos. GBMF3850 and GBMF9073 to Cornell University.
Sample preparation was facilitated in part by the Cornell NanoScale Facility,
a member of the National Nanotechnology Coordinated Infrastructure (NNCI),
which is supported by the National Science Foundation (Grant No.
NNCI-2025233). K.M.S. acknowledges the support of Air Force Office of
Scientific Research Grant No. FA9550-21-1-0168 and National Science Foundation
DMR-2104427.

##  Appendix A Scanning SQUID Susecptometry

Our SQUID susceptometer features a counter-wound field coil pickup loop pair
that allows us to null the SQUID response to the field coil excitation. To
perform local magnetic susceptibility measurements, the SQUID is retracted as
far as possible from the sample surface (about
$30\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$) and the susceptometer is
balanced by adjusting the current flowing in the two arms of the field coil
until the flux coupled into the SQUID pickup loops is minimized. Any remaining
SQUID signal after balancing is recorded as an offset signal which represents
the response of the SQUID in to the field coil in the absence of the sample.
This offset is then subtracted from the height sweep data used to characterize
the susceptometer geometry in Fig. 1b and from all height sweep data used to
extract the penetration depth and superfluid density (Figs. 2-3). No offset
was subtracted from imaging measurements of the local susceptibility (Fig 2a,
Fig 4a) as these images are used to characterize the spatial variations in the
diamagentic response of the sample.

To avoid spurious signals from vortex motion during the penetration depth
measurements, the susceptibility data were collected conducted under near zero
field conditions. A superconducting solenoid was used to adjust the out-of-
plane magnetic field threading the sample until no vortices were observed
within the $$140\text{\,}\times$$140\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{/}$^{2}$ field of view after warming and cooling the sample and SQUID
through their respective critical temperatures. This procedure places a rough
upper bound of
$\frac{\Phi_{0}}{$140\text{\,}\times$$140\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{/}$^{2}}\approx$100\text{\,}\mathrm{nT}\text{/}$$ on the magnetic flux
density in the sample environment. Before and after conducting the temperature
sweeps used to acquire the magnetic susceptibility data, the SQUID was used to
acquire an image confirming that vortices were absent from the measurement
area.

In the main text, we utilized an analytical expression (Eq. 1) relating the
magnetic penetration depth to the geometry of our magnetic susceptibility
measurement. Although the data is well-fit by this model, it is strictly only
valid deep in the Pearl limit where, $\lambda>>d$. Muon spin relaxation
measurements of the zero-temperature penetration depth in Sr2RuO4 have found
that at the lowest temperatures,
$\lambda(T)\approx$190\text{\,}\mathrm{nm}\text{/}$$ [29], which is only a few
times our film thickness of $35\text{\,}\mathrm{nm}\text{/}$. In the thin-film
samples that we measure here, $T_{c}$ is substantially suppressed from the
values observed in the best single crystals. Using calculations of the zero
temperature superfluid density in the presence of disorder, described below,
we estimate that $\lambda(0)$ in our samples
($T_{c}\approx$900\text{\,}\mathrm{mK}\text{/}$$,
$\Gamma_{N}\approx$0.7\text{\,}\mathrm{K}\text{/}$$) is a factor of $\sim 1.7$
larger than $\lambda(0)$ in the best single crystals. We therefore estimate
that at all temperatures measured, $\lambda(T)\gtrsim 10d$.

To extract the temperature dependence of $\lambda(T)$ from the magnetic
susceptibility measurements, we observe that the expression Eq. 1 may be
separated into the product of a temperature-dependent part and a temperature-
independent part that depends only on the geometric configuration of the SQUID
and sample,

$\frac{M(z,T)}{M_{FC}}=-\frac{ad}{2\lambda^{2}(T)}\left[1-\frac{2z}{\sqrt{a^{2}+4z^{2}}}\right]=\frac{A_{geo}(z)}{\lambda^{2}(T)}.$
(4)

Where $A_{geo}(z)=-\frac{ad}{2}\left[1-\frac{2z}{\sqrt{a^{2}+4z^{2}}}\right]$,
encodes the sample-susceptometer geometry. From this expression, the
temperature dependence of the penetration depth,
$\delta\lambda(T)/\lambda_{o}$ may be written,

$\frac{\delta\lambda(T)}{\lambda_{o}}=\frac{\lambda(T)-\lambda_{o}}{\lambda_{o}}=\sqrt{\frac{M(z_{o},T=0)}{M(z_{o},T)}}-1,$
(5)

allowing $\delta\lambda(T)/\lambda_{o}$ to be directly calculated from
measurements of $M(z_{o},T)$.

##  Appendix B Superfluid Density Calculations

In this work we calculate the superfluid density in the weak-coupling limit of
BCS theory.

To calculate the superfluid density in the presence of disorder, we turn to
the self-consistent t-matrix approximation [20, 25]. Within the approximation,
impurities are treated as isotropic point scatters. For simplicity, we again
model the system as a single circular Fermi surface, and choose the simplest
separable pairing potential $V_{\textbf{k},\textbf{k'}}$ with a d-wave form
factor,

$V_{\textbf{k},\textbf{k'}}=V_{0}\Omega_{\textbf{k}}\Omega_{\textbf{k'}}.$ (6)

Where,

$\Omega_{\textbf{k}}\propto\cos{k_{x}a}-\cos{k_{y}a},$ (7)

with $a$ the lattice spacing. Given this Fermi surface and pairing potential,
we find solutions to the gap equation,

$\Delta_{\textbf{k}}=2\pi
T\sum_{\omega_{n}>0}^{\omega_{0}}\left<V_{\textbf{k}\textbf{k}^{\prime}}\frac{\Delta_{\textbf{k}^{\prime}}}{\left(\tilde{\omega}_{n}^{2}+\Delta_{\textbf{k}^{\prime}}^{2}\right)}\right>_{FS}.$
(8)

Here, $\Delta_{\textbf{k}}=\psi(T)\Omega_{\textbf{k}}$ is the superconducting
gap, with $\psi(T)$ the temperature-dependent gap amplitude. $\omega_{n}=2\pi
T\left(n+\frac{1}{2}\right)$ are the fermionic Matsubara frequencies and
$\left<\cdot\cdot\cdot\right>_{FS}$ denotes an average over the Fermi surface.
Adding disorder to the system renormalizes the Matsubara frequencies:

$\tilde{\omega}_{n}=\omega_{n}+\pi\Gamma\frac{\left<N_{\textbf{k}}(\tilde{\omega}_{n})\right>_{FS}}{c^{2}+\left<N_{\textbf{k}}(\tilde{\omega}_{n})\right>_{FS}}.$
(9)

Where $c$ is the cotangent of the scattering phase shift and $\Gamma$
parameterizes the density of scattering sites. The limits $c<<1$ and $c>>1$
correspond to strong (unitary) and weak (Born) scattering respectively. For
this work, we observe that impurities appear to act as strong pair-breaking
scatterers in Sr2RuO4 [14, 23] and take $c=0$. Within weak-coupling BCS
theory, the gap equation may be solved by finding the $\tilde{\omega}_{n}$ and
$\psi(T)$ that satisfy,

$\ln{\left(\frac{T_{c0}}{T}\right)}=2\pi
T\sum_{\omega_{n}>0}^{\infty}\left(\frac{1}{\omega_{n}}-\left<\frac{\Omega_{k}^{2}}{\left(\tilde{\omega}_{n}^{2}+\psi^{2}\Omega_{k}^{2}\right)^{1/2}}\right>_{FS}\right).$
(10)

Equations 10 and 9 may be solved self-consistently at a range of temperatures
to obtain $\tilde{\omega}_{n}$ and $\psi(T)$. Once $\tilde{\omega}_{n}$ and
$\psi(T)$ are known, the superfluid density in the presence of disorder may be
calculated,

$\frac{\rho_{s}(T)}{\rho_{s00}}=2\pi
T\sum_{\omega_{n}>0}^{\infty}\left<\left<\frac{\Delta_{\textbf{k}}^{2}}{\left(\tilde{\omega}_{n}^{2}+\Delta_{\textbf{k}}^{2}\right)^{3/2}}\right>\right>_{FS}.$
(11)

$\left<\left<\cdot\cdot\cdot\right>\right>_{FS}$ denotes a velocity-weighted
average over the Fermi surface, and in the case of a circular Fermi surface,

$\left<\left<A(\phi)\right>\right>_{FS}=\frac{\int_{0}^{2\pi}A(\phi)v_{F,x}^{2}}{\int_{0}^{2\pi}v_{F,x}^{2}}.$
(12)

Examples of calculated $\psi(T)$ and $\rho_{s}(T)$ under different impurity
scattering strengths $\Gamma_{N}$ are presented in the Supplementary
Information.

##  Appendix C Sample Growth

The Sr2RuO4 thin film was grown in a Veeco Gen10 molecular-beam epitaxy (MBE)
system on a
$(\mathrm{LaAlO_{3}})_{\mathrm{0.29}}-(\mathrm{SrAl}_{0.5}\mathrm{Ta}_{0.5}\mathrm{O}_{3})_{0.71}$
substrate from CrysTec GmbH. The substrate used for the growth was screened to
have a miscut of less than
$0.05\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$, which is important to
reduce the formation of out-of-phase boundaries. The films were grown at a
substrate temperate of $810\text{\,}\mathrm{\SIUnitSymbolCelsius}\text{/}$ as
measured using an optical pyrometer operating at
$1550\text{\,}\mathrm{nm}\text{/}$. Elemental strontium (99.99% purity) and
elemental ruthenium (99.99% purity) evaporated from a low-temperature effusion
cell and a Telemark electron beam evaporator, respectively, were used for
growing the Sr2RuO4 film. The films were grown with a strontium flux of
$2.6\text{\times}{10}^{13}\text{\,}\mathrm{atoms}\text{/}\text{(}{\mathrm{cm}}^{2}\text{\,}\mathrm{s}\text{)}$
and a ruthenium flux of
$1.8\text{\times}{10}^{13}\text{\,}\mathrm{atoms}\text{/}\text{(}{\mathrm{cm}}^{2}\text{\,}\mathrm{s}\text{)}$
in a background of distilled ozone ($\sim$80% $\mathrm{O_{3}}$ \+ 20%
$\mathrm{O_{2}}$ made from oxygen gas with 99.994% purity.) The background
oxidant pressure during growth was
$3\text{\times}{10}^{-6}\text{\,}\mathrm{Torr}\text{/}$. At the end of the
growth the strontium and ruthenium shutters were closed simultaneously, and
the sample was cooled down to below
$250\text{\,}\mathrm{\SIUnitSymbolCelsius}\text{/}$ in a background pressure
of distilled ozone of $1\text{\times}{10}^{-6}\text{\,}\mathrm{Torr}\text{/}$.
Further details of the adsorption-controlled growth conditions for the growth
of Sr2RuO4 thin films by MBE can be found elsewhere [12].

##  Appendix D Fabrication Procedure

The Sr2RuO4 films grown by MBE were subsequently patterned into rigs for
scanning SQUID microscopy using standard photolithography, sputter deposition,
and ion milling techniques. First, a Pt meander structure used for navigation
on the sample was defined using photolithography. Next,
$25\text{\,}\mathrm{nm}\text{/}$ of platinum, with
$5\text{\,}\mathrm{nm}\text{/}$ titanium adhesion layer, was sputtered onto
the Sr2RuO4 film with an AJA sputtering tool, followed by a standard lift-off
processes. A second photolithography step was used to define the ring
geometry, followed by ion milling with an AJA ion mill to remove the excess
Sr2RuO4 film.

## References

  * Maeno  _et al._ [1994] Y. Maeno, H. Hashimoto, K. Yoshida, S. Nishizaki, T. Fujita, J. Bednorz, and F. Lichtenberg, Superconductivity in a layered perovskite without copper, nature 372, 532 (1994).
  * Ishida  _et al._ [1998] K. Ishida, H. Mukuda, Y. Kitaoka, K. Asayama, Z. Mao, Y. Mori, and Y. Maeno, Spin-triplet superconductivity in Sr2RuO4 identified by 17O Knight shift, Nature 396, 658 (1998).
  * Luke  _et al._ [1998] G. M. Luke, Y. Fudamoto, K. Kojima, M. Larkin, J. Merrin, B. Nachumi, Y. Uemura, Y. Maeno, Z. Mao, Y. Mori, _et al._ , Time-reversal symmetry-breaking superconductivity in Sr2RuO4, Nature 394, 558 (1998).
  * Nelson  _et al._ [2004] K. Nelson, Z. Mao, Y. Maeno, and Y. Liu, Odd-parity superconductivity in Sr2RuO4, Science 306, 1151 (2004).
  * Xia  _et al._ [2006] J. Xia, Y. Maeno, P. T. Beyersdorf, M. Fejer, and A. Kapitulnik, High resolution polar kerr effect measurements of Sr2RuO4: Evidence for broken time-reversal symmetry in the superconducting state, Physical review letters 97, 167002 (2006).
  * Kidwingira  _et al._ [2006] F. Kidwingira, J. Strand, D. Van Harlingen, and Y. Maeno, Dynamical superconducting order parameter domains in Sr2RuO4, Science 314, 1267 (2006).
  * Hicks  _et al._ [2014] C. W. Hicks, D. O. Brodsky, E. A. Yelland, A. S. Gibbs, J. A. Bruin, M. E. Barber, S. D. Edkins, K. Nishimura, S. Yonezawa, Y. Maeno, _et al._ , Strong increase of $t_{c}$ of Sr2RuO4 under both tensile and compressive strain, Science 344, 283 (2014).
  * Steppke  _et al._ [2017] A. Steppke, L. Zhao, M. E. Barber, T. Scaffidi, F. Jerzembeck, H. Rosner, A. S. Gibbs, Y. Maeno, S. H. Simon, A. P. Mackenzie, _et al._ , Strong peak in $t_{c}$ of Sr2RuO4 under uniaxial pressure, Science 355, eaaf9398 (2017).
  * Pustogow  _et al._ [2019] A. Pustogow, Y. Luo, A. Chronister, Y.-S. Su, D. Sokolov, F. Jerzembeck, A. P. Mackenzie, C. W. Hicks, N. Kikugawa, S. Raghu, _et al._ , Constraints on the superconducting order parameter in Sr2RuO4 from oxygen-17 nuclear magnetic resonance, Nature 574, 72 (2019).
  * Ishida  _et al._ [2020] K. Ishida, M. Manago, K. Kinjo, and Y. Maeno, Reduction of the 17O Knight shift in the superconducting state and the heat-up effect by nmr pulses on Sr2RuO4, Journal of the Physical Society of Japan 89, 034712 (2020).
  * Li  _et al._ [2021] Y.-S. Li, N. Kikugawa, D. A. Sokolov, F. Jerzembeck, A. S. Gibbs, Y. Maeno, C. W. Hicks, J. Schmalian, M. Nicklas, and A. P. Mackenzie, High-sensitivity heat-capacity measurements on Sr2RuO4 under uniaxial pressure, Proceedings of the National Academy of Sciences 118, e2020492118 (2021).
  * Nair  _et al._ [2018] H. P. Nair, J. P. Ruf, N. J. Schreiber, L. Miao, M. L. Grandon, D. J. Baek, B. H. Goodge, J. P. Ruff, L. F. Kourkoutis, K. M. Shen, _et al._ , Demystifying the growth of superconducting Sr2RuO4 thin films, APL Materials 6, 101108 (2018).
  * Tsuei and Kirtley [2000] C. Tsuei and J. Kirtley, Pairing symmetry in cuprate superconductors, Reviews of Modern Physics 72, 969 (2000).
  * Mackenzie  _et al._ [1998] A. Mackenzie, R. Haselwimmer, A. Tyler, G. Lonzarich, Y. Mori, S. Nishizaki, and Y. Maeno, Extremely strong dependence of superconductivity on disorder in Sr2RuO4, Physical review letters 80, 161 (1998).
  * Bonalde  _et al._ [2000] I. Bonalde, B. D. Yanoff, M. Salamon, D. Van Harlingen, E. Chia, Z. Mao, and Y. Maeno, Temperature dependence of the penetration depth in Sr2RuO4: evidence for nodes in the gap function, Physical review letters 85, 4775 (2000).
  * Gardner  _et al._ [2001] B. W. Gardner, J. C. Wynn, P. G. Björnsson, E. W. Straver, K. A. Moler, J. R. Kirtley, and M. B. Ketchen, Scanning superconducting quantum interference device susceptometry, Review of Scientific Instruments 72, 2361 (2001).
  * Huber  _et al._ [2008] M. E. Huber, N. C. Koshnick, H. Bluhm, L. J. Archuleta, T. Azua, P. G. Björnsson, B. W. Gardner, S. T. Halloran, E. A. Lucero, and K. A. Moler, Gradiometric micro-SQUID susceptometer for scanning measurements of mesoscopic samples, Review of Scientific Instruments 79, 053704 (2008).
  * Mueller  _et al._ [2023] E. Mueller, Y. Iguchi, F. Jerzembeck, J. O. Rodriguez, M. Romanelli, E. Abarca-Morales, A. Markou, N. Kikugawa, D. A. Sokolov, G. Oh, _et al._ , Superconducting penetration depth through a van hove singularity: Sr2RuO4 under uniaxial stress, arXiv preprint arXiv:2312.05130 (2023).
  * Kosztin and Leggett [1997] I. Kosztin and A. J. Leggett, Nonlocal effects on the magnetic penetration depth in d-wave superconductors, Physical review letters 79, 135 (1997).
  * Hirschfeld and Goldenfeld [1993] P. J. Hirschfeld and N. Goldenfeld, Effect of strong scattering on the low-temperature penetration depth of a d-wave superconductor, Physical Review B 48, 4219 (1993).
  * Abrikosov and Gor’kov [1960] A. A. Abrikosov and L. P. Gor’kov, Contribution to the theory of superconducting alloys with paramagnetic impurities, Zhur. Eksptl’. i Teoret. Fiz. 39 (1960).
  * Akima  _et al._ [1999] T. Akima, S. NishiZaki, and Y. Maeno, Intrinsic superconducting parameters of Sr2RuO4, Journal of the Physical Society of Japan 68, 694 (1999).
  * Ruf  _et al._ [2024] J. P. Ruf, H. M. Noad, R. Grasset, L. Miao, E. Zhakina, P. H. McGuinness, H. P. Nair, N. J. Schreiber, N. Kikugawa, D. Sokolov, _et al._ , Controllable suppression of the unconventional superconductivity in bulk and thin-film Sr2RuO4 via high-energy electron irradiation, arXiv preprint arXiv:2402.19454 (2024).
  * Fang  _et al._ [2021] Y. Fang, H. P. Nair, L. Miao, B. Goodge, N. J. Schreiber, J. P. Ruf, L. F. Kourkoutis, K. M. Shen, D. G. Schlom, and B. J. Ramshaw, Quantum oscillations and quasiparticle properties of thin film ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$, Phys. Rev. B 104, 045152 (2021).
  * Lee-Hone  _et al._ [2017] N. Lee-Hone, J. Dodge, and D. Broun, Disorder and superfluid density in overdoped cuprate superconductors, Physical Review B 96, 024501 (2017).
  * Kirtley  _et al._ [2003] J. Kirtley, C. Tsuei, V. Kogan, J. Clem, H. Raffy, and Z. Li, Fluxoid dynamics in superconducting thin film rings, Physical Review B 68, 214505 (2003).
  * Burganov  _et al._ [2016] B. Burganov, C. Adamo, A. Mulder, M. Uchida, P. King, J. Harter, D. Shai, A. Gibbs, A. Mackenzie, R. Uecker, _et al._ , Strain control of fermiology and many-body interactions in two-dimensional ruthenates, Physical review letters 116, 197003 (2016).
  * Liu  _et al._ [2018] Y.-C. Liu, W.-S. Wang, F.-C. Zhang, and Q.-H. Wang, Superconductivity in Sr2RuO4 thin films under biaxial strain, Physical Review B 97, 224522 (2018).
  * Riseman  _et al._ [1998] T. Riseman, P. Kealey, E. Forgan, A. Mackenzie, L. Galvin, A. Tyler, S. Lee, C. Ager, D. M. Paul, C. Aegerter, _et al._ , Observation of a square flux-line lattice in the unconventional superconductor Sr2RuO4, Nature 396, 242 (1998).
  * Low  _et al._ [2021] D. Low, G. Ferguson, A. Jarjour, B. T. Schaefer, M. D. Bachmann, P. J. Moll, and K. C. Nowack, Scanning squid microscopy in a cryogen-free dilution refrigerator, Review of Scientific Instruments 92 (2021).
  * Huttema  _et al._ [2009] W. Huttema, J. Bobowski, P. Turner, R. Liang, W. Hardy, D. Bonn, and D. Broun, Stability of nodal quasiparticles in underdoped YBa2Cu3O6+y probed by penetration depth and microwave spectroscopy, Physical Review B 80, 104509 (2009).


# The DKU-DukeECE-Lenovo System for the Diarization Task
of the 2021 VoxCeleb Speaker Recognition Challenge

###### Abstract

This report describes the submission of the DKU-DukeECE-Lenovo team to the
VoxCeleb Speaker Recognition Challenge (VoxSRC) 2021 track 4. Our system
includes a voice activity detection (VAD) model, a speaker embedding model,
two clustering-based speaker diarization systems with different similarity
measurements, two different overlapped speech detection (OSD) models, and a
target-speaker voice activity detection (TS-VAD) model. Our final submission,
consisting of 5 independent systems, achieves a DER of 5.07% on the challenge
test set.

Index Terms: Speaker Diarization, Target-Speaker Voice Activity Detection

##  1 Introduction

As the speaker embedding becomes more and more robust, the conventional
diarization system also achieves good performance since the speaker confusion
has been significantly reduced. To further improve the performance and reduce
diarization error rate (DER), many researches focus on overlapped speech
detection (OSD) to reduce the missed speaker error, including speech
separation [1], overlap detection [2], end-to-end neural speaker diarization
(EEND) [3] and target-speaker voice activity detection (TS-VAD) [4].

We also explore many different OSD models in this challenge. First, a ResNet-
based model with LSTM back-end is employed for overlap detection, where the
output is 1 for overlapped speech and 0 otherwise. Second, an x-vector- and
ResNet-based TS-VAD model is used to refine the output from the conventional
diarization systems. Finally, we propose a 2-speaker TS-VAD model for overlap
detection, where a pair of speaker embeddings are fed to the TS-VAD model, and
the overlapped speech regions between these two speakers are detected.
Compared with the original TS-VAD, this method is not restricted to the number
of speakers.

##  2 Dataset Description

Since our TS-VAD model takes a long time for inference, we only use the last
46 recordings in the test dataset as our validation dataset, referred to as
VAL46. The development dataset with the remaining data in the test dataset is
used as our development dataset, which contains 402 recordings and is referred
to as DEV402. The detailed dataset used in this challenge for each model
include:

  * •

Voice activity detection (VAD): AMI [5], ICSI [6], ISL (LDC2004S05), NIST
(LDC2004S09), SPINE1&2 (LDC2000S87, LDC2000S96, LDC2001S04, LDC2001S06,
LDC2001S08), AISHELL-4 [7], DIHARD II [8] and DIHARD III [9] are the mixed
training set. DEV402 and VAL46 is used for fine-tuning and validation,
respectively.

  * •

Speaker embedding: We use Voxceleb 1 & 2 [10] as the training set.

  * •

Agglomerative hierarchical clustering (AHC): We directly tune the parameters
on DEV402.

  * •

LSTM-based similarity measurement with spectral clustering: We use the same
dataset as VAD does.

  * •

Overlap detection: We use DEV402 for training and VAL46 for validation.

  * •

2-speaker TS-VAD & TS-VAD: These models are first trained on the data
simulated by Librispeech [11]. Then we transfer the model to VoxConverse [12]
dataset with the data simulated by DEV402. Finally, we fine-tune the model on
DEV402 and validate it on VAL46.

  * •

Data augmentation: We perform data augmentation with MUSAN [13] and RIRs [14]
corpus.

Figure 1: The data simulation strategy for TS-VAD.

##  3 Detailed Model Configuration

The model architectures of VAD, overlap detection, and TS-VAD are very
similar, including a ResNet34 [15], a statistical pooling layer, two BiLSTM
[16] layers and two fully-connected layers with sigmoid function. For
ResNet34, the channel width is {32, 64, 128, 256} with kernel size of 3. Each
BiLSTM layer contains 256 units per direction with a dropout rate of 0.1. The
two fully connected layers contain 128 and 1 unit, respectively.

###  3.1 VAD

We use a ResNet34 as the front-end model to extract the frame-level feature
map. Next, a statistical pooling layer is employed on the feature map every
$S$ frames. Finally, two BiLSTM layers and two fully-connected layers with a
sigmoid function produce the posterior probability of speech. In our
experiments, the input is 8s chunked wav, and the acoustic feature is 80-dim
log Mel-filterbank energies with a frame length of 25ms and a frameshift of
10ms.

The model is first trained on the mixed training set for 100 epochs with a
learning rate of 0.0001 and then fine-tuned on DEV402 for 100 epochs with a
learning rate of 0.00001 with binary cross-entropy (BCE) loss and Adam
optimizer. We train four models with different $S=\\{1,2,4,8\\}$ to obtain the
VAD outputs with different resolutions. Finally, we average the outputs from
these four models in frame level. The threshold for speech decision is set to
0.6. Table 1 shows the false alarm (FA) and miss detection (MISS) on VAL46.

Table 1: The false alarm (FA), miss detection (MISS) and accuracy of the VAD model | FA [%] | MISS [%] |  Accuracy [%]
---|---|---|---
S=1 | 2.57 | 1.56 | 96.33
S=2 | 2.97 | 1.30 | 96.17
S=4 | 3.22 | 1.32 | 95.93
S=8 | 3.40 | 1.59 | 95.54
Fusion | 2.37 | 1.66 | 96.43

###  3.2 Speaker Embedding

The ResNet34 structure is employed as the front-end pattern extractor, which
learns a frame-level representation from the input acoustic feature. A global
statistic pooling (GSP) layer projects the variable length input to the fixed-
length vector. Next, a 128-dim fully connected layer is adopted as the speaker
embedding layer. The ArcFace [17] (s=32,m=0.2) is used as a classifier. The
detailed configuration of the neural network is the same as [18]. Table 2
shows the EER on the Voxceleb1 test set. The input is 2$\sim$4s chunked wav,
and the acoustic feature is 80-dim log Mel-filterbank energies with a frame
length of 25ms and a frameshift of 10ms.

Table 2: The EER of the speaker embedding model Training data | EER (%)
---|---
Voxceleb 2 | 1.23

###  3.3 Clustering-based System

####  3.3.1 AHC

We use a similar AHC method as [1] without speech separation. First, speaker
embeddings are extracted from the uniformly segmented speech with a length of
1.28s and shift of 0.32s, and two consecutive segments are merged into a
longer segment if the distance is greater than a segment threshold. The
pairwise similarity is measured by cosine distance. Next, we perform a plain
AHC on the similarity matrix with a relatively high stop threshold to obtain
the clusters with high confidence. These clusters are split into “long
clusters” and “short clusters” by the total duration in each cluster, and the
central embedding of each cluster is the mean of all speaker embeddings in the
cluster. Later, each short cluster is assigned to the closest long cluster by
the cosine distance of central embedding. Finally, if a short cluster is too
different from all long clusters, which means that the distance between them
is lower than a speaker threshold, we treat it as a new speaker.

All of these parameters are directly tuned on the DEV402 by grid search. In
our experiments, the segment threshold is 0.54, the stop threshold is 0.62,
the duration for classifying long and short clusters is 6s, and the speaker
threshold is 0.2.

####  3.3.2 LSTM-based Similarity Measurement with Spectral Clustering

We use the same LSTM-based model as [19]. The model consists of two BiLSTM and
two fully connected layers with a sigmoid function. Speaker embeddings are
also extracted from the uniformly segmented speech with a length of 1.28s and
shift of 0.64s. Next, the speaker embedding sequence
$[\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{n}]$ is concatenated with
repeated $\mathbf{x}_{i}$ as the input and produce the i-th row of the
affinity matrix $\mathbf{S}$:

$\mathbf{S}_{i}=[\mathbf{S}_{i,1},\mathbf{S}_{i,2},...,\mathbf{S}_{i,n}]=f(\left[\begin{matrix}\mathbf{x}_{i}\\\
\mathbf{x}_{1}\end{matrix}\right],\left[\begin{matrix}\mathbf{x}_{i}\\\
\mathbf{x}_{2}\end{matrix}\right],...,\left[\begin{matrix}\mathbf{x}_{i}\\\
\mathbf{x}_{n}\end{matrix}\right]),$ (1)

where $f$ is the LSTM-based neural network, n is set to 64 in our experiments.
More details can be found in [20].

The model is trained on the mixed training set for 100 epochs and fine-tuned
on DEV402 for 100 epochs. The model is optimized with BCE loss and Adam
optimizer with a learning rate of 0.001. After obtaining the affinity matrix
$\mathbf{S}$, we employ spectral clustering (SC) to get the final diarization
results.

###  3.4 Overlap Detection

The overlap detection model and training process are the same as that of the
VAD model except for the training data, labels, and resolutions. We train the
overlap detection model on DEV402, and we only average the outputs from two
models with $S=\\{1,2\\}$ since the overlapped regions are much shorter than
the speech regions in the VAD task. The label is 1 for overlapped speech and 0
otherwise. After an overlapped region is detected, we replace the label with
two closest speakers near this region. The threshold for overlap decision is
set to 0.8. The input is 8s chunked wav, and the acoustic feature is 80-dim
log Mel-filterbank energies with a frame length of 25ms and a frameshift of
10ms.

Table 3:  The performance of different speaker diarization systems. ∗Since the VAD model used in the 1st submissions is not well trained, the improvement of DER compared with other systems may come from both overlap detection (OD) and VAD on the test set.  Model | Submission | VAL46 | VoxSRC21 test set |
---|---|---|---|---
| DER[%] | JER[%] | DER[%] | JER[%] |
Baseline | - | - | - | 17.99 | 38.72 |
AHC + Oracle VAD | - | 2.61 | 21.93 | - | - |
LSTM-SC + Oracle VAD | - | 3.16 | 28.04 | - | - |
1 AHC |  1∗ | 4.42 | 26.42 |  6.77∗ |  27.66∗ |
2 LSTM-SC | - | 4.97 | 31.04 | - | - |
3 AHC + OD | - | 3.98 | 26.08 | - | - |
4 LSTM-SC + OD | - | 4.58 | 30.70 | - | - |
5 AHC + 2-spk TS-VAD as OD | 4 | 3.96 | 25.82 | 5.45 | 27.55 |
6 LSTM-SC + 2-spk TS-VAD as OD | - | 4.56 | 30.51 | - | - |
7 System 3 + TS-VAD | - | 4.53 | 28.39 | - | - |
8 System 5 + TS-VAD | - | 4.51 | 28.33 | - | - |
Fusion (3+4+7) | 2 | 3.93 | 25.68 | 5.36 | 29.10 |
Fusion (3+4+5+7) | 3 | 4.02 | 27.11 | 5.82 | 29.78 |
Fusion (3+4+5+6+8) | 5 | 4.10 | 30.97 | 5.07 | 29.16 |

###  3.5 TS-VAD

####  3.5.1 Data Simulation

We simulate two 900-hour pre-training datasets. One is simulated from
LibriSpeech, and another is simulated from DEV402. To obtain a simulated
dataset that is similar to the VoxConverse dataset, we first extract the label
from the DEV402, which is a matrix of size $T\times N$, where $T$ is the
number of frames and $N$ is the number of speakers. Then, for the label of
each speaker, we fill the active regions with a single speaker’s speech from
LibriSpeech or DEV402. Finally, we take the sum of all speeches as one
simulated recording. Note that the non-speech regions in the labels are first
removed. Figure 1 shows an example of the process of data simulation.

####  3.5.2 Training details

TS-VAD has achieved an excellent performance on CHIME6 [4] and DIHARD III [21]
challenge. Unlike the previous method using i-vector, we employ ResNet-based
x-vector as the target-speaker embedding.

The TS-VAD model is also similar to the VAD model except that the feature maps
produced by ResNet need to be concatenated with a target speaker embedding.
The concatenated features are then fed to the BiLSTM layers and fully
connected layers. Since TS-VAD training needs large computing resources, we
only train a model with $S=2$ for statistical pooling. The number of target
speakers embedding $N$ is set to 8. The parameters of front ResNet34 are
initialized from our speaker embedding model.

The model is first trained on the simulated LibriSpeech for 10 epochs with
front ResNet34 frozen, and then it is trained for another 10 epochs with all
parameters. Next, we transfer this model to VoxConverse data by training on
the simulated DEV402 for 10 epochs. Finally, we fine-tune the model on DEV402
for 200 epochs and validate on VAL46. The learning rate is 0.0001 when
training on simulated data and 0.00001 during the fine-tuning stage. The model
is optimized by BCE loss and Adam optimizer. The input is 16s chunked wav, and
the acoustic feature is 80-dim log Mel-filterbank energies with a frame length
of 25ms and a frameshift of 10ms.

####  3.5.3 Inference

For inference, the non-speech regions are first removed by VAD, and the wavs
are split into 16s chunks. Next, speaker embeddings are extracted given the
results from a clustering-based method. We only consider those speaker
embeddings with 16s or longer speech. If the number of speakers is less than
8, we use zero-vectors as the fake embeddings. If it is greater than 8, we
discard the speaker embeddings with shorter speech, but their labels are kept
in the final results. The threshold for speaker decision is set to 0.8.

###  3.6 2-speaker TS-VAD for Overlap Detection

The training data, model configuration, and training process are the same as
the TS-VAD in Sec. 3.5 except that the number of target speaker $N$ is 2. For
each recording, we select at most 5 speakers with the longest speech for
inference, which provides at most ${5\choose 2}=10$ pairs of target speaker
embeddings. After obtaining the speaker decision of each pair of speakers by a
threshold of 0.8, we update the results with the detected overlapped speech
regions. This 2-speaker TS-VAD method can be applied to any data without
considering the number of speakers. In this challenge, we only consider
modifying the overlapped speech regions, but the single speaker region can
also be iteratively refined according to the output from this 2-speaker TS-VAD
model. And we will do more experiments later.

###  3.7 Data Augmentation

We perform online data augmentation [22] with MUSAN and RIRs corpus. For
background additive noise, we use ambient noise, music, television, and babble
noise. For reverberation, we perform convolution with 40,000 simulated room
impulse responses from small and medium rooms. The data augmentation is
employed for all models which take acoustic features as input.

###  3.8 System Fusion

To further improve the performance and robustness, we fuse our systems by
DOVER-Lap [23].

##  4 Experimental Results

Table 3 shows the results on both VAL46 and the challenge test set. For the
clustering-based system, the AHC method achieves a DER of $4.42\%$ on VAL46
and $6.77\%$ on the test set, which is our first submission. We employed our
best VAD model for all systems on VAL46, but our first submission includes a
poor VAD model, and it may not correctly reveal the improvements brought by
OSD.

For the 2nd submission, we fused systems 3, 4, and 7 using DOVER-Lap with
custom weight tuned on VAL46, and we obtained a DER of $4.41\%$ on VAL46 and
$5.36\%$ on the test set.

For the 3rd submission, we fuse systems 3, 4, 5, and 7 using DOVER-Lap with
custom weight tuned on VAL46. However, the DER on the test set goes lower. We
did not know if it was system 5 that reduced the performance, so we directly
submitted system 5, and it shows a DER of $5.45\%$. Therefore, the reason may
be that we tuned the weights so aggressively, and the fused system is over-
fitted on VAL46.

Finally, we fused systems 3, 4, 5, 6, and 8 with rank-based weighting and
achieve a DER of $5.07\%$ on the test set, which ranked 1st at the VoxSRC
2021.

##  5 Conclusions

In this report, we describe our system for the VoxSRC 2021. To achieve better
performance, we mainly focus on the overlapped speech detection. We employ
overlap detection and TS-VAD to reduce the missed speaker error. In addition,
we propose a 2-speaker TS-VAD framework to detect the overlapped speech
between each pair of speakers. Our experiments show that detecting the
overlapped speech regions can significantly improve performance.

## References

  * [1] X. Xiao, N. Kanda, Z. Chen, T. Zhou, T. Yoshioka, S. Chen, Y. Zhao, G. Liu, Y. Wu, J. Wu _et al._ , “Microsoft speaker diarization system for the voxceleb speaker recognition challenge 2020,” in _ICASSP_. IEEE, 2021, pp. 5824–5828.
  * [2] K. Boakye, B. Trueba-Hornero, O. Vinyals, and G. Friedland, “Overlapped speech detection for improved speaker diarization in multiparty meetings,” in _ICASSP_. IEEE, 2008, pp. 4353–4356.
  * [3] Y. Fujita, S. Watanabe, S. Horiguchi, Y. Xue, and K. Nagamatsu, “End-to-end neural diarization: Reformulating speaker diarization as simple multi-label classification,” _arXiv preprint arXiv:2003.02966_ , 2020.
  * [4] I. Medennikov, M. Korenevsky, T. Prisyach, Y. Khokhlov, M. Korenevskaya, I. Sorokin, T. Timofeeva, A. Mitrofanov, A. Andrusenko, I. Podluzhny, A. Laptev, and A. Romanenko, “Target-Speaker Voice Activity Detection: A Novel Approach for Multi-Speaker Diarization in a Dinner Party Scenario,” in _INTERSPEECH_ , 2020, pp. 274–278.
  * [5] J. Carletta, S. Ashby, S. Bourban, M. Flynn, M. Guillemot, T. Hain, J. Kadlec, V. Karaiskos, W. Kraaij, M. Kronenthal _et al._ , “The ami meeting corpus: A pre-announcement,” in _International workshop on machine learning for multimodal interaction_. Springer, 2005, pp. 28–39.
  * [6] A. Janin, D. Baron, J. Edwards, D. Ellis, D. Gelbart, N. Morgan, B. Peskin, T. Pfau, E. Shriberg, A. Stolcke _et al._ , “The icsi meeting corpus,” in _ICASSP_ , vol. 1. IEEE, 2003, pp. I–I.
  * [7] Y. Fu, L. Cheng, S. Lv, Y. Jv, Y. Kong, Z. Chen, Y. Hu, L. Xie, J. Wu, H. Bu _et al._ , “Aishell-4: An open source dataset for speech enhancement, separation, recognition and speaker diarization in conference scenario,” _arXiv preprint arXiv:2104.03603_ , 2021.
  * [8] N. Ryant, K. Church, C. Cieri, A. Cristia, J. Du, S. Ganapathy, and M. Liberman, “The second dihard diarization challenge: Dataset, task, and baselines,” _arXiv preprint arXiv:1906.07839_ , 2019.
  * [9] N. Ryant, P. Singh, V. Krishnamohan, R. Varma, K. Church, C. Cieri, J. Du, S. Ganapathy, and M. Liberman, “The third dihard diarization challenge,” _arXiv preprint arXiv:2012.01477_ , 2020.
  * [10] A. Nagrani, J. S. Chung, W. Xie, and A. Zisserman, “Voxceleb: Large-scale speaker verification in the wild,” _Computer Speech & Language_, vol. 60, p. 101027, 2020.
  * [11] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: an asr corpus based on public domain audio books,” in _ICASSP_. IEEE, 2015, pp. 5206–5210.
  * [12] J. S. Chung, J. Huh, A. Nagrani, T. Afouras, and A. Zisserman, “Spot the Conversation: Speaker Diarisation in the Wild,” in _INTERSPEECH_ , 2020, pp. 299–303.
  * [13] D. Snyder, G. Chen, and D. Povey, “MUSAN: A Music, Speech, and Noise Corpus,” _arXiv:1510.08484_ , 2015.
  * [14] T. Ko, V. Peddinti, D. Povey, M. L. Seltzer, and S. Khudanpur, “A study on data augmentation of reverberant speech for robust speech recognition,” in _ICASSP_. IEEE, 2017, pp. 5220–5224.
  * [15] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _CVPR_ , 2016, pp. 770–778.
  * [16] M. Schuster and K. K. Paliwal, “Bidirectional recurrent neural networks,” _IEEE transactions on Signal Processing_ , vol. 45, no. 11, pp. 2673–2681, 1997.
  * [17] J. Deng, J. Guo, N. Xue, and S. Zafeiriou, “Arcface: Additive angular margin loss for deep face recognition,” in _CVPR_ , 2019.
  * [18] X. Qin, M. Li, H. Bu, R. K. Das, W. Rao, S. Narayanan, and H. Li, “The ffsvc 2020 evaluation plan,” _arXiv preprint arXiv:2002.00387_ , 2020.
  * [19] Q. Lin, R. Yin, M. Li, H. Bredin, and C. Barras, “LSTM Based Similarity Measurement with Spectral Clustering for Speaker Diarization,” in _INTERSPEECH_ , 2019, pp. 366–370.
  * [20] Q. Lin, W. Cai, L. Yang, J. Wang, J. Zhang, and M. Li, “DIHARD II is Still Hard: Experimental Results and Discussions from the DKU-LENOVO Team,” in _Odyssey_ , 2020, pp. 102–109.
  * [21] Y. Wang, M. He, S. Niu, L. Sun, T. Gao, X. Fang, J. Pan, J. Du, and C.-H. Lee, “Ustc-nelslip system description for dihard-iii challenge,” _arXiv preprint arXiv:2103.10661_ , 2021.
  * [22] W. Cai, J. Chen, J. Zhang, and M. Li, “On-the-fly data loader and utterance-level aggregation for speaker and language recognition,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 1038–1051, 2020.
  * [23] D. Raj, L. P. Garcia-Perera, Z. Huang, S. Watanabe, D. Povey, A. Stolcke, and S. Khudanpur, “Dover-lap: A method for combining overlap-aware diarization outputs,” in _SLT_. IEEE, 2021, pp. 881–888.