\section{Introduction}

\textbf{A major issue for the implementation of large scale superconducting quantum circuits  is the interaction with interfacial two-level system defects (TLS) that leads to qubit relaxation and impedes qubit operation in certain frequency ranges that also drift in time.  
Another major challenge comes from non-equilibrium quasiparticles (QPs) that result in qubit dephasing and relaxation. }
\textbf{
In this work we show that such QPs can also serve as a source of TLS. 
Using spectral and temporal mapping of TLS-induced fluctuations in frequency tunable resonators, we identify a subset of the general TLS population that are highly coherent TLS with a low reconfiguration temperature $\sim$ 300 mK, and a non-uniform density of states.
These properties can be understood if these TLS are formed by QPs trapped in shallow subgap states formed by spatial fluctutations of the superconducting order parameter $\Delta$. 
Magnetic field measurements of one such TLS reveals a link to superconductivity. Our results imply that trapped QPs can induce qubit relaxation.}



It is becoming increasingly evident that the most coherence limiting TLS reside outside the qubit junctions
 on the surface of metals and dielectrics surrounding them \cite{klimov2018,lisenfeld2019,burnett2014,degraaf2017,bilmes2019} and their slowly fluctuating dynamics poses a significant challenge for quantum computation \cite{supremacy,klimov2018,lisenfeld2019,schlor2019, burnett2019}.
A wide range of techniques has been developed to study and mitigate such TLS, by developing qubits
and resonators into probes of a wider range of material properties \cite{schneider2019,degraaf2017,degraaf2018,grabovskij2012,lisenfeld2019,geaney2019,degraaf2015, lisenfeld2015,lisenfeld2015b}.
However, charged surface TLS and paramagnetic impurities also result in a stochastic and locally varying backdrop for QPs in the superconductor itself. It is thus important to consider the implications on the ever-present excess number of non-equilibrium QPs in superconducting quantum devices \cite{Henrique2019, jin2015, devoretPRL} which must be eliminated in order to increase qubit dephasing times. It is also likely that trapped QPs
are responsible for the non-equilibrium relaxation of transmons recently observed \cite{devoretPRL}. 
Non-equlibrium QPs can be generated by electromagnetic radiation \cite{devisser2012} and from rare high energy cosmic particles impinging the sample, the latter inducing correlated errors in all qubits of a surface code architecture \cite{moore2012}. 

Here we reveal the implications of a spatially fluctuating $\Delta$ on such QPs. 
Scanning tunneling spectroscopy studies typically find spatial variations in $\Delta$ of the order of $\delta\Delta = 10-20\%$  in moderately
disordered superconductors \cite{carbillet2019,lemarie2013,Liao2019}.
$\delta\Delta$ gets smaller in very clean superconductors, however,
even in exceptionally clean films with negligible intrinsic magnetic
disorder an everpresent surface spin density
of the order of $\sim5\cdot10^{17}$ m$^{-2}$ \cite{degraaf2017, anton2013, saveskul2019} results in both flux noise \cite{degraaf2017, quintana2017, anton2013, kumar2016, paladino2014,muller2017} and spatially
non-uniform gap suppression \cite{saveskul2019}. In thin films the gap is also non-uniformely suppressed
by the Altshuler-Aronov effect  due to impurity scattering \cite{Altshurer_Aronov, carbillet2019}, as well as thickness variations \cite{ivry2014}, particularly relevant for Al \cite{aldelta}. The clean, low resistance NbN films that we chose to study in this work are characterized by a low carrier density \cite{Wong2017}. In such materials it is very likely that the DOS varies significantly, resulting in gap variations. 

Crucially, 
we find that the optimal fluctuation able to trap a QP only has a few bound states. As a result a trapped QP behaves similar to `textbook' TLS, as shown in Fig. \ref{fig:sketch}. We hereafter refer to this type of quasiparticle TLS as `qTLS', which coexist with a large number of conventional glassy TLS.  The spatial extent of the gap fluctuation sets the frequency of the qTLS, which means that the density of states (DOS) of these objects is very non-uniform (in contrast to conventional TLS in glasses). We predict a maximum in the qTLS DOS around 6-10 GHz, the frequency range of most qubits (transmons).
Relatively shallow traps also imply that even a modest temperature can reshuffle the QPs. 
We expect the underlying qualitative physics of charge defects formed by trapped QPs to be applicable to most materials, including
superconductors in the dirty limit, such as Al used for qubits, but the relevant frequency ranges may differ somewhat due to different material parameters.

\begin{figure*}[t!]
\includegraphics[scale=0.6]{fig1.png}
\caption{qTLS formation. a) Illustration of typical fluctuations of the order parameter $\Delta$ and the formation of wells with multiple bound states in which quasiparticles can be trapped, effectively forming a TLS (`qTLS'). The scale of the trap determines the qTLS frequency and escape rate. b) qTLS are separated by large distances compared to length-scales at which fluctuations in $\Delta$ occur, fluctuations which may originate from impurity scattering or surface spins. Importantly the qTLS are separated in distance from the superconductor oxide surface where conventional TLS defects, adsorbants etc. reside, and form the usual bath of TLS.  c) Numerically evaluated typical elongated fluctuation of the order parameter capable of hosting multiple bound states for quasiparticles, forming a qTLS. Here $\tilde{\xi}=\sqrt{\Delta_0/E} \xi_0$, the scaling assumes $(k_F \xi_0)^{-1}= 0.14$, $\Delta_0/E = 40$, appropriate for our NbN samples.  }
\label{fig:sketch} 
\end{figure*}

This model is supported by a series of experimental observations from detailed mapping of parameter fluctuations in high-Q frequency tunable superconducting NbN resonators \cite{sumedh2019}. 
Such resonators allow us to map individual TLS in energy and trace them in a broad range of external parameters, such as magnetic or electric field \cite{lisenfeld2019}.
Through analysis of a large number of individual TLS we find:
(i) most qTLS are highly coherent, some with detected linewidths down to $\lesssim 70$ kHz, the resonator linewidth. (ii) The qTLS landscape irreversibly reshuffles at very modest temperatures of $\sim 300$ mK,  (iii) The qTLS DOS appears strongly suppressed at lower frequencies. And (iv) tracking a qTLS in magnetic field reveals a qTLS energy that scale in a way similar to the expected scaling of $\Delta$.
In particular, the presence of a low energy scale  $\lesssim1\text{K},$ is inconsistent with conventional TLS theory in which all scales
are set by chemical energies \cite{phillips}. The low energy scale is instead natural for TLS formed by QPs trapped in regions of locally smaller gap. 
We start by outlining the main features of the model that explains these observations.

\addtocontents{toc}{\protect\setcounter{tocdepth}{-1}}
\section*{Results}
To understand the emergence of qTLS from a fluctuating $\Delta = \Delta_0 + \delta\Delta({\bf r})$
we consider a general model that assumes gaussian
fluctuations of the gap, $\left\langle \delta\Delta({\bf r})\delta\Delta({\bf r'})\right\rangle =g\delta({\bf r}-{\bf r'})$, where $g\simeq \delta \Delta\Delta_0\xi_0^2$ is the strength of fluctuations, and $\xi_0$ the BCS coherence length. 
The exact origin of these fluctuations is not relevant
for the conclusions. The density of subgap states is determined by
 rare fluctuations of the gap that reduce the QP 
energy sufficiently below the gap edge in a uniform material. 
The key point is that such fluctuations  can create traps with a shape and depth that allow them to effectively trap QPs.
As the QPs near the gap edge move with momenta close to the Fermi momentum $p_{F}$, it follows that trapped QPs will have similar momenta, which can point in an arbitrary direction.
This freedom results in the optimal trap being very anisotropic, with a shape elongated along the direction
of $p_{F}$. In Fig. \ref{fig:sketch}c we show the typical shape of such a trap obtained through numerical simulations.
In general, the probability to find a fluctuation of depth E scales exponentially with the area $A$ of the fluctuation, $P\sim \exp{(-E^2A)}$. Efficient traps thus have small area, however, a too small dimension along the direction of $p_F$ will forbid the formation of a bound state due to the large kinetic energy of the QP. Therefore optimal traps do not favor isotropic fluctuations. 
The wave function of the QP oscillates quickly along the direction of $p_F$. The
ground and excited states in the trap differ in the number of oscillations:
since the total number of oscillations is large it is not surprising
that each trap typically contains more than one bound state.


\begin{figure*}[t!]
\includegraphics[scale=0.63]{fig2.png} \caption{Spectral \& temporal fluctuations of the internal quality factor and
spectral reconfiguration of TLS. (a) Initial measurement taken at
10 mK and $\langle n\rangle\approx40$. For each pixel the resonator
line-shape is measured and fitted to extract the internal quality
factor $Q_{i}$. 
The measurement is repeated as the resonator is tuned in frequency (the typical frequency resolution is $\approx 50$ kHz), this measurement
vs detuning is then repeated over time. The color scale shows the normalised
variation $\delta Q_{i}=(Q_{i}-\langle Q_{i}\rangle)/\langle Q_{i}\rangle$. 
(b) same measurement performed after
(a) with the temperature
set to 300 mK. (c) the same measurement was repeated again after returning
to 10 mK. The red arrows indicate a broad TLS jumping into the measurement window. Dashed boxes indicates two strong
TLS. For repeated experiments see \cite{suppl}. d) The time-averaged data from panels (a)-(c) respectively.}
\label{fig:spectraltemporal} 
\end{figure*}


Due to the presence of excited bound states, a QP in such a
trap forms an effective TLS with typical energy splitting $\hbar\omega_{\text{qTLS}}=\zeta\left|\widetilde{E}\right|$
where $\widetilde{E}$ is the energy of the trapped QP as measured
from the gap edge and  $\zeta\approx0.7$ is found from theory \cite{suppl}. 
One consequence of these traps is that, in contrast to conventional TLS, the DOS of qTLS is not constant.
A large $\omega_{qTLS}$ implies very deep traps that are exponentially
rare. We find the density
\begin{equation}
\rho_{\text{qTLS}}(\omega)\simeq \nu(\widetilde{E}) = \frac{1}{\xi_{0}^{2}\Delta_{0}}\text{exp}\left[-\frac{\eta}{g}\frac{\xi_{0}^{2}\Delta_{0}^{2}}{\sqrt{k_{F}\xi_{0}}}\left(\frac{\hbar\omega}{\zeta\Delta_{0}}\right)^{5/4}\right],
\end{equation}
where $\eta\approx 10$. Likewise, small $\omega_{qTLS}$ are suppressed since
QPs in shallow traps anihilate each other efficiently, leading to 
\begin{equation}
\rho_{\text{qTLS}}(\omega) \leq \nu_0(\widetilde{E}) = \frac{\sqrt{2k_{F}\xi_{0}}}{\xi_{0}^{2}\Delta_{0}[ \text{ln}(\tau_{\text{exp}} \Delta_0/\hbar) ]^2}\left({\frac{\hbar\omega}{\zeta\Delta_{0}}}\right)^{3/4},
\end{equation}
in the absence of QP generation. Here $\tau_{\text{exp}}$ is an experimental timescale.


Due to the spatial asymmetry, the transition between the lowest and 
first excited state in each trap is expected to have a significant electric dipole
moment, resulting in a strong interaction with quantum devices.
In this picture the TLS is formed in a single well, in contrast to the more conventional double-well picture.

In similar films  $\delta\Delta$ constitute 10\% of the total gap \cite{carbillet2019, lemarie2013, Liao2019}, which gives in our devices $\delta\Delta\approx2.5$ K. 
At $T=10$ mK tunneling out of such a local minima is exponentially suppressed, however, at temperatures approaching fractions of $\delta\Delta$ there is a finite probability that the QP escapes and either finds another minima or recombines, resulting in a typical escape temperature of the order of $ 200-300$ mK; each event resulting in a different TLS landscape.


The unusual shape of the
traps and the presence of multiple bound states in each trap distinguish clean superconductors from the dirty limit  \cite{Larkin_Ovchinnikov, Lamacraft_Simons, Silva_Ioffe}. In the latter a QP is scattered many times by defects
while moving in a particular direction inside the trap, making
the formation of elongated traps impossible. 
The physics of  qTLS traps is  similar to the well studied phenomena of the appearance of states with negative energy
in disordered conductors \cite{Halperin_Lax66,Cardy1978}, with two important differences: electrons at the bottom of the band in a disordered conductor have much smaller momenta, and the wave functions do not have zeros. 






\begin{figure*}[t!]
\includegraphics[scale=0.54]{fig3.png} \caption{qTLS Properties. a) Histogram of extracted TLS linewidths (FWHM) obtained
at 10 mK across a series of cooldowns and different resonators in
the frequency range 3-7 GHz.  b)  Expected scaling of the qTLS density of states (blue solid region) bounded by  $\nu(\widetilde{E})$ and $\nu_0(\widetilde{E})$ (see text for details). Not to scale. c) Number of strongly coupled TLS as a function of frequency. The data
points are obtained by counting the number of TLS present in a spectral
and temporal map. Data taken on two different samples (circles and triangles) in seven separate cooldowns (different colors). Duplicate markers at the same frequency are before/after thermal reconfiguration. Multiple markers from the same resonator have been offset somewhat in frequency for better visibility. 
Error bars indicate miscounting the number of TLS by one, e.g. if one
TLS happens to be present just outside the measurement range. The solid line is the average of observations, bounded by the standard deviation. The inset shows the same data plotted against the corresponding trap depth from $\Delta_0$. All data obtained at 10 mK and $\langle N\rangle<100$.}
\label{fig:linewidths} 
\end{figure*}



We now turn to the experimental work where our observations can be explained by this model.
To be able to study individual TLS and spectral and temporal fluctuations
of superconducting resonator parameters we use kinetic inductance
tunable superconducting resonators with centre frequencies $f_0$ in the
range 3-7 GHz, described in detail in Ref. \cite{sumedh2019}. The resonators
are cooled down to a base temperature of 10 mK in a well filtered
dilution refrigerator. A small dc current ($<1$ mA) is applied to change the kinetic inductance and tune $f_0$.
We measure the transmitted microwave signal, $S_{21}$, using a vector network analyser and 
extract the resonator frequency and quality factors Q. In some of the
experiments we also utilise a gate electrode mounted in the lid of
the sample enclosure, such that we can affect the TLS energy through
the applied gate voltage $V_g$: $E_{{\rm TLS}}=\sqrt{\epsilon_{0}^{2}+\gamma^2 V_{g}^{2}}$ \cite{bilmes2019,lisenfeld2019}. Here $\epsilon_{0}$ is the TLS minimum energy and $\gamma$ the TLS voltage coupling strength.

In Fig. \ref{fig:spectraltemporal}a we show the relative variation
of the internal quality factor $\delta Q_{i}=(Q_{i}-\langle Q_{i}\rangle)/\langle Q_{i}\rangle$, 
as we tune the resonator frequency. The measurement is continously
repeated over several days to produce a spectral and temporal map
of the variations in loss. This reveals large drops in $\delta Q_{i}$: at detunings of 
$3.5$ and $14$ MHz (from 6.995 GHz) TLS dips remain in the same place whereas at
other frequencies both switching-like behaviour and drift can be observed.
We then raise the temperature of the cryostat to 300 mK and repeat
the same measurement (Fig. \ref{fig:spectraltemporal}b), and return to 10 mK (Fig. \ref{fig:spectraltemporal}c). We observe that a strong TLS 
remains stable over more than 60 hours in each individual measurement,
but thermal cycling to a mere 300 mK makes another TLS  appear at a different frequency.
This is quite remarkable as the temperature is increased to a value
smaller than the energy level splitting of the TLS itself (7 GHz $\approx$
350 mK), and certainly not expected from conventional TLS-glass physics
in which all scales are set by chemical energies \cite{phillips}.
Separate high power measurements show that these dips are saturable, and not caused by e.g. local variations
in background transmission. We have repeated the same experiment but visiting multiple temperatures
on the way to 300 mK (see  \cite{suppl} for additional data) which confirms that reconfiguration only
happens near 300 mK.

\begin{figure*}[t!]
\includegraphics[scale=0.68]{fig4.png} \caption{A single strongly coupled qTLS in magnetic field. a) Transmission
spectrum of the resonator strongly coupled to a qTLS at selected magnetic
fields, as a function of applied electrostatic gate voltage, showing
the coherent coupling of a single qTLS to the resonator. b) The resonance
frequency as a function of magnetic field (solid blue line) together
with the fitted qTLS minimum energy $\epsilon_{0}$ (red markers).
Other extracted parameters are qTLS-resonator coupling strength $g_0=270$
kHz, damping rate $\gamma_{{\rm qTLS}}=200$ KHz and $\gamma=10$ MHz/V, consistent with $p_0\approx 1$ $\rm e\AA$ in our geometry.}
\label{fig:tlsbfield} 
\end{figure*}

Occasionally we observe very dissipative TLS with linewidths in excess of several MHz in our measurement window (red arrow in Fig. \ref{fig:spectraltemporal}c). Here we concentrate instead on the much more coherent TLS that we consistently observe, the narrowest with
line-widths down to $\sim70$ kHz \cite{suppl}, limited by the readout resonator line-width. In Fig. \ref{fig:linewidths}a we show the histogram of
linewidths of the observed TLS, which peaks at values corresponding to coherence times in excess
of several microseconds, much longer than what has been reported for individually resolved TLS in qubits \cite{lisenfeld2019,lisenfeld2015}. The energy leakage caused by  such highly coherent TLS could pose a significant challenge for quantum error correction \cite{Brown2019-1}.
Naturally, qTLS interact with conventional TLS and slow fluctuators that result in parameter fluctuations.

We now turn to the qTLS DOS. Fig. \ref{fig:linewidths}b shows the expected frequency dependence of the qTLS DOS, given by the bounds of eq. (1) and (2).
This is in good agreement with the number of strongly coupled qTLS
that we observe in resonators with different $f_0$, shown in Fig. \ref{fig:linewidths}c; data obtained across several cooldowns and devices. A similar dependence was obtained by tuning the TLS energy by an applied electrostatic field, see \cite{suppl}. 
In contrast, for conventional non-interacting TLS in amorphous materials we expect
a DOS $P(E)={\rm {const}}$ \cite{phillips}, while for interacting
TLS we expect a weak pseudogap given by $P(E)=E^{\mu}$ \cite{faoro2015, moshe2020}.  $\mu\approx0.2$ is typically found experimentally \cite{burnett2014,degraaf2018}, which would for $f_{0}=3-6$ GHz only result in a $\sim10\%$
variation in the TLS DOS. 
As expected, we find through the temperature dependent shift of $f_0$ which samples all TLS, including the more numerous conventional glassy TLS,  that the intrinsic loss tangent fall within $F\tan\delta_{i}=0.7\pm0.1 \cdot10^{-5}$,
with no observed dependence on $f_0$ \cite{suppl}. 



The data in Fig. \ref{fig:linewidths}b has not been
scaled to the surface area occupied by each resonator, which would yield an even stronger non-linear
dependence ($A\propto f_{0}$). If we make the assumption
that the observed qTLS uniformely occupy the surface of the superconductor we get for $f_0=4.9$ GHz an observed density of qTLS of $0.004$ GHz$^{-1}$ ${\rm \mu}$m$^{-2}$, or $1$ GHz$^{-1}$ per $(\sim3000\xi_0)^2$.
In our devices the distance between electrodes is small, resulting
in significant electric fields away from the immediate vicinity
of metal edges. In fact, the maximum observed TLS coupling strengths of
up to $100-200$ kHz are obtained for electric dipoles of size $p_0\sim 1$
${\rm e\lyxmathsym{\AA}}$ (the expected scale of a qTLS dipole moment \cite{suppl} and commonly found also for conventional TLS \cite{lisenfeld2019, sarabi2016}) across a significant part of the device
surface. Near metal edges we would expect coupling strengths well exceeding
1 MHz \cite{suppl}. This we do not observe, implying instead that the qTLS are detected across
most of the device surface.
The above density is much smaller than the total density of TLS we extract
from $F\tan\delta_{i}$. We find a surface density
of $\rho_{{\rm \Delta,A}}\approx50$ GHz$^{-1}$ ${\rm \mu}$m$^{-2}$
of weakly coupled TLS which are mainly responsible for the background
loss and noise \cite{suppl} previously studied in detail \cite{burnett2014, degraaf2018}. These conventional TLS coexist with the qTLS and constitute the much broader background of loss and fluctuations. However,
since the qTLS are located at the metal surface they interact
weakly with the conventional glassy TLS bath which is mainly situated
some distance away, e.g. in the oxide on the superconductor and adsorbants
on top, resulting in higher qTLS coherence. 



Finally, to determine if there indeed is a link between the qTLS and
superconductivity we measure a strongly coupled qTLS 
in an applied
magnetic field $B$ parallel to the superconducting film, shown in Fig.
\ref{fig:tlsbfield}a. As expected, we observe a quadratic suppression
of $f_0$ owing to the kinetic inductance 
increasing as $\Delta$ is suppressed. At the same time
we track the anti-crossing with the qTLS and extract the minimum qTLS
energy $\epsilon_0(B)$. As shown in Fig. \ref{fig:tlsbfield}b
$\epsilon_0$ is also suppressed in magnetic field on a
scale comparable to the weak (few \%) suppression of $\Delta$. 
We also note that the conventional TLS bath remains unaffected up to moderate fields \cite{degraafAPL2018}.


High-resolution spectral mapping of TLS in qubits has only been performed in the range $ 5-6$ GHz \cite{klimov2018}.  Below $3$ GHz qubits
in other studies were limited by other loss mechanisms or statistics are insufficient \cite{nguyen2019,quintana2017,yan2015}.
Only Ref. \cite{quintanathesis} hints at a lower TLS density in Al qubits below $\simeq 3$ GHz,
suggesting that the mechanism
responsible for the formation of qTLS could indeed be quite common, and motivating the need to explore lower frequencies.
The critical parameter for the DOS, the strength of fluctuations $g$, is to the best of our knowledge not known for thin Al films. However, it may even be enhanced due to thickness fluctuations contributing to $\delta \Delta$ \cite{ivry2014,aldelta} and the upper frequency for qTLS is bounded by a smaller $\Delta_0$.  
We  note that in the dirty limit the details of the theory developed here is not strictly valid. However, the general mechanism would still apply, but the optimal shape of fluctuations is likely to be different.



We now turn to the implications that these qTLS states have for a range of superconducting technologies and for the coherence of qubits. A typical sample of $\sim 1$ cm$^{2}$ size is hit by a high energy
cosmic ray once in every $\sim10$ s, creating a shower of phonons that quickly relax to low frequencies where their relaxation time becomes long. These
phonons produce QPs well above $\Delta$. 
Eventually these QPs fall into the traps discussed here where they remain
essentially forever at low temperatures. Such localization of QPs prevents their recombination and equilibration. Those close to a qubit can result in its excitation or relaxation, and hence the QP localization results in a
long period of qubit performance degradation. As this degradation affects all qubits on a
sample it becomes a catastrophic event for quantum computation, a conclusion supported by recent data
on charge noise in qubits \cite{mcdermott}
 and by the studies of QP bursts in granular Al
resonators \cite{Henrique2019} and kinetic inductance detectors \cite{karatsu2019}. 

Identifying and eliminating strongly coupled TLS is crucially important for scaling-up of superconducting quantum computing.  
All the mechanisms that lead to gap fluctuations are difficult to control in thin films used for resonators and qubits, but our results suggest these must be carefully addressed in order to achieve this goal.  


\section*{Acknowledgements}
We thank R. McDermott and M. Gershenson for fruitful discussions. Samples were fabricated in the nanofabrication facilities of the department of microtechnology and nanoscience at Chalmers university. This work was supported by the UK Department of Business, Energy and Industrial Strategy (BEIS), the EU Horizon 2020 research and innovation programme (grant agreement 766714/HiTIMe), the Swedish Research Council (VR) (grant agreements 2016-04828 and 2019-05480), EU H2020 European Microkelvin Platform
(grant agreement 824109), and Chalmers Area of Advance NANO/2018.
J.J.B acknowledges financial support from the Industrial Strategy Challenge Fund Metrology Fellowship as part of BEIS.

\let\oldaddcontentsline\addcontentslin
\renewcommand{\addcontentsline}[3]{
