# Ultrastrong Magnon-Photon Coupling Achieved by Magnetic Films in Contact
with Superconducting Resonators

Alberto Ghirri<EMAIL_ADDRESS>Istituto Nanoscienze - CNR, Centro
S3, via G. Campi 213/A, 41125, Modena, Italy  Claudio Bonizzoni  Dipartimento
di Scienze Fisiche, Informatiche e Matematiche Università di Modena e Reggio
Emilia, via G. Campi 213/A, 41125, Modena, Italy  Istituto Nanoscienze - CNR,
Centro S3, via G. Campi 213/A, 41125, Modena, Italy  Maksut Maksutoglu
Institute of Nanotechnology, Gebze Technical University, 41400, Gebze,
Kocaeli, Turkey  Alberto Mercurio  Dipartimento di Scienze Matematiche e
Informatiche, Scienze Fisiche e Scienze della Terra, Università di Messina,
I-98166 Messina, Italy  Omar Di Stefano  Dipartimento di Scienze Matematiche e
Informatiche, Scienze Fisiche e Scienze della Terra, Università di Messina,
I-98166 Messina, Italy  Salvatore Savasta  Dipartimento di Scienze Matematiche
e Informatiche, Scienze Fisiche e Scienze della Terra, Università di Messina,
I-98166 Messina, Italy  Marco Affronte  Dipartimento di Scienze Fisiche,
Informatiche e Matematiche Università di Modena e Reggio Emilia, via G. Campi
213/A, 41125, Modena, Italy  Istituto Nanoscienze - CNR, Centro S3, via G.
Campi 213/A, 41125, Modena, Italy

###### Abstract

Coherent coupling between spin wave excitations (magnons) and microwave
photons in a cavity may disclose new paths to unconventional phenomena as well
as for novel applications. Here, we present a systematic investigation on YIG
(Yttrium Iron Garnet) films on top of coplanar waveguide resonators made of
superconducting YBCO. We first show that spin wave excitations with frequency
higher than the Kittel mode can be excited by putting in direct contact a 5
$\mu$m thick YIG film with the YBCO coplanar resonator (cavity frequency
$\omega_{c}/2\pi=8.65$ GHz). With this configuration, we obtain very large
values of the collective coupling strength $\lambda/2\pi\approx 2$ GHz and
cooperativity $C=5\times 10^{4}$. Transmission spectra are analyzed by a
modified Hopfield model for which we provide an exact solution that allows us
to well reproduce spectra by introducing a limited number of free parameters.
It turns out that the coupling of the dominant magnon mode with photons
exceeds 0.2 times the cavity frequency, thus demonstrating the achievement of
the ultrastrong coupling regime with this architecture. Our analysis also
shows a vanishing contribution of the diamagnetic term which is a peculiarity
of pure spin systems.

##  I Introduction

The interplay between magnetic excitations and electromagnetic radiation has
recently assumed a pivotal role in many fields of research such as magnonics,
spintronics, magneto-opto-mechanics and information processing for its
potentialities in the development of hybrid systems and devices [1, 2, 3, 4].
The control of magnon-photon coupling and cooperativity is one of the keys for
enabling the exploitation of unique functionalities related to the coherent
dynamics in these systems. Novel applications, including memory devices [5],
coherent spin pumping [6], haloscopes for axion detection [7], microwave-
optical transducers [8] and coherent microwave sources [9] have been already
tested. An open issue is the realization of all-on-chip devices for their
efficient integration in microwave circuits [4].

In the prototypical case of a ferromagnetic sample embedded in a microwave
resonator, spin waves couple with resonant electromagnetic modes [10] and the
system can be modelled by combining Maxwell and Landau-Lifshitz-Gilbert (LLG)
equations [11, 12]. This classical description works remarkably well in the
case of the ferrimagnetic Yttrium Iron Garnet (YIG) [3], that has been studied
in detail showing a combination of several optimal features, including the
exceptionally low damping of magnetization precession [13].

Figure 1: Sketch of the YBCO/sapphire CPW resonator with the YIG/GGG film
positioned above.

In the quantum regime, spin waves are collective (bosonic) excitations
(magnons) that may coherently interact with cavity photons [14, 15]. The Dicke
model and its evolution in the Hopfield version, which also includes the
diamagnetic term, have been developed to describe a large variety of hybrid
light-matter quantum systems [16, FornDiazRM19]. In the strong coupling
regime, the rotating wave approximation (RWA) and the (Jaynes) Tavis-Cummings
Hamiltonian have been mainly used to interpret experimental data [3]. However,
there is a general trend now to push these studies beyond conventional
coupling regimes and more specifically to reach coupling strengths, $\lambda$,
being a non-negligible fraction of the cavity frequency, $\omega_{c}$. Under
these conditions novel effects related to processes that do not conserve the
number of excitations in the system have been predicted [16, 54]. For
$\lambda/\omega_{c}\geq$0.1, we refer to the ultrastrong coupling (USC)
regime, in which counter-rotating terms, neglected in RWA, must be considered.
In this context, the description of the fundamental interaction between the
electromagnetic field and the magnetic system, including the diamagnetic term
which plays a key role in the superradiant phase transition [17, 18, 19, 20],
still need to be clarified and tested on real magnetic materials.

Figure 2: Calculated profile of the oscillating field (magnetic component) for
the fundamental mode of the bare CPW resonator. (a) Top view of the
distribution $h_{ac}(x,y,z=0)$ showing the magnetic antinode in the middle of
the resonator. (b) Profile of $\mathbf{h_{ac}}$ (phase = 0) in the plane
perpendicular to the $x$-direction. The $x$ component of $\mathbf{h_{ac}}$ is
vanishingly small. (c) Mean value of $h_{ac}(z)$ plotted as a function of the
$z$ height. Figure 3: Broadband characterization of the YIG film. (a)
Ag/alumina microstrip, case $\\#$A. Top: field-swept transmission spectra
(temperature $T=50$ K), taken at different continuous wave frequencies, as
indicated. Bottom: plot of the excitation frequencies as a function of the
measured resonance field; the dashed line shows the fit with Eq. 1. (b)
YBCO/sapphire CPW line, case $\\#$B ($T=30$ K). Center: Spectral map showing
$\partial S_{21}(\omega,H_{0})/\partial{H_{0}}$. Top: Plot of $S_{21}$ and
$\partial S_{21}/\partial{H_{0}}$ as function of $H_{0}$ for
$\omega/2\pi=8.65$ GHz. Right: Plot of $S_{21}$ and $\partial
S_{21}/\partial{H_{0}}$ as function of $\omega$ for $H_{0}=0.21$ T. (c) Same
as (b) with YIG in contact with the CPW line, case $\\#$C.

Few experimental reports show magnetic systems achieving the USC regime with
microwave resonators. The USC regime was first reported in the case of mm-size
YIG crystals in 3D cavities [21, 22, 23, 24, 25] and few other magnetic
materials [26, 27, 28, 29]. In view of the realization of scalable
architectures, small magnets coupled to superconducting planar resonant
geometries have been recently designed and observed to achieve the
(ultra-)strong coupling regime [30, 31, 32, 33, 34, 35, 36, 37]. In spite of
these encouraging results, the optimization of these hybrid systems for the
USC is still largely unexplored. For instance, the geometry of the hybrid
superconductor-magnet system needs to be optimized in order to maximize their
mutual coupling, while magnetic materials may be preferably insulating in
order to minimize the damping. Moreover, the vicinity to a superconducting
material, which has to be resilient to magnetic fields [38] to allow for
ferromagnetic resonance (FMR) experiments, may significantly affect the
spectrum of magnetic excitations through several mechanisms [33, 39].

In this article, we address the problem of reaching the USC regime in coupled
superconductor/ferrimagnet hybrid architectures. By means of a series of
systematic experiments carried out with YIG films and planar transmission line
resonators, we demonstrate the achievement of high coupling rates by
positioning the magnetic film in direct contact with the superconducting
resonator. In our geometry, the excitation of spin waves take place at the
superconductor/ferrimagnet interface, where the amplitude of the microwave
field is maximum. The optimized magnon-photon coupling results in collective
coupling strengths as large as 0.2 times the cavity frequency. Data analysis,
carried out with a modified Hopfield model for which we provide an exact
solution, also evidences vanishingly small diamagnetic coupling for magnon
excitations in YIG.

The article is organised as follows: we start presenting the experimental
methods in Section II whist in Section III we report numerical simulations of
the profile of the oscillating field which allow us to estimate volume and
coupling strength involved in the spin-photon process. Experimental results
are presented in progressive way by considering the broadband characterization
first in Section IV, then transmission spectra on resonators in Section V.
Modelling and data analysis are reported in Section VI, then discussion and
concluding remarks in the latest sections.

##  II Experimental methods

We prepared planar resonators with different sizes and geometries, including
meanders with multiple resonance frequencies and inverse anapole resonators
with focused electromagnetic radiation [40], to study the coupling with
different YIG samples. Here we report on two types of planar devices:
microstrip and coplanar waveguide (CPW) resonator (Fig. 1). For each of them,
the corresponding broadband transmission line, having the same dimensions and
geometry (except for the presence of the input and output coupling gaps
defining the half-wavelength resonator), has been fabricated to investigate
the spin wave excitation spectrum in the different cases. Microstrip lines
having $500~{}\mu$m wide central strips were obtained by wet etching of
$\mathrm{Ag/Al_{2}O_{3}}$ films [41]. Superconducting coplanar waveguide (CPW)
lines (Fig. 1) were fabricated from commercial YBa2Cu3O7 (YBCO) films
deposited on sapphire. Etching was carried out by Ar plasma in a reactive ion
etching (RIE) chamber. The central conductor had length $l\approx 6$ mm, width
$w=(17\pm 1)~{}\mu$m and separation $s=(14\pm 1)~{}\mu$m from the lateral
ground planes. Further details are given in [42].

The experiments were carried out at low temperature in applied magnetic field
($\mathbf{H_{0}}$), which was oriented in the plane of the film along $x$
(Fig. 1). Frequency-swept spectra were acquired using a Vector Network
Analyzer with typical incident power $P_{inc}=-8$ dBm. We studied YIG films
grown by liquid-phase epitaxy on $500~{}\mu$m thick gadolinium gallium garnet
(GGG) substrates. Unless specified, the film has rectangular shape with
thickness of $5~{}\mu$m and in-plane sizes of $4\times 3~{}\mathrm{mm}^{2}$.
The sample was mounted in diverse ways on the aforementioned microwave
transmission lines and resonators. We anticipate that different spin-wave
spectra and, ultimately, different magnon-photon coupling strengths were
obtained in different cases. In the following we thus consider three
experimental configurations: YIG film glued on metallic microstrip (case
$\\#$A); YIG film glued on YBCO CPW line (case $\\#$B); YIG film in contact
with the YBCO CPW line (case $\\#$C). In $\\#$A and $\\#$B the estimated film-
resonator separation was on the order of $10~{}\mu$m. In case $\\#$C we
obtained a good contact between YBCO and YIG surfaces thanks to a
polytetrafluoroethylene (PTFE) screw that gently pushed the GGG substrate
against the YBCO surface from the backside.

##  III Electromagnetic simulations

Finite-element electromagnetic simulations were carried out using commercial
software to evaluate amplitude and distribution of the resonator fields. The
maps calculated for the CPW resonator (Fig. 2(a,b)) show the expected in-plane
and out-of-plane profiles of the microwave magnetic field, $\mathbf{h_{ac}}$.
The resonator field is mostly confined around the central line between the
lateral ground planes, in a region of approximate width $w+2s=45~{}\mathrm{\mu
m}$ [42]. To evaluate the $z$ dependence of the resonator field, we calculated
the in-plane averaged value of $h_{ac}$ as a function of the $z$ height (Fig.
2(c)). The decay of $h_{ac}(z)$ is quasi-exponential and for
$z=10~{}\mu\mathrm{m}$ it results in a reduction of a factor $\approx 2$ with
respect to the maximum value.

We can exploit these electromagnetic simulations to estimate the amplitude of
the vacuum fluctuation $b_{vac}$ and, consequently, of the spin-photon
coupling strength ($g_{s}$) expected with this particular CPW resonator. By
rescaling the calculated magnetic field to the single photon power level [42],
we obtain $b_{vac}\approx 3$ nT. This value is in perfect agreement with the
value derived from $b_{vac}\approx\mu_{0}\omega_{c}/(4w)\sqrt{h/Z_{0}}$ [43,
30], being $h=6.626\times 10^{-34}$ J s the Planck constant and
$Z_{0}=58~{}\mathrm{\Omega}$ the nominal impedance of the CPW resonator. The
spin-photon coupling strength then results $g_{s}=\gamma b_{vac}/4=21$ Hz,
where $\gamma$=28.02 GHz/T is the electron’s gyromagnetic ratio. For
comparison, the value of the spin-photon coupling expected for the fundamental
mode of the microstrip resonator is less than 1 Hz.

Figure 4: Coupling of the YIG film with different planar resonators
($P_{inc}=-8$ dBm). (a) Microstrip resonator, case $\\#$A ($T=50$ K). (b) CPW
resonator, case $\\#$B ($T=30$ K). (c) YIG in contact with the YBCO CPW
resonator, case $\\#$C ($T=30$ K). The blue dashed lines display
$\omega_{FMR}(H_{0})$ (Eq. 1) while the white dash-dot lines indicate the
resonators’ frequencies. The dotted lines correspond to the frequency of a
broad mode due to the metal box hosting the CPW lines.

##  IV Broadband transmission spectroscopy characterization

We first present a series of transmission ($S_{21}$) spectra acquired with
broadband transmission lines by sweeping the microwave frequency ($\omega$) at
steady values of the external magnetic field $H_{0}$ (Fig 3). In the case of
the microstrip ($\\#$A), the main resonance dip follows the Kittel relation
(Fig. 3(a)), which, for in-plane magnetized films, reduces to [44]:

$\displaystyle\omega_{FMR}=\sqrt{\omega_{H}(\omega_{H}+\omega_{M})},$ (1)

where $\omega_{H}/2\pi=\gamma\mu_{0}H_{0}$ and
$\omega_{M}/2\pi=\gamma\mu_{0}M_{s}$, being $\mu_{0}=4\pi\times 10^{-7}$ H/m
the vacuum permeability and $\mu_{0}M_{s}=0.245$ T the saturation
magnetization of YIG as reported in the literature [13, 42].

For case $\\#$B, the numerical derivative of the signal $\partial
S_{21}(\omega,H_{0})/\partial{H_{0}}$ is shown in Fig. 3(b). Additional spin
wave resonance modes are well visible in this case. Such modes are commonly
observed for YIG films [45, 46] but their frequency specifically depends by
the profile of the microwave field $\mathbf{h_{ac}}$. In particular, narrow
CPW lines efficiently excite travelling spin waves with finite wavevector
$0<k\leq 2\pi/s$ [47]. Given the in-plane magnetization of the film
($\mathbf{H_{0}}||x$) and the negligible $x$-component of $\mathbf{h_{ac}}$
(Fig. 2), the dispersion of the Damon-Eshbach modes follows the characteristic
expression [12]:

$\displaystyle\omega_{DE}=\sqrt{\left(\omega_{H}+\frac{\omega_{M}}{2}\right)^{2}-\left(\frac{\omega_{M}}{2}\right)^{2}e^{-2kd}},$
(2)

where $k=k_{s}=2\pi/s=4.5\times 10^{5}~{}\mathrm{rad~{}m^{-1}}$ (dash-dot line
in Fig. 3(b)) [48]. We note that $\omega_{DE}$ is near the maximum limit of
Eq. 2, $\omega_{H}(H_{0})+\omega_{M}/2$, whilst the coupling between the
magnetization of YIG and GGG [49] is not evident from the measured spectra
[42].

Fig. 3(c) displays the spectrum obtained in case $\\#$C. Here the transmission
spectrum shows several additional spin wave modes whose origin can be related
to the special boundary conditions imposed by the contact of the magnetic YIG
with the superconducting YBCO planes, as observed in similar experiments [50,
51]. Such modes are located in a wide absorption band between
$\omega_{DE}(H_{0})$ and $\omega_{H}(H_{0})+\omega_{M}$ and their typical
linewidth is $10\lesssim\kappa_{m}\lesssim 40$ MHz.

##  V Enhanced Spin Photon Coupling and Cooperativity

Once we had identified the main features of the spin wave excitation spectrum,
we performed experiments in a similar configuration, with the YIG film on
half-wavelength resonators having dimensions and geometries equivalent to the
ones of the broadband lines used in Fig. 3 [42]. The $S_{21}$ map taken with
the microstrip resonator (case $\\#$A) is shown in (Fig. 4(a)). The lower
branch has an oblique asymptote with the FMR mode (Eq. 1), while the upper
branch is barely visible owing to multiple resonances with spin wave modes.
This situation is common with other spin ensemble [52, 53] or YIG [3] systems
in a resonator mode. From Fig. 4(a) the collective coupling strength can be
roughly estimated to be of order of hundreds of MHz.

The YBCO CPW resonator has, with respect to the microstrip, a smaller mode
volume that results in a larger spin-photon coupling $g_{s}$ (Sect. III). A
large anticrossing, having well defined polariton branches, is obtained in
this case ($\\#$B, Fig. 4(b)). Even larger couplings were obtained when the
YIG film was put in direct contact with the YBCO resonator (case $\\#$C): with
this geometry the splitting raises to $2\lambda/2\pi\approx 4$ GHz (Fig.
4(c)), a factor $\approx 2$ larger than case $\\#$B. Within the anticrossing
gap, no additional modes are visible except for $\omega_{FMR}$. This suggests
that this mode is weakly coupled to the CPW resonator. Additional spin wave
modes are observed above $\approx 10$ GHz due to the presence of a spurious
box mode (Fig. 4(b,c)).

It is worth to estimate the cooperativity, that is commonly defined as
$C=4\lambda^{2}/(\kappa_{m}\kappa_{c})$ [54]. From our spectra we take
$\kappa_{m}/2\pi=40$ MHz and $\kappa_{c}/2\pi=8$ MHz as decay rates of magnon
(Sect. IV) and uncoupled cavity [42] modes. With these numbers, we achieve
$C=5\times 10^{4}$. Considering the normalized parameter
$U=\sqrt{(4\lambda^{2}/\kappa_{m}\kappa_{c})(\lambda/\omega_{c})}=108$, we
infer that our results ranks among the best performing physical platforms for
the reaching of the USC regime [54].

##  VI System Hamiltonian

We model our system by considering a quantized single-mode electromagnetic
field (with $\omega_{c}$ cavity frequency) interacting with an ensemble of
magnetic moments. We consider collective operators for the spin ensemble, the
quantization of both spin excitations and the electromagnetic field which
allows us to introduce the respective bosonic operators $\hat{a}$ and
$\hat{b}$. Due to the vanishing orbital angular momentum of Fe3+ in YIG [12],
we expect a prominent Zeeman interaction of the type $\hat{\mathcal{H}}_{\rm
Z}=-g_{e}\hat{\mathbf{\sigma}}\cdot\mu_{B}\hat{\mathbf{h}}$ for a single spin.
Here
$\commutator{\hat{\sigma}_{j}}{\hat{\sigma}_{k}}=i\epsilon_{jkl}\hat{\sigma}_{l}$
are the Pauli operators, $\hat{\mathbf{h}}$ is the magnetic field component of
the cavity resonator while $\mu_{B}$ is the Bohr magneton and $g_{e}\approx 2$
in the case of a simple electron. However, to not exclude the possibility of
having orbital angular momentum contributions in our hybrid system, we also
consider this degree of freedom including a diamagnetic term, which comes from
the usual minimal coupling replacement.

The total Hamiltonian then reads ($\hbar=1$)

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$
$\displaystyle\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\omega_{b}}{2}\sum_{j=1}^{N}\sigma_{z}^{(j)}+\frac{\lambda}{2\sqrt{N}}\sum_{j=1}^{N}\sigma_{x}^{(j)}\left(\hat{a}+\hat{a}^{\dagger}\right)$
(3) $\displaystyle+\beta\left(\hat{a}+\hat{a}^{\dagger}\right)^{2}~{},$

where $\hat{a}$ ($\hat{a}^{\dagger}$) is the photon annihilation (creation)
operator, $\omega_{c}$ is the cavity resonance frequency, $\omega_{b}$ is the
resonance frequency of a single spin, $\lambda$ is the collective light-matter
coupling, and $\beta$ is the coefficient of the diamagnetic term.

By using the collective spin operators
$\hat{J}_{z}\equiv(1/2)\sum_{j=1}^{N}\hat{\sigma}_{z}^{(j)}$ and
$\hat{J}_{x}=\hat{J}_{+}+\hat{J}_{-}\equiv(1/2)\sum_{j=1}^{N}\hat{\sigma}_{x}^{(j)}$,
we can apply the Holstein-Primakoff transformations [55]

$\hat{J}_{z}\to\hat{b}^{\dagger}\hat{b}-\frac{N}{2}~{},\quad\hat{J}_{+}\to\hat{b}^{\dagger}\sqrt{N-\hat{b}^{\dagger}\hat{b}}~{},\quad\hat{J}_{-}=\hat{J}_{+}^{\dagger}~{},$
(4)

where $\hat{b}$ and $\hat{b}^{\dagger}$ are the magnon annihilation and
creation operators respectively, which obey to the standard bosonic
commutation relations. In the thermodynamic limit (i.e. $N\to\infty$) we can
approximate $\hat{J}_{+}\approx\sqrt{N}\hat{b}^{\dagger}$. Then, by applying
the Holstein-Primakoff transformation, we obtain:

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$
$\displaystyle\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{b}\hat{b}^{\dagger}\hat{b}+\lambda\left(\hat{b}+\hat{b}^{\dagger}\right)\left(\hat{a}+\hat{a}^{\dagger}\right)$
(5) $\displaystyle+\beta\left(\hat{a}+\hat{a}^{\dagger}\right)^{2}~{},$

that is the well-known Hopfield Hamiltonian [56].

We consider the dependence of the magnon resonance frequency $\omega_{b}$ to
the external magnetic field $H_{0}$ as described by
$\omega_{b}=\sqrt{\omega_{H}(\omega_{H}+\omega_{M})}+\Delta$, leaving as the
sole free parameter the energy shift $\Delta$ characterizing high frequency
magnons for the next step of our investigation. In our analysis we also leave
as free parameters the cavity frequency $\omega_{c}$, the collective coupling
$\lambda$ and the factor of the diamagnetic term $\beta$.

The Hamiltonian in Eq. 3 of the main text can be expressed in terms of two
non-interacting harmonic oscillators
$\hat{\mathcal{H}}=\Omega_{-}\hat{P}_{-}^{\dagger}\hat{P}_{-}+\Omega_{+}\hat{P}_{+}^{\dagger}\hat{P}_{+}$,
where $\hat{P}_{\pm}$ are the polariton operators, which are linear
combinations of light and matter operators
$\hat{P}_{\mu}=c_{1}^{(\mu)}\hat{a}+c_{2}^{(\mu)}\hat{b}+c_{3}^{(\mu)}\hat{a}^{\dagger}+c_{4}^{(\mu)}\hat{b}^{\dagger}$,
with $\mu=\pm$. To fit our parameters, we need first to find the polariton
frequencies, and, being a proper bosonic excitation of the system, the
operator $\hat{P}_{\mu}$ fulfills the equation of motion of the harmonic
oscillator $[\hat{P}_{\mu},\hat{\mathcal{H}}]=\Omega_{\mu}\hat{P}_{\mu}$.
Since the polariton operator $\hat{P}_{\mu}$ is a linear combination of the
light and matter operators, we need to calculate first the commutator of the
latter with the Hamiltonian

$\displaystyle\commutator{\hat{a}}{\hat{\mathcal{H}}}$ $\displaystyle=$
$\displaystyle\omega_{c}\hat{a}+\lambda\left(\hat{b}+\hat{b}^{\dagger}\right)$
$\displaystyle\commutator{\hat{b}}{\hat{\mathcal{H}}}$ $\displaystyle=$
$\displaystyle\omega_{b}\hat{b}+\lambda\left(\hat{a}+\hat{a}^{\dagger}\right)$
$\displaystyle\commutator{\hat{a}^{\dagger}}{\hat{\mathcal{H}}}$
$\displaystyle=$
$\displaystyle-\omega_{c}\hat{a}^{\dagger}-\lambda\left(\hat{b}+\hat{b}^{\dagger}\right)$
$\displaystyle\commutator{\hat{b}^{\dagger}}{\hat{\mathcal{H}}}$
$\displaystyle=$
$\displaystyle-\omega_{b}\hat{b}^{\dagger}-\lambda\left(\hat{a}+\hat{a}^{\dagger}\right)~{},$

and the polariton frequencies $\Omega_{\mu}$ are obtained by finding the
positive eigenvalues of the following Hopfield matrix

$\mathcal{M}=\begin{pmatrix}\omega_{c}+2\beta&\lambda&-2\beta&-\lambda\\\
\lambda&\omega_{b}&-\lambda&0\\\ 2\beta&\lambda&-\omega_{c}-2\beta&-\lambda\\\
\lambda&0&-\lambda&-\omega_{b}\end{pmatrix}~{}.$ (6)

leading to:

$\Omega_{\pm}=\frac{1}{\sqrt{2}}\sqrt{\tilde{\omega}_{c}^{2}+\omega_{b}^{2}\pm\sqrt{\left(\tilde{\omega}_{c}^{2}-\omega_{b}^{2}\right)^{2}+16\omega_{c}\omega_{b}\lambda^{2}}}~{},$
(7)

where $\tilde{\omega}_{c}=\sqrt{\omega_{c}(\omega_{c}+4\beta)}$.

The above equation fits well the peaks of the $S_{21}$ spectrum (case $\\#$C).
Figure 5 shows the best fit result obtained with these parameters:
$\omega_{c}/2\pi=8.65$ GHz, $\Delta/2\pi=2.05$ GHz, and $\lambda/2\pi=2.002$
GHz. For what concerns the diamagnetic parameter $\beta$, the only nontrivial
result (i.e. nonzero result from fit) has been obtained by assuming a
dependence on the magnon frequency: $\beta=\alpha/\sqrt{\omega_{b}}$, which
can be justified by the fact that the presence of this term is dominant at low
frequencies. With this assumption, we obtained $\alpha/\sqrt{2\pi}=3\times
10^{-3}\ \text{GHz}^{\frac{3}{2}}$, corresponding to a value of the
diamagnetic coefficient $\beta/2\pi\sim 10^{-3}$ GHz on resonance condition.
The ratio between the collective coupling and the cavity frequency is
approximately 0.23, fulfilling the criterion $\lambda/\omega_{c}>0.1$ for USC.
The fit confirms that the influence of the diamagnetic term is almost
negligible leading us to conclude that the system couples to the resonator
mainly through the spins. Notice that its value is more than two orders of
magnitude smaller than the standard diamagnetic term for electric dipolar
interactions (Appendix D).

##  VII Discussion

The dispersion characteristics of the spin wave spectrum in an infinite
ferromagnetic film have been reported by Kalinikos and Slavin [57]. Taking
into account both dipole-dipole and exchange interactions, theory predicts the
excitation of perpendicular standing spin waves (PSSW) modes even in the long
wavenumber ($k_{y}d\ll 1$) limit. These modes are due to the broken
translational invariance along the film thickness. Remarkably, in these
conditions, the lowest mode shows a quasi-uniform profile with a dispersion
equation very similar to the Damon-Eshbach dipolar surface mode (Eq. 2).
Conversely, higher order PSSW modes display a nonuniform magnetization profile
along $z$ [57, 58].

Broadband spectroscopy data in Fig. 3 evidence that the spin wave spectrum is
influenced by two main factors: (i) the distribution of the exciting
electromagnetic field and (ii) the boundary conditions at the ferrimagnet-
superconductor interface. The effect of (i) emerges from the comparison
between Fig. 3(a) and (b). In the case of the wide microstrip line, the modes
calculated with Eq. 1 and Eq. 2, the latter by considering $k=2\pi/w^{\prime}$
with $w^{\prime}=500~{}\mu\mathrm{m}$, are very near. They both follow the
measured dispersion of lowest resonance mode in Fig. 3(a), showing that the
magnetostatic approximation is valid in this case. Conversely, in the case of
the narrow CPW line, the modes calculated with Eq. 1 and Eq. 2 ($k=2\pi/w$)
are different. Since the main absorption line is found for
$\omega(H_{0})\approx\omega_{DE}(H_{0})$ in Fig. 3(b), according to the
dipole-exchange theory [57, 58], we expect this mode shows a homogeneous
magnetization profile along the film thickness.

The effect of (ii) is clearly evident from the direct comparison between Fig.
3(b) and (c), which shows the appearance of additional modes at
$\omega>\omega_{DE}$ for superconductor and ferrimagnet in direct contact. The
interplay between a magnetically ordered film in the vicinity of a
superconductor is an interesting and open issue. In the first instance, we can
assume that the Meissner effect (perfect diamagnetism) imposes the expulsion
of the oscillating field at the interface. Intuitively this can be visualized
as a superconducting plane reflecting the image of the magnetic excitations in
YIG [59, 60]. Analysis of the portion of the spectrum can be attempted by
simulations as suggested in [60], yet the dispersion law may depend to a large
extent on the specific materials and geometry of the problem. We cannot
exclude other mechanisms, like those related to vortex configuration and
dynamics [39], anisotropy-induced surface pinning and other interfacial
effects [61], although these were reported for thin ferromagnetic films.
Particular effects can also be induced on the superconductor and these may
depend on the type of superconducting wave function and on its coherence
length. Whilst our experiments represent a case study involving insulating YIG
and the high critical temperature YBCO superconductor that merits further
attention, this issue goes beyond the scope of our research that is focused on
the ultrastrong coupling regime which, instead, is achieved also in different
conditions and geometries.

In the experiments with YIG film and resonators (Fig. 4) the profile of the
resonator field (Fig. 2) may introduce additional quantization of the spin
wave modes [62]. Consistently with Fig. 3(c)), we note that the fitted value
of the frequency shift $\Delta/2\pi=2.05$ GHz is
$\omega_{DE}<\Delta<\omega_{M}+\omega_{H}$, thus within the band of spin wave
resonance modes observed for YIG and YBCO CPW line in direct contact. This
frequency shift can be captured by finite-element electromagnetic simulations
[25], in which YIG film and resonator are respectively modelled as gyrotropic
medium and perfect electric conductor (Appendix A).

Figure 5: Best fit of the transmission spectrum obtained in case $\\#$C with
the YIG film pressed on top of the superconducting YBCO CPW. The obtained
parameters are: $\omega_{c}/2\pi=8.65$ GHz, $\Delta/2\pi=2.05$ GHz,
$\lambda/2\pi=2.002$ GHz, $\alpha/\sqrt{2\pi}=3\times 10^{-3}\
\text{GHz}^{\frac{3}{2}}$.

The collective coupling strength $\lambda/2\pi=2.002$ GHz obtained by fitting
the experimental data with Eq. 7 (Fig. 5) can be spelled out as
$\lambda=g_{s}\sqrt{2s_{\mathrm{Fe}}N}$, where $s_{\mathrm{Fe}}=5/2$ is the
ground state spin of YIG and $N=1.8\times 10^{15}$ is the total number of
spins. The latter resulted much higher than the mean number of photons in the
resonator [63, 42]. The spin-photon coupling is $g_{s}=21$ Hz as derived in
Sect. III. Considering the spin density $\rho=2\times
10^{28}~{}\mathrm{m}^{-3}$ [3] and the effective area of
$4~{}\mathrm{mm}\times 45~{}\mu\mathrm{m}$ in which the microwave field
overlaps the YIG film (Fig. 2(a,b)), we can estimate an upper bound
$d^{\prime}\approx 500~{}\mathrm{nm}\ll d$ for the thickness of the portion of
YIG film coupled with the CPW resonator. This suggests that the magnon modes
that are effectively coupled to the cavity one are located in close vicinity
with the superconducting resonator. This observation is supported by similar
experiments we performed with 20 $\mu$m thick YIG film, which provide results
quite similar to those obtained with the 5 $\mu$m YIG film, thus confirming
that the coupling is confined within few $\mu$m (Appendix B). The ratio
between the collective couplings in configuration $\\#$C and $\\#$B is
$g_{s,C}/g_{s,B}\sqrt{N_{C}/N_{B}}\approx 2$ (Fig. 4(b,c)). From finite-
element simulations, we estimated $g_{s,C}/g_{s,B}\approx 2$ due to the
exponential decay of $b_{vac}$ with $z$ distance (Fig. 2(c)). We therefore
expect $N_{C}\approx N_{B}$. These observations indicate that in our
experiments the achievement of the USC regime takes place mainly as a
consequence of the optimized magnon-photon coupling strength.

The Hopfield model we used to analyze our spectra allows us to overcome the
RWA approximation generally used in previous works with magnetic systems, thus
providing a quantum description of the problem that can be applied to safely
explore the USC regime. We stress that the analysis of the spectrum reported
in Sect. VI is quite robust since it relies on a minimal set of free
parameters. Obviously, by introducing more degrees of freedom in the
Hamiltonian the spectrum can also be fitted well. For instance, one may wonder
whether more magnetic modes couple to cavity photons simultaneously. Results
of simulations with three (and more) modes are reported in Appendix C. It
always results that the best fit is obtained with one dominant mode, ultra-
strongly coupled to the resonator’s one, plus additional modes coupled more
weakly. As a matter of fact, on the basis of the sole fitting of the
polaritonic branches, one cannot exclude that additional magnetic excitations
are involved in our and in similar experiments reported in the literature.
However, the representation with only one effective mode allows putting
stringent bounds to the other terms of the Hamiltonian. Specifically, this is
the case for the diamagnetic term, that in our analysis results vanishingly
small. This result is consistent to what is expected for a pure spin
interaction, although the issue is still debated in the literature. For
instance, we just mention that conclusions reported in Ref. [35] lead to a
finite diamagnetic contribution yet, as those authors concluded, this may
arise from surface plasmonic modes or by the different nature of magnetic
material (permalloy). In our case, the absence of this diamagnetic term may be
relevant for the observation of superradiant phase transition expected for
$\lambda/\omega_{c}>0.5$ (see discussion in Appendix D), making pure spin
systems interesting and unique in this perspective.

As concerns a possible route for applications, the coexistence of
superconductivity and magnetism is not trivial since the presence of a
magnetic field can be detrimental to the superconductor. In our experiments,
the use of high critical temperature YBCO resonators, resilient to high
magnetic field [38], indicates a good option for the realization of this kind
of hybrid device. Our experiments also show that the contact between
superconductor and magnet, and in particular the vanishing gap between the
two, is critical to enhancing the coupling between magnetic and microwave
modes (Fig. 4(b,c)). The direct growth of YIG on top of superconducting oxide
is not straightforward but there can be different options to overcome this
technical issue [64, 65]. The availability of commercial YIG films of
excellent quality allows the easy implementation of our experiment in
different geometries.

##  VIII Conclusions

In summary, our experimental and theoretical results show the achievement of
the ultrastrong coupling with a YIG film positioned in direct contact with a
superconducting CPW resonator. The obtained collective coupling strength
$\lambda/2\pi=2.002$ GHz is improved by at least one order of magnitude with
respect to previous reports involving YIG films in 3D cavities [66] or bulk
YIG crystals in planar resonators [30, 31]. The estimated $\lambda/\omega_{c}$
ratio of 0.23 and the cooperativity $C=5\times 10^{4}$ are among the largest
reported so far for magnetic systems yet we believe that there are still
margins to further increase them for instance by using planar resonators with
more confined mode volumes and higher quality factors, or by means of YIG
films detached from the GGG substrate [64, 65] to reduce magnetic losses [67,
68]. The very small diamagnetic coefficient (compared to standard electric
dipole interactions) observed in these systems makes them suitable for
exploring superradiant phase transitions [17, 18, 19, 20].

###### Acknowledgements.

We thank Prof. Bulat Rameev for useful discussions. This work was partially
supported by European Community through FET Open SUPERGALAX project (grant
agreement No. 863313) and by NATO Science for Peace and Security Programme
(NATO SPS Project No. G5859). MM acknowledges TUBITAK-BIDEB for the support
under the 2219 scholarship program. SS acknowledges the Army Research Office
(ARO) (Grant No. W911NF1910065).

##  Appendix A Additional finite-element electromagnetic simulations

Figure A.1: Transmission spectral maps obtained by finite-element
electromagnetic simulations. (a) The YIG film is lifted of 10 $\mu$m with
respect to the surface of the CPW resonator. (b) The YIG film is in contact
with the conducting surfaces of the resonator. Blue dashed lines indicate
$\omega_{FMR}(H_{0})$ while the dotted lines show $\omega_{c}$ in the two
cases.

We carried out electromagnetic simulations with a commercial software (CST
Microwave Studio) to evaluate the scattering parameters of the coupled system
composed by CPW resonator and YIG film. The superconductor was modelled as a
perfect electric conductor, whilst the magnetic film was included on the lower
face of the GGG substrate [69] with a thickness $d=5~{}\mu$m.

Finite-element simulations were carried out by assuming that the precession of
the magnetization in the ferrite can be described by the gyrotropic model, in
which the permeability is modelled as a nonsymmetric Polder tensor with
characteristic frequency dependence [12]. The magnetic dispersion of YIG was
defined by directly introducing the Larmor frequency,
$\omega_{H}/2\pi=\gamma\mu_{0}H_{0}$, and the gyrotropic frequency,
$\omega_{M}/2\pi=\gamma\mu_{0}M_{s}$, as input parameters of the simulator,
with $\gamma$=28.02 GHz/T, $\mu_{0}=4\pi\times 10^{-7}$ H/m and
$\mu_{0}M_{s}=0.245$ T.

The simulation of frequency spectra was repeated for increasing values of the
external magnetic field ($\mathbf{H_{0}}=H_{0}\hat{x}$) to obtain the the
spectral maps shown in Fig. A.1. The evolution of the coupled YIG-resonator
modes displays the appearance of polaritonic branches. The simulations carried
out with the YIG film lifted of 10 $\mu$m (panel (a)) or in contact (panel
(b)) with the CPW resonator, essentially reproduce the experimental trend
shown in Fig. 4(b,c), showing a good correspondence with the measured
splittings of the polaritonic branches. In agreement with the measured
spectra, the spectral maps in Fig. A.1 show that in both cases the lower
polariton branch converge towards $\omega_{FMR}(H_{0})$ at low frequency. The
increase of the spltting from panel (a) to (b) determines a shift of the
anticrossing towards lower magnetic field, which is compatible with the
frequency shift $\Delta$ introduced in Sect. VI.

For the sake of completeness, we mention that finite element simulations fail
to mimic the polariton lineshapes shown in [42]. A more detailed model is
probably required to reproduce broadening effects, this however goes beyond
the scope of our work.

##  Appendix B Effects of YIG films with different thickness

Figure B.2 shows a direct comparison between transmission data taken with the
same CPW resonator and films of different thickness and size. In these
experiments, the sample was held on the resonator with a copper spring, this
resulted in a lower coupling with respect to data reported in the main body of
the article. Panel (a) shows the spectral map taken with the YIG/GGG film
having a thickness of $5~{}\mu$m and an area of $\approx 4\times
3~{}\mathrm{mm}^{2}$. Panel (b) shows data acquired on a YIG/GGG film having a
thickness of $20~{}\mu$m and area of $\approx 5\times 1~{}\mathrm{mm}^{2}$.
The splitting indicated by the blue arrow, corresponding to
$2\lambda/2\pi\approx 3.6$ GHz, is comparable in the two cases.

Figure B.2: Comparison between data obtained with YIG/GGG films of different
thickness ($T=30$ K): (a) $5~{}\mu$m, (b) $20~{}\mu$m. The red arrows indicate
the splitting $2\lambda/2\pi\approx 3.6$ GHz. The blue dashed lines display
$\omega_{FMR}(H_{0})$. Both measurements were carried out by means of the same
YBCO CPW resonator, horizontal dash-dot lines show the frequency of the
fundamental mode in the two cases.

##  Appendix C Multimode fit

Figure C.3: Best fit of the transmission spectrum obtained in case $\\#$C with
the YIG film pressed on top of the superconducting YBCO CPW by considering
three magnonic modes instead of only one, are $\omega_{c}/2\pi=8.64$ GHz,
$\Delta_{0}/2\pi=0.23$ GHz, $\Delta_{1}/2\pi=1.92$ GHz, $\Delta_{2}/2\pi=3.17$
GHz, $\lambda_{0}/2\pi=0.192$ GHz, $\lambda_{1}/2\pi=1.93$ GHz,
$\lambda_{2}/2\pi=0.2$ GHz, $\alpha/\sqrt{2\pi}=0.035\
\text{GHz}^{\frac{3}{2}}$. Among the three modes, only one is ultrastrongly
coupled with the cavity field. Again, the factor $\alpha$ of the diamagnetic
term is very small.

Despite the fitting analysis presented in the main body of the article
involving one magnonic mode is statistically significant on its own, we extend
the model by introducing more magnonic modes. In addition to the results shown
in the main text, we present here the fit results considering three magnonic
modes (Fig. C.3). The method can be straightforwardly extended from the case
of single mode, and the parameters obtained from the fit are
$\omega_{c}/2\pi=8.64$ GHz, $\Delta_{0}/2\pi=0.23$ GHz, $\Delta_{1}/2\pi=1.92$
GHz, $\Delta_{2}/2\pi=3.17$ GHz, $\lambda_{0}/2\pi=0.192$ GHz,
$\lambda_{1}/2\pi=1.93$ GHz, $\lambda_{2}/2\pi=0.2$ GHz,
$\alpha/\sqrt{2\pi}=0.035\ \text{GHz}^{\frac{3}{2}}$. It is worth noting that,
among the three modes, only one results to be ultrastrongly coupled with the
field of the resonator, corroborating the validity of considering only one
magnonic mode. The obtained diamagnetic factor $\alpha$ is slightly larger
than the single mode case.

##  Appendix D On the diamagnetic term

Figure D.4: Left: polariton eigenfrequencies considering only one magnonic
mode and using the same parameters of the fit in the main text, that are:
$\omega_{c}/2\pi=8.65$ GHz, $\Delta/2\pi=2.05$ GHz, $\lambda/2\pi=2.002$ GHz,
$\alpha/\sqrt{2\pi}=3\times 10^{-3}\ \text{GHz}^{\frac{3}{2}}$. Here the
fitted term of the diamagnetic term is relatively small. However, the coupling
$\lambda$ is still too small to show the superradiant phase transition. Right:
Polariton eigenfrequencies with the same parameters as before, except for
$\lambda/2\pi=3$ GHz, which is large enough to achieve the superradiant phase
transition. The continuous blue lines correspond to the case of $\beta=0$ (no
diamagnetic term), clearly showing a critical point in the region of small
magnetic fields. The red dashed lines correspond to
$\beta=\lambda^{2}/\omega_{b}$ which avoids the superradiant phase transition.

It is useful to compare the obtained diamagnetic terms with the standard one
$\beta_{\rm std}=\lambda^{2}/\omega_{b}$, which comes from the minimal
coupling replacement [17, 70] for electric dipolar interactions. Notice that
the expression for $\beta_{\rm std}$ is fixed by gauge-invariance requirements
[70]. This constraint, preventing superradiance phase transitions, does not
hold in the presence of magnetic interactions [17, 18, 19, 20]. We calculate
the ratio $\mathcal{B}\equiv\beta/\beta_{\rm
std}=\alpha\sqrt{\omega_{b}}/\lambda^{2}$, which is $\mathcal{B}\approx 0.002$
and $\mathcal{B}\approx 0.027$ for the single mode and three modes cases,
respectively. These results pave the way for the transition to a superradiant
phase [17], as the obtained values are extremely low. To clarify the role of
the diamagnetic term, Fig. D.4 shows the polariton eigenfrequencies in three
different conditions. The left panel shows the eigenfrequencies using the same
parameters of the fit obtained in the main text, that are:
$\omega_{c}/2\pi=8.65$ GHz, $\Delta/2\pi=2.05$ GHz, $\lambda/2\pi=2.002$ GHz,
$\alpha/\sqrt{2\pi}=3\times 10^{-3}\ \text{GHz}^{\frac{3}{2}}$. The coupling
strength is still too small to achieve a superradiant phase transition, but a
slightly higher value would be enough to achieve it. Indeed, the continuous
blue lines of the right panel show the polariton eigenfrequencies in the
absence of the diamagnetic term and a larger coupling ($\beta=0$ and
$\lambda/2\pi=3$ GHz). With these parameters, we can see a critical point in
the region of small magnetic fields. On the other hand, as shown by the red
dashed lines, the superradiant phase transition is forbidden by including the
standard diamagnetic term of the Hopfield model
($\beta=\lambda^{2}/\omega_{b}$).

## References

  * Pirro  _et al._ [2021] P. Pirro, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Advances in coherent magnonics, Nature Reviews Materials 6, 1114 (2021).
  * Lachance-Quirion  _et al._ [2019] D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Hybrid quantum systems based on magnonics, Applied Physics Express 12, 070101 (2019).
  * Zare Rameshti  _et al._ [2022] B. Zare Rameshti, S. Viola Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C.-M. Hu, H. X. Tang, G. E. Bauer, and Y. M. Blanter, Cavity magnonics, Physics Reports 979, 1 (2022).
  * Li  _et al._ [2020] Y. Li, W. Zhang, V. Tyberkevych, W.-K. Kwok, A. Hoffmann, and V. Novosad, Hybrid magnonics: Physics, circuits, and applications for coherent information processing, Journal of Applied Physics 128, 130902 (2020).
  * Zhang  _et al._ [2015] X. Zhang, C.-L. Zou, N. Zhu, F. Marquardt, L. Jiang, and H. X. Tang, Magnon dark modes and gradient memory, Nature Communications 6, 8914 (2015).
  * Bai  _et al._ [2017] L. Bai, M. Harder, P. Hyde, Z. Zhang, C.-M. Hu, Y. P. Chen, and J. Q. Xiao, Cavity mediated manipulation of distant spin currents using a cavity-magnon-polariton, Phys. Rev. Lett. 118, 217201 (2017).
  * Crescini  _et al._ [2018] N. Crescini, D. Alesini, C. Braggio, G. Carugno, D. Di Gioacchino, C. S. Gallo, U. Gambardella, C. Gatti, G. Iannone, G. Lamanna, C. Ligi, A. Lombardi, A. Ortolan, S. Pagano, R. Pengo, G. Ruoso, C. C. Speake, and L. Taffarello, Operation of a ferromagnetic axion haloscope at $$m_a=58$\backslash$,$\backslash$upmu $\backslash$mathrm {ev}$$, The European Physical Journal C 78, 703 (2018).
  * Hisatomi  _et al._ [2016] R. Hisatomi, A. Osada, Y. Tabuchi, T. Ishikawa, A. Noguchi, R. Yamazaki, K. Usami, and Y. Nakamura, Bidirectional conversion between microwave and light via ferromagnetic magnons, Phys. Rev. B 93, 174427 (2016).
  * Yao  _et al._ [2023] B. Yao, Y. S. Gui, J. W. Rao, Y. H. Zhang, W. Lu, and C.-M. Hu, Coherent microwave emission of gain-driven polaritons, Phys. Rev. Lett. 130, 146702 (2023).
  * Roberts  _et al._ [2004] R. W. Roberts, B. A. Auld, and R. R. Schell, Magnetodynamic Mode Ferrite Amplifier, Journal of Applied Physics 33, 1267 (2004).
  * Auld [1963] B. A. Auld, Coupling of electromagnetic and magnetostatic modes in ferrite loaded cavity resonators, Journal of Applied Physics 34, 1629 (1963).
  * Gurevich and Melkov [1996] A. Gurevich and G. Melkov,  _Magnetization Oscillations and Waves (1st ed.)_ (CRC Press., 1996).
  * Maier-Flaig  _et al._ [2017] H. Maier-Flaig, S. Klingler, C. Dubs, O. Surzhenko, R. Gross, M. Weiler, H. Huebl, and S. T. B. Goennenwein, Temperature-dependent magnetic damping of yttrium iron garnet spheres, Phys. Rev. B 95, 214423 (2017).
  * Imamoğlu [2009] A. Imamoğlu, Cavity qed based on collective magnetic dipole coupling: Spin ensembles as hybrid two-level systems, Phys. Rev. Lett. 102, 083602 (2009).
  * Soykal and Flatté [2010] O. O. Soykal and M. E. Flatté, Strong field interactions between a nanomagnet and a photonic cavity, Phys. Rev. Lett. 104, 077202 (2010).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nature Reviews Physics 1, 19 (2019).
  * Nataf and Ciuti [2010] P. Nataf and C. Ciuti, No-go theorem for superradiant quantum phase transitions in cavity qed and counter-example in circuit qed, Nature Communications 1, 72 (2010).
  * Mazza and Georges [2019] G. Mazza and A. Georges, Superradiant quantum materials, Phys. Rev. Lett. 122, 017401 (2019).
  * Andolina  _et al._ [2020] G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, Theory of photon condensation in a spatially varying electromagnetic field, Phys. Rev. B 102, 125137 (2020).
  * Román-Roche  _et al._ [2021] J. Román-Roche, F. Luis, and D. Zueco, Photon condensation and enhanced magnetism in cavity qed, Phys. Rev. Lett. 127, 167201 (2021).
  * Goryachev  _et al._ [2014] M. Goryachev, W. G. Farr, D. L. Creedon, Y. Fan, M. Kostylev, and M. E. Tobar, High-cooperativity cavity qed with magnons at microwave frequencies, Phys. Rev. Applied 2, 054002 (2014).
  * Bourhill  _et al._ [2016] J. Bourhill, N. Kostylev, M. Goryachev, D. L. Creedon, and M. E. Tobar, Ultrahigh cooperativity interactions between magnons and resonant photons in a yig sphere, Phys. Rev. B 93, 144420 (2016).
  * Kostylev  _et al._ [2016] N. Kostylev, M. Goryachev, and M. E. Tobar, Superstrong coupling of a microwave cavity to yttrium iron garnet magnons, Applied Physics Letters 108, 062402 (2016).
  * Zare Rameshti  _et al._ [2015] B. Zare Rameshti, Y. Cao, and G. E. W. Bauer, Magnetic spheres in microwave cavities, Phys. Rev. B 91, 214430 (2015).
  * Bourcin  _et al._ [2023] G. Bourcin, J. Bourhill, V. Vlaminck, and V. Castel, Strong to ultrastrong coherent coupling measurements in a yig/cavity system at room temperature, Phys. Rev. B 107, 214423 (2023).
  * Flower  _et al._ [2019] G. Flower, M. Goryachev, J. Bourhill, and M. E. Tobar, Experimental implementations of cavity-magnon systems: from ultra strong coupling to applications in precision measurement, New Journal of Physics 21, 095004 (2019).
  * Liensberger  _et al._ [2019] L. Liensberger, A. Kamra, H. Maier-Flaig, S. Geprägs, A. Erb, S. T. B. Goennenwein, R. Gross, W. Belzig, H. Huebl, and M. Weiler, Exchange-enhanced ultrastrong magnon-magnon coupling in a compensated ferrimagnet, Phys. Rev. Lett. 123, 117204 (2019).
  * Białek  _et al._ [2021] M. Białek, J. Zhang, H. Yu, and J.-P. Ansermet, Strong coupling of antiferromagnetic resonance with subterahertz cavity fields, Phys. Rev. Applied 15, 044018 (2021).
  * Everts  _et al._ [2020] J. R. Everts, G. G. G. King, N. J. Lambert, S. Kocsis, S. Rogge, and J. J. Longdell, Ultrastrong coupling between a microwave resonator and antiferromagnetic resonances of rare-earth ion spins, Phys. Rev. B 101, 214414 (2020).
  * Huebl  _et al._ [2013] H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, High cooperativity in coupled microwave resonator ferrimagnetic insulator hybrids, Phys. Rev. Lett. 111, 127003 (2013).
  * Morris  _et al._ [2017] R. G. E. Morris, A. F. van Loo, S. Kosen, and A. D. Karenowska, Strong coupling of magnons in a YIG sphere to photons in a planar superconducting resonator in the quantum limit, Scientific Reports 7, 11511 (2017).
  * Hou and Liu [2019] J. T. Hou and L. Liu, Strong coupling between microwave photons and nanomagnet magnons, Phys. Rev. Lett. 123, 107702 (2019).
  * Golovchanskiy  _et al._ [2021a] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, M. Weides, V. V. Ryazanov, A. A. Golubov, A. V. Ustinov, and M. Y. Kupriyanov, Ultrastrong photon-to-magnon coupling in multilayered heterostructures involving superconducting coherence via ferromagnetic layers, Science Advances 7, eabe8638 (2021a).
  * Li  _et al._ [2022] Y. Li, V. G. Yefremenko, M. Lisovenko, C. Trevillian, T. Polakovic, T. W. Cecil, P. S. Barry, J. Pearson, R. Divan, V. Tyberkevych, C. L. Chang, U. Welp, W.-K. Kwok, and V. Novosad, Coherent coupling of two remote magnonic resonators mediated by superconducting circuits, Phys. Rev. Lett. 128, 047701 (2022).
  * Golovchanskiy  _et al._ [2021b] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, A. A. Golubov, M. Y. Kupriyanov, V. V. Ryazanov, and A. V. Ustinov, Approaching deep-strong on-chip photon-to-magnon coupling, Phys. Rev. Applied 16, 034029 (2021b).
  * Macêdo  _et al._ [2021] R. Macêdo, R. C. Holland, P. G. Baity, L. J. McLellan, K. L. Livesey, R. L. Stamps, M. P. Weides, and D. A. Bozhko, Electromagnetic approach to cavity spintronics, Phys. Rev. Appl. 15, 024065 (2021).
  * Martínez-Losa del Rincón  _et al._ [2023] S. Martínez-Losa del Rincón, I. Gimeno, J. Pérez-Bailón, V. Rollano, F. Luis, D. Zueco, and M. J. Martínez-Pérez, Measuring the magnon-photon coupling in shaped ferromagnets: Tuning of the resonance frequency, Phys. Rev. Appl. 19, 014002 (2023).
  * Ghirri  _et al._ [2015] A. Ghirri, C. Bonizzoni, D. Gerace, S. Sanna, A. Cassinese, and M. Affronte, Yba2cu3o7 microwave resonators for strong collective coupling with spin ensembles, Applied Physics Letters 106, 184101 (2015).
  * Niedzielski  _et al._ [2023] B. Niedzielski, C. L. Jia, and J. Berakdar, Magnon-fluxon interaction in coupled superconductor/ferromagnet hybrid periodic structures, Phys. Rev. Appl. 19, 024073 (2023).
  * Bonizzoni  _et al._ [2022] C. Bonizzoni, M. Maksutoglu, A. Ghirri, J. van Tol, B. Rameev, and M. Affronte, Coupling Sub-nanoliter BDPA Organic Radical Spin Ensembles with YBCO Inverse Anapole Resonators, Applied Magnetic Resonance (2022).
  * Ghirri  _et al._ [2020] A. Ghirri, C. Herrero, S. Mazerat, T. Mallah, O. Moze, and M. Affronte, Coupling nanostructured csnicr prussian blue analogue to resonant microwave fields, Advanced Quantum Technologies 3, 1900101 (2020).
  * [42] See Supplemental Material at [URL will be inserted by publisher] for experimental details, additional transmission spectroscopy data and simulations .
  * Tosi  _et al._ [2014] G. Tosi, F. A. Mohiyaddin, H. Huebl, and A. Morello, Circuit-quantum electrodynamics with direct magnetic coupling to single-atom spin qubits in isotopically enriched 28si, AIP Advances 4, 087122 (2014).
  * Kittel [1948] C. Kittel, On the theory of ferromagnetic resonance absorption, Phys. Rev. 73, 155 (1948).
  * Kittel [2004] C. Kittel,  _Introduction to Solid State Physics_, 8th ed. (Wiley, 2004).
  * Kajiwara  _et al._ [2010] Y. Kajiwara, K. Harii, S. Takahashi, J. Ohe, K. Uchida, M. Mizuguchi, H. Umezawa, H. Kawai, K. Ando, K. Takanashi, S. Maekawa, and E. Saitoh, Transmission of electrical signals by spin-wave interconversion in a magnetic insulator, Nature 464, 262 (2010).
  * Maksymov and Kostylev [2015] I. S. Maksymov and M. Kostylev, Broadband stripline ferromagnetic resonance spectroscopy of ferromagnetic films, multilayers and nanostructures, Physica E: Low-dimensional Systems and Nanostructures 69, 253 (2015).
  * Kennewell  _et al._ [2007] K. J. Kennewell, M. Kostylev, and R. L. Stamps, Calculation of spin wave mode response induced by a coplanar microwave line, Journal of Applied Physics 101, 09D107 (2007).
  * Wang  _et al._ [2020] L. Wang, Z. Lu, X. Zhao, W. Zhang, Y. Chen, Y. Tian, S. Yan, L. Bai, and M. Harder, Magnetization coupling in a yig/ggg structure, Phys. Rev. B 102, 144428 (2020).
  * Tsutsumi  _et al._ [1996] M. Tsutsumi, T. Fukusako, and S. Yoshida, Propagation characteristics of the magnetostatic surface wave in the ybco-yig film-layered structure, IEEE Transactions on Microwave Theory and Techniques 44, 1410 (1996).
  * Oates  _et al._ [1997] D. Oates, A. Pique, K. Harshavardhan, J. Moses, F. Yang, and G. Dionne, Tunable ybco resonators on yig substrates, IEEE Transactions on Applied Superconductivity 7, 2338 (1997).
  * Diniz  _et al._ [2011] I. Diniz, S. Portolan, R. Ferreira, J. M. Gérard, P. Bertet, and A. Auffèves, Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories, Phys. Rev. A 84, 063810 (2011).
  * Ghirri  _et al._ [2016] A. Ghirri, C. Bonizzoni, F. Troiani, N. Buccheri, L. Beverina, A. Cassinese, and M. Affronte, Coherently coupling distinct spin ensembles through a high-${T}_{c}$ superconducting resonator, Phys. Rev. A 93, 063855 (2016).
  * Forn-Díaz  _et al._ [2019b] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong coupling regimes of light-matter interaction, Rev. Mod. Phys. 91, 025005 (2019b).
  * Holstein and Primakoff [1940] T. Holstein and H. Primakoff, Field dependence of the intrinsic domain magnetization of a ferromagnet, Phys. Rev. 58, 1098 (1940).
  * Hopfield [1958] J. J. Hopfield, Theory of the contribution of excitons to the complex dielectric constant of crystals, Phys. Rev. 112, 1555 (1958).
  * Kalinikos and Slavin [1986] B. A. Kalinikos and A. N. Slavin, Theory of dipole-exchange spin wave spectrum for ferromagnetic films with mixed exchange boundary conditions, Journal of Physics C: Solid State Physics 19, 7013 (1986).
  * Demokritov and Slavin [2021] S. O. Demokritov and A. N. Slavin, Spin waves, in  _Handbook of Magnetism and Magnetic Materials_, edited by J. M. D. Coey and S. S. Parkin (Springer International Publishing, Cham, 2021) pp. 281–346.
  * Golovchanskiy  _et al._ [2018a] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, V. V. Bolginov, V. V. Ryazanov, A. A. Golubov, and A. V. Ustinov, Ferromagnet/superconductor hybridization for magnonic applications, Advanced Functional Materials 28, 1802375 (2018a).
  * Golovchanskiy  _et al._ [2018b] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, V. V. Ryazanov, A. A. Golubov, and A. V. Ustinov, Modified dispersion law for spin waves coupled to a superconductor, Journal of Applied Physics 124, 233903 (2018b).
  * Lee  _et al._ [2023] B. H. Lee, T. Fakhrul, C. A. Ross, and G. S. D. Beach, Large anomalous frequency shift in perpendicular standing spin wave modes in biyig films induced by thin metallic overlayers, Phys. Rev. Lett. 130, 126703 (2023).
  * Bayer  _et al._ [2006] C. Bayer, J. Jorzick, S. O. Demokritov, A. N. Slavin, K. Y. Guslienko, D. V. Berkov, N. L. Gorn, M. P. Kostylev, and B. Hillebrands, Spin-wave excitations in finite rectangular elements, in  _Spin Dynamics in Confined Magnetic Structures III_, edited by B. Hillebrands and A. Thiaville (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 57–103.
  * Sage  _et al._ [2011] J. M. Sage, V. Bolkhovsky, W. D. Oliver, B. Turek, and P. B. Welander, Study of loss in superconducting coplanar waveguide resonators, Journal of Applied Physics 109, 063915 (2011).
  * Trempler  _et al._ [2020] P. Trempler, R. Dreyer, P. Geyer, C. Hauser, G. Woltersdorf, and G. Schmidt, Integration and characterization of micron-sized yig structures with very low gilbert damping on arbitrary substrates, Applied Physics Letters 117, 232401 (2020).
  * Kosen  _et al._ [2019] S. Kosen, A. F. van Loo, D. A. Bozhko, L. Mihalceanu, and A. D. Karenowska, Microwave magnon damping in yig films at millikelvin temperatures, APL Materials 7, 101120 (2019).
  * Zhang  _et al._ [2016] X. Zhang, C. Zou, L. Jiang, and H. X. Tang, Superstrong coupling of thin film magnetostatic waves with microwave cavity, Journal of Applied Physics 119, 023905 (2016).
  * Danilov  _et al._ [1989] V. V. Danilov, D. L. Lyfar’, Y. V. Lyubon’ko, A. Y. Nechiporuk, and S. M. Ryabchenko, Low-temperature ferromagnetic resonance in epitaxial garnet films on paramagnetic substrates, Soviet Physics Journal 32, 276 (1989).
  * Mihalceanu  _et al._ [2018] L. Mihalceanu, V. I. Vasyuchka, D. A. Bozhko, T. Langner, A. Y. Nechiporuk, V. F. Romanyuk, B. Hillebrands, and A. A. Serga, Temperature-dependent relaxation of dipole-exchange magnons in yttrium iron garnet films, Phys. Rev. B 97, 214405 (2018).
  * Connelly  _et al._ [2021] D. A. Connelly, H. R. O. Aquino, M. Robbins, G. H. Bernstein, A. Orlov, W. Porod, and J. Chisum, Complex permittivity of gadolinium gallium garnet from 8.2 to 12.4 ghz, IEEE Magnetics Letters 12, 1 (2021).
  * Garziano  _et al._ [2020] L. Garziano, A. Settineri, O. Di Stefano, S. Savasta, and F. Nori, Gauge invariance of the dicke and hopfield models, Phys. Rev. A 102, 023718 (2020).

##  Appendix E Supplemental material

###  E.1 Experimental details

Figure S1: Schematic representation of (a) microstrip resonator and (b) CPW
resonator. In (a) the regions covered with Ag are black while the uncovered
alumina substrate is white. In (b) YBCO regions are black while sapphire is
white. The images on the right show blow-ups of launcher and central line of
the CPW resonator as indicated by the different colors.

Microstrip lines were fabricated by optical lithography and wet etching of Ag
films (thickness $\approx 3~{}\mu\mathrm{m}$), which were thermally evaporated
onto alumina $14\times 10\times 0.63~{}\mathrm{mm}^{3}$ substrates [41]. The
central strip is $w^{\prime}=500~{}\mu$m wide, while two
$\delta^{\prime}=300~{}\mu$m coupling gaps were made to define the
corresponding resonator (Fig. S1(a)). Superconducting coplanar waveguide (CPW)
broadband lines and resonators (Fig. S1(b)) were fabricated from commercial
YBa2Cu3O7 (YBCO) films (thickness 330 nm) deposited on a sapphire substrate
and diced into $8\times 5\times 0.43~{}\mathrm{mm^{3}}$ blocks. Etching was
carried out by Ar plasma in a reactive ion etching (RIE) chamber. The central
conductor has typical width $w=(17\pm 1)~{}\mu$m and is separated by $s=(14\pm
1)~{}\mu$m from the lateral ground planes (Fig. S2(c)). In the CPW resonator,
the coupling gaps have width $\delta=140~{}\mu$m.

Figure S2: Photographs of (a) Ag/alumina microstrip resonator and (b)
YBCO/sapphire CPW transmission line, both with the YIG/GGG film positioned on
top. (c) Blow-up of the central region of the CPW line. The lateral widths are
$w=(17\pm 1)~{}\mu$m and $s=(14\pm 1)~{}\mu$m.

We studied $5~{}\mu$m thick yttrium iron garnet (YIG) films which were grown
by liquid-phase epitaxy on a gadolinium gallium garnet (GGG) substrate with
(111) crystallographic orientation (Matesy GmbH). Reflection ($S_{11}$) and
transmission ($S_{21}$) spectra were acquired in the 0.1-18 GHz range by means
of a Vector Network Analyzer (VNA). The experiments were carried out in the
10-50 K temperature ($T$) range. Planar transmission lines and resonators were
installed in a cryomagnetic set-up (Quantum Design PPMS), having variable
temperature control and external magnetic field up to 7 T, by means of a
cryogenic insert wired with a pair of silver-plated stainless steel coaxial
cables connecting the low temperature stage to the external circuitry. The
microstrip transmission line and resonator (Fig. S2(a)) were mounted in a
brass box and connected to the coaxial line by means of two metallic pins
positioned on the microstrip. Coplanar waveguide (CPW) transmission line and
resonator (Fig. S2(b)) were installed in a copper box and wire bonded (Al
wire) to a printed-circuit board. Gold (thickness 200 nm) pads were patterned
on the edges of the YBCO resonator to facilitate the bonding.

Figure S3: Transmission spectra of the YIG film acquired by the microstrip
broadband line at different temperatures. (a) Comparison between fixed-
frequency field-swept spectra taken at 10 K and 50 K. (b,c) Spectral maps
measured by sweeping the frequency at progressively increasing magnetic field.
The temperature is (b) 10 K and (c) 50 K.

In all the measurements, the incident power supplied to the planar
transmission line is $P_{inc}\approx-8$ dBm. Spectra taken with $P_{inc}$ in
the range between -28 dBm and 2 dBm showed no variations with the microwave
power. The number of photons can be calculated as [63]

$\displaystyle n=\frac{P_{inc}Q_{L}10^{-IL/20}}{\pi hf_{0}^{2}},$ (8)

where $f_{0}=\omega_{0}/2\pi$ is the fundamental mode frequency, $Q_{L}$ is
the loaded quality factor, $IL$ is the insertion loss and $h$ is the Planck
constant. Depending on the resonator’s parameters (see below), the estimated
number of photons results $n_{mic}\approx 7\times 10^{8}$ for the microstrip
resonator and $n_{CPW}\approx 3\times 10^{10}$ for the CPW resonators. In all
the experiments the number of photons resulted much smaller than the estimated
number of spins.

###  E.2 Additional transmission spectroscopy data

The YIG film has been initially characterized by means of the broadband
microstrip line, which allowed the acquisition of ferromagnetic resonance
spectra at different excitation frequencies in the 4.5-16.5 GHz range (Fig.
S3). Fixed-frequency field-swept transmission spectra taken at 10 and 50 K
show similar features below $\approx 13.5$ GHz while little differences in
absorption intensity and line position can be observed for higher frequencies
(Fig. S3(a)). Spectra taken with the microstrip line are characterized by
maximum absorption amplitude in correspondence to the Kittel mode, while
additional higher modes have lower amplitude. The measured spectral maps (Fig.
S3(b,c)) show that the dependence on the lowest mode can be fitted with the
Kittel formula (Eq. 1 in the main text) by considering the saturation
magnetization $\mu_{0}M_{s}=0.245$ T both for 10 and 50 K.

Figure S4: Transmission spectra of the YIG film acquired by the CPW broadband
line. (a) Comparison between frequency spectra taken for cases $\\#B$ and
$\\#C$. (b) Direct comparison between the spectra at 0.21 T after background
subtraction. Dashed lines show the integral calculated between 8 and 12 GHz.

Fig. S4(a) shows broadband transmission-vs-frequency spectra acquired by YBCO
CPW lines (case $\\#B$ and $\\#C$, as defined in the main text) for zero and
0.21 T applied magnetic field. In panel (b) the spectra taken at 0.21 T are
directly compared to evidence the increased depth and width of the $S_{21}$
spectrum when YIG and YBCO CPW line are put in direct contact. To quantify
this point, we calculate the total absorption related to the transitions in
Fig. 3(b) and (c) in the main text, which can be obtained by integrating
$S_{21}(\omega)$ with respect to the frequency. For $H_{0}=0.21$ T, the ratio
of the integrals calculated between 8 and 12 GHz (Fig. S4) gives
$I_{C}/I_{B}=3.49$. The total number of magnon excitations in this frequency
range thus results much larger when the YIG film is positioned in direct
contact with the YBCO resonator.

Figure S5: Transmission spectral maps obtained with the same YIG film and
different resonator. The background transmission curve has been subtracted
from the data. (a) Metallic microstrip resonator ($T=50$ K, case $\\#$A).
Inset: transmission-vs-frequency spectrum taken at $H_{0}=0.45$ T. (b)
Superconducting CPW resonator ($T=30$ K, case $\\#$B). (c) CPW resonator with
the YIG film pressed against the superconductor ($T=30$ K, case $\\#$C). In
(b) and (c) the insets show the transmission-vs-frequency spectrum taken at
$H_{0}=2$ T.  Figure S6: Evolution of the derivative spectra (case $\\#$C)
plotted as a function of frequency for fixed applied magnetic fields. Green
squares and blue circles indicate respectively the position of the lower and
upper polaritonic branches, yellow triangles shows the evolution of the FMR
mode.

Figure S5 shows $S_{21}(\omega,H_{0})$ transmission spectral maps acquired
with microstrip and CPW resonators corresponding to the same dataset reported
in Fig. 4 in the main text. Since resonator and YIG are already coupled in
zero magnetic field, we exploited the resilience of YBCO in high magnetic
field [38] to characterize the uncoupled CPW resonators. From transmission-vs-
frequency spectra acquired at $H_{0}=2$ T, we obtain the resonator frequencies
$\omega_{0,B}/2\pi=9.1$ GHz and $\omega_{0,C}/2\pi=8.5$ GHz, respectively for
case $\\#$B and $\\#$C (Fig. S5(b) and (c)). The discrepancy between the
latter and $\omega_{c}/2\pi=8.65$ GHz obtained from the fit (Fig. 5 in the
main text), indicates that at 2 T resonator and magnon modes are still weakly
coupled, due to their large collective coupling strength. Since CPW resonators
are fabricated from the same lithographic mask, the differences between
$\omega_{0,B}$ and $\omega_{0,C}$ originate from specific experimental
details, including sample installation in diverse ways. The insertion loss is
$IL_{B}=32$ dB and $IL_{C}=30$ dB. In both cases, the loaded quality factor is
$Q_{L}\approx 1000$, which is dominated by losses introduced by the GGG
substrate (permittivity $\epsilon_{r}=11.99$ and $\tan{\delta}=5.2\times
10^{-3}$ at room temperature [69]). The corresponding decay rate of the
resonator ($\kappa_{0}=\omega_{0}/Q_{L}$) amounts to $\kappa_{0,B}/(2\pi)=9$
MHz and $\kappa_{0,C}/(2\pi)=8$ MHz. For comparison, the microstrip resonator
(case $\\#$A) displays $\omega_{0,A}/2\pi=9.5$ GHz, $IL_{A}=42$ dB and
$Q_{L}\approx 120$ (Inset in Fig. S5(a)), resulting in the decay rate
$\kappa_{0,A}/(2\pi)=79$ MHz.

Figure S7: Transmission spectral data taken with the YIG film in contact with
the CPW resonator. $S_{21}$ (a) and $\partial S_{21}/\partial H_{0}$ (b)
spectra plotted as a function of the magnetic field for different frequencies,
as indicated. (c) Spectral map acquired for $H_{0}$ spanning between 0.38 to
-0.38 T.

Figure S6 shows the evolution of derivative spectra corresponding to the
dataset in Fig. 4(c) of the main text and plotted as a function of the
frequency for fixed fields. The lower polaritonic branch typically shows
narrower and more marked derivative peaks with respect to the upper one. This
behavior is evident also in the spectra plotted in Fig. S7(a,b) as a function
of the magnetic field. Transmission $S_{21}$ spectra taken at positive and
negative field shows a quite similar behavior (panel (c)).

The coupling between YIG film and superconducting CPW resonator was studied by
means of spectral maps acquired at $T=10,~{}30$ and 50 K (Fig. S8). In this
range of temperatures, the coupling strength remains essentially constant,
while it is possible to observe a broadening of the polaritonic modes for
decreasing temperatures. Since the quality factor of the YBCO resonator
typically increases from 50 to 10 K [38], this behavior can be ascribed to the
effect of the GGG substrate, which is known to exhibit a paramagnetic behavior
below 70 K that is reported to increase the damping of magnetization
precession in the YIG film grown on its surface [67, 68, 65]. We just mention
here that narrow linewidths, lower than those measured at 300 K, were reported
for temperatures as low as mK with substrate-free YIG films obtained by
detaching the GGG in different manners [65, 64].

Figure S8: Comparison between derivative maps taken at different temperatures,
as indicated. The YIG/GGG film had an area of $4\times 2~{}\mathrm{mm}^{2}$
and a thickness of $5~{}\mu$m. The red dashed lines show the Kittel equation
($\mu_{0}M_{s}=0.245$ T) while the white dash-dot lines indicate the frequency
of the resonator determined at $H_{0}=2$ T. Cyan arrows have the same length
and indicate the splitting of the polaritonic branches. Figure S9:
Transmission spectral map measured with the GGG substrate (volume $4\times
2\times 0.5~{}\mathrm{mm}^{3}$) positioned in contact with the CPW resonator.
The red dashed line shows the Kittel equation calculated with
$\mu_{0}M_{s}=0.245$ T. The upper panel shows the $\partial S_{21}/\partial
H_{0}$-vs-$H_{0}$ spectrum taken at 9 GHz.

Additional effects reported for YIG films grown on paramagnetic GGG are the
deviation of the resonance field from that of a free thin YIG plate [67] and
the coupling between the magnetizations of ferrimagnetic YIG and paramagnetic
GGG [49]. Fig. S9 shows the transmission map of a YIG/GGG film measured with
the GGG side in contact with the CPW resonator. In contrast with Fig. S8(b),
only a weak absorption line is observed in this case. From such spectra, as
well as from data in Figs. S3 and S9 or in Fig. 3 of the main text, the
presence of additional coupled modes [49] was not observed. This likely
follows from the different conditions in our experiment, particularly lower
microwave power and narrower (and superconducting) CPW line with respect to
Ref. [49]. Thus, the magnetization coupling between YIG and GGG cannot be
inferred from our data. However, the presence of the paramagnetic GGG
substrate probably influences the value of the saturation magnetization
$\mu_{0}M_{s}=0.245$ T entering in the Kittel equation (Eq. 1 in the main
text), that well reproduces the field dependence of the FMR line for
temperatures between 50 and 10 K (Fig. S3). We stress that little changes in
the choice of $\mu_{0}M_{s}$ by no means affect the interpretation of the
results reported in our work.

###  E.3 Additional finite-element electromagnetic simulations

Figure S10: Electromagnetic simulation of the bare CPW resonator. (a)
Comparison between experimental ($T=30$ K) and simulated spectra. (b) Spectral
map showing the normalized value of $h_{ac}$ around the central conductor. (c)
$x$-dependence of the $y$-averaged value of $h_{ac}$ calculated for $z$
heights of 0 and $10~{}\mu$m. (d) Calculated dependence of the microwave
magnetic field along the $z$-axis. The solid line show the fit with Eq. 9.

Preliminary simulations of the electromagnetic field profile with commercial
codes (CST Microwave Studio) have been used as a guide to evaluate and the
distribution of the microwave field and to optimize the effective volume of
interaction. The CPW resonator was modelled with realistic geometry,
dimensions and materials parameters while the superconducting film was
approximated as a perfect electric conductor (PEC).

We initially calculated the transmission spectrum of the bare resonator . To
obtain a substantial correspondence with the measured spectrum (Fig. S10(a)),
we modelled the sapphire substrate (volume $8\times 5\times
0.43~{}\mathrm{mm}^{3}$) with the permittivity $\epsilon_{r}=10.74$, while
losses where included by forcing $\tan\delta$ to values higher than those
expected for sapphire. The frequency of the fundamental mode results in good
agreement with the value obtained from
$\omega_{c,bare}/2\pi=c/2\sqrt{\epsilon_{eff}}l$ being $l=6.02$ mm the length
of the central strip (Fig. S1(b)) and $\epsilon_{eff}=5.77$ the effective
permittivity calculated for the CPW line.

The simulated distribution of the magnetic component of the fundamental mode
of the resonator ($h_{ac}$) shows that the field is localized along $y$ in a
region of approximate width $w+2s=45~{}\mu\mathrm{m}$ around the central
conductor (panel (b)). The value of $h_{ac}$ shows a pronounced decay with
increasing $z$ distance from the upper surface of the resonator ($z=0$), as it
can be seen in the profiles calculated for $z=0$ and $10~{}\mu$m and plotted
as a function of $x$ (Fig. S10(c)). To quantify this trend, we calculated the
average of $\mathbf{h_{ac}}$ in regions centered in the middle of the CPW
resonator (area $4~{}\mathrm{mm}\times 45~{}\mu\mathrm{m}$) with progressively
increasing $z$ height. The obtained points ($\bar{h}_{ac}(z)$) show a quasi
exponential dependence from $z$ (Fig. S10(d)) that roughly follows

$\displaystyle\bar{h}_{ac}(z)=h_{max}\exp(-z/\eta)$ (9)

being $h_{max}=\bar{h}_{ac}(z=0)$ and $\eta=13.5~{}\mu$m. Since
$\bar{h}_{ac}(10~{}\mathrm{\mu m})/h_{max}=0.45$ and by assuming that the
normalized spin-photon coupling $g_{s}(z)/g_{s}(z=0)$ scales as
$\bar{h}_{ac}(z)/h_{max}$, we obtain the ratio
$g_{s}(z=0)/g_{s}(\mathrm{z=10~{}\mu m})=2.2$. This value is in good agreement
with the ratio of the anticrossing splittings in configurations $\\#$C and
$\\#$B derived from Fig. 4 in the main text.

The absolute value of $h_{max}$ obtained from the simulation can be used to
estimate the amplitude of the zero-point vacuum fluctuation, $b_{vac}$. Being
$b_{max}=\mu_{0}h_{max}=26$ mT the value obtained for $P_{inc}=0.5$ W, we
obtain $b_{vac}=b_{max}\sqrt{P_{vac}/P_{inc}}\approx 3$ nT by considering the
threshold power $P_{vac}\approx 5$ fW for single-photon operation (Eq. 8). On
the other hand, the vacuum magnetic fluctuation can be calculated from [43]

$\displaystyle
b_{vac}\approx\frac{\mu_{0}\omega_{c}}{4w}\sqrt{\frac{h}{Z_{0}}}$ (10)

with $h=6.626\times 10^{-34}$ J s and $Z_{0}=58~{}\mathrm{\Omega}$. From Eq.
10 we get $b_{vac}\approx 3$ nT in excellent agreement with the value
resulting from finite-element simulations.

As a next step we modeled the effect of a GGG substrate having volume of
$4\times 3\times 0.5~{}\mathrm{mm}^{3}$ (Fig. S10). By using as permittivity
$\epsilon_{r}=11.99$ and dielectric loss tangent $\tan{\delta}=5.2\times
10^{-3}$ for GGG [69], the simulated fundamental mode frequency
$\omega_{c,sim}/2\pi=8.62$ GHz results in good agreement with
$\omega_{c}/2\pi=8.65$ GHz obtained from the fit of the experimental data
(Fig. 5 in the main text). With respect to the bare resonator, the presence of
the GGG substrate determines the shift of the fundamental mode towards lower
frequency and the broadening of resonant peaks. The latter gives rise to a
higher decay rate $\kappa_{c}$. The distribution of
$\mathbf{h_{ac}}=h_{x}\hat{x}+h_{y}\hat{y}+h_{z}\hat{z}$ simulated with the
GGG substrate above the resonator is in line with that of the bare resonator
(Fig. 2 in the main text): the magnetic antinode is located in the middle of
the resonator, the microwave field is localized around the central conductor
(panel (c)) and $h_{x}$ is negligible as expected for a quasi-transverse
electromagnetic (TEM) mode (panel (d)).

Figure S11: Electromagnetic simulations of the CPW resonator with the GGG
substrate positioned on top. (a) Representation of the model used for the
simulations. Colors indicate different materials: PEC, grey; Au, yellow;
sapphire and GGG, cyan. (b) Simulated reflection $S_{11}$ and transmission
$S_{21}$ spectra showing the comparison between bare resonator and GGG cases.
(c) Vertical section of the simulated $\mathbf{h_{ac}}$ field in the middle of
the resonator (phase $\phi=0$). (d) Blow-up of the CPW line showing the
simulated root-mean square cartesian components of $\mathbf{h_{ac}}$. Letters
in (a,c,d) correspond to $\delta=140~{}\mu$m, $w=17~{}\mu$m and $s=14~{}\mu$m.

###  E.4 Modelling of spin systems in a cavity

Magnetic excitations, i.e. spin waves in magnetically ordered materials, can
be described by considering the precession of the macroscopic magnetic vector
M around the effective local magnetic field $\mathbf{H_{eff}}$ through the
Landau-Lifshizt-Gilbert (LLG) equation, that also accounts for damping by a
phenomenological term characterised by the Gilbert parameter. Exchange and
dipolar interactions between different local values of M can be included in
the LLG equation as effective local fields, $\mathbf{H_{exc}}$ and
$\mathbf{H_{dip}}$ respectively, that contribute to $\mathbf{H_{eff}}$ along
with the externally applied static field $\mathbf{H_{0}}$ [12]. Both exchange
and dipolar interactions are naturally included in a microscopic description
of a regular lattice of interacting spins, the relative spectrum of excitation
shows spin wave modes having characteristic dispersion between energy
($\hbar\omega$) and propagation wavevector ($k$) [3]. For an isotropic
infinite lattice, the low lying states assume parabolic dispersion law
$\omega\approx k^{2}$ with a characteristic stationary ($k$=0) ground state
(Kittel mode). Yet, for a real specimen with a specific shape, the finite size
and the boundary conditions, as well as the way the antenna excites them, can
make the spectrum of magnetic excitations much richer with the appearance of
characteristic bands.

The collective coupling $\lambda$ between magnetic and microwave modes is the
key parameter that needs to be compared, in the first instance, with the
damping of both the magnetic system (i.e. the dissipation rate $\kappa_{m}$)
and the photon losses ($\kappa_{c}$) of the cavity. Since $\lambda$ determines
the rate at which the two systems exchange energy and information at
resonance, when $\lambda\geq\kappa_{m},\kappa_{c}$ we have coherent dynamics
of the two systems, and the so-called strong coupling regime. The dipolar
coupling of single spins to the oscillating magnetic component of radiation
$\mathbf{b_{vac}}$ has the relative coupling strength $g_{s}=\frac{1}{4}\gamma
b_{vac}$ which turns out to be very weak, typically on the order of 0.1 to 50
Hz for $b_{vac}\sim$ nT. For spin ensembles, either paramagnets or
magnetically ordered systems, the collective coupling strength
$\lambda=g_{s}\sqrt{N_{s}}$ is enhanced by a factor scaling as $\sqrt{N_{s}}$
[14], being $N_{s}$ the number of spin transitions in the ensemble interacting
with the electromagnetic radiation. For instance, with $N_{s}\sim
10^{12}-10^{14}$ we can get $\lambda=g_{s}\sqrt{N_{s}}=1$ to 10 MHz. For YIG,
$N_{s}=2s_{\mathrm{Fe}}N$, where $N$ is the number of spins and
$s_{\mathrm{Fe}}=5/2$ is the ground state spin of Fe3+. In a YIG sphere,
$\kappa_{m}$ is of the order of MHz, thus for experiments with cavities having
high quality factor, that is small $\kappa_{c}$, the strong coupling regime is
achieved (see references in the main text).