# Flux noise in disordered spin systems

José Alberto Nava Aquino  Rogério de Sousa  Department of Physics and
Astronomy, University of Victoria, Victoria, British Columbia V8W 2Y2, Canada
Centre for Advanced Materials and Related Technology, University of Victoria,
Victoria, British Columbia V8W 2Y2, Canada

###### Abstract

Impurity spins randomly distributed at the surfaces and interfaces of
superconducting wires are known to cause flux noise in Superconducting Quantum
Interference Devices (SQUIDs), providing a dominant mechanism for decoherence
in flux-tunable superconducting qubits. While flux noise is well characterized
experimentally, the microscopic model underlying spin dynamics remains a great
puzzle. The main problem is that first-principles theories based on
integration of the quantum Heisenberg equations of motion for interacting
spins are too computationally expensive to capture spin diffusion over large
length scales, hindering comparisons between microscopic models and
experimental data. In contrast, third principles approaches lump spin dynamics
into a single phenomenological spin-diffusion operator $D\nabla^{2}$, that is
not able to describe the quantum noise regime and connect to microscopic
models and different disorder scenarios such as spin clusters. Here we propose
an intermediate “second principles” method to describe general spin
dissipation and flux noise in the quantum regime. It leads to the
interpretation that flux noise arises from the density of paramagnon
excitations at the edge of the superconducting wire, with paramagnon-
paramagnon interactions leading to spin diffusion, and interactions between
paramagnons and other degrees of freedom such as phonons, electrons, and two-
level systems leading to spin energy relaxation. At high frequency $\omega$ we
obtain an upper bound for flux noise, showing that the (super)Ohmic noise
observed in experiments is not originating from interacting spin impurities.
We apply the method to Heisenberg models in two dimensional square lattices
with a random distribution of vacancies, with nearest-neighbor spins coupled
by constant exchange. Explicit numerical calculations of flux noise show that
it follows the observed power law $A/\omega^{\alpha}$, with amplitude $A$ and
exponent $\alpha$ depending on temperature and inhomogeneities such as spatial
confinement and disorder. These results are compared to experiments in niobium
and aluminum devices. The method establishes a connection between flux noise
experiments and microscopic Hamiltonians with the goal of identifying relevant
microscopic mechanisms and guiding strategies for reducing flux noise.

##  I Introduction

While progress in experimental realization of quantum computers based on
superconducting wires and Josephson junctions has been remarkable [1], the
noise level in current devices greatly reduces their capacity to solve
problems, and washes out their “quantum advantage”. One key issue is the trade
off between scalability and flux noise. Qubit frequency tunability is
essential to circumvent the frequency crowding problem faced by
superconducting circuits with more than 100 qubits. This requires the addition
of Superconducting Quantum Interference Devices (SQUIDs) to the circuit,
increasing their sensitivity to flux noise [2, 3]. A similar issue plagues
SQUID qubits, in that additional qubit interconnection increases the impact of
flux noise [4].

The origin of flux noise in superconducting devices remains unknown, although
there is consensus that it arises from the dynamics of spin centers (magnetic
impurities) near the superconducting wires [5, 6, 7, 8, 9, 10, 11] (see Fig.
1). This conclusion is supported by experiments showing noise amplitude
following the Curie susceptibility law ($\tilde{\chi}(\omega=0)\propto 1/T$,
where $T$ is temperature) [7]. However, there are conflicting opinions about
the magnetic order of the spins causing noise. While some authors claimed
proximity to a spin glass phase [7, 12], recent experiments were able to
measure a Curie-Weiss susceptibility $\tilde{\chi}(\omega=0)\propto
1/(T-T_{{\rm CW}})$ which rules out the spin glass scenario. Instead, [11]
measured $T_{{\rm CW}}\approx-10~{}{\rm mK}<0$ indicating proximity to an
_antiferromagnetic_ phase, while [13] measured $T_{CW}\approx+5~{}{\rm mK}>0$
indicating proximity to a _ferromagnetic_ phase.

Figure 1: Origin of flux noise. Section of a superconducting wire with spin
impurities randomly distributed at its surface. The flux produced by each spin
is given by $\hat{\Phi}_{i}=-\bm{F}(\bm{R}_{i})\cdot\bm{\hat{s}}_{i}$, where
$\bm{\hat{s}}_{i}$ is the spin operator of an impurity located at
$\bm{r}=\bm{R}_{i}$. The “flux vector” $\bm{F}(\bm{r})$ points along the
magnetic field produced by the current density $\bm{J}_{SC}$ shown in the
figure. Spin impurity dynamics leads to background flux noise that limits
coherence times for SQUID-based and flux-tunable qubits.

Measurements of flux noise [14, 15, 12, 11] are quite puzzling. They reveal
approximate $\tilde{S}_{\Phi}(\omega)\propto 1/\omega^{\alpha}$ frequency
dependence over several decades of frequency, and show quantum-noise asymmetry
$\tilde{S}_{\Phi}^{-}(\omega)=\tilde{S}_{\Phi}(\omega)-\tilde{S}_{\Phi}(-\omega)$
due to spontaneous emission [11]. To our knowledge all theories of flux noise
available in the literature are semiclassical so they are not able to describe
this asymmetry [8, 12, 16]. Moreover, flux noise was shown to become either
Ohmic $\propto\omega$ [17] or super-Ohmic $\propto\omega^{3}$ above $4$ GHz
[18, 11], and a key open question is whether or not this high frequency
contribution is also due to impurity spins.

The main mechanism for spin energy relaxation in ferromagnetic metals, the so
called Gilbert damping, relies on magnetic excitations (magnons) decaying into
electron-hole pair excitations in the metal. In a superconductor these are
exponentially suppressed at temperatures much lower than the superconducting
energy gap, making the Gilbert damping constant exponentially small [19]. In
addition, the weak magnetic fields in superconducting devices implies the
spin-orbit and hyperfine spin-flip rates are close to zero [5]. The only
remaining mechanism for spin energy decay in superconductors is due to the
interaction between each impurity spin with nearby amorphous two-level systems
[5]. Such an interaction leads to wide distributions of single-spin-flip rates
$\Gamma_{i}$ for different impurity spins $i$ [20].

It is known in nuclear magnetic resonance experiments that in the presence of
spin-spin interaction, spins with $\Gamma_{i}=0$ can relax by diffusing their
nonequilibrium magnetization towards sites with $\Gamma_{i}>0$ [21]. However,
no theory to date has been able to capture the coexistence of spin-spin
interaction with wide distributions of $\Gamma_{i}$.

When $\Gamma_{i}=0$ for all spins $i$, the total spin magnetization is
conserved, and the spin fluctuations due to spin-spin exchange interaction
necessarily obey a spin diffusion equation at long wavelengths [22, 23]. A
recent “pump and probe” experiment [13] measured the flux time correlation
function $\langle\hat{\Phi}(t)\hat{\Phi}(0)\rangle$ in SQUIDs, and showed that
it behaved similar to Brownian motion:
$\langle\Phi(t)\Phi(0)\rangle/\langle\Phi^{2}\rangle\propto 1-{\rm
const.}\sqrt{t}$ in the $1-1000$ $\mu$s time range. It provides evidence that
the dynamics of flux $\Phi$ in a superconducting device is described by the
phenomena of diffusion in this time range.

The usual theory for spin dynamics in a disordered spin system is based on the
assumption that the system is in a spin glass phase, that is uniform and
translational invariant [24]. These theories are not satisfactory for modeling
flux noise in superconducting devices for two reasons. First, there is
evidence that the impurity spins are in the paramagnetic (non-spin glass)
phase [11, 13]. Second, inhomogeneity and lack of translation invariance play
a crucial role. For example, the flux produced by spins located close to the
superconducting wire edge is much larger than the flux produced by spins away
from the edges [6, 25]. There is also the desire to know what is the impact of
nonhomogeneous spin distributions, such as impurity spin clustering [15, 16].

Describing spin diffusion from “first principles”, i.e. by integrating the
Heisenberg equations of motion for a model of interacting spins, is a well
known challenge of theoretical physics [23]. The standard method is what we
call “third-principles approach”: It assumes the spins can be described by a
continuous magnetization density $\bm{M}(\bm{r},t)$ that satisfies the
phenomenological equation $\frac{\partial\bm{M}}{\partial
t}=D\nabla^{2}\bm{M}$, with $D$ the spin diffusion constant. With all physical
properties lumped into a phenomenological constant $D$, the third-principles
theory can not establish a connection to microscopic model spin Hamiltonians,
spin spatial distributions, and the impact of inhomogeneity and clusters.

A serious shortcoming of the third principles approach is that it requires the
assumption of a hard boundary condition such as $\bm{M}(\bm{r},t)=\bm{0}$ at
wire edges and the surface of spin clusters in order to ensure total spin
conservation across the boundaries [8, 13]. A hard boundary condition like
this is unjustifiable, and in fact is known to be violated in magnetic systems
due to the appearence of confined surface or edge magnons [26]. Developing a
theory of spin dynamics that properly accounts for the boundary effects is of
crucial importance because flux noise is known to be dominated by spins at the
edge of the wire where the supercurrent is maximum [6, 25].

The purpose of this article is to propose a “second principles” theoretical
framework for spin dissipation (diffusion plus relaxation) that includes
quantum noise and is more “microscopic” than the usual third-principles
approach. The goal is to establish a connection between flux noise measured in
experiments and microscopic spin Hamiltonians, without the prohibitive
computational cost associated to the first-principles approach. To do this we
assume spin dissipation according to a random walk model governed by the
parameters of the spin Hamiltonian such as the microscopic exchange
interaction between each pair of spins.

Below we describe general theoretical results, and then present explicit
numerical calculations of flux noise for the Heisenberg model with nearest
neighbour interactions in the paramagnetic phase ($T>T_{c}^{{\rm mag}}$). Our
calculations are done in a finite spin lattice with a random distribution of
vacancies, showing explicit predictions for spatial confinement (wire edges)
and disorder due to random distribution of vacancies across the wire’s surface
as well as wide distributions of individual spin-flip rates $\Gamma_{i}$.

##  II Model for flux noise and linear response theory

We start by describing the impact of wire currents on impurity spins and how
it leads to a general expression for flux noise. The magnetic moment of an
impurity spin is given by $-g\mu_{B}\hat{\bm{s}}_{i}$, where $g\approx 2$ is
the $g$-factor, $\mu_{B}$ is the Bohr magneton, and $\hat{\bm{s}}_{i}$ is a
dimensionless spin operator for an impurity located at position $\bm{R}_{i}$.
It couples to the superconducting wire current density by producing a flux
[25],

$\hat{\Phi}=-\sum_{i}x_{i}\bm{F}(\bm{R}_{i})\cdot\hat{\bm{s}}_{i},$ (1)

where the sum goes over all sites $\bm{R}_{i}$ of a virtual square lattice
containing $N$ sites. The variable $x_{i}=1$ when there is a spin at the
virtual site, and $x_{i}=0$ otherwise, and the spin density is
$\sigma=\sum_{i}x_{i}/N=N_{s}/N$, where $N_{s}$ is the number of spins. The
flux vector $\bm{F}(\bm{r})$ is directly proportional to $\bm{B}_{I}(\bm{r})$,
the magnetic field produced by the wire’s current density:
$\bm{F}(\bm{r})=g\mu_{B}\bm{B}_{I}(\bm{r})/I$, where $I$ is the total current
flowing through the wire.

The wire’s current in turn affect spins by imprinting an _external_ local
field

$\bm{h}_{i}=-g\mu_{B}\bm{B}_{I}(\bm{R}_{i})=-I\bm{F}(\bm{R}_{i}),$ (2)

that couples to the spin according to ${\cal
H}_{c}=-\sum_{i}x_{i}\bm{h}_{i}\cdot\hat{\bm{s}}_{i}$, so $\bm{h}_{i}$ has
dimensions of energy. When the local field $\bm{h}_{i}$ is time dependent, the
spins respond according to

$\langle\hat{s}_{i}^{a}(t)\rangle_{h\neq
0}=\langle\hat{s}_{i}^{a}(t)\rangle_{h=0}+\sum_{j,b}\int_{-\infty}^{\infty}dt^{\prime}\chi_{ij}^{ab}(t-t^{\prime})h_{j}^{b}(t^{\prime}),$
(3)

where $a,b=x,y,z$ and the dynamical susceptibility is given by the linear
response formula,

$\chi^{ab}_{ij}(t-t^{\prime})=x_{i}x_{j}\frac{i}{\hbar}\theta(t-t^{\prime})\langle[\hat{s}_{i}^{a}(t),\hat{s}_{j}^{b}(t^{\prime})]\rangle.$
(4)

Defining spin noise as

$\tilde{S}^{ab}_{ij}(\omega)=x_{i}x_{j}\int_{-\infty}^{\infty}dte^{i\omega
t}\langle[\hat{s}_{i}^{a}(t)-\langle\hat{s}_{i}^{a}\rangle][\hat{s}_{j}^{b}(0)-\langle\hat{s}_{j}^{b}\rangle]\rangle,$
(5)

and using Eq. (4) we obtain the general relationship between susceptibility
and spin noise,

$\tilde{\chi}^{ab}_{ij}(\omega)=\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}d\omega^{\prime}\frac{1-e^{-\frac{\hbar\omega^{\prime}}{k_{B}T}}}{\omega^{\prime}-\omega-i\eta}\tilde{S}^{ab}_{ij}(\omega^{\prime}),$
(6)

where $\eta\rightarrow 0^{+}$. Taking the imaginary part and using the fact
that $\tilde{S}^{ab}_{ij}(\omega)=\tilde{S}^{ba}_{ji}(\omega)^{*}$, we get the
fluctuation-dissipation theorem for spins,

$\tilde{S}^{ab}_{ij}(\omega)=\frac{2\hbar}{1-e^{-\hbar\omega/k_{B}T}}\frac{1}{2i}[\tilde{\chi}^{ab}_{ij}(\omega)-\tilde{\chi}^{ba}_{ji}(\omega)^{*}].$
(7)

The flux noise is then given by

$\displaystyle\tilde{S}_{\Phi}(\omega)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dt\textrm{e}^{i\omega
t}\left\langle\delta\hat{\Phi}(t)\delta\hat{\Phi}(0)\right\rangle$ (8)
$\displaystyle=$
$\displaystyle\sum_{i,j,a,b}F^{a}(\bm{R}_{i})\tilde{S}^{ab}_{ij}(\omega)F^{b}(\bm{R}_{j}),$

where $\delta\hat{\Phi}(t)=\hat{\Phi}(t)-\langle\hat{\Phi}\rangle$ denotes
flux fluctuation.

##  III Theory of spin dynamics

Our goal is to compute flux noise for the general Heisenberg quantum spin
Hamiltonian,

${\cal
H}=-\frac{1}{2}\sum_{i,j}x_{i}x_{j}J_{ij}\hat{\bm{s}}_{i}\cdot\hat{\bm{s}}_{j}-\sum_{i}x_{i}\bm{h}_{i}\cdot\hat{\bm{s}}_{i},$
(9)

plus spin energy relaxation due to other degrees of freedom such as phonons,
electron-hole excitations, and two-level system defects [20]. Here $J_{ij}$ is
the exchange interaction between spins $i$ and $j$, that can be ferromagnetic
($J_{ij}>0$) or antiferromagnetic ($J_{ij}<0$), and $\bm{h}_{i}$ is the local
field defined in Eq. (2).

###  III.1 Static mean-field theory

In mean-field theory we neglect higher order fluctuations by approximating
$\langle\hat{s}_{i}^{a}\hat{s}_{j}^{b}\rangle\approx\langle\hat{s}_{i}^{a}\rangle\langle\hat{s}_{j}^{b}\rangle$.
We simplify the notation by writing
$\langle\hat{\bm{s}}_{i}\rangle=\bm{s}_{i}$, i.e. the spin vector without a
hat denotes the average of the spin operator (a real vector). The mean-field
approximation is exactly the same as the “classical spin model” used by many
authors, e.g. [16]. An additional approximation in mean-field theory is to
assume the system’s entropy can be written as a sum of single-spin entropies
[27]:

$\langle{\cal
S}\rangle=k_{B}\sum_{i}x_{i}\left[\ln{2}-2|\bm{s}_{i}|^{2}-\frac{4}{3}(|\bm{s}_{i}|^{2})^{2}+\mathcal{O}((|\bm{s}_{i}|^{2})^{3})\right].$
(10)

This expression is specific to spin-1/2 impurities. Note that cutting the
expansion to fourth order affects the result only when $T\ll T_{c}^{{\rm
mag}}$, where $T_{c}^{{\rm mag}}$ is a critical temperature for a phase where
$s_{i}>0$. The free energy is thus given by:

$\displaystyle{\cal F}$ $\displaystyle=$ $\displaystyle\langle{\cal
H}\rangle-T\langle{\cal
S}\rangle=-\frac{1}{2}\sum_{i,j}x_{i}x_{j}J_{ij}\bm{s}_{i}\cdot\bm{s}_{j}-\sum_{i}x_{i}\bm{h}_{i}\cdot\bm{s}_{i}$
(11)
$\displaystyle-(k_{B}T)\sum_{i}x_{i}\left[\ln{2}-2|\bm{s}_{i}|^{2}-\frac{4}{3}(|\bm{s}_{i}|^{2})^{2}\right].$

Usually thermal equilibrium is realized by the set of $\bm{s}_{i}$ that leads
to the global minimum of the free energy. For example, take $x_{i}=1$ for all
$i$, $J_{ij}=J>0$ for nearest neighbors and zero otherwise, and a lattice with
periodic boundary conditions (b.c.). In this case the global minimum of ${\cal
F}$ is realized by $\bm{s}_{i}=\bm{s}^{{\rm eq}}$ for all $i$ (the
ferromagnetic homogeneous state), with free energy given by

$\frac{{\cal F}}{N}=\left(2k_{B}T-\frac{zJ}{2}\right)(s^{{\rm
eq}})^{2}+\frac{4}{3}k_{B}T(s^{{\rm eq}})^{4}-k_{B}T\ln{2},$ (12)

where $z$ is the number of nearest neighbours for each site of the lattice.
From this expression we see that a global minimum with $s^{{\rm eq}}>0$
appears only when the first term changes sign, leading to critical temperature
$k_{B}T_{c}^{{\rm mag}}=k_{B}T_{{\rm CW}}=zJ/4$. The same calculation can be
done for $J<0$ when the lattice can be partitioned into two sublattices with
one being n.n. to the other. For this case the global minimum of ${\cal F}$ is
realized by $\bm{s}_{i}=+\bm{s}^{{\rm eq}}$ for one sublattice and
$\bm{s}_{i}=-\bm{s}^{{\rm eq}}$ for the other (antiferromagnetic homogeneous
state). This leads to $k_{B}T_{c}=z|J|/4$. However, $k_{B}T_{{\rm
CW}}=-z|J|/4$ because the magnetic susceptibility does not have a singularity
at $T_{c}$ (it’s the staggered susceptibility that is singular at $T_{c}$).
These are the well-known mean-field results for phase transitions in the
spin-$1/2$ Heisenberg model [27].

###  III.2 Dynamical mean-field theory

Based on fundamental theories of spin dynamics [27] we propose the following
generalized equation of motion for the spins,

$\frac{d\bm{s}_{i}}{dt}=\frac{1}{\hbar}\bm{s}_{i}\times\bm{H}_{i}-\sum_{j}D_{ij}\bm{H}_{j}-\Gamma_{i}(\bm{s}_{i}-\bm{s}_{i}^{{\rm
inst~{}eq}}),$ (13)

expected to be valid for frequencies smaller than a cut-off $\Omega_{c}$ to be
discussed later. In addition to the usual spin precession, this includes a
discrete version of the _intra-spin dissipation_ operator $D_{ij}$, together
with isotropic _spin energy relaxation_ $\Gamma_{i}$. The $\Gamma_{i}$ drives
$\bm{s}_{i}$ towards its “instantaneous equilibrium” value

$\bm{s}_{i}^{{\rm inst~{}eq}}=\bm{s}_{i}^{{\rm
eq}}+\sum_{j}\tilde{\chi}_{ij}(0)\delta\bm{h}_{j}(t),$ (14)

which is time-dependent due to local field dynamics,
$\bm{h}_{i}(t)=\bm{h}_{i}^{{\rm eq}}+\delta\bm{h}_{i}(t)$, where
$\bm{h}_{i}^{{\rm eq}}$ is the static part. The quantities $\bm{s}_{i}^{{\rm
eq}}$ do not depend on time, they are thermal equilibrium averages calculated
assuming $\delta\bm{h}_{i}(t)=0$, i.e. they only depend on $\bm{h}_{i}^{{\rm
eq}}$ and other static free energy parameters;
$\tilde{\chi}_{ij}^{ab}(0)=\partial(s_{i}^{{\rm eq}})^{a}/\partial h_{j}^{b}$
is the $\omega=0$ susceptibility, assumed isotropic ($\propto\delta_{ab}$) to
be consistent with our Hamiltonian (9). We call Eq. (14) the “instantaneous
approximation”, because it assumes the other degrees of freedom causing spin
energy relaxation relax much faster than the spins themselves, so that the
spin system remains in thermal equilibrium with the other non-spin degrees of
freedom at all times. Note how Eq. (14) introduces the $\omega=0$
susceptibility _self-consistently_ into the equation of motion (13).

The _internal_ spin field is defined as

$\displaystyle\bm{H}_{i}$ $\displaystyle=$
$\displaystyle-\frac{\partial\mathcal{F}}{\partial\bm{s}_{i}}$
$\displaystyle=$ $\displaystyle
x_{i}\left\\{\sum_{j}x_{j}J_{ij}\bm{s}_{j}+\bm{h}_{i}-4k_{B}T\left[1+\frac{4}{3}s_{i}^{2}\right]\bm{s}_{i}\right\\},$

and the thermal equilibrium spin averages $\bm{s}_{i}^{{\rm eq}}$ are
determined by imposing time independence, $\delta\bm{h}_{i}(t)=\bm{0}$ and
$\frac{d\bm{s}_{i}}{dt}=\bm{0}$ for all $i$. This implies $\bm{s}_{i}^{{\rm
eq}}$ must be found by solving the system of equations

$\frac{1}{\hbar}\bm{s}_{i}^{{\rm eq}}\times\bm{H}_{i}^{{\rm
eq}}-\sum_{j}D_{ij}\bm{H}_{j}^{{\rm eq}}=\bm{0},$ (16)

where $\bm{H}_{i}^{{\rm eq}}$ is Eq. (III.2) with $\bm{s}_{i}=\bm{s}_{i}^{{\rm
eq}}$ and $\bm{h}_{i}=\bm{h}_{i}^{{\rm eq}}$. Note how Eq. (16) is always
satisfied for $\bm{H}_{i}^{{\rm eq}}=\bm{0}$, a smooth local minimum of the
free energy. However, other solutions with $\bm{H}_{i}^{{\rm eq}}\neq\bm{0}$
may arise in the presence of site-dependent local fields $\bm{h}_{i}^{{\rm
eq}}$.

The three terms in the right hand side of Eq. (13) correspond to reactive
dynamics, intra-spin-system dissipation (e.g. diffusion), and spin energy
relaxation due to other degrees of freedom, respectively. The reactive term is
non-dissipative, it does not change sign under time reversal $t\rightarrow-t$
so it has the same symmetry as the left hand side. The second and third terms
on the right hand side do change sign under time reversal, leading to an
irreversible approach to thermal equilibrium (the arrow of time). These terms
must be added ad hoc to the linearized equations of motion so that the zeroth
(attainment of thermal equilibrium) and second (entropy always increases) laws
of thermodynamics are obeyed. That is, the system is able to reach thermal
equilibrium, and the free energy always decreases as a function of time when
the system is in contact with a thermal reservoir.

A few notes about the microscopic origin of spin dissipation $D_{ij}$ are
warranted. The normal modes of Eq. (13) are called _magnons_ and _paramagnons_
, to be defined below. The reactive terms of Eq. (13) describe the dynamics of
noninteracting (para)magnons obtained by the mean-field approximation. Exactly
the same results are obtained from different methods, e.g. using a Holstein-
Primakoff transformation to convert spin operators into Bosonic
creation/destruction operators; transforming Hamiltonian (9) and keeping
contributions that are quadratic in these Bosonic operators leads to the same
magnon modes as Eq. (13) with $D_{ij}=\Gamma_{i}=0$ [26]. However, the higher
order terms that are neglected in this quadratic approximation can be
interpreted as describing (para)magnon-(para)magnon interactions. The
introduction of $D_{ij}\neq 0$ serves to account for these interactions
phenomenologically.

##  IV Specification of dissipation matrix $D_{ij}$ in the presence of
confinement and disorder

To go beyond the third-principles assumption of long wavelength spin diffusion
we need to come up with a specification for $D_{ij}$ that respects several
physical constraints. To do this, we take inspiration from random walk models
in a lattice. The key idea is that exchange interaction $J_{ij}$ is the main
driver for each random walk step, a spin “flip-flop”. A sequence of many flip-
flops will lead to diffusion. The constraint of total spin conservation
motivates our postulation of the following spin dissipation matrix:

$D_{ij}=\frac{d_{0}(T)}{\hbar\bar{J}_{c}}\left(x_{i}x_{j}|J_{ij}|-\delta_{ij}\sum_{k}x_{i}x_{k}|J_{ik}|\right).$
(17)

Here $d_{0}(T)$ is a function of temperature to be determined by fitting the
theory to experiments (note $d_{0}(T)$ is dimensionless). This is introduced
to account for critical behaviour of the spin diffusion constant near
$T_{c}^{{\rm mag}}$ [28]. The quantity

$\bar{J}_{c}=\frac{1}{N_{c}}\sum_{j,k\in{\rm cluster~{}}c}x_{j}x_{k}|J_{jk}|$
(18)

is the average exchange times coordination number for $c^{th}$ cluster, the
cluster that contains spin $i$. Such a cluster is defined as the set of all
spins $j$ such that either $J_{ij}\neq 0$ or there exists a set of sites
$k_{1},k_{2},\ldots,k_{n}$ such that
$J_{ik_{1}}J_{k_{1}k_{2}}J_{k_{2}k_{3}}\cdots J_{k_{n}j}\neq 0$. $N_{c}$ is
the number of spins in the isolated cluster $c$. The following motivates this
choice:

  1. 1.

When the local external field $\bm{h}_{i}$ and spin relaxation rate
$\Gamma_{i}$ are both zero, Eq. (17) preserves total spin. For each isolated
spin cluster, $\sum_{i,j\in{\rm cluster}}D_{ij}\bm{H}_{j}=\bm{0}$, therefore
summing Eq. (13) over all spins in a cluster leads to
$\frac{d}{dt}(\sum_{i\in{\rm cluster}}\bm{s}_{i})=\bm{0}$, so the total spin
in each cluster is a constant of the motion. Therefore, we do not need to
assume hard boundary conditions [8, 13] to describe confined systems such as
spin clusters and wire edges.

  2. 2.

This choice for $D_{ij}$ gives rise to diffusion in the long wavelength
regime. E.g. for the homogeneous nearest-neighbor model in the square lattice
with $J_{ij}=J_{i,i+\bm{v}}=J$ for $\bm{v}=\pm a\hat{x},\pm a\hat{y}$, we get

$-\sum_{i,j}\delta(\bm{r}-\bm{R}_{i})D_{ij}\bm{H}_{j}=\frac{d_{0}(T)a^{2}k_{B}T}{\hbar}\nabla^{2}\bm{M}(\bm{r}),$
(19)

where we assumed high temperature ($k_{B}T\gg J$ and $(\bm{s}_{i})^{2}\ll 1$)
and took the continuum limit by defining the magnetization density
$\bm{M}(\bm{r})=-\sum_{i}\bm{s}_{i}\delta(\bm{r}-\bm{R}_{i})$. Equation (19)
may be compared to experiments that show spin diffusion constant increasing
with temperature [12].

  3. 3.

Consider the time derivative of the free energy in each isolated spin cluster
when $\bm{h}_{i}$ is independent of time:

$\displaystyle\frac{d\mathcal{F}}{dt}$ $\displaystyle=$
$\displaystyle\sum_{i}\frac{\partial\mathcal{F}}{\partial\bm{s}_{i}}\cdot\frac{d\bm{s}_{i}}{dt}=\frac{d_{0}(T)}{\hbar\bar{J}_{c}}\sum_{i}\bm{H}_{i}\cdot\sum_{j}x_{i}x_{j}|J_{ij}|(\bm{H}_{j}-\bm{H}_{i})+\sum_{i}\Gamma_{i}\bm{H}_{i}\cdot(\bm{s}_{i}-\bm{s}_{i}^{{\rm
eq}})$ (20) $\displaystyle=$
$\displaystyle-\frac{d_{0}(T)}{\hbar\bar{J}_{c}}\sum_{i<j}x_{i}x_{j}|J_{ij}|(\bm{H}_{i}-\bm{H}_{j})^{2}-\sum_{i}\Gamma_{i}\left(\mathcal{F}(\bm{s}_{i})-\mathcal{F}(\bm{s}_{i}^{{\rm
eq}})+\mathcal{O}\left[(\delta\bm{s}_{i})^{2}\right]\right).$

The first term on the RHS is always negative, showing that our choice for
$D_{ij}$ always tends to decrease the free energy as time increases (i.e. it
obeys the 2nd law of thermodynamics). This justifies our use of modulus of
$J_{ij}$ in Eq. (17).

The second term on the RHS of Eq. (20) is negative provided that the deviation
out of equilibrium is small and $\mathcal{F}(\bm{s}_{i}^{{\rm eq}})$ is a
local minimum of the free energy. Therefore, the coupling to other nonspin
degrees of freedom $\Gamma_{i}$ pushes the system towards a local minimum of
the free energy, without subjecting to spin conservation.

##  V Calculation of dynamical susceptibility

For small deviations from equilibrium we write $\bm{h}_{i}=\bm{h}_{i}^{{\rm
eq}}+\delta\bm{h}_{i}(t)$, and $\bm{s}_{i}=\bm{s}_{i}^{{\rm
eq}}+\delta\bm{s}_{i}(t)$, where both $\delta\bm{h}_{i}(t)$ and
$\delta\bm{s}_{i}(t)$ are small time-dependent perturbations.

We plug these into the equation of motion (13), drop non-linear terms such as
$(\delta\bm{s}_{i})^{2},(\delta\bm{h}_{i}(t))^{2}$, and use Eq. (16) to
simplify:

$\displaystyle\frac{d}{dt}(\delta\bm{s}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar}\bm{s}_{i}^{{\rm
eq}}\times\left\\{(\delta\bm{h}_{i})+\sum_{j}x_{j}J_{ij}(\delta\bm{s}_{j})-4k_{B}T\left[1+\frac{4}{3}(s_{i}^{{\rm
eq}})^{2}\right](\delta\bm{s}_{i})-\frac{32}{3}k_{B}T[\bm{s}_{i}^{{\rm
eq}}\cdot(\delta\bm{s}_{i})]\bm{s}_{i}^{{\rm eq}}\right\\}$ (21)
$\displaystyle-\sum_{j}D_{ij}x_{j}\left\\{(\delta\bm{h}_{j})+\sum_{k}x_{k}J_{jk}(\delta\bm{s}_{k})-4k_{B}T\left[1+\frac{4}{3}(s_{j}^{{\rm
eq}})^{2}\right](\delta\bm{s}_{j})-\frac{32}{3}k_{B}T[\bm{s}_{j}^{{\rm
eq}}\cdot(\delta\bm{s}_{j})]\bm{s}_{j}^{{\rm eq}}\right\\}$
$\displaystyle+\frac{1}{\hbar}(\delta\bm{s}_{i})\times\bm{H}_{i}^{{\rm
eq}}-\Gamma_{i}(\delta\bm{s}_{i})+\Gamma_{i}\sum_{j}\tilde{\chi}_{ij}(0)(\delta\bm{h}_{j}).$

We assume $\bm{s}_{i}^{{\rm eq}}=s_{i}^{{\rm eq}}\hat{\bm{z}}$ and
$\bm{H}^{{\rm eq}}_{i}=H^{{\rm eq}}_{i}\bm{\hat{z}}$ and break this down into
two equations, one for $\delta s_{i}^{z}$ obtained by dot product with
$\bm{\hat{z}}$ on both sides of Eq. (21), and the other for $\delta
s_{i}^{+}=\delta s^{x}_{i}+i\delta s^{y}_{i}$ obtained by dot product with
$(\bm{\hat{x}}+i\bm{\hat{y}})$. Taking the time Fourier transform we get two
decoupled equations:

$\displaystyle(\omega\bm{{\sf I}}-\bm{{\sf P}})\cdot\delta\tilde{\bm{s}}^{z}$
$\displaystyle=$ $\displaystyle
i\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\delta\tilde{\bm{h}}^{z},$ (22a) $\displaystyle(\omega\bm{{\sf
I}}-\bm{{\sf M}})\cdot\delta\tilde{\bm{s}}^{+}$ $\displaystyle=$
$\displaystyle\left[i\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)-\frac{1}{\hbar}\bm{{\sf s}}^{{\rm
eq}}\right]\cdot\delta\tilde{\bm{h}}^{+},$ (22b)

where $\delta\tilde{\bm{s}}^{z},\delta\tilde{\bm{s}}^{+}$ and
$\delta\tilde{\bm{h}}^{z},\delta\tilde{\bm{h}}^{+}$ are $N$-component column
vectors, and $\bm{{\sf
I}},\bm{{\sf\Gamma}},\bm{{\sf\tilde{\chi}}}_{0},\bm{{\sf D}},\bm{{\sf
s}}^{{\rm eq}}$ are $N\times N$ matrices. They are defined by $[\bm{{\sf
I}}]_{ij}=x_{i}\delta_{ij}$,
$[\bm{{\sf\Gamma}}]_{ij}=x_{i}\Gamma_{i}\delta_{ij}$,
$[\bm{{\sf\tilde{\chi}}}_{0}]_{ij}=\tilde{\chi}_{ij}(0)$, $[\bm{{\sf
D}}]_{ij}=x_{i}x_{j}D_{ij}$, and $[\bm{{\sf s}}^{{\rm eq}}]_{ij}=x_{i}s^{{\rm
eq}}_{i}\delta_{ij}$.

The matrices $\bm{{\sf P}}$ and $\bm{{\sf M}}$ are the _paramagnon_ and
_magnon_ matrices, respectively. They are given by

$\displaystyle\bm{{\sf P}}$ $\displaystyle=$
$\displaystyle-i\left\\{\bm{{\sf\Gamma}}+\bm{{\sf D}}\cdot\bm{{\sf
J}}-4k_{B}T\bm{{\sf D}}\cdot\left[\bm{{\sf I}}+4\left(\bm{{\sf s}}^{{\rm
eq}}\right)^{2}\right]\right\\},$ (23a) $\displaystyle\bm{{\sf M}}$
$\displaystyle=$ $\displaystyle-i\left\\{\bm{{\sf\Gamma}}+\bm{{\sf
D}}\cdot\bm{{\sf J}}-4k_{B}T\bm{{\sf D}}\cdot\left[\bm{{\sf
I}}+\frac{4}{3}\left(\bm{{\sf s}}^{{\rm eq}}\right)^{2}\right]\right\\}$
$\displaystyle+$ $\displaystyle\frac{1}{\hbar}$ $\displaystyle\left\\{\bm{{\sf
H}}^{{\rm eq}}-\bm{{\sf s}}^{{\rm eq}}\cdot\left[\bm{{\sf
J}}-4k_{B}T\left(\bm{{\sf I}}+\frac{4}{3}\left(\bm{{\sf s}}^{{\rm
eq}}\right)^{2}\right)\right]\right\\},$ (23b)

where $[\bm{{\sf J}}]_{ij}=x_{i}x_{j}J_{ij}$, and $[\bm{{\sf H}}^{{\rm
eq}}]_{ij}=H^{{\rm eq}}_{i}\delta_{ij}$.

The eigenvalues of $\bm{{\sf P}}$ and $\bm{{\sf M}}$ are paramagnon and magnon
frequencies, respectively. While paramagnons cause spin fluctuations along
$\bm{s}_{i}^{{\rm eq}}$, the magnons cause fluctuations perpendicular to
$\bm{s}_{i}^{{\rm eq}}$. In the paramagnetic phase ($T>T_{c}^{{\rm mag}}$)
with zero external fields, $\bm{{\sf s}}^{{\rm eq}}=\bm{{\sf H}}^{{\rm
eq}}=\bm{{\sf 0}}$ and the matrices $\bm{{\sf P}}$ and $\bm{{\sf M}}$ become
identical, signaling the presence of isotropic spin fluctuations (i.e.
paramagnons are three-fold degenerate). From now on we shall focus our
discussion on this paramagnetic regime, so we only need to consider the
spectrum of $\bm{{\sf P}}$.

Diagonalize $\bm{{\sf P}}$ with a transformation $\bm{{\sf U}}$ such that

$\bm{{\sf U}}^{-1}\cdot\bm{{\sf P}}\cdot\bm{{\sf U}}=\bm{{\sf
P}}_{d}=-i\sum_{m}\gamma_{m}\bm{\hat{e}}_{m}\otimes\bm{\hat{e}}_{m}^{T},$ (24)

where $m$ labels the paramagnon mode with frequency $-i\gamma_{m}$, with
$\bm{\hat{e}}_{m}$ unit column vectors,
$\bm{\hat{e}}_{m}^{T}=(0,\ldots,0,1,0,\ldots,0)$, etc. All elements of
$\bm{{\sf P}}$ are pure complex, therefore the $\bm{{\sf U}}$ and $\bm{{\sf
U}}^{-1}$ can be chosen to have real elements. Take the complex conjugate of
Eq. (24) and use $\bm{{\sf P}}^{*}=-\bm{{\sf P}}$ to see that
$\gamma^{*}_{m}=\gamma_{m}$.

Apply $\bm{{\sf U}}^{-1}$ on both sides of Eq. (22a),

$\left(\omega-\bm{{\sf P}}_{d}\right)\cdot\bm{{\sf
U}}^{-1}\cdot\delta\tilde{\bm{s}}^{z}=i\bm{{\sf
U}}^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\delta\tilde{\bm{h}}^{z},$ (25)

and invert the diagonal matrix to get an exact expression for the
susceptibility:

$(\delta\tilde{\bm{s}}^{z})=i\sum_{m}\bm{{\sf
U}}\cdot\frac{\bm{\hat{e}}_{m}\otimes\bm{\hat{e}}_{m}^{T}}{\omega+i\gamma_{m}}\cdot\bm{{\sf
U}}^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\delta\tilde{\bm{h}}^{z}.$ (26)

Transformation (24) implies that column vector $\bm{v}_{m}=\bm{{\sf
U}}\cdot\bm{\hat{e}}_{m}$ is the right eigenvector of $\bm{{\sf P}}$
associated to the eigenvalue $-i\gamma_{m}$. Similarly, the line vector
$\bm{v}_{m}^{-1}=\bm{\hat{e}}_{m}^{T}\cdot\bm{{\sf U}}^{-1}$ is the left
eigenvector of $\bm{{\sf P}}$ associated to the same eigenvalue. We have
$\bm{v}^{-1}_{m}\cdot\bm{v}_{m^{\prime}}=\delta_{mm^{\prime}}$, but the set
$\\{\bm{v}_{m}\\}$ is not mutually orthogonal because $[\bm{{\sf P}},\bm{{\sf
P}}^{{\dagger}}]\neq\bm{{\sf 0}}$ ($\bm{{\sf P}}$ is not normal).

Using these results, Eqs. (22a) and (26) imply the following exact expression
for the isotropic dynamical susceptibility:

$\displaystyle\bm{{\sf\tilde{\chi}}}^{ab}(\omega)$ $\displaystyle=$
$\displaystyle i\left(\omega\bm{{\sf I}}-\bm{{\sf
P}}\right)^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\delta_{ab}$ (27a) $\displaystyle=$ $\displaystyle
i\sum_{m}\frac{\bm{v}_{m}\otimes\left[\bm{v}_{m}^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\right]}{\omega+i\gamma_{m}}\delta_{ab},$ (27b)

valid for $T>T_{c}^{{\rm mag}}$ with $a,b=x,y,z$. The dynamical susceptibility
has poles at the paramagnon frequencies $\omega=-i\gamma_{m}$, and these
contribute to dissipation and noise. The paramagnons are said to be purely
dissipative because $\gamma_{m}$ is real, leading to $e^{-i\omega
t}=e^{-\gamma_{m}t}$ for the decay of the spin excitations. For $T>T_{c}^{{\rm
mag}}$ we have $\gamma_{m}\geq 0$ because the paramagnetic phase is stable.

Conservation of total spin follows from $\bm{1}^{T}\cdot\bm{{\sf
D}}=\bm{0}^{T}$, where $\bm{1}$ is the column vector with all $N$ components
equal to $1$. When in addition $\bm{{\sf\Gamma}}=\bm{{\sf 0}}$ all
contributions to the paramagnon matrix (23a) have $\bm{{\sf D}}$ on the left,
so it follows that the $m=0$ mode defined by
$\bm{\hat{v}}^{-1}_{0}=\bm{1}^{T}/\sqrt{N}$ is a left eigenvector of $\bm{{\sf
P}}$ associated to $\gamma_{0}=0$. This is true for general $J_{ij}$. Because
$\bm{v}_{0}^{-1}\cdot\bm{{\sf D}}=\bm{0}^{T}$, $m=0$ does not contribute to
the sum in the dynamical susceptibility Eq. (27b). This occurs as a direct
consequence of the conservation law
$\frac{d}{dt}(\bm{1}^{T}\cdot\delta\bm{s}^{z})=\bm{0}$, so we say $m=0$ is the
nondissipative Goldstone paramagnon.

At high frequency $\omega\gg{\rm Max}_{m}\\{\gamma_{m}\\}$, Eq. (27b) leads to

$\bm{{\sf\tilde{\chi}}}^{ab}(\omega)\approx\frac{i\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)}{\omega}\delta_{ab},$ (28)

because $\sum_{m}\bm{v}_{m}\otimes\bm{v}_{m}^{-1}=\bm{{\sf I}}$. As Eq. (13)
has an upper frequency cut-off, Eq. (28) should be taken as an upper bound on
the modulus of the susceptibility.

Now consider the opposite limit, $\omega\rightarrow 0$. Since
$\bm{{\sf\tilde{\chi}}}^{ab}(0)\equiv\bm{{\sf\tilde{\chi}}}_{0}\delta_{ab}$,
setting $\omega=0$ in Eq. (27a) leads to $\left(\bm{{\sf
P}}+i\bm{{\sf\Gamma}}\right)\cdot\bm{{\sf\tilde{\chi}}}_{0}=i\bm{{\sf D}}$,
and using Eq. (23a) we get

$\bm{{\sf\tilde{\chi}}}^{ab}(\omega=0)=\bm{{\sf\tilde{\chi}}}_{0}\delta_{ab}=\left(4k_{B}T\bm{{\sf
I}}-\bm{{\sf J}}\right)^{-1}\delta_{ab}.$ (29)

This is the generalized Curie-Weiss susceptibility for a nonhomogeneous spin
system (valid for $T>T_{c}^{{\rm mag}}$).

Exact analytical results can be obtained for the special case of a
translation-invariant spin system. If the system is close to being
translation-invariant, e.g. only a few vacancies are present so that
$\sigma=\langle x_{i}\rangle_{i}\lesssim 1$, and $\Gamma_{i}$ does not depend
appreciably on $i$, a homogeneous approximation (HA) can be proposed. The HA
replaces $x_{i}$, $\Gamma_{i}$, $J_{ij}$, and $D_{ij}$ by their spatial
averages $\sigma$, $\bar{\Gamma}$, $\bar{J}_{ij}$, and $\bar{D}_{ij}$ (the
latter two depending only on $\bm{R}_{j}-\bm{R}_{i}$), making the problem
analytically solvable. In this case Appendix A shows that the exact
susceptibility can be obtained by spatial Fourier transformation,

$\tilde{\chi}^{ab}(\bm{q},\omega)=\frac{i\left(\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right)\delta_{ab}}{\omega+i\left\\{-\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]+\bar{\Gamma}\right\\}},$
(30)

where
$\tilde{D}(\bm{q})=\sum_{\bm{v}}\bar{D}_{i,i+\bm{v}}e^{-i\bm{q}\cdot\bm{v}}$
and
$\tilde{J}(\bm{q})=\sum_{\bm{v}}\bar{J}_{i,i+\bm{v}}e^{-i\bm{q}\cdot\bm{v}}$.
The paramagnon modes are labelled by $m=\bm{q}\in 1^{{\rm st}}$ Brillouin
zone, each with frequency eigenvalue

$\gamma_{m}=\gamma_{\bm{q}}=-\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]+\bar{\Gamma},$
(31)

and associated right and left eigenvectors $\bm{v}_{m}=\bm{e}_{\bm{q}}$,
$\bm{v}_{m}^{-1}=\bm{e}_{\bm{q}}^{{\dagger}}/N_{s}$, respectively, where
$\bm{e}^{{\dagger}}_{\bm{q}}=\left(e^{-i\bm{q}\cdot\bm{R}_{0}},\ldots,e^{-i\bm{q}\cdot\bm{R}_{N-1}}\right)$,
and $N_{s}$ is the number of occupied sites forming a translation-invariant
lattice.

When $\omega\rightarrow 0$, Eq. (30) leads to

$\tilde{\chi}^{ab}(\bm{q},0)=\frac{1}{4k_{B}}\frac{1}{T-T_{{\rm
CW}}(\bm{q})}\delta_{ab},$ (32)

where $T_{{\rm CW}}(\bm{q})=\tilde{J}(\bm{q})/(4k_{B})$ is the Curie-Weiss
temperature in Fourier space.

Without translation invariance, e.g. in the presence of spin clusters, Eq.
(29) shows that the $\omega=0$ susceptibility may have several different
temperature poles $T_{{\rm CW}}(0)$, each associated with different clusters
having different $J_{ij}$ or number of neighbours. However, $T_{{\rm CW}}$
only depends on spin-spin interaction, it does not depend on relaxation
parameters $\Gamma_{i}$.

##  VI Flux noise and paramagnon density

The explicit expression for flux noise is obtained by plugging Eq. (27b) into
Eq. (7) and using Eq. (8)

$\tilde{S}_{\Phi}(\omega)=\frac{2\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\sum_{m,a}\frac{\bm{F}^{aT}\cdot\bm{v}_{m}\left[\bm{v}_{m}^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\bm{F}^{a}\right]}{\omega^{2}+\gamma_{m}^{2}},$ (33)

where $\bm{F}^{aT}=\left(F^{a}(\bm{R}_{0}),\ldots,F^{a}(\bm{R}_{N-1})\right)$
represents the $a$-component of the flux vector for all spins. A convenient
way to interpret this expression is to write it in terms of a density of
Lorentzian contributions

$\tilde{S}_{\Phi}(\omega)=\frac{2\pi\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\int_{-\infty}^{\infty}d\gamma\frac{\gamma/\pi}{\omega^{2}+\gamma^{2}}\rho_{\Phi}(\gamma),$
(34)

where $\rho_{\Phi}(\gamma)$ is the paramagnon flux density, defined as

$\rho_{\Phi}(\gamma)=\frac{1}{\gamma}\sum_{m,a}(\bm{F}^{aT}\cdot\bm{v}_{m})\left[\bm{v}_{m}^{-1}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\bm{F}^{a}\right]\delta(\gamma-\gamma_{m}),$ (35)

where $\delta(x)$ is the Dirac delta function.

When the spin system is translation-invariant, the paramagnon flux density is
given by

$\rho_{\Phi}(\gamma)=\frac{1}{N_{s}\gamma}\sum_{\bm{q}}\left|\tilde{\bm{F}}(\bm{q})\right|^{2}\left[\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right]\delta(\gamma-\gamma_{\bm{q}}),$
(36)

where $\tilde{\bm{F}}(\bm{q})=\sum_{j}\bm{F}_{j}e^{-i\bm{q}\cdot\bm{R}_{j}}$
is the flux vector in Fourier space.

Finally, in the high frequency limit $\omega\gg{\rm Max}_{m}\\{\gamma_{m}\\}$
Eq. (33) implies

$\tilde{S}_{\Phi}(\omega)=\frac{2\hbar}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\frac{1}{\omega}\sum_{a}\bm{F}^{aT}\cdot\left(\bm{{\sf\Gamma}}\cdot\bm{{\sf\tilde{\chi}}}_{0}-\bm{{\sf
D}}\right)\cdot\bm{F}^{a}.$ (37)

Since Eq. (13) has an upper frequency cut-off $\Omega_{c}$, the actual
$\tilde{S}_{\Phi}(\omega)$ is expected to be less than Eq. (37) at
$\omega>\Omega_{c}$. In this regime Eq. (37) provides an upper bound on flux
noise.

##  VII Paramagnon flux density in the presence of spin diffusion from an
infinite plane of spins

Before we display numerical calculations with our explicit expression for
$D_{ij}$ shown in Eq. (17), it is of value to consider a simplified model for
$D_{ij}$ based on the third-principles theory. Assume the spin system is an
infinite square lattice with no vacancies (homogeneous). Moreover, assume the
superconducting wire is also infinite with flux vector given by the “edge
model”shown in Fig. 1:

$\bm{F(\bm{r})}=F_{0}(\delta_{x,-W/2}-\delta_{x,W/2})\hat{\bm{z}},$ (38)

where $W$ is the wire width and $\delta_{x,\pm W/2}$ are Kronecker delta
functions. Equation (38) gives a good description of thin-film wires where it
is shown that $\bm{F}_{i}$ is sharply peaked at the wire edges [25].

The simplest model for $D_{ij}$ is to emulate the third-principles theory. To
do this, assume the phenomenological paramagnon relaxation rates are given by
$\gamma_{\bm{q}}=D(T)q^{2}+\bar{\Gamma}$, with $D(T)$ the spin diffusion
constant. Inspection of Eq. (31) shows that this $\gamma_{\bm{q}}$ is obtained
by a choice of $D_{ij}$ that has the following Fourier representation at low
$q$:

$\tilde{D}(\bm{q})=-\frac{D(T)q^{2}}{4\left(k_{B}T-J\right)}.$ (39)

Plug these into Eq. (36) along with $\tilde{\chi}(\bm{q},0)\approx
1/[4(k_{B}T-J)]$ and
$|\tilde{\bm{F}}(\bm{q})|^{2}=4F_{0}^{2}N_{sy}^{2}\sin^{2}{(q_{x}W/2)}\delta_{q_{y},0}$
to get

$\displaystyle\rho_{\Phi}(\gamma)$ $\displaystyle=$
$\displaystyle\frac{4F_{0}^{2}N_{sy}^{2}}{N_{s}\gamma}\int_{-\pi/a_{0}}^{\pi/a_{0}}\frac{dq}{\frac{2\pi}{N_{sx}a_{0}}}\sin^{2}{\left(\frac{q_{x}W}{2}\right)}\frac{\bar{\Gamma}+Dq^{2}}{4(k_{B}T-J)}\delta\left(\gamma-
Dq^{2}-\bar{\Gamma}\right)$ (40) $\displaystyle=$
$\displaystyle\frac{F_{0}^{2}a_{0}N_{sy}}{2\pi(k_{B}T-J)\sqrt{D}}\sin^{2}{\left[\sqrt{\frac{(\gamma-\gamma_{{\rm
min}})W^{2}}{4D}}\right]}\frac{\theta(\gamma_{{\rm
max}}-\gamma)\theta(\gamma-\gamma_{{\rm min}})}{\left(\gamma-\gamma_{{\rm
min}}\right)^{1/2}},$

where $\theta(x)$ is the Heaviside step function, with $\gamma_{{\rm
min}}=\bar{\Gamma}$ and $\gamma_{{\rm max}}=D(\pi/a_{0})^{2}+\bar{\Gamma}$
defining the region where $\rho_{\Phi}(\gamma)$ is nonzero. The sine squared
represents interference between the two edges of the wire; in most cases this
averages out to $1/2$ either because of small fluctuations in wire shape or
frequency resolution during integration over $\gamma$. Apart from this,
$\rho_{\Phi}(\gamma)$ follows a power law in frequency

$\rho_{\Phi}(\gamma)=\frac{C(T)\theta(\gamma_{{\rm
max}}-\gamma)\theta(\gamma-\gamma_{{\rm min}})}{\gamma^{\alpha}},$ (41)

with exponent $\alpha=1/2$ and amplitude $C(T)$.

Plug Eq. (41) into Eq. (34) to get

$\tilde{S}_{\Phi}(\omega)=\frac{2\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}C(T)\frac{b_{\alpha}(\omega)}{\omega^{\alpha}},$
(42)

where

$b_{\alpha}(\omega)=\int_{\gamma_{{\rm min}}/\omega}^{\gamma_{{\rm
max}}/\omega}dx\frac{x^{1-\alpha}}{1+x^{2}}.$ (43)

Thus, when $\gamma_{{\rm min}}\ll\omega\ll\gamma_{{\rm max}}$,
$b_{\alpha}(\omega)\approx\pi/[2\sin{(\pi\alpha/2)}]$ and the flux noise
scales as a power law in frequency with the same exponent $\alpha$,
$\tilde{S}_{\Phi}(\omega)\propto 1/\omega^{\alpha}$.

We emphasize that the $\alpha=1/2$ obtained in this section is a direct
consequence of assuming long-wavelength diffusion in an infinite and
homogeneous (spin density $\sigma=1$) lattice of spins. The next section shows
explicit numerical calculations of the impact of spatial confinement and
disorder using Eq. (17) for $D_{ij}$.

(a)

(b)

Figure 2: Explicit calculations of paramagnon flux density
$\rho_{\Phi}(\gamma)$ for a confined system of spins ($20\times 20$ virtual
lattice with open b.c. along $x$), using the edge model for the flux vector
(Eq. (38) and Fig. 1). Results are shown for $k_{B}T/|J|=12$ and
$\Gamma_{i}=0$ for (a) $J>0$ (FM) and (b) $J<0$ (AFM), for spin densities
$\sigma=1,0.75,0.5$. The cases with $\sigma<1$ were averaged over 512 random
instances, each containing a different distribution of vacancies in the
virtual lattice with the same $N_{s}=\sigma N$. These plots demonstrate that
$\rho_{\Phi}(\gamma)$ remains a power law in $\gamma$, even in the presence of
confinement (open b.c.) and disorder ($\sigma<1$). However, the exponent
$\alpha>1/2$, in contrast to $\alpha=1/2$ in the absence of confinement and
disorder.

##  VIII Numerical evaluation of the paramagnon flux density: Heisenberg model
in the 2d square lattice with a random distribution of vacancies

For numerical calculations with our proposed $D_{ij}$ in Eq. (17), we consider
a $20\times 20$ ($N=400$) “virtual” square lattice. We randomly populate
$N_{s}\leq N$ sites with spins, yielding spin density $\sigma=N_{s}/N$. The
remaining unoccupied virtual sites are called vacancies, see Fig. 1. All
calculations below are done with open boundary condition (b.c.) along $x$, and
periodic b.c. along $y$, describing spins confined within the region of the SC
wire.

All calculations assume the nearest-neighbour (n.n.) Heisenberg model in a
square lattice with $J_{ij}=J\sum_{\bm{v}}\delta_{i,j+\bm{v}}$ for $\bm{v}=\pm
a_{0}\bm{\hat{x}},\pm a_{0}\bm{\hat{y}}$. Divide matrix (23a) by
$d_{0}(T)|J|/\hbar$, so that the eigenvalues $\gamma_{m}$ are expressed in
units of $d_{0}(T)|J|/\hbar$. That way the problem is now specified by
$k_{B}T/|J|$, $\hbar\Gamma_{i}/(d_{0}(T)|J|)$, and
$\hbar\omega/(d_{0}(T)|J|)$. As a result, when $\Gamma_{i}=0$, the calculated
exponent $\alpha$ is independent of the choice of $d_{0}(T)$ (although
$\gamma_{{\rm min}},\gamma_{{\rm max}}$ do depend on $d_{0}(T)$).

Figure 2 shows explicit calculations of $\rho_{\Phi}(\gamma)$ using Eq. (35)
with numerical calculations of the eigenvalues $\gamma_{m}$ of $\bm{{\sf P}}$,
with the $\delta(x)$ function approximated by a Gaussian with standard
deviation $0.1\gamma_{max}$. The flux vector was given by the “edge model”,
see Eq. (38) and Fig. 1. Results for $k_{B}T/|J|=12$ and $\Gamma_{i}=0$ are
shown for both $J>0$ (FM) and $J<0$ (AFM), using spin densities
$\sigma=1,0.75,0.5$. The latter two are averaged over 512 random instances,
each containing a different distribution of vacancies in the virtual lattice
with the same $N_{s}=\sigma N$. These results demonstrate that
$\rho_{\Phi}(\gamma)$ remains a power law in $\gamma$, even in the presence of
confinement (open b.c.) and spatial disorder ($\sigma<1$). However, the
corresponding exponent $\alpha>1/2$, contrasting to the case without
confinement/disorder (Section VII).

Figure 3 shows explicit calculations of the frequency exponent $\alpha$
appearing in both $\rho_{\Phi}(\gamma)$ and $\tilde{S}_{\Phi}(\omega)$ as a
function of $T$ and $\sigma$, with all other parameters like in Fig. 2. It is
seen that $\alpha$ _decreases_ with $T$ for the FM model, and has the opposite
behavior for the AFM model. At high $T$, both FM/AFM models lead to
$\alpha\approx 0.7$ for $\sigma=1$; this demonstrates the importance of
confinement. As $\sigma$ decreases from $1$, $\alpha$ further deviates from
its infinite/homogeneous value of $1/2$. This demonstrates the impact of
disorder. Note how the dependence of $\alpha$ on $\sigma$ is nonmonotonic.

Interestingly, several experiments with niobium devices measure $\alpha=0.7$
[29, 4].

(a)

(b)

Figure 3: Explicit calculations of flux noise frequency exponent $\alpha$ as a
function of temperature for $\Gamma_{i}=0$ in the n.n. Heisenberg model with
confinement and disorder, with all other parameters as in Fig. 2. (a) FM case
with $J>0$, (b) AFM case with $J<0$. Note how $\alpha$ decreases (increases)
with $T$ for the FM (AFM) cases. In all cases $\alpha>1/2$ demonstrating the
relevance of confinement and disorder. The dependence of $\alpha$ on $\sigma$
is nonmonotonic.

(a)

(b)

Figure 4: Noise frequency exponent $\alpha$ in the presence of both spin-spin
interaction $J$ and $\Gamma_{i}=\Gamma_{{\rm max}}\exp{\left(-\lambda_{{\rm
max}}r_{i}\right)}$, with $r_{i}\in[0,1]$ a random number for each site $i$.
The spin density is $\sigma=1$ for all calculations (no vacancies), so the
only source of disorder is the variations in $\Gamma_{i}$; all other
parameters are the same as in Fig. 3. The parameters $\Gamma^{\prime}_{{\rm
max}}$ are renormalized relaxation cut-offs, assumed to be constant (see
text). Results show that the introduction of $\Gamma_{i}$ to the interacting
spin system increases the values of $\alpha$ in the low to intermediate $T$
range, with the impact of $\Gamma_{i}$ washed out at high $T$.

##  IX Disorder due to wide distribution of relaxation rates $\Gamma_{i}$ and
spin-spin interactions

Another model for disorder is to assume an arbitrary distribution of
relaxation rates $\Gamma_{i}$. When $J_{ij}=D_{ij}=0$, $\rho_{\Phi}(\gamma)$
and $\tilde{S}_{\Phi}(\omega)$ can be computed exactly. The eigenvalues of
$\bm{{\sf P}}$ are simply $\gamma_{i}=\Gamma_{i}$, and the eigenvectors are
the unit column vectors $\hat{\bm{e}}_{i}$. Using Eq. (35) the paramagnon flux
density becomes

$\displaystyle\rho_{\Phi}(\gamma)$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma}\sum_{i}\left|\bm{F}_{i}\right|^{2}\Gamma_{i}\tilde{\chi}_{0}\delta\left(\gamma-\Gamma_{i}\right)$
(44) $\displaystyle=$
$\displaystyle\frac{1}{4k_{B}T}\sum_{i}\left|\bm{F}_{i}\right|^{2}\int d\Gamma
p(\Gamma)\frac{\Gamma}{\gamma}\delta\left(\gamma-\Gamma\right)$
$\displaystyle=$
$\displaystyle\frac{1}{4k_{B}T}\sum_{i}\left|\bm{F}_{i}\right|^{2}p(\gamma),$

where in the second line we plugged Eq. (29) and took an average using
$p(\Gamma)$, the probability density for rates $\Gamma_{i}$. The justification
for wide distributions of spin-flip rates $\Gamma_{i}$ for spin impurities is
given in [5, 20]. Each spin interacts with one or more amorphous two-level
systems (TLSs), leading to the cross-relaxation rate (joint spin flip and TLS
switch):

$\Gamma_{i}=\Gamma_{{\rm max}}(T)e^{-\lambda},$ (45)

where $\Gamma_{{\rm max}}(T)$ is a cut-off for $\Gamma_{i}$, and $\lambda$ is
a random variable uniformly distributed in the interval $[0,\lambda_{{\rm
max}}]$; it models the barrier for TLS switch, see Fig. 3(a) of [5]. Such a
model is described by the probability density

$p(\Gamma)=\frac{1}{\lambda_{{\rm
max}}}\frac{1}{\left|\frac{d\Gamma}{d\lambda}\right|}=\frac{1}{\lambda_{{\rm
max}}}\frac{1}{\Gamma},$ (46)

for $\Gamma_{{\rm min}}(T)<\Gamma<\Gamma_{{\rm max}}(T)$, and zero otherwise,
where $\Gamma_{{\rm min}}=\Gamma_{{\rm max}}e^{-\lambda_{{\rm max}}}$. Plug
this into Eq. (44) and use Eq. (34) to get the flux noise

$\tilde{S}_{\Phi}(\omega)=\frac{\pi\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\frac{\sum_{i}\left|\bm{F}_{i}\right|^{2}}{4k_{B}T\lambda_{{\rm
max}}}\frac{1}{\omega},$ (47)

for $\Gamma_{{\rm min}}<\omega<\Gamma_{{\rm max}}$, with constant
$\tilde{S}_{\Phi}(\omega)=\tilde{S}_{\Phi}(\Gamma_{{\rm min}})$ for
$\omega<\Gamma_{{\rm min}}$, and $\tilde{S}_{\Phi}(\omega)=0$ for
$\omega>\Gamma_{{\rm max}}$. From now on we will refer to this model as the
_spin $1/f$ model_, since it is the spin equivalent of the well-known $1/f$
noise model in semiconducting devices.

Now consider the impact of nonzero spin-spin interaction $J_{ij}$ and
dissipation $D_{ij}$. Figure 4 shows numerical calculations using the second-
principles theory, for $\sigma=1$ and other parameters as in Fig. 3, plus the
spin $1/f$ model with different choices of $\Gamma^{\prime}_{{\rm
max}}=\hbar\Gamma_{{\rm max}}(T)/(d_{0}(T)|J|)$. In Eq. (47) we see that the
$T$ dependence in $\Gamma_{{\rm max}}$ does not affect the noise provided that
$\Gamma_{{\rm min}}<\omega<\Gamma_{{\rm max}}$; this occurs because low
frequency noise is independent on the cut-off for $\Gamma_{i}$. The same
argument applies in the presence of $J_{ij}$, provided $\Gamma^{\prime}_{{\rm
max}}$ is a sufficiently large constant in our calculations. We simulated the
spin $1/f$ model by choosing $\Gamma_{i}=\Gamma_{{\rm
max}}\exp{\left(-\lambda_{{\rm max}}r_{i}\right)}$ with $\lambda_{{\rm
max}}=20$ and $r_{i}\in[0,1]$ a random number generated for each of the 400
sites.

The calculations should be compared to the case of $J_{ij}=0$, that has
$\alpha=1$ (Eq. (47)). As shown in Fig. 4, the addition of a wide distribution
of $\Gamma_{i}$ to the interacting spin system increases $\alpha$ in the low
to intermediate $T$ range. At high $T$, the impact of $\Gamma_{i}$ is washed
out.

##  X Simple expressions for flux noise for comparison to experiments

The previous sections showed that the paramagnon flux density scales as a
power law in $\gamma$, $\rho_{\Phi}(\gamma)=C/\gamma^{\alpha}$, for
$\gamma_{{\rm min}}\ll\gamma\ll\gamma_{{\rm max}}$, and
$\rho_{\Phi}(\gamma)\approx 0$ outside this range. The exponent $\alpha$ and
amplitude $C$ depend on temperature, confinement, and disorder.

As discussed above Eq. (30), when $\sigma\lesssim 1$ and $\Gamma_{i}$ is
nearly uniform we can use the homogeneous approximation (HA) described in
Appendix A to obtain simpler analytical expressions. The HA replaces $x_{i}$,
$\Gamma_{i}$, $J_{ij}$, and $D_{ij}$ by their spatial averages, $\sigma$,
$\bar{\Gamma}$, $\bar{J}_{ij}$, and $\bar{D}_{ij}$, making the problem
analytically solvable. In particular, $\rho_{\Phi}(\gamma)$ is approximated by
the translation-invariant Eq. (36). Consider its integral over all $\gamma$:

$\displaystyle\int d\gamma\rho_{\Phi}(\gamma)$ $\displaystyle=$
$\displaystyle\frac{1}{N_{s}}\sum_{\bm{q}}\left|\tilde{\bm{F}}(\bm{q})\right|^{2}\frac{1}{\gamma_{\bm{q}}}\left[\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right]$
(48) $\displaystyle=$
$\displaystyle\frac{1}{N_{s}}\sum_{\bm{q}}\frac{\left|\tilde{\bm{F}}(\bm{q})\right|^{2}}{4k_{B}T-\tilde{J}(\bm{q})},$
$\displaystyle=$
$\displaystyle\frac{\sigma\sum_{i}\left|\bm{F}_{i}\right|^{2}}{4\left(k_{B}T-\sigma\bar{J}\right)},$

where in the second line we used Eqs. (31) and (32), and in the third line we
assumed the typical scale for variations in $\bm{F}(\bm{r})$ is much larger
than the lattice spacing $a_{0}$, so that $\tilde{J}(\bm{q})$ can be
approximated by $\tilde{J}(0)=4\sigma\bar{J}$. Now equate this to $\int
d\gamma C/\gamma^{\alpha}\approx C\gamma_{{\rm max}}^{1-\alpha}/(1-\alpha)$
(valid for $\alpha<1$) in order to obtain an expression for amplitude $C$;
plug this into Eq. (42) to get

$\tilde{S}_{\Phi}(\omega)=\frac{\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\frac{\sigma(1-\alpha)\sum_{i}\left|\bm{F}_{i}\right|^{2}}{2\gamma_{{\rm
max}}^{1-\alpha}\left(k_{B}T-\sigma\bar{J}\right)}\frac{b_{\alpha}(\omega)}{|\omega|^{\alpha}}.$
(49)

When $\alpha\lesssim 1$, the temperature dependence of $\gamma_{{\rm
max}}\propto d_{0}(T)T$ is washed out, leading to a simple relation valid for
$\hbar\omega\ll k_{B}T$:

$\frac{\tilde{S}_{\Phi}(\omega,T)}{\tilde{S}_{\Phi}(\omega,T\gg\sigma\bar{J}/k_{B})}\approx\frac{k_{B}T}{k_{B}T-\sigma\bar{J}}.$
(50)

Note how within HA and for $\alpha\lesssim 1$ this noise amplitude ratio is
independent of $d_{0}(T)$ and other details such as values of $\Gamma_{i}$.

Figure 5 shows this dependence for the FM ($\bar{J}>0$) and AFM ($\bar{J}<0$)
cases, and compares Eq. (50) to numerical evaluation with $T=10\bar{J}/k_{B}$
and $\hbar\omega/|\bar{J}|=0.1$ chosen as the high temperature and frequency,
respectively, and $d_{0}(T)=1$, $\Gamma_{i}=0$. It shows that the amplitude of
$\tilde{S}_{\Phi}(\omega,T)$ for a given $\omega$ _decreases_ with increasing
$T$ for the FM model. For the AFM model it instead _increases_ with $T$.
Extrapolating to lower temperatures (including $T<0$ for AFM) allows
determination of the value of $\sigma\bar{J}$ that models the spins.

The singularity in $\tilde{S}_{\Phi}(\omega)$ for the FM model as $T$ gets
close to $T_{c}^{{\rm mag}}=\sigma\bar{J}/k_{B}$ is due to the formation of
clusters of spins with nonzero magnetization (short range order) [15, 25]. As
seen in Fig. 5, numerical results deviate from Eq. (50) demonstrating the
effect of confinement, disorder and the value of $\alpha$ ($\alpha<1$ for most
T, see Fig. 3) on the $T$ dependence of spin clusters.

Figure 5: Temperature dependence of flux noise amplitude. Due to the formation
of magnetized spin clusters, the noise amplitude for the FM model has a
singularity at $T=T_{c}^{{\rm mag}}=\sigma|\bar{J}|/k_{B}$. As a result, the
FM $\tilde{S}_{\Phi}(\omega)$ decreases with increasing $T$. For the AFM model
($\bar{J}<0$), $\tilde{S}_{\Phi}(\omega)$ is instead slightly increasing with
$T$. The solid curves are based on Eq. (50), with dashed ones representing
simple extrapolation to determine $\sigma|\bar{J}|$ for the AFM model. Data
points are numerical evaluation with the second principles theory, using
$T=10\bar{J}/k_{B}$ and $\hbar\omega/|\bar{J}|=0.1$ as the high $T$ and low
$\omega$, respectively, and $d_{0}(T)=1$, $\Gamma_{i}=0$.

Similar considerations apply to the antisymmetric flux noise

$\displaystyle\tilde{S}^{-}_{\Phi}(\omega)$ $\displaystyle=$
$\displaystyle\tilde{S}_{\Phi}(\omega)-\tilde{S}_{\Phi}(-\omega)$ (51)
$\displaystyle=$
$\displaystyle\frac{\sigma(1-\alpha)\sum_{i}\left|\bm{F}_{i}\right|^{2}}{2\gamma_{{\rm
max}}^{1-\alpha}\left(k_{B}T-\sigma\bar{J}\right)}\hbar|\omega|^{1-\alpha}b_{\alpha}(\omega).$

This is seen to have weak frequency dependence when $\alpha\lesssim 1$.
Plotting $1/\tilde{S}^{-}_{\Phi}(\omega)$ as a function of $T$ yields a
straight line that extrapolates to zero at $T=\sigma\bar{J}<0$. In Quintana et
al. [11] (inset of Fig. 3) this procedure reveals $\sigma\bar{J}=-10$ mK.

##  XI Discussion and conclusions

Several experiments in SQUID-based devices [29, 15, 12, 18, 11, 4] have
concluded that flux noise follows the empirical law

$\tilde{S}_{\Phi}(\omega)=\frac{A}{\omega^{\alpha}},$ (52)

where amplitude $A$ and exponent $\alpha<1$ are both temperature-dependent.
Here we developed a theory of flux noise due to interacting spins that is able
to calculate $A$ and $\alpha$ for realistic model impurity spin systems with
disorder due to vacancies and wide distributions of spin relaxation rates, and
for spins confined in bounded regions such as SC wires.

To achieve this we needed to develop a method that does not rely on the
“third-principles” diffusion operator $D\nabla^{2}$. Our “second-principles”
method instead assumes lattice sites are coupled by a dissipation matrix
$D_{ij}$ given by Eq. (17). We showed that this choice obeys fundamental
principles such as total spin conservation and the second law of
thermodynamics. Our prescription for $D_{ij}$ depends on spin Hamiltonian
parameters such as exchange interaction, establishing a direct connection
between flux noise and model spin Hamiltonians.

A central concept is the interpretation of flux noise in terms of (para)magnon
excitations. Flux noise is shown to be directly related to the density of edge
(para)magnons, the ones that lead to fluctuations at the superconducting wire
edges, where spin flips cause the largest flux changes to the device. While
$D_{ij}$ accounts for the interactions between (para)magnons, the rates
$\Gamma_{i}$ describe interactions between spins and other degrees of freedom
such as phonons, electron gas excitations, and two-level system defects.

Section VII shows that choosing $D_{ij}$ consistent with $D\nabla^{2}$ (third-
principles theory) and assuming an infinite, translation-invariant spin system
leads to temperature-independent noise exponent $\alpha=1/2$, in contradiction
to experiments. In contrast, Sections VIII and IX show that numerical
calculations with the “second principles” prescription for $D_{ij}$ (Eq. (17))
makes $\alpha$ temperature-dependent and in the range observed in experiments.
Both $D_{ij}$ and $\Gamma_{i}$ seem to be required to explain experiments.

We can separate experiments in two groups, the ones using niobium and the ones
using aluminum devices. Measurements in niobium devices show both amplitude
$A$ and exponent $\alpha$ decreasing with increasing $T$ (Fig. 3 of [15]).
According to our Eq. (50) and Fig. 5 this $A(T)$ requires a ferromagnetic
model ($J>0$). Our calculated $\alpha(T)$ for $\Gamma_{i}=0$ ranges from $0.9$
at low $T$ to $0.7$ at high $T$ for a confined system (spins only on top of
the wire) with a small number of vacancies ($\sigma\lesssim 1$), but the
exponent can approach $0.5$ when the spin system extends beyond the wire (Eq.
(40)). In comparison, experimental measurement shows $\alpha(T)$ going from
$0.8$ to $0.4$ with increasing $T$ [15]. Therefore, a n.n. Heisenberg model
with $J>0$ and $\sigma\lesssim 1$ and low $\Gamma_{i}\ll J/\hbar$ provides a
reasonable model.

Much less data exists for aluminum devices, but one experiment definitely
shows that their $A(T)$ and $\alpha(T)$ are qualitatively different from
niobium. In [11], $A(T)$ is shown to increase with $T$, and based on our Fig.
5 an antiferromagnetic model ($J<0$) is required. The authors reached the same
conclusion by measuring asymmetric noise $\tilde{S}_{\Phi}^{-}(\omega)$ and
extrapolating to $T<0$ to get $T_{{\rm CW}}=-10$ mK. According to our
interpretation this implies $\sigma\bar{J}/k_{B}=-10$ mK. The measured
exponent $\alpha=0.96-1.05$ can not be explained by our theory with
$\Gamma_{i}=0$ (Fig. 3b). However, Fig. 4b shows that introducing a wide
distribution of $\Gamma_{i}$’s (the spin $1/f$ model) makes $\alpha\approx 1$
at low temperatures. This indicates the necessity of a model with both
$\Gamma_{i}$ and $J$ nonzero in aluminum devices. Measurements of $\alpha(T)$
over a wide temperature range are not yet available to confirm this scenario.

In addition to low frequency flux noise, experiments also measure Ohmic
($\propto\omega$) [17] or super-Ohmic ($\propto\omega^{3}$) [18, 11] flux
noise in the GHz range. We now argue that this behaviour can not arise from
interacting spins alone.

Our theory is fundamentally based on the assumption of “hydrodynamics”, i.e.
that spin degrees of freedom can be described by the classical equation of
motion (13) [30]. As a result, it overestimates the noise for
$\omega>\Omega_{c}$, where $\Omega_{c}$ is a high-frequency cut-off. The cut-
off $\Omega_{c}$ can be estimated from exact calculations of the moments of
the noise spectrum at $T\rightarrow\infty$. Calculations for the 3d Heisenberg
model [22] suggests $\hbar\Omega_{c}\sim 10J/\hbar$ for our 2d case. When
$\omega>\Omega_{c}$, $\tilde{S}_{ij}(\omega)$ drops off faster than
$1/\omega^{2}$, so that our Eq. (33) becomes an upper bound for flux noise. In
Eq. (37) this upper bound was shown to be $\propto 1/\omega$ when quantum
noise is included. As the $1/\omega$ upper bound holds for all interacting
spin models, this allows us to conclude that the high frequency Ohmic [17] or
super-Ohmic [18, 11] flux noise observed in SQUIDs _can not originate from a
model of interacting impurity spins_. A likely source is the normal resistance
due to excited quasiparticles, either by thermal or nonequilibrium sources
such as cosmic rays [31].

In conclusion, we developed a “second principles” theory of flux noise due to
interacting spins that is able to account for the confinement and disorder
present in realistic impurity spins systems on superconducting devices. The
theoretical framework allows explicit prediction of the amplitude and exponent
of flux noise due to different wire geometries and spin disorder scenarios,
such as random vacancies and wide distributions of spin-flip rates due to
interactions with amorphous TLSs. Comparing numerical results to experiments
allowed us to specify different spin Hamiltonians for niobium and aluminum
devices. Generalizations of the theory that include time dependent external
currents and fields can be used to design of optimal control strategies that
reduce the impact of flux noise on quantum devices.

###### Acknowledgements.

This work was supported by NSERC (Canada) through its Discovery program (Grant
number RGPIN-2020-04328). The authors thank M. Amin, R. Harris, P. Kovtun, and
T. Lanting for useful discussions.

##  Appendix A Exact solution of the homogeneous case and homogeneous
approximation

When the spin system is translation-invariant, the “second principles” method
yields exact analytic expressions for the spin noise for general $D_{ij}$
using spatial Fourier transforms. This is the case when the system has
periodic b.c., and the vacancies are organized in a regular sublattice of the
full virtual lattice. The homogeneous case (no vacancies, $\sigma=1$) with
periodic b.c. is the most relevant example.

When the system is not translation invariant, we may take spatial averages
over parameters $x_{i}$, $D_{ij}$, $J_{ij}$, and $\Gamma_{i}$ in order to
force its equation of motion to become exactly solvable. In this case the
resulting analytic solution is called homogeneous approximation (HA).

In the paramagnetic phase ($T>T_{c}^{{\rm mag}}$), the EOM Eq. (21) becomes

$\displaystyle\frac{d}{dt}(\delta s_{i}^{a})$ $\displaystyle=$
$\displaystyle-\sum_{j}D_{ij}x_{j}\left[(\delta
h_{j}^{a})+\sum_{k}x_{k}J_{jk}(\delta s_{k}^{a})\right.$
$\displaystyle\left.-4k_{B}T(\delta s_{j}^{a})\right]-\Gamma_{i}(\delta
s_{i}^{a})+\Gamma_{i}\sum_{j}\tilde{\chi}_{ij}(0)(\delta h_{j}^{a}),$

for $a=x,y,z$. If this is not translation invariant, replace
$x_{i},D_{ij},\Gamma_{i}$ by their average values:

$\displaystyle\bar{x}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i}x_{i}=\sigma,$ (54a)
$\displaystyle\bar{D}_{\bm{v}}$ $\displaystyle=$
$\displaystyle\frac{1}{N\sigma^{2}}\sum_{i}x_{i}x_{i+\bm{v}}D_{i,i+\bm{v}},$
(54b) $\displaystyle\bar{J}_{\bm{v}}$ $\displaystyle=$
$\displaystyle\frac{1}{N\sigma^{2}}\sum_{i}x_{i}x_{i+\bm{v}}J_{i,i+\bm{v}},$
(54c) $\displaystyle\bar{\Gamma}$ $\displaystyle=$
$\displaystyle\frac{1}{N\sigma}\sum_{i}x_{i}\Gamma_{i},$ (54d)
$\displaystyle\bar{\chi}_{\bm{v}}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{N\sigma^{2}}\sum_{i}\tilde{\chi}_{i,i+\bm{v}}(0).$
(54e)

Equation (21) becomes

$\displaystyle\frac{d}{dt}(\delta s_{i}^{a})$ $\displaystyle=$
$\displaystyle-\sum_{\bm{v}}\sigma\bar{D}_{\bm{v}}\left[(\delta
h_{i+\bm{v}}^{a})+\sum_{\bm{v}^{\prime}}\sigma\bar{J}_{\bm{v}^{\prime}}(\delta
s_{i+\bm{v}+\bm{v}^{\prime}}^{a})\right.$ $\displaystyle\left.-4k_{B}T(\delta
s_{i+\bm{v}}^{a})\right]-\bar{\Gamma}(\delta
s_{i}^{a})+\bar{\Gamma}\sum_{\bm{v}}\bar{\chi}_{\bm{v}}(0)(\delta
h_{i+\bm{v}}^{a}).$

Take the Fourier transform in both sides

$\delta\tilde{s}^{a}_{\bm{q}}(\omega)=\int
dt\sum_{j}e^{-i\left(\bm{q}\cdot\bm{R}_{j}-\omega t\right)}(\delta
s_{j}^{a}),$ (56)

to obtain the dynamical susceptibility,

$\tilde{\chi}^{ab}(\bm{q},\omega)=\frac{\left[\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right]\delta_{ab}}{-i\omega+\left\\{-\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]+\bar{\Gamma}\right\\}},$
(57)

where

$\displaystyle\tilde{D}(\bm{q})$ $\displaystyle=$
$\displaystyle\sigma\sum_{\bm{v}}\bar{D}_{\bm{v}}e^{-i\bm{q}\cdot\bm{v}},$
(58a) $\displaystyle\tilde{J}(\bm{q})$ $\displaystyle=$
$\displaystyle\sigma\sum_{\bm{v}}\bar{J}_{\bm{v}}e^{-i\bm{q}\cdot\bm{v}}.$
(58b)

This result implies the zero-frequency susceptibility,

$\tilde{\chi}^{ab}(\bm{q},\omega=0)=\frac{1}{4k_{B}}\frac{\delta_{ab}}{T-T_{{\rm
CW}}(\bm{q})},$ (59)

with Curie-Weiss temperature $T_{{\rm
CW}}(\bm{q})=\frac{1}{4k_{B}}\tilde{J}(\bm{q})$.

Using the fluctuation-dissipation theorem the spin noise in Fourier space
becomes exactly equal to

$\displaystyle\tilde{S}^{ab}(\bm{q},\omega)$ $\displaystyle=$
$\displaystyle\frac{2\hbar\omega}{1-e^{-\hbar\omega/k_{B}T}}$
$\displaystyle\times\frac{\left[\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right]\delta_{ab}}{\omega^{2}+\left\\{-\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]+\bar{\Gamma}\right\\}^{2}}.$

In the notation of Eq. (27b), the paramagnon modes are labelled by
$m=\bm{q}\in 1^{{\rm st}}$ Brillouin zone, each with frequency eigenvalue

$\gamma_{m}=\gamma_{\bm{q}}=-\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]+\bar{\Gamma},$
(61)

and right and left eigenvectors $\bm{v}_{m}=\bm{e}_{\bm{q}}$,
$\bm{v}_{m}^{-1}=\bm{e}_{\bm{q}}^{{\dagger}}/N_{s}$, respectively, where
$\bm{e}^{{\dagger}}_{\bm{q}}=\left(e^{-i\bm{q}\cdot\bm{R}_{0}},\ldots,e^{-i\bm{q}\cdot\bm{R}_{N-1}}\right)$.

Note that we assumed the presence of $N_{s}$ $\bm{R}_{i}$’s forming a
translation-invariant lattice, so there are $N_{s}$ $\bm{q}$’s in the $1^{{\rm
st}}$ Brillouin zone. As a consequence Eq. (A) does not have $\sigma$
appearing explicitly in the numerator.

Similar to what was done in Section VI, we can define a mode density to
interpret Eq. (A),

$\displaystyle\tilde{S}^{ab}(\bm{q},\omega)=\frac{2\pi\hbar\omega}{1-e^{-\frac{\hbar\omega}{k_{B}T}}}\int
d\gamma\frac{\gamma/\pi}{\omega^{2}+\gamma^{2}}\rho^{ab}(\gamma,\bm{q}).$ (62)

The $\rho^{ab}(\gamma,\bm{q})$ is called paramagnon wavevector density (in
contrast to the paramagnon flux density defined by Eq. (35)). It’s given by

$\rho^{ab}(\gamma,\bm{q})=\frac{1}{\gamma}\left[\bar{\Gamma}\tilde{\chi}(\bm{q},0)-\tilde{D}(\bm{q})\right]\delta(\gamma-\gamma_{\bm{q}})\delta_{ab}.$
(63)

For the n.n. Heisenberg model in the 2d square lattice we get

$\displaystyle\tilde{D}(\bm{q})$ $\displaystyle=$
$\displaystyle-\frac{d_{0}(T)}{\hbar}\left[\sin^{2}{\left(\frac{q_{x}a_{0}}{2}\right)}+\sin^{2}{\left(\frac{q_{y}a_{0}}{2}\right)}\right],$
(64a) $\displaystyle\tilde{J}(\bm{q})$ $\displaystyle=$ $\displaystyle 2\sigma
J\left[\cos{(q_{x}a_{0})}+\cos{(q_{y}a_{0})}\right].$ (64b)

We remark that $\tilde{D}(\bm{q})$ does not depend on $\sigma$ because from
Eq. (18) $\bar{J}_{c}=4\sigma|J|$, so that $\sigma$ cancels out in the
definition of $\tilde{D}(\bm{q})$. These results lead to the effective
diffusion constant in the HA,

$\displaystyle D_{{\rm hom}}(\bm{q})$ $\displaystyle\equiv$
$\displaystyle-\frac{\tilde{D}(\bm{q})\left[4k_{B}T-\tilde{J}(\bm{q})\right]}{q^{2}}$
$\displaystyle=$
$\displaystyle\frac{d_{0}(T)a_{0}^{2}}{\hbar}\left[k_{B}T-\sigma
J\left(1-\frac{q^{2}a_{0}^{2}}{4}\right)\right]+{\cal O}(q^{4}).$

Finally, in the HA approximation the flux noise with edge flux vector Eq. (38)
is given by

$\tilde{S}_{\Phi}(\omega)=4F_{0}^{2}\frac{N_{sy}}{N_{sx}}\sum_{q}\sin^{2}{\left(\frac{qW}{2}\right)}\tilde{S}^{zz}(q\bm{\hat{x}},\omega),$
(66)

where $q=\frac{2\pi}{N_{sx}a_{0}}\left(n-\frac{N_{sx}}{2}\right)$ with
$n=0,1,\ldots,N_{sx}-1$, where $N_{s}=N_{sx}N_{sy}$ is the number of occupied
sites.

## References

  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandão, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Hutchings  _et al._ [2017] M. D. Hutchings, J. B. Hertzberg, Y. Liu, N. T. Bronn, G. A. Keefe, M. Brink, J. M. Chow, and B. L. T. Plourde, Tunable Superconducting Qubits with Flux-Independent Coherence, Phys. Rev. Appl. 8, 044003 (2017).
  * Chávez-Garcia  _et al._ [2022] J. M. Chávez-Garcia, F. Solgun, J. B. Hertzberg, O. Jinka, M. Brink, and B. Abdo, Weakly Flux-Tunable Superconducting Qubit, Phys. Rev. Appl. 18, 034057 (2022).
  * Zaborniak and de Sousa [2021] T. Zaborniak and R. de Sousa, Benchmarking Hamiltonian Noise in the D-Wave Quantum Annealer, IEEE Trans. Quantum Eng. 2, 3100206 (2021).
  * de Sousa [2007] R. de Sousa, Dangling-bond spin relaxation and magnetic $1/f$ noise from the amorphous-semiconductor/oxide interface: Theory, Phys. Rev. B 76, 245306 (2007).
  * Koch  _et al._ [2007] R. H. Koch, D. P. DiVincenzo, and J. Clarke, Model for 1/f Flux Noise in SQUIDs and Qubits, Phys. Rev. Lett. 98, 267003 (2007).
  * Sendelbach  _et al._ [2008] S. Sendelbach, D. Hover, A. Kittel, M. Mück, J. M. Martinis, and R. McDermott, Magnetism in SQUIDs at Millikelvin Temperatures, Phys. Rev. Lett. 100, 227006 (2008).
  * Faoro and Ioffe [2008] L. Faoro and L. B. Ioffe, Microscopic origin of low-frequency flux noise in Josephson circuits, Phys. Rev. Lett. 100, 227005 (2008).
  * Kumar  _et al._ [2016] P. Kumar, S. Sendelbach, M. A. Beck, J. W. Freeland, Z. Wang, H. Wang, C. C. Yu, R. Q. Wu, D. P. Pappas, and R. McDermott, Origin and Reduction of $1/f$ Magnetic Flux Noise in Superconducting Devices, Phys. Rev. Appl. 6, 041001(R) (2016).
  * de Graaf  _et al._ [2017] S. E. de Graaf, A. A. Adamyan, T. Lindström, D. Erts, S. E. Kubatkin, A. Y. Tzalenchuk, and A. V. Danilov, Direct Identification of Dilute Surface Spins on Al2O3: Origin of Flux Noise in Quantum Circuits, Phys. Rev. Lett. 118, 057703 (2017).
  * Quintana  _et al._ [2017] C. M. Quintana, Y. Chen, D. Sank, A. G. Petukhov, T. C. White, D. Kafri, B. Chiaro, A. Megrant, R. Barends, B. Campbell, Z. Chen, A. Dunsworth, A. G. Fowler, R. Graff, E. Jeffrey, J. Kelly, E. Lucero, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, P. Roushan, A. Shabani, V. N. Smelyanskiy, A. Vainsencher, J. Wenner, H. Neven, and J. M. Martinis, Observation of Classical-Quantum Crossover of $1/f$ Flux Noise and Its Paramagnetic Temperature Dependence, Phys. Rev. Lett. 118, 057702 (2017).
  * Lanting  _et al._ [2014] T. Lanting, M. H. Amin, A. J. Berkley, C. Rich, S.-F. Chen, S. LaForest, and R. de Sousa, Evidence for temperature-dependent spin diffusion as a mechanism of intrinsic flux noise in SQUIDs, Phys. Rev. B 89, 014503 (2014).
  * [13] T. Lanting, M. H. Amin, C. Baron, M. Babcock, J. Boschee, S. Boixo, V. N. Smelyanskiy, M. Foygel, and A. G. Petukhov, Probing Environmental Spin Polarization with Superconducting Flux Qubits, arXiv:2003.14244 .
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nat. Phys. 7, 565 (2011).
  * Anton  _et al._ [2013] S. M. Anton, J. S. Birenbaum, S. R. O’Kelley, V. Bolkhovsky, D. A. Braje, G. Fitch, M. Neeley, G. C. Hilton, H.-M. Cho, K. D. Irwin, F. C. Wellstood, W. D. Oliver, A. Shnirman, and J. Clarke, Magnetic Flux Noise in dc SQUIDs: Temperature and Geometry Dependence, Phys. Rev. Lett. 110, 147002 (2013).
  * Atalaya  _et al._ [2014] J. Atalaya, J. Clarke, G. Schön, and A. Shnirman, Flux $1/f^{\alpha}$ noise in two-dimensional Heisenberg spin glasses: Effects of weak anisotropic interactions, Phys. Rev. B 90, 014206 (2014).
  * Lanting  _et al._ [2011] T. Lanting, M. H. S. Amin, M. W. Johnson, F. Altomare, A. J. Berkley, S. Gildert, R. Harris, J. Johansson, P. Bunyk, E. Ladizinsky, E. Tolkacheva, and D. V. Averin, Probing high-frequency noise with macroscopic resonant tunneling, Phys. Rev. B 83, 180502(R) (2011).
  * Yan  _et al._ [2016] F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, The flux qubit revisited to enhance coherence and reproducibility, Nature Comm. 7, 12964 (2016).
  * Bell  _et al._ [2008] C. Bell, S. Milikisyants, M. Huber, and J. Aarts, Spin Dynamics in a Superconductor-Ferromagnet Proximity System, Phys. Rev. Lett. 100, 047002 (2008).
  * Belli  _et al._ [2020] M. Belli, M. Fanciulli, and R. de Sousa, Probing two-level systems with electron spin inversion recovery of defects at the Si/SiO2 interface, Phys. Rev. Res. 2, 033507 (2020).
  * Bloembergen [1949] N. Bloembergen, On the interaction of nuclear spins in a crystalline lattice, Physica 15, 386 (1949).
  * de Gennes [1958] P. G. de Gennes, Inelastic magnetic scattering of neutrons at high temperatures, J. Phys. Chem. Solids 4, 223 (1958).
  * Bennett and Martin [1965] H. S. Bennett and P. C. Martin, Spin Diffusion in the Heisenberg Paramagnet, Phys. Rev. 138, A608 (1965).
  * Hertz and Klemm [1983] J. A. Hertz and R. A. Klemm, Spin-glass dynamics with conserved magnetization, Phys. Rev. B 28, 3849 (1983).
  * LaForest and de Sousa [2015] S. LaForest and R. de Sousa, Flux-vector model of spin noise in superconducting circuits: Electron versus nuclear spins and role of phase transition, Phys. Rev. B 92, 054502 (2015).
  * Beairsto  _et al._ [2021] S. Beairsto, M. Cazayous, R. S. Fishman, and R. de Sousa, Confined Magnons, Phys. Rev. B 104, 134415 (2021).
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, Cambridge, U.K., 1995).
  * Halperin and Hohenberg [1967] B. Halperin and P. Hohenberg, Generalization of scaling laws to dynamical properties of a system near its critical point, Phys. Rev. Lett. 19, 700 (1967).
  * Wellstood  _et al._ [1987] F. C. Wellstood, C. Urbina, and J. Clarke, Low frequency noise in dc superconducting quantum interference devices below 1 K, Appl. Phys. Lett. 50, 772 (1987).
  * Kadanoff and Martin [1963] L. P. Kadanoff and P. C. Martin, Hydrodynamic equations and correlation functions, Ann. Phys. 24, 419 (1963).
  * Wilen  _et al._ [2021] C. D. Wilen, S. Abdullah, N. A. Kurinsky, C. Stanford, L. Cardani, G. d’Imperio, C. Tomei, L. Faoro, L. B. Ioffe, C. H. Liu, A. Opremcak1, B. G. Christensen, J. L. DuBois, and R. McDermott, Correlated charge noise and relaxation errors in superconducting qubits, Nature 594, 369 (2021).