# Macroscopic Magnetic Monopoles in a 3D-Printed Mechano-Magnet

H. A. Teixeira  Laboratory of Spintronics and nanomagnetism (LabSpiN),
Departamento de Física, Universidade Federal de Viçosa, Viçosa, 36570-900,
Minas Gerais, Brazil  M. F. Bernardo  Laboratory of Spintronics and
nanomagnetism (LabSpiN), Departamento de Física, Universidade Federal de
Viçosa, Viçosa, 36570-900, Minas Gerais, Brazil  M. D. Saccone
<EMAIL_ADDRESS>Theoretical Division (T4) and Center for Nonlinear Studies,
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA  F.
Caravelli<EMAIL_ADDRESS>Theoretical Division (T4), Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  C. Nisoli<EMAIL_ADDRESS>Theoretical Division (T4), Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, USA  C. I. L. Araujo<EMAIL_ADDRESS>Laboratory of Spintronics
and nanomagnetism (LabSpiN), Departamento de Física, Universidade Federal de
Viçosa, Viçosa, 36570-900, Minas Gerais, Brazil

###### Abstract

The notion of magnetic monopoles has puzzled physicists since the introduction
of Maxwell’s Equations and famously Dirac had hypothesized them in the context
of quantum mechanics. While they have proved experimentally elusive as
elementary particles, the concept has come to describe excitations or
topological defects in various material systems, from liquid crystals, to Hall
systems, skyrmion lattices, and Bose-Einstein condensate. Perhaps the most
versatile manifestation of magnetic monopoles as quasiparticles in matter has
been in so-called spin ice materials. There, they represent violations of the
ice rule, carry a magnetic charge, and can move freely unbound. We have built
a mechano-magnet realized via 3D-printing, that consists of mechanical rotors
on which macroscopic magnets can pivot. By controlling the relative height of
the rotors we can achieve different regimes for magnetic monopoles, including
the free monopole state. We then explore their driven dynamics under field. In
the future, integration of our proof of principle in an elastic matrix can
lead to novel macroscopic mechano-magnetic materials, to explore unusual
piezomagnetism and magnetostriction, with applications to actuators and soft-
robotics.

##  I Introduction

While electrical dipoles can be broken into separable, opposite charges,
breaking a magnetic dipole does not separate its north and south poles. It
results instead in two smaller dipoles. The possibility of elementary,
quantum, magnetic monopoles has been explored theoretically by Dirac Dirac
(1931), but no such particle was ever found Preskill (1984). The concept,
however, has proved useful in describing emergent quasi-particles in liquid
crystals Chuang  _et al._ (1991), skyrmion lattices Milde  _et al._ (2013),
Bose-Einstein condensates Ray  _et al._ (2014), and in the context of the
anomalous Hall effect Fang  _et al._ (2003). The notion of a magnetic monopole
is perhaps the closest to the classical intuition in a special class of
frustrated binary magnets called spin ices, either crystal grown Ramirez  _et
al._ (1999); Bramwell and Gingras (2001); Fennell  _et al._ (2009); Ryzhkin
(2005); Castelnovo  _et al._ (2008); Morris  _et al._ (2009); Bramwell  _et
al._ (2009); Bovo  _et al._ (2013) or artificially fabricated at the
nanoscale, the so-called “artificial spin ices” Wang  _et al._ (2006); Skjærvø
_et al._ (2019); Mól  _et al._ (2009); Mengotti  _et al._ (2010); Ladak  _et
al._ (2011); Nascimento  _et al._ (2012); Perrin  _et al._ (2016); Farhan  _et
al._ (2019); Goryca  _et al._ (2021); Ostman and et al. (2018); Caravelli
(2020). These magnetic materials are made of interacting moments that can be
described as binary variables arranged along the edges of a lattice. The
moments are atomic in rare earth pyrochlores, and have sizes of tens or
hundreds of nanometers in nanofabricated artificial spin ices.

We show here for the first time that magnetic monopoles can be realized at the
macroscale, in a mechano-magnet made of dipoles pivoting on mechanical rotors,
shown in Fig. 1, a method pioneered before for kagome spin ice Mellado  _et
al._ (2012) and non-degenerate square ice Gonçalves  _et al._ (2020). We
develop the approach for a square lattice, and incorporate the relative offset
Möller and Moessner (2006) known to achieve degeneracy and thus free monopoles
Chern  _et al._ (2014); Perrin  _et al._ (2016); Farhan  _et al._ (2019).

Figure 1: The system. a) Design of the 3D-printed rotor, conceived to act as a
dipole in a macroscopic artificial spin ice system. Each dipole is composed of
a main axis of 15mm length, able to rotate in the polar direction, and two
neodymium magnets of 5mm diameter. The height of the axis can be changed in
three offsets, in multiples of $5mm$. b) 84 rotors arrayed along the geometry
of a square spin ice, fit at the center of an Helmholtz coil of $20cm$ radius.
The inset shows the evolution of magnetic field in function of time in the
demagnetization protocol. c) The observed vertex configurations. $T_{1},T_{2}$
obey the ice rule (2-in/2-out), of which $T_{2}$ is magnetized along the blue
line. The $T_{3}$ configuration violates the ice rule and carries a charge
$\pm 2$.

To understand our system, consider the model degenerate square spin ice of
Fig. 2. The low energy state is degenerate, and corresponds to disordered
configurations of moments that obey the so-called ice-rule Bernal and Fowler
(1933); Pauling (1935) in which two moments point in, two point out of each
vertex. If we define as charge of a vertex the number of moments pointing in
minus the number of moments pointing out, then ice-rule vertices have zero
magnetic charge. But flipping one moment creates a pair of opposite $\pm 2$
charges, or magnetic monopoles Ryzhkin (2005); Castelnovo  _et al._ (2008),
that break the ice rule. Further flips can separate them without more
violations of the ice rule and thus without an energy cost Castelnovo  _et
al._ (2008). Note that these monopoles are in fact north and south poles of a
long floppy dipole of sort, a magnetized line of moments (e.g. in grey in Fig.
2) often called “Dirac string”. However, because the manifold is disordered,
other strings of head-to-toe moments starting in a monopole and ending in
another could be equally chosen as Dirac strings (Fig. 2d), which are
therefore not uniquely defined. This leads to the interpretation of the low
energy manifold in terms of free monopole excitations of an ice rule obeying
spin vacuum.

Analogous magnetic monopoles had been proposed in pyrochlore spin ices where
they have also motivated a search for their current, or magnetricity Bramwell
_et al._ (2009). They were recently revealed in square artificial spin ice ice
at the nanoscale Perrin  _et al._ (2016); Farhan  _et al._ (2019). Early
realizations Wang  _et al._ (2006) possessed an ordered, antiferromagnetic
ground state because moments converging perpendicularly in a vertex interact
more strongly than those converging collinearly, thus lifting the degeneracy
of the ice rule. However, degeneracy can be regained by offsetting horizontal
and vertical magnets in the $z$ direction Möller and Moessner (2006); Chern
_et al._ (2014); Perrin  _et al._ (2016); Farhan  _et al._ (2019).

There are obviously enormous differences in system, geometry, dynamics, scale,
and interaction, between atomic magnetic moments in pyrochlore spin ice,
superparamagnetic nanoislands at thermal equilibrium, and macroscopic magnets
pivoting on a rotor. Yet, we show here that their behavior can be
conceptualized in similar ways. We implement the square ice geometry with
magnetic rotors and an offset perpendicular to the plane. This allows us to
explore in a macroscopic system the three different structural phases of
square ice–the antiferromagnetic state, the degenerate state, and the line
state–and the driven dynamics of their monopoles.

Figure 2: Schematics of magnetic monopole excitations as violations of the ice
rule. a) An example of ice rule configuration for Square Spin Ice: moments are
represented by arrows; in each vertex, two moments point in, two point out.
The number of such configuration scales exponentially with the number of
moments, leading to a finite Pauling entropy per moment. b) Flipping one
moment (yellow) in an ice rule configuration breaks the ice rule, and creates
a monopole pair: a $q=+2$ monopole corresponding to three moments pointing in
and one pointing out (blue circle), and a $q=-2$ monopole corresponding to
three moments pointing out and one pointing in (red circle). c) Further flips
of moments can separate the monopoles, leaving behind a so-called ”Dirac
String” (yellow spins), which is however not univocally determined. d) A
different possible Dirac string is marked in yellow, the previous one in light
yellow. Note that flipping any entire Dirac string annihilates the monopole
pair. Figure 3: AC Annealing of the mechano-magnetic spin ice. a) Snapshot of
experimental setup. c) Vertex configurations observed after the AC
demagnetization protocol. b) Fractions of vertex types obtained after the
demagnetization protocol as a function of the offset show $T_{1}$ vertices
being promoted at low values of $h/a$, $T_{2}$ vertices being promoted at high
values of $h/a$. c) Fractions of vertex types obtained via a mechanical
dynamical simulation of a much larger system agree qualitatively well with the
experimental results. d-f) Real space configurations after AC demagnetization
for offsets $h_{1}=5mm$ (d), $h_{2}=10mm$ (e) and $h_{3}=15mm$ (f).

##  II System

We have utilized a 3D printer to fabricate rotors equipped with neodymium
spherical magnets of 5mm diameter that are allowed to spin in the polar angle
direction, as depicted in Fig. 1a. The rotors are arranged pointing along the
edges of a square lattice of lattice constant $a=27mm$.

The pivoting rotors can be set at four different heights, differing by $5mm$:
$h_{0}=0,h_{1}=5mm$, $h_{2}=10mm$, and $h_{3}=15mm$. The rotors converging in
a vertex form four possible configurations, shown in Fig. 1c, and called
$T_{1},T_{2},T_{3},$ Wang  _et al._ (2006). $T_{1},T_{2}$ vertices obey the
ice rule, but $T_{1}$ vertices have zero net magnetic moment. $T_{3}$ vertices
corresponds to violations of the ice rule in the form of $q=\pm 3\mp 1=\pm 2$
monopoles. A fourth set of configurations, $T_{4}$ vertices, or monopoles of
charge $q=\pm 4$, correspond to all spin pointing in or all out, are highly
energetic and never observed.

We follow the general recipe of Ref. Möller and Moessner (2006) and offset the
vertically aligned moments with respect to the horizontally aligned ones, to
achieve different interaction strengths among perpendicular and collateral
rotors ($J_{1},J_{2}$ in Fig. 1c). For an offset of zero, $T_{1}$ vertex
configurations have the lowest energy, and thus the lowest energy of the array
corresponds to an antiferromagnetic tessellation of $T_{1}$ vertices. For very
high offset, the magnetized $T_{2}$ vertices are expected to be lowest in
energy, and the collective lowest energy state is expected to be a so-called
line state, as it corresponds to magnetization lines flowing through the
system Perrin  _et al._ (2016); King  _et al._ (2021); Nisoli (2021).

At some intermediate offset, by reasoning within a simplified nearest neighbor
model for simple binary moments, one expects that the ground state of the
system corresponds to the degenerate six-vertex model, and that its
excitations correspond to free monopoles Möller and Moessner (2006). While our
mechanical system is more complex than such model of binary spins, our results
confirm the same picture, as we shall show.

##  III Annealing via AC Demagnetization

For each of the four offsets, we prepared the system via an AC demagnetization
protocol Ke  _et al._ (2008); Rodrigues  _et al._ (2013). The external
magnetic field generated by a Helmholtz coil with radius $R=20cm$ decreases in
oscillatory steps of amplitude $B_{x}=$ 5 Oe starting from a field of 40 Oe.
Results obtained after the demagnetization processes are presented in Fig. 3.
In all cases ice-rule vertices are promoted, while monopoles are suppressed.

Figure 4: Driving Macroscopic monopoles by Field Inversion in an anti-
ferromagnetic spin ice (offset $h_{1}$). a) Magnetic moment measured in units
of rotors vs. applied field. b) Vertex statistics vs. applied field. Note the
proliferation of $T_{1}$ antiferromagnetic vertices. c) Real-space snapshots
of moment inversion. Yellow moments have flipped with respect to the previous
frame. Figure 5: Driving Macroscopic monopoles by Field Inversion in a line
state spin ice (offset $h_{3}$). a) Magnetic moment measured in units of
rotors vs. applied field. b) Vertex statistics vs. applied field. Note the
absence of $T_{1}$ configurations, as $T_{2}$s transform into monopoles
($T_{3}$) and then into $T_{2}$ configuration of opposite magnetic moment.

For zero offset we obtain a completely antiferromagnetic configuration, the
expected tessellation of $T_{1}$ vertices of zero net magnetic moment. At the
offset $h_{1}$, the energy difference between the $T_{2},T_{3}$ and the
$T_{1}$ vertex configurations is reduced, and we obtain antiferromagnetic
domains of $T_{1}$ vertices, as shown in the schematics of Fig. 3d, in which
the rotors are represented by arrows pointing in the direction of their
projection in the plane. Note that, unlike the binary variables of standard
spin ices, here the rotors have a continuous degree of freedom: their pivoting
angle. The domains are separated by domain walls, which we draw according to
the convention of Ref. Nisoli (2020) by joining together the magnetic moment
of $T_{2}$ vertices. These domains are also “Faraday Lines” Nisoli (2020)
entering negative monopoles and exiting positive ones. The obtained
configuration has a small net magnetization per moment
$Mx_{h1}={\color[rgb]{0,0,0}0.095}$, defined by assigning to each arrow a
moment $\pm 1$.

For the largest offset $h_{3}$ the lowest energy configuration corresponds to
$T_{2}$ vertices. Indeed, Fig. 3b,d shows that most vertices are in the
magnetized $T_{2}$ configuration. The Faraday lines no longer describe domain
walls, but a line state that covers almost the entire sample, which has now a
total magnetic moment $M_{x}=0.42$, $M_{y}=0$, as also seen previously for
nanomagnetic realizations Perrin  _et al._ (2016).

At some intermediate offset, the $T_{1}$ and $T_{2}$ configurations are
predicted to become degenerate. For ideal dipoles, this critical height is
estimated at $h/a\simeq 0.419$, but the finite length of the dipole lowers
that ratio. In our case, from $h_{c}=a-l/\sqrt{2}$, using the length of the
rotor $l=1.2mm$ and lattice constant $a=27mm$ return the offset $h\simeq
10mm$. corresponding to our $h_{2}$.

The results of the annealing at offset $h_{2}$ is shown in Fig. 3c. The sample
is fairly demagnetized, with $M_{x}=-0.047$, $M_{y}=0.143$. Note that the
ratio of relative occurrence among $T_{2}$ and $T_{1}$ vertices is $\simeq
1.66$, quite close to the ratio among their multiplicities $2=4/2$ (there are
four types of $T_{2}$ vertices, two types of $T_{1}$ vertices) and similar to
that observed in literature for nanomagnetic systems either demagnetized via
AC demagnetization or in thermal equilibrium Mól  _et al._ (2010); Perrin  _et
al._ (2016); Farhan  _et al._ (2019). Interestingly, the observed frequency of
vertex configurations corresponds to that found in a thermal state for
$T/J=1.1$, as it can be verified via a Metropolis Monte Carlo simulation of
nearest neighbors coupling $J$. This is not completely surprising, as it has
been shown before that AC demagnetization can lead to ensembles similar to the
thermal ones Nisoli  _et al._ (2007, 2010); Perrin  _et al._ (2019).

In Supplementary Materials, we show how to simulate our system within a driven
mechanical model of rotors. Our numerical results agree well with the
experimental findings, as shown in Fig. 3c.

Finally, we note that the geometric Gauss’s law is fulfilled in all samples.
If we call $\phi$ the flux of the magnetic moments inside the system, defined
as number of boundary spins pointing in minus those pointing out, than $\phi$
must be equal the net charge in the system, or $\phi=Q_{\text{tot}}$. That
number is zero for the offset of zero and $5mm$, the antiferromagnetic case,
and $-8$ in the other two cases.

Figure 6: Driving Macroscopic monopoles by Field Inversion in a degenerate
spin ice state (offset $h_{2}$). a) Magnetic moment measured in units of
rotors vs. applied field. b) Vertex statistics vs. applied field.

##  IV Monopoles Driven by Field

We explore monopole dynamics by driving them with an in-plane external
magnetic field in the $[1,0]$ direction. We saturate the material under a
field of $\sim 35$ Oe and then invert the field. When field is inverted
linearly, the behavior of the system proceeds by short bursts of collective
rearrangement of spins. Each time the system changes, we pause the field
inversion and collect data at that field. Therefore, in the data shown in Fig.
4, Fig. 5, and Fig. 6 the system remains unchanged as the field decreases
between the various recorded data points, which are then shown as real space
images.

In Fig. 4 we report results for the offset $h_{1}$ (antiferromagnetic spin
ice). For this case, the interaction is strong and a field of 35 Oe cannot
fully polarize it. Under inversion, at low negative field, $T_{2}$ vertices
transform first into monopoles, and then into $T_{1}$ vertices, generating
anti-ferromagnetic domains whose domain walls remain magnetized. Then as the
negative field grows, the kinetics becomes largely dominated by the movement
of domain walls. This effect can be seen in the curves for vertex
configurations (Fig. 4b) and also in the magnetization curve (Fig. 4a). The
latter is not as sharp as the other two cases (shown below), reflecting the
complexity of the vertex transformations involved in the inversion of the
overall moment. The $T_{2}$ vertices change their orientation via elaborate
set of “reactions”. Typically, for a vertex, the reactions pathway is
$T_{2}\to T_{3}\to T_{1}\to T_{3}\to T_{2}$ and cannot always be completed.
This leaves and abundance of magnetized monopoles and causes an asymmetry
between the initial state polarized to the right and the final state polarized
to the left. This asymmetry is absent in the other two cases.

In Fig. 5 we show the case of the largest offset $h_{3}$, whose low energy
manifold is a line state. At saturation horizontal spins are all pointing in
the same direction. Vertical spins form ferromagnetic rows of alternating
direction, due to their antiferromagnetic interactions. Because of the large
offset, and unlike in the $h_{1}$ case, $T_{1}$ vertices have now higher
energy than $T_{2}$. Thus $T_{1}$ vertices do not appear as an intermediate
step in the pathway of moment inversion. The moment inversion proceeds as
$T_{2}\to T_{3}\to T_{2}$, that is via creation of monopoles ($T_{3}$) and
their propagation through the system, which leaves in its wake a line of
$T_{2}$ vertices properly aligned along the field. Sometimes, entire rows of
moments flip. This can always be interpreted as a monopole “entering” the
system from the right boundary and being pushed out at the left boundary. The
monopoles, driven by the field, appear to propagate ballistically. Expectedly,
the positive monopoles mostly move along the direction of the field, the
negative ones in the opposite direction. The lack of the complex vertex-
reactions of the previous case makes the magnetization profile sharper (Fig.
5a). This allows for a reasonable, if approximated, determination of a
coercive field at around 10 Oe. The statistics of vertex configurations (Fig.
5b) looks therefore considerably flatter than in the other cases, with a small
maximum in monopole presence around the coercive field. At high negative
field, the monopoles are eliminated by pushing them out of the boundaries.

The case of offset $h_{2}$, or degenerate square ice, is intermediate between
the two others (Fig. 5). In the absence of a field, the $T_{1}$ and $T_{2}$
vertex configurations are energetically degenerate. However, in presence of a
field $T_{2}$ vertices aligned along the field possess lower energy than the
$T_{1}$ vertices, because of the Zeeman contribution from the field Gilbert
_et al._ (2015); Goryca  _et al._ (2021). During reversal, the energetic
advantage of the newly misaligned $T_{2}$ vertices is lost as they then oppose
the external field. This might explain why the $T_{1}$ configurations reach a
maximum at the coercive field (again at circa 10 Oe) where their energy is
lower than most $T_{2}$ vertices. The initial slope in relative magnetization
is large as the same ballistic monopole motion observed in the ferromagnetic
system can quickly reverse lines of spins, even while $T_{1}$ vertices allow
for more complex motion. There is a secondary leap in magnetization around 20
Oe as the population of $T_{1}$ vertices decrease while $T_{3}$ increase. This
is the latter portion of the vertex reactions described for the
antiferromagnetic spin ice, a $T_{1}\to T_{3}\to T_{2}$ process. The monopoles
may move more freely and, with a lesser slope than the initial reversal,
gradually convert away all the $T_{1}$ vertices while increasing the net
magnetic moment.

Finally in the panels $a$ of Figs. 4-6 we also show in red the simulated lines
of magnetization after an annealing that is sinusoidal with a time- decreasing
envelope. While there is no perfect agreement, as the real system is more
complex to simulate mechanically, the simulations show similar qualitative
differences among the three cases seen in experiments.

##  V Discussion

We have provided proof of principle that subtle magnetism of degenerate
frustrated manifolds and magnetic monopoles can be realized at the macroscale,
where macroscopic magnetic monopoles can be driven by external field. We have
found that models and ideas developed at the atomic or nanoscopic scale are
applicable to a much different macroscopic system. In particular, we have
reproduced the three theoretically predicted regimes already realized at the
nanoscale Wang  _et al._ (2006); Perrin  _et al._ (2016); Farhan  _et al._
(2019) and in the superconductive qubits of a quantum annealer King  _et al._
(2021), and similar monopole dynamics, but at the macroscale.

Because the constrained disorder of the degenerate square ice (realized by us
at offset $h_{2}$), this is another classical and macroscopic example of a
seemingly topological phase Henley (2010); Castelnovo and Chamon (2007), due
to a kinetically constrained dynamics. One key difference is that, however,
the real degrees of freedom of this system are continuous, and we did observe
in simulations and experiments fast cascade transitions which reabsorb
monopole pairs. Such fast cascades will be the subject of future experiments
and numerical analysis. It is in fact typically the case that rapid cascades
in system with a threshold dynamics leads to self-organized criticality and
punctuated equilibria Bak and Sneppen (1993). Avalanches were seen in
antiferromagnetic square ice Bingham  _et al._ (2021). Our system is an ideal
platform to study monopole-related avalanches at the macroscale.

Our work also show another instance of apparently thermal ensembles achieved
in an athermal system via external driving D’Anna  _et al._ (2003). In the
future it would be interesting to investigate relaxation procedures based on
mechanical rather than magnetic drives, such as vibro-fluidization. Finally,
in the context of the recent introduction of ideas from frustration and spin
ice in the field of mechanical metamaterals Meeussen  _et al._ (2019); Pisanty
_et al._ (2021); Merrigan  _et al._ (2021), our work can open new venues for
macroscopic magneto-mechanical systems, of unusual magnetostriction or
piezomagnetism.

## Acknowledgments

The authors thank CNPq, CAPES and FAPEMIG (Brazilian agencies) for financial
support. The work of F. Caravelli, C. Nisoli and M. Saccone was carried out
under the auspices of the NNSA of the U.S. DoE at LANL under Contract No. DE-
AC52-06NA25396, financed via LDRD-ER grant PRD20190195.

## References

  * Dirac (1931) P. A. M. Dirac, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 133, 60 (1931).
  * Preskill (1984) J. Preskill, Annual Review of Nuclear and Particle Science 34, 461 (1984).
  * Chuang  _et al._ (1991) I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science 251, 1336 (1991).
  * Milde  _et al._ (2013) P. Milde, D. Köhler, J. Seidel, L. Eng, A. Bauer, A. Chacon, J. Kindervater, S. Mühlbauer, C. Pfleiderer, S. Buhrandt, _et al._ , Science 340, 1076 (2013).
  * Ray  _et al._ (2014) M. W. Ray, E. Ruokokoski, S. Kandel, M. Möttönen, and D. Hall, Nature 505, 657 (2014).
  * Fang  _et al._ (2003) Z. Fang, N. Nagaosa, K. S. Takahashi, A. Asamitsu, R. Mathieu, T. Ogasawara, H. Yamada, M. Kawasaki, Y. Tokura, and K. Terakura, Science 302, 92 (2003).
  * Ramirez  _et al._ (1999) A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry, Nature 399, 333 (1999).
  * Bramwell and Gingras (2001) S. T. Bramwell and M. J. Gingras, Science 294, 1495 (2001).
  * Fennell  _et al._ (2009) T. Fennell, P. Deen, A. Wildes, K. Schmalzl, D. Prabhakaran, A. Boothroyd, R. Aldus, D. McMorrow, and S. Bramwell, Science 326, 415 (2009).
  * Ryzhkin (2005) I. Ryzhkin, Journal of Experimental and Theoretical Physics 101, 481 (2005).
  * Castelnovo  _et al._ (2008) C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).
  * Morris  _et al._ (2009) D. J. P. Morris, D. Tennant, S. Grigera, B. Klemke, C. Castelnovo, R. Moessner, C. Czternasty, M. Meissner, K. Rule, J.-U. Hoffmann, _et al._ , Science 326, 411 (2009).
  * Bramwell  _et al._ (2009) S. T. Bramwell, S. Giblin, S. Calder, R. Aldus, D. Prabhakaran, and T. Fennell, Nature 461, 956 (2009).
  * Bovo  _et al._ (2013) L. Bovo, J. Bloxsom, D. Prabhakaran, G. Aeppli, and S. Bramwell, Nature communications 4, 1 (2013).
  * Wang  _et al._ (2006) R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, Nature 439, 303 (2006).
  * Skjærvø  _et al._ (2019) S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Nature Reviews Physics 2, 1 (2019).
  * Mól  _et al._ (2009) L. Mól, R. Silva, R. Silva, A. Pereira, W. Moura-Melo, and B. Costa, Journal of Applied Physics 106, 063913 (2009).
  * Mengotti  _et al._ (2010) E. Mengotti, L. J. Heyderman, A. F. Rodríguez, F. Nolting, R. V. Hügli, and H.-B. Braun, Nat. Phys. 7, 68 (2010).
  * Ladak  _et al._ (2011) S. Ladak, D. Read, W. Branford, and L. Cohen, New Journal of Physics 13, 063032 (2011).
  * Nascimento  _et al._ (2012) F. S. Nascimento, L. A. S. Mol, W. A. Moura-Melo, and A. R. Pereira, New Journal of Physics 14, 115019 (2012).
  * Perrin  _et al._ (2016) Y. Perrin, B. Canals, and N. Rougemaille, Nature 540, 410 (2016).
  * 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).
  * Goryca  _et al._ (2021) M. Goryca, X. Zhang, J. Li, A. Balk, J. Watts, C. Leighton, C. Nisoli, P. Schiffer, and S. Crooker, Physical Review X 11, 011042 (2021).
  * Ostman and et al. (2018) E. Ostman and et al., Nature Phys. 14, 375–379 (2018).
  * Caravelli (2020) F. Caravelli, Eur. Phys. Lett. 130 (2020).
  * Mellado  _et al._ (2012) P. Mellado, A. Concha, and L. Mahadevan, Physical review letters 109, 257203 (2012).
  * Gonçalves  _et al._ (2020) R. Gonçalves, A. Gomes, R. Loreto, F. Nascimento, W. Moura-Melo, A. Pereira, and C. de Araujo, Journal of Magnetism and Magnetic Materials 502, 166471 (2020).
  * Möller and Moessner (2006) G. Möller and R. Moessner, Phys. Rev. Lett. 96, 237202 (2006).
  * Chern  _et al._ (2014) G.-W. Chern, C. Reichhardt, and C. Nisoli, Applied Physics Letters 104, 013101 (2014).
  * Bernal and Fowler (1933) J. Bernal and R. Fowler, The Journal of Chemical Physics 1, 515 (1933).
  * Pauling (1935) L. Pauling, Journal of the American Chemical Society 57, 2680 (1935).
  * King  _et al._ (2021) A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, and A. Lopez-Bezanilla, Science 373, 576 (2021).
  * Nisoli (2021) C. Nisoli, AIP Advances 11, 105313 (2021).
  * Ke  _et al._ (2008) X. Ke, J. Li, C. Nisoli, P. E. Lammert, W. McConville, R. Wang, V. H. Crespi, and P. Schiffer, Physical review letters 101, 037205 (2008).
  * Rodrigues  _et al._ (2013) J. Rodrigues, L. Mól, W. Moura-Melo, and A. Pereira, Applied Physics Letters 103, 092403 (2013).
  * Nisoli (2020) C. Nisoli, EPL (Europhysics Letters) 132, 47005 (2020).
  * Mól  _et al._ (2010) L. Mól, W. Moura-Melo, and A. Pereira, Physical Review B 82, 054434 (2010).
  * Nisoli  _et al._ (2007) C. Nisoli, R. Wang, J. Li, W. McConville, P. Lammert, P. Schiffer, and V. Crespi, Phys. Rev. Lett. 98, 217203 (2007).
  * Nisoli  _et al._ (2010) C. Nisoli, J. Li, X. Ke, D. Garand, P. Schiffer, and V. H. Crespi, Phys. Rev. Lett. 105, 047205 (2010).
  * Perrin  _et al._ (2019) Y. Perrin, B. Canals, and N. Rougemaille, Physical Review B 99, 224434 (2019).
  * Gilbert  _et al._ (2015) I. Gilbert, G.-W. Chern, B. Fore, Y. Lao, S. Zhang, C. Nisoli, and P. Schiffer, Physical Review B 92, 104417 (2015).
  * Henley (2010) C. L. Henley, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).
  * Castelnovo and Chamon (2007) C. Castelnovo and C. Chamon, Physical Review B 76, 174416 (2007).
  * Bak and Sneppen (1993) P. Bak and K. Sneppen, Physical Review Letters 71, 4083 (1993).
  * Bingham  _et al._ (2021) N. Bingham, J. Park, A. Simon, W. Zhu, J. Watts, X. Zhang, J. Batley, K. Dahmen, C. Leighton, and P. Schiffer, Bulletin of the American Physical Society 66 (2021).
  * D’Anna  _et al._ (2003) G. D’Anna, P. Mayor, A. Barrat, V. Loreto, and F. Nori, Nature 424, 909 (2003).
  * Meeussen  _et al._ (2019) A. S. Meeussen, E. C. Oguz, Y. Shokef, and M. van Hecke, Nature Physics 16, 307 (2019).
  * Pisanty  _et al._ (2021) B. Pisanty, E. C. Oğuz, C. Nisoli, and Y. Shokef, SciPost Physics 10, 136 (2021).
  * Merrigan  _et al._ (2021) C. Merrigan, C. Nisoli, and Y. Shokef, Physical Review Research 3, 023174 (2021).
  * Note (1) Experimental video available at https://www.youtube.com/watch?v=oxnfRsieK2A.
  * Note (2) We are happy to share the code upon request. Please write to<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
  * Note (3) Simulation video available at https://www.youtube.com/watch?v=sgVFy4ghFvs at the degeneracy point.

## Appendix

###  V.1 Experiment setup

The rotors presented in Figure 1 were developed by a conventional 3d printing
process using Polylactic acid material. The main halter total length of 12mm
is formed by two neodymium magnetic spheres with 5mm diameter each. This main
part is connected to the external arms by a conic groove, allowing rotation
with minimal friction. Three different distances among grooves were used to
perform the z-axis offsets h=5,10 and 15mm discussed in the paper, with the
halters aligned in y direction. The rotors were disposed in a square lattice,
at a distance of L=27mm, and a total length for the apparatus of 20cm. This
size was used in order to ensure the uniformity of the external field along
the experiment. The external field was applied by a commercial Helmholtz coil.

To anneal the sample, we employed an oscillatory demagnetization protocol,
starting with a current of 5A (the maximum allowed in the coil), corresponding
to an external field of 35Oe; we then decreased the magnetic field by steps of
5Oe per second, and took a picture of the apparatus. The reversal processes
were performed starting with system saturation, with a maximum field of 35Oe
along the x-axis. We then mapped the population density observed in each
picture to a snapshot, as in the paper. A video of the experiment is available
at a video repository. 111Experimental video available at
https://www.youtube.com/watch?v=oxnfRsieK2A.

###  V.2 Molecular dynamics simulations, annealing and monopoles movement

We performed molecular dynamics simulations, by numerically implementing the
force applied on each magnetic charge (approximating every magnetic rotor as a
pair of charges) and projecting onto the axis of rotation.222We are happy to
share the code upon request. Please write to<EMAIL_ADDRESS>or
<EMAIL_ADDRESS>We performed a direct numerical integration of the
equations via a Verlet method (e.g. the second order direct Euler
discretization), with an integration time step $\delta t=10^{-2}$, in the
absence of noise sources (zero temperature). The parameters used are those we
inferred to best match the experimental results, as we describe below. The
equations of motion for each rod were those derived in Mellado  _et al._
(2012) with the addition of an external field coupling

$\displaystyle I\frac{d^{2}\alpha_{i}}{dt^{2}}$ $\displaystyle=$
$\displaystyle T_{i}-\eta\frac{d\alpha_{i}}{dt}+\tau_{i}$ (1) $\displaystyle
T_{i}$ $\displaystyle=$
$\displaystyle\left(\frac{a\mu_{0}}{4\pi}\right)\sum_{j\neq
i}q_{i}q_{j}\hat{r}_{i}^{cm}\frac{\vec{r}_{ij}}{|\vec{r}_{ij}|^{3}}$ (2)
$\displaystyle\tau_{i}$ $\displaystyle=$
$\displaystyle\mu_{B}m(t)p_{i}^{\perp}\cdot(\vec{s}_{i}(\theta_{i})\times\hat{h});$
(3)

In the equations above, the total torque is comprised of torques between
nearest neighbor rotors, $T_{i}$, the kinetic friction coefficient is $\eta$,
and external field torque $\tau_{i}$. $\vec{r}_{ij}$ is the vector joining the
magnetic charges $q_{i}$ and $q_{j}$, while $\hat{r}_{i}^{cm}$ is unit vector
pointing towards the charge $q_{i}$ from the center of mass of a rotor. The
quantity $\tau_{i}$ is the torque force applied to the dipole, where
$p_{i}^{\perp}$ is the vector associated with the dipole rotation axis,
$\hat{h}$ is the direction of the external field, and $\vec{s}_{i}$ is the
vector associated to the dipole; $m(t)$ is the time depedent magnitude of the
field, and $\mu_{B}$ is the effective coupling between the dipole and external
field. We obtained a reasonable fit with the experimental results by choosing
$\eta=1$, $a=1$ (mm). We did not introduce static friction, which should
explain the difference between the numerical and experimental results in the
main text. $I$ was chosen so that in absence of forcing, the relaxation time,
controlled by $\eta/I\approx 1$ is such that the dipoles relax within a
timescale of $\approx 1$ second (a number of steps equal to $1/dt$), which is
consistent with independent experimental observations. We have allowed for a
staggered height of the monopoles on a checkerboard lattice (e.g. an island
has height 0 and the neighbors have height $h$). Thus, the distance between
two islands is given by $r_{ij}=\sqrt{d_{ij}^{2}+(h_{i}-h_{j})^{2}}$, where
$d_{ij}$ is the distance on the plane.

Figure 7: Numerical annealing of the macroscopic dipoles. The magnetization as
a function of the external fields (in units of $\mu_{B}=1$). We see that the
shape of the annealing resembles the experimental results. The magnetization
is $M$ is the one parallel to the external field $\vec{H}=m(t)\hat{H}$.

First, we annealed the sample with a decaying, oscillating field applied
diagonally as to affect both subsets of spins equally. We attempted to place
the system in the ground state (at $h=0$ this corresponds to a Type I
vertices). This can be seen in Fig. 7.

Figure 8: A snapshot of a simulation, for $h=0.48mm$(the degeneracy point) at
a fixed external magnetic field. Red and blue represent the negative and
positive T3 vertices, e.g. the monopoles. The lengths of the arrows represent
the projection of the plane of the simulated rotors.

A video of the annealing is available at an online repository for L=20, while
a snapshot is shown in Fig. 8.333Simulation video available at
https://www.youtube.com/watch?v=sgVFy4ghFvs at the degeneracy point. The value
of $q$ has been chosen to be $1.5$, in order to obtain a crossover of the
populations such that for $h\approx 0.55$ the populations of Type I and Type
II are identical, which is consistent with experimental results.

Figure 9: Population versus the field strength $m(t)$, $(m_{0}=1)$ for two
values of $h$. On the top we have the populations for $h=0$, where we see that
Type I and Type II cross in population as a function of the external field. On
the bottom, we plot the populations for $h=0.8$, where we see that Type II are
stable as a function of $h$. This is also consistent with the experimental
observations.

Once all the experimental parameters have been found, we can use molecular
dynamics to estimate more precisely the location of the degenerate point, as
function of $h$, in which the energies of the the Type I and Type II vertices
are equal. Similarly to the experiments, we have a crossover between Type I
and Type II at around $h=44$ (mm), with the presence of Type III vertices
increasing as a function of $h$. During the annealing procedure, the vertices
being excited differ depending on the height parameter. This can be seen in
Fig. 9, where different annealing curves for the vertex populations are
observed. For $h=0$, Type I and Type II populations cross as we swipe the
external field diagonally across the lattice. For larger values of $h$,
however, Type II vertices become more stable, indicating a crossover in
energies. These numerical simulations are also consistent with experimental
observations.

Figure 10: Tracked monopoles in the simulation of a $20\times 20$ lattice,
after an annealing procedure. We note that for $h=0$ there is the formation of
a higher number of monopoles. For $h=0.5$, there is a lower density of
monopoles, but these can move ballistically.

We studied numerically the evolution of the systems under a hysteresis loop,
comparable to the experiments depicted in Figures 4-6. In this case, the field
was applied parallel to one set of spins and perpendicular to the other to
match the experimental protocol while the height offset was set to $h=0.45,$
and 1 respectively to match the experimental vertex type populations. The
average in plane magnetization, as plotted in Figures 4-6, is a good fit to
the experimental values. A subjective feature relevant to our understanding of
these systems is that the degenerate system evolves in leaps at a time when a
threshold of external field is cleared. Since this feature is present in both
the experiment and simulations, it is relevant to consider as a new dynamical
feature of the system’s collective motion.

Finally, we tracked the movements of the Type III vertices, e.g. the
monopoles, comparing their movement both at $h=0$ and for $h=0.5$, e.g. near
the overlapping points. We tracked every single monopole in a simulation. If
no monopoles were created or annihilated, we could simply consider the closest
monopole in a successive frame to be the same monopole. However, if the
distance between closest monopole was greater than a lattice spacing, the
monopole was considered annihilated as no monopole traveled farther in a
frame. Finally, the remaining monopoles with no previous time pairs were
considered new monopoles. We plotted them over time in Fig. 10. Each colored
curve represents the position of a monopole. The first observation is that we
can in fact track the monopoles in continuous space rather than discrete: this
is due to the fact that as a particular dipole flips, the monopole
continuously moves from one site to another. For $h=0$, then, we note that the
monopoles are localized, e.g. they are created and destroyed without
significant motion from site to site and thus have a finite lifetime. In terms
of movement, for $h=0$ we observe only a handful of monopoles moving
ballistically. On the other hand, for $h=0.5$, while there is a lower density
of monopoles, these can move over longer ballistic trajectories, implying that
at this particular point the monopoles are less bound, as we should expect.