# From Vertices to Vortices in magnetic nanoislands

Michael Saccone  Theoretical Division (T4), Los Alamos National Laboratory,
Los Alamos, New Mexico 87545, USA  Center for Nonlinear Studies (T-Division),
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA  Jack C.
Gartside  Blackett Laboratory, Imperial College London, London SW7 2AZ, United
Kingdom, Imperial College London, London SW7 2AZ, United Kingdom  Kilian D.
Stenning  Blackett Laboratory, Imperial College London, London SW7 2AZ, United
Kingdom, Imperial College London, London SW7 2AZ, United Kingdom  Will R.
Branford  Blackett Laboratory, Imperial College London, London SW7 2AZ, United
Kingdom, Imperial College London, London SW7 2AZ, United Kingdom  London
Centre for Nanotechnology, University College London, London WC1H 0AH, United
Kingdom  Francesco Caravelli  Theoretical Division (T4), Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

Recent studies in magnetic nanolithography show that a variety of complex
magnetic states emerge as a function of a single magnetic island’s aspect
ratio. We propose a model which, in addition to fitting experiments, predicts
magnetic states with continuous symmetry at particular aspect ratios and
reveals a duality between vortex and vertex states. Our model opens new means
of engineering novel types of artificial spin systems, and their application
to complex magnetic textures in devices and computing.

##  I Introduction

Remarkable advances in nanofabrication have led to the possibility of
engineering designer artificial magnetic systems. In particular, Artificial
Spin Ices (ASI) are arrays of magnetically frustrated nanoislands with low-
energy state degeneracy which emulate ice modelsWang  _et al._ (2006); Skjærvø
_et al._ (2020). Currently, the literature on ASI is expanding from the study
of novel phases of matter in many-body physics Farhan  _et al._ (2019), to
diverse applications including reconfigurable magnonics Gliga  _et al._
(2020), neuromorphic Saccone  _et al._ (2022) and reservoir computing Gartside
_et al._ (2022). The magnetic islands (MIs) of the original Artificial Square
Ice Wang  _et al._ (2006) had an aspect ratio typically above four, meaning
they were four times longer than wide, causing the island to magnetize in a
collinear state along the longitudinal direction (a ‘macrospin’ state).
Islands with lower aspect ratios can also support the emergence of a vortex
domain or even a double-vortex domain, both with zero net magnetization.
Understanding the emergence of vortex domains in single islands typically
requires the use of costly micromagnetic simulations Gliga  _et al._ (2020),
which numerically integrate the Landau-Lifschitz-Gilbert equation. The
emergence of vortices occurs when the aspect ratio is around two to three has
been recently used in order to amplify the number of configurations that
arrays of nanoislands can have, in particular granting an advantage for
reservoir computing Gartside  _et al._ (2022). Increasingly numerous island
states introduce an enriched microstate space and emergent metamaterial memory
dynamicsKeim  _et al._ (2019); Merrigan  _et al._ (2022). We developed a
simple yet feature-rich model to explore the emergence of vortices and guide
the ASI field in developing islands that encourage the formation of complex
and dynamic magnetic textures.
Artificial spin ices can be used to engineer phases of matter since at large
aspect ratios, the interaction between fully magnetized (macrospin)
nanoislands can be approximated using the dipolar interaction model: collinear
islands behave ferromagnetically, while orthogonal ones antiferromagnetically.
Since each island in an ASI has two possible magnetic states, one can map a
single vertex in a spin ice as an Ising model, with $s_{u}$, $s_{r}$, $s_{d}$
and $s_{l}$ the Ising states of the up, right, down and left island
respectively, with $s_{*}=\pm 1$. If we call $J_{\perp}$ the coupling for
perpendicular spins, and $J_{\|}$ for collinear ones, the Hamiltonian of a
coordination-four spin ice vertex can be written in the form
$H_{vert}=\frac{1}{2}\Big{(}J_{\perp}(s_{r}-s_{l})(s_{u}-s_{d})-J_{\|}(s_{l}s_{r}+s_{d}s_{u})\Big{)}.$

These interaction patterns can be repeated in any ASI which is a decimation of
a square ice. However, as this model fails to accomodate the existence of
multiple states in an island, here we argue that a similar but continuous
model of a vertex can be used to study the complex features of a single
island.

##  II  Model

We employ a generalized bipartite mean field model Gallo and Contucci (2008)
to multiple mean field domains. We partition our system in magnetic domains
composed of large numbers of Ising spins. For a magnetic system such a spin
ice, one can either choose the magnetic domains to be those of a fully
magnetized island, e.g. a macrospin, or those internal to an island. We assume
that our system is partitioned into $K$ subsystems composed of $N$ spins,
based on the symmetry of the problem. The first approximation is that each
subsystem is composed of perfectly aligned Ising spins, an approximation which
is justified for island’s aspect ratios above one. The spins in domain $A$
interact with the spins in the domain $B$, with a certain exchange interaction
$J_{AB}$, and among themselves with a certain effective interaction $J_{AA}$.
The graph which describes the interation between the sets is described by an
adjacency matrix ${\mathcal{A}}_{AB}$, $A=1\cdots K$. Each subset of spins is
of size $N_{1}\cdots N_{K}$. After a brief calculation we can see that, in the
spirit of the bipartite case, the Hamiltonian for a larger set of mean field
spins is given by
$H=\sum_{A=1,B=1}^{K}J_{AB}{\mathcal{A}}_{AB}\bigg{(}\sum_{i\in
A}^{N_{A}}s_{i}\bigg{)}\bigg{(}\sum_{j\in
B}^{N_{B}}s_{j}\bigg{)}+\sum_{A=1}^{K}h_{A}\bigg{(}\sum_{i\in
A}^{N_{A}}s_{i}\bigg{)}$, where $\mathcal{A}$ has only zeros and ones, while
$J_{AB}$ are coupling constants. As in the case of the bipartite mean field
model, we introduce for every set of spins a mean field parameter
$m_{A}=\frac{1}{N}\sum_{i\in A}s_{i}$, and a Lagrange multiplier $\lambda_{A}$
to enforce the constraints. Rescaling all the couplings by $J_{AB}\rightarrow
J_{AB}/\sqrt{N_{A}N_{B}}$, we then obtained a mean field partition function
via Steepest Descent. After a brief calculation, we obtain the emerging mean
field equations, given by $m_{A}=\tanh(\sum_{B}\beta Q_{AB}m_{B}+\beta h_{A})$
in terms of effective couplings $Q_{AB}=J_{AB}\mathcal{A}_{AB}$.

It is interesting at this point to note that the fixed points of these
equations are, for sufficiently low temperature, the minima of the Hamiltonian

$\displaystyle\mathcal{H}=\frac{1}{2}\vec{m}^{t}\mathbf{Q}\vec{m},$ (1)

but whether or not the minima may be reached depends on the type of dynamics.
As we show in the Supp. Mat., this model is able to encode the ground states
of various spin ice models. One important comment is that if we study spin ice
vertices, the effective matrix $\mathbf{Q}$ has the same structure as the
Hamiltonian of a vertex. The key difference is that while the vector $\vec{s}$
is composed of discrete states, the vector $\vec{m}$ is composed of continuous
magnetic states $-1\leq m_{i}\leq 1$. As such, the minimum of the Hamiltonian
of eqn. (1) can be studied using its eigenvectors and eigenvalues. When the
$m_{i}$ represent subdomains, the key advantage over micromagnetic simulations
is that stable states become a tractable eigenvector problem.

Figure 1: (a) Magnetic states of a single island and their interpretation. The
single vortex and anti-ferromagnetic order are interpreted as real 1 vortex
and 2-vortex states. Such nomenclature is consistent with the corresponding
relaxed states, with these as initial conditions, obtained from micromagnetic
simulations. (b) Ground state phase diagram from eqn. (1), with corresponding
eigenvectors and eigenvalues. The dashed yellow line represents the behavior
of the effective parameters, as a function of the aspect ratio, that fit the
experiments. (c) MFM image of isolated islands at various aspect ratios, with
identified states. The black bar is $5\mu$m. (d) Comparison between the
experimental populations of island states after AC demagnetization as a
function of the aspect ratio, fit by the populations obtained from the mean-
field theory, with relaxation given by Metropolis dynamics.

Magnetic island and eigenvalues of $4\times 4$ structured matrices. Let us now
restrict our discussion to the behavior of a single MI, which is the building
block of all artificial spin ices.We consider a partitioning of the island as
in Fig. 1 (a). We assume that the island can be divided in four subdomains
$m_{1},\cdots,m_{4}$ assumed to be oriented to the right when positive. By
symmetry considerations, we can then introduce three couplings $J_{d}$ between
the domain $1$ and $4$ and $3$ and $2$ across the diagonal, $J_{f}$ for the
domains $1$ and $2$ and $3$ and $4$, and then $J_{a}$ for $1$ and $3$, and $2$
and $4$. Since our effective spins are continuous, one can also introduce a
diagonal component $q$. While this interaction is important for the dynamics,
since by symmetry considerations $q$ is identical for all the spins, it is not
relevant in our considerations below. The interaction matrix is then given by

$\displaystyle\footnotesize{\mathbf{Q}}=\left(\begin{array}[]{cccc}0&J_{f}&J_{d}&J_{a}\\\
J_{f}&0&J_{a}&J_{d}\\\ J_{d}&J_{a}&0&J_{f}\\\ J_{a}&J_{d}&J_{f}&0\\\
\end{array}\right)$ (6)

What is remarkable is that for this matrix, the eigenvectors do not depend on
the particular values of $J_{a},J_{d},J_{f}$, although clearly the eigenvalues
do. The eigenvalues $\eta_{i}$ and eigenvectors of this matrix are
$\eta_{1}=J_{a}-J_{d}-J_{f}$, $\eta_{2}=-J_{a}+J_{d}-J_{f}$,
$\eta_{3}=-J_{a}-J_{d}+J_{f}$, $\eta_{4}=J_{a}+J_{d}+J_{f}$, with
corresponding eigenvectors, $\vec{v}_{1}=(1,-1,-1,1)$,
$\vec{v}_{2}=(-1,1,-1,1)$ $\vec{v}_{3}=(-1,-1,1,1)$, $\vec{v}_{4}=(1,1,1,1)$.
Clearly, also the states with opposite signs are valid eigenvectors.
As it turns out, if we look at the graphics of Fig.1 (a), these correspond
respectively to: 1) an antiferromagnetic order state, which we interpret as a
2-vortex state; 2) a domain-wall non-magnetized state; 3) a 1-vortex state 4)
and a fully magnetized state. To show that this nomenclature is justified, we
initialized a micromagnetic simulation in these states, and shown that they
relax to the corresponding real vortex states (see Supp. Mat.). If one
considers this a squashed square ice vertex, there is a duality between
“vertex type” states and the single island eigenvectors. In fact Type I
vertices, correspond to the 2-vortex state. Type II vertices may be
constructed from the 1-vortex and fully magnetized states. Type IV vertices,
highly energetic magnetic charge states, are non-magnetized state equivalent
of the domain wall type. The remaining Type III vertices are superpositions of
these eigenvectors, fittingly as they exist as excitations from the spin ice
ground state. This simple model incorporates immediately, only via symmetry
considerations, the experimentally observed magnetic states of a single MI and
those relevant to square ice studies. Interestingly, only four degrees of
freedom analytically capture what previously required micromagnetics, despite
only crudely approximating the shape of vortices. Finer structures would
emerge if more degrees of freedom were allowed through the introduction of a
subdivided grid, but these alone present a compelling correspondence to
nature. Whether the underlying magnetic interactions are due to exchange,
dipolar interactions, or a combination of the two, is irrelevant in this model
insofar as the appropriate behaviors are captured. Despite the fact that the
entire field of ASI requires this feature, it is often given for granted.
Focusing on a single island, what is then important is which eigenvalue is the
minimum one as a function of the parameters $J_{f},J_{a},J_{d}$. Since we know
that one of the states has to be fully magnetized, we can safely set $J_{f}<0$
and study the minimum eigenvalue problem as a function of $J_{a}$ and $J_{d}$
only. The result is shown in the diagram of Fig.1 (b). In the top right region
$J_{a}>0,J_{d}>0$, we have as the ground state the 1-vortex state. In the top
left region $J_{a}<J_{f}$, $J_{d}>J_{f}$ we have the two vortex states. The
bottom left of the diagram is given by fully magnetized states, while the
bottom right by a non-magnetized state, which is however rarely observed. The
interface lines between two regions $i,j$ imply that we have eigenvalue
degeneracy in the ground state. Thus, on those particular lines $J_{a}=\pm
J_{f}$, or $J_{d}=\pm J_{f}$, we have that the magnetic ground state can be
written in a continuous superposition of the form
$\vec{v}_{ij}(\theta)=\frac{\cos\theta}{2}\vec{v}_{i}+\frac{\sin\theta}{2}\vec{v}_{j}$
(U(1) symmetry). At the interface between three regions $i,j,k$, the
degeneracy becomes triple, and the corresponding eigenvector at that
particular point can be written in terms of Euler angles as
$\vec{v}_{ijk}(\theta,\phi)=\frac{\cos\theta\cos\phi}{2}\vec{v}_{i}+\frac{\sin\theta\cos\phi}{2}\vec{v}_{j}+\frac{\sin\phi}{2}\vec{v}_{k}$
(O(3) symmetry). The dashed black color lines of Fig.1 (b) are due to the fact
that the ground state becomes an excited state, which is then degenerate. The
semi-magnetized states shown in Fig. 1 (a) are never observed in experiments.
One way to explain this phenomenon is that these can only be obtained as a
superposition of all the eigenvectors $v_{1},\cdots,v_{4}$. However, in the
ground state phase diagram, there is no point of contact between the four
eigenvectors, so these states should only be seen at higher temperatures.

##  III Experiments

Let us now use the model above to intepret the experimental data on islands
with various aspect ratio. Permalloy nanomagnets of thickness 25 nm were
fabricated via electron beam lithography with $1\mu$m inter-island spacing
such that dipolar interaction between neighbouring islands is negligible; in
total, we studied 1350 nanoislands with aspect ratios of 0.5-3.0. These were
effectively annealed via an AC demagnetisation protocol, and the resulting
states were imaged via magnetic force microscopy (MFM). We then counted the
different occurrence of magnetic texture states (macrospin, vortex and double-
vortex). A subset of the islands imaged are shown in Fig. 1 c). As we can see,
as a function of the aspect ratio we only observe states given by a macrospin
and single and double vortices Ehrmann and Blachowicz (2019). Thus, in the
phase diagram of Fig. 1 (b) we must be in a region around $J_{a}\approx-J_{f}$
and $J_{d}\approx J_{f}$. Because the nanomagnets are athermal at room
temperature (the energy barrier between states is much higher than the thermal
energy at 300 K), the system can be trapped in a metastable state. To
reproduce these results we introduce a form of dynamics for the Hamiltonian of
eqn. (1). Note that for $r=1$ the majority of the population is in the
1-vortex state, and for $r\approx 2$ we are predominantly in the macrospin
state. At intermediate aspect ratios, we observe an increasing population of
2-vortex states.
As a result, we hypothesize that parametrically, we are moving along a curve
$J_{a}=-E_{a}J_{f}(r+r_{a})$, $J_{d}\propto E_{d}\frac{J_{f}^{2}}{J_{a}^{2}}$.
If $E_{0}=J_{f}$, then at $J_{a}=J_{f}$ we are at the triple point, which
allows the coexistence of all states. Such function is shown in Fig. 1 (b),
which is the yellow dashed curve. A sample of the experimental results is
shown in Fig. 1 (c), and the vertex populations for different aspect ratios
are shown in Fig. 1 (d), averaged over 1350 experiments, with light lines. The
theoretically obtained populations after a Metropolis annealing are shown in
the same figure (averaged over 100 samples per point), and are given by hard
dashed lines. We see a good match between the experimental and the theoretical
results, with fitting parameters $E_{d}=0.6$, $E_{a}=0.6$, $r_{a}=1.2$, all of
order one.

##  IV Conclusions

We provided a simple model of the emergent states of magnetic nanoislands, as
a function of its aspect ratio, showing that these are the eigenvectors of a
$4\times 4$ matrix representing the couplings. This seemingly arbitrary coarse
graining of a nanoisland not only explained previously seen vortices, but
pointed to a new, parametrically robust, two vortex state that we then
discovered experimentally. Simultaneously, this justifies the lack of domain
wall states and opens the possibility of observing degenerate continuous
ground states. By simple symmetry considerations, we attained a mapping
between ASI vertex states and MIs states that generated new insight into why
magnetic domains form. As shown in Supp. Mat., our model also reproduces the
low energy regime of known spin ices. While micromagnetic simulations are
sufficiently fast for understanding properties of 1-10 nanomagnets before an
experiment, this mean-field model excels at handling 10-1000 nanomagnets and
grasping their non-binary behaviors. Emergent states grant higher memory
capacity, analogies to more complex models in statistical physicsSchiffer and
Nisoli (2021), and the potential to embed a wider class of computational
problems within patterned nanomagnet arrays. Additionally, this is the first
time a triple-point of macrospin, vortex and double vortex magnetic textures
has been observed. Our model also enables dynamical analysis to explain
observed emergent phenomena such as avalanchesChern  _et al._ (2014) and
coupling to magnon modesLendinez and Jungfleisch (2019); Urazhdin  _et al._
(2014).

Acknowledgements. The work of F.C. and M.S. was carried out under the NNSA of
the U.S., DoE at LANL, Contract No. DE-AC52-06NA25396 (LDRD grant -
PRD20190195). W.R.B. and J.C.G. were supported by the Leverhulme Trust
(RPG-2017-257) and the Imperial College London President’s Excellence Fund for
Frontier Research. J.C.G. was supported by the Royal Academy of Engineering
under the Research Fellowship programme. We would like to thank C. Nisoli and
P. Schiffer for various comments.

## References

  * Wang  _et al._ (2006) R. Wang, C. Nisoli, R. Freitas, J. Li, W. McConville, B. Cooley, M. Lund, N. Samarth, C. Leighton, V. Crespi, _et al._ , Nature 439, 303 (2006).
  * Skjærvø  _et al._ (2020) S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Nature Reviews Physics 2, 13 (2020).
  * Farhan  _et al._ (2019) A. Farhan, M. Saccone, C. F. Petersen, S. Dhuey, R. V. Chopdekar, Y.-L. Huang, N. Kent, Z. Chen, M. J. Alava, T. Lippert, _et al._ , Science advances 5, eaav6380 (2019).
  * Gliga  _et al._ (2020) S. Gliga, E. Iacocca, and O. G. Heinonen, APL Materials 8, 040911 (2020).
  * Saccone  _et al._ (2022) M. Saccone, F. Caravelli, K. Hofhuis, S. Parchenko, Y. A. Birkhölzer, S. Dhuey, A. Kleibert, S. Van Dijken, C. Nisoli, and A. Farhan, Nature Physics 18, 517 (2022).
  * Gartside  _et al._ (2022) J. C. Gartside, K. D. Stenning, A. Vanstone, H. H. Holder, D. M. Arroo, T. Dion, F. Caravelli, H. Kurebayashi, and W. R. Branford, Nature Nanotechnology 17, 460 (2022).
  * Keim  _et al._ (2019) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Reviews of Modern Physics 91, 035002 (2019).
  * Merrigan  _et al._ (2022) C. Merrigan, D. Shohat, C. Sirote, Y. Lahini, C. Nisoli, and Y. Shokef, arXiv preprint arXiv:2204.04000 (2022).
  * Urazhdin  _et al._ (2014) S. Urazhdin, V. E. Demidov, H. Ulrichs, T. Kendziorczyk, T. Kuhn, J. Leuthold, G. Wilde, and S. O. Demokritov, Nature Nanotechnology 9, 509 (2014).
  * Gallo and Contucci (2008) I. Gallo and P. Contucci, Mathematical Physics Electronic Journal 14 (2008).
  * Ehrmann and Blachowicz (2019) A. Ehrmann and T. Blachowicz, Journal of Magnetism and Magnetic Materials 475, 727 (2019).
  * Schiffer and Nisoli (2021) P. Schiffer and C. Nisoli, Applied Physics Letters 118, 110501 (2021).
  * Chern  _et al._ (2014) G.-W. Chern, C. Reichhardt, and C. Reichhardt, New Journal of Physics 16, 063051 (2014).
  * Lendinez and Jungfleisch (2019) S. Lendinez and M. Jungfleisch, Journal of Physics: Condensed Matter 32, 013001 (2019).
  * Note (1) Https://github.com/ww1g11/JuMag.jl.
  * Note (2) Https://math.nist.gov/oommf/.
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Glauber (1963) R. J. Glauber, Journal of Mathematical Physics 4, 294 (1963).
  * Peretto (1984) P. Peretto, Biological Cybernetics 50, 51 (1984).
  * Note (3) The change of sign compared to Peretto is because we are considering an anti-ferromagnetic system.
  * Caravelli (2020) F. Caravelli, Eur.Phys. Lett. 130 (2020).
  * Caravelli and Nisoli (2020) F. Caravelli and C. Nisoli, New Journal of Physics 22 (2020).

##  Appendix A Micromagnetic Simulation

Figure 2: Vector fields representing the steady state magnetic moments
attained from micromagnetic simulation in OOMMF on a elliptical system of
spins with aspect ratio 2. Different seeds of random numbers converge to the
a) fully-magnetized, b) 1-vortex, and c) 2-vortex states.

###  A.1 Evolution from initial conditions

We evolved from the states of the paper, e.g. respectively the fully-
magnetized, 1-vortex and the 2-vortex states, and relaxed a stadium shaped
magnetic island. We have used the package
JuMag.jl111https://github.com/ww1g11/JuMag.jl, written in Julia, in which it
is easier to initialize the system. In Fig. 3 we can see the initial states
(blue lines) versus the relaxed magnetic states. Our simulations were
performed on a grid of $100\times 50$ points, with $2nm^{2}$ cell size, using
the parameters for Permalloy. As we can see, these initial states relax into
the fully-magnetized, 1-vortex and 2-vortex states respectively (the red
lines), showing that our interpretation of the paper is correct from a
micromagnetic perspective. We have used random initial conditions for an oval
shape matrix of aspect ratio 2, for the parameter of permalloy, using the
software OOMMF 222https://math.nist.gov/oommf/, obtaining similar stable
relaxed states. These are shown in Fig. 11.

Figure 3: Evolution of the initial states (blue) from our main text into the
relaxed magnetic states for a stadium shaped island with aspect ratio 2 (50nm
$\times$ 100nm). We see how the system evolves from these states into the
fully magnetized, 1-vortex and 2-vortices states respectively. Result obtained
with the JuMag.jl package, with the parameters of permalloy.

##  Appendix B Further details on experiments

Let us now use the model above to intepret the experimental data on islands
with various aspect ratio. Permalloy nanomagnets of thickness 25 nm were
fabricated via electron beam lithography with $1\mu$m inter-island spacing
such that dipolar interaction between neighbouring islands is negligible and
islands behave as isolated. 1350 nanoislands were fabricated with lengths of
460-715 nm, widths of 225-275 nm and aspect ratios of 0.5-3.0. The islands
were effectively annealed via an AC demagnetisation protocol, initially
saturating islands at 65 mT (well-above their coercive fields) and then
oscillating down to 0 field in 0.5 mT steps. The resulting island states were
imaged via magnetic force microscopy (MFM) and the occurrence of different
magnetic texture states (macrospin, vortex and double-vortex) were counted. A
subset of the islands imaged are shown in Fig. 1 of the main text.

##  Appendix C Derivation of the mean field equations

The underlying hypothesis for the mean field equations for artificial spin ice
is that all nanoislands are composed of Curie-Weiß domains. We derive mean
field quantities representing the average magnetic moment of collections of
those domains in equilibrium and equations governing the time evolution of the
mean field moments.

####  C.0.1 Generalized case: arbitrary graph

We want to solve a generalization of the Curie-Weiß model. We have $k=1\cdots
K$, and the graph which describes the interation between the sets is described
by an adjacency matrix ${\mathcal{A}}_{AB}$, $A=1\cdots K$. Each subset of
spins is of size $N_{1}\cdots N_{K}$.

Figure 4: The Weiß domains in a magnetic nanoisland which we use in our mean
field approximation. These are aligned along the longest direction of the
island.

Here we would like to write down an Hamiltonian which describes this system,
with an external field $h_{A}$ homogeneous on every set, but nonhomogeneous
otherwise.

One way to solve it is the following. First, introduce $M_{A}=\sum_{i\in
A}s_{i}$ and $M_{B}=\sum_{j\in B}s_{j}$, and $1$ smartly into the partition
function, as

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\sum_{\\{s\\}}\int_{-N_{A}}^{N_{A}}dM_{A}\int_{-N_{B}}^{NB}dM_{B}\delta(M_{A}-\sum_{i\in
A}s_{i})\delta(M_{B}-\sum_{j\in
B}s_{j})e^{\beta(J_{A}(M_{A})^{2}+J_{AB}M_{A}M_{B}+J_{B}(M_{B}^{2})}$
$\displaystyle=$
$\displaystyle\sum_{\\{s\\}}\int_{-N_{A}}^{N_{A}}dM_{A}\int_{-N_{B}}^{NB}dM_{B}\int_{-i\infty}^{i\infty}d\lambda\int_{-i\infty}^{i\infty}d\eta
e^{\lambda(M_{A}-\sum_{i\in A}}s_{i})e^{\eta(M_{B}-\sum_{j\in
B}s_{j})}e^{\beta(J_{A}(M_{A}^{2}+J_{AB}M_{A}M_{B}+J_{B}M_{B}^{2})}$
$\displaystyle=$
$\displaystyle\int_{-N_{A}}^{N_{A}}dM_{A}\int_{-N_{B}}^{NB}dM_{B}\int_{-i\infty}^{i\infty}d\lambda\int_{-i\infty}^{i\infty}d\eta\
e^{\lambda M_{A}+\eta
M_{B}+\beta(J_{A}M_{A}^{2}+J_{AB}M_{A}M_{B}+J_{B}M_{B}^{2})}$ $\displaystyle\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdot\sum_{\\{s_{i}\in A\\}}e^{-\sum_{i\in
A}s_{i}\lambda}\sum_{\\{s_{j}\in B\\}}e^{-\sum_{j\in B}s_{j}\eta}$
$\displaystyle=$
$\displaystyle\int_{-N_{A}}^{N_{A}}dM_{A}\int_{-N_{B}}^{NB}dM_{B}\int_{-i\infty}^{i\infty}d\lambda\int_{-i\infty}^{i\infty}d\eta\
e^{\lambda M_{A}+\eta
M_{B}+\beta(J_{A}M_{A}^{2}+J_{AB}M_{A}M_{B}+J_{B}M_{B}^{2})}$ $\displaystyle\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \cdot e^{N_{A}\log 2\cosh(\lambda)+N_{B}\log
2\cosh(\eta)}\ $

We now rescale $M_{A}=N_{A}m_{a}$ and $M_{B}=N_{B}m_{B}$, and perform the
rescaling $J_{A}=\tilde{J}_{A}/2N_{A}$, $J_{B}=\tilde{J}_{B}/2N_{B}$, and
$J_{AB}=\tilde{J}_{AB}/\sqrt{N_{A}N_{B}}$. We then have

$\displaystyle Z$ $\displaystyle=$ $\displaystyle\frac{1}{\xi
N^{2}}\int_{-1}^{1}dm_{A}\int_{-1}^{1}dm_{B}\int_{-i\infty}^{i\infty}d\lambda\int_{-i\infty}^{i\infty}d\eta\
e^{N\lambda m_{A}+\xi N\eta
m_{B}+N\beta(\frac{\tilde{J}_{A}}{2}m_{A}^{2}+\frac{\tilde{J}_{AB}}{\sqrt{\xi}}m_{A}m_{B}+\xi\frac{\tilde{J}_{B}}{2}m_{B}^{2})}$
(7) $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdot e^{N\log
2\cosh(\lambda)+\xi N\log 2\cosh(\eta)}.$

We can now perform the saddle point approximation, writing $Z\approx
e^{Nf(m_{a},m_{b},\lambda,\eta)}$, and we obtain the mean field equations:

$\displaystyle\partial_{m_{A}}f=0$ $\displaystyle\rightarrow$
$\displaystyle\lambda+\beta\tilde{J}_{A}m_{A}+\beta\frac{\tilde{J}_{AB}}{\sqrt{\xi}}m_{B}=0$
(8) $\displaystyle\partial_{m_{B}}f=0$ $\displaystyle\rightarrow$
$\displaystyle\eta+\beta\xi\tilde{J}_{B}m_{B}+\beta\frac{\tilde{J}_{AB}}{\sqrt{\xi}}m_{A}=0$
(9) $\displaystyle\partial_{\lambda}f=0$ $\displaystyle\rightarrow$
$\displaystyle m_{A}+\tanh(\lambda)=0$ (10) $\displaystyle\partial_{\eta}f=0$
$\displaystyle\rightarrow$ $\displaystyle m_{B}+\tanh(\eta)=0$ (11)

which can be rewritten, in terms of $m_{A}$ and $m_{B}$, as the equations

$\displaystyle m_{A}$ $\displaystyle=$
$\displaystyle\tanh\left[\beta(\tilde{J}_{A}m_{A}+\frac{\tilde{J}_{AB}}{\sqrt{\xi}}m_{B})\right].$
$\displaystyle m_{B}$ $\displaystyle=$
$\displaystyle\tanh\left[\beta(\xi\tilde{J}_{B}m_{B}+\frac{\tilde{J}_{AB}}{\sqrt{\xi}}m_{A})\right].$
(12)

If we define $\theta_{A}=\beta\tilde{J}_{A}$,
$\theta_{B}=\beta\xi\tilde{J}_{B}$, and
$\theta_{AB}=\frac{\beta\tilde{J}_{AB}}{\sqrt{\xi}}$, the equations read

$\displaystyle m_{A}$ $\displaystyle=$
$\displaystyle\tanh\left(\theta_{A}m_{A}+\theta_{AB}m_{B}\right).$
$\displaystyle m_{B}$ $\displaystyle=$
$\displaystyle\tanh\left(\theta_{B}m_{B}+\theta_{AB}m_{A}\right).$ (13)

It is clear from the derivation above that we can extend the same procedure to
an arbitrary set of parameters for a Hamiltonian of the form

$H=\sum_{A=1,B=1}^{K}J_{AB}{\mathcal{A}}_{AB}(\sum_{i\in A}s_{i})(\sum_{j\in
B}s_{j})+\sum_{A=1}^{K}h_{A}(\sum_{i\in A}s_{i})$ (14)

where $\mathcal{A}$ has only zeros and ones, while $J_{AB}$ are coupling
constants.

As in the case of the bipartite mean field model, we introduce for every set
of spins a mean field parameter $m_{A}=\frac{1}{N}\sum_{i\in A}s_{i}$, and a
Lagrange multiplier $\lambda_{A}$ to enforce the constraints. Also, we rescale
all the couplings by $J_{AB}\rightarrow J_{AB}/\sqrt{N_{A}N_{B}}$. It is not
hard to see that the partition function can be written as

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\int_{-1}^{1}[dM]\int_{-i\infty}^{\infty}[d\lambda]e^{\beta(\sum_{A=1,B=1}^{K}\sqrt{N_{A}N_{B}}J_{AB}{\mathcal{A}}_{AB}m_{A}m_{B}}$
(15) $\displaystyle\ \ \ \ \ \ \ \cdot
e^{\sum_{A=1}^{K}N_{A}h_{A}+\sum_{A=1}^{K}N_{A}\log 2\cosh(\lambda_{A})}$

Now, assuming that $N_{A}=\xi_{A}N$ and taking the thermodynamic limit, we get
a set of $2K$ mean field equations of the form:

$\displaystyle m_{A}+\tanh(\lambda_{A})=0$
$\displaystyle\lambda_{A}+\beta\sum_{B}{\mathcal{A}}_{AB}Q_{AB}m_{B}+h_{A}=0$
(16)

where $Q_{AB}=\frac{\tilde{J}_{AB}}{\xi_{A}\xi_{B}}$. As a result, we have

$m_{A}=\tanh(\sum_{B}{\mathcal{A}}_{AB}\beta Q_{AB}m_{B}+\beta h_{A}).$ (17)

In the continuous limit of an infinite graph the index $A$ can become a
continuous variable. We have thus the result that

$m(x)=\tanh\beta\left(\int dyQ(x,y)m(y)+h(x)\right),$ (18)

which is a generalization of a single mean field to an infinite set of order
parameters. It is not hard to see that it can be obtained from a continuous
field theory. Now define the following $M(y)=\text{atanh}(m(y))$. Thus we have
the following nonlinear integral equation of the second kind:

$M(x)=\beta h(x)+\beta\int dyQ(x,y)\tanh(M(y)),$ (19)

whose solution can be obtained via iterations.

Additionally, we can see that the effective Hamiltonian of eqn. (15) is
composed of a mean field energy Hamiltonian given by

$H=\frac{1}{2}\sum_{A=1,B=1}^{K}Q_{AB}m_{A}m_{B}+\sum_{A=1}^{K}h_{A}m_{A}.$
(20)

We will use this Hamiltonian as an effective energy for our islands’
interaction.

##  Appendix D The case of two islands: bipartite static mean field equations

The case of two nano-islands features both the Curie-Weiß temperature of the
single islands and the interaction between the islands themselves determining
the state of the magnetic moments. Specifically, this is a scenario in which
two nanoislands can have both ferro- and antiferro-magnetic interactions, yet
behave individually as mean field ferromagnets, as shown in Fig. 4. In the
standard treatment of the system with discrete spins: the two islands tend to
align head-to-head or head-to-tail at low temperatures. This is essentially
right, but as we show below, at the level of mean field, the phase diagram of
the two cases is different.

Figure 5: The coupling between the nanoislands: for orthogonal islands it is
anti-ferromagnetic, while for parallel it is ferromagnetic.

As stated previously,

$\displaystyle m_{A}$ $\displaystyle=$
$\displaystyle\tanh\left(\theta_{A}m_{A}+\theta_{AB}m_{B}\right),$
$\displaystyle m_{B}$ $\displaystyle=$
$\displaystyle\tanh\left(\theta_{B}m_{B}+\theta_{AB}m_{A}\right),$ (21)

A numerical solution to these equations determines three solutions over
different parameter regimes. In phase 1, $m_{A}=m_{B}=0$, the demagnetized
state. In phase 2, $m_{A}=m_{B}=1$, the fully ferromagnetic case. We can see
however that there is a third phase, Phase 3, which can be seen in Fig. 6
(green) at which the system is in an intermediate state $m_{A}=m_{B}=1/2$.
While in the top left diagram this is easy to see, this phase shrinks and lies
at the boundary between phase $A$ and $B$ for even largers $\theta_{ab}$, and
it disappears only when the diagram is in phase 1.

The diagram of what parameters correspond to these phases for
$\theta_{ab}=0.05,0.5,1$ is presented in Fig. 6.

Figure 6: Phase diagram for the bipartite mean field model as a function of
$\theta_{a}$ and $\theta_{b}$, for $\theta_{ab}\in[0,1]$. Yellow is a phase
$3$, green is a phase $2$, while blue is a phase $1$. It is interesting to
note that between the phase $1$ and $3$ there is almost everywhere a
interfacial phase $2$, aside along the line $\theta_{a}=\theta_{b}$, in which
there is a direct transition $1$ to $3$. We see that for
$\theta_{ab}\rightarrow 1$, the phase $1$ shrinks, leaving space for a phase
$3$.

For $\theta_{AB}<0$ the interfacial region with a phase $2$ disappears, and
there is a transition between the phase $1$ and phase $2$ only. See for
instance Fig. 6.

As a result, we see that although the macroscopic notion of head-to-tail and
head-to-head behavior of the two islands is essentially correct, at the mean
field level one can have interfacial mixed boundaries between the fully
magnetized and non-magnetized case as the temperature is lowered, and
depending on the interaction strength (distance between the island).

##  Appendix E Mean field spin ice: Dynamics

We consider three models of dynamics for the spin ice effective mean field
model: Model A, Metropolis, and Glauber dynamics Hohenberg and Halperin
(1977). Both in the model A and the Metropolis moves, we enforce $-1\leq
m_{a}\leq 1$. We will compare our results again the Glauber Glauber (1963)
dynamics for the discrete spins, with an acceptance probability for a spin
flip given

$\displaystyle p_{A\rightarrow B}=\frac{e^{-\frac{\Delta E}{\kappa
T}}}{1+e^{-\frac{\Delta E}{\kappa T}}}$ (22)

with $\Delta E=2\sum_{b}J_{ab}s_{b}$.

###  E.1 Model A dynamics

The first model we consider is the so called Model A relaxational dynamics,
defined by Hohenberg and Halperin (1977)

$\displaystyle\frac{\partial m_{A}}{\partial_{t}}=-\gamma\frac{\delta
H}{\delta m_{A}}+\xi_{A}(t),$ (23)

with $\xi_{A}(t)$ being a stochastic term with the properties

$\displaystyle\langle\xi_{A}(t)\rangle$ $\displaystyle=$ $\displaystyle 0$
(24) $\displaystyle\langle\xi_{A}(t)\xi_{B}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle 2T\gamma\delta_{AB}\delta(t-t^{\prime}).$ (25)

above, $\gamma$ is the relaxation rate.

The solution to such a linear stochastic differential equation is well known:

$\displaystyle\vec{s}(t)=e^{-\gamma\bm{Q}t}\vec{s}_{0}+e^{\gamma\bm{Q}t}\int_{0}^{t}ds\
e^{\gamma\bm{Q}s}\vec{\xi}(s),$ (26)

and thus

$\displaystyle\langle\vec{s}(t)\rangle$ $\displaystyle=$ $\displaystyle
e^{-\gamma\bm{Q}t}\vec{s}_{0}+e^{\gamma\bm{Q}t}\int_{0}^{t}ds\
e^{\gamma\bm{Q}s}\langle\vec{\xi}(s)\rangle,$ (27) $\displaystyle=$
$\displaystyle e^{-\gamma\bm{Q}t}\vec{s}_{0}.$ (28)

If we consider that $-1\leq s_{i}(t)\leq 1$, the asymptotic state is the one
determined by the maximum positive eigenvalue. If the maximum eigenvalue is
negative, the asymptotic state is unmagnetized.

In the case of a four island vertex, these dynamics agree with simple
energetic considerations. Considering four islands where $s_{i}=1$ when they
point up and to the right,

$\mathbf{Q}=\begin{pmatrix}q&-J_{1}&J_{2}&J_{1}\\\ -J_{1}&q&J_{1}&J_{2}\\\
J_{2}&J_{1}&q&-J_{1}\\\ J_{1}&J_{2}&-J_{1}&q\\\ \end{pmatrix}$ (29)

The eigenvalues and corresponding eigenvectors of this matrix are

$\mathbf{v}_{1}=\begin{pmatrix}1\\\ -1\\\ -1\\\
1\end{pmatrix},\mathbf{v}_{2}=\begin{pmatrix}1\\\ 0\\\ 1\\\
0\end{pmatrix},\mathbf{v}_{3}=\begin{pmatrix}0\\\ 1\\\ 0\\\
1\end{pmatrix},\mathbf{v}_{4}=\begin{pmatrix}-1\\\ -1\\\ 1\\\ 1\end{pmatrix}.$
(30) $\lambda_{1}=q+2J_{1}-J_{2},$ $\lambda_{2}=\lambda_{3}=q+J_{2},$
$\lambda_{4}=q-2J_{1}-J_{2}.$ (31)

If we only consider positive nonzero parameters, the system asymptotically
approaches $\pm\mathbf{v}_{1}$ when $J_{1}>J_{2}$ and a superposition of
$\pm\mathbf{v}_{2}$ and $\pm\mathbf{v}_{3}$ when $J_{2}>J_{1}$. The former
case corresponds to both possible Type I vertices while the fully saturated
superpositions of the latter are the four Type II vertices, precisely the
local low energy states of square ASI. When $J_{1}=J_{2}$, the eigenvectors
abruptly change:

$\mathbf{v}_{1}=\begin{pmatrix}1\\\ -1\\\ 0\\\
0\end{pmatrix},\mathbf{v}_{2}=\begin{pmatrix}1\\\ 0\\\ 1\\\
0\end{pmatrix},\mathbf{v}_{3}=\begin{pmatrix}0\\\ 1\\\ 0\\\
1\end{pmatrix},\mathbf{v}_{4}=\begin{pmatrix}-1\\\ -1\\\ 1\\\ 1\end{pmatrix}.$
(32) $\lambda_{1}=\lambda_{2}=\lambda_{3}=q+J_{1},$
$\lambda_{4}=q-2J_{1}-J_{2}.$ (33)

For all positive and nonzero parameters, all six ice rule obeying vertices are
possible as average asymptotes.

Figure 7: Asymptotic vertex types for model A dynamics. In both plots, $q=10$
and $J_{2}=1$. Left: $J_{1}$ is varied while $\xi=.1$. Right: $\xi$ is varied
while $J_{1}=1$.

To confirm that this is relevant without noise averaging, the stochastic
differential equations were Euler integrated and the asymptotic states were
recorded in terms of vertex types (Fig. 7) . $J_{1}$ was set to a range of
parameters in one case, while the square root of the variance of the noise
$\xi$ was modified in another case. In all scenarios, $q=10$ and $J_{2}=1$.
When $J_{1}$ varied, $\xi=.1$, and when $\xi$ varied, $J_{1}=1$. In the former
case, Type 1 vertices dominate for high $J_{1}$, Type 2 dominate and coexist
with Type 3 vertices when $J_{1}$ is low. In the latter case, Type I and III
occupy 25% of asymptotic states and Type II are 50 % of the states. This
gradually transitions to the ergodic arrangements of states at high $\xi$,
with an intermediate region where the average spin squared dips to around 50%.

###  E.2 Metropolis moves

The Metropolis dynamics, employed in the paper to reproduce the populations of
the single vertex, is given a two step move. The first, is a proposed move
given by a diffusion of the mean field value, as in eqn.(25). We pick randomly
a variable $m_{A}$, and propose a diffusive move

$\displaystyle m_{A}\rightarrow m_{A}^{\prime}=m_{A}+\xi$ (34)

with $\langle\xi\rangle=0$, and $\langle\xi(t)^{2}\rangle=2\kappa T\gamma$.
Then, the accept the move with probability given by the Metropolis acceptance
rate,

$\displaystyle r=\text{min}(1,e^{-\frac{H-H^{\prime}}{\kappa T}}),$ (35)

which enforces detailed balance. This is the model we simulated in the main
text, with $T=J/2$ and $\sigma=0.05$.

###  E.3 Mean field Glauber dynamics

Given the Hamiltonian

$H=\sum_{AB}\mathcal{A}_{AB}J_{AB}\sum_{a\in A}s_{a}\sum_{b\in
B}s_{b}+\sum_{A}h_{A}(\sum_{i\in A}s_{i})$ (36)

In the following we will assume $J_{AB}\equiv J$. We consider the following
Glauber dynamics Glauber (1963) for each individual spin. We have

$\displaystyle\frac{w_{a}(\vec{s})}{w_{a}(s_{a})}$ $\displaystyle=$
$\displaystyle\frac{P_{eq}(s_{a})}{P_{eq}(s)}=e^{-2\beta
Js_{a}\sum_{b}A_{ab}s_{b}-2\beta h_{a}s_{a}}$ (37) $\displaystyle=$
$\displaystyle\frac{1-s_{a}\tanh(\beta J\sum_{b\in N(a)}A_{ab}s_{b}+\beta
h_{a}s_{a})}{1+s_{a}\tanh(\beta J\sum_{b\in N(a)}A_{ab}s_{b}+\beta
h_{a}s_{a})}.$

As it is known, a rate which satisfies the detail balance is given by

$w_{i}(s)=\frac{1}{2}\left(1-s_{i}\tanh(\beta\sum_{j\in N(i)}A_{ij}s_{j}+\beta
h_{i}s_{i})\right)$ (38)

The dynamics for each spin is given by

$s_{i}(t+\Delta t)=\begin{cases}s_{i}(t)&p_{+}=1-w(s)\Delta t\\\
-s_{i}(t)&p_{-}=w(s)\Delta t\\\ \end{cases}$ (39)

Thus,

$\frac{ds_{i}}{dt}=-2\langle s_{i}w(s)\rangle=-\langle
s_{i}\rangle+\langle\tanh(\beta\sum_{j\in\mathcal{N}(i)}J_{ij}s_{j}+\beta
h_{i}s_{i})\rangle$ (40)

If we now replace the Hamiltonian for the spin ice, we note that can write
coarse grained differential equations by replacing $S_{j}^{A}=\langle
s_{j}\rangle$, and defining $m_{A}(t)=\langle\frac{1}{N}\sum_{j\in
A}s_{j}(t)\rangle$. We obtain as a result the following differential
equations:

$\dot{\mathbf{m}}=-\mathbf{m}+\tanh(\beta\mathbf{Q}\mathbf{m}+\beta\mathbf{h}),$
(41)

which generalizes the mean field equations we obtained for the static case to
the dynamical regime. For a larger set of mean field parameter we can use also
these equations to search for solutions of the static mean field equations.

It is known that for the case of synchronous Glauber dynamics the asymptotic
state distribution is described by

$P_{eq}(s)=\mathcal{N}e^{-\beta\tilde{H}(s)},$ (42)

where $\tilde{H}=\frac{1}{\beta}\sum_{A}\log 2\cosh(\beta\Delta_{a})$ is
Peretto’s pseudo-Hamiltonian Peretto (1984) and $\mathcal{N}$ a normalization
constant, with $\Delta_{a}=J\sum_{B}A_{ab}m_{b}$.333The change of sign
compared to Peretto is because we are considering an anti-ferromagnetic
system. For large values of $\beta$ (i.e. low temperature), we have

$\tilde{H}\underbrace{\approx}_{\beta\gg 1}J\sum_{a}|\sum_{b}A_{ab}m_{b}|$
(43)

which is different from the original Hamiltonian, in particular it is linear
in the spins, but depends on the absolute value for an arbitrary coupling
matrix $J_{ij}$. Solving for the ground state of the system above can be done
via a linear program.

For the case of spin ices this is sufficient to recover the right ground
state, as the ground state of a spin ice system is vertex-frustrated. We saw
at the beginning of the paper that we can write

$H=J\sum_{v}(\sum_{\beta=1}B_{v,\beta}s_{\beta})^{2}.$ (44)

It is not surprising then to see that the ground state (not the excitations)
of the system above is equivalent to the one of

$\displaystyle\tilde{H}$ $\displaystyle=$ $\displaystyle
J\sum_{\alpha}|\sum_{v}B_{\alpha,v}\sum_{\beta=1}B_{v,\beta}s_{\beta}|$ (45)
$\displaystyle\leq$ $\displaystyle\sum_{v}\sum_{\alpha}J|B_{\alpha
v}|\sum_{\beta=1}B_{v,\beta}s_{\beta}|$ $\displaystyle\leq$
$\displaystyle\sum_{v}\tilde{J}_{v}|\sum_{\beta=1}B_{v,\beta}s_{\beta}|$
$\displaystyle\leq$
$\displaystyle\tilde{J}_{max}\sum_{v}|\sum_{\beta=1}B_{v,\beta}s_{\beta}|$

which is the pseudo-Hamiltonian for the spin ice model from the mean field
dynamics. In this sense, as long as we only care about the ground state of the
model, the mean field dynamics should provide the right asymptotic ground
states; this statement depends on the fact that one can have ground states in
which $\sum_{\beta}B_{v\beta}s_{\beta}=0,\forall v$, which are exactly the ice
states. This fact is shown in the next section numerically.

##  Appendix F Mean field dynamics: single spin

The mean field dynamics of an Ising system driven by Glauber dynamics reduce
to

$\dot{m}(t)=-m(t)+\tanh[\beta Q(m(t)+\beta h(t))]$ (46)

for a single spin. We can analyze this simplified system in a variety of
parameter regimes to qualitatively understand the limits of the systems
responsiveness to external field and the ”quality” of the system. The quality
is an appropriate measure of the similarity of outputs given similar inputs
(kernel quality) minus the similarity of outputs with dissimilar inputs
(generalization capability). A natural measure of this emerges from a high
temperature analysis while the slow dynamics have simple enough behavior that
one may comment on quality by inspection.

###  F.1 Slow Dynamics

When considering this system as a reservoir, it is pertinent to analyze the
response of $m(t)$ to $h(t)$. We begin in the case of slow dynamics where
$m(t)$ is close to a fixed point and $h(t)$ evolves sufficiently gradually. In
this case, the derivative of the mean moment approaches zero, yielding a self
consistency equation for the magnetic moment:

$m=\tanh[\beta Qm+\beta h].$ (47)

This takes the form of a static mean field Ising system. First, let us
consider the ferromagnetic case ($Q>0$). At zero field there is a paramagnetic
to ferromagnetic transition when the slope of the hyperbolic tangent evaluated
at zero grows equal to one, resulting in a pitchfork bifurcation when $\beta
Q=1$. A single stable solution, $m=0$, transitions to an unstable $m=0$ and
two stable nonzero solutions resulting from the new intersection of the
hyperbolic tangent with the line. The nonzero solutions for $m$ are stable and
grow in value with decreasing $\beta$, corresponding to the system settling
into a two-fold degenerate ground state. When $h$ is small but not zero, the
values of the solutions are shifted in the direction of the field. With a
large enough $h$, the unstable and smaller stable points converge to a single
unstable fixed point in a saddle node bifurcation. The physical interpretation
of this convergence to a single moment in the direction of the field is that
the coercivity of the system has been overcome and a hysteritic moment
opposing the external field must flip.

Figure 8: Stable fixed points plotted at different temperatures with $Q=1$.
The paramagnetic to ferromagnetic transition is apparent at $k_{B}T=1$ as the
solutions split into two curves. Magnetic susceptibility and coercivity
increases as the temperature is further cooled.

One can infer that limiting this system to slow dynamics constrains it to act
like a simple hystersis function. Any applied field can be mapped to an output
field if the system’s previous state is known. This directly defies the echo
state property and constrains the parameter space of the Ising systems used in
reservoir computing. If the temperature is too low with slow dynamics, a
nonzero coercivity makes for ambiguous computation. This may otherwise be
described in the language of the quality of the system. When the coercivity is
non-zero, the kernel quality decays as similar inputs may produce different
outputs due to the influence of the system’s history. Further, whenever the
system is close to saturation, the generalization capability falls as
different fields in this range all produce a similarly high mean moment.
Generally speaking, a slow reservoir computes well when coercivity is low
(intermediate temperature) and the system is not near saturation (low field).

The antiferromagnetic case ($Q<0$) behaves almost precisely like the
paramagnetic regime of the ferromagnetic case. When the temperature is further
reduced, the shape of the hyperbolic tangent function approaches a step
function, leading to a sharper and sharper response of the magnetic moment.
This lack of hystersis makes it a more reliable component of a reservoir.

###  F.2 High Temperature Limit

When the temperature is sufficiently high, i.e. where $Q<<1$ and $Qh(t)<<1$,
we may linearize the equations of motion and attain exact solutions.
Approximating to first order,

$\dot{m}(t)=-\lambda m(t)+h(t).$ (48)

where the characteristic timescale $\lambda=1-Q$. This equation equivalently
describes the dynamics of a simple RC circuit. An ansatz to this equation
ansatz may be written as a homogenous and inhomogenous solution:
$m(t)=m_{h}(t)+m_{i}(t)$. The homogenous solution, $m_{h}(t)=A\exp(-\lambda
t)$, decays only when $\lambda>0$. An exponential growth of this term would be
unphysical and suppressed by higher order terms in the expansion. The analysis
within this section is therefore only valid when $Q<1$, consistent with the
previous assumptions.

The particular solution depends on the timeseries $h(t)$. Assuming this
function is periodic with frequency $\omega_{h}$, we write the Fourier series
decomposition:

$h(t)=\sum_{n}C_{n}\exp(in\omega_{h}t).$ (49)

Matching coefficients for the inhomogenous solution, we find

$m_{i}(t)=\sum_{n}\frac{1}{in\omega+\lambda}C_{n}\exp(in\omega t).$ (50)

The coefficients are reduced by a factor proportional to the number of the
harmonic, indicating the low pass filtering property which becomes most
prominent when $n\omega$ and $\lambda$ are of the same order, implying that
$\lambda$ is the natural cutoff frequency of the system.

To understand the quality of this system, we calculate how the asympotic and
therefore imhomogenous solutions with different driving fields correlate over
time. Defining a correlator $\langle m(t)m^{{}^{\prime}*}(t)\rangle$ where
$m^{{}^{\prime}}(t)$ solves $\dot{m}^{{}^{\prime}}(t)=-\lambda
m^{{}^{\prime}}(t)+h^{{}^{\prime}}(t)$ and
$h^{{}^{\prime}}(t)=\sum_{n}C_{n}^{{}^{\prime}}\exp(in\omega_{h}t)$.
$\omega_{h}$ is the same as the other signal because most reservoir computers
fix input lengths and any harmonics that aren’t integer multiples of each
other will have vanishing overlaps. The average is conducted forward in time:

$\langle
m(t)m^{{}^{\prime}*}(t)\rangle=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}m(t)m^{{}^{\prime}*}(t)dt$
(51)

Choosing a starting time arbitrarily late enough such that the homogenous
solution approaches zero, only the inhomogenous solution is nonzero.

$\langle
m(t)m^{{}^{\prime}*}(t)\rangle=\lim_{T\to\infty}\frac{1}{T}\sum_{nn^{{}^{\prime}}}\int_{0}^{T}\frac{C_{n}C_{n^{\prime}}^{*}}{n^{2}\omega^{2}+\lambda^{2}}\exp(i(n-n^{{}^{\prime}})\omega
t)dt$ (52)

Over a single period, the integrand will always go to zero unless the argument
of the exponential goes to zero. This leaves only the terms where
$n=n^{{}^{\prime}}$ and the overlap simplifies to:

$\langle
m(t)m^{{}^{\prime}*}(t)\rangle=\sum_{n}\frac{C_{n}C_{n}^{*}}{n^{2}\omega^{2}+\lambda^{2}}$
(53)

It follows then that at high frequencies this model gives a quality dependency
of the type $1/\omega^{2}$.

##  Appendix G Vertices and Vortices

###  G.1 Vertices

The famed degeneracy of artificial spin ice is often analzyed through the lens
of energetics. That is, extensively numerous states of the system possess the
same energy. The states are categorized at vertices, natural building blocks
of the system. We analyze an isolated vertex from the square ice system to
observe how ice properties may be achieved dynamically and what the
transitions between states imply. To model this system, we label all four
$m_{i}=1$ states as pointing up and to the right, indexed counterclockwise
from the rightmost magnet. The external field is applied equally to all
magnets along the same direction. The interactions between magnets are
parameterized in terms of $J_{1}$, the interaction between spins oriented
perpendicular to one another which is typically stronger in planar square ice,
and $J_{2}$, the typically weaker interaction of collinear macrospins. The
interaction matrix is thus

$\mathbf{Q}=\begin{pmatrix}q&-J_{1}&J_{2}&J_{1}\\\ -J_{1}&q&J_{1}&J_{2}\\\
J_{2}&J_{1}&q&-J_{1}\\\ J_{1}&J_{2}&-J_{1}&q\\\ \end{pmatrix}$ (54)

and the field vector is

$\mathbf{h}(t)=h(t)\begin{pmatrix}1\\\ 1\\\ 1\\\ 1\end{pmatrix}.$ (55)

###  G.2 Low magnetization fixed point

Considering the case a small argument to the $\tanh$ function and a static
field $h(t)=h$, one finds a fixed point at $m_{i}=\frac{h}{q+J_{2}-1/\beta}$
for all $i$. This intuitively implies that the fixed point is slightly offset
from 0 by the external field, shifting the entirety of the phase space and
therefore allowing the system to switch between basins of attraction. There is
a critical temperature where $1/\beta=q+J_{2}$ that maximizes the sensitivity
to external field, implying that this is an ideal temperature for reservoir
computing as only small amounts of driving are required to alter the entire
system’s state. This will be explored further in the numerical analysis of the
phase space.

The eigenvectors of the Jacobian about the saddle node are all associated with
real eigenvalues that change sign depending on the system parameters. Positive
eigenvalue directions are evolved towards and typically manifest as the system
approaching a high $m_{i}$ state of a energetically favorable vertex type.
Specifically, the eigenvectors and their corresponding eigenvalues are the
following:

$\mathbf{v}_{1}=\begin{pmatrix}1\\\ -1\\\ -1\\\
1\end{pmatrix},\mathbf{v}_{2}=\begin{pmatrix}1\\\ 0\\\ 1\\\
0\end{pmatrix},\mathbf{v}_{3}=\begin{pmatrix}0\\\ 1\\\ 0\\\
1\end{pmatrix},\mathbf{v}_{4}=\begin{pmatrix}-1\\\ -1\\\ 1\\\ 1\end{pmatrix}.$
(56)

The corresponding eigenvalues are

$\lambda_{1}=\beta(q+2J_{1}-J_{2})-1,$
$\lambda_{2}=\lambda_{3}=\beta(q+J_{2})-1,$
$\lambda_{4}=\beta(q-2J_{1}-J_{2})-1.$ (57)

The first eigenvector is precisely the direction of a Type I vertex, the
second and third comprise collinear spin ordering that together lead to Type
II states, and the fourth direction is a type IV vertex. In all cases, beta
may be reduced low enough to make the eigenvalues negative, effectively
raising the temperature high enough to ”melt” the individual spins into
paramagnets and make the fixed point stable. On the other extreme, the self
energy $q$ dominates the interactions and the fixed point is unstable. When
parameters are similar to experiments, that is, $\beta q>>1$, and the rest of
the parameters are on the same order of magnitude, the relative strength of
$J_{1}$ and $J_{2}$ makes the Type I and Type II associated eigenvalues
respectively more positive, leading towards a faster and more dominant
evolution of the system towards those vertex types.

###  G.3 Basins of attraction

To characterize the phase space, we initialized the system across evenly
spaced points on the physically relevant interval [-1,1]. We then evolved the
system over time until the system converged and recorded the closest Ising
values of $m_{i}$ and average net moment, $\langle m_{i}^{2}\rangle$. We then
categorize the Ising states into vertex types and record the average fraction
of vertex types as parameters are varied 11).

The temperatures between $10^{-1}$ and $10^{1}$ represent the expected
behavior for a spin ice. The Type I and II vertices occur proportional to the
1:2 ratio of microstates they comprise while all other states are near zero.
Below $10^{-1}$, a new type III state emerges due to the ability to partially
magnetize individual spins. Finding a low energy configuration through
reducing the moments’ magnitudes is not found in nanoscale experiments,
potentially due to the inability to appropriately cool the systems, but does
occur in macroscopic dipole spin ice.

Figure 9: Fraction of the system evolving towards vertex types as a function
of temperature for $J_{1}=100$, $J_{2}=100$, $q=10$, and $h=0$. Both the
fraction and average magnetization parameter, $\langle m_{i}^{2}\rangle$, are
plotted on the same scale. The critical temperature of maximum fluctuation is
plotted at the goldenrod line.  Figure 10: Vertex types at a range of field
strengths for $J_{1}=100$, $J_{2}=100$, $q=10$, and three separate
temperatures.  Figure 11: Fraction of the system evolving towards vertex types
as a function of temperature for $J_{1}=100$, $J_{2}=100$, $q=10$, and $h=0$.
Both the fraction and average magnetization parameter, $\langle
m_{i}^{2}\rangle$, are plotted on the same scale. The critical temperature of
maximum fluctuation is plotted at the goldenrod line.

###  G.4 Vortices

Figure 12: Interactions for a single island: we have diagonal, horizontal and
vertical interactions.

Recent experiments have shown the potential for real spin islands to develop
vortex states: spin textures no longer dominated by a single direction but
rather spins that wrap around the center of the nanomagnet. A simple model to
capture the emergence of a vortex in this mean field format is one in which
the island is divided into four quadrants, each with an individual mean field
spin. By symmetry, there are three interaction strengths between these
quadrants: $J_{f}$ between collinear spins with ferromagnetic dipolar
interactions, $J_{a}$ between parallel, non-collinear spins with
antiferromagnetic dipolar interactions, and $J_{d}$ between spins across the
diagonal from one another. In the convention that all spins are positive when
pointing to the right, the interaction matrix is

$\mathbf{Q}=\begin{pmatrix}q&J_{f}&J_{d}&J_{a}\\\ J_{f}&q&J_{a}&J_{d}\\\
J_{d}&J_{a}&q&J_{f}\\\ J_{a}&J_{d}&J_{f}&q\\\ \end{pmatrix}$ (58)

Applying the same fixed point analysis as before, we find that all
$m_{i}=\frac{h}{q+J_{d}+J_{f}-J_{a}}$ and the following eigenvectors of the
Jacobian:

$\displaystyle\mathbf{v}_{mag}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}1\\\ 1\\\ 1\\\
1\end{pmatrix},\mathbf{v}_{1vort}=\frac{1}{2}\begin{pmatrix}-1\\\ -1\\\ 1\\\
1\end{pmatrix},$ (59) $\displaystyle\mathbf{v}_{nomag}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}-1\\\ 1\\\ -1\\\
1\end{pmatrix},\mathbf{v}_{2vort}=\frac{1}{2}\begin{pmatrix}1\\\ -1\\\ -1\\\
1\end{pmatrix},$ (60)

where the factor $\frac{1}{2}$ is for normalization. The corresponding
eigenvalues are

$\displaystyle\lambda_{mag}$ $\displaystyle=$
$\displaystyle\beta(q+J_{a}+J_{f}+J_{d})-1,$ $\displaystyle\lambda_{1vort}$
$\displaystyle=$ $\displaystyle\beta(q-J_{a}+J_{f}-J_{d})-1,$
$\displaystyle\lambda_{nomag}$ $\displaystyle=$
$\displaystyle\beta(q-J_{a}-J_{f}+J_{d})-1,$ $\displaystyle\lambda_{2vort}$
$\displaystyle=$ $\displaystyle\beta(q+J_{a}-J_{f}-J_{d})-1.$ (61)

####  G.4.1 Statics

We note that the Hamiltonian can be written in the form

$\displaystyle H_{island}=\frac{1}{2}\vec{s}^{t}\ \mathbf{Q}\vec{s}$ (62)

with $\vec{s}=(s_{1},s_{2},s_{3},s_{4})$, and that thus the eigenvectors are
the same. Let us call $\rho_{k}$ the corresponding eigenvalues. We have
$\lambda_{k}=\beta\rho_{k}-1$. The semi-magnetized states can be written then
as superpositions of the magnetized, unstructured unmagnetized, and the 1 and
2 vortices solutions.

This also implies that, if we call $c_{k}=\vec{s}\cdot\vec{v}_{k}$, we can
rewrite

$\displaystyle H_{island}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\rho_{mag}c_{mag}^{2}+\rho_{nomag}c_{nomag}^{2}$
(63) $\displaystyle+$
$\displaystyle\rho_{1vort}c_{1vort}^{2}+\rho_{2vort}c_{2vort}^{2})$ (64)

Minimizing the Hamiltonian is thus associated to the minimum eigenvalue as a
function of $J_{a},J_{f},J_{d}$. Since $q$ is a shift, we can set it to zero
in this case. Moreover, we have the freedom of setting $J_{f}=-1$, thus
reducing the phase diagram to $J_{d}$ and $J_{a}$ in units of $J_{f}$. This is
equivalent to the minimum eigenvalue problem for the matrix

$\mathbf{\tilde{Q}}=\begin{pmatrix}0&-1&-\tilde{J}_{d}&-\tilde{J}_{a}\\\
-1&0&-\tilde{J}_{a}&-\tilde{J}_{d}\\\ -\tilde{J}_{d}&-\tilde{J}_{a}&0&-1\\\
-\tilde{J}_{a}&-\tilde{J}_{d}&-1&0\\\ \end{pmatrix}$ (65)

with eigenvalues

$\displaystyle\tilde{\rho}_{mag}$ $\displaystyle=$
$\displaystyle-\tilde{J}_{a}-\tilde{J}_{d}-1$ (66)
$\displaystyle\tilde{\rho}_{nomag}$ $\displaystyle=$
$\displaystyle\tilde{J}_{a}-\tilde{J}_{d}+1$ (67)
$\displaystyle\tilde{\rho}_{1vort}$ $\displaystyle=$
$\displaystyle\tilde{J}_{a}+\tilde{J}_{d}-1$ (68)
$\displaystyle\tilde{\rho}_{2vort}$ $\displaystyle=$
$\displaystyle-\tilde{J}_{a}+\tilde{J}_{d}+1$ (69)

####  G.4.2 Dynamics

For the dynamics, however, $q$ becomes important as it affects the role that
that temperature has. We can largely draw the same conclusions with respect to
the $q$ dependence: larger values of $q$ drive all eigenvalues towards
positivity which in turn makes all modes unstable, leading to an evolution
towards ergodic states. The 4th state is the first to become stable and
therefore cease to exist in an ensemble of systems. This state is effectively
the island splitting into two macrospins aligned energetically unfavorably
with one another. The remaining three states (1 being a macrospin, 2 being a
vortex, and 3 being two vortices) emerge when two interaction strengths
combined are larger than the third. This allows for coexistence of these three
states across a wide variety of interaction ranges while the precise size of
their basins of attraction are tunable via the balance of interaction
strengths (Fig. 5).

One may additionally observe a point of maximum sensitivity of the fixed point
to external field: $q=J_{a}-J_{d}-J_{f}$. This is especially relevant to the
exploitation of vortex formation in FMR facilitated reservoir computing. The
prediction of double vortices with appreciable populations should allow
another form of system output which in term increases the amount of processing
the system is capable of.

###  G.5 Square and Hexagonal Ices: Glauber vs. Mean field dynamics

Let us now check that this model reproduces also the ground states of the
Hexagonal Artificial Spin Ice (HeX), and the Artificial Square Ice (ASI).

The first thing to notice is that, given the eigenvalues of the effective
Ising interaction matrix $J_{ij}$, the ordering occurs for temperatures below
the maximum eigenvalues $T_{o}=\Lambda_{max}$. Above the temperature $T_{o}$,
the system is in a paramagnetic phase and the spins converge exponentially to
$m=0$. At temperatures below $T_{o}$, we observe an interesting phenomenon,
which we describe in Fig. 13 for HeX and ASI. We observe that as time evolves,
below the ordering temperatures the Ising models some of the spins will
converge towards values $+1$ and $-1$. Some of the spins instead converge
towards the value zero. The asymptotic states depends on the initial
condition. These spins are those that can be freely flipped without leaving
the ice manifold, and we can thus assign to them random positive and negative
spin values. Given this arrangement, in Fig. 13 we show the final
configuration both for ASI and HeX, which can be promptly checked to be in the
ice manifold.

Figure 13: Relaxation of each spin of the Artificial Square Ice (ASI) and the
Hexagonal Spin Ice (HeX) in the mean field approximation below the ordering
temperatures $T_{o}$, and the corresponding final states. It can be seen that
the system converges to the known correct ground state for the two systems.

At the temperature $T_{o}$, the mean field equations still converge to $m=0$,
but the regime is not exponential anymore, but a power law which can be
extrapolate to be roughly $\sqrt{t^{-1}}$, which is the mean field dynamical
exponent. A well known example of possible application of artificial square
ice is the possibility of having monopole like charges in spin ices without a
string tension. Such degeneracy can be obtained via mediated spin interactions
Caravelli (2020) or staggered island heights Farhan  _et al._ (2019).

Figure 14: Mean field magnetization of each spin of an ASI at $T=T_{o}$. We
observe that the average spin $m_{av}(t)=\frac{1}{N}\sum_{i}|m_{i}(t)|$
relaxes as a power law $t^{\alpha}$ , with exponent roughly $\alpha\approx
0.5$ consistent with the mean field approximation.

###  G.6 Embedding Logic gates

Let us also consider another application of this model, e.g. logic gates
embedded in magnetic nanoislands. As discussed in Caravelli and Nisoli (2020),
the asymptotic states are not necessarily in the spin ice manifold, but we
find interesting to see whether the mean field dynamics can be implemented to
obtain better results, compared to the discrete Glauber approach, to the
ground designed ground state of these models.

The gate functionality can be obtained Caravelli and Nisoli (2020) by studying
a simple spin Hamiltonian of the form

$H_{3}=\sigma_{3}\left(J(\sigma_{1}+\sigma_{2})+h\right),$ (70)

which we assume that the spins $\sigma_{1}$ and $\sigma_{2}$ are inputs, and
$\sigma_{3}$ is our output and free to fluctuate at a certain temperature $T$.
If $2|J|>h$ the ground state of the system are those of a logic gates (N)AND
and (N)OR, which is why this embedding is successful in reproducing a wide
range of logic gates.

The hyerarchy problem can be solved by modulating couplings $J$ among islands.
For the model without horizontal (input) spin interactions, we choose
$\epsilon$ such that $|{J_{k}}/({2+\epsilon})|<|h_{k}|$. As in Caravelli and
Nisoli (2020), we used the coupling rescaling between layers given by

$|J_{k-1}|\leq\frac{|J_{k}|}{(2+\epsilon)}.$ (71)

We performed numerical simulations of randomly chosen gates via a Glauber
Glauber (1963) spin dynamics with exponential annealing for a system whose
couplings scale as in Eq. (71), versus the mean field dynamics results. The
probability of a spin flip in the Glauber model is associated to the energy
change $\Delta E$, and is given by $p(k)=(1+e^{\Delta E/T_{k}(t)})^{-1}.$

As in previous papers Caravelli and Nisoli (2020), the control parameter for
the fidelity of the circuit we introduce the overlap between the spin system
and the equivalent logic circuit, not only for final output, but for all the
intermediate layers. For each gate, we consider the metric given by the gate
overlap as the quantity,

$\mathcal{Q}=\frac{\left[\vec{S}\cdot\vec{L}+(1-\vec{S})\cdot(1-\vec{L})\right]}{N},$
(72)

where $N$ is the total number of gates. Full functionality corresponds to
$\mathcal{Q}=1$, and completely random system have $Q=0$. In Fig. 15 we
compare the mean field and Glauber dynamics results, showing that the results
are consistent with previous studies.

Figure 15: Average corresponding energy and overlap for random gates in the
mean field approximation. We oberve that in the mean field approximation, the
average overlap is close to $1$ (full functionality of the gates) and higher
than the Glauber result for trees of depth $L=5,7,9$. Each line is averaged
over 100 simulations for the Glauber dynamics with discrete spins, and $100$
initial conditions for the mean field corresponding equations. The annealing
rate is identical for both systems.


# Observation of Traveling Breathers and Their Scattering in a Two-Fluid
System

Yifeng Mao<EMAIL_ADDRESS>Department of Applied Mathematics,
University of Colorado, Boulder, CO 80309, USA  Sathyanarayanan Chandramouli
Department of Mathematics, Florida State University, Tallahassee, FL 32306,
USA.  Wenqian Xu  Department of Applied Mathematics, University of Colorado,
Boulder, CO 80309, USA  Mark Hoefer  Department of Applied Mathematics,
University of Colorado, Boulder, CO 80309, USA

###### Abstract

The observation of traveling breathers (TBs) with large-amplitude oscillatory
tails realizes an almost 50-year-old theoretical prediction [1] and
generalizes the notion of a breather. Two strongly nonlinear TB families are
created in a core-annular flow by interacting a soliton and a nonlinear
periodic (cnoidal) carrier. Bright and dark TBs are observed to move faster or
slower, respectively, than the carrier while imparting a phase shift.
Agreement with model equations is achieved. Scattering of the TBs are observed
to be physically elastic. The observed TBs generalize to many continuum and
discrete systems.

Fifty years ago, special, localized two-soliton bound state solutions of the
sine-Gordon equation were called “breathers” because they incorporate two time
scales: one associated with propagation and the other associated with internal
oscillation [2]. The term breather has since been generalized to solutions of
other completely integrable systems such as the modified Korteweg-de Vries
(mKdV) and the focusing nonlinear Schrödinger (NLS) equations [3]. However,
the strict notion of localized breather solutions appears to be limited to
integrable equations. Numerical simulations [4, 5, 6] and mathematical
analysis [7, 8] demonstrate that breathers in some non-integrable equations
necessarily display small oscillatory tails. Consequently, they have been
referred to as quasi-breathers [3] or nonlocal solitary waves [9], although,
in some instances, these oscillatory tails may be vanishingly small [10].

Small amplitude breathers can be approximated by bright soliton solutions of
the focusing NLS equation. Consequently, breathers have been observed in many
situations such as water waves [11, 12], internal waves [13], nonlinear
optical [14, 15, 16] and matter [17, 18] waves, rogue waves in these systems
[19, 20, 21, 22, 23], magnetic materials [24], and discrete systems [25, 26].

Breathers with non-vanishing oscillatory tails have been studied extensively
in lattice systems where they are commonly referred to as traveling intrinsic
localized modes or traveling breathers [27, 25, 28] with experimental
observations of their existence in [29, 30]. In the continuum setting, special
solutions, originally called solitary dislocations, were derived for the KdV
equation as the nonlinear superposition of a soliton and a cnoidal traveling
periodic wave in the seminal work of Kuznetsov and Mikhailov [1]. This was
achieved using inverse scattering theory (IST) by introducing a discrete
eigenvalue, representing the soliton, in either the finite or semi-infinite
gap of the cnoidal wave spectrum for the Schrödinger operator associated with
KdV. This solution is a traveling breather with large amplitude oscillatory
tails that can be of either elevation (bright) or depression (dark) type
depending on the eigenvalue location in the semi-infinite gap or finite gap,
respectively. The cnoidal wave can also be considered a soliton lattice and so
this traveling breather can be viewed as a soliton-soliton lattice
interaction, generalizing the localized two-soliton bound state breathers in,
e.g., sine-Gordon and mKdV [3]. The scattering theory for KdV traveling
breathers indicates that they interact with each other elastically,
experiencing only a spatial shift [1]. Despite much ensuing analysis [31, 32,
33, 34, 35, 36, 37, 38, 39, 40] and corresponding physical interest, the
observation and scattering of strongly nonlinear traveling breathers in
continuum mechanics is lacking.

In this Letter, we create and scatter strongly nonlinear elevation wave or
bright (BB) and depression wave or dark (DB) traveling breathers at the
interface between two viscous fluids by the overtaking interaction of a
carrier wave with a soliton or vice-versa. The traveling breather oscillatory
tails—the carrier background—are large amplitude. The results experimentally
prove that this class of breathers can be interpreted as a nonlinear
superposition between a soliton and a cnoidal-like wave, realizing the
spectral interpretation of a traveling breather from IST developed almost 50
years ago [1, 40].

Additionally, we measure traveling breather properties during free propagation
and scattering, finding that they retain their solitonic character, i.e.,
interactions are physically elastic. Unimodal and bimodal geometries of
internal oscillation are observed. Qualitative or quantitative features of
individual traveling breathers and their scattering agree with predictions
from KdV traveling breather theory or numerical simulations of a viscous two-
fluid nonlinear model, respectively. Increasing traveling breather speeds
imply decreasing (BB) or increasing (DB) carrier phase shifts. Slower (faster)
traveling breathers exhibit a negative (positive) spatial shift post-
interaction. The observation of traveling breathers in this model continuum
system and their characterization are relevant to the physics of nonlinear
dispersive media more broadly.

The experiments are conducted in a tall acrylic column of $5\text{cm}\times
5\text{cm}\times 180\text{cm}$, described in [41]. The column consists of a
pressure-driven viscous core fluid with a free interface to a miscible,
heavier viscous reservoir fluid with a small core-to-reservoir fluid viscosity
ratio. We precisely control the injection of the buoyant, interior fluid in
the low Reynolds number regime, which is convectively unstable [42, 43], so
that straight conduits are established with constant injection. Time-varying
injection results in a conduit that exhibits an azimuthally symmetric
interface with cross-sectional area $A(z,t)$, where $z>0$ is the vertical
spatial coordinate in the camera view that is slightly above the injection
site. We identify $t=0$ as the time at the initiation of interfacial imaging.
At the injection site, the cross-sectional area satisfies the Hagen-Poiseuille
flow law $Q\propto A^{2}$ so it can be precisely controlled by temporally
varying the injection rate $Q$ [44, 45].

Viscous core-annular interfacial waves are modeled by the strongly nonlinear
conduit equation expressing mass conservation [45, 46]

$A_{t}+\left(A^{2}\right)_{z}-\left(A^{2}\left(A^{-1}A_{t}\right)_{z}\right)_{z}=0,$
(1)

given in non-dimensional form where $z=\tilde{z}/L$, $t=\tilde{t}/T$,
$A=\tilde{A}/A_{0}$ for dimensional quantities $\tilde{z},\tilde{t},\tilde{A}$
and $L$, $T$, $A_{0}$ are the fitted vertical length, time, and cross-
sectional area, respectively, obtained from separate measurements of linear
dispersive waves on the background injection rate $Q_{0}\propto A_{0}^{2}$
[47]. The conduit equation (1) and variants of it are models of magma flow
[48] and channelized water flow in glaciers [49]. The conduit equation is
apparently not integrable and can be reduced to the KdV equation in the long-
wavelength and small-amplitude regime [47, 50]. Cnoidal-like waves in the form
$A(z,t)=g(\theta)$, $\theta=kz-\omega t$, $g(\theta+2\pi)=g(\theta)$ and
solitons are generated by evaluating numerically computed traveling wave
solutions of the conduit equation [51] at the injection site, which determines
the injection rate time-series $Q(t)\propto g(-\omega t)^{2}$. The
characterization and reliable generation of soliton and cnoidal-like waves
have, individually, been reported in [45, 52, 41].

Figure 1: Bright traveling breather (BB) in a viscous fluid conduit formed
from a soliton overtaking a cnoidal-like periodic traveling wave. (a)
$90^{\circ}$ clockwise rotated time-lapse. (b) Space-time contour of the
cross-sectional area. (c) Background cnoidal-like wave and BB measurements at
different heights. (d) Spatial and temporal profiles (black solid) and the
cnoidal-like wave solution of the conduit equation (1) corresponding to
measured parameters (dashed).

Figure 1 displays the spatio-temporal evolution of a soliton overtaking a
cnoidal-like wave. At $t=0$, we have prepared the conduit as a nonlinear
cnoidal-like wave (the carrier) that abruptly terminates at the injection site
to a constant flow rate for $t>0$. The carrier propagates upward for positive
$z$ and experiences a modulation region as it transitions from periodic to a
constant background (Fig. 1a for $t>20$s, Fig. 1b) [47]. This modulated region
consists of a dispersive shock wave (DSW) [53], a constant region whose value
$A_{\rm min}$ coincides with the minimum of the adjacent carrier, and a
carrier modulation [53, 54]. The averaged carrier’s amplitude ($a_{\rm cn}$),
mean ($\overline{A}_{\rm cn}$), and dimensional angular frequency
($\tilde{\omega}=\omega/T$) are reported in Fig. 1c. The measured angular
wavenumber ($\tilde{k}=k/L$) is $\tilde{k}=2.35\pm 0.01\text{cm}^{-1}$. After
terminating the carrier, a large-amplitude soliton is injected, which is
transmitted through the DSW [55] and the constant region $A_{\rm min}$.
Because the soliton speed on the constant region $A_{\rm min}$ exceeds the
carrier’s phase speed ($\tilde{v}_{\rm ph}=\tilde{\omega}/\tilde{k}$), the
soliton overtakes the carrier and forms a new coherent structure. The
amplitude of the coherent structure extracted from experiment ($a_{\rm
BB}(z)=\max_{t}A(z,t)-A_{\rm min}$) oscillates with propagation. In Fig. 1c,
we plot $a_{\rm BB}(z)$, its envelope, and the carrier’s properties extracted
from each spatial slice of Fig. 1b to characterize propagation up the conduit.
The amplitude, mean, and frequency of the carrier and the coherent structure’s
amplitude envelope exhibit small fluctuations with height $z$. We conclude
that the coherent structure is a BB consisting of the nonlinear superposition
of a soliton and a cnoidal-like wave.

Next, we extract the BB’s maximum value for $\tilde{z}\gtrsim 40$cm in Fig. 1b
that follows a zig-zag path with constant speed regularly interspersed between
regions of rapid acceleration/deceleration. The BB speed across the entire
zig-zag path is $\tilde{v}_{\rm BB}=0.63\pm 0.01$cm/s [47], which exceeds the
measured carrier phase speed $\tilde{v}_{\rm ph}=0.319\pm 0.005$cm/s. This is
consistent with KdV BB solutions [1, 39]. A more thorough analysis of the BB
trajectory [47] shows that it can also be viewed as a soliton undergoing a
sequence of phase shifts while interacting with a soliton lattice that
composes the cnoidal-like wave [1, 56].

In Fig. 1d, we extract experimental time and spatial slices of $A(z,t)$. The
measured BB’s phase shift in space and time $\Delta\theta\in(-\pi,\pi]$ is
defined as the difference between the left and right carrier phases. We find
that $\Delta\theta_{\rm BBz}=0.72\pi$ in space and $\Delta\theta_{\rm
BBt}=0.75\pi$ in time are comparable, as expected. This BB with a positive
phase shift is consistent with large amplitude BB solutions of the KdV
equation [39].

Figure 2: Absorption of a solitary wave by a larger-mean, larger-amplitude
cnoidal-like periodic traveling wave resulting in a dark traveling breather.
(a) $90^{\circ}$ clockwise rotated time-lapse images. (b) Space-time contour.
(c) Background cnoidal-like wave and DB measurements. (d) Spatial and temporal
profiles (black solid) and the cnoidal-like wave solution of (1) corresponding
to measured parameters (dashed).

To generate a DB, the cnoidal-like carrier is set to overtake the soliton. So
that the carrier fully absorbs the soliton, it is essential to increase the
mean flow rate of the carrier relative to the soliton. This approach is
reported in Fig. 2a,b where, now, soliton-cnoidal-like wave interaction
results in a depression defect, a DB. For each fixed $z$, we define the DB
amplitude $a_{\rm DB}(z)$ as the difference between the maximum of the carrier
and the minimum of the DB upper envelope. The measured averaged carrier
parameters are given in Fig. 2c and $\tilde{k}=1.95\pm 0.03\text{cm}^{-1}$. In
contrast to the BB case, the carrier frequency shown in Fig. 2c increases up
the conduit. The carrier with a larger mean than the soliton leads to more
carrier fluctuations than in the BB case. Near the top of the conduit, we
observe a 5% increase in the carrier mean and a corresponding change in other
DB properties. Generally, the injected interior fluid is suspended at the top
after rising through the exterior fluid. We attribute DB modulations for
$\tilde{z}>110$cm to the slow diffusion of interior fluid from the top, which
lowers the exterior fluid’s density $\rho_{e}$. The observed 5% mean increase
can be explained by a correspondingly small 2.5% decrease in the exterior to
interior density difference via the Hagen-Poiseuille law in which
$\overline{A}\propto(\rho_{e}-\rho_{i})^{-1/2}$ [46]. Note that a similar mean
increase near the top is observed for the BB in Fig. 1c. Despite this,
propagation of the DB is robust over a large portion of the conduit
$\tilde{z}\in[80,110]$cm where carrier and DB fluctuations are modest.

The measured DB speed is $\tilde{v}_{\rm DB}=0.27\pm 0.01$cm/s over
$\tilde{z}\in[80,110]\text{cm}$ while the carrier phase speed is larger
$\tilde{v}_{ph}=0.45\pm 0.01$cm/s, a characteristic feature of KdV DBs [1,
39]. Figure 2d shows spatial, temporal slices and the determination of the
carrier phase shift. The phase shifts are $\Delta\theta_{\rm DBz}=0.24\pi$ in
space and $\Delta\theta_{\rm DBt}=0.25\pi$ in time.

Figure 3: Experiment (light gray) compared with simulation of the conduit
equation (1) (black) with initial conditions from experiment, and the KdV
traveling breather solution (blue). (a) BB from Fig. 1d. (b) DB from Fig. 2d.
The top four panels in (a,b) present the evolution over one period.

While qualitative features of observed BBs and DBs agree with KdV theory, KdV
is a quantitative model only for small amplitude $a\ll 1$, long waves $k\ll 1$
subject to the dominant balance $a\sim k^{2}$, where $a$ and $k$ are the
nondimensional amplitude and wavenumber of the carrier, respectively. Since
$a\gtrsim 1$ and $k\approx 0.4$ in our experiments, the dominant balance and
smallness conditions are not satisfied. Consequently, we perform numerical
simulations of the strongly nonlinear conduit equation (1) with periodic
boundary conditions and an initial condition extracted from the time sliced
experimental sections in Figs. 1d and 2d. Figure 3 depicts a comparison of the
simulations, experiment, and a KdV traveling breather with fitted carrier and
phase shift in the comoving reference frame [47]. In Fig. 3a, an oscillation
period of the BB ($\tilde{T}_{\rm BB}=2\pi/(\tilde{v}_{\rm
BB}\tilde{k}-\tilde{\omega})\approx 9$s) is shown in the top four panels with
the, much smaller than observation, KdV BB solution. Viewed as a soliton-
soliton lattice interaction, the BB exhibits the unimodal interaction geometry
according to the Lax categories of two-soliton interactions [57, 58, 59]. For
longer times, the simulation and experiment agree very well, demonstrating the
stability, robustness, and accuracy of conduit BBs relative to observation.
Small carrier discrepancies are due to unphysical periodic boundary conditions
and the acknowledged diffusion issues. The DB in Fig. 3b has period of
oscillation $\tilde{T}_{\rm DB}=2\pi/(\tilde{\omega}-\tilde{v}_{\rm
DB}\tilde{k})\approx 17$s, depicted in the top four panels. The dimensional
time scale $T$ from eq. (1) used here is 5% smaller than the nominal, measured
$T$ extracted from separate linear wave measurements [47]. We attribute this
small discrepancy to interior fluid diffusion and limitations of the conduit
equation as a model [41]. The KdV DB solution differs from observation, which
exhibits a bimodal interaction according to the Lax categories [59]. The DB
simulation closely tracks experiment subject to some carrier discrepancies.

Figure 4: Observation of traveling breather scattering on a fixed carrier
background. (a,c,e) Space-time contour revealing conduit dynamics. (b,d,f)
Traveling breather trajectories (blue and red solid) in the $z$-$t$ plane with
linear fits (dashed) and spatial shifts $\Delta\tilde{z}_{j}$. Subfigures
provide representative spatial profiles at specific times. (a,b) BB-BB
scattering. (c,d) DB-DB scattering. (e,f) BB-DB scattering.

The reliable creation of single BBs and DBs in viscous fluid conduits enables
the investigation of traveling breather scattering. We create multiple
traveling breathers on a fixed carrier background by nonlinear superposing the
periodic wavetrain with two solitons of differing amplitude. Figure 4a depicts
the scattering of two BBs. The BB created from the larger-amplitude soliton
travels at a faster speed and overtakes the smaller amplitude BB. Both BBs are
observed to propagate coherently before and after interaction, maintaining
essentially the same shape. Figure 4b tracks BB peak trajectories. We
separately fit the trajectories before and after interaction with linear
functions to derive BB speeds $\tilde{v}_{j}$ and spatial shifts
$\Delta\tilde{z}_{j}$ ($j=1,2$) occurring as a result of the interaction (see
Tab. 1). For each BB, the speeds pre and post interaction are the same, within
very small error tolerances, indicating that BB-BB scattering is physically
elastic. The faster (slower) BB experiences a positive (negative) spatial
shift.

Figure 4c presents DB scattering. The larger amplitude DB originating from the
larger amplitude soliton moves faster than the smaller amplitude, slower DB.
Trajectories in Fig. 4d trace the positions of the DB upper envelope minima.
Again, DB spatial shifts are observed and DB speeds pre and post interaction
are conserved (Tab. 1). The faster (slower) DB undergoes a positive (negative)
spatial shift.

Since all BBs (DBs) are faster (slower) than the carrier, a BB must overtake a
DB for them to interact as shown in Fig. 4e. The carrier must have finite
spatial extent in order to concurrently overtake and be overtaken by a soliton
to create a BB and DB. The speed measurements of BB, DB trajectories (Tab. 1)
depicted in Fig. 4f again indicate that the scattering is physically elastic.
The DB spatial shift is negative but the positive BB spatial shift is an order
of magnitude smaller, hence the BB is scarcely affected by the DB. In all
cases of observed traveling breather scattering, we find that they are
physically elastic and faster (slower) traveling breathers experience a
positive (negative) spatial shift due to the interaction. This observation is
consistent with KdV traveling breather scattering characterized in [1, 38].

In Fig. 4, we launch traveling breather pairs on the same carrier background.
By considering each traveling breather independently, we can infer their
nonlinear dispersion relation. In Fig. 4a, individual BB phase shifts after
interaction are $\Delta\theta_{\rm BB1}=-1.08\pi$ and $\Delta\theta_{\rm
BB2}=-0.85\pi$. Note that it is standard for a phase shift below $-\pi$ to be
shifted to $(0,\pi)$ by adding a $2\pi$ period. However, by maintaining the
phase shift sign, we can deduce that an increasing BB speed corresponds to an
increasing BB amplitude and decreasing BB phase shift. This is consistent with
the KdV BB dispersion relation [39, 47]. In Fig. 4c, we measure the DB phase
shifts post interaction to be $\Delta\theta_{\rm DB1}=0.61\pi$ and
$\Delta\theta_{\rm DB2}=0.52\pi$. An increasing DB speed corresponds to an
increasing DB amplitude and phase shift. For KdV DBs, there are two branches
of the DB nonlinear dispersion relation, slow and fast associated with
positive and negative phase shifts, respectively [39, 47]. Our experimental
observations are qualitatively consistent with the nonlinear dispersion
relation’s slow branch of KdV DBs.

| Traveling Breather 1 | Traveling Breather 2
---|---|---
$\tilde{v}_{1}$ before (cm/s) |  $\tilde{v}_{1}$ after (cm/s) |  $\Delta\tilde{z}_{1}$ (cm) |  $\tilde{v}_{2}$ before (cm/s) |  $\tilde{v}_{2}$ after (cm/s) |  $\Delta\tilde{z}_{2}$ (cm)
BB-BB | $0.65\pm 0.01$ | $0.64\pm 0.01$ | 1.54 | $0.55\pm 0.01$ | $0.53\pm 0.01$ | -2.69
DB-DB | $0.25\pm 0.01$ | $0.24\pm 0.01$ | 1.75 | $0.20\pm 0.01$ | $0.19\pm 0.01$ | -1.97
BB-DB | $0.61\pm 0.01$ | $0.62\pm 0.01$ | 0.15 | $0.22\pm 0.01$ | $0.21\pm 0.01$ | -1.22
Table 1: Speed and spatial shift measurements of traveling breather
scattering.

This work experimentally verifies that traveling breathers in a continuum
system can be generated from the interaction of solitons and periodic cnoidal-
like traveling waves. While continuum traveling breather theory has been
developed for integrable systems, our observations and numerical simulations
of the two-fluid system suggest that traveling breathers exist in non-
integrable systems as well. Traveling breathers appear to be natural,
practical and common phenomena that occur in continuum environments
accompanied by a possibly large amplitude oscillatory background,
distinguishing them from localized breathers. While spectral data analysis
techniques have been developed to extract this type of information from
complex signals using IST-based nonlinear Fourier analysis [60, 61], we have
physically realized this nonlinear superposition principle with a simple
method. The scattering of traveling breathers is observed to be physically
elastic (energy conserving) while exhibiting a spatial shift. This work lays
the foundation for exploring traveling breathers in other continuum systems
such as optics, fluids, condensed matter and anywhere nonlinear waves occur.

The authors would like to thank the Isaac Newton Institute for Mathematical
Sciences, Cambridge, for support and hospitality during the programme
Dispersive hydrodynamics: mathematics, simulation and experiments, with
applications in nonlinear waves, where work on this paper was discussed. This
work was supported by NSF DMS-1816934.

## References

  * Kuznetsov and Mikhailov [1975] E. A. Kuznetsov and A. V. Mikhailov, Stability of stationary waves in nonlinear weakly dispersive media, Sov. Phys. JETP 40, 855 (1975), translation from Russian of E.A. Kuznetsov, A.V. Mikhailov, ZhETF, 67, 1717 (1974).
  * Ablowitz  _et al._ [1973] M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, Method for solving the sine-gordon equation, Physical Review Letters 30, 1262 (1973).
  * Ablowitz and Segur [1981] M. J. Ablowitz and H. Segur,  _Solitons and the inverse scattering transform_ (SIAM, 1981).
  * Kudryavtsev [1975] A. Kudryavtsev,  _Solitonlike solutions for a Higgs scalar field_ , Tech. Rep. (Institute of Theoretical and Experimental Physics, 1975).
  * Ablowitz  _et al._ [1979] M. J. Ablowitz, M. D. Kruskal, and J. Ladik, Solitary wave collisions, SIAM Journal on Applied Mathematics 36, 428 (1979).
  * Campbell  _et al._ [1983] D. K. Campbell, J. F. Schonfeld, and C. A. Wingate, Resonance structure in kink-antikink interactions in $\varphi$4 theory, Physica D: Nonlinear Phenomena 9, 1 (1983).
  * Segur and Kruskal [1987] H. Segur and M. D. Kruskal, Nonexistence of small-amplitude breather solutions in phi 4 theory, Physical review letters 58, 747 (1987).
  * Soffer and Weinstein [1999] A. Soffer and M. I. Weinstein, Resonances, radiation damping and instability in hamiltonian nonlinear wave equations, Inventiones mathematicae 136, 9 (1999).
  * Boyd [1990] J. P. Boyd, A numerical calculation of a weakly non-local solitary wave: the $\phi$ 4 breather, Nonlinearity 3, 177 (1990).
  * Kalisch  _et al._ [2022] H. Kalisch, M. A. Alejo, A. J. Corcho, and D. Pilod, Breather Solutions to the Cubic Whitham Equation (2022), arxiv:2201.12074 [nlin, physics:physics] .
  * Pethiyagoda  _et al._ [2017] R. Pethiyagoda, S. W. McCue, and T. J. Moroney, Spectrograms of ship wakes: Identifying linear and nonlinear wave signals, J. Fluid Mech. 811, 189 (2017).
  * Chabchoub  _et al._ [2019] A. Chabchoub, K. Mozumi, N. Hoffmann, A. V. Babanin, A. Toffoli, J. N. Steer, T. S. van den Bremer, N. Akhmediev, M. Onorato, and T. Waseda, Directional soliton and breather beams, Proc. Natl. Acad. Sci. 116, 9759 (2019).
  * Grimshaw  _et al._ [2017] R. Grimshaw, J. C. B. da Silva, and J. M. Magalhaes, Modelling and observations of oceanic nonlinear internal wave packets affected by the Earth’s rotation, Ocean Modelling 116, 146 (2017).
  * Mollenauer  _et al._ [1980] L. F. Mollenauer, R. H. Stolen, and J. P. Gordon, Experimental Observation of Picosecond Pulse Narrowing and Solitons in Optical Fibers, Phys. Rev. Lett. 45, 1095 (1980).
  * Mandelik  _et al._ [2003] D. Mandelik, H. S. Eisenberg, Y. Silberberg, R. Morandotti, and J. S. Aitchison, Observation of Mutually Trapped Multiband Optical Breathers in Waveguide Arrays, Phys. Rev. Lett. 90, 253902 (2003).
  * Boyd [2013] R. W. Boyd,  _Nonlinear Optics_ (Academic Press, 2013).
  * Di Carli  _et al._ [2019] A. Di Carli, C. D. Colquhoun, G. Henderson, S. Flannigan, G.-L. Oppo, A. J. Daley, S. Kuhr, and E. Haller, Excitation Modes of Bright Matter-Wave Solitons, Phys. Rev. Lett. 123, 123602 (2019).
  * Luo  _et al._ [2020] D. Luo, Y. Jin, J. H. V. Nguyen, B. A. Malomed, O. V. Marchukov, V. A. Yurovsky, V. Dunjko, M. Olshanii, and R. G. Hulet, Creation and Characterization of Matter-Wave Breathers, Phys. Rev. Lett. 125, 183902 (2020).
  * Kharif  _et al._ [2009] C. Kharif, E. N. Pelinovsky, and A. Slunyaev,  _Rogue Waves in the Ocean_ , Advances in Geophysical and Environmental Mechanics and Mathematics (Springer, Berlin, 2009).
  * Kibler  _et al._ [2010] B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev, and J. M. Dudley, The Peregrine soliton in nonlinear fibre optics, Nature Phys 6, 790 (2010).
  * Chabchoub  _et al._ [2011] A. Chabchoub, N. Hoffmann, and N. Akhmediev, Rogue wave observation in a water wave tank, Physical Review Letters 106, 204502 (2011).
  * Närhi  _et al._ [2016] M. Närhi, B. Wetzel, C. Billet, S. Toenger, T. Sylvestre, J.-M. Merolla, R. Morandotti, F. Dias, G. Genty, and J. M. Dudley, Real-time measurements of spontaneous breathers and rogue wave events in optical fibre modulation instability, Nat Commun 7, 13675 (2016).
  * Tikan  _et al._ [2022] A. Tikan, F. Bonnefoy, G. Roberti, G. El, A. Tovbis, G. Ducrozet, A. Cazaubiel, G. Prabhudesai, G. Michel, F. Copie, E. Falcon, S. Randoux, and P. Suret, Prediction and manipulation of hydrodynamic rogue waves via nonlinear spectral engineering, Phys. Rev. Fluids 7, 054401 (2022).
  * Wu [2011] M. Wu, Nonlinear spin waves in magnetic film feedback rings, in  _Solid State Physics_ , Vol. 62, edited by R. E. Camley and R. L. Stamps (Academic Press, Burlington, 2011) pp. 163–224.
  * Flach and Gorbach [2008] S. Flach and A. V. Gorbach, Discrete breathers — Advances in theory and applications, Phys. Rep. 467, 1 (2008).
  * Boechler  _et al._ [2010] N. Boechler, G. Theocharis, S. Job, P. G. Kevrekidis, M. A. Porter, and C. Daraio, Discrete breathers in one-dimensional diatomic granular crystals, Physical review letters 104, 244302 (2010).
  * Iooss and James [2005] G. Iooss and G. James, Localized waves in nonlinear oscillator chains, Chaos: An Interdisciplinary Journal of Nonlinear Science 15, 015113 (2005).
  * James [2018] G. James, Travelling breathers and solitary waves in strongly nonlinear lattices, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20170138 (2018).
  * Sato and Sievers [2007] M. Sato and A. J. Sievers, Driven localized excitations in the acoustic spectrum of small nonlinear macroscopic and microscopic lattices, Physical review letters 98, 214101 (2007).
  * Kim  _et al._ [2015] E. Kim, F. Li, C. Chong, G. Theocharis, J. Yang, and P. G. Kevrekidis, Highly nonlinear wave propagation in elastic woodpile periodic structures, Physical review letters 114, 118002 (2015).
  * Belokolos  _et al._ [1994] E. D. Belokolos, A. I. Bobenko, V. Z. Enolski, A. R. Its, and V. B. Matveev,  _Algebro-Geometric Approach to Nonlinear Integrable Equations_ (Springer, New York, 1994).
  * Gesztesy and Svirsky [1995] F. Gesztesy and R. Svirsky,  _M-KdV Solitons on the Background of Quasi-Periodic Finite-Gap Solutions_ (American Mathematical Society, Providence, UNITED STATES, 1995).
  * Hu  _et al._ [2012] X.-R. Hu, S.-Y. Lou, and Y. Chen, Explicit solutions from eigenfunction symmetry of the korteweg–de vries equation, Physical Review E 85, 056607 (2012).
  * Lou [2015] S. Lou, Consistent riccati expansion for integrable systems, Studies in Applied Mathematics 134, 372 (2015).
  * Takahashi [2016] D. A. Takahashi, Integrable model for density-modulated quantum condensates: Solitons passing through a soliton lattice, Phys. Rev. E 93, 062224 (2016).
  * Kuznetsov [2017] E. A. Kuznetsov, Fermi–pasta–ulam recurrence and modulation instability, JETP letters 105, 125 (2017).
  * Nakayashiki [2020] A. Nakayashiki, One Step Degeneration of Trigonal Curves and Mixing of Solitons and Quasi-Periodic Solutions of the KP Equation, in  _Geom. Methods Phys. XXXVIII_ , Trends in Mathematics, edited by P. Kielanowski, A. Odzijewicz, and E. Previato (Springer International Publishing, Cham, 2020) pp. 163–186.
  * Bertola  _et al._ [2022] M. Bertola, R. Jenkins, and A. Tovbis, Partial degeneration of finite gap solutions to the korteweg-de vries equation: soliton gas and scattering on elliptic background, arXiv preprint arXiv:2210.01350 (2022).
  * Hoefer  _et al._ [2023] M. A. Hoefer, A. Mucalica, and D. E. Pelinovsky, Kdv breathers on a cnoidal background, J. Phys. A (2023), to appear.
  * Krichever [1975] I. M. Krichever, Potentials with zero coefficient of reflection on a background of finite-zone potentials, Functional Analysis and Its Applications 9, 161 (1975).
  * Mao and Hoefer [2023] Y. Mao and M. A. Hoefer, Experimental investigations of linear and nonlinear periodic travelling waves in a viscous fluid conduit, Journal of Fluid Mechanics 954, A14 (2023).
  * Martin  _et al._ [2009] J. Martin, N. Rakotomalala, D. Salin, L. Talon, _et al._ , Convective/absolute instability in miscible core-annular flow. part 1: Experiments, Journal of fluid mechanics 618, 305 (2009).
  * Selvam  _et al._ [2009] B. Selvam, L. Talon, L. Lesshafft, and E. Meiburg, Convective/absolute instability in miscible core-annular flow. part 2. numerical simulations and nonlinear global modes, Journal of Fluid Mechanics 618, 323 (2009).
  * Scott  _et al._ [1986] D. R. Scott, D. J. Stevenson, and J. A. Whitehead, Observations of solitary waves in a viscously deformable pipe, Nature 319, 759 (1986).
  * Olson and Christensen [1986] P. Olson and U. Christensen, Solitary wave propagation in a fluid conduit within a viscous matrix, Journal of Geophysical Research: Solid Earth 91, 6367 (1986).
  * Lowman and Hoefer [2013] N. K. Lowman and M. A. Hoefer, Dispersive hydrodynamics in viscous fluid conduits, Physical Review E 88, 023016 (2013).
  * [47] See Supplemental Material at — for additional mathematical and experimental details.
  * Barcilon and Richter [1986] V. Barcilon and F. M. Richter, Nonlinear waves in compacting media, Journal of Fluid Mechanics 164, 429 (1986).
  * Stubblefield  _et al._ [2020] A. G. Stubblefield, M. Spiegelman, and T. T. Creyts, Solitary waves in power-law deformable conduits with laminar or turbulent fluid flow, Journal of Fluid Mechanics 886 (2020).
  * Harris  _et al._ [2006] S. E. Harris, P. A. Clarkson, _et al._ , Painlevé analysis and similarity reductions for the magma equation, SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 2, 068 (2006).
  * Maiden and Hoefer [2016] M. D. Maiden and M. A. Hoefer, Modulations of viscous fluid conduit periodic waves, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472, 20160533 (2016).
  * Helfrich and Whitehead [1990] K. R. Helfrich and J. A. Whitehead, Solitary waves on conduits of buoyant fluid in a more viscous fluid, Geophysical & Astrophysical Fluid Dynamics 51, 35 (1990).
  * Maiden  _et al._ [2016] M. D. Maiden, N. K. Lowman, D. V. Anderson, M. E. Schubert, and M. A. Hoefer, Observation of dispersive shock waves, solitons, and their interactions in viscous fluid conduits, Physical review letters 116, 174501 (2016).
  * Gavrilyuk and Shyue [2021] S. Gavrilyuk and K.-M. Shyue, Singular solutions of the bbm equation: analytical and numerical study, Nonlinearity 35, 388 (2021).
  * Maiden  _et al._ [2018] M. D. Maiden, D. V. Anderson, N. A. Franco, G. A. El, and M. A. Hoefer, Solitonic dispersive hydrodynamics: theory and observation, Physical review letters 120, 144101 (2018).
  * Whitham [1984] G. B. Whitham, Comments on periodic waves and solitons, IMA J. Appl. Math. 32, 353 (1984).
  * Lax [1968] P. D. Lax, Integrals of nonlinear equations of evolution and solitary waves, Commun. Pure Appl. Math. 21, 467 (1968).
  * Craig  _et al._ [2006] W. Craig, P. Guyenne, J. Hammack, D. Henderson, and C. Sulem, Solitary water wave interactions, Phys. Fluids 18, 057106 (2006).
  * Lowman  _et al._ [2014] N. K. Lowman, M. A. Hoefer, and G. A. El, Interactions of large amplitude solitary waves in viscous fluid conduits, Journal of fluid mechanics 750, 372 (2014).
  * Costa  _et al._ [2014] A. Costa, A. R. Osborne, D. T. Resio, S. Alessio, E. Chrivì, E. Saggese, K. Bellomo, and C. E. Long, Soliton Turbulence in Shallow Water Ocean Surface Waves, Phys. Rev. Lett. 113, 108501 (2014).
  * Brühl  _et al._ [2022] M. Brühl, P. J. Prins, S. Ujvary, I. Barranco, S. Wahls, and P. L. F. Liu, Comparative analysis of bore propagation over long distances using conventional linear and KdV-based nonlinear Fourier transform, Wave Motion 111, 102905 (2022).