\section{Center-of-mass motion}

The motion of the center of mass $(X,Y)$ of a magnetic bubble is conveniently represented by a complex variable $r_1 = (X-iY)/2$. The general motion of a bubble with inertial mass $\mathcal M$ and a gyrotropic coefficient $\mathcal G$ in a parabolic potential with stiffness $\mathcal K$ is a superposition of two spiral motions: 
\begin{equation}
r_1(t) = \sum_{m = \pm 1} a_{m} e^{i\omega_{m} t - \Gamma_{m} t}. 
\label{eq:r1-t}
\end{equation}
For weak dissipation, $\Gamma_{m} \ll \omega_{m}$, 
\begin{equation}
\omega_{\pm 1} \approx -\frac{\mathcal G}{2 \mathcal M} 
	\pm \sqrt{\left( \frac{\mathcal G}{2 \mathcal M} \right)^2 + \frac{\mathcal K}{\mathcal M}},
\quad
\frac{\Gamma_{+1}}{\Gamma_{-1}} \approx \left| \frac{\omega_{+1}}{\omega_{-1}} \right|.
\end{equation}
The top panel of Fig.~\ref{fig:COM} shows the best fit to Eq.~(\ref{eq:r1-t}) with $\omega_{+1}/2\pi = 0.97$ GHz, $\Gamma_{+1} = 0.25$ ns$^{-1}$, $\omega_{-1}/2\pi = -4.27$ GHz, and $\Gamma_{-1} = 1.0$ ns$^{-1}$. 

\begin{figure}[b]
\includegraphics[width=0.95\columnwidth]{figS1}
\caption{Top panel: The displacement of the bubble's center of mass $(X,Y)$ as a function of time (points) and the best fit to the two-mode approximation (\ref{eq:r1-t}) (lines). Bottom panel: The difference between the data and the best fit (points) and the best fit to a superposition of two spiral motions with frequencies $\omega_{+1} \pm \omega_0$ (lines). The vertical scales in the two panels differ by a factor of 30.}
\label{fig:COM}
\end{figure}

A small systematic deviation between the best-fit line and the data is plotted in the bottom panel of Fig.~\ref{fig:COM}. Much of it can be accounted for by a superposition of two spiral motions with frequencies $\omega_{+1} \pm \omega_0$, where $\omega_0/2\pi = 3.15$ GHz is the eigenfrequency of the bubble's breathing mode. This is likely the effect of an anharmonic coupling between the two modes. 

We thus conclude that any unidentified modes contribute no more than 0.1 nm in amplitude, or about 1 per cent, to the center-of-mass motion. 

\section{Locating the domain wall}

We define the domain wall as a line of points in the plane of the film where out-of-plane magnetization $M_z(x,y)$ vanishes. Because our numerical simulations were done on a discrete lattice, we needed an algorithm to extract a continuous line from discrete data points. The location of the domain wall was determined in two steps: we first identified a discrete set of points with $M_z=0$; we then fit their positions to a line, Eq.~(5) in the main text.

\begin{figure}
\includegraphics[width = 0.6\columnwidth]{figS2}
\caption{Tracing a domain wall (solid line).  Circles denote the centers of lattice cells in the simulation. Squares are points at cell boundaries (dashed lines).  Crosses are locations of the domain wall obtained by linear extrapolation. Open red symbols mean $M_z>0$, filled blue symbols mean $M_z<0$.}
\label{find_boundary}
\end{figure}

The first step is illustrated in Fig.~\ref{find_boundary}. The magnetization at the boundary of a cell is determined by interpolation with the magnetization at the center of the current cell and its nearest neighbors as follows:
\begin{eqnarray*}
\mathbf M_{i,j+\frac{1}{2}} & = & \left( \mathbf M_{i,j} + \mathbf M_{i,j+1} \right)/2, \\
\mathbf M_{i+\frac{1}{2},j} & = & \left( \mathbf M_{i,j} + \mathbf M_{i+1,j} \right)/2, \\
\mathbf M_{i+\frac{1}{2},j+\frac{1}{2}} & = & \left( \mathbf M_{i,j} + \mathbf M_{i,j+1} + \mathbf M_{i+1,j} + \mathbf M_{i+1,j+1} \right)/4.
\end{eqnarray*}
If the wall intersects the cell then the sign of $M_z$ at one or more points at its boundary differs from that at the center of the cell. At the crudest level, we could use the centers of intersected cells as a proxy for the location of the domain wall. To refine this result, we used a linear interpolation for $M_z$ between the center of a cell and the points on its boundary to find the locations with $M_z = 0$. The refinement reduced the discretization noise by a factor of 70, which was particularly important for modes with higher azimuthal numbers, which had small amplitudes. 

\section{Local interactions} 

Tension $\sigma$ and alignment stiffness $\kappa$ come from the local portion of the energy functional. Consider a  ferromagnetic film of thickness $t$ with two domains of magnetization, $\mathbf m = (0,0,1)$ and $\mathbf m = (0,0,-1)$, separated by a straight domain wall of length $\ell$ along the $x$-axis.  The magnetization field can be parametrized as 
\[
\mathbf m(y) = (\sin{\theta}\cos{\psi}, \sin{\theta}\sin{\psi}, \cos{\theta})
\]
with $\theta = \theta(y)$ and $\psi = \operatorname{const}$. In the state of lowest energy, in-plane magnetization points along the wall, $\psi = 0$ (Bloch wall). The energy contains three terms: exchange 
\[
U_\mathrm{exchange} = 
A t \int d^2r \, (\partial_i \mathbf m) \cdot (\partial_i \mathbf m) 
	= A t\ell \int_{-\infty}^\infty dy \left( \frac{d\theta}{dy} \right)^2, 
\]
easy-axis anisotropy  
\[
U_\mathrm{anisotropy} = 
-K t \int d^2r \,  m_z^2 = -K t\ell  \int_{-\infty}^\infty dy\, \cos^2{\theta},
\] 
and the energy of the magnetic field, which can be evaluated as the Coulomb interaction of magnetic charges at the top and bottom surfaces of the film with area densities $\mp M_z(\mathbf r)$: 
\[
U_\mathrm{dipolar} =
	\frac{\mu_0}{4\pi} \int d^2 r_1 \, d^2 r_2 \, 
V(\mathbf r_1 - \mathbf r_2) M_z(\mathbf r_1) M_z(\mathbf r_2) 
\]
with the interaction kernel 
\[
V(\mathbf r) = \frac{1}{r} - \frac{1}{\sqrt{r^2+t^2}} \equiv \mathbb{V}(r).
\]
For future convenience, we extend the definition of $\mathbb V(x)$ to negative values of the argument so that $\mathbb V(-x) = \mathbb V(x)$: 
\begin{equation}
\mathbb{V}(x) = \frac{1}{|x|} - \frac{1}{\sqrt{x^2+t^2}} 
= \lim_{\epsilon \to 0} \frac{1}{\sqrt{x^2+\epsilon^2}} - \frac{1}{\sqrt{x^2+t^2}}.
\label{eq:Vbb}
\end{equation}
The interaction kernel $V(\mathbf r)$ has a peak with a characteristic width of the order of the film thickness $t$. If magnetization varies slowly on that length scale, we may approximate 
$V(\mathbf r) \approx 2 \pi t \delta(\mathbf r)$. Then the dipolar energy assumes a local form:
\begin{equation}
U_\mathrm{dipolar} \approx t \int d^2r \, \frac{\mu_0 M_z^2}{2}
	= \frac{\mu_0M^2 t\ell}{2} \int_{-\infty}^\infty dy \, \cos^2{\theta}.
\label{eq:dipolar-local}
\end{equation}
Put another way, we make a local approximation for the magnetic field inside the film, $H_z(\mathbf r) \approx -M_z(\mathbf r)$, whose energy density $\mu_0 H^2/2 \approx \mu_0 M_z^2/2$. 

The energy cost of a domain wall is thus
\begin{equation}
U = \ell t \int_{-\infty}^{\infty} dy 
	\left[
		A {\theta'}^2 + (K - \mu_0 M^2/2) \sin^2{\theta} 
	\right].
\label{eq:energy-local}
\end{equation}
Minimization of this energy yields the domain-wall profile 
\begin{equation}
\cos{\theta(y)} = \tanh{\left( 
		\frac{y}{\lambda} \sqrt{Q-1} 
	\right)},
\end{equation}
where $\lambda = \sqrt{2A/\mu_0 M^2}$ is the exchange length and $Q = 2K/\mu_0 M^2 > 1$ is the ``quality factor.'' The minimized energy per unit length gives line tension 
\begin{equation}
\sigma = U/\ell = 2 \mu_0 M^2 t \lambda \sqrt{Q - 1}. 
\end{equation}

When in-plane magnetization deviates from the direction of the wall, it creates bulk magnetic charges with volume density $-\nabla \cdot \mathbf M = - d M_y/dy$. This induces an additional magnetic field $H_y \approx - M_y$ and thus generates additional energy density $\mu_0 H_y^2/2 \approx (\mu_0 M^2/2) \sin^2{\theta} \sin^2{\psi}$. The domain-wall energy per unit length increases to
\begin{equation}
U(\psi)/\ell = 2 \mu_0 M^2 t \lambda \sqrt{Q - 1 + \sin^2{\psi}} \sim \sigma + \kappa \psi^2/2
\end{equation}
for small angles $\psi$.  The strength of coupling between the azimuthal angle and the direction of the wall is
\begin{equation}
\kappa = \frac{1}{\ell}\left. \frac{d^2 U}{d \psi^2} \right|_{\psi=0} = \frac{\sigma}{Q-1}. 
\end{equation}

\section{Long-range interaction}

Equation (\ref{eq:dipolar-local}) captures the local part of dipolar energy, whose density is determined by the local value of magnetization. Inhomogeneities in magnetization produce a stray magnetic field that gives rise to a nonlocal component, 
\begin{equation}
U_\mathrm{stray} = \frac{\mu_0}{4\pi} \int d^2r_1 \, d^2r_2 \, 
	G(\mathbf r_1 - \mathbf r_2) \, 
	\nabla M_z(\mathbf r_1) \cdot \nabla M_z(\mathbf r_2),
\label{eq:U-gradients}
\end{equation}
where $\mathbf r = (x,y)$ and $\nabla = (\partial_x, \partial_y)$ are two-dimensional vectors. 
The kernel is defined as a solution to the equation $\nabla^2 G(\mathbf r) = -V(\mathbf r)$ and has the explicit form 
\begin{equation}
G(\mathbf r) =	\sqrt{r^2+t^2} - r - t \operatorname{arsinh}(t/r) \equiv \mathbb{G}(r).
\label{eq:G-exact}
\end{equation}
As we did for $\mathbb V(x)$, we extend the definition of $\mathbb{G}(x)$ to $x < 0$ so that $\mathbb{G}(-x) = \mathbb{G}(x)$. Note that 
\begin{equation}
\mathbb G''(x) + \mathbb G'(x)/x = - \mathbb V(x).
\label{eq:Gbb-Vbb}
\end{equation}


For an infinitely sharp domain wall, Eq.~(\ref{eq:U-gradients}) reduces to a double line integral,
\begin{equation}
U_\mathrm{stray}[\mathbf r] 
	= \frac{\mu_0 M^2}{\pi} \int d \mathbf r_1 \cdot d \mathbf r_2 \, G(\mathbf r_1 - \mathbf r_2).
\label{eq:U-nonlocal-exact}
\end{equation}
The asymptotic form of the kernel (\ref{eq:G-exact}) is $G(\mathbf r) \sim -t^2/2r$ for $r \gg t$. Hence 
\begin{equation}
U_\mathrm{stray}[\mathbf r] 
	\approx -\frac{\sigma_d}{2} \int \frac{d \mathbf r_1 \cdot d \mathbf r_2}{|\mathbf r_1 - \mathbf r_2|}
\label{eq:U-nonlocal-approx}
\end{equation}
with ``dipolar tension'' 
\begin{equation}
\sigma_d = \frac{\mu_0 M^2t^2}{\pi}. 
\label{eq:sigma-d}
\end{equation}

The simplified version of the stray-field interaction (\ref{eq:U-nonlocal-approx}) is logarithmically divergent at short distances. One way to handle the divergence is to impose a short-distance cutoff, $|\mathbf r_1 - \mathbf r_2| > b$, where $b$ is a length scale of the order of the film thickness $t$. The energy of a straight wall due to its stray field, computed with the exact expression (\ref{eq:U-nonlocal-exact}), is
\[
U_\mathrm{stray} \sim -\sigma_d \ell 
	\left[
		\ln{(2\ell/t)} + 1/2
	\right],
\]
whereas the simplified version (\ref{eq:U-nonlocal-approx}) with a cutoff yields 
\begin{equation}
U_\mathrm{stray} \sim -\sigma_d \ell 
	\left[
		\ln{(\ell/b)} - 1
	\right].
\label{eq:U-straight-cutoff}
\end{equation}
The two expressions agree if we choose the cutoff to be $b = (t/2) e^{-3/2}$.


Local contributions to the energy of the domain wall $\sigma \ell$ can be absorbed into the stray-field energy (\ref{eq:U-straight-cutoff}) at the expense of renormalizing the cutoff parameter, 
\begin{equation}
b = (t/2) e^{-3/2+\sigma/\sigma_d}.
\end{equation}


\section{Energy of transverse fluctuations}

We compute the energy of a nearly straight domain wall, $y(x) \approx 0$. The local part can be written as
\begin{equation}
U_\mathrm{loc}[y] = \int_0^\ell dx \frac{\sigma y'^2}{2} 
= \ell \int_0^\infty \frac{dk}{2\pi} \sigma k^2 |y_k|^2,
\label{eq:U-flucts-local}
\end{equation}
where $y_k = \ell^{-1/2} \int_0^\ell dx \, y(x) e^{-ikx}$ is the spatial Fourier transform of $y(x)$. 

The energy of the stray field (\ref{eq:U-nonlocal-exact}) can be expanded to $\mathcal O(y^2)$ with the aid of the following results: 
\[
d\mathbf r_1 \cdot d \mathbf r_2 = dx_1 \, dx_2 
	\left[
		1 + y'(x_1)  y'(x_2)
	\right],
\]
\[
G(\mathbf r) = \mathbb G(\sqrt{x^2+y^2}) 
	= \mathbb G(x) 
		+ \mathbb G'(x) \, y^2/2x
		+ \mathcal O(y^4).
\] 
The nonlocal portion of the energy of transverse fluctuations thus consists of two parts,
\begin{eqnarray*}
U_\mathrm{stray}[y] &=& \frac{\mu_0 M^2}{\pi} \int dx_1 \, dx_2 \, \mathbb G(x_1-x_2) y'(x_1) y'(x_2)
\\
	&+& \frac{\mu_0 M^2}{\pi} \int dx_1 \, dx_2 \, \frac{\mathbb G'(x_1-x_2)}{2(x_1-x_2)} (y_1-y_2)^2.
\end{eqnarray*}
After an integration by parts and a Fourier transform, we obtain 
\[
U_\mathrm{stray}[y] = \ell \int_0^\infty \frac{dk}{2\pi} A_k |y_k|^2,
\]
where 
\[
A_k = \frac{2 \mu_0 M^2 }{\pi} \int_{-\infty}^\infty dx \, [1-\cos{(kx)}] 
	\left[
		\mathbb G''(x) + \mathbb G'(x)/x
	\right].
\]
With the aid of Eqs.~(\ref{eq:Vbb}) and (\ref{eq:Gbb-Vbb}), we obtain
\[
A_k = -\frac{4\mu_0 M^2}{\pi} \left[
		K_0(kt) + \ln{(kt/2)} + C
	\right],
\]
where $K_0(x)$ is a modified Bessel function and $C = 0.577\ldots$ is the Euler constant. In the infrared limit, $kt \to 0$,
\[
A_k \sim \sigma_d k^2 
	\left[ 
		\ln{(kt/2)} + C - 1
	\right].
\]
Upon adding the local term (\ref{eq:U-flucts-local}), we obtain
\begin{equation}
U[y] \sim \ell \int_0^\infty \frac{dk}{2\pi} \sigma_d k^2 \ln{(ka)} |y_k|^2
\label{eq:U-nonlocal-asymptotic}
\end{equation}
with the short-distance scale 
\begin{equation}
a = (t/2)e^{C - 1 + \sigma/\sigma_d}. 
\label{eq:a}
\end{equation}

\end{document}


