# Anomalous Loss Reduction Below Two-Level System Saturation in Aluminum
Superconducting Resonators

###### Abstract

Superconducting resonators are widely used in many applications such as qubit
readout for quantum computing, and kinetic inductance detectors. These
resonators are susceptible to numerous loss and noise mechanisms, especially
the dissipation due to two-level systems (TLS) which become the dominant
source of loss in the few-photon and low temperature regime. In this study,
capacitively-coupled aluminum half-wavelength coplanar waveguide resonators
are investigated. Surprisingly, the loss of the resonators was observed to
decrease with a lowering temperature at low excitation powers and temperatures
below the TLS saturation. This behavior is attributed to the reduction of the
TLS resonant response bandwidth with decreasing temperature and power to below
the detuning between the TLS and the resonant photon frequency in a discrete
ensemble of TLS. When response bandwidths of TLS are smaller than their
detunings from the resonance, the resonant response and thus the loss is
reduced. At higher excitation powers, the loss follows a logarithmic power
dependence, consistent with predictions from the generalized tunneling model
(GTM). A model combining the discrete TLS ensemble with the GTM is proposed
and matches the temperature and power dependence of the measured internal loss
of the resonator with reasonable parameters.

###### Abstract

The Supplemental Material discusses the experimental details and measurement
setup, reviews the standard and generalized tunneling models of two-level
systems, reviews the master equation formalism for the fluctuators in the
generalized tunneling model, and presents an overview of the nonequilibrium
quasiparticle and phonon model used to describe the nonequilibrium
quasiparticle loss in the resonator. The electric field coupled to the TLS
dipoles is estimated from the peak electric field in the resonator obtained in
a finite element solver.

###### keywords:

Two-level system, Superconducting resonators

Tamin Tai Jingnan Cai Steven M. Anlage**

Dr. Tamin Tai, Jingnan Cai, Prof. Steven M. Anlage
Quantum Materials Center, Department of Physics, University of Maryland,
College Park, MD 20742-4111, USA
Email Address<EMAIL_ADDRESS>
##  1 Introduction

Two-dimensional (2D) planar high internal quality factor ($Q_{i}$)
superconducting resonators have been widely fabricated and investigated in
recent times for applications such as single photon detectors,[1] kinetic
inductance detectors,[2] and quantum buses in quantum computing technology[3].
Tremendous progress has been made in terms of design, fabrication and
measurement techniques, which has led to orders of magnitude increase in
coherence time and improved quantum fidelity of the quantum gates [3, 4, 5].
In microwave measurements, although all qubits are operated at an excitation
frequency well below the superconducting gap energy, microwave photons can be
absorbed by quasiparticles, which in turn interact with the phonon bath,
creating non-equilibrium distributions of both quasiparticles and phonons[6,
7, 8, 9]. This process affects the population of quasiparticles, in addition
to pair-breaking process induced by cosmic rays,[10] higher order microwave
harmonics, and stray infrared radiation [11, 12, 13]. These non-equilibrium
quasiparticles are one limiting factor on superconducting resonator $Q_{i}$
and qubit coherence, which can reduce both the qubit relaxation time
$(T_{1}^{Qubit})$ and the coherence time $(T_{2}^{Qubit})$[14].

Another comparable loss mechanism due to two-level systems (TLS) is also
ubiquitous in 2D superconducting resonators [15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28, 29, 30]. Despite the elusive microscopic origin of the
TLS (some recent works suggesting hydrogen impurities in alumina as one
candidate for TLS [31, 32]), TLS can be simply modeled as electric dipoles
that couple to the microwave electric field. In general, TLS are abundant in
amorphous solids and can also exist in the local defects of crystalline
materials. They are found in three kinds of interfaces in the superconducting
resonators: the metal-vacuum interface due to surface oxide or contaminants;
the metal-dielectric substrate interface due to residual resist chemicals and
buried adsorbates; and the dielectric substrate-vacuum interface with
hydroxide dangling bonds, processing residuals, and adsorbates[33]. To address
these issues, different kinds of geometry of coplanar waveguide (CPW)
structure have been proposed and fabricated, with more care given to the
surface treatment to alleviate the TLS losses [34]. For example, a trenched
structure in the CPW helps to mitigate the metal-dielectric TLS interaction
with the resonator fields[35, 36]. These efforts have improved the 2D
resonator intrinsic quality factor to more than 1 million in recent
realizations of high-$Q_{i}$ resonators [35, 37, 38, 39, 36]. Nevertheless,
TLS still exist even in extremely high $Q_{i}$ 3D superconducting radio
frequency cavities used in particle accelerator applications[40]. Recently,
other sources of TLS loss have been proposed based on quasiparticles trapped
near the surface of a superconductor[41].

Clearly TLS loss is a universal issue in superconducting resonators. However,
at microwave frequencies, this loss was long thought to be constant under low
microwave power and low temperature below TLS saturation [16, 42, 43, 18, 17,
28]. Measurements in this regime were limited due to the constraints of noise
levels in both electronic equipment and the thermal environment. Therefore,
experimental investigation of TLS at low temperatures and microwave excitation
are important, and would assist the superconducting quantum information
community to understand its effect on operating quantum devices.

Figure 1: An SEM image of the aluminum CPW resonator on a sapphire substrate.
(b) (c) Zoom-in SEM images of the left and right capacitive couplers. (d) AFM
image highlighting the tapered center conductor with a 1 $\mu m$ wide center
trace near the center of the resonator, and (e) AFM topography image
highlighting the 5 $\mu m$ wide capacitive coupler from (b) or (c). Note that
the AFM probe scanning direction is 45 degrees with respect to the center-line
direction to reduce AFM scanning artefacts. (f) Line scan profile of AFM image
to show thickness of the center line in (d). (g) Line scan profile of the
capacitive coupler. Both line scans show an Al film thickness of 70 nm.

We have designed a 2D half wavelength resonator with a tapering geometry that
gradually shrinks the signal line width $w$ from 50 $\mu m$ down to 1 $\mu m$
at the center where many three-junction flux qubits could be hosted and
strongly coupled for the study of the collective behavior of quantum meta-
materials. Analogous to cavity quantum electrodynamics, qubits serve as
artificial meta-atoms with mutual coupling [44, 45, 46, 47, 48] and can be
read out through the dispersive frequency shift of the cavity[49, 50, 51].
Theoretical publications discussing the physics of qubit arrays coupled to the
harmonic cavities predict a number of novel collective behaviors of these
meta-atoms [52, 53, 54]. In this paper, we report our finding on the TLS loss
in the low power and low temperature limit of this particular design of
capacitively-coupled half-wavelength resonator, without the qubits. The
technique of very low power microwave measurement with low noise to enhance
the signal-to-noise ratio (SNR) is critical for measuring this TLS behavior.

##  2 Experimental Methods

Aluminum (Al) half-wavelength $(\lambda/2)$ CPW resonators on sapphire
substrates were designed with a center line width $w=50$ $\mu m$ and spacing
$s=30$ $\mu m$ (the distance between center conductor line and ground plane as
illustrated in Fig. 1(b) to maintain the characteristic impedance near 50
$\Omega$ in the meander part. At the center of the resonator a tapering
structure narrows the center line width down to $w=1\>\mu m$ and spacing to
$s=12\>\mu m$, which gradually increases the characteristic impedance to 100
$\Omega$ at the resonator center. Fig. 1 (a) shows a perspective view of the
resonator in a diced chip with a designed fundamental frequency around 3.6
GHz. The entire resonator is surrounded by many 10 $\mu m$ by 10 $\mu m$
vortex moats. The resonator is symmetric and capacitively coupled through 5
$\mu m$ gaps (Fig. 1 (b) and 1 (c)) in the center conductor. A topographic
image of the narrowed resonator center section is shown in Fig. 1 (d) with a
critical dimension around $w=$ 1 $\mu m$ in width. Line cuts shown in the AFM
images in Fig. 1 (d), (e) show that the Al film is 70 nm thick.

This CPW resonator was fabricated using standard photo-lithography procedures.
First, a 70 nm thick Al film was deposited on a 3-inch diameter sapphire wafer
using thermal evaporation technology with a background pressure of $\sim
3\times 10^{-7}$ mbar. Then a thin SHIPLEY1813 photo-resist was coated on top
of the film and exposed to UV through the designed photomask. The UV exposed
wafer was developed and then wet etched by commercial Transene Aluminum
Etchant. The remaining photoresist was stripped off by acetone and the entire
wafer was cleaned by methanol and isopropanol. Finally, the wafer was coated
in a protective photo-resist and then diced into many chips. After dicing, the
protective photo-resist was removed and the chip was mounted on a printed
circuit board bolted inside a copper box. Several lumps of indium were pressed
between the on-chip ground planes and the copper box ground to achieve a
continuous ground contact, which mitigates parasitic resonant microwave modes
due to uneven electrical grounding. The indium lumps also secured the chip in
the center of the printed circuit board. The on-chip transmission line is
wire-bonded to the center conductor of the transmission line on the printed
circuit board by gold wires. Finally, the copper box is capped by a copper lid
to eliminate stray light illumination.

The device was placed in a closed Cryoperm cylinder in a BlueFors (BF-XLD 400)
croygen-free dilution refrigerator (base temperature 10 mK) to minimize any
stray DC magnetic field, and the shield was thermally anchored to the mixing
chamber plate. The microwave excitation was attenuated by a series of
attenuators in the input line at different cooling stages in the dilution
fridge before going into the resonator to reduce the noise. The transmitted
signal was amplified twice through a cryogenic amplifier and a room
temperature amplifier before being measured by a Keysight N5242A vector
network analyzer (VNA). The low power measurements were performed using the
smallest intermediate frequency bandwidth (1 Hz) of the VNA, with a 400 kHz
span across the resonance, following 5 averages to reduce the random noise. A
thru calibration of the setup was performed in a separate cool down to
determine the overall loss/gain in the transmission lines leading to the
resonator. Further details of the experimental setup for the high SNR
measurement at very low microwave power can be found in section I of the
Supplemental Material.

##  3 Experimental Data

The measured transmitted signal $(S_{21}(f))$ has a fundamental $(\lambda/2)$
resonance peak around $f=$ 3.644 GHz at the fridge base temperature when
sweeping the frequency, $f$. The complex $S_{21}(f)$ signal is fitted to an
equivalent circuit model of a two-port resonator capacitively coupled to
external microwave excitation [55, 9].

$\displaystyle
S_{21}(f)=|S_{21,in}||S_{21,out}|\left(\frac{Q_{L}/Q_{c}}{1+2jQ_{L}(\frac{f}{f_{0}}-1)}e^{j\phi}\right)+C_{0}$
(1)

where $|S_{21,in}|$ and $|S_{21,out}|$ are the net loss or gain in the
transmission of the input and output line, respectively. $Q_{L}$ is the loaded
quality factor. $Q_{c}$ is the coupling quality factor representing the
dissipation to the external circuit, $j=\sqrt{-1}$, $f_{0}$ is the resonance
frequency of the half-wavelength $(\lambda/2)$ CPW resonator, $\phi$ is the
phase and $C_{0}$ is an offset in the complex $S_{21}$ plane due to background
contributions.[55] The internal quality factor, $Q_{i}$, inversely
proportional to the the internal loss, $\delta=Q_{i}^{-1}$, is extracted from
the identity $1/Q_{L}\equiv 1/Q_{i}+1/Q_{c}$. The absorbed power $P_{ab}$ of
the resonator is characterized by the average number of circulating microwave
photons in the cavity on resonance, which can be estimated using the
approximation [9, 56] $\langle
n\rangle=\frac{2Q_{L}^{2}P_{in}}{Q_{c}\hbar\omega_{0}^{2}}$ for a two-port
device, where $\hbar$ is the reduced Planck constant, and $\omega_{0}=2\pi
f_{0}$ is the angular frequency of the resonance.

Figure 2:  (a) Temperature dependent first harmonic resonant frequency shift
$\Delta f/f_{0}(6\text{mK})$, with $\Delta f=f_{0}-f_{0}(6\text{mK})$ of the
$\lambda/2$ Aluminum co-planar waveguide resonator on sapphire substrate
measured at different excitation powers (average photon numbers). Here
$f_{0}(6\text{mK})$ is the resonance frequency measured at the base
temperature for each excitation power . (b) Temperature dependent loss
(inverse of intrinsic quality factor, $Q_{i}^{-1}$) at its first harmonic
frequency of a Aluminum co-planar waveguide resonator on sapphire substrate
measured at different circulating photon numbers $\langle n\rangle$. Some of
the error bars are smaller than the data point such as those for the high
power and temperature measurements.

Figure 2(a) illustrates the temperature dependence of the fractional resonant
frequency shift from the resonance frequency at lowest temperature,
$(f_{0}(T)-f_{0}(6\text{mK}))/f_{0}(6\text{mK})$, for different circulating
microwave photon numbers inside the CPW resonator, where 6 mK is the measured
fridge base temperature. The resonance frequencies start at their maxima at
the fridge base temperature and then show local minima around 60 mK. This
phenomenon seems to be independent of the average circulating photon number
and can be explained by the standard tunneling model of TLS [43]. Upon further
increasing the temperature above 150 mK, the resonance frequencies quickly
decrease due to the thermal quasiparticles, which increases the real and
imaginary parts of the surface impedance of the superconducting resonator. The
inset focuses on the low temperature regime and shows a very small power
dependence that is qualitatively similar to the strong field correction to the
frequency shift in STM proposed by Gao, which predicts smaller frequency
shifts for higher power [57].

The temperature dependence of the measured internal loss is shown in Fig.
2(b). For high power measurements ($\langle n\rangle>10^{6}$), the loss is
constant at low temperatures (below 150 mK) which is expected for the typical
non-interacting TLS. At higher temperatures, the loss increases due to thermal
quasiparticles. For low power measurements ($\langle n\rangle<10^{6}$),
starting from the minimum temperature, the loss has an unusual increase at low
temperatures, from the base temperature to a peak at 40 mK. The loss then
drops with increasing temperature following the equilibrium value of the
population difference in TLS [20, 58]. Similar to the high power measurements,
the loss rises again above 150 mK due to thermal quasiparticles. The observed
loss decrease with decreasing temperature from 40 mK to 10 mK has not been
explicitly acknowledged and discussed in prior work of microwave
superconducting resonators until recently [59]. Indications of an upturn in
$Q_{i}(T)$ has been attributed to poor SNR and therefore treated as not
statistically significant[58, 60].

##  4 Modeling

###  4.1 Frequency shifts

The power and temperature dependent frequency shifts have primarily been
attributed to TLS and the dynamics of quasiparticles. These two mechanisms
could overlap and become difficult to distinguish in the operation of many
superconducting devices, including resonators and qubits, [14]. A simple model
that combines both quasiparticles and TLS contribution in one equation
describes the resonance frequency $\Delta f$ data in Fig. 2,[25, 34, 61]

$\displaystyle\frac{f_{0}(T)-f_{0}(0)}{f_{0}(0)}=$
$\displaystyle\frac{\delta_{0}}{\pi}\left(Re\left[\Psi(\frac{1}{2}+\frac{\hbar\omega}{2\pi
jk_{B}T})\right]-log(\frac{\hbar\omega}{2\pi k_{B}T})\right)$
$\displaystyle-\frac{\alpha}{2}\left(\frac{n_{qp}}{2N_{0}\Delta_{S0}}\left[1+\sqrt{\frac{2\Delta_{S0}}{\pi
k_{B}T}}\exp(\zeta)I_{0}(\zeta)\right]\right)$ (2)

where $\zeta=\frac{hf_{0}}{2k_{B}T}$, $f_{0}$ is the resonance frequency as a
function of the temperature, $\delta_{0}$ is the zero temperature and zero
power loss tangent from the TLS, $\Psi(\cdot)$ is the digamma function,
$\alpha=L_{kinetic}/L_{total}$ is the kinetic inductance fraction of the CPW
resonator, $N_{0}$ is the single spin density of states, $\Delta_{S0}$ is the
aluminum superconducting gap at zero temperature, and $I_{0}(\cdot)$ is the
0th order modified Bessel function of the first kind. The first term in Eq. 2
represents the frequency shift caused by the TLS mechanism [20, 23] and the
second term is the frequency shift due to quasiparticles using the Bardeen-
Cooper-Schrieffer (BCS) model for $k_{B}T,\ hf_{0}\ll\Delta_{S0}$, and written
explicitly in terms of quasiparticle number density $n_{qp}$ [61], including
both thermal and non-equilibrium quasiparticles. However, the model with only
thermal quasiparticle $n_{th}=2N_{0}\sqrt{2\pi
k_{B}T\Delta_{S0}}\exp\left(-\frac{\Delta_{S0}}{k_{B}T}\right)$ (valid for
$T\ll T_{c}$) seems to match the measurement sufficiently well, where
$N_{0}=10^{47}$ $J^{-1}\>m^{-3}\approx 1.74\times 10^{4}\mu eV^{-1}\mu m^{-3}$
is the single spin density of states at the Fermi level [7, 13].

Figure 3:  Temperature dependent fundamental ($\lambda/2$) mode resonant
frequency $f_{0}(T)$ of the Al CPW resonator on sapphire substrates at an
external microwave excitation creating around one circulating photon. The
inset highlights the low temperature regime where the frequency shift is
dominated by the TLS mechanism. The dots are experimental data and solid line
is the model fit to Eq. (2).

The fit to the frequency shift data is shown in Fig. 3, and the extracted
fitting parameters indicate that the aluminum superconducting gap at zero
temperature is $\Delta_{S0}\sim$170 $\mu eV$, a value close to the BCS gap
approximation which is $1.76k_{B}T_{c}$ with transition temperature
$T_{c}=1.12\>K$. The values of the other fitting parameters are $\alpha\approx
0.014$, and $\delta_{0}=9.6\times 10^{-6}$. The values of $\alpha$ and
$\delta_{0}$ are consistent with other results on a variety of similar
superconducting resonators[23, 27, 62, 34].

###  4.2 Internal loss

Since the temperature dependent internal loss is dominated by the well-known
thermal quasiparticles above 150 mK, this analysis focuses only on the low
temperature data. The power dependence of the loss $Q_{i}^{-1}(T)$ is shown in
Fig. 4 at different temperatures below the onset of thermal quasiparticle
effects. Clearly, the loss shows a gradual power dependence above the low-
power saturation, similar to previous experimental observations[63, 27, 64,
37], and is not consistent with STM shown as the dashed curves.

To account for the slower power dependence, many improvements on the STM have
been proposed, such as introducing more than one species of TLS in the
dielectrics, [65, 66, 67, 68, 69] and accounting for the nonuniform field
distribution in the resonator [24]. In addition, there is another approach
that generalizes the STM to include the strong interactions among TLS [70, 71,
72], resulting in the generalized tunneling model (GTM) that can produce the
logarithmic power dependence shown as the black dotted line in Fig. 4.
However, none of the existing models predicts a strong temperature dependence
of loss below the TLS saturation. To interpret this unusual loss reduction in
our aluminum resonators at low power and low temperature, we go beyond the
assumption of a uniform distribution of TLS and invoke the discrete TLS
contribution to the loss at low temperatures. A simple modification that sums
over the discrete and detuned TLS responses near the resonance as in Eq. (15)
is proposed. When combined with GTM, this model reproduces the full power and
temperature dependence of the loss data: the gradual power dependence at high
power as well as the observed anomalous temperature dependence of loss for
$\langle n\rangle<10^{2}$ and $T<50$ mK. To lay the foundations of the
proposed model, several key concepts of STM and GTM are first reviewed.

Figure 4:  Internal loss $Q_{i}^{-1}$ as a function of power (measured by
photon number $\langle n\rangle$) at different temperatures for an Aluminum
resonator on a sapphire substrate. The scatter plots are experimental data
points, and the dashed lines are the fitting curves from the STM given in Eq.
(9). The large deviation from STM power dependence above TLS saturation power
is typical for strongly interacting TLS, which are better treated with GTM.
The black dotted line is the logarithmic power dependence at high excitation
power predicted by GTM.

####  4.2.1 Conventional model for TLS loss

The TLS formalism is based on a simple model for a single TLS that can be
described by the Hamiltonian,
$H_{TLS}=\frac{1}{2}\begin{pmatrix}-\Delta&\Delta_{0}\\\
\Delta_{0}&\Delta\end{pmatrix}$ where $\Delta$ is the asymmetry of the double
well potential and $\Delta_{0}$ is the tunneling barrier energy between the
potential wells[20]. A typical resonator hosts an ensemble of TLS with
different values of $\Delta$ and $\Delta_{0}$ with their distribution function
given as

$\displaystyle
P(\Delta,\Delta_{0})=(1+\mu)\left(\frac{\Delta}{\Delta_{max}}\right)^{\mu}\frac{P_{0}}{\Delta_{0}}\quad\begin{cases}\mu=0&\text{no
TLS interactions \cite[cite]{[\@@bibref{}{Phillips87}{}{}]}}\\\ \mu\approx
0.3&\text{strong TLS interactions
\cite[cite]{[\@@bibref{}{cuevas_density_1989,faoro_interacting_2015}{}{}]}}\end{cases}$
(3)

where $P_{0}\approx 10^{44}J^{-1}m^{-3}$ is the density of states for TLS,
$\Delta_{max}\sim\varepsilon_{max}\sim 100K$ [72], and the exponent $\mu$
approximately characterizes the logarithmic correction to the density of
states of TLS due to strong TLS interactions.

The dynamics of a single TLS can be described by the linearized Bloch
equations of the pseudospin $\vec{S}(t)$ (see Supplementary Material Section
II for details) that is characterized by the following four rates: Rabi
frequency $\Omega\propto|\vec{E}|$, the frequency of the driving field
$\omega$, the splitting between the two eigenenergies of TLS
$\varepsilon=\sqrt{\Delta^{2}+\Delta_{0}^{2}}$, and the longitudinal and
transverse relaxation rates of the TLS $\Gamma_{1,2}$ which are defined as

$\displaystyle\Gamma_{1}=\left(\frac{\Delta_{0}}{\varepsilon}\right)^{2}\left[\frac{\gamma_{L}^{2}}{v_{L}^{5}}+\frac{2\gamma_{T}^{2}}{v_{T}^{5}}\right]\frac{\varepsilon^{3}}{2\pi\rho\hbar^{4}}$
$\displaystyle\coth(\frac{\varepsilon}{2k_{B}T})=\left(\frac{\Delta_{0}}{\varepsilon}\right)^{2}\Gamma_{1}^{max}\quad\text{\cite[cite]{[\@@bibref{}{Jack72,Hunk76,Bla77,Gao_phdthesis}{}{}]}}$
(4) $\displaystyle\Gamma_{2}=\Gamma_{2}^{ph}+\Gamma_{ds}$
$\displaystyle\text{where}\quad\Gamma_{2}^{ph}=\Gamma_{1}/2\quad\text{and}\quad$
$\displaystyle\Gamma_{ds}\sim
10^{-3}(k_{B}T/\varepsilon_{max})^{\mu}k_{B}T/\hbar\quad\text{\cite[cite]{[\@@bibref{}{jankowiak_origin_1993,faoro_interacting_2015}{}{}]}}$
(5)

Eq.(4) describes the longitudinal relaxation rate dominated by the phonon
process where $\gamma_{L}$ and $\gamma_{T}$ are the longitudinal and
transverse deformation potentials, respectively, $v_{L}$ and $v_{T}$ are the
longitudinal and transverse sound velocities, $\rho$ is the mass density, and
$\Gamma_{1}^{max}$ is the maximum $\Gamma_{1}$ for the TLS with energy
splitting $\varepsilon$, when $\Delta_{0}=\varepsilon$. Eq.(5) defines the
transverse relaxation rate where $\Gamma_{ds}$ is the dephasing rate of the
resonant TLS energy level $\varepsilon$, caused by its interactions with
thermally activated TLS whose $\varepsilon\lesssim k_{B}T$, and is valid for
low temperature measurement ($T<1$ K) [74]. We note that $\Gamma_{ds}\sim
10^{6}\ \text{Hz}$ dominates over $\Gamma_{1}\sim 10^{3}\ \text{Hz}$ in the
typical cryogenic measurement of amorphous dielectrics [72]. Therefore,
$\Gamma_{2}$ is often approximated as $\Gamma_{ds}$ and scales as
$T^{1+\mu}=T^{1.3}$.

In STM, the resonant dielectric response of a single TLS is expressed as [15,
17, 71]:

$\chi_{res}=\frac{m(\omega-\varepsilon/\hbar-i\Gamma_{2})}{(\omega-\varepsilon/\hbar)^{2}+\Gamma_{2}^{2}(1+\Omega^{2}\Gamma_{1}^{-1}\Gamma_{2}^{-1})}$
(6)

where $m=\tanh(\varepsilon/(2k_{B}T))/2$ is the equilibrium value of $\langle
S_{z}^{0}\rangle$. The single TLS loss corresponds to the imaginary part of
the response function in Eq.(6) which is in the form of a Lorentzian in
$\varepsilon/\hbar$ centered at $\omega$ with a width:

$\displaystyle w=\Gamma_{2}\sqrt{1+\kappa}$ (7)

where $\kappa=\Omega^{2}\Gamma_{1}^{-1}\Gamma_{2}^{-1}$. For a typical TLS
with $\varepsilon/h\approx$ 5 GHz and at reasonably low temperatures and
powers, the width of its response $w\approx\Gamma_{2}\sim 1\
\text{MHz}\ll\omega$. Due to this sharp Lorentzian response function, the
total loss is dominated by the resonant TLS whose energies
$\varepsilon\sim\hbar\omega$.

The total dielectric loss is simply the integral of the single TLS
contribution Eq.(6) over the distribution of the TLS Eq.(3) [17, 20, 57].

$\delta_{TLS}=\frac{1}{\epsilon_{r}\epsilon_{0}}\int\int\int
P(\varepsilon,\Delta_{0})\left(\frac{\Delta_{0}d_{0}}{\varepsilon}\right)^{2}\frac{\cos^{2}{\theta}}{\hbar}\frac{m\Gamma_{2}}{\Gamma_{2}^{2}(1+\kappa)+(\varepsilon/\hbar-\omega)^{2}}d\varepsilon
d\Delta_{0}d\theta$ (8)

where $\epsilon_{r}\epsilon_{0}$ is the permittivity of the host dielectric
material, $P(\varepsilon,\Delta_{0})$ is the distribution function of coherent
TLS, obtained from Eq. (3) with a change of variable from $\Delta$ to
$\varepsilon$, $d_{0}$ is the maximum transition electric dipole moment of the
TLS with energy splitting $\varepsilon$, $\vec{E}$ is the applied microwave
electric field on the TLS dipole, and $\theta$ is the angle between the
applied electric field and the TLS dipole moment.

Evaluating this integral leads to the famous STM prediction of TLS loss [43]

$\delta_{TLS}=\frac{\pi
P_{0}d_{0}^{2}}{3\epsilon_{r}\epsilon_{0}}\frac{\tanh{\frac{\hbar\omega}{2k_{B}T}}}{\sqrt{1+(\Omega/\Omega^{c})^{2}}}$
(9)

where $\Omega^{c}\propto\sqrt{\Gamma_{1}^{max}\Gamma_{2}}$ is the critical
Rabi frequency that characterizes the saturation of TLS. The loss is expected
to have an inverse square root dependence on power after the TLS saturation,
$\delta_{TLS}\sim\Omega^{-1}\propto P_{ab}^{-0.5}$ for $\Omega\gg\Omega^{c}$,
which is much faster than observed in the data in Fig. 4. In fact, if one fits
the data with a general power law [27] where the square root in the
denominator of Eq.(9) is replaced by a fitting parameter, the resulting
exponent is around -0.15, indeed a slower power dependence than predicted in
STM.

####  4.2.2 Effect of fluctuators on TLS loss

The dephasing rate introduced in Eq.(5) describes the resonant TLS energy
drift resulting from an average of weak interactions among TLS [20, 72], which
cannot incorporate stochastic strong interactions, such as those from
fluctuators [75, 76, 71, 70, 72, 77]. Fluctuators can be modeled as incoherent
TLS whose $\Gamma_{2}^{ph}\geq\varepsilon$ as opposed to the coherent TLS
introduced above in STM[72]. If strongly coupled with the coherent resonant
TLS, the fluctuators can move the latter in and out of resonance with a jump
rate $\gamma$ and effectively create a random telegraphic noise on the energy
level $\varepsilon\rightarrow\varepsilon+\xi(t)$. This random noise cannot be
treated as a contribution to the spectral diffusion as in Eq. (5), since they
are not Gaussian-distributed. The exact solution to the Bloch equation will
depend on the relationship between
$\gamma,\Omega,\Gamma_{1}^{max},\text{and}\Gamma_{2}$. $\Gamma_{1}^{max}$ is
abbreviated to $\Gamma_{1}$ for clarity in the following discussion, which
mainly focuses on the interaction between fluctuators and one resonant TLS.
Thus, the distribution of values of $\Gamma_{1}$ for an ensemble of TLS is not
invoked until the last step of integration to calculate the loss, and is not
relevant to the fluctuators-induced effect.

When the jump rate $\gamma$ is slow compared to the dynamics of the resonant
TLS characterized by the rates $\Omega,\Gamma_{1},\Gamma_{2}$, the stationary
solution similar in form to Eq.(6) can still be used with the substitution
$\varepsilon\rightarrow\varepsilon+\xi$. After averaging over the distribution
of the fluctuator jumps $\xi$, the response of a single TLS coupled to
low-$\gamma$ fluctuators is obtained (See Supplementary Material Section III
A):

$\displaystyle\chi_{res}=m\frac{\omega-\varepsilon/\hbar-i(\Gamma_{2}+\Gamma_{f}/\sqrt{1+\kappa})}{(\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f})^{2}+(\omega-\varepsilon/\hbar)^{2}}$
(10)

which has the same form as Eq.(6) but with the width of the Lorentzian widened
by $\Gamma_{f}\propto(k_{B}T/\hbar)\\\
(k_{B}T/\varepsilon_{max})^{\mu}\lesssim\Gamma_{2}$ due to the weakly-coupled
low-$\gamma$ fluctuators. For a continuous distribution of TLS such as
$P(\varepsilon,\Delta_{0})$, the total internal loss is calculated by
integrating Eq.(10) over the distribution function.

$\displaystyle\delta_{TLS}=\frac{1}{\epsilon_{r}\epsilon_{0}}\int\int\int
P(\varepsilon,\Delta_{0})\left(\frac{\Delta_{0}d_{0}}{\varepsilon}\right)^{2}\frac{\cos^{2}{\theta}}{\hbar}\frac{m}{\sqrt{1+\kappa}}\frac{(\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f})}{(\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f})^{2}+(\varepsilon/\hbar-\omega)^{2}}d\varepsilon
d\Delta_{0}d\theta$ (11)

Clearly, the last fraction in the integral is a Lorentzian which evaluates to
a constant after integration over $\varepsilon$, resulting in the same
prediction for internal loss as the STM [72].

On the other hand, when the dynamics of the resonant TLS is dominated by a
fast jump rate, $\gamma\gtrsim\Omega,\Gamma_{1},\Gamma_{2}$, a probabilistic
description of the resonant TLS must be adopted. Instead of directly solving
the linearized Bloch equations Eq.(S2,S3), the master equation or the
evolution equation of the probability distribution $\rho(\vec{S})$ of the
Bloch vector $\vec{S}=(\langle S_{x}^{1}\rangle,\langle
S_{y}^{1}\rangle,\langle S_{z}^{0}\rangle)$ is introduced: [71],

$\displaystyle\frac{\partial\rho}{\partial
t}+\frac{d}{d\vec{S}}(\frac{d\vec{S}}{dt}\rho)=\gamma[\delta(S_{z}^{0}-m)\delta(S_{x}^{1})\delta(S_{y}^{1})-\rho]$
(12)

where $d\vec{S}/dt$ is given in Eq.(S2, S3) with a time independent
$\varepsilon$ where the random jumps $\xi(t)$ are dropped since the fast jumps
are averaged out over a long time. (See Supplementary Material Section III B).
The TLS loss is then extracted by solving for the average $y$ component of
$\vec{S}$, $\overline{\langle S_{y}^{1}\rangle}=\int\rho\langle
S_{y}^{1}\rangle d\vec{S}$, and integrating the solution over a probability
distribution of fluctuator jump rates $P(\gamma)\propto 1/\gamma$ [71, 72],
and a distribution of resonant TLS energies given in Eq.(3).The loss has a
logarithmic dependence on power:

$\displaystyle\delta=m\delta_{0}\text{arcsinh}\left(\frac{\gamma}{\Omega}\right)\Biggr{|}_{\text{max}(\Omega,\Gamma_{2})}^{\gamma_{\text{max}}}\xrightarrow{\gamma_{max}\gg\Omega\gg\Gamma_{1},\Gamma_{2}}m\delta_{0}\ln(\frac{\gamma_{max}}{\Omega})$
(13)

where $\gamma_{max}$ is the maximal jump rate from the fluctuators. This
expression explains the high power limit of the data in Fig. 4 where the
losses from different temperatures converge to a linear trend in the linear-
log plot. This high $\gamma$ fluctuator loss will saturate to a constant value
$\sim m\ln(\gamma_{max}/\Gamma_{2})$ once $\Omega\lesssim\Gamma_{2}$ (See
Supplementary Material Section III B). Thus, it will not affect the low power
behavior of the TLS loss.

More complicated is the case of intermediate jump rates where
$\gamma\sim\Omega,\Gamma_{1},\Gamma_{2}$. A similar master equation as in
Eq.(12) needs to be solved with $\varepsilon\rightarrow\varepsilon+\xi_{k}$
for each different fluctuator state $k$. The jumps are no longer ignored since
their rates are close to the other dynamics ($\Omega,\Gamma_{2},\Gamma_{1}$)
in the system. After obtaining the average solution $\overline{\langle
S_{y}^{1}\rangle}$ for the TLS with energy levels $\varepsilon_{k}$, the same
recipe for the loss calculation can be applied, namely, integrating the
average solution over $P(\gamma)$ and then integrating over the distribution
of the coherent TLS $P(\varepsilon,\Delta_{0})$ . The loss is then (See
Supplementary Material Section III C)

$\displaystyle\delta=\delta_{0}\int\frac{m\Gamma_{2}}{\Gamma_{2}^{2}+(\varepsilon-\omega)^{2}}\left(\frac{1}{1+\kappa(\varepsilon)}\ln{\frac{\gamma_{h}}{\gamma_{l}}}+\frac{\kappa(\varepsilon)(1-n)}{(1+\kappa(\varepsilon))(1+n\kappa(\varepsilon))}\ln{\frac{1+\kappa(\varepsilon)+\gamma_{h}/\Gamma_{1}(1+n\kappa(\varepsilon))}{1+\kappa(\varepsilon)+\gamma_{l}/\Gamma_{1}(1+n\kappa(\varepsilon))}}\right)d\varepsilon$
(14)

where
$\kappa(\varepsilon)=\Xi/\Gamma_{1}=\Omega^{2}\Gamma_{2}/[(\varepsilon/\hbar-\omega)^{2}+\Gamma_{2}^{2}]/\Gamma_{1}$,
and $n$ is the probability that a given TLS is resonant, and is typically
small for a system of many ($\sim 10$) fluctuators (see supplementary material
for [71]), and thus ignored in the final model. $\gamma_{h,l}$ are the upper
and lower bounds of the jump rates and are defined such that
$\gamma_{h}\gtrsim\Xi+\Gamma_{1},\sqrt{\Gamma_{2}^{2}+(\varepsilon/\hbar-\omega)}\gtrsim\gamma_{l}$.
These limits translate to a range for power $\Omega$ where the power
dependence of the loss is dominated by this model:
$\Gamma_{2}\gtrsim\Omega\gtrsim\sqrt{\Gamma_{1}\Gamma_{2}}$. Within this
range, the loss from intermediate $\gamma$ fluctuators is approximately
$\delta_{0}\ln[(\Omega^{2}+\Gamma_{2}^{2})/(2\Omega^{2})]$, a faster
logarithmic power dependence than Eq. (13). At higher powers, the loss becomes
constant $\sim m\ln(2)$. At lower power, the loss saturates to another
constant $m\ln(\Gamma_{2}/\Gamma_{1})$.
In summary, the three different fluctuation rates correspond to three
different power ranges for the power dependence of the loss. In the high power
limit $\Omega\gtrsim\Gamma_{2}$, the effect of fluctuators that induce large
$\gamma$ dominates and leads to a logarithmic power dependence; in the
intermediate power regime
$\Gamma_{2}\gtrsim\Omega\gtrsim\sqrt{\Gamma_{1}\Gamma_{2}}$, the fluctuators
with intermediate $\gamma$ give rise to a faster logarithmic power dependence,
but meanwhile the saturation of TLS just as in STM has a comparable or even
stronger power dependence and overlap in the same power regime; and finally in
the low power limit $\Omega<\sqrt{\Gamma_{1}\Gamma_{2}}$, the typical TLS
saturation in STM is recovered as the contributions from all three different
types of fluctuators become constant in power. The above description
qualitatively matches our experimental observation in Fig. 4.

####  4.2.3 Fit to the internal loss measurement

Although the power dependence of our data as in Fig. 4 agrees with the effect
of fluctuators in the GTM, the original model does not reproduce the observed
temperature dependence. The GTM predicts the same temperature dependence of
the TLS loss in the low power limit as in STM [72] shown as the orange dashed
curve in Fig. 5 (b), which clearly deviates from the extracted low power loss
of TLS. To reconcile this difference, we propose a simple modification to the
TLS model to account for the discrete coherent TLS near the resonance.
Consider the discrete form of the integral in the STM loss expression, Eq.
(11),

$\delta_{TLS}=\frac{P_{0}d_{0}^{2}\Delta\varepsilon}{3\hbar\epsilon_{r}\epsilon_{0}}\sum_{n}\tanh(\frac{\varepsilon_{n}}{2k_{B}T})\frac{\Gamma_{2}+\Gamma_{f}/\sqrt{1+\kappa}}{(\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f})^{2}+(\varepsilon_{n}/\hbar-\omega)^{2}}$
(15)

where the index $n$ denotes the coherent TLS near the resonance and
$\Delta\varepsilon$ is the average energy spacing in the TLS spectrum. We
believe that Eq.(15) is justified since the number of coherent TLS inside the
resonator bandwidth is $\sim 1$ for a TLS volume around $100\mu m^{3}$ [78],
and many previous works have observed individual TLS in microwave resonators
[79, 80, 81, 82]. For the TLS exactly on resonance, $\varepsilon=\hbar\omega$,
its loss $\delta_{TLS}\sim\Gamma_{2}^{-1}\propto T^{-1.3}$ at low power, and
is the classic result for the single TLS model in STM[83]. However, this
stands in clear contrast to the observed reduction in loss at low temperature
in Fig. 2 (b) and Fig. 5.

It is thus required that none of the TLS is exactly on resonance
($\nu=\varepsilon_{0}/\hbar-\omega\neq 0$ where $\varepsilon_{0}$ stands for
the energy level of the coherent TLS closest to resonance), a reasonable
assumption given the sparse TLS distribution in the frequency spectrum for a
small volume of TLS-inhabiting dielectrics. Mathematically, the width of the
Lorentzian in the summation $w=\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f}$ dictates
the transition from the low temperature reduced loss to the high temperature
equilibrium result. For small $w$, a discrete sum will deviate from the
integral since the Lorentzian is under-sampled. While for a Lorentzian with
large $w$, a discrete sum with the same sampling rate will approximate the
integral better. Specifically, at low powers ($\kappa\ll 1$),
$w=\Gamma_{2}+\Gamma_{f}$ increases with the temperature and
$w=\Gamma_{2}+\Gamma_{f}\sim\nu$ marks the transition temperature between the
two regimes. For low temperatures ($w\ll\nu$), the Lorentzian term becomes
roughly proportional to $w=\Gamma_{2}+\Gamma_{f}$ which gives the almost
linear temperature dependence of loss. For higher power, $w$ increases with
$\kappa$, which pushes the transition temperature lower and suppresses the low
temperature reduction in loss. And eventually at high powers ($\kappa\gg 1$)
such that $w>\nu$ for all temperatures, the equilibrium temperature dependence
$m=(1/2)\tanh{(\hbar\omega/(2k_{B}T))}$ in STM is recovered in the entire
temperature range. The same discrete summation can be applied to Eq. (14) for
intermediate $\gamma$ fluctuators. On the other hand, Eq. (13) for high
$\gamma$ fluctuators is only modified with the substitution
$\Gamma_{2}\rightarrow\sqrt{\Gamma_{2}^{2}+\nu^{2}}$ due to the sparse TLS
assumption (See Supplementary Material Section III B). The final model that
combines all three contributions is able to reproduce the full temperature
($T=8-110$ mK) and power ($\langle n\rangle=10^{-1}-10^{8}$) dependence of the
loss shown as the solid curves in Fig. 5 (a).

Figure 5:  a) The least squares fit of the GTM with discrete summation to the
full power and temperature dependence of the measured internal loss over the
power and temperature range displayed. b) Plot of $\delta_{TLS}^{0}(T)$
extracted from the average of the low power loss below TLS saturation in
Fig.4. The orange dashed curve is the temperature dependence of STM loss below
saturation power $\propto\tanh{(\varepsilon/(2k_{B}T))}$. The purple dash-
dotted (light blue densely dotted) curve is from the discrete summation of
individual TLS contributions for low (intermediate)-$\gamma$ fluctuators at
zero applied power. The green dotted curve is the temperature dependent low
power limit of the TLS loss induced by high $\gamma$ fluctuators. The blue
solid curve is the sum of contributions from the low, intermediate, and high
$\gamma$ fluctuators. c) Comparison of the temperature dependent rates
determined from the least squares fit.

The fit shows reasonable agreement with the data, with $\text{RMSE}=0.014$.
The fit has a total of 11 fitting parameters for the entire temperature and
power range, which has less degrees of freedom than fitting the data from
different temperatures individually. The different contributions to the loss
below TLS saturation power is plotted in Fig.5 (b) which illustrates that the
discrete TLS coupled to low and intermediate $\gamma$ fluctuates are
responsible for the loss reduction. Different rates in the model determined
from the fit are summarized in Fig. 5 (c). The numerical values for
$\Gamma_{1,2}$ and $\gamma_{max,min}$ are typical for TLS in amorphous
materials [72]. The rates also satisfy the following assumptions in the model:
$\Gamma_{2}\gtrsim\Gamma_{f}$, $\gamma_{max}\gg\Gamma_{2}$. In addition, the
low temperature loss reduction occurs around $40$ mK as expected, when
$\Gamma_{2}+\Gamma_{f}<\nu$, the width of the response is smaller than the
detuning between TLS and the resonance. The other quantities extracted from
the fit are listed below: the volume of TLS-inhabiting dielectrics $18.5\ \mu
m^{3}$, the intrinsic TLS loss $\delta^{TLS}_{0}=5.5\times 10^{-6}$, the other
loss $\delta_{other}=1.29\times 10^{-5}$, and the minimum fluctuator rate
$\gamma_{min}=2.9\times 10^{-4}\text{ Hz}$, and the exponent for the energy
dependent density of states $\mu=0.2$.

##  5 Discussion

The discrete and detuned TLS formalism will not affect the high
$\gamma$-fluctuator contribution to internal loss, since the width of
Lorentzian in the calculation of loss of high $\gamma$ fluctuators $w$ is
widened by $\gamma$ such that $w\sim\gamma>\Delta\varepsilon/\hbar$ (See
Supplementary Material Section III B), which is indicated by the almost flat
region in the green dotted curve at low temperature in Fig.5 (b). However, the
loss from intermediate $\gamma$ fluctuators could be subject to the low
coherent TLS density but to a lesser degree than that from the low $\gamma$
fluctuators, since although the bandwidth of their response $\sim\Gamma_{2}$
is the same (See Supplementary Material Section III C), there are many
intermediate-$\gamma$-fluctuator-induced sublevels for one TLS in one Rabi
cycle which effectively increases the density of available TLS energy levels.
In order to avoid over fitting, this effect was not included in the model
where the same density of states for TLS are assumed for those coupled to
intermediate $\gamma$ fluctuators and the low $\gamma$ fluctuators. Thus, the
same $\Delta\varepsilon$ value is shared for the two different contributions.
This simplification could lead to an underestimation of the loss in the
intermediate power region, as illustrated by the deviation between the fit and
data from $\langle n\rangle=10^{2}\text{ to }10^{6}$.

The treatment above is largely classical where the TLS are treated as dipoles
under classical field. A quantum mechanical approach that studies the Jaynes-
Cummings model of a single TLS strongly coupled to a photon predicts a linear
temperature dependence of the loss similar to our observation [83]. However,
it should be noted that the photon frequency in our measurement 3.64 GHz
corresponds to weak photon-TLS coupling, since the Rabi frequency from the
effective field of a single photon is much weaker than the relaxation rates
$\Gamma_{1,2}$. Additionally, the loss from strongly coupled TLS is predicted
to show saturation in power at $\langle n\rangle\sim 1$, clearly lower than
the observed saturation in the data at $\langle n\rangle\sim 10$ which
corresponds to the weak coupling regime and reproduces the classical result
[83].

Although the fluctuators significantly affect the TLS internal loss, they
should have limited effect on the frequency shifts [72]. The proposed discrete
and detuned TLS formalism would not modify the STM frequency shift prediction
either, because unlike the Lorentzian response function that governs the
internal loss, the response function for frequency shift does not have a
resonant shape and is not sensitive to the reduced sampling from the discrete
TLS.

Ever since the importance of TLS interactions in amorphous solids was
recognized by Yu and Leggett [84], there have been numerous experimental works
demonstrating evidence of TLS interactions,[75, 70, 85] and theoretical works
treating the interacting TLS beyond STM [86, 87, 88], with a recent example by
Burin and Maksymov where they used a similar Master equation formalism [28].
However, the spectral diffusion is assumed to follow a Gaussian process,
unlike the fluctuators introduced by Faoro and Ioffe [71, 72], and the loss is
predicted to have a power dependence faster than STM by a logarithmic factor,
contrary to our observation.

At higher temperatures (above 150 mK), the quasiparticle effects become
important, which corresponds to the upturn in loss in Fig. 2. The
quasiparticles loss is related to its density $n_{qp}$ as,

$\displaystyle\frac{1}{Q_{qp}}$
$\displaystyle=\frac{2\alpha}{\pi}\frac{\sinh(\zeta)K_{0}(\zeta)n_{qp}}{N_{0}\sqrt{2\pi
k_{B}T\Delta_{S0}}}$ $\displaystyle\text{where}\quad n_{qp}$
$\displaystyle=n_{th}+n_{noneq}=2N_{0}\sqrt{2\pi
k_{B}T\Delta_{S0}}\exp\left(-\frac{\Delta_{S0}}{k_{B}T}\right)+n_{noneq}$ (16)

where $n_{noneq}$ is the non-equilibrium quasiparticle density. Similar to the
fit for frequency shift, the model with only $n_{th}$ matches our data with
the same set of fitting parameters $\Delta_{0}=170\mu eV$ and $\alpha=0.014$.
A calculation of the increased quasiparticle density including both thermal
and non-equilibrium quasiparticles at high photon numbers in the half
wavelength resonator based on Mattis-Bardeen equations [89, 90] can be found
in Supplementary Material Section V. However, the results lack any strong
temperature or power dependence below 100 mK. Note that this calculation
includes the dynamics of the non-equilibrium quasiparticle finite lifetime due
to recombination and trapping, with and without photon illumination [91, 92,
93, 94]. (Refer to Fig. S5 in the Supplementary Material).

##  6 Conclusion

We have designed and fabricated capacitively-coupled half wavelength
superconducting aluminum microwave resonators with minimum critical dimension
of $1\ \mu m$ in the center conducting line of the CPW. The temperature and
power dependence of the resonator $Q_{i}$ deviate from the classical standard
tunneling model results. At high applied powers, the internal loss shows
logarithmic power dependence, a signature of the generalized tunneling model
with fluctuators. At powers below TLS saturation, the internal loss decreases
from 50 mK down to the fridge base temperature. We attribute this behavior to
the detuning between TLS and the resonance frequency in a discrete TLS
ensemble. Upon cooling, the single TLS response bandwidth, proportional to
$\Gamma_{2}\propto T^{1.3}$, decreases. When the bandwidth drops below the
detuning between TLS and the resonance frequency defined by the CPW resonator,
the resonant TLS response decreases and contributes less to the internal loss.
The generalized tunneling model is revisited and modified with the discrete
TLS formalism resulting in a comprehensive fit to the measured loss in the
entire low temperature and low power range, with a reasonable set of
parameters.

Acknowledgements

This work is funded by the US Department of Energy through Grant # DESC0018788
(sample preparation and measurements), and the National Science Foundation
through Grant # NSF DMR-2004386. We acknowledge use of facilities at the
Maryland Quantum Materials Center and the Maryland NanoCenter. We thank Rangga
Budoyo and Braden Larsen for assistance with the non-equilibrium
superconductor simulation. We thank Dr. Ben Palmer and Yizhou Huang at
Laboratory for Physical Sciences for assistance with the thermal evaporator
and helpful suggestions.

## References

  * [1] J. Zmuidzinas,  _Annual Review of Condensed Matter Physics_ 2012, _3_ , 1 169.
  * [2] J. J. A. Baselmans, S. J. C. Yates,  _AIP Conference Proceedings_ 2009, _1185_ , 1 160.
  * [3] P. Krantz, M. Kjaergaard, F. Yan, T. Orlando, S. Gustavsson, W. Oliver,  _Applied Physics Reviews_ 2019, _6_ 021318\.
  * [4] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, W. D. Oliver,  _Annual Review of Condensed Matter Physics_ 2020, _11_ , 1 369.
  * [5] M. D. Hutchings, et al,  _Phys. Rev. Appl._ 2017, _8_ 1\.
  * [6] J.-J. Chang, D. J. Scalapino,  _Physical Review B_ 1977, _15_ , 5 2651.
  * [7] D. J. Goldie, S. Withington,  _Superconductor Science and Technology_ 2013, _26_ 015004\.
  * [8] P. J. de Visser, D. J. Goldie, P. Diener, S. Withington, J. J. A. Baselmans, T. M. Klapwijk,  _Physical Review Letters_ 2014, _112_ , 4 047004.
  * [9] R. P. Budoyo,  Ph.D. thesis, University of Maryland, College Park, USA, 2015,  URL https://doi.org/10.13016/M2TD9N.
  * [10] A. P. Vepsäläinen, A. H. Karamlou, J. L. Orrell, A. S. Dogra, B. Loer, F. Vasconcelos, D. K. Kim, A. J. Melville, B. M. Niedzielski, J. L. Yoder, S. Gustavsson, J. A. Formaggio, B. A. VanDevender, W. D. Oliver,  _Nature_ 2020, _584_ , 7822 551.
  * [11] P. J. de Visser, J. J. A. Baselmans, P. Diener, S. J. C. Yates, A. Endo, T. M. Klapwijk,  _Physical Review Letters_ 2011, _106_ , 16 167004.
  * [12] R. Barends, J. Wenner, M. Lenander, Y. Chen, R. C. Bialczak, J. Kelly, E. Lucero, P. O’Malley, M. Mariantoni, D. Sank, H. Wang, T. C. White, Y. Yin, J. Zhao, A. N. Cleland, J. M. Martinis, J. J. A. Baselmans,  _Applied Physics Letters_ 2011, _99_ , 11 113507.
  * [13] R. P. Budoyo, J. B. Hertzberg, C. J. Ballard, K. D. Voigt, J. R. A. Z. Kim, C. J. Lobb, F. C. Wellstood,  _Phys. Rev. B_ 2016, _93_ 024514\.
  * [14] N. P. de Leon, K. M. Itoh, D. Kim, K. K. Mehta, T. E. Northup, H. Paik, B. S. Palmer, N. Samarth, S. Sangtawesin, D. W. Steuerman,  _Science_ 2021, _372_ , 6539 eabb2823.
  * [15] P. W. Anderson, B. I. Halperin, c. M. Varma,  _The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics_ 1972, _25_ , 1 1.
  * [16] W. A. Phillips,  _Journal of Low Temperature Physics_ 1972, _7_ , 3 351\.
  * [17] M. Von Schickfus, S. Hunklinger,  _Physics Letters A_ 1977, _64_ , 1 144.
  * [18] J. L. Black, B. I. Halperin,  _Physical Review B_ 1977, _16_ , 6 2879.
  * [19] J. L. Black,  _Physical Review B_ 1978, _17_ , 6 2740.
  * [20] W. A. Phillips,  _Reports on Progress in Physics_ 1987, _50_ , 12 1657\.
  * [21] J. M. Martinis, K. B. Cooper, R. McDermott, M. Steffen, M. Ansmann, K. D. Osborn, K. Cicak, S. Oh, D. P. Pappas, R. W. Simmonds, C. C. Yu,  _Physical Review Letters_ 2005, _95_ , 21 210503.
  * [22] J. Gao, J. Zmuidzinas, B. A. Mazin, H. G. LeDuc, P. K. Day,  _Applied Physics Letters_ 2007, _90_ , 10 102507.
  * [23] J. Gao, M. Daal, A. Vayonakis, S. Kumar, J. Zmuidzinas, B. Sadoulet, B. A. Mazin, P. K. Day, H. G. Leduc,  _Applied Physics Letters_ 2008, _92_ , 15 152505.
  * [24] A. D. O’Connell, M. Ansmann, R. C. Bialczak, M. Hofheinz, N. Katz, E. Lucero, C. McKenney, M. Neeley, H. Wang, E. M. Weig, A. N. Cleland, J. M. Martinis,  _Applied Physics Letters_ 2008, _92_ , 11 112903.
  * [25] S. Kumar, J. Gao, J. Zmuidzinas, B. A. Mazin, H. G. LeDuc, P. K. Day,  _Applied Physics Letters_ 2008, _92_ , 12 123503.
  * [26] R. Barends, H. L. Hortensius, T. Zijlstra, J. J. A. Baselmans, S. J. C. Yates, J. R. Gao, T. M. Klapwijk,  _Applied Physics Letters_ 2008, _92_ , 22 223502.
  * [27] P. Macha, S. H. W. v. d. Ploeg, G. Oelsner, E. Il’ichev, H.-G. Meyer, S. Wünsch, M. Siegel,  _Applied Physics Letters_ 2010, _96_ , 6 062503.
  * [28] A. L. Burin, A. O. Maksymov,  _Physical Review B_ 2018, _97_ , 21 214208.
  * [29] C. Müller, J. H. Cole, J. Lisenfeld,  _Reports on Progress in Physics_ 2019, _82_ , 12 124501\.
  * [30] C. C. Yu, H. M. Carruzzo,  Two-level systems and the tunneling model: A critical view, 2021,  URL https://arxiv.org/abs/2101.02787.
  * [31] A. M. Holder, K. D. Osborn, C. J. Lobb, C. B. Musgrave,  _Physical Review Letters_ 2013, _111_ , 6 065901.
  * [32] L. Gordon, H. Abu-Farsakh, A. Janotti, C. G. V. de Walle,  _Scientific Reports_ 2014, _4_ , 1 7590.
  * [33] W. D. Oliver, P. B. Welander,  _MRS Bulletin_ 2013, _38_ , 10 816.
  * [34] A. Bruno, G. de Lange, S. Asaad, K. L. van der Enden, N. K. Langford, L. DiCarlo,  _Appl. Phys, Lett._ 2015, _106_ 182601\.
  * [35] G. Calusine, A. Melville, W. Woods, R. Das, C. Stull, V. Bolkhovsky, D. H. D. Braje1, D. K. Kim, X. Miloshi, D. Rosenberg, A. Sevi, J. L. Yoder, E. Dauler, W. D. Oliver,  _Appl. Phys. Lett._ 2018, _112_ , 6 1.
  * [36] A. Melville, G. Calusine, W. Woods, K. Serniak, E. Golden, B. M. Niedzielski, D. K. Kim, A. Sevi, J. L. Yoder, E. A. Dauler, W. D. Oliver,  _Applied Physics Letters_ 2020, _117_ , 12 124004.
  * [37] J. M. Sage, V. Bolkhovsky, W. D. Oliver, B. Turek, P. B. Welander,  _J. Appl. Phys._ 2011, _109_ , 12 063915.
  * [38] B. Chiaro, A. Megrant, A. Dunsworth, Z. Chen, R. Barends, B. Campbell, Y. Chen, A. Fowler, I. C. Hoi, E. Jeffrey, J. Kelly, J. Mutus, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White1, J. M. Martinis,  _Supercond. Sci. Technol._ 2016, _29_ , 12 104006.
  * [39] C. Richardson, A. Alexander, C. G. Weddle, B. Arey, M. Olszta,  _J. Appl. Phys._ 2020, _127_ 235302\.
  * [40] A. Romanenko, R. Pilipenko, S. Zorzetti, D. Frolov, M. Awida, S. Belomestnykh, S. Posen, A. Grassellino,  _Phys. Rev. Appl._ 2020, _13_ 034032\.
  * [41] S. E. de Graaf, L. Faoro, L. B. Ioffe, S. Mahashabde, J. J. Burnett, T. Lindström, S. E. Kubatkin, A. V. Danilov, A. Y. Tzalenchuk,  _Science Advances_ 2020, _6_ , 51 eabc5055.
  * [42] J. Jäckle,  _Zeitschrift für Physik A Hadrons and nuclei_ 1972, _257_ , 3 212.
  * [43] S. Hunklinger, W. Arnold,  _3 - Ultrasonic Properties of Glasses at Low Temperatures_ , volume 12, 155–215,  Academic Press,  ISBN 0893-388X, 1976.
  * [44] C. Du, H. Chen, S. Li,  _Physical Review B_ 2006, _74_ , 11 113105.
  * [45] A. L. Rakhmanov, A. M. Zagoskin, S. Savel’ev, F. Nori,  _Physical Review B_ 2008, _77_ , 14 144507.
  * [46] P. Jung, A. V. Ustinov, S. M. Anlage,  _Superconductor Science and Technology_ 2014, _27_ , 7 073001.
  * [47] A. M. Zagoskin, D. Felbacq, E. Rousseau,  _EPJ Quantum Technology_ 2016, _3:2_ , 1.
  * [48] N. Lazarides, G. P. Tsironis,  _Physics Reports_ 2018, _752_ 1\.
  * [49] J. Q. You, F. Nori,  _Nature_ 2011, _474_ , 7353 589.
  * [50] P. Macha, G. Oelsner, J.-M. Reiner, M. Marthaler, S. André, G. Schön, U. Hübner, H.-G. Meyer, E. Il’ichev, A. V. Ustinov,  _Nat Commun_ 2014, _5:5146_.
  * [51] K. V. Shulga, E. Il’ichev, M. V. Fistul, I. S. Besedin, S. Butz, O. V. Astafiev, U. Hübner, A. V. Ustinov,  _Nature Communications_ 2018, _9_ , 1 150.
  * [52] Z. L. Xiang, S. Ashhab, J. Q. You, F. Nori,  _Rev. Mod. Phys_ 2013, _85_ , 2 623.
  * [53] S. I. Mukhin, M. V. Fistul,  _Superconductor Science and Technology_ 2013, _26_ , 8 084003.
  * [54] P. A. Volkov, M. V. Fistul,  _Physical Review B_ 2014, _89_ , 5 054507.
  * [55] P. J. Petersan, S. M. Anlage,  _Journal of Applied Physics_ 1998, _84_ , 6 3392.
  * [56] S. J. Weber, K. W. Murch, D. H. Slichter, R. Vijay, I. Siddiqi,  _Applied Physics Letters_ 2011, _98_ , 17 172510.
  * [57] J. Gao,  PhD dissertation, Caltech, Department of Physics, Caltech, 2008,  URL https://resolver.caltech.edu/CaltechETD:etd-06092008-235549.
  * [58] C. R. H. McRae, H. Wang, J. Gao, M. R. Vissers, T. Brecht, A. Dunsworth, D. P. Pappas, J. Mutus,  _Review of Scientific Instruments_ 2020, _91_ , 9 091101\.
  * [59] K. D. Crowley, R. A. McLellan, A. Dutta, N. Shumiya, A. P. M. Place, X. H. Le, Y. Gang, T. Madhavan, N. Khedkar, Y. C. Feng, E. A. Umbarkar, X. Gui, L. V. H. Rodgers, Y. Jia, M. M. Feldman, S. A. Lyon, M. Liu, R. J. Cava, A. A. Houck, N. P. de Leon,  Disentangling Losses in Tantalum Superconducting Circuits, 2023,  URL http://arxiv.org/abs/2301.07848,  ArXiv:2301.07848 [cond-mat, physics:quant-ph].
  * [60] J. Burnett, L. Faoro, T. Lindstrom,  _Superconductivity Science and Technology_ 2016, _29_ 044008\.
  * [61] J. Gao, J. Zmuidzinas, A. Vayonakis, P. Day, B. Mazin, H. Leduc,  _J Low Temp Phys_ 2008, _151_ , 1 557.
  * [62] D. P. Pappas, M. R. Vissers, D. S. Wisbey, J. S. Kline, J. Gao,  _IEEE Transactions on Applied Superconductivity_ 2011, _21_ , 3 871.
  * [63] H. Wang, M. Hofheinz, J. Wenner, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, M. Weides, A. N. Cleland, J. M. Martinis,  _Applied Physics Letters_ 2009, _95_ , 23 233508.
  * [64] M. S. Khalil, F. C. Wellstood, K. D. Osborn,  _IEEE Transactions on Applied Superconductivity_ 2011, _21_ , 3 879.
  * [65] M. Schechter, P. C. E. Stamp,  _Journal of Physics: Condensed Matter_ 2008, _20_ , 24 244136.
  * [66] M. S. Khalil, M. J. A. Stoutimore, S. Gladchenko, A. M. Holder, C. B. Musgrave, A. C. Kozen, G. Rubloff, Y. Q. Liu, R. G. Gordon, J. H. Yum, S. K. Banerjee, C. J. Lobb, K. D. Osborn,  _Applied Physics Letters_ 2013, _103_ , 16 162601.
  * [67] M. Schechter, P. C. E. Stamp,  _Physical Review B_ 2013, _88_ , 17 174202, publisher: American Physical Society.
  * [68] M. Schechter, P. Nalbach, A. L. Burin,  _New Journal of Physics_ 2018, _20_ , 6 063048, publisher: IOP Publishing.
  * [69] L. Yu, S. Matityahu, Y. J. Rosen, C.-C. Hung, A. Maksymov, A. L. Burin, M. Schechter, K. D. Osborn,  _Scientific Reports_ 2022, _12_ , 1 16960, number: 1 Publisher: Nature Publishing Group.
  * [70] J. Burnett, L. Faoro, I. Wisby, V. L. Gurtovoi, A. V. Chernykh, G. M. Mikhailov, V. A. Tulin, R. Shaikhaidarov, V. Antonov, P. J. Meeson, A. Y. Tzalenchuk, T. Lindström,  _Nature Communications_ 2014, _5_ 4119\.
  * [71] L. Faoro, L. B. Ioffe,  _Physical Review Letters_ 2012, _109_ , 15 157005.
  * [72] L. Faoro, L. B. Ioffe,  _Physical Review B_ 2015, _91_ , 1 014201.
  * [73] E. Cuevas, R. Chicón, M. Ortuño,  _Physica B: Condensed Matter_ 1989, _160_ , 3 293.
  * [74] R. Jankowiak, G. J. Small,  _Chemical Physics Letters_ 1993, _207_ , 4 436.
  * [75] J. Gao, M. Daal, J. M. Martinis, A. Vayonakis, J. Zmuidzinas, B. Sadoulet, B. A. Mazin, P. K. Day, H. G. Leduc,  _Applied Physics Letters_ 2008, _92_ , 21 212504.
  * [76] S. Kumar, J. Gao, J. Zmuidzinas, B. A. Mazin, H. G. LeDuc, P. K. Day,  _Applied Physics Letters_ 2008, _92_ , 12 123503.
  * [77] S. E. de Graaf, S. Mahashabde, S. E. Kubatkin, A. Y. Tzalenchuk, A. V. Danilov,  _Physical Review B_ 2021, _103_ , 17 174103.
  * [78] B. Sarabi,  Ph.D. thesis, University of Maryland, College Park, 2014,  URL https://drum.lib.umd.edu/handle/1903/16186,  Accepted: 2015-02-06T06:41:40Z.
  * [79] J. Lisenfeld, C. Müller, J. H. Cole, P. Bushev, A. Lukashenko, A. Shnirman, A. V. Ustinov,  _Physical Review Letters_ 2010, _105_ , 23 230504.
  * [80] G. J. Grabovskij, T. Peichl, J. Lisenfeld, G. Weiss, A. V. Ustinov,  _Science_ 2012, _338_ , 6104 232.
  * [81] B. Sarabi, A. N. Ramanayaka, A. L. Burin, F. C. Wellstood, K. D. Osborn,  _Physical Review Letters_ 2016, _116_ , 16 167002.
  * [82] C.-C. Hung, L. Yu, N. Foroozani, S. Fritz, D. Gerthsen, K. D. Osborn,  _Phys. Rev. Applied_ 2022, _17_ 034025\.
  * [83] M. Bhattacharya, K. D. Osborn, A. Mizel,  _Physical Review B_ 2011, _84_ , 10 104517.
  * [84] C. Yu, A. Leggett,  _Comments Cond. Mat. Phys_ 1988, _14_ , 4 231.
  * [85] J. Lisenfeld, G. J. Grabovskij, C. Müller, J. H. Cole, G. Weiss, A. V. Ustinov,  _Nature Communications_ 2015, _6_ , 1 6182.
  * [86] S. N. Coppersmith,  _Physical Review Letters_ 1991, _67_ , 17 2315.
  * [87] A. L. Burin, Y. Kagan,  _Physica B: Condensed Matter_ 1994, _194-196_ 393\.
  * [88] D. C. Vural, A. J. Leggett,  _Journal of Non-Crystalline Solids_ 2011, _357_ , 19 3528.
  * [89] D. C. Mattis, J. Bardeen,  _Physical Review_ 1958, _111_ , 2 412.
  * [90] J. P. Turneaure, J. Halbritter, H. A. Schwettman,  _Journal of Superconductivity_ 1991, _4_ , 5 341.
  * [91] A. Rothwarf, B. N. Taylor,  _Physical Review Letters_ 1967, _19_ , 1 27.
  * [92] W. H. Parker,  _Physical Review B_ 1975, _12_ 3667\.
  * [93] T. Guruswamy, D. J. Goldie, S. Withington,  _Superconductor Science and Technology_ 2015, _28_ , 5 054002.
  * [94] L. Grünhaupt, N. Maleeva, S. T. Skacel, M. Calvo, F. Levy-Bertrand, A. V. Ustinov, H. Rotzinger, A. Monfardini, G. Catelani, I. M. Pop,  _Physical Review Letters_ 2018, _121_ , 11 117001.

Supplemental Material for Anomalous Loss Reduction Below Two-Level System
Saturation in Aluminum Superconducting Resonators

##  S1 Experimental Details

Fig. S1 schematically shows the setup of the cryogenic transmission
measurement with the VNA to characterize the resonator response. A variety of
attenuators (produced by XMA) are used on each cryogenic stage to thermalize
the center conductors of the coaxial cables. The total attenuation in the
input line is -66 dB. Both the input line and output line on either side of
the device are filtered by commercial microwave low-pass filters. The input
line has a Marki Microwave low pass filter (FLP-1460) with 3-dB cutoff
frequency at 14.6 GHz and the output line has another Marki Microwave low pass
filter (FLP-1250) with 3-dB cutoff frequency at 12.5 GHz. The output line goes
through the cryogenic isolator (QUINSTAR Technology QCI-G0301201AM) with
working frequency band 3-12 GHz, and then the signal is amplified by 36 dB
using a commercial high-electron mobility transistor amplifier (Low Noise
Factory LNF-LNC0.3_14A with typical noise temperature 4.2 K) at the 4K stage,
and then the transmitted signal is further amplified by 37 dB using another
room temperature amplifier (HEMT) (Low Noise Factory LNF-LNC2_6A with typical
noise temperature 50 K at ambient temperature).

Fig. S1: Schematic of the microwave measurement setup for the study of the
aluminum $\lambda/2$ resonator. The VNA at room temperature sends a signal
from port 1 to the cryostat. The signal is attenuated at each stage of the
cryostat before passing through the low-pass filter (LPF) and entering the
device under test (DUT). The DUT is surrounded by a Cryoperm magnetic shield.
The output signal also passes through a low-pass filter before going through
0-dB attenuators that thermalize the coaxial cable center conductor. The
signal passes through an isolator and is amplified at the 4 K stage and at
room temperature, before entering the VNA in port 2.

##  S2 Conventional model for TLS loss

The TLS formalism is based on a simple model for a single TLS that can be
described by the Hamiltonian,
$H_{TLS}=\frac{1}{2}\begin{pmatrix}-\Delta&\Delta_{0}\\\
\Delta_{0}&\Delta\end{pmatrix}$ , where $\Delta$ is the asymmetry of the
double well potential and $\Delta_{0}$ is the tunneling barrier energy between
the potential wells [20, 9]. Here
$\varepsilon=\sqrt{\Delta^{2}+\Delta_{0}^{2}}$ is the energy splitting of the
two eigenvalues of $H_{TLS}$. A typical resonator hosts an ensemble of TLS
with different values of $\Delta$ and $\Delta_{0}$. In the standard tunneling
model (STM), a uniform distribution of the ensemble in $\Delta$ and log
uniform distribution in $\Delta_{0}$ is assumed [20]. Whereas in the
generalized tunneling model (GTM), the distribution takes on an additional
weak dependence of $\Delta$ due to a logarithmic correction introduced by the
strong TLS interactions [73, 72]. The distribution function can be
mathematically expressed as

$\displaystyle
P(\Delta,\Delta_{0})=(1+\mu)\left(\frac{\Delta}{\Delta_{max}}\right)^{\mu}\frac{P_{0}}{\Delta_{0}}\quad\begin{cases}\mu=0&\text{no
TLS interactions}\\\ \mu\approx 0.3&\text{strong TLS interactions}\end{cases}$
(S1)

where $P_{0}\approx 10^{44}J^{-1}m^{-3}$ is the density of states for TLS,
$\Delta_{max}\sim\varepsilon_{max}\sim 100K$ [72], and the exponent $\mu$
approximately characterizes the logarithmic correction to the density of
states of TLS due to strong TLS interactions.

The dynamics of a single TLS can be described by the linearized Bloch
equations of the pseudospin $\vec{S}(t)=\vec{S}^{0}(t)+\vec{S}^{1}(t)$, where
$\vec{S}^{0}$ is the solution to the homogeneous system without the external
field, and $\vec{S}^{1}$ is the linear solution with frequency $\omega$ of the
driving field. Under the rotating wave approximation, the Bloch equations
become: [20, 57, 71]

$\displaystyle i\frac{d\langle S^{+}\rangle}{dt}=\Omega\langle
S^{0}_{z}\rangle-(\omega-\varepsilon/\hbar+i\Gamma_{2})\langle S^{+}\rangle$
(S2) $\displaystyle\frac{d\langle
S^{0}_{z}\rangle}{dt}=\Omega\imaginary\langle
S^{+}\rangle-\Gamma_{1}\left(\langle S^{0}_{z}\rangle-m)\right)$ (S3)

with the Rabi frequency $\Omega$ characterizing the absorbed power,
$\Omega\propto\sqrt{P_{ab}}$, and is defined as

$\Omega=\frac{2d_{0}\Delta_{0}}{\hbar\varepsilon}|\vec{E}|$ (S4)

where $S^{+}=S^{1}_{x}+iS^{1}_{y}$, $\Gamma_{1}$ and $\Gamma_{2}$ are the two
phenomenological rates that describe the longitudinal ($S_{z}$) and transverse
($S_{x,y}$) relaxations, and $m=\tanh(\varepsilon/(2k_{B}T))/2$ is the
equilibrium value of the $\langle S_{z}^{0}\rangle$. In STM, the dielectric
response of a single TLS can be obtained from the stationary solution to
$\langle S^{+}\rangle$ [15, 17, 71]:

$\langle
S^{+}\rangle=\frac{m\Omega(\omega-\varepsilon/\hbar-i\Gamma_{2})}{(\omega-\varepsilon/\hbar)^{2}+\Gamma_{2}^{2}(1+\Omega^{2}\Gamma_{1}^{-1}\Gamma_{2}^{-1})}$
(S5)

The single TLS loss corresponds to the imaginary part of the response function
in Eq.(S5) which is in the form of a Lorentzian in $\varepsilon/\hbar$
centered at $\omega$ with a width:

$\displaystyle w=\Gamma_{2}\sqrt{1+\kappa}$ (S6)

where $\kappa=\Omega^{2}\Gamma_{1}^{-1}\Gamma_{2}^{-1}$. The total dielectric
loss is simply the integral of the single TLS contribution Eq.(S5) over the
distribution of the TLS Eq.(S1) [17, 20, 57].

$\delta_{TLS}=\frac{1}{\epsilon_{r}\epsilon_{0}}\int\int\int
P(\varepsilon,\Delta_{0})\left(\frac{\Delta_{0}d_{0}}{\varepsilon}\right)^{2}\frac{\cos^{2}{\theta}}{\hbar}\frac{m\Gamma_{2}}{\Gamma_{2}^{2}(1+\kappa)+(\varepsilon/\hbar-\omega)^{2}}d\varepsilon
d\Delta_{0}d\theta$ (S7)

where $\epsilon_{r}\epsilon_{0}$ is the permittivity of the host dielectric
material, $P(\varepsilon,\Delta_{0})$ is obtained from Eq. (S1) with a change
of variable from $\Delta$ to $\varepsilon$, $d_{0}$ is the maximum transition
electric dipole moment of the TLS with energy splitting $\varepsilon$,
$\vec{E}$ is the applied microwave electric field on the TLS dipole, and
$\theta$ is the angle between the applied electric field and the TLS dipole
moment.

For a typical TLS with $\varepsilon/h\approx$ 5 GHz and at reasonably low
temperatures and powers, the width of its response $w\approx\Gamma_{2}\sim 1\
\text{MHz}\ll\omega$. Due to this sharp Lorentzian response function, the
total loss is dominated by the resonant TLS whose energies
$\varepsilon\sim\hbar\omega$. Before analytically evaluating the integral in
Eq.(S7), the expressions for $\Gamma_{1,2}$ need to be introduced.

The longitudinal relaxation rate $(\Gamma_{1})$ of a single TLS is dominated
by the phonon process: [42, 43, 18, 57]

$\displaystyle\Gamma_{1}=\left(\frac{\Delta_{0}}{\varepsilon}\right)^{2}\left[\frac{\gamma_{L}^{2}}{v_{L}^{5}}+\frac{2\gamma_{T}^{2}}{v_{T}^{5}}\right]\frac{\varepsilon^{3}}{2\pi\rho\hbar^{4}}\coth(\frac{\varepsilon}{2k_{B}T})=\left(\frac{\Delta_{0}}{\varepsilon}\right)^{2}\Gamma_{1}^{max}$
(S8)

where $\gamma_{L}$ and $\gamma_{T}$ are the longitudinal and transverse
deformation potentials, respectively, $v_{L}$ and $v_{T}$ are the longitudinal
and transverse sound velocities, $\rho$ is the mass density, and
$\Gamma_{1}^{max}$ is the maximum $\Gamma_{1}$ for the TLS with energy
splitting $\varepsilon$, when $\Delta_{0}=\varepsilon$.

The transverse relaxation rate $(\Gamma_{2})$ is defined as

$\displaystyle\Gamma_{2}=\Gamma_{2}^{ph}+\Gamma_{ds}$
$\displaystyle\text{where}\quad\Gamma_{2}^{ph}=\Gamma_{1}/2\quad\text{and}\quad\Gamma_{ds}\sim$
$\displaystyle
10^{-3}(k_{B}T/\varepsilon_{max})^{\mu}k_{B}T/\hbar\quad\text{\cite[cite]{[\@@bibref{}{jankowiak_origin_1993,faoro_interacting_2015}{}{}]}}$
(S9)

$\Gamma_{ds}$ is the dephasing rate of the resonant TLS energy level
$\varepsilon$, caused by its interactions with thermally activated TLS whose
$\varepsilon\lesssim k_{B}T$, and is valid for low temperature measurement
($T<1$ K) [74].

After substituting the expressions for $\Gamma_{1}$, $\Gamma_{2}$ and $\Omega$
from Eqs. (S8, S9, S4) into the integral for the loss Eq. (S7) and using
$\mu=0$, the famous STM prediction of TLS loss is obtained, [43]

$\delta_{TLS}=\frac{\pi
P_{0}d_{0}^{2}}{3\epsilon_{r}\epsilon_{0}}\frac{\tanh{\frac{\hbar\omega}{2k_{B}T}}}{\sqrt{1+(\Omega/\Omega^{c})^{2}}}$
(S10)

where $\Omega^{c}\propto\sqrt{\Gamma_{1}^{max}\Gamma_{2}}$ is the critical
Rabi frequency that characterizes the saturation of TLS. The loss is expected
to have an inverse square root dependence on power after the TLS saturation,
$\delta_{TLS}\sim\Omega\propto P_{ab}^{-0.5}$ for $\Omega\gg\Omega^{c}$, which
is much faster than observed in the data in Fig. 4. This leads us to look at
the generalized tunneling model with the inclusion of fluctuators.

##  S3 Master equation formalism for fluctuators

The dynamics of a single TLS is described in the Bloch equations [20, 57, 71]:

$\displaystyle i\frac{d\langle S^{+}\rangle}{dt}=\Omega\langle
S^{0}_{z}\rangle-(\omega-(\varepsilon+\xi(t))/\hbar+i\Gamma_{2})\langle
S^{+}\rangle$ $\displaystyle\frac{d\langle
S^{0}_{z}\rangle}{dt}=\Omega\imaginary\langle
S^{+}\rangle-\Gamma_{1}\left(\langle S^{0}_{z}\rangle-m)\right)$

The solution to the Bloch equations depend on the relation between the
transition rates induced by fluctuators $\gamma$ and the other inherent
dynamics: $\Omega,\Gamma_{1,2},\omega-\varepsilon/\hbar$. We structure this
section according to the three different regimes of the rate $\gamma$, namely
high $\gamma>\Omega,\Gamma_{2}$, intermediate $\gamma\approx\Omega\text{ or
}\Gamma_{2}$, and low $\gamma<\Gamma_{1}$.

###  S3.1 Low $\gamma$ fluctuators

For low $\gamma$ fluctuators, the single TLS response is given by the
stationary solution to Eq.(S2,S3), since the TLS dynamics is much faster than
the jump rates $\gamma$. We can further separate the analysis for two
different types: the weakly-coupled fluctuator which shifts the energy level
of TLS by a small energy $\xi$ and cannot bring a TLS in or out of resonance,
and the strongly-coupled fluctuator whose $\xi$ is large enough to achieve
this. The weakly-coupled fluctuators induce detunings that follow a Lorentzian
distribution with width
$\Gamma_{f}\propto(k_{B}T/\hbar)(k_{B}T/\varepsilon_{max})^{\mu}\lesssim\Gamma_{2}$[72]:

$\displaystyle
P(\xi(t))=\frac{1}{\pi}\frac{\Gamma_{f}}{\Gamma_{f}^{2}+\xi(t)^{2}}$ (S11)

The single TLS imaginary response under a single fluctuator can be obtained
from Eq.( S5) by replacing $\varepsilon$ with $\varepsilon+\xi(t)$ [72]:

$\imaginary\langle
S^{+}\rangle=\frac{-m\Omega\Gamma_{2}}{(\omega-(\varepsilon+\xi(t))/\hbar)^{2}+\Gamma_{2}^{2}(1+\kappa)}$
(S12)

where $\kappa=\Omega^{2}/(\Gamma_{1}\Gamma_{2})$. Integrating Eq.(S12) over
the distribution of detuning $P(\xi(t))$, one can obtain the single TLS loss
[72]:

$\delta(\varepsilon)\propto\frac{m}{\sqrt{1+\kappa}}\frac{\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f}}{(\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f})^{2}+(\omega-\varepsilon/\hbar)^{2}}$
(S13)

The second fraction is in the form of a Lorentzian with width
$\Gamma_{2}\sqrt{1+\kappa}+\Gamma_{f}$, and the first term is the same as the
STM formula. If a continuous distribution function for TLS
$P(\varepsilon,\Delta_{0})$ is assumed, one can obtain the total internal loss
by integrating Eq.(S13) over $P(\varepsilon,\Delta_{0})$, resulting in the
same internal loss as in STM.

This distribution $P(\xi(t))$ may not apply to the strongly-coupled
fluctuators which can completely shift the TLS in and out of resonance, since
there are typically only a few strong fluctuators in surface dielectrics [72].
They contribute to the imaginary part of a single TLS response as a slow
random telegraph noise (i.e. $\xi(t)$ in Eq.(S12) ). The internal loss will
still be in the form of the STM loss Eq.(S5). Since our loss measurement is
obtained by a slow frequency sweep (typically around 15 minutes) on the
network analyzer through the resonance, the effect of the few strongly-coupled
fluctuators is likely averaged out.

Three methods were used to estimate the average effect of the low $\gamma$
fluctuators on the internal loss. Intuitively, we could estimate the total
loss by summing up the contribution from all the TLS near the resonance. The
key assumption is the low density of the resonant TLS which justifies the
discrete treatment. An estimation based on [78] yields $\sim 1$ TLS in the
bandwidth of the resonance for a volume of TLS-inhabiting dielectrics around
30 $\mu m^{3}$, which supports our assumption. The first method is thus a
summation of all the STM single TLS loss in the form of Eq.(S13) near the
resonance with
$\varepsilon=...,\hbar(\omega+\nu)-\Delta\varepsilon,\hbar(\omega+\nu),\hbar(\omega+\nu)+\Delta\varepsilon,...$
given in Eq.(15). A finite detuning $\hbar\nu<\Delta\varepsilon$ is included
such that the TLS closest to resonance will not give a diverging response as
$T\rightarrow 0$, and instead contributes to the loss as
$m\Gamma_{2}\Delta\varepsilon/(\hbar\nu^{2})\propto T^{1+\mu}$ at low
temperature. This detuning is a consequence of a sparse TLS distribution where
$\Delta\varepsilon/\hbar\gtrsim\Gamma_{2}$. Thus, any given TLS is rarely on
resonance due to the low density of states. We should emphasize that this
treatment is independent from the master equation formalism for the
fluctuators, and can be applied to STM with low TLS density where the use of
distribution function $P(\varepsilon,\Delta_{0})$ is inappropriate.

The second method assumes that the detuning $\nu$ is not much larger than
$\Gamma_{2}$ so that the distribution $P(\xi(t))$ for the weakly-coupled
fluctuators can still be applied to the case $\xi\sim\hbar\nu$. The
probability of a TLS with a detuning $\nu$ to be on resonance under the
influence of the low $\gamma$ fluctuators can be calculated by integrating
$P(\xi(t))$:

$P_{\text{on
resonance}}=\left.\frac{1}{\pi}\arctan(\frac{\xi}{\hbar\Gamma_{f}})\right\rvert_{\hbar(\nu-\Gamma_{2})}^{\hbar(\nu+\Gamma_{2})}$
(S14)

which also leads to a monotonic increasing temperature dependence below 50mK
as in the data.

The third method estimates the number of strongly-coupled fluctuators capable
of moving the coherent TLS in and out of resonance. Since the energy drift
caused by fluctuators is directly related to the interaction energy, which is
dipole-dipole like, $U(r)\sim r^{-3}$, we could find a volume around the near-
resonant TLS which hosts fluctuators that interact strongly enough. For
$\Gamma_{2}<\nu$, the TLS is detuned, and the criterion for being a strong
fluctuator is that the induced energy shift can bring the TLS in resonance
such that $\hbar(\nu-\Gamma_{2})<U(r)<\hbar(\nu+\Gamma_{2})$. This volume is a
spherical shell with inner and outer radius defined by the bounds on the
interaction energy $U(r)$. For $\Gamma_{2}>\nu$, the TLS is already in
resonant, and the strong fluctuators are those that move the TLS out of
resonance, with $U(r)>\hbar(\Gamma_{2}-\nu)$. This volume is a sphere with
radius defined by $U(r)=\hbar(\Gamma_{2}-\nu)$. The total number of strong
fluctuators can then be estimated as:

$N_{\text{strong
fluctuators}}=(1+\mu)P_{0}\frac{k_{B}T}{\varepsilon_{max}}^{\mu}V_{\text{strong
fluctuators}}k_{B}T\propto\quad\begin{cases}\frac{2\Gamma_{2}^{2}}{\nu^{2}-\Gamma_{2}^{2}}\quad(\Gamma_{2}<\nu)\\\
\frac{\Gamma_{2}}{\Gamma_{2}-\nu}\quad(\Gamma_{2}>\nu)\end{cases}$ (S15)

where the density of states for TLS from Eq.(S1) is used, and the energy range
is set to the thermal energy $k_{B}T$. The two expressions Eq.(S14,S15) are
plotted below as a function of $\Gamma_{2}\propto T^{1+\mu}$.

Fig. S2:  The blue curve is the probability for a TLS detuned with $\nu$ to be
resonant under the influence of weakly-coupled fluctuators. The yellow curve
is the normalized estimated number of strongly-coupled fluctuators. They are
plotted as a function of $\Gamma_{2}\propto T^{1.3}$, which also represents
their temperature dependence. It can be concluded that for both the strong and
weak fluctuators, as temperature decreases ($\Gamma_{2}/\nu\rightarrow 0$) the
TLS are increasingly detuned.

The x axis of Fig. S2, $\Gamma_{2}$ is normalized by the detuning $\nu$. When
$\Gamma_{2}$ increases to $\nu$, the probability for TLS to be resonant under
the weak fluctuators (Blue curve in Fig. S2) approaches a constant
$2/\pi\arctan(\Gamma_{2}/\Gamma_{f})$, which is denoted by the horizontal line
where $\Gamma_{f}=0.8\Gamma_{2}$. The interpretation of the number of strong
fluctuators (Yellow curve in Fig.S2) is more nuanced. When $\Gamma_{2}<\nu$,
the coherent TLS itself is detuned, the strong fluctuators shift the TLS in
resonance, and therefore their number is positively correlated with the TLS
response. However, for $\Gamma_{2}>\nu$, the strong fluctuators shift the
already resonant TLS out of resonance, and thus their number is negatively
correlated with the TLS response. Consequently, as $\Gamma_{2}$ or temperature
increases, the resonant response becomes stronger. The divergence near
$\Gamma_{2}\sim\nu$ represents the failure of this estimation in the range
where any arbitrarily small shift can move the TLS in and out of resonance,
and the strong fluctuator volume goes to infinity. In reality, this volume is
bounded by the host material for TLS. For $\Gamma_{2}\gg\nu$, the case of TLS
exactly on resonance is recovered, and the number of strong fluctuators is
clearly independent of temperature. In both the probability and the number
estimates, the detuned TLS becomes resonant more often as $\Gamma_{2}$ or
temperature increases, which qualitatively explains the increase in loss in
temperature at low temperatures in the measurement.

It should be noted that if applied to TLS exactly on resonance, the above
arguments will lead to the conclusion that low $\gamma$ fluctuators present no
temperature dependence for both the strong and weak low $\gamma$ fluctuators,
which forms the basis of the claim that no anomalous temperature dependence on
the loss is expected in the original GTM work [72].

###  S3.2 High $\gamma$ fluctuators

For fluctuators that induce fast jump rates on the resonant TLS,
$\gamma\gtrsim\Omega,\Gamma_{1,2}$, their effects on the dynamics are
described by the master equation [71]:

$\displaystyle\frac{\partial\rho}{\partial
t}+\frac{d}{d\vec{S}}(\frac{d\vec{S}}{dt}\rho)=\gamma[\delta(\langle
S_{z}^{0}\rangle-m)\delta(\langle S_{x}^{1}\rangle)\delta(\langle
S_{y}^{1}\rangle)-\rho]$ (S16)

where $\vec{S}=(\langle S_{x}^{1}\rangle,\langle S_{y}^{1}\rangle,\langle
S_{z}^{0}\rangle)$ and $\rho$ is the probability for a TLS to be in the state
$\vec{S}$. The case for high $\gamma$ fluctuators is particularly simple,
where the fast random telegraph noise is averaged out. Thus, $d\vec{S}/dt$ is
solved from Eq.(S2,S3) with $\xi(t)=0$. The average y-component response of
the single TLS with energy splitting $\varepsilon$ can then be solved from the
master equation Eq.(S16) as:

$\overline{\langle S_{y}^{1}\rangle}=\int\rho\langle S_{y}^{1}\rangle
d\vec{S}=\frac{m\Omega}{\sqrt{1+\frac{\Omega^{2}}{(\gamma+\Gamma_{1})(\gamma+\Gamma_{2})}}}\frac{\sqrt{(\gamma+\Gamma_{2})^{2}+\Omega^{2}\frac{\gamma+\Gamma_{2}}{\gamma+\Gamma_{1}}}}{(\gamma+\Gamma_{2})^{2}+\Omega^{2}\frac{\gamma+\Gamma_{2}}{\gamma+\Gamma_{1}}+(\omega-\varepsilon/\hbar)^{2}}$
(S17)

where again the second fraction is in the form of a Lorentzian with width $w$
given as

$w=\sqrt{(\gamma+\Gamma_{2})^{2}+\Omega^{2}\frac{\gamma+\Gamma_{2}}{\gamma+\Gamma_{1}}}$
(S18)

This width is guaranteed to be larger than $\nu$ for large $\gamma$.
Therefore, we could still proceed with integrating the Lorentzian (the second
fraction) in Eq.(S17) over $\varepsilon$, which leads to a constant
independent of temperature and power. After another integration over the
probability distribution of $\gamma$, we can express the total internal loss
as

$\delta=m\delta_{0}\text{arcsinh}\left(\frac{\gamma}{\Omega}\right)\Biggr{|}_{\text{max}(\Omega,\sqrt{\Gamma_{2}^{2}+\nu^{2}})}^{\gamma_{\text{max}}}$
(S19)

The function $\text{arcsinh}(x)$ is reduced to $\ln(2x)$ for $x\gg 1$ giving
the observed logarithmic power dependence. And naturally for
$\Omega\rightarrow 0$, $\delta\rightarrow
m\delta_{0}\ln(\gamma_{\text{max}}/\sqrt{\Gamma_{2}^{2}+\nu^{2}})$, which is a
constant limit at low power much like the TLS saturation in STM.

The rate $\sqrt{\Gamma_{2}^{2}+\nu^{2}}$ is calculated from the characteristic
rate of the TLS dynamics with $\xi(t)=0$ and characterizes the boundary
between fast and intermediate jump rates $\gamma$.

###  S3.3 Intermediate $\gamma$ fluctuators

When $\gamma$ is comparable to $\Omega$ and/or $\Gamma_{1,2}$, several jumps
could be expected before the TLS relaxes to its ground state. Therefore, the
interaction among the fluctuators need to be accounted for, which complicates
the modelling compared to the low and high $\gamma$ fluctuators. The dynamical
effect is still described by a similar master equation Eq.(S16) for high
$\gamma$ fluctuators with $\rho$ replaced by $\rho_{k}$ for each different
state $k$ of the fluctuators. Moreover, the fast random telegraph noise is not
averaged out, and the dynamics of the transverse component can be treated as
stationary in a rotating wave approximation ($d\langle S^{+}\rangle/dt=0$).
The resulting master equation to be solved is then [71]

$\displaystyle\frac{\partial\rho_{k}}{\partial
t}+\frac{d}{dz}[(\Gamma_{1}(m-z)-\Xi_{k}z)\rho_{k}]=\gamma_{kn}\rho_{n}$ (S20)
$\displaystyle
y_{k}=\frac{-\Omega\Gamma_{2}}{\Gamma_{2}^{2}+(\varepsilon_{k}/\hbar-\omega)^{2}}z_{k}$
(S21)

where $\langle S_{z}^{0}\rangle$, $\langle S_{y}^{1}\rangle$ are abbreviated
as $z$ and $y$, and
$\Xi_{k}=\Omega^{2}\Gamma_{2}/[(\varepsilon_{k}-\omega)^{2}+\Gamma_{2}^{2}]$,
and $\gamma_{kn}$ is the transition rate between the fluctuator state $k$ and
$n$. $\Gamma_{1}$ and $m$ are not subscripted with $k$ since their
$\varepsilon_{k}$ dependence contains a slowly varying function in the form
$\tanh[\varepsilon_{k}/(2k_{B}T)]$ that should not vary much near the
resonance. $z_{k}$ can be solved from Eq.(S20) and substituted in Eq.(S21) to
obtain $y_{k}$

$\displaystyle
y_{k}=\frac{-\Omega\Gamma_{2}}{\Gamma_{2}^{2}+(\varepsilon_{k}-\omega)^{2}}\frac{n_{k}}{1+\frac{1+n_{k}\gamma/\Gamma_{1}}{\gamma+\Gamma_{1}}\Xi_{k}}$
(S22)

where $n_{k}$ is the probability for a TLS to have energy $\varepsilon_{k}$,
$\gamma=\sum_{n}(\gamma_{kn}+\gamma_{nk})$ is the totoal effective jump rate
for fluctuator state $k$. The single TLS y-component response with the energy
level $\varepsilon_{k}$ can then be calculated from integrating Eq.(S22) over
the probability distribution $P(\gamma)\propto 1/\gamma$

$\displaystyle\overline{y_{k}}\propto\frac{-m\Omega\Gamma_{2}}{\Gamma_{2}^{2}+(\varepsilon_{k}-\omega)^{2}}\left(\frac{1}{1+\kappa_{k}}\ln{\frac{\gamma_{h}}{\gamma_{l}}}+\frac{\kappa_{k}(1-n_{k})}{(1+\kappa_{k})(1+n_{k}\kappa_{k})}\ln{\frac{1+\kappa_{k}+\gamma_{h}/\Gamma_{1}(1+n_{k}\kappa_{k})}{1+\kappa_{k}+\gamma_{l}/\Gamma_{1}(1+n_{k}\kappa_{k})}}\right)$
(S23)

where $\kappa_{k}=\Xi_{k}/\Gamma_{1}$ is reduced to
$\kappa=\Omega^{2}/(\Gamma_{1}\Gamma_{2})$ when the TLS is on resonance
($\varepsilon_{k}=\hbar\omega$). The upper and lower bounds for the
intermediate fluctuators are obtained from the TLS dynamics as in Sec.S3.2
such that
$\gamma_{h}\gtrsim\Xi_{k}+\Gamma_{1},\sqrt{\Gamma_{2}^{2}+(\varepsilon_{k}/\hbar-\omega)^{2}}\gtrsim\gamma_{l}$.

The total loss can then be obtained from summing up the contribution from all
the different TLS energy levels just as in Eq.(15).

$\displaystyle\delta=\frac{\delta_{0}}{\pi}\sum_{k}\frac{m\Gamma_{2}}{\Gamma_{2}^{2}+(\varepsilon_{k}-\omega)^{2}}\left(\frac{1}{1+\kappa_{k}}\ln{\frac{\gamma_{h}}{\gamma_{l}}}+\frac{\kappa_{k}(1-n_{k})}{(1+\kappa_{k})(1+n_{k}\kappa_{k})}\ln{\frac{1+\kappa_{k}+\gamma_{h}/\Gamma_{1}(1+n_{k}\kappa_{k})}{1+\kappa_{k}+\gamma_{l}/\Gamma_{1}(1+n_{k}\kappa_{k})}}\right)$
(S24)

If one assumes that the energy drift induced by the individual fluctuator
follows the same Lorentzian distribution $P(\xi)$ for the weak fluctuators in
Sec.S3.1, the total drift will follow a Lorentzian distribution with a wider
width $P(\xi)\propto N\Gamma_{f}/((N\Gamma_{f})^{2}+\xi^{2})$ since there are
$N$ fluctuators affecting the TLS in the case of intermediate fluctuators. If
we assume that this width $N\Gamma_{f}$ is comparable to the detuning $\nu$ ,
then the broadening of the TLS energy levels will ensure the existence of
resonant TLS. Therefore, one can apply the continuous distribution of TLS
$P(\varepsilon,\Delta_{0})$. The final result of the internal loss after the
integration is then

$\displaystyle\delta\propto\tanh(\frac{\hbar\omega}{2k_{B}T})\left(\frac{1}{1+\kappa}\ln{\frac{\gamma_{h}}{\gamma_{l}}}+\frac{\kappa(1-n)}{(1+\kappa)(1+n\kappa)}\ln{\frac{1+\kappa+\gamma_{h}/\Gamma_{1}(1+n\kappa)}{1+\kappa+\gamma_{l}/\Gamma_{1}(1+n\kappa)}}\right)$
(S25)

The loss takes on a logarithmic power dependence in the intermediate power
range
$\left(\sqrt{\Gamma_{1}\Gamma_{2}}<\Omega<\sqrt{\Gamma_{2}^{2}+\nu^{2}}\right)$,
$\delta\sim m\ln((\Omega^{2}+\Gamma_{2}^{2})/(2\Omega^{2}))$. This loss
saturates to a constant $\sim m\ln(\sqrt{\Gamma_{2}^{2}+\nu^{2}}/\Gamma_{1})$
for lower powers, and to a small constant $\sim m\ln(2)$ for higher powers.

In summary, the power dependence of the loss can be separated into three
different regimes where in the high power limit $\Omega\gtrsim\Gamma_{2}$, the
effect of fluctuators that induce large $\gamma$ dominates and leads to a
logarithmic power dependence; in the intermediate power regime
$\Gamma_{2}\gtrsim\Omega\gtrsim\sqrt{\Gamma_{1}\Gamma_{2}}$, the fluctuators
with intermediate $\gamma$ give rise to a faster logarithmic power dependence,
but meanwhile the saturation of TLS just as in STM has a comparable if not
stronger power dependence and overlap in the same power regime; and finally in
the low power limit $\Omega<\sqrt{\Gamma_{1}\Gamma_{2}}$, the typical TLS
saturation in STM is recovered as the fluctuator contributions from all three
different regimes become constant in power. The above description matches our
experimental observation in Fig. 4.

The final fit to the data can be obtained from summing up all three
contributions from different regimes of fluctuators: the discrete summation
for the low and intermediate $\gamma$ fluctuators, Eq.(15) and Eq.(S24)
respectively, and the logarithmic power dependent loss from the high $\gamma$
fluctuators Eq.(S19), and is shown in Fig. 5.

##  S4 Alternative models for fitting TLS loss

Since none of the existing TLS model showed a strong temperature dependence of
loss below TLS saturation, the alternatives need to be composed of the
discrete TLS formalism to account for the temperature dependence and some
other TLS models that produce slower power dependence.

###  S4.1 Two TLS model

The two TLS or , by extension, multiple TLS model is a popular alternative in
explaining the slower power dependence than STM. The multiple species of TLS
could be attributed to nonuniform field distribution in the resonator [64],
different dipole moments of the TLS [69], and different microscopic origins of
the TLS [67]. The exact model for TLS loss is simply a sum of two different
STM loss contributions:

$\displaystyle\delta_{TLS}=\sum_{n}\delta_{0}^{TLS,1}\frac{m_{n}\Gamma_{2}}{\Gamma_{2}^{2}(1+\kappa)+(\varepsilon_{n}/\hbar-\omega)^{2}}+\delta_{0}^{TLS,2}\frac{m_{n}\Gamma_{2}}{\Gamma_{2}^{2}(1+\eta\kappa)+(\varepsilon_{n}/\hbar-\omega)^{2}}$
(S26)

where $m_{n}=\tanh(\varepsilon_{n}/(2k_{B}T))$ and its subscript explicitly
denotes its dependence on $\varepsilon_{n}$. The difference between the two
species of the TLS is encoded in their intrinsic losses
$\delta_{0}^{TLS,1},\delta_{0}^{TLS,2}$ which accounts for different
participation ratios, and the factor $\eta$ which could be due to different
field strength experienced by the TLS dipoles, different dipole moments, or
different relaxation rates of the two species of TLS. The resulting fit in
Fig.S3 shows a systematic deviation between the data and the fit, with
$\text{RMSE}=0.021$.

Fig. S3: The discrete two TLS model fit to the loss data.

There are 9 fitting parameters in this model. The resulting relaxation rates
from the fit are $\Gamma_{2}=4.7\times 10^{5}\text{ to }2\times 10^{7}Hz$ and
$\Gamma_{1}=6.8\times 10^{3}\text{ to }10^{4}Hz$. The losses from the
different contributions are determined as $\delta_{other}=1.29\times 10^{-5}$,
$\delta^{TLS,1}_{0}=6.9\times 10^{-6}$, and $\delta^{TLS,2}_{0}=6.8\times
10^{-6}$. The other fitting parameters are volume of the TLS-inhabiting
dielectrics $2\mu m^{3}$, the asymmetry ratio $\eta=18$, and the exponent of
the power dependent density of states of TLS $\mu=0.2$. We should note that
this fit might be improved by including more species of TLS. However, this
multiple-TLS model lacks a physical motivation and is not considered in this
work.

###  S4.2 Power law fit

Another well-known phenomenological fitting method introduces a free power law
dependence to the saturation power term $1+\kappa$ [27]. The model when
combined with the discrete formalism becomes

$\displaystyle\delta_{TLS}=\delta_{0}\sum_{n}\frac{m_{n}\Gamma_{2}}{\Gamma_{2}^{2}(1+\kappa)^{\beta}+(\varepsilon_{n}/\hbar-\omega)^{2}}$
(S27)

where $\beta$ is the exponent for the free power law dependence, and $\beta=1$
corresponds to STM prediction.

Fig. S4: The discrete power law fit to the loss data

The discrete free power law has a large deviation from the data, with
$\text{RMSE}=0.036$. The 8 fitting parameters are the power law exponent
$\beta=0.38$, the exponent in the power dependent density of states of TLS
$\mu=0.2$, the volume of TLS-inhabiting dielectrics $2.6\mu m^{3}$, the loss
from other mechanisms $\delta_{other}=1.29\times 10^{-5}$, the intrinsic TLS
loss $\delta_{0}=8\times 10^{-6}$, and the two relaxation rates are
$\Gamma_{1}=700\text{ to }10^{3}Hz$ and $\Gamma_{2}=3.3\times 10^{4}\text{ to
}1.4\times 10^{6}Hz$.

Moreover, neither of the alternative models provide a clear interpretation of
the physics in the TLS loss, since the phenomenological fitting parameters
could have different microscopic origins, i.e. nonuniform field distribution,
or different dipole moments or relaxation rates among several distinct
ensembles of TLS. Not only does the discrete form of GTM in this work generate
better fits to the measured loss, but it also provides a clear physical
explanation for the observed loss: the low spectral density of coherent TLS
and their interactions with the fluctuators are responsible for the observed
temperature and power dependence.

##  S5 Non-Equilibrium Quasiparticle Treatment and $n_{qp}$ and $Q_{qp}$
Estimates

The aim of this section is to examine the non-equilibrium quasiparticle
density by simulating the physics using the non-equilibrium energy-dependent
distribution of quasiparticles. This model (scattering model) was established
by considering the coupled quasiparticle and phonon systems.[6] The details of
this model are described in references [7] [8] and [9] and employ numerical
methods to discretize the distribution of quasiparticles $f(E)$ at energy $E$
and the phonon distribution $n(\Omega)$ at energy $\Omega$ by solving the
kinetic equations in steady state for a given flux of microwave photons and
higher frequency radiation.

Fig. S5:  The calculated quasiparticle distribution $f(E)$ (a) and phonon
distribution $n(\Omega)$ (b) as a function of normalized energy for different
circulating numbers of photons in the half wavelength resonator. Note that in
the plot of $f(E)$, a double y axis is used due to the different scales of
thermal and non-equilibrium distributions.

Fig. S5 (a) and Fig. S5 (b) show the calculated quasiparticle distribution and
phonon distribution, respectively, for three selected circulated numbers of
photons in the half wavelength resonator. This calculation is performed by
assuming the fridge base temperature $T_{b}=10\>mK$, the resonator drive
frequency 3.6442 GHz ($\hbar\omega=23\>\mu eV$), superconducting energy gap
$\Delta_{s0}=188$ $\mu eV$, $Q_{i}=10^{5}$, $Q_{c}=1.5\times 10^{6}$ which we
obtain from fitting the resonance, resonator center conductor volume
$8.6\times 10^{-14}$ $m^{3}$ and effective temperature $T_{eff}=189\>mK$ due
to stray light illumination and radiation which creates enhanced phonon
generation, as described by the Parker model [92]. Table SI lists all of the
parameter values used in this model. The black dashed line in Fig. S5(a)
indicates the thermal distribution of quasiparticles without any microwave
excitation at $T_{b}=10$ $mK$. At low circulating photon numbers ($<n>=1$ or
325 ), the quasiparticle distribution is enhanced significantly above the
$T_{b}$ thermal distribution, and the phonon distribution is also enhanced. At
high circulating photon numbers, jumps appear in the electron and phonon
distributions every $\hbar\omega$ because microwave drive at high power
significantly affects the distributions. Big step jumps at $E=3\Delta_{s0}$
and $\Omega=2\Delta_{s0}$ are due to pair breaking and recombination
processes.

From the quasiparticle distribution $f(E)$, one can calculate the
quasiparticle density, $n_{qp}$, by numerical integration over all energy,
with its density of states $\rho(E)$

$n_{qp}=4N_{0}\int_{\Delta_{s0}}^{\infty}f(E)\rho(E)dE\quad\text{with}\quad\rho(E)=\frac{E}{\sqrt{E^{2}-\Delta_{s0}^{2}}}$
(S28)

The calculated $n_{qp}$ at different circulated photon numbers $\langle
n\rangle$ are shown in Fig. S6(a). Because this calculation is performed at
the 10 mK fridge base temperature, the entire term $n_{qp}$ (Eq. (16) of the
main text) can be regarded as the contribution of non-equilibrium quasi-
particle density since the number of thermal quasiparticles is extremely
small. The calculated $n_{qp}$ is around $50\;\mu m^{-3}$ when $\langle
n\rangle\;<10^{6}$ . In addition, the calculated result of $n_{qp}$ as a
function of fridge/bath temperature ($T_{b}$) at $\langle n\rangle=1$ is shown
in Fig. S6(b) as circle dots. When the bath temperature is below 150 $mK$,
$n_{qp}$ remains constant. This is consistent to the assumptions made in the
frequency shift fit at low temperature where quasiparticle density is a
constant. The continuous line of Fig. S6(b) is from Eq. (16) of the main text.

For the quality factor due to the quasiparticles, $Q_{qp}$ is defined as

$\frac{1}{Q_{qp}}=\alpha\frac{\sigma_{1}}{\sigma_{2}}$ (S29)

where the kinetic inductance ratio $\alpha=0.0185$ was obtained from the
frequency shift fit in the main text. The real and imaginary parts of the
complex conductivity $\sigma=\sigma_{1}-j\sigma_{2}$ can be expressed
respectively by the Mattis-Bardeen formula [89] given by

$\displaystyle\frac{\sigma_{1}}{\sigma_{N}}(\omega)=\frac{2}{\hbar\omega}\int_{\Delta_{s0}}^{\infty}\left[f(E)-f(E+\hbar\omega)\right]g_{1}(E)dE$
(S30)
$\displaystyle\frac{\sigma_{2}}{\sigma_{N}}(\omega)=\frac{2}{\hbar\omega}\int_{\Delta_{s0}-\hbar\omega}^{\infty}\left[f(E)-2f(E+\hbar\omega)\right]g_{2}(E)dE$
(S31)

where $\sigma_{N}$ is the normal-state conductivity,
$g_{1}(E)=h_{1}(E,E+\hbar\omega)\rho(E)$ and
$g_{2}(E)=h_{1}(E,E+\hbar\omega)\frac{E}{\sqrt{\Delta_{s0}^{2}-E^{2}}}$ with
$h_{1}(E,E^{\prime})=(1+\frac{\Delta_{s0}^{2}}{E*E^{\prime}})\rho(E^{\prime}).$
Here we use the non-equilibrium distribution function $f(E)$ discussed above.

The calculated result of $Q_{qp}$ as a function of average circulated photon
numbers is shown in Fig. S7. In the regime of low circulating photon numbers,
$Q_{qp}$ remains a constant and then gradually increases as the circulating
photon numbers in the half wavelength resonator increase. Overall, the
$Q_{qp}$ is on the order of $10^{7}$, and this result verifies our assumption
that $Q_{qp}>Q_{TLS0}$ in the $Q_{i}(T)$ fitting. In other words, the loss in
the half wave length resonator at low temperatures and low circulating photon
numbers is dominated by the TLS losses.

Fig. S6:  (a) The calculated quasiparticle density as a function of circulating photon numbers in the half wavelength resonator. Here $T_{b}=10\>mK$ and $T_{eff}=189\>mK$. (b) The calculated quasiparticle density as a function of fridge temperature at $<n>=1$. Note the dots are calculated from the non-equilibrium model, and the solid line is from Eq. (16) of the main text.  Fig. S7:  The calculated quality factor of quasi-particle as a function of circulating photon numbers $<n>$ in the half wavelength resonator. Here $T_{b}=10\>mK$ and $T_{eff}=189\>mK$, and other parameters are listed in Table SI. Superconducting gap ($\Delta_{s0}$) = $188\>\mu eV$ |  RF photon energy ($\hbar\omega$) = $23\>\mu eV$
---|---
|
$Q_{i}=1.0\times 10^{5}$ |  $Q_{c}=1.5\times 10^{6}$
|
Quasiparticle-phonon time =438 ns [8] |  Characteristic phonon time = 0.26 ns [8]
|
Phonon escape time =0.17 ns [8] |  Resonator volume= $8.6\times 10^{-14}\;m^{3}$
|
Base Temperature $T_{b}$=10 mK |  Phonon effective temperature ($T_{eff}$)=189 mK
Table SI:  Parameters used for our non equilibrium quasi-particle calculation.

##  S6 CST Microwave Simulation of CPW Electric Fields

A model of the CPW microwave resonator was constructed in CST Microwave
Studio. The model structure (Fig. S8) represents the entire CPW resonator and
coupling capacitors, both of which reproduce the geometrical structure in our
superconducting chip. The superconductor is modeled as a perfect electric
conductor for the purpose of E field calculation. The results below are
obtained from the finite element frequency domain solver in CST.

Fig. S8: CPW resonator model in CST microwave studio. (a) The geometry of the
CPW resonator model in CST. The coupling capacitors and the entire resonator
structure reproduce those in the experiment. (b) The top view of the E field
strength on the substrate-vacuum interface at the fundamental resonance of
3.647 GHz, which has a node in the center of the resonator, a typical standing
wave pattern in a half wavelength resonator. (c) The close-up top view of the
E field vector plot on the substrate-vacuum interface near the coupling
capacitor. The maximum E field is found around the corner of the center strip
that is part of the coupling capacitor.The E field that contributes to the TLS
model is estimated from the average of the field along the resonator side of
the coupler about $1.16\times 10^{8}Vm^{-1}$, highlighted in the yellow-orange
line. (d) The side view of the E field vector plot on the cross section
through the coupler at the center strip corner. The maximum is at the corner
of the center strip on the substrate-vacuum interface.

The CST simulation shows that on resonance at 3.647 GHz, the electric field
can attain a maximum of $6\times 10^{8}$ V/m on the substrate-vacuum interface
at the corner of the center strip that is part of the coupling capacitor on
the input side. This calculation was done under a 0.5 watt excitation level.
An average electric field of the adjacent area is estimated to be $1.16\times
10^{8}$ $V/m$. By scaling this power down to that required to achieve one
circulating photon in the resonator, we estimate the average electric field of
the region nearby the coupling capacitor of the resonator to be $0.2$ V/m. In
Fig. 4 and 5, the Rabi frequency $\Omega$ on the x-axis is estimated from this
electric field and a dipole moment estimated from [81, 82].