Quantum electrodynamics with magnetic textures

The spin Hamiltonian of a ferromagnet can be written as:

$$\begin{eqnarray}&&{H}_{S}={H}_{\mathrm{ex}}+{H}_{\mathrm{dip}}+{H}_{{\rm{K}}}+{H}_{\mathrm{DMI}}+{H}_{{\rm{Z}}}.\end{eqnarray} \tag{ 1 }$$

Here, ${H}_{\mathrm{ex}}=-{J}_{{ij}}{\vec{S}}_{i}\,{\vec{S}}_{j}$ is the Heisenberg exchange energy, H_dip is the dipolar term that accounts for the total magnetostatic energy of the system (shape anisotropy), H_K is the magnetocrystalline anisotropy that sets some preferred axis for the magnetization and ${H}_{\mathrm{DMI}}={D}_{{ij}}{\vec{S}}_{i}\times {\vec{S}}_{j}$ is the Dzyaloshinskii–Moriya interaction (DMI) term, with D_ij the strength of the asymmetric interaction. Finally, ${H}_{Z}={g}_{{\rm{e}}}{\mu }_{{\rm{B}}}\overrightarrow{B}\left({{\bf{r}}}_{i}\right)\overrightarrow{{S}_{i}}$ is the Zeeman coupling between spins and the externally applied magnetic field $\vec{B}({{\bf{r}}}_{i})$. ${\vec{S}}_{i}={({S}_{i}^{x},{S}_{i}^{y},{S}_{i}^{z})}^{t}$ are spin-angular momentum operators, g_e is the electron g-factor and μ_B is the Bohr magneton. The spin operators satisfy $[{S}_{i}^{\alpha },{S}_{j}^{\beta }]$ = ${\rm{i}}{\delta }_{{ij}}\,{\epsilon }_{\alpha \beta \gamma }{S}^{\gamma }$ with α, β, γ = x, y, z. The exchange length λ_ex defines the characteristic length scale of the magnetic textures under study.This amounts to λ_ex ∼ 5 nm typically, being much smaller than the dimensions of the nanodiscs where they are stabilized. Under these conditions, quantum fluctuations can be neglected and the equations of motion for the average magnetic moment, i.e. ${\vec{\mu }}_{j}={g}_{{\rm{e}}}{\mu }_{{\rm{B}}}\langle {\vec{S}}_{j}\rangle$, can be casted in the Landau–Lifshitz-Gilbert (LLG) equation [27, 28]:

$$\begin{eqnarray}&&\frac{\partial {\vec{\mu }}_{j}}{\partial t}=\frac{{\gamma }_{e}}{1+{\alpha }_{\mathrm{LLG}}^{2}}\left\{{\vec{\mu }}_{j}\times {\vec{B}}_{\mathrm{eff}}+\frac{{\alpha }_{\mathrm{LLG}}}{{\mu }_{{\rm{B}}}}\left[{\vec{\mu }}_{j}\times ({\vec{\mu }}_{j}\times {\vec{B}}_{\mathrm{eff}})\right]\right\}.\end{eqnarray} \tag{ 2 }$$

Here, ${\vec{\mu }}_{j}=({\mu }_{j}^{x},{\mu }_{j}^{y},{\mu }_{j}^{z})$, γ_e is the electron gyromagnetic ratio, α_LLG is the dimensionless Gilbert damping of the material, and ${\vec{B}}_{\mathrm{eff}}$ is an effective field that accounts for all terms considered in the spin Hamiltonian (1).

Depending on the magnetic object under study (material, size and shape) and the initial conditions (or the magnetic history), relaxing equation (2) leads to the stabilisation of different magnetic textures. These might have very interesting properties such as topological protection. Two archetypal textures are the magnetic vortex and the skyrmion. The former is the ground magnetic state in flat micro- or nanoscopic magnetic discs with negligible magnetocrystalline anisotropy (i.e. H_K = 0 and H_DMI = 0 in the spin Hamiltonian (1)). The minimisation of surface magnetic charges makes spins lie preferentially parallel to the borders of the disc, leading to the characteristic vortex spin curling in clockwise or counterclockwise fashion, as shown in figure 1(a). In the vortex core, spins turn out-of-plane pointing up or down defining the vortex polarity [29]. This yields four possible states that are, in principle, degenerate. The lowest energy excitation peculiar to vortices is the gyrotropic mode, in which the vortex core gyrates around its equilibrium position. It is predominantly excited by means of an in-plane oscillating magnetic field.

Zoom In Zoom Out Reset image size

Skyrmion stabilization is a little more tricky, since it requires the presence of terms favouring non-collinear spin ordering such as the DMI (i.e. ${H}_{\mathrm{DMI}}\ne 0$ in the spin Hamiltonian (1)). Néel-like skyrmions can be typically stabilized in confined geometries such as flat micro- or nanodiscs by means of the interfacial-DMI. This interaction arises in multilayers where ultra-thin ferromagnets are combined with materials having large spin–orbit coupling [24]. Additionally, a preferred out-of-plane magnetic ordering is required. The latter is usually achieved by the application of external magnetic fields, or by using materials with perpendicular magnetic anisotropy (i.e. ${H}_{{\rm{K}}}\ne 0$ in the spin Hamiltonian (1)). Skyrmions are characterised by a central core pointing in the opposite direction to the surrounding magnetization, so that spins can be projected once onto the unit sphere (see figure 1(b)). Using perpendicular (out-of-plane) oscillating fields it is possible to excite the breathing mode peculiar to the skyrmion, in which it conserves radial symmetry [30].

The objective of this work is to study the coupling between vortex or skyrmion excitations and photons. In particular, photons generated by currents in superconducting circuits, either propagating in transmission lines or stationary photons living in cavities. The total Hamiltonian, including the magnetic nanodisc, the electromagnetic field and their interaction is decomposed as:

$$\begin{eqnarray}&&{H}_{T}={H}_{S}+{H}_{Q}+{H}_{I}.\end{eqnarray} \tag{ 3 }$$

Here, H_S is the spin Hamiltonian (see equation (1)) and H_Q is the photonic Hamiltonian. In the case of one dimensional (1D) transmission lines it reads [31];

$$\begin{eqnarray}&&{H}_{Q}={\hslash }\int {\rm{d}}\omega \,\omega \,{r}_{\omega }^{\dagger }{r}_{\omega }+{\hslash }\int {\rm{d}}\omega \,\omega \,{l}_{\omega }^{\dagger }{l}_{\omega },\end{eqnarray} \tag{ 4 }$$

where $[{r}_{\omega },{r}_{\omega ^{\prime} }^{\dagger }]=\delta (\omega -\omega ^{\prime} )$ are right-moving modes while $[{l}_{\omega },{l}_{\omega ^{\prime} }^{\dagger }]=\delta (\omega -\omega ^{\prime} )$ are the left-moving ones. Here, we are considering 1D fields where the k index stands for the photon wavenumber for the transverse mode solutions of the wave equations. Owing to the fact that light propagates in 1D (z with our choice, see figure 2), we know that ω_k = vk with v the propagation velocity. In transmission lines $v=\sqrt{1/{LC}}$ with L (C) being the effective inductance (capacitance) per unit of length.

Zoom In Zoom Out Reset image size

H_I is the Zeeman coupling between the magnetic moments and the magnetic field $\vec{B}({\bf{r}})$ generated by the circuit:

$$\begin{eqnarray}&&{H}_{I}=-{g}_{{\rm{e}}}{\mu }_{{\rm{B}}}\displaystyle \sum _{i}{\vec{{ \mathcal S }}}_{i}\,\vec{B}({{\bf{r}}}_{i}).\end{eqnarray} \tag{ 5 }$$

Finally, the magnetic field at point r = (x, y, z) is quantized in the Coulomb gauge as [32, ch 10]:

$$\begin{eqnarray}&&\vec{B}({\bf{r}})=\sqrt{\displaystyle \frac{{\hslash }}{2{{vc}}^{2}{\epsilon }_{0}}}\left(\begin{array}{c}{{ \mathcal B }}_{x}(x,y)\\ {{ \mathcal B }}_{y}(x,y)\\ 0\end{array}\right)\int {\rm{d}}\omega \sqrt{\omega }\left({r}_{\omega }^{\dagger }{{\rm{e}}}^{{\rm{i}}\omega z/v}+{l}_{\omega }^{\dagger }{{\rm{e}}}^{-{\rm{i}}\omega z/v}+{\rm{h}}.{\rm{c}}.\right).\end{eqnarray} \tag{ 6 }$$

Here, h.c. means hermitean conjugate, ${{ \mathcal B }}_{x}$ (${{ \mathcal B }}_{y}$) are the transverse x (y) components of the field density having units of inverse length, c is the speed of light in vacuum and ₀ is the vacuum permittivity. We emphasise that we are mainly interested in discussing the single photon coupling. Thus, these components must be understood as the root mean square value of the field (${\vec{b}}_{0}$) generated by the vacuum current fluctuations in the circuit (i_rms). In section 3 we explain how to compute them. Experimentally, the single photon regime can be achieved at very low temperatures (i.e. few mK) and upon severe attenuation of the input/output signal lines.

Quantum systems coupled to 1D bosonic fields, as in equations (5) and (6),are of high interest. They are named waveguide quantum electrodynamics (QED) to differentiate them from cavity-QED setups, where the quantum system is coupled to the discrete stationary modes of a cavity. 1D waveguides are used to enhance the light and matter coupling in order to, e.g. mediate effective interactions between spatially separated quantum systems or to generate or manipulate quantum states of light as single or N-photon wavepackets. Here, we use this waveguide QED setup to perform broadband spectroscopy of magnetic textures.

An ideal experimental setup for resolving the excitation spectrum of magnetic nanodiscs will use the field fluctuations in a waveguide b₀ as the perturbation field, as sketched in figure 2(a). Owing to the continuum spectrum of the waveguide (using superconducting circuits, the operation range can be safely assumed to lie between few MHz up to 15–20 GHz) we can access all characteristic modes in a single experiment. In particular, the excitation spectrum can be obtained by measuring the transmission through the superconducting waveguide. To see how it works, we study the dynamics of a nanodisc coupled to a waveguide through equation (5). We highlight that we will consider two locations, i.e. position #1 and position #2 in figure 2(a). The former will be used to excite modes susceptible to in-plane magnetic fields. On the other hand, position #2 will be used to excite modes susceptible to an out-of-plane excitation. In both cases, the rms field is negligible along the z direction (${b}_{0}^{z}\sim 0$).

The field fluctuations on the waveguide can be estimated numerically in the following way. We assume that a total current i_rms is circulating through the superconducting transmission line with given width w and thickness. The spatial distribution of i_rms is calculated using the software 3D-MLSI that allows one to solve the London equations in thin film superconducting wires [33]. We then use the same software to compute the spatial distribution of ${\vec{b}}_{0}$ generated by i_rms. As an example, we use i_rms = 11 nA and assume a Nb central transmission line with a constriction of width w = 500 nm and thickness 150 nm (see section 3.1 for more details about the constriction). Figures 2(b) and (c) show the resulting density plots of b₀^x and ${b}_{0}^{y}$ along the z = 0 and y = 0 planes corresponding to positions #1 and #2, respectively. As it can be seen, the excitation magnetic field $\vec{B}({\bf{r}})={\vec{b}}_{0}$ is non-homogeneous along the surface of the disc located at both locations #1 and #2. However, to simplify calculations we will approximate $\vec{B}({\bf{r}})$ to the one passing through the disc centre (core of the topological solution), i.e. $\vec{B}({\bf{r}})\cong \vec{B}({{\bf{r}}}_{c})$, with r_c the centre of the nanodisc. As we will see in section 3.1, this is a good approximation for the first resonant modes of both vortices and skyrmions. For a disc located at position #1 (#2), $\vec{B}({{\bf{r}}}_{c})$ has only x (y) component. Accordingly, $\vec{B}({\bf{r}})\cong {B}^{x}({{\bf{r}}}_{c})={b}_{0}^{x}({{\bf{r}}}_{c})$ at position #1 whereas $\vec{B}({\bf{r}})\cong {B}^{y}({{\bf{r}}}_{c})={b}_{0}^{y}({{\bf{r}}}_{c})$ at position #2. Consequently, we can write a simplified total coupling Hamiltonian

$$\begin{eqnarray}&&{H}_{T}={H}_{S}+{H}_{Q}+\displaystyle \sum _{j}{S}_{j}^{\alpha }\,\int {\rm{d}}\omega \,{\lambda }_{\alpha }(\omega )\,\left({{\rm{e}}}^{{\rm{i}}\omega z/v}{r}_{\omega }^{\dagger }+{{\rm{e}}}^{-{\rm{i}}\omega z/v}{l}_{\omega }^{\dagger }+{\rm{h}}.{\rm{c}}.\right),\end{eqnarray} \tag{ 7 }$$

with ${\lambda }_{\alpha }(\omega )=\sqrt{\tfrac{{\hslash }}{2{{vc}}^{2}{\epsilon }_{0}}}{{ \mathcal B }}_{\alpha }({x}_{c},{y}_{c})\sqrt{{\omega }_{k}}$ and α = x and y for positions #1 and #2, respectively.

So far, we have just written the Zeeman coupling (5) in a convenient way by using the quantized field in the waveguide, see equation (6). This is rather general, but we are interested in the coupling to magnetic textures and, in particular, in resolving their excitation spectrum in a transmission experiment. This facilitates the job, since it is known that using linear response theory³ and input–output theory it is possible to find the scattering matrix. For the transmission function this results in [34]:

$$\begin{eqnarray}&&T(\omega )=1-\displaystyle \sum _{i}\displaystyle \frac{{\rm{\Gamma }}({\omega }_{n})}{{\rm{\Gamma }}({\omega }_{n})+{\rm{\Delta }}{\omega }_{n}+{\rm{i}}({\omega }_{n}-\omega )},\end{eqnarray} \tag{ 8 }$$

while the reflection is given by R(ω) = T(ω) − 1. The transmitted (reflected) measured power is given by $| T(\omega ){| }^{2}$ ($| R(\omega ){| }^{2}$). Here, ω_n are the resonant frequencies of the different normal modes (n) of the magnetic texture whereas Δω_n are their corresponding linewidths (being proportional to α_LLG). Γ (ω_n) is the emission rate of the texture in the waveguide which depends on λ_α (ω). Its exact calculation is based, in principle, in knowing the analytic solution for the magnetic textures, which is not possible in general. However, it is always possible to compute it numerically as we will explain in section 3.

In this section we compute $| T(\omega )|$ numerically for a model soft magnetic disc in the absence of both DMI and magnetocrystalline anisotropy. This is to say, we restrict ourselves to the study of the normal modes peculiar to magnetic vortices. Simulations are performed using the GPU-accelerated code MuMax³ [35]. We consider a disc with radius r = 500 nm and thickness t = 80 nm. Finally, for a more general discussion, we assume prototype material parameters.

The saturation magnetization of conventional ferromagnets, e.g. Fe or Ni₈₀Fe₂₀ (Py), amounts to M_s ∼ 1.7 MA m⁻¹ and M_s ∼ 0.8 MA m⁻¹, respectively [36]. Other interesting candidate is YIG with a very reduced M_s ∼ 0.14 MA m⁻¹ [37]. Finally, much attention is currently being pay to ferromagnetic composites like Heusler alloys or Co_xFe_1−x. For example, NiMnSb exhibits M_s ∼ 0.85 MA m⁻¹ whereas the saturation magnetization of Co_xFe_1−x amounts to M_s ∼ 2.4 MA m⁻¹ for x = 0.25 [38, 39]. Typical values of the exchange stiffness constant are A(Fe) = 21 pJ m⁻¹, A(Py) = 13 pJ m⁻¹, A(YIG) = 1.9 pJ m⁻¹, A(NiMnSb) = 9.0 pJ m⁻¹ and A(Co_0.25Fe_0.75) = 26 pJ m⁻¹. Finally, very low α_LLG values can be currently achieved using low damping ferromagnets like the above mentioned Heusler and Co_xFe_1−x alloys reaching 10⁻³ and 5 × 10⁻⁴, respectively [38, 39]. Typical ferromagnetic metals such as Fe and Py exhibit larger α_LLG of the order of 2 × 10⁻³ and 8 × 10⁻³, respectively [36]. Currently, record low damping values are reported for the insulating ferrimagnet YIG, having α_LLG ∼ 5 × 10⁻⁵, although at the cost of a largely reduced saturation magnetization [37]. With this considerations in mind we set M_s = 1 MA m⁻¹, A_ex = 15 pJ m⁻¹ and α_LLG = 10⁻³. These values yield a total exchange length ${\lambda }_{\mathrm{ex}}=\sqrt{2{A}_{\mathrm{ex}}/{\mu }_{0}{M}_{{\rm{s}}}}\sim 5\,\mathrm{nm}$, typical of conventional ferromagnets.

The excitation spectrum of the prototype ferromagnetic disc is obtained by applying a perturvative field ${B}^{\alpha }(\tau )=A\,\mathrm{sinc}({\omega }_{\mathrm{cutoff}}\tau )$ with amplitude A = 1 mT and ${\omega }_{\mathrm{cutoff}}/2\pi =50\,\mathrm{GHz}$. Calculating the Fast Fourier Transform (FFT) of the resulting spatially-averaged time-varying magnetization along the field direction results in the spectrum plotted in figure 3(a) (bottom panel). The black curve results when applying an in-plane driving field (i.e. α = x, position #1) whereas the red one corresponds to an out-of-plane excitation (i.e. α = y, position #2). Each peak corresponds to a normal mode of the ferromagnetic disc labelled g1, g2, a1, a2, a3 and r. Fitting each resonance peak to a Lorentzian function allows estimating ω_n and Δω_n with n = {g1, g2, a1, a2, a3 and r}. These values, together with the estimated Γ(ω_n) computed in section 3.1, are plugged into equation (8) resulting in the transmission curve plotted in figure 3(a) (top panel).

Zoom In Zoom Out Reset image size

We calculate now the spatial distribution of the amplitude of the FFT of the spatially-resolved time-varying magnetization. The resulting plots can be seen in figure 3(b) for the upper and bottom disc surfaces and the transverse cut. This representation helps to visualise the magnetization profile of each mode that we describe in the following. The lowest energy mode (g1) corresponds to the vortex core translation around the central position. The gyration sense is only given by the vortex polarity regardless the direction of the in-plane magnetization which plays no role. This is the archetypal gyrotropic mode where the vortex core is only minimally distorted through the disc thickness, as it can be seen in the transverse cut shown in figure 3(b). At slightly larger energy, we observe the second gyrotropic mode (g2) where the vortex gyrates in the upper/bottom disc surface with a π phase shift. In this way, the vortex core itself flexes through the disc thickness yielding one node in the centre. This node can be easily appreciated in the transverse cut shown in figure 3(b). Several azimuthal modes (a1, a2 and a3) can be observed as well. Here, the vortex core also gyrates but there is an additional magnetization spiral that shows up in the upper/bottom surface plots in figure 3(b). This spiral breaks radial symmetry and exhibits a increasing number of nodes along the radial direction for increasing mode number. Importantly, these modes are not homogeneous trough the disc thickness but curl at opposite directions from top to the bottom surfaces [40]. Finally, in the case of an out-of-plane excitation field one finds a dominant mode conserving radial symmetry (r). As shown in figure 3(b), this mode is neither homogeneous through the disc thickness.

The energy of each mode can be slightly tuned by means of an external DC magnetic field B_DC. This can be seen in figure 3(c), where the energy spectrum is plotted against B_DC. Curves are calculated for an out-of-plane field applied along the same direction as the vortex core magnetization. The gyrotropic mode increases linearly with the applied magnetic field whereas higher energy modes exhibit a non-monotonous behaviour. The presence of other low-amplitude resonant modes can be appreciated as well. This is the case of, e.g. modes visible between 12 and 14 GHz that will not be studied here.

To obtain the coupling g_n we proceed as follows. We fix our attention to a half-wavelength CPW resonator, where a central superconducting line of a few microns is surround by two ground planes. The geometry of this experiment is similar to that sketched in figure 2(a) for the case of a transmission line. In the present case, two capacitors are used to confine the photons inside the resonator yielding different quantized modes with frequency ω_C = ω_n. First, we need to compute the rms of the zero point current fluctuations flowing through the central conductor:

$$\begin{eqnarray}&&{i}_{\mathrm{rms}}={\omega }_{{\rm{C}}}\sqrt{\displaystyle \frac{{\hslash }\pi }{2{Z}_{0}}},\end{eqnarray} \tag{ 13 }$$

with Z₀ = 50 Ω the typical impedance of the resonator.

The spatial distribution of i_rms is calculated using the software 3D-MLSI as described in section 2.3. This current generates a field $\vec{B}({\bf{r}})$ that is used to compute b₀ at the centre of the nanodisc. Besides, $\vec{B}({\bf{r}})$ also serves to excite the topological solution. Therefore, we use a perturvative field B^α(τ) = ${b}_{0}^{\alpha }$ (r_c) sin (ωτ) and compute the resulting maximum response in magnetization ${\rm{\Delta }}{M}_{n}^{\alpha }$ for each mode n. Here, α = x for position #1 whereas α = y for position #2. Approximating the non-homogeneous rms field to the one produced at the centre of the nanodisc is a good approximation for the gyrotropic, first azimuthal modes and the breathing mode peculiar to skyrmions. This is so since the magnetization modulation concentrates mostly on the disc's central region (see figure 3(b)). This is not the case of higher-order azimuthal and radial modes. In these cases, however, our approach underestimates the resulting ΔM_n by a factor ∼1.5 only, as checked numerically. The coupling is finally calculated using equation (11). Finally, it remains to compute the rate Γ (ω_n) appearing in the transmission equation (8). This can be done using the relation between the waveguide-resonator coupling [46]:

$$\begin{eqnarray}&&{g}_{n}=\sqrt{\displaystyle \frac{{\rm{\Gamma }}({\omega }_{n}){\omega }_{n}}{\pi }}.\end{eqnarray} \tag{ 14 }$$

In the simulations, we assume that a constriction of width w is patterned in the central conductor to increase the strength of the resulting b₀ field. The reader is referred to figure 2 for the definition of w and to [47] and [19] for a discussion on the increase of the coupling strength using nanoconstrictions. In our calculations, w = 500 nm (in a 150 nm thick Nb central conductor) for the magnetic vortex whereas w = 50 nm (in a 50 nm thick Nb conductor) for the skyrmion.

In the case of the vortex, we consider a disc with radius r = 500 nm and thickness t = 80 nm and prototype material parameters, i.e. M_s = 1 MA m⁻¹ and A_ex = 15 pJ m⁻¹. We let the damping vary between 10⁻⁴ < α_LLG < 10⁻². We highlight that coupling to each normal mode is calculated assuming a different CPW resonator with characteristic frequency ω_C = ω_n, with n = {g1, g2, a1, a2, a3, r and s}. In figure 4 we plot the resulting strong coupling condition given by equation (12) for the different resonant modes and for different values of the damping parameter. The coupling condition depends linearly on α_LLG, increasing for decreasing damping parameter. Interestingly, although ΔM_n decreases for high order modes, all of them exhibit similar values of the coupling condition, increasing slightly for increasing frequency. This is due to the fact that the reduced ΔM_n is compensated by an increased strength of the zero-point current fluctuations at larger frequencies, see equation (13). Additionally, Δω_n also increases for increasing mode number as Δω_n ∝ α_LLGω_n. Finally, we highlight that the coupling strength is independent of the vortex polarity (up/down) and circulation (clockwise/counterclockwise).

Zoom In Zoom Out Reset image size

For the disc size analysed here, using materials with low α_LLG ∼ 10⁻³ will assure reaching the strong coupling regime for all resonant modes. This is well within the current state-of-the-art using low-damping ferromagnets like Heusler alloys or CoFe. Using the latter material, even larger coupling factors are expected stemming from its large saturation magnetisation. Moreover, the reported value of α_LLG ∼ 5 × 10⁻⁴ was measured for a very thin 10 nm thick Co_0.25Fe_0.75 film at room-temperature [39]. This means that even lower damping parameters could result in the proposed ferromagnetic disc (t = 80 nm) and working at very low temperatures (mK range).

These numbers can be compared with the expected strong coupling condition for an uniformly magnetised thin-film. For this purpose, we assume a 5 nm thick 1 μm × 13 μm magnet, which yield the same volume as that of the disc used for vortex stabilisation. We also assume the same material parameters, i.e. M_s = 1 MA m⁻¹, A_ex = 15 pJ m⁻¹ and α_LLG = 10⁻³. The thin-film lies on top of a w = 1 μm constriction on a Nb conduction line (with thickness 150 nm). The resonator produces a relatively homogeneous magnetic field b₀ along the x direction (as for position #1, see figure 1(a)). Finally, an external bias magnetic field B_DC is applied along the z direction (long side of the magnet). Under these circumstances the uniform Kittel mode will be excited in the thin-film, the frequency of which is related to B_DC as ${\omega }_{\mathrm{Kittel}}\cong {\gamma }_{e}\sqrt{{B}_{\mathrm{DC}}({B}_{\mathrm{DC}}+{\mu }_{0}{M}_{{\rm{s}}})}$ [48]. Numerical simulations yield $4g/{\rm{\Delta }}{\omega }_{\mathrm{Kittel}}\sim 3$, depending very smoothly on ω_Kittel. This is of the same order of magnitude as $4g/{\rm{\Delta }}{\omega }_{n}\cong \{0.95,0.85,1.0,1.1,2.1,2.3\}$ obtained for n = {g1, g2, a1, a2, a3 and r}, respectively, when assuming α_LLG = 10⁻³. We highlight that excitation of the Kittel mode involves precession of all magnetic moments at unison. In contrast, the gyrotropic and first azimuthal modes cause a sizeable magnetization modulation in a very small region of the nanomagnet (disc's centre).

Regarding the skyrmion, we will focus on thin-film nanodiscs with perpendicular anisotropy and interfacial-DMI. Under these circumstances, a Néel skyrmion can be stabilized in the absence of an external bias magnetic field [49–52]. We will assume a thin multilayer disc with r = 50 nm and t = 1 nm. Material parameters given in the literature for systems allowing skyrmion stabilisation can vary significantly. For example, M_s = 0.6 MA m⁻¹, A_ex = 15 pJ m⁻¹ and perpendicular magnetic anisotropy constant K_u = 0.8 MJ m⁻³ or M_s = 1.4 MA m⁻¹, A_ex = 27 pJ m⁻¹ and K_u = 1.4 MJ m⁻³ are used for Pt/Co bilayers [49, 50]. Regarding the interfacial DMI strength D_i, it usually varies between 2 mJ m⁻² < D_i < 6 mJ m⁻². In our simulations we use M_s = 1 MA m⁻¹, A_ex = 15 pJ m⁻¹, K_u = 1 MJ m⁻³ and D_i = 3 mJ m⁻². Finally, in order to keep consistency with numerical results obtained for the vortex, we let the damping vary between 10⁻⁴ < α_LLG < 10⁻².

Reaching strong coupling with the skyrmion is much more difficult. The strong coupling condition for the breathing mode peculiar to skyrmions lies almost one order of magnitude below that for the vortex (see figure 4). This is mostly due to the smallness of the disc used in the simulations. Ultra-thin films are required to obtain sizeable perpendicular anisotropy constants and interfacial DMI, limiting enormously the maximum thickness of the structures with skyrmionic ground states (t ∼ 1 nm). In addition to that, ultra-thin films usually yield damping parameters in the α_LLG ∼ 0.1 range. With these considerations in mind we can conclude that reaching the strong coupling condition using skyrmions in state-of-the-art thin film multilayers would require strongly increasing the intensity of the zero-point field fluctuations in the cavity. Increasing b₀ can be achieved by patterning even narrower constrictions or by using low-impedance cavities [47]. A different possibility would imply using low-damping materials for skyrmion stabilisation. For instance, promising values of ∼3 × 10⁻⁴ have been recently reported in bismuth doped YIG showing large perpendicular anisotropy [53].